models of the Alpine Fault seismic cycle: stress transfer in the mid-crust

We examine the expected elastic and inelastic strain accumulation and release across the transpressive Alpine Fault of New Zealand using a ﬁnite element model, which incorporates a frictional fault and material behaviour in the crust based on geophysical measurements and rock mechanics experiments. A zone of localized ductile creeping is predicted beneath the part of the fault that slips seismically. Localized creeping results from (a) thermal weakening along the fault ramp and (b) enhanced stresses due to interaction with the shallow, brittle fault above. The creeping zone controls the shorter-wavelength feature seen in the GPS-determined surface velocities. Using the model we show that: (1) seismic faulting loads the creeping shear zone elastically, transferring stress to the mid-crust, and enhancing creep localization in the lower crust; (2) seismic loading and additional thermal weakening, due to long-term advection and exhumation of rocks along the Alpine Fault, can explain the shallow depth and localization of creeping below the Alpine Fault several hundred years after an earthquake and (3) model surface velocities are in good agreement with present-day velocities determined from campaign GPS measurements. Model strain rates show a transient response after faulting that lasts for the ﬁrst ca. 20 per cent of the interseismic period ( ca. 100 yr), with a gradual broadening in the surface velocity signal and almost constant strain rates thereafter until the next faulting event occurs. and post-seismic effects. We show that the stress transfer between


I N T RO D U C T I O N
The manner in which stresses build up and are released on and around interplate faults is crucial to estimating their importance to short-term and long-term deformation, and to seismic hazard. A particular question is how strain becomes partitioned into elastic and inelastic components during the seismic cycle. While deformation in the upper crust focuses along brittle faults, with a transfer between elastic and permanent brittle strain during seismic events, the influence of ductile localization in the lower ductile crust during the seismic cycle is still unclear. A related question is whether ductile shear zones load stress on to overlying faults during the interseismic period, or whether these shear zones merely serve to passively relieve stresses transferred from faults in the post-seismic period (e.g. Gilbert et al. 1994).
Here, we examine these questions using a comparison between a natural example, the Alpine Fault of New Zealand, and simplified numerical models incorporating static effects of fault behaviour. We show that stresses transferred to the mid-crust during an Alpine Fault type earthquake may exert a first-order effect on localization in underlying ductile crust. The results suggest that ductile localization due to overlying faults may account for a large proportion of the strain observed in exhumed mylonite zones.  Inset (b) shows plate boundary through New Zealand, with convergence rate (based on NUVEL-1A, in mm yr −1 ) between Pacific and Australian plates indicated by arrows. Relief map (a) shows location of Alpine Fault, New Zealand Southern Alps, and transects through the Alps. L1 = line 1, L2 = line 2, SIGHT transects. Box outlined by dashed lines and marked 'G' is approximate area of geodetic campaign data used to construct interseismic velocity field normal to the Alpine Fault shown in Fig. 2. Profile T is Franz Josef thermal section used in modelling by Shi et al. (1996). Panels below relief map show (c) averaged topography normal to Alpine Fault from box G (Beavan et al. 1999;Beavan et al. 2004) (V.E. = vertical exaggeration) and (d) a simplified interpretation of the collision zone based on seismicity (Leitner et al. 2001, dots with error bars), seismic reflection and refraction from the SIGHT experiment, and inversion of the geodetic data, from Sutherland et al. (2000) based on Davey et al. (1998). AF = Alpine Fault; LVZ = low velocity zone  described in text. Short grey lines indicate dipping reflectors imaged by Kleffman et al. (1998). Fault seismic movement is sparse, because no large earthquakes have occurred within historic times. However, geological estimates place the last three Alpine Fault ruptures at ca. 1717 AD, 1620 AD and 1430 AD, with the extent of rupture estimated to have been several hundred kilometres in each case, giving estimated magnitudes of M w > 7 and possibly over M w 8 (Yetton 1998;Berryman et al. 1998;Yetton et al. 1998).
Plots of earthquake locations with depth in the central South Island indicate that the seismogenic zone extends to 12 ± 2 km depth beneath much of the South Island, shallowing slightly beneath the high part of the Southern Alps, which may indicate ele-vated isotherms there (Leitner et al. 2001). The stress field derived from focal data is fairly uniform with maximum principal compressive stress directions of 110 • -120 • , approximately coincident with maximum principal shortening (or contractional strain) rates derived from GPS data. The seismicity rate of the Alpine Fault is comparable to locked sections of the San Andreas Fault (Leitner et al. 2001).
A seismic low-velocity anomaly and wide-angle reflectors have been imaged along the downdip trend of the surface fault in two geophysical transects across the central Southern Alps (Lines 1 and 2, Fig. 1; Stern et al. 2001;Scherwath et al. 2003;Van Avendonk et al. 2004). If these represent a zone of enhanced shearing and fluid release as suggested by Stern et al. (2001) they would indicate that the Alpine Fault is a major crustal structure that continues as a localized shear zone to depths of more than 30 km. This interpretation is supported by the coincidence of the low-velocity zone with a highconductivity region (interpreted to represent interconnected fluids) imaged by magnetotelluric observations beneath the central Southern Alps (Wannamaker et al. 2002). The seismic velocity data also show extremely strong lower crustal reflectors that have been interpreted to represent old oceanic crust below the greywacke .
First-order elastic strength beneath the Southern Alps can be estimated from the ratio between P-wave (Vp) and S-wave (Vs) seismic velocities. 3-D models constructed for the Southern Alps region by Eberhart-Phillips & Bannister (2002) show fairly uniform crustal Vp from 5 to 25 km depth averaging from 5.5 to 6.5 km s −1 , typical of greywacke schist. Like the observations reported above, their analysis shows the low-velocity region coinciding with the downdip extension of the Alpine Fault to at least 14 km depth. The imaged crustal root is deepest ca. 30 km east of the Main Divide and is asymmetric with a sharper gradient on the northwestern side (Eberhart-Phillips & Bannister 2002;Henrys et al. 2004).
A principal control on long-term Alpine Fault dynamics is the rapid uplift and exhumation of material along the fault and its downdip extension, which has exposed material from a depth of over 20 km within the past few million years (e.g. Cooper 1980;Little et al. 2002a;Koons et al. 2003). Studies of high-strain mylonites exhumed from these depths along the hanging wall of the Alpine Fault confirm the presence of a crustalscale shear zone inferred from the seismic data, by showing that strain is localized to depths of 20 km or more along a narrow mylonitic shear zone less than 4 km wide (Sibson et al. 1979;Little et al. 2005). Modelling has shown how positive feedback between exhumation and advection of material up the Alpine Fault can promote thermal weakening of the fault zone, which is thought to be a key control on tectonics of the Southern Alps (e.g. Koons 1987Koons , 1990Beaumont et al. 1992;Shi et al. 1996;Batt & Braun 1999;Gerbault et al. 2003). The Alpine Fault region can, therefore, be considered to be thermally weakened compared to surrounding material.

GPS data and elastic inversions
The present-day relative velocity field across the Southern Alps can be derived from Global Positioning System (GPS) repeat campaign and continuous measurements (Beavan et al. 1999(Beavan et al. , 2004Ellis et al. 2006, Fig. 1) (Fig. 2). The interseismic velocity field from these measurements has been projected onto a cross-section near SIGHT line 2 by Beavan et al. (2004) (location shown in Fig. 1). The resulting velocity profiles and derived strain-rate estimates show how strain is currently accumulating across the central Southern Alps and Alpine Fault (Fig. 2). Fault-parallel shear strain is concentrated in a narrow zone to the east of the Alpine Fault, while fault-normal velocities indicate more distributed deformation over a zone of ca. 100 km east of the fault.
Inversion of the interseismic surface velocities using elastic dislocation models suggests that part of the observed signal is due to displacement along a localized dipping shear zone in the crust beneath the Alpine Fault, with the remainder of the surface motion originating from a deeper source (Beavan et al. 1999(Beavan et al. , 2004Moore et al. 2002). The deeper source was found necessary to produce the  Beavan et al. (2004). Short grey lines indicate dipping reflectors imaged by Kleffman et al. (1998), and numbers are velocities in mm yr −1 for strike-parallel and strike-normal components of motion relative to Alpine Fault, respectively. (b) Profiles of horizontal velocity data averaged from box G (Fig. 1a), showing strike-parallel and strike-normal velocities with 1σ uncertainty estimates (grey) projected onto a profile normal to the Alpine Fault trace. The black lines are the best-fit dislocation model predictions, using the dislocation model shown in (a). (c) Profile of vertical velocity data with 1σ uncertainty estimates (grey) projected onto a profile normal to the Alpine Fault trace compared with prediction from dislocation model (black line). long-wavelength fall-off in velocities away from the Alpine Fault, and was modelled by Beavan et al. (1999) as a shear zone dipping opposite to the Alpine Fault below ca. 30 km depth, suggestive of a 'mantle subduction'-type model (e.g. Wellman 1979;Beaumont et al. 1996;Gerbault et al. 2003). However, it can equally well be fit by a zone of distributed shortening over a width of <100 km (e.g. Moore et al. 2002;Ellis et al. 2006). Ellis et al. (2006) show that variations in crustal elastic strength do not significantly affect the elastic inversions. However, because of the simplifying assumption of discrete creeping patches embedded in an elastic half-space used in those inversions, it is not possible to estimate how much of the current strain that is accumulating is recoverable and elastic (to be released in the next earthquake) versus permanent (due to inelastic processes). Here, we examine the proportion of strain accumulating in the central Southern Alps that is permanent (inelastic), and the interplay between seismic faulting and post-seismic effects. We show that the stress transfer between brittle and ductile crust occurs in a localized fashion, and that most elastic stress accumulates and is released in the upper crust.

M O D E L S E T -U P
We conduct a numerical finite element analysis using the Abaqus software (Abaqus/Standard 2004). Following Liu & Bird (2006) and because we are interested in short-term processes, we consider the present-day structure, crustal thickness, and thermal structure of the Southern Alps rather than evolving the collision zone over time. The model set-up is shown in Fig. 3. Surface topography corresponds to an average of a swath extending 25 km on either side of SIGHT line 2 (Fig. 1). The model ends at 40 km depth. Elastic moduli and density are derived from 3-D inversions of local earthquake and active source data from profile Y = −30 km in Eberhart-Phillips & Bannister (2002), and gravity is applied to the model.
The model approximates the central Southern Alps region by a 2-D cross-section in which out-of-plane motion is allowed, consistent with geophysical observations that deformation and material properties along-strike change over longer wavelengths compared to variations across-strike (e.g. Eberhart-Phillips & Bannister 2002). This approximation allows us to use a dense, structured finite element mesh (average mesh spacing = 1 km horizontally and vertically). We do not consider deformation below the crust in any detail, but instead assume that the mantle lithosphere behaves as a strong, cold frictional-elastic body, consistent with observations that mantle beneath the Southern Alps has been thickened into a cold root ca. 100 km wide (Stern et al. 2000). No adjustments due to erosion or flexure are considered in the model. Vergne et al. (2001) showed how these effects may be neglected for short-term deformation because we assume, to first order, a stationary topography.

Properties of the Alpine Fault
The Alpine Fault is represented by a frictional contact surface dipping at 45 • down to 20 km depth (Fig. 3, Table 1). The geometry roughly coincides with the present-day geometry of the Alpine Fault of New Zealand for this region (Stern et al. 2001;. In the experiments, we impose a uniform frictional coefficient μ = 0.17 along the fault from the surface to the base of the model crust, even though the actual depth extent of frictional slip is a result of the model. A coefficient of 0.17 represents a weak fault that may have undergone strain weakening and be subject to fluid overpressure, and is consistent with inferences of low fault strength from heat-flow observations along major faults (e.g. Lachenbruch & Sass 1980) and from modelling (e.g. Liu & Bird 2002).
For simplicity, we assume constant hydrostatic fluid pressure along the fault, although we acknowledge that fluid pressure variations may play an important role in the seismic cycle (e.g. Miller et al. 1999;Bosl & Nur 2000;Hobbs et al. 2004). Only static frictional behaviour is considered, so that the frictional-ductile transition along the fault in the model corresponds to the maximum rupture depth in the synoptic shear zone model of Scholz (1988).

Properties of the crust and mantle
The crust and mantle lithosphere deform by a combination of linear elasticity, Coulomb frictional plasticity, and thermally activated power-law creep ( Table 2). The transition from frictional to creep behaviour is determined naturally by the model, according to the applied geotherm, ambient pressure and strain rates, and assumed creep parameter values described below. The maximum Coulomb shear stress, τ y , along any optimally oriented plane is computed for an internal angle of friction φ = 30 • (i.e. approximately Byerlee's law conditions;Byerlee 1978;Sibson 1977) and a cohesion C = 1 MPa. The effective mean stress, p = p − p f , where p is total (mean) stress and p f is the fluid pressure (here assumed to remain at hydrostatic, as discussed above), is used in the equation: where maximum shear stress τ y is half of the maximum differential stress, (i.e. τ y = (1/2)(σ 1 − σ 3 )), and indices 1 and 3 refer to maximum and minimum principal stresses, respectively (e.g. Mandl 1988). No strain softening as a result of fault damage is applied in the model rocks surrounding the fault tip. Below the seismogenic zone depth, ductile creep will occur at rates determined by ambient temperatures, differential stresses, dominant mineral assemblages and pore fluids (e.g. Kohlstedt et al. 1995). Power-law ductile creep follows the constitutive equation: where τ v andε are the second invariants of the deviatoric stress and strain-rate tensors, respectively, Q is activation energy, n is the power-law exponent, R the universal gas constant, and T is temperature (K). (Individual deviatoric strain-rate tensor components are related to equivalent deviatoric stress tensor components using the effective viscosity derived from 2). For a 3-D model, B is related to the pre-exponential constant, A, by B = A −1 n . We choose crustal creep parameters to represent a quartz-dominated rheology, consistent with geological and geophysical interpretations that the upper and mid-crust of the Southern Alps is composed of quartzofeldspathic schist. A is determined (as are Q and n) from laboratory creep experiments on wet synthetic quartzite (Paterson & Luan 1990).

Crust-mantle interface
As mentioned above, the model only extends to a depth of 40 km, because the purpose of the experiments is to focus on the behaviour of the Alpine Fault in the crust. We define the crust-mantle interface at a transitional density of 2900 kg m −3 , where densities are derived from Eberhart-Phillips & Bannister (2002). Since according to this definition continental crustal thickness is everywhere less than 40 km a small amount of mantle lithosphere is incorporated at the base of the model. The density boundary is chosen so that our 'mantle' incorporates a dense layer at the base of the crust, which is thought to represent oceanic crust (Smith et al. 1995;Davey et al. 1998). Model crust is distinguished from mantle lithosphere by allowing thermally activated ductile creep within it; as the model only extends to a depth of 40 km, no thermally activated creep is allowed in the mantle lithosphere, only frictional and elastic behaviour. When discussing model results, we use the terms 'upper, mid, and lower' crust, where 'upper' corresponds to the area deforming predominantly by frictional mechanisms, 'lower' to that deforming by thermally activated creep, and 'mid' to the transitional region between them.

Boundary conditions
Velocity boundary conditions are applied at the left and right edges of the model, representing the oblique plate convergence velocity  (2002), with a maximum resolution of 1 × 1 km. Boundary between crust and anomalous oceanic crust/mantle is defined by density = 2900 kg m −3 (dashed line on Fig. 2b). (e) Prescribed temperature profile 1 (used in Models 1 and 2), based on thermal profile 1 from Shi et al. (1996) with no frictional shear heating. (f) Prescribed temperature profile 2 (used in Model 3), modified from Shi et al. (1996) to incorporate changes suggested by Batt & Braun (1999) and Leitner et al. (2001). from NUVEL-1A (DeMets et al. 1990Beavan et al. 1999;Beavan & Haines 2001). At the base, we consider two end-members for the style of deformation: the first end-member has distributed mantle lithosphere deformation 100-km wide (e.g. Molnar et al. 1999;Stern et al. 2000;Little et al. 2002b), whereas the second assumes a more localized shear zone 20 km wide in the mantle lithosphere beneath the downdip extension of the Alpine Fault (e.g. Wellman 1979;Beaumont et al. 1996;Batt & Braun 1999;Gerbault et al. 2003). Liu & Bird (2006) conclude that the present geophysical data are consistent with either end-member, but their model predictions slightly favour localized mantle deformation. Because we prescribe mantle kinematics at a depth of 40 km, we cannot properly address the behaviour of mantle lithosphere in our models. In particular, the assumptions used here create abnormally large deviatoric stresses within modelled parts of the mantle. However, while the choice of boundary condition for the lower boundary has an effect on the broad-wavelength signal seen at the surface (e.g. Ellis et al. 2006), model sensitivity studies have shown that it does not affect mechanisms of stress transfer along the crustal fault to any great extent.

Thermal structure
A key requirement to model the present-day structure of the Southern Alps is to impose a realistic temperature profile with depth. Heat-flow measurements are few, and are likely to be influenced by near-surface fluid circulation (Shi et al. 1996). Heat-flow studies, fission-track analyses, fluid inclusion studies, and numerical modelling all show the importance of high-erosional exhumation rates in producing a heat-flow anomaly near the Alpine Fault (Allis & Shi 1995;Koons 1987;Shi et al. 1996;Tippett & Kamp 1993;Beaumont et al. 1994Beaumont et al. , 1996Batt & Braun 1999;Gerbault et al. 2003). There is less agreement between different models about the shape and magnitude of the thermal anomaly. For example, Shi et al. (1996) and Gerbault et al. (2003) show a maximum upward deflection of isotherms to the east of the surface trace of the Alpine Fault, matched by a downward perturbation to the west and at greater depths, whereas the models of Batt & Braun (1999) and Koons (1987) show a larger upward perturbation of isotherms focused on the Alpine Fault itself. However, Leitner et al. (2001) in a study of the maximum seismogenic depth for three transects across the Southern Alps, found no evidence of a downward displacement in the seismogenic zone as would be expected if isotherms dipped to the west of the Alpine Fault as predicted by Shi et al. (1996). The seismogenic zone below most of the South Island extends to ca. 12 km depth (Leitner et al. 2001). If the seismogenic depth approximately indicates the depth to the 350 • C isotherm (e.g. Scholz 1990), this suggests higher near-surface geothermal gradients than those used by Shi et al. (1996), more like the values used by Koons (1987). The analysis of Leitner et al. (2001) also suggests shallowing of the lower boundary of the seismogenic zone by ca. 4 km a few km east of the Alpine Fault in the central part of the Southern Alps, consistent with a higher geotherm due to advective heat transport of the rapidly exhuming hanging wall block.
To account for the uncertainties in thermal models, we take two 'end-member' configurations. In Models 1 and 2, we impose a thermal structure (temperature profile 1) based on Model 1 from the Franz Josef thermal profile of Shi et al. (1996) without frictional heating (Fig. 3e). In Model 3, we use a thermal profile (temperature profile 2) that is based on the observations from Leitner et al. (2001) where the 350 • C isotherm is assumed to occur at a depth of 12 km,  Table 2. Material properties used in models 1-3; other properties described in Fig. 3.

S TA RT -U P P H A S E
A start-up phase is used to attain gravitational stability and to simulate the long-term build-up of differential stress in the crust (Fig. 4).
The three panels (a-c) in Fig. 4 show results for Models 1-3. Model 1 uses temperature profile 1 and assumes that the mantle lithosphere deforms over a broad region beneath the collision zone; Model 2 uses temperature profile 1 but prescribes the mantle lithosphere to deform over a shorter length-scale; and Model 3 is similar to Model 2 but uses temperature profile 2. Transpressive boundary conditions corresponding approximately to those derived from NUVEL-1A are applied at the right hand boundary at rates of 10 mm yr −1 (normal to strike) and 35 mm yr −1 (parallel to strike). The start-up phase is run in two parts: first, the boundary conditions are applied over a time period of 200 ka (i.e. giving 2 km of convergence normal to the contact surface). The deformation from this phase is then discarded, while the stress field is used as an initial condition in an 'undeformed' model and the boundary conditions then applied for a time period of 50 ka (i.e. 500 m convergence normal to the contact surface). This twostage set-up allows stresses to reach frictional yield in the upper crust and mantle, and steady ductile flow in the lower crust, without running the model to large amounts of deformation.
During the start-up phase, the contact surface (which may be regarded as a 'fault') has a uniform coefficient of friction of 0.17 prescribed along it. The contact algorithm in Abaqus determines slip as a function of depth and time-that is, slip along the fault is solved for and is not prescribed. In all three models shown in Fig. 4, the fault slips steadily along its upper portion (the top 10-15 km) once frictional yield along the contact surface is exceeded. Frictional slip tapers off below 10-15 km because the mid-and lower crust are ductile and weaker than the frictional fault due to effects of increasing pressure and temperature. As a result, at mid-crustal depths and below, ductile flow around (but not on) the contact surface takes over from the frictional motion along the fault, as described in similar models by Cattin & Avouac (2000), Chéry et al. (2001), and Rolandone & Jaupart (2002). Note that while the model includes a contact surface throughout the crust, it is only an active feature for frictional displacement, and the location of the contact surface below this displacement (i.e. below 10-15 km) has no specific effect on the creep. Thus the thermal and frictional properties of each model determine to what depth the contact surface becomes the fault. In the following discussion 'fault' means the surface that has experienced frictional displacement.
The stress field at the end of the start-up phase in Models 1-3 shows the following characteristics (Fig. 4): (1) a differential stress maximum in the mid-crust (ca. 8-10 km) with stresses decreasing in the lower crust due to ductile flow, as indicated in the plots to the right of each panel in Fig. 4; (2) a decrease in the differential stress within a few kilometres of the up-thrown side of the fault, resulting from frictional yielding along it and (3) high stresses in the mantle lithosphere beneath the fault, resulting from the basal boundary condition applied there, which are not considered realistic (as discussed in the previous section). Apart from the different stress distributions in the upper mantle, the resulting stress fields for Models 1 and 2 are very similar. Model 3 has lower stresses in the crust due to the higher temperatures of temperature profile 2. As a result, the high-strength layer and frictional yielding in the crust are restricted to the top 10 km. In all three models, frictional slip along the fault takes up a high proportion of the prescribed far-field velocities at depths where the fault is actively slipping (ca. top 10 km of the fault plane, Fig. 4d). This is especially the case for Model 3, where 65 per cent of strike-normal and 80 per cent of strike-parallel plate boundary motion is accommodated along the fault, although frictional slip ends at a shallower depth than Models 1 and 2. The higher amount of frictional slip in Model 3 cf. the other models is caused by the additional thermal weakness beneath the fault resulting from use of temperature profile 2, as discussed in Section 5.9.
All three models show differential stresses >100 MPa building up around the lower termination of the actively slipping fault (downthrown side of fault, ca. 8-11 km depth, Figs 4a-c). This stress build-up results from the termination of slip over this depth interval, where frictional faulting gives way to thermally activated ductile creep (White 1996;Ellis & Stöckhert 2004a,b). The slipping fault loads the surrounding mid-crust elastically, by an amount that depends on the elastic strength of the medium and the depth range over which slip tapers (ca. 10-15 km, Models 1 and 2; and ca. 8-10 km, Model 3). The stress build-up is limited by the ability of the midcrust to relax differential stress by ductile flow, so that by the end of the start-up phase, a steady state is reached. Directions of principal compressive stress (not shown) are consistent with observations from seismic focal mechanisms (Leitner et al. 2001). However, principal stress directions are rotated within a few km of the fault plane as described in Ellis & Stöckhert (2004a).
Despite the differences in slip and stress patterns, the overall form of Models 1-3 is very similar. Although we continue to show results from all three models in the rest of the paper, we particularly focus on Model 3, because: (1) the depth at which frictional slip terminates in Model 3 most closely matches that predicted from inversion of GPS data in the Southern Alps; (2) the temperature profile better matches that deduced from the observations of Leitner et al. (2001) and (3) the strong, frictional part of the crust in Model 3 is slightly thinner than the other two models, with no significant downward deflection in the maximum depth of the frictional zone to the east of the fault, consistent with the observations from Leitner et al. (2001).

R E S U LT S F O R T H E S E I S M I C C YC L E
To simulate the effects of a seismic cycle, the fault is locked for a specified period of 500 yr after the start-up phase described above. Then, following this 'interseismic' interval, the fault friction coefficient is reduced back down to its initial value (μ = 0.17) for a short period (1 day), simulating the static response to the faulting step of the seismic cycle. No dynamic or rate-state effects are modelled. Note that to ensure that we have reached a dynamic steady-state where each earthquake cycle produces identical deformation, we run the models through a number of complete seismic cycles until a cyclically constant pattern emerges, and then analyse one cycle in detail.
To keep things simple, we use the same frictional coefficient during the seismic phase as we did during the set-up phase. It may be argued that a slightly lower frictional strength for the fault should be used during the seismic interval, so that the longterm average of the fault strength is the same as for the set-up phase. However, the stress drop over the upper part of the fault is a small fraction of the total stress along the fault, so the difference in behaviour using a slightly lower friction coefficient is not significant. Fig. 5 shows the change in differential stress in the models when the contact surface friction drops suddenly, simulating the static effect of an earthquake. As expected, the upper portion of the fault shows a small stress drop of a few MPa. This is matched by a larger increase in stress at the lower end of the slipping part of the fault (depths ca. 7-10 km; Figs 5b, d and f, solid line). The frictional slip along the fault of several metres (Figs 5c, e and g) transfers stress from the upper frictional crust to the mid-crust. Stress changes result from elastic loading and unloading of the crust near the fault plane, and their magnitude depends on the bulk elastic strength and the interplay between elastic and viscous rheologies in the mid-crust. The stress changes for seismic and interseismic steps are almost exact mirror images of each other, indicating that the models have attained a kind of cyclical 'steady state' where short-term stress changes during the seismic cycle oscillate around the stress fields shown in Fig. 4. The stress drops during faulting in the top 10 km (Models 1, 2) or 7 km (Model 3) are consistent with small stress drops inferred from seismic observations (e.g. Bouchon 1997;McGarr & Fletcher 2002), but the stress increases at ca. 15 km (Models 1, 2) or 8-10 km (Model 3) depth show that the stress behaviour is different in the semi-brittle mid-crust (e.g. Ellis & Stöckhert 2004a,b;Rolandone et al. 2004).
The elevated stresses at the lower termination of the fault, enhanced by the sudden slip every 500 yr, causes localized ductile flow in the mid-and lower crust beneath the fault (e.g. for Model 3; Fig. 6a). Strain rates attain their maximum just after the seismic step is completed, when mid-crustal differential stresses along the fault plane are highest. The elevated strain rates are directly caused by termination of fault slip, and there is no component of dynamic ductile localization (as defined by Montesi & Zuber 2002) since the model contains no dynamic localization mechanism. Instead, focused ductile shearing occurs due to the imposed (as defined by Montesi & Zuber 2002) and localized stress transfer from the fault. After the entire interseismic part of the cycle, the accumulated creep strain (Figs 6b-d) beneath the fault is highly localized in a lobe dipping beneath the fault and down to the base of the crust, and ca. 20-100 times higher than background rates.
The focusing of strain beneath the fault in the models has an effect on the displacements accumulating at the surface during the seismic cycle (e.g. Model 3; Figs 6e-g). Enhanced strain rates beneath the fault during the interseismic period cause a surface displacement pattern similar to that observed accumulating today from GPS measurements (Fig. 7). In particular, localized creep in the mid-lower crust causes a short-wavelength fall off in both strikenormal and strike-parallel components above the (locked) fault. By adding this signal to that produced during the seismic interval, it is possible to predict the total displacement due to one entire seismic cycle (Figs 6e-g bold lines). These show, as expected, a discontinuity at the surface trace of the fault, but also distributed deformation (contributed by the interseismic part of the signal). For the particular thermal profile and boundary conditions used here, Model 3 predicts that about 50 per cent of normal and 80 per cent of strike-slip plate motion in the crust will be accommodated on the fault, with the rest distributed on either side over a zone ca. 100 km wide.

D I S C U S S I O N
Any comparison between the model results and the true behaviour of the Alpine Fault of New Zealand must necessarily be at a simple level, because of the assumptions and approximations used in constructing the models. Chief among these are: (1) the assumption that out-of-plane deformation is uniform along strike, when differences in long-term exhumation rates  along-strike in the central Southern Alps have been documented (e.g. Little et al. 2002aLittle et al. , 2005; (2) dynamic strain localization has been ignored, yet might be expected-for example, due to weakening of the material by progressive development of mylonitic fabric and fluid release in the downdip extension of the Alpine Fault, which would enhance localization at mid-and lower-crustal depths; (3) the temperature profiles are not well constrained; (4) the resolution of the model, limited to 1 km, will tend to broaden any localization effects; (5) heterogeneities in elastic, frictional and ductile strength on scales less than 1 km are under-represented due to limitations in material data set and model resolution; (6) the models allow no true evolution in time of topography or rheology; (7) dynamic stress effects are neglected during the seismic step; (8) the imposed seismic cycle is regular in time and space and (9) the presence of secondary weak faults within the brittle crust (representing, for example, faults along the Main Divide of the Southern Alps; Cox & Findlay 1995) is not considered.
Bearing these shortcomings in mind, the models are nevertheless expected to provide some insight into the behaviour of real plate boundary faults such as the Alpine Fault of New Zealand, as discussed below.

Boundary effects and comparison of model results to GPS observations in the New Zealand Southern Alps
Previous work has shown how the surface interseismic deformation can be fit adequately with an elastic dislocation model that has far fewer parameters than we use in the models here (e.g. Fig. 2). The aim here is to show that a more physically consistent model can also adequately fit the data, as well as providing additional insight into the interaction between elastic and inelastic processes. Models 1-3 are designed to investigate near-fault effects only -the limited horizontal and vertical dimensions of the models precludes investigation of broader-scale features associated with mantle deforming over wavelengths of greater than about 100 km. For these reasons, nothing would be gained by a detailed, statistical comparison between the models investigated and observations of GPS deformation. However, qualitatively, the horizontal velocities predicted from the models for a late stage in the interseismic step are similar to those seen in the contemporary GPS signal above the Alpine Fault (Fig. 7). The elevated creeping lobe beneath the 'locked' fault seen in Models 1-3 (Figs 6b-d)   et al. (1999,2004) (Fig. 2b). Based on comparison between the models and the geodetic results, we suggest that the creeping patch in the geodetic inversions represents a mylonitic shear zone, extending downdip along the Alpine Fault to the base of the quartz-dominated crust and deforming in a narrow, focused lobe (e.g. . Although this mylonite zone is likely to have experienced dynamic localization with increasing strain, the model results shown here suggest that localization in a narrow shear zone can be achieved simply by a combination of thermal weakening due to exhumation along the fault plane, and stress transfer between the upper and midcrust along the Alpine Fault (i.e. 'imposed' localization, Montesi & Zuber 2002). This stress transfer occurs seismically and is due to the fall-off of frictional slip at mid-crustal levels, which loads the midcrust elastically. The extra stress is then relieved by post-seismic localized creeping in the mid-lower crust extending downwards beneath the lower end of the seismically slipping fault. This localized creeping is aided by thermal weakening of ductile crust with increasing depth. Rapid exhumation of the hanging wall block along the Alpine Fault, leading to an upward deflection in isotherms (and thermal weakening) near it, has also enhanced this localization.

Comparison between Models 1, 2 and 3
All three models that we investigated show similar qualitative behaviour, as discussed in Section 5.1. The main qualitative difference is the deeper slip distribution and smaller magnitude fault stress cycling in Models 1 and 2 versus Model 3, owing to the difference in geothermal gradient (Fig. 5). The lower gradient used in Models 1 and 2 deepens the transition zone at which frictional slip along the fault plane gives way to thermally activated ductile creep compared to the hotter case considered in Model 3. Correspondingly, frictional slip falls off over a larger depth interval in Models 1 and 2, causing a smaller stress increase during the seismic period (Figs 5b, d cf. f). It might be expected that the 'stronger' ductile behaviour in Models 1 and 2 would cause greater amounts of plate boundary deformation to be taken up along the fault plane, giving rise to larger amounts of frictional slip. However, the temperature profile used in Model 3 (Fig. 3f) also has elevated temperatures beneath the fault at mid-crustal depths. The effect is to further localize creep there, helping to load the overlying fault. This point is discussed further in Section 5.9.
A comparison between Models 1 and 2-3 demonstrates that the behaviour of the lower lithosphere has little influence on fault slip in the upper crust. The predicted vertical velocity for Model 1 is slightly lower than for Models 2 and 3 (Fig. 7). Model 2 has slightly more frictional slip (Fig. 5 e vs. c) which suggests that localized mantle deformation is enhancing ductile creep in the lower crust and causing further loading of the fault plane, but the effect is rather small.

'Decoupling' versus 'detachment' between crust and mantle lithosphere
The steady-state behaviour illustrated during the model set-up (Fig. 4) is very similar for Models 1, 2 and 3. In all cases, the weak fault in combination with thermal weakening caused by erosional exhumation along the fault controls crustal deformation. All three models (Fig. 4) show a low strength in the lower crust, depending on the assumed rheology and geotherm. A weak, ductile rheology in the lower crust is supported by the cut-off in seismicity observed at mid-crustal depths beneath the Southern Alps (Leitner et al. 2001).
The weak strength of the base of the crust and the observation that Models 1 and 2 have a similar behaviour despite different prescribed deformation in underlying mantle lithosphere, suggests that the crust is decoupled from the mantle there. However, this does not necessarily indicate that wholescale detachment between crust and mantle lithosphere is occurring. Despite the low strength, the plot  of creep strain (Figs 6b-d) indicates little relative shearing between them, except in the region immediately above the deforming mantle lithosphere, which acts to transfer localized crustal strain to the more distributed mantle deformation prescribed by the boundary conditions. The difference between 'decoupling' (the inability of one layer to transmit stress to another) and 'detachment' (a strong shear zone developing between layers) was investigated by . They found that crustal localization can prevent detachment from lower layers, even though little stress is transmitted across the crustmantle interface. If both crust and mantle lithosphere experience frictional and viscous weakening early on in the collision, before substantial thickening occurs at the plate boundary, they will continue to deform at a coincident location.

Effect of rheology on frictional slip termination distance
The depth at which frictional slip ceases along the model contact surface depends on its frictional strength, and the rheology of the surrounding medium. The steeper the geothermal gradient, the shallower the transition from frictional slip to ductile creep in surrounding regions. The depth range over which frictional slip decays is also narrower for a steeper geothermal gradient, causing higher elastic stress loading after faulting (Fig. 5f). The difference in slip decay can be seen by comparing Models 1 and 2 versus Model 3 (Figs 5c, e and f), which has a steeper geothermal gradient near the fault (Fig. 3). Note that the depth at which frictional slip ceases in the models is likely to be greater than the true seismogenic depth, as no distinction is made in the models between regions of unstable and conditionally stable slip (cf. Scholz 1990).
In cases where fault slip tapers more suddenly than that shown here, for example due to asperities or changes in fault geometry, even larger stresses may build up at the fault termination (e.g. Ellis & Stöckhert 2004a). Even for the comparatively 'smooth' fault case examined in Models 1-3, the seismically induced perturbation of the stress field beneath the fault tip, with stress cycling of up to ca. ±15 MPa, causes significant localized deformation along the extension of the fault in the ductile lower crust.

Where does elastic strain accumulate?
A critical question for hazard estimation along major faults such as the New Zealand Alpine Fault, is how much of the strain accumulating during an interseismic period will be released during the next major earthquake. The models can be used to investigate this question by separately plotting the interseismic strain that is accumulating elastically versus that accumulating inelastically (e.g. Fig. 8).
Most of the elastic strain accumulates in the brittle crust in a narrow region around and above the fault, above the transition from frictional slip along the fault plane, to aseismic creep below. By comparison, inelastic straining (predominantly creep) is focused around the downdip continuation of the fault, in a narrow lobe extending down to the base of the ductile crust.
This result shows that most lower-crustal deformation is permanent and localized, whereas upper crust deformation is transient and elastic, explaining why elastic half-space inversions (which assume a localized slipping patch in the mid and lower crust) work quite well at modelling interseismic deformation (e.g. Beavan et al. 1999Beavan et al. , 2004. These simpler models can, therefore, provide quite good estimates of how much elastic strain is accumulating within the crust. We caution though, that the models here do not consider deformation integrated over thousands of years, during which time permanent frictional strain in the upper crust (insignificant in Fig. 8) will accumulate and may eventually localize onto a system of secondary faults.

Predictions of permanent displacement at the surface
Plots of surface displacement during each component of the seismic cycle (Fig. 6e-g) can be used to estimate how much of the displacement currently accumulating there will be turned into permanent displacement over multiple seismic cycles. The sum of the interseismic and seismic components shows that about half to threequarters (50-75 per cent) of the permanent surface displacement in Model 3 is taken up along the fault plane, with the rest distributed, mostly to the right of the fault, over a wavelength of about 50-100 km. The proportion of slip occurring along the fault in Model 3 is in rough agreement with GPS and geological estimates that 60-75 per cent of total plate motion is taken up along the Alpine Fault (Norris & Cooper 2001;Beavan et al. 1999).

Change in interseismic signal with time after last major earthquake
Transient post-seismic surface displacements following major earthquakes have been documented for many tectonic settings (e.g. Denali 2002 (Pollitz 2006); Landers 1992; (Pollitz et al. 2000); Loma Prieta 1989 (Pollitz et al. 1998;Segall et al. 2000)). Postseismic effects have variously been attributed to afterslip along the lower part of the fault plane (e.g. Shen et al. 1994;Savage & Svarc 1997;Owen et al. 2002;Hudnut et al. 2002); viscoelastic relaxation in the ductile lower crust or mantle (Deng et al. 1998;Pollitz et al. 2000Pollitz et al. , 2001Freed & Bürgmann 2004); or poroelastic rebound (Peltzer et al. 1998;Fialko 2004). It is, therefore, reasonable to ask whether the models investigated here predict a transient post-seismic response that would be detectable by GPS geodesy, bearing in mind the following: (1) Frictional afterslip in the models is included in the seismic interval, since we do not distinguish between unsteady slip over seconds, and steady slip over a longer time periods of days to weeks following an earthquake; (2) Viscoelastic relaxation in the mantle is not modelled, as we only consider the upper 40 km of the system; (3) Poroelastic effects are not modelled. Any change in the predicted model surface velocities over the interseismic period will, therefore, only indicate effects of viscoelastic relaxation in the lower crust.
Model 3, shows subtle, but consistent, changes in the interseismic velocity field at the surface as a function of time (Fig. 9). This is consistent with the enhanced strain rates below the locked fault decaying post-seismically over time as stresses are relieved there. The decay time depends on the relaxation of stresses in the mid-crust, which in turn depends on the assumed geotherm and rheology. The velocity signal gradually broadens in time as less deformation is taken up along the downdip extension of the fault.
The small magnitude by which the surface velocities change over the interseismic period in the models is much less than has been observed after large earthquakes in other tectonic settings, particularly the western US, where large post-seismic effects have been attributed to a weak and low-viscosity upper mantle (e.g. Freed & Bürgmann 2004). Some of the discrepancy in magnitude may be caused by the various model approximations listed above, as well as the fact that the models do not distinguish between seismic slip and frictional afterslip on the fault plane. In addition, the mantle lithosphere beneath the Southern Alps has been thickened and is relatively cool (Stern et al. 2000), whereas the mantle lithosphere beneath the western US is thought to be thin and warm, and is likely to have a much lower viscosity (e.g. Dixon et al. 2004). The one earthquake in the central South Island for which we have measurements of coseismic and post-seismic deformation, the 1994 Arthur's Pass earthquake, did not experience any measurable post-seismic decay (Ellis & Beavan 2001). For these reasons we are not convinced that a large post-seismic signal caused by mantle relaxation would follow a major Alpine Fault earthquake. However, owing to the approximations used in the models, we cannot rule out either mantle relaxation or significant afterslip in addition to the effects seen in Fig. 9.
It is conceivable that the changes over the interseismic interval seen in Model 3 (e.g. maximum change of ca. 10 per cent, Fig. 9) could be measured geodetically, if the signal were continuously ob-strike-parallel velocity (mm yr served in the first 100 yr or so after a major earthquake. Later into the interseismic period, the rate of change in the surface velocity field slows down and it is doubtful that it could be measured. In the Southern Alps, the last large Alpine Fault earthquake is thought to have occurred several hundred years ago (e.g. Wells et al. 1999).
On the basis of comparison with Model 3, we, therefore, expect that the interseismic velocity field is now steady in time, although it may be influenced by smaller seismic events occurring on secondary faults. In Model 3, the near-fault signal late in the interseismic period is about 10 per cent less than its averaged value for the entire interseismic period.

Controls on initiation and style of crustal localization
The models explored here illustrate that frictional weakening along major faults can enhance deformation in both upper (frictional) and lower (ductile) continental crust. This view is supported by observations from the Southern Alps, where evidence of anomalous resistivities, seismic reflectivity and inferences from rocks exhumed at the surface all support the theory that localized ductile flow occurs directly downdip from the surface trace of the Alpine Fault (Wannamaker et al. 2002;Stern et al. 2001;Little et al. 2002a;Norris 2004). We have shown that imposed localization (Montesi & Zuber 2002) in the lower crust, together with thermal weakening related to rapid exhumation of the hanging wall block along a major reverse fault, is sufficient to explain many of these observations. There are two alternative views on the causes of crustal localization. Firstly, localization may arise spontaneously within ductile lower crust (e.g. Montesi & Zuber 2002;Regenauer-Lieb et al. 2004). This lower-crustal localization and associated ductile weakening may then trigger development of a weak frictional fault above it. A second view is that underthrusting of mantle lithosphere along a shear zone (e.g. the 'antifault' proposed by Wellman 1979) triggers the crustal shear zone above it (e.g. Beaumont et al. 1996;Little et al. 2002a).
The question of whether ductile or frictional localization arises first is most likely a 'chicken or egg' scenario. A combination of the two may lead to rapid localization of crustal deformation at plate boundaries. However, many plate boundaries, including the Alpine Fault, have evolved from earlier configurations (Sutherland et al. 2000), where a pre-existing weakness probably initiates localized deformation (e.g. Gerbault et al. 2003). Rapid exhumation along faults may then further contribute to localization of deformation at the plate boundary (e.g. Koons 1987;Beaumont et al. 1992;Batt & Braun 1999).
The stress concentration at the lower termination of the frictional part of the fault and the non-linear ductile rheology act to focus ductile creep in a narrow dipping lobe beneath it. The lobe is narrowest directly below the fault termination, in agreement with observations from some exhumed fault zones (e.g. Sibson 1977;West & Hubbard 1997;Hanmer 1988) but its width does not increase significantly even near the base of the crust (Fig. 6), because of the decrease in viscosity with increasing temperature. Model 3 indicates a width of about 3-5 km throughout the lower crust. This is comparable, though slightly larger, than suggested by  based on geologic observations, who estimate that most deformation on the downdip extension of the Alpine Fault takes place along a 1-2 km wide mylonitic shear zone through the middle and lower crusts. The discrepancy may be partly due to the limited resolution of the grid in Model 3, but also can be explained in terms of dynamic localization (not modelled here) due to effects such as grain-size reduction and fluid release. Additionally,  state that the 1-2 km width derived in their study is a minimum estimate, and that the mylonite zone may have been truncated during uplift. A more recent study by Little et al. (2005) suggests that the shear zone may be up to 4 km wide.

Imposed weakening in lower crust: does the fault load the lower crust, or vice versa?
Elevated strain rates beneath the frictionally slipping part of the fault in Models 1-3 can cause strain-rate weakening (a lowering in the effective viscosity) in the surrounding region. Over time, the uplift and exhumation of mid-lower crustal rocks along the fault has also caused an elevation in temperatures near it ( Fig. 3; Shi et al. 1996), which further reduces the viscosity near the downdip extension of the fault. Therefore, although the models include no spontaneous weakening of the lower crust, an 'imposed weakening' occurs beneath the fault. Strain is concentrated there, and the system reaches an approximate steady state where stress changes in interseismic and seismic periods roughly balance. Motion along the fault clearly 'stress loads' the mid-crust and induces elevated creeping there (Figs 5b, d and f). We may then ask: Does ductile straining beneath the fault contribute to 'stress loading' of the fault during the interseismic period?
To determine this, an additional 'interseismic' step was run for Model 3, with ductile creep suppressed. This allowed us to estimate stress loading along the fault due to creeping below it (i.e. the nominal case), by comparison with a case with no creep. For the case with creeping, an average stress increase of ca. 1.5 MPa occurs during the interseismic period along the upper 10 km of the fault, giving an average stressing rate of 3 kPa yr −1 . When ductile creep is suppressed, the fault stress increases at an average rate of 1.2 kPa yr −1 . This indicates that in Model 3, the fault loading is enhanced by a factor of ca. 2.5× due to elevated creeping below it. The effect of thermal weakening on stress loading along the fault can also be seen by comparing frictional slip magnitudes in Model 3 versus Models 1 and 2 (Figs 5c, e and f). The higher slips in Model 3 correspond to the greater perturbation in the geotherm used in that model. If ongoing creep caused ductile strain weakening beneath the fault (e.g. due to development of a mylonitic fabric, and/or release of fluids), fault loading from below would be further enhanced. Even without these additional factors, the models show that a weak fault and thermally weakened shear zone in the lower crust can act to reinforce and load each other over the seismic cycle.

C O N C L U S I O N S
Using numerical models, we have simulated first-order effects of the New Zealand Alpine Fault seismic cycle on deformation in the surrounding crust. In particular, we have shown that: (1) Stress transfer from upper crust to mid-crust due to seismic faulting can promote localization beneath a major, frictionally weak fault such as the Alpine Fault of New Zealand.
(2) This stress transfer may help to explain localized deformation beneath the Alpine Fault, as predicted by geophysical and geological studies (3) A comparison between model surface velocities throughout the seismic cycle is in reasonable agreement with present-day horizontal velocities measured by GPS, and with slip rates estimated along the Alpine Fault from geological studies. The model can be used to visualize the different parts of an entire seismic cycle, and thus help to tie observations at different timescales together.
(4) Predicted surface velocities initially fall off rapidly, and then more gradually, over the interseismic period. Model predictions, if applicable to the Alpine Fault seismic cycle, suggest that for times greater than ca. 100 yr after a large earthquake, it will not be possible to resolve exactly which part of the cycle we are in from surface velocities.
Further tests are needed to clarify additional effects due to ductile strain weakening in the lower crust, fluids, and mantle lithosphere dynamics. However, the results demonstrate that even without dynamic ductile weakening, the low frictional strength of a fault such as the Alpine Fault of New Zealand, coupled with thermal weakening due to rapid exhumation up the fault ramp, can account for the near-fault interseismic velocity field presently observed at the surface.

A C K N O W L E D G M E N T S
This study was funded by the Foundation for Research, Science, and Technology under the GNS Science programme 'Geological Hazards and Society'. Some of the figures were created using the GMT public domain software. Discussions with Laura Wallace, Steve Cohen, Rupert Sutherland, Jean Chéry, Martha Savage, Tim Little, Richard Norris, Martin Reyners, Joan Gomberg, Paul Bodin, and many other colleagues helped improve the ideas behind the manuscript. Careful and constructive reviews by Roland Bürgmann and an anonymous reviewer, as well as informal reviews by Fred Davey, Stuart Henrys, and Rupert Sutherland are appreciated. The updated vertical velocity data were reproduced with the permission of Peter Molnar, Brad Hager, Tom Herring and Paul Denys.