Compression of the entire continental lithosphere is considered using twodimensional numerical models, in order to study the influence of the lithospheric mantle on the geometry of continental collision in its initial stages. The models are based on the central section of New Zealand Southern Alps, where continental collision has occurred along the Alpine Fault since about 7 Ma. They incorporate brittleelastic ductile rheology, heat transfer, surface processes, and fault localisation. The models are compared to the surface relief, the GPS convergence velocity, the measured electrical conductivity, and the geometry of the crustal root imaged from seismic velocity measurements.

The crustal deformation is characterized by localized uplift at the plate boundary (Alpine Fault) and by two secondary zones of faulting. One is located about 60 to 80 km east of the Alpine Fault, at the start of upper crust bending (or tilting), and the other is located about 100 to 130 km east of the Alpine Fault as a result from shear strain propagating to the surface through the ductile lower crust.

The observed asymmetric shape of the crustal root is best reproduced for a mantle lithosphere strength of the order of 200 MPa, and an intermediate rate of strainsoftening. A lower strength of the mantle lithosphere can produce symmetric thickening, but provides an amplitude of the crustal root too small when compared to observations.

The observed 20 km offset between the maximum in surface relief and the crustal root was not satisfactorily reproduced. This offset is most likely due to the three dimensionality of oblique collision in New Zealand Southern Alps.

INTRODUCTION

The Australian/Pacific plate boundary through the South Island of New Zealand has progressively formed since about 45 Ma ago, transforming an initial Eocene passive margin into the dextral strikeslip Alpine fault, with about 450 km of offset [START_REF] Wellman | The geology of the west coast from Abut Head to Milford sound, Part 1[END_REF][START_REF] Carter | Cenozoic history of southern New Zealand: an accord between geological observations and platetectonic predictions[END_REF][START_REF] Sutherland | Plate boundary deformation in South Island, New Zealand, is related to inherited lithospheric structure[END_REF]. Oblique convergence commenced about 7 Ma ago, and has resulted in the uplift of the Southern Alps [START_REF] Wellman | An uplift map for the South Island of New Zealand, and a model for uplift of the Southern Alps[END_REF], De Mets et al., 1990[START_REF] Beavan | Contemporary horizontal velocity and strainrate fields of the PacificAustralian plate boundary zone through New Zealand[END_REF], Figure 1).

In regions of oblique relative plate motion and where both shear and compressive deformation occur together, twodimensional (2D) models (for example, [START_REF] Willett | Orography and orogeny: The effects of erosion on the structure of mountain belts[END_REF][START_REF] Batt | The tectonic evolution of the Southern Alps, New Zealand: insights from fully thermally coupled dynamical modelling[END_REF]) can help understand which deformation structures are due solely to the compressive component of the system. [START_REF] Koons | Twosided orogen; collision and erosion from the sandbox to the Southern Alps[END_REF] and [START_REF] Beaumont | The continental collision zone South Island, New Zealand: comparison of geodynamical models and observations[END_REF] first applied such models to the Southern Alps of New Zealand. They identified that the parameters that control the largescale geometry of the orogen are a combination of crustal rheological properties, surface processes, and boundary conditions applied at the base of the crust. However, these models assume that the Pacific crust is driven by lithospheric mantle against the Australian lithosphere, so that an imposed basal velocity discontinuity (commonly referred to as the 'S point') controls the location of deformation. Modelling of the interaction between the crust and the mantle lithosphere is biased by this prescribed discontinuous boundary condition. During the mountain building process, deflection of the Moho should also be constrained by the ability of the underlying mantle lithosphere to deform. In this paper we incorporate the mantle lithosphere when examining the growth of the Southern Alps.

The lithospheric mantle can undergo tectonic shortening in three principal ways: homogeneous (or volumetric) shortening by distributed thickening, folding (or buckling) by subperiodical bending of competent lithospheric layers, or intracontinental subduction by asymmetric underthrusting of one block under the other along a major shear zone. While the occurrence of symmetric or asymmetric mantle thickening is an old debate [START_REF] Houseman | Convective instability of a thickened boundary layer and its relevance for the thermal evolution of continental convergence belts[END_REF]Mattauer, 1986), it is only more recently that largescale folding is understood to be a primary response of the lithosphere to induced compressional stress (e.g. [START_REF] Martinod | Periodic instabilities during compression or extension of the lithosphere 1. Deformation modes from an analytical perturbation method[END_REF][START_REF] Burov | Continental lithosphere folding in Central Asia (II): constraints from gravity and topography[END_REF][START_REF] Gerbault | Do faults trigger folding in the lithosphere?[END_REF][START_REF] Cloetingh | Lithosphere folding; primary response to compression? (from Central Asia to Paris Basin)[END_REF]. Lithospheric folding develops due to its strength contrast with low viscosity asthenosphere, and evolves as a single megafold where rheological weaknesses concentrate [START_REF] Burg | Lithospheric scale folding; numerical modelling and application to the Himalayan syntaxes[END_REF][START_REF] Cloetingh | Lithosphere folding; primary response to compression? (from Central Asia to Paris Basin)[END_REF].

At the timescale of 10 to 20 Ma, symmetric downwelling of the mantle lithosphere may also result from a convective or a RayleighTaylor instability, due to cold and dense thickening mantle lithosphere (e.g. [START_REF] Houseman | Convective instability of a thickened boundary layer and its relevance for the thermal evolution of continental convergence belts[END_REF]. This mechanism has been proposed for the Southern Alps of New Zealand, in order to explain recorded Pwave delays [START_REF] Molnar | Continuous deformation versus faulting through the continental lithosphere of New Zealand[END_REF][START_REF] Stern | Teleseismic P wave delays and modes of shortening the mantle lithosphere beneath South Island, New Zealand[END_REF]. According to models of large scale mantle compression from [START_REF] Pysklywec | Modeling the behaviour of the continental mantle lithosphere during plate convergence[END_REF], a RayleighTaylor type instability or asymmetric intracontinental subduction can evolve into dripping or slab breakoff of the mantle lithosphere, depending on the rate of convergence, the density distribution and the strength of the mantle lithosphere.

In contrast to the approach of [START_REF] Pysklywec | Modeling the behaviour of the continental mantle lithosphere during plate convergence[END_REF], we investigate the interaction between crust and mantle lithosphere at the New Zealand collisional plate boundary, by considering only positive density gradient with depth and incorporating preexisting heterogeneities based on known rheological discontinuities (terranes of different tectonic history, Figure 1a). For example the Alpine Fault is a major strikeslip feature that has been in existence for about 45 Ma prior to the actual compressional episode [START_REF] Sutherland | Plate boundary deformation in South Island, New Zealand, is related to inherited lithospheric structure[END_REF], and we model it as an initial vertical zone of weakness. Since in addition, we ask how the strength contrasts determine the asymmetry of the lithospheric mantle and its effects on crustal deformation, we assume an initial state of the plate boundary between two opposing approaches. One end member approach assumes a complete stress and deformation discontinuity across a narrow zone (the 'S point' models are related to this approach), and the other end member accounts for lateral continuity of stresses across the plate boundary (allowing for lithospheric scale folding).

The mechanical properties of the lithosphere (on a geological timescale) are examined and compared to several available data sets, in particular the image of the crustal root deduced from seismic velocity modeling [START_REF] Scherwath | Crustal structure from seismic data in the vicinity of the Alpine Fault[END_REF], Figure 1b) and contemporary GPS measurements. Instead of presenting an exhaustive collection of numerical models, this study demonstrates how a well constrained numerical model can bring insight into the development of a young continental collision zone such as New Zealand Southern Alps.

NUMERICAL MODELLING 1 Numerical approach

General features

We use the mixed finite elementsfinite differences code Parovoz [START_REF] Poliakov | Diapirism and topography[END_REF]. It is based on the Fast Lagrangian Analysis Continuum method (FLAC, Cundall & Board, 1988), which incorporates an explicit timemarching scheme, and allows the use of a wide range of constitutive laws such as brittleelasticductile rheology derived by rock experimentalists (for example [START_REF] Brace | Limits on lithospheric stress imposed by laboratory experiments[END_REF]Kohlstedt, 1980, Ranalli, 1995). It handles initiation and propagation of nonpredefined faults (shear bands). The present version of Parovoz includes diffusive surface processes, heat advection and conduction including initial age dependent temperature field. Parovoz has proved efficient in modeling many tectonic features, such as salt diapirism (Poliakov et al., 1993), lithospheric folding [START_REF] Gerbault | Do faults trigger folding in the lithosphere?[END_REF], lithospheric rifting [START_REF] Buck | Abyssal hills formed by stretching oceanic lithosphere[END_REF][START_REF] Burov | Erosion and rheology controls on synrift and postrift evolution: verifying old and new ideas using a fully coupled numerical model[END_REF]. More details are given in the Appendix concerning the numerical method, calculation of the geotherm, and rheological laws.

The lithosphere is modeled as a medium of 300 by 30 quadrilateral elements, with a total length of 600 km and depth of 60 km (Figure 2). Both lateral borders are free to slip vertically, the left border is fixed in the horizontal direction and a velocity V x =4.5 10 10 m/s (15 mm/yr) is applied on the right border. The lithosphere floats on the asthenosphere within the gravity field; hydrostatic boundary conditions are applied at the bottom of the model, with an underlying density of 3250 kg/m 3 . The surface is stress free.

An initial temperature distribution is calculated, which includes a thermomechanical age dependency, analogous to that of oceanic lithospheres, and incorporates crustal radioactive heat (Burov & Diament, 1995, see Appendix). A thermomechanical age equal to 100 Ma is assigned here, based on the last episode of extension that affected the New Zealand shelf [START_REF] Laird | Cretaceous continental rifts: New Zealand region[END_REF] and was likely to have rejuvenated most of the lithosphere.

The model lithospheric plate is divided into areas of different mechanical properties (Table 1). The initial crustal thickness, or the depth to the Moho, is set at 25 km; this choice is constrained by several observations from seismological studies, indicating a 'far field' crustal thickness of 22 to 26 km, [START_REF] Scherwath | Crustal structure from seismic data in the vicinity of the Alpine Fault[END_REF][START_REF] Godfrey | Crustal structure and thermal anomalies of the Dunedin Region, South Island, New Zealand[END_REF]. The density distribution is 2650 kg/m 3 , 2900 kg/m 3 and 3200 kg/m 3 respectively for the first 19 km, from 19 to 25 km, and from 25 to 60 km depth.

Rheological setup

Elasticviscousbrittle behavior is modeled with a pressuredependent Coulomb criterion for brittle failure, and a temperaturedependent creep powerlaw for ductile behavior (see Table 1 and Appendix). For each element of the model and at any timestep, if the shear stress provided by the creep law reaches the MohrCoulomb plastic yield, then nonassociative failure occurs. In this sense, the brittleductile transition is selfconsistently defined in the modeling.

Above the Moho, a vertical division separates the Australian crust, to the left (west) from the Pacific crust to the right (east). The 'Australian crust' and 'Pacific crust' rheologies are chosen to fit geological and geophysical data. The Australian crust is composed of Paleozoic plutonic rocks (Western Province) that are found in northwest South Island, in Fiordland, and offshore on the Challenger and Campbell Plateaus (Figure 1a). There is no evidence for internal deformation in the Australian lithosphere, apart from a broad flexural down warp towards the Alpine Fault [START_REF] Sircompe | The South Westland Basin; seismic stratigraphy, basin geometry and evolution of a foreland basin within the Southern Alps collision zone, New Zealand[END_REF]. Powerlaw creep parameters for the Australian crust are thus chosen so that the lower crust remains elastobrittle (Table 1). East of the Alpine Fault, the Eastern Province rocks are composed of lowgrade metasediments (greywacke to schist) that were accreted in Mesozoic times. These rocks are younger than the Western province rocks, and are likely to be mechanically weaker. A dominant quartzite rheology (Table 1) best matches the observed geology, the seismological brittleductile transition (BDT) at 12 km [START_REF] Leitner | A focused look at the Alpine Fault, New Zealand, focal mechanisms and stress Observations[END_REF], the convergence rate, and constraints on the geotherm.

In between the Australian and Pacific crust, a zone of weakness is inserted (4 elements wide and located at X = [0; 10] km), and simulates the Alpine Fault (Figure 2). Inclusion of a preexisting fault zone derives from the observation that Neogene plate boundary deformation along the Alpine Fault is most likely an inherited Eocene and older structure [START_REF] Sutherland | Plate boundary deformation in South Island, New Zealand, is related to inherited lithospheric structure[END_REF]. If this preexisting fault zone is not present, then a folding instability develops in the Pacific crust. The strength (or viscosity) contrast between the brittle crust and the lower ductile crust only needs to be of about 3 orders of magnitude for folding to develop after less than 5% of shortening [START_REF] Martinod | Periodic instabilities during compression or extension of the lithosphere 1. Deformation modes from an analytical perturbation method[END_REF][START_REF] Burov | The effective elastic thickness of continental lithosphere: what does it real ly mean?[END_REF][START_REF] Gerbault | Do faults trigger folding in the lithosphere?[END_REF].

Below the Moho, dry olivine is assumed to control the strength of the mantle lithosphere (see creep law parameters in Table 1). Strainsoftening is included, with details in the Appendix: its effects are discussed later in the paper. A zone of weakness must also be inserted in the mantle, or else periodic undulations develop, due to periodic folding of the 'strong' mantle lithosphere. Preliminary numerical tests show that both assumptions of a preexisting weak zone or a thermal perturbation lead to similar results. Since a thermal anomaly would better correspond to extra heat present at depth due to sustained shear along the trend of the Alpine Fault, a perturbation of 100 o C is inserted below the crustal weak zone, at X = [10; +20] km and Y = [40;60] km.

Results after 7 Ma of shortening are presented in Figures 3 and4 for the reference model. In these figures, the origin of the horizontal axis (X=0) is located at the initial preexisting crustal weak zone. Comparisons with electrical conductivity, geodesy, seismic Moho and gravity anomaly are shown in Figures 5,6, and 7. The following sections discuss deformation in the crust, the role of lithospheric mantle strength on the geometry of the orogen, and the relative position of the maximum topography and the crustal root, compared to observations across the Southern Alps.

Deformation in the crust

Surface uplift

Surface deformation is displayed through time in Figure 3. Maximum elevation develops above the initial weak zone (at X=0), evolves 10 km towards the west, and reaches about 3000 meters within 7 Ma (Figure 3a). A secondary maximum develops close to X=100 km (this position reduces as time progresses).

The whole lithosphere is affected by a component of regional surface uplift (Figure 3b), due to homogeneous thickening (or volumetric or pureshear thickening). Homogeneous thickening of an elastic or viscous plate is a fundamental way to accommodate compression, complementary to localized shear and mountain growth [START_REF] Burov | Continental lithosphere folding in Central Asia (II): constraints from gravity and topography[END_REF][START_REF] Burg | Lithospheric scale folding; numerical modelling and application to the Himalayan syntaxes[END_REF][START_REF] Gerbault | At what stress level is the Indian Ocean lithosphere buckling ?[END_REF]. This regional or mean surface uplift (different from rock exhumation) is around 0.15 mm/yr, with a maximum of 5 mm/yr (3000 meters in 7 Ma, Figure 3a, 3b), with respect to the initial zero level. In reality the mean density of the lithosphere probably increases with increasing distance away from the South Island (as the lithospheric rock composition becomes more similar to oceanic lithosphere). Our model cannot fit real far field surface elevations, because such lateral density contrasts are not taken into account.

The maximum erosion rate is located about 10 km west of the maximum of the topography (Figure 3c). It reaches a quasistable value of about 5 mm/yr after 2 Ma (Figure 3d). The sedimentation rate is lower than the erosion rate, partly because of the open mass flow condition at the borders of the model, and partly because of the numerical approximation of sedimentation (see Appendix).

Sediment budget studies (for example, [START_REF] Griffiths | Distribution of mean annual precipitation across some steepland regions of New Zealand[END_REF][START_REF] Adams | Age of metamorphism and uplift of the Haast schist Group at Haast Pass, Lake Wanaka and Lake Hawea[END_REF] and cooling ages [START_REF] Tippett | Fission track analysis of the late Cenozoic vertical kinematics of continental Pacific crust, South Island, New Zealand[END_REF], both indicate that uplift rates increase exponentially up to about 10 mm/yr at the Alpine Fault. There is also an indication for the initiation of rapid erosion in the Southern Alps beginning at 23 Ma [START_REF] Tippett | Fission track analysis of the late Cenozoic vertical kinematics of continental Pacific crust, South Island, New Zealand[END_REF]. [START_REF] Koons | The topographic evolution of collisional mountain belts: a numerical look at the Southern Alps, New Zealand[END_REF][START_REF] Koons | Twosided orogen; collision and erosion from the sandbox to the Southern Alps[END_REF] studied the role of asymmetric erosion and recognized that rapid removal of material from the west coast controls the asymmetry of deformation on either side of the Main Divide, a point developed also by [START_REF] Beaumont | Erosional control of active compressional orogens[END_REF][START_REF] Beaumont | The continental collision zone South Island, New Zealand: comparison of geodynamical models and observations[END_REF].

The factor of two difference in the erosion rate between our models and the observations is due to surface processes being modeled with a numerical algorithm that smoothes surface deformation. This is also the reason why fluvial transport is not incorporated in the modeling, as developed by [START_REF] Koons | Twosided orogen; collision and erosion from the sandbox to the Southern Alps[END_REF] and [START_REF] Beaumont | Erosional control of active compressional orogens[END_REF][START_REF] Beaumont | The continental collision zone South Island, New Zealand: comparison of geodynamical models and observations[END_REF], because it is 'erased' by this numerical effect.

Temperature, strain and stress distribution

Figure 4b displays the deflection of initially horizontal isotherms after 7 Ma of shortening. Isotherms are uplifted in the upper crustal area just east of the modeled 'Alpine Fault zone' (see vertical dashed line referring to its modeled position on Figure 4). They are deflected downwards at midcrustal depths. This geometry is very similar to results obtained from previous thermal modeling [START_REF] Shi | Thermal modelling of the Southern Alps[END_REF], dashed lines on Figure 4b), which was compared to heat flow measurements. This supports the fact that advection of heat is the main mode of heat transfer as a response to 7 Ma of compression.

Figure 4c displays distribution of the shear stress σ s = ((σ xx σ yy ) 2 /2 + σ xy 2 ) 1/2 . One can identify three strong layers: the Australian lithosphere, the Pacific upper crust, and the Pacific lithospheric mantle. Note the geometry of weak zones that develop in the center of the model, at the plate boundary.

Figure 4d displays the accumulated total shear strain ε s = ((ε xx ε yy ) 2 /2 + ε xy 2 ) 1/2 , with time, while Figure 4e displays the instantaneous shear strain rate and instantaneous velocity vectors. Deformation accumulates in a 'crustal root', which progressively increases in width and depth. Lowstrength Pacific material 'piles up' against the Australian lithosphere, connecting (i) the 'Alpine Fault zone' in the upper crust, emerging to the surface at X~ 30 km, (ii) a horizontal shear zone in the ductile lower crust, extending towards the east, and (iii) a 'westdipping' shear zone propagating downwards into the lithospheric mantle.

Shear strain localization the Alpine Fault

The 'Alpine Fault' shear zone develops at the emplacement of the initial vertical weak zone in the crust. This zone progressively inclines and propagates into the Australian crust and towards the surface (compare Figures 4e with inclined shear zone, and 4f, which shows rheological layers and vertical initial plate boundary zone). Its geometry is determined by two conditions: (1) the preexisting vertical weakzone, (2) the development of a shear zone inclined between 30 o and 45 o from the maximum principal strain and principal stress orientation (nonassociated plasticity with friction angle φ=30 o and dilatation angle ψ=0°, [START_REF] Vermeer | The orientation of shear bands in biaxial tests[END_REF][START_REF] Gerbault | Prediction of faulting from the theories of elasticity and plasticity, what are the limits ?[END_REF].

The 'Alpine Fault' shear zone is wider than the mapped fault: the numerical model produces shear bands that occupy a width of about 5 numerical elements. Since the present model deals with 2 km wide elements, the modelled Alpine Fault zone is thus 10 km wide. Because the initial vertical weak zone is not in the 'correct' orientation to accommodate thrust faulting inclined at 30° to 45°, the active shear zone is superimposed only on part of the initial weak zone: this superimposition (Figure 4d, 4f) occurs around 10 15 km depth, where the shear strength is highest with depth, but lowest in horizontal direction. This geometry is logical if one considers how a fault would develop when using smallest energy patterns.

The initial vertical weak zone in the model and after 7 Ma does not represent the present day transform shear zone. 2D modelling cannot address the issue of strainpartitioning. However, 3D modelling (Gerbault et al., in press) also indicates a progressive inclination with time of the vertical plate boundary, and along which transcurrent shear occurs. This is consistent with indications of nonpartitioning in the Southern Alps (e.g. [START_REF] Braun | Three dimensional numerical experiments of strain partitioning at oblique plate boundaries: implications for contrasting tectonic styles in California and South Island, New Zealand[END_REF].

Shear strain localization general features

The modelled Alpine Fault shear zone curves at depth to merge with a ductile shear zone in the lower crust (Figures 4d,4e). At this level, deformation propagates 'backwards', to the east, within the approximately 10 km thick low viscosity lower crust, and allows the upper crust to detach from the lithospheric mantle. In section 3, we illustrate how the eastern tip of this lower crust ductile shear zone can be connected via uppercrustal faulting to the surface, at a characteristic distance ranging at 100 to 130 km from the Alpine Fault shear zone.

In contrast to the development of localized shear strain at the Alpine Fault zone, diffuse faulting in the brittle upper crust is represented numerically by shear bands `jumping' in both dip directions in time and space (Figure 4e). This model behavior may correspond to the abundant reactivated faults observed within the Torlesse greywackes east of the Alpine Fault [START_REF] Oliver | Sheet J36D and part sheet J35 Clearwater[END_REF][START_REF] Cox | The Main Divide fault Zone and its role in the formation of the Southern Alps, New Zealand[END_REF], Long et al., 2001).

Figure 5a displays the shear strain using ellipses: each ellipse represents how an initially perfect circle is deformed through time. Darkest ellipses are those that undergo the most shear strain (same color distribution as Figure 4d). Figure 5b shows shearstrain ellipses representation for a model identical to the reference model, but using a more refined grid (400x40 elements). A higher resolution allows us to reproduce a narrower 'Alpine Fault zone' and examine lower crustal deformation in more detail. The high resolution models could take more than a week to run, thus preventing us from establishing them as a reference case.

In the crustal root, progressive rotation of the highstrain ellipses demonstrates how the Pacific crust diverges in response to indentation by the Australian plate (the sharp western border of high strain ellipses coincides with rheological discontinuities inserted in the model). The Alpine Fault zone consists of multiple ellipse directions, associated with thickening of the crust, and fault shearing (Figure 5).

Inclined ellipses located at the base of the lower crust indicate that mantle lithosphere drives material from beneath the crustal root. This shows that although a uniform velocity is applied at the east end of the model, the mantle and the crust move with different rates with respect to each other from about 130 km east of the plate boundary (Figure 5).

In the lower crust at about X ~ 25 km, 60 km east of the plate boundary, the shear orientation of ellipses shifts, with locally a horizontal orientation corresponding to less horizontal compression. Above, brittle deformation accumulates in the upper crust, in a zone located about 60 km east from the Alpine Fault zone: it is the result of uplift of the upper crust, as it is tilted (or flexed), and dragged with the ductile lower crust (Figures 4d,5). Tilting and flexure of the upper crust is consistent with the fanshape of geological strain markers [START_REF] Little | Kinematics of oblique continental collision inferred from ductile microstructures and strain in midcrustal Alpine Schist, central South Island[END_REF]).

The modeled strain, observed faults and electrical conductivity in the Southern Alps

The localisation of high strain about 60 to 80 km southeast of the Alpine Fault (Figure 5) is in very good agreement with a series of active structures, the Lake HeronForest Creek Faults and the Irishman Creek Fault, which are structurally similar faults that dip either east or west [START_REF] Beanland | Field guide to sites of active earth deformation[END_REF][START_REF] Oliver | Sheet J36D and part sheet J35 Clearwater[END_REF][START_REF] Woodward | Seismic tests for the South Island Transect[END_REF]. This characteristic distance also coincides with the eastern edge of the Bouguer gravity anomaly.

Electrical resistivity models derived from a magnetotelluric profile experiment [START_REF] Scherwath | Crustal structure from seismic data in the vicinity of the Alpine Fault[END_REF] show a Ushaped high conductivity structure under the Southern Alps [START_REF] Wannamaker | Fluid generation and pathways beneath an active compressional orogen, the New Zealand Southern Alps, inferred from magnetotelluric data[END_REF]. The highest conductivities (40100 Ohm/m) lie at 25 to 30 km depth below the surface (Figure 5) and extend upward as narrownearvertical conductors near the Alpine Fault trace and 6070 km east of the Alpine Fault. This latter branch of the conductor appears to coincide with fault zones (the Forest CreekLake Heron Faults), that would allow deep crustal fluids to migrate to the surface [START_REF] Templeton | stable isotopic evidence for mixing between metamorphic fluids and surfacederived waters during recent uplift in the Southern Alps of New Zealand[END_REF]. Among the deformation mechanisms that are expected to operate in lower crustal ductile rocks, diffusion creep may dominate when the volume of fluids is small (<0.5 vol %), and under relatively low stress conditions [START_REF] Tullis | Diffusion creep in feldspar aggregates: experimental evidence[END_REF]. [START_REF] Wannamaker | Fluid generation and pathways beneath an active compressional orogen, the New Zealand Southern Alps, inferred from magnetotelluric data[END_REF] discuss the cause of such a high conductive signal and, while the presence of graphite cannot be dismissed, these authors favor the presence of interconnected fluids resulting from prograde metamorphic reactions occurring in the crustal root.

In the present numerical models, we have not accounted for the presence of fluids being produced and evolving with time. However, the superposition of the approximate shape of the high conductor on the strain ellipses, in Figure 5, indicates a plausible link of conductivity as an image of the high shear zones in the crust. The change in ellipse orientations at about X= 25 km, is interpreted as a zone of less compression, due to clockwise bending of the upper crust to the west of this zone, and in which interconnected fluids may accumulate.

Figure 5c represents stress lines, each line showing the local orientation of the most compressive stress. This figure suggests that the conductivity signal is caused by inner crust deformation mechanisms, because it is located above and to the west of the change in stress orientation corresponding to extensional flexure of the upper mantle. The superposition of earthquakes along our 2D model (from [START_REF] Leitner | A focused look at the Alpine Fault, New Zealand, focal mechanisms and stress Observations[END_REF] shows good agreement with areas of elevated compression and also seem to coincide with local stress orientation changes.

Comparison to contemporary plate motion, and superficial width of the orogen

Figure 6 displays a comparison between the observed GPS convergence velocity and the result from the model. This modeled convergence velocity is evaluated after about 7 Ma of shortening, for one timestep, which is typically of the order of 30 years. We see that the reference model has a generally higher convergence rate than observations (thick light grey line on Figure 6). This is because the reference model is defined to be extending for 400 km east of the plate boundary, and to have a convergence rate of 15 mm/yr (these values were used firstly, to avoid numerical border effects on the resulting deformation, and secondly, to match plate tectonic models that indicate about 100 km of convergence within 7 Ma). Present convergent rates based on GPS observations [START_REF] Beavan | Contemporary horizontal velocity and strainrate fields of the PacificAustralian plate boundary zone through New Zealand[END_REF], however, record a velocity of about 8 mm/yr, as measured across the South Island's width of 200 km, and for a time period of 10 years. While our aim in this paper is not to explain this mismatch of convergence velocity from geophysical data, we also plotted the resulting velocity for an alternative model, in which the distance from the border of the model to the plate boundary has been reduced to 200 km, and a convergence rate of 8 mm/yr imposed. After 7 Ma, the crust and mantle deformation from this alternative model are similar to those from the reference model, except that the crustal root is less deep, since the shortening rate is slower. Comparison of the surface velocity with the GPS data are now more consistent (thin dark grey line on Figure 6).

Two interesting features are noticeable on Figure 6: first, we note a similar stabilisation of convergence velocity, in the data and in the models, at about 80 to 100 km from the Alpine Fault: this feature can be related to the area where the crust detaches from the mantle lithosphere. Second, the breadth of increasing convergence velocity at the Alpine Fault is broader in the data than in the model: this is due to elastic strain accumulating from the last earthquake cycle. This velocity increase, over time, would become sharper in GPS measurements, and thus more similar to the models.

About 100 km east of the Alpine Fault in the central Southern Alps, the topography increases to about 1500 meters. Many strikeslip and reverse faults coexist in this area. Previous analytic, analog and numerical models also reproduce this zone of secondary surface uplift, which represents conjugate brittle shear with respect to the Alpine fault, i.e. the other edge of the doublevergent wedge characteristic of orogenic process [START_REF] Koons | Twosided orogen; collision and erosion from the sandbox to the Southern Alps[END_REF][START_REF] Koons | Geodetic analysis of model oblique collision and comparison to the Southern Alps of New Zealand[END_REF][START_REF] Beaumont | The continental collision zone South Island, New Zealand: comparison of geodynamical models and observations[END_REF]. Our models emphasise that the amplitude of this secondary zone of uplift depends on a number of parameters. Firstly, erosion and sedimentation processes are known to help localise deformation [START_REF] Beaumont | The continental collision zone South Island, New Zealand: comparison of geodynamical models and observations[END_REF]. If surface processes are prescribed to be slow, deformation is less asymmetric and diffuses over a broader area, to the east in the case of the Southern Alps. Secondly, if the lower crust is assumed to be relatively stronger than in the reference model, then ductile shear deformation is confined to a smaller depth range, and tends to propagate more to the east, in order to occupy the same volume. If the Moho is chosen initially deeper (say at 30 km), then the distance at which this secondary uplift occurs increases. Finally, a strong mantle lithosphere also contributes to producing wider surface deformation. This is discussed in following section (3).

3 Deformation in the lithospheric mantle

Observations of deformation at and under the Moho in the Southern Alps

Seismic velocity models [START_REF] Davey | Preliminary results from a geophysical study across a modern, continentcontinent collisional plate boundarythe Southern Alps[END_REF], Scherwath et al., 2000, EberhartPhillips & Bannister, 2001) provide an image of a crustal root about 100 km wide, with a maximum depth of 40 to 45 km, and increasing southward between Transect 1 and Transect 2 (Figure 1a,b). The deepest part of the crust lies about 20 km east from the Main Divide (itself about 20 km east from the Alpine Fault). These thickness and location of the crustal root indicate that the orogen is not in isostatic equilibrium, and is corroborated by isostatic gravity anomalies [START_REF] Woodward | The crustal structure of the Southern Alps, New Zealand, as determined by gravity[END_REF]). In addition, teleseismic Pwaves delays indicate a zone of highvelocity in the mantle below the Southern Alps [START_REF] Stern | Teleseismic P wave delays and modes of shortening the mantle lithosphere beneath South Island, New Zealand[END_REF]Kohler & EberhartPhillips, submitted). The observation of up to 1 second advance in Pwave arrivals (recorded from three Western Pacific earthquakes) is satisfied with a highspeed body located at about 120 km depth and 100 km wide, with a maximum speed advance of ~7% [START_REF] Stern | Teleseismic P wave delays and modes of shortening the mantle lithosphere beneath South Island, New Zealand[END_REF]. [START_REF] Stern | Teleseismic P wave delays and modes of shortening the mantle lithosphere beneath South Island, New Zealand[END_REF] argue that this highspeed zone is due to highdensity mantle, which would also explain observed gravity anomalies.

In the following we compare our 'reference' model to the observed surface and Moho in the Southern Alps, and then discuss the role of lithospheric mantle rheology, with support of additional models that were modified from the reference model. How these results compare with the highdensity mantle model developed by [START_REF] Stern | Teleseismic P wave delays and modes of shortening the mantle lithosphere beneath South Island, New Zealand[END_REF] is addressed in the conclusions.

Comparison to the seismic Moho and gravity

Figure 7 displays the surface relief (a) and the Moho (b) across the Southern Alps, superimposed on the same horizons obtained from the reference model (thin gray line). The shape of the crustal root (Figure 7b) is satisfactorily reproduced. It is difficult to match these horizons in the far field because we do not account for progressive change in lithospheric rock composition, away from the plate boundary. The modeled topography is consistently higher than observed since uniform uplift acts in our model by raising the entire surface horizon. Since level 0 of topography corresponds to a regional value, corresponding to far field isostatic equilibrium that the model cannot take into account, we empirically subtracted 500 meters on its representation Figure 7a.

The model shows no significant offset between maximum topography and the deepest portion of the crustal root. In the present reference model, it is even slightly to the east of the crustal root; in alternative models (not shown) with a slightly weaker Australian upper crust, we obtain a maximum topography located about 10 km west of the crustal root. The cause for this absence of offset in the model, when compared to observations, is further discussed in section 3.6.

Figure 7c displays the comparison between observed and modeled gravity anomaly, the calculation of which is detailed in the Appendix. The density contrast (300 kg/m 3 ) at the Moho contributes about 50 % of the gravity anomaly , the density contrast between the upper and the lower crusts (250 kg/m 3 ) contributes about 30 %, and the density contrast at the base of the model (50 kg/m 3 ) contributes about 15 %, of the total gravity anomaly. 100 mGals were added to the overall modeled anomaly, which is again justified by far field isostatic equilibrium and the reference geoid level not being taken into account in our models. The model and observed data do not match on the eastern side of the crustal root because of the presence of a thick Cenozoic sedimentary basin beneath Canterbury plains and extending offshore [START_REF] Field | Cretaceous and Cenozoic sedimentary basins and geological evolution of the Canterbury Region, South Island, New Zealand[END_REF]), which has not been taken into account in the initial setup of our model. The mismatch in detail to the west of the crustal root is most likely due to density variations in the Australian crust.

Seismological data indicate a seismogenic zone in the Pacific plate at a depth from 25 km to 60 km towards the west, which is interpreted as old oceanic crust and mantle lithosphere being tilted down towards the Australian plate [START_REF] Reyners | Subcrustal Earthquakes in the central south Island, New Zealand, and the root of the Southern Alps[END_REF]. This seismogenic zone (high shear strain in a brittle domain) is most likely to be associated with local mineralogical heterogeneities rather than a layer of significantly higher strength than the surrounding ductile crust or lithospheric mantle.

Remarks on the rheology of the lithospheric mantle

Lateral density contrasts between the Australian and Pacific plates can obviously be a first order cause of the direction of 'mantle subduction' in the Southern Alps of New Zealand. Although the relative position of initial weaknesses in the mantle and in the crust is also crucial, no argument justifies an asymmetric positioning. But the mantle's strength and rate of strain softening are also first order controls on the geometry of mantle lithosphere thickening, as recently illustrated by [START_REF] Pysklywec | Modeling the behaviour of the continental mantle lithosphere during plate convergence[END_REF].

There are several ways in which the strength of the mantle lithosphere can vary. For example, a warmer geotherm tends to induce lower viscosity, due to the temperature dependent powerlaw constitutive behaviour. Parameters of the powerlaw may vary when extrapolated to real mantle depths, and other mechanisms of deformation than dislocation creep may even dominate (see [START_REF] Ranalli | Rheology of the Earth[END_REF]. The choice of dry versus wet olivine creep parameters modifies the depth of the brittleductile transition in the mantle lithosphere (e.g. Rannalli, 1995;Pysklywec et al., in press), and thus the thickness of its 'strong' layer: consequently the lengthscale of mantle deflection is modified. With the large uncertainties in a realistic constitutive law for the mantle lithosphere, and after having tested about 300 models, we choose to demonstrate the control of strength and strainsoftening in a simple manner, and regardless of which mechanism controls this strength in detail: since maximum strength is bounded by both friction and cohesion, we prescribe zero friction and vary the cohesion (see table 1). However, this cohesion should not be considered as a physical cohesive value for rocks.

The effect of lithospheric mantle's strength and the symmetry of its deformation

From a mechanical point of view, Moho topography (i.e. wavelength and amplitude) is dependent on the ability of mantle lithosphere to deform. The stronger the lithospheric mantle (either viscous or brittle), the easier it is to deflect because it can support large internal stress contrasts responding to lateral loading. This is similar to the mechanism of lithospheric folding (e.g. [START_REF] Gerbault | Do faults trigger folding in the lithosphere?[END_REF], in which the strong mantle lithosphere (dark layer in Figure 4c) can buckle 'plastically' with a wavelength proportional to about 46 times its thickness, and an amplitude proportional to its strength ratio with the surrounding layers. In contrast, a mantle with relatively low viscosity deforms by uniform viscous thickening. In this case the amplitude of deformation is diminished, although deformation may be laterally more diffuse. The small amplitude of mantle deformation then less affects crustal deformation.

Models were produced for a mantle lithosphere strength of 500 MPa (not shown). As expected, broad flexure of the mantle lithosphere leads to broad surface deformation. The secondary maximum of topography, mentioned previously, tends to develop more easily, as it coincides with flexural uplift of the mantle lithosphere.

Figure 8 displays two models with an initial mantle strength of 50 MPa, reducing to 20 MPa. In the first model (Figure 8a), the 'Australian' mantle underthrusts the 'Pacific' mantle after 7 Ma of convergence. In the second model (Figure 8b), with a greater rate of softening and a lower coefficient of erosion (see figure caption for values), two conjugate shear zones develop so that mantle deformation remains symmetric. These models show that the combination of several parameters control the mode of mantle lithosphere thickening when it has a low strength (20 MPa is the value of cohesion); this is understandable when regarding the unpredictability of preferential movement along two conjugate thrust faults, which results from local equilibrium in mass transfer along vertical columns. However, none of our models with initial mantle strength of 50 MPa allow for a large width and amplitude of crustal thickening, because there is insufficient build up of stress to support significant mantle lithosphere deflection. Models with an initial mantle strength of 50 MPa produce a crustal root about 50 km wide and less than 35 km deep, in contrast to models with an initial strength greater than 200 MPa, in which 'Pacific plate' subduction systematically develops.

In our models, (a)symmetry of the crustal root develops associated with (a)symmetry of entire lithospheric thickening.The symmetrical deformation of Figure 8b is in agreement with indications from P wave delays for symmetric lithospheric thickening. However, the significantly smaller shape of the Moho renders this model unsatisfactory. More complex constitutive laws, such as temperature dependent strain softening [START_REF] Zhang | Simple shear deformation of olivine aggregates[END_REF], may be required in order to reproduce asymmetric crustal root and symmetric mantle lithosphere thickening, as indicated by seismological analyses in the Southern Alps.

The effect of strain softening

If strain softening is too rapid, deformation remains localized over a narrower region than indicated by observations from the Southern Alps. An extreme situation is that deformation remains concentrated in the initial vertical zones of weaknesses. In contrast, if strain softening is imposed too slowly, then lithospheric thickening is too wide. In absence of strainsoftening, largescale folding of the Pacific mantle develops.

Models of initial mantle strength equal to 50 MPa (Figure 8) illustrate the importance of strain softening on the geometry of mantle deformation. To further illustrate its effect on the width of the orogen for a stronger mantle lithosphere, Figure 9 depicts a case where strainsoftening is imposed more slowly than in the reference model (reduction from 200 MPa to 50 MPa after 200% of shearstrain rather than 100%). Broad surface deformation occurs, with development of another maximum of topography located about 120 km east of the Alpine Fault zone. Shearstrain propagates to the east into the ductile lower crust, and emerges to the surface via thrust faults: the secondary maximum of topography becomes as high as the primary one after 6 Ma of shortening.

Discussion with respect to the 'S point' models

Uniform lithospheric thickening can be approximated by considering simple elastic compressibility: results of models for a purely elastic lithosphere provide after 5 Ma and 75 km of shortening, surface uplift y t = 750 m (uplift rate of 0.15mm/yr) and bottom subsidence y b =2000 m (thickening of 4mm/yr). By assuming a linear dependency of the vertical displacements with depth, D(y) = y t + (y t y b )y/h o , the depth at which material diverges upwards or downwards (D(y)=0) is found to be located at the midcrustal depth S=16.4 km.

Plastic and viscous behaviors trigger heterogeneous deformation that lead to continuous readjustment of the 'S point' position, through time and space: the S point evolves within the crustal root, and would be best represented by a diffuse zone, rather than a point. Our numerical models show that the present geometry of upper crustal deformation in the Southern Alps is not very sensitive to the strength of the mantle lithosphere, because it is less than 200 MPa (and undergoes sufficiently rapid strainsoftening). In this sense, Beaumont et al.'s (1996) assumption of 'S point' discontinuity at the base of the crust is a realistic approximation, as long as the models' interest is not in the exact shape of the base of the crust and deeper deformation.

In the present numerical models, the maximum of the topography and the deepest part of the crustal root remain on the same vertical axis (see Figure 7). The observed ~20 km offset was never satisfactorily reproduced, despite a large number of parametrical tests, including a variation of rheological properties, a varying coefficient of erosion between west and east of the Main Divide, and a tuning of strain softening coefficients. An offset of about 10 km at best was obtained when the friction angle for the Australian crust is reduced to 10 o , rather than 30 o . The reason for the model not producing significant offset of surface and crustal root deformation may be that only diffusive surface processes were taken into account in our models. When Beaumont et al. (1996, model 9) account for strain softening and asymmetric erosion processes, a significant offset is obtained (see also [START_REF] Willett | Orography and orogeny: The effects of erosion on the structure of mountain belts[END_REF]. However, we suggest that the achievement of this offset by [START_REF] Beaumont | The continental collision zone South Island, New Zealand: comparison of geodynamical models and observations[END_REF] is probably due to the overconstraint imposed by their basal velocity discontinuity.

The observation of increasing offset of the Bouguer gravity anomaly with respect to the plate boundary, towards the southeast, where compression and surface topography decrease, indicates that this offset may find its cause in the three dimensionality of the orogen. Recent 3D numerical modeling (Gerbault et al., in press) supports the occurrence of southward ductile lateral flow within the lower crust and along the plate boundary. Crustal thickening tends to develop towards the southeast, showing that the overall deformation results from a tendency to equilibrate buoyancy forces in both the vertical section and laterally, according to pressure variations. This thickening results from the geometry of rheological discontinuities with respect to the orientation of the relative plate velocity vector, and coincides with the shape of the Bouguer anomaly.

CONCLUSION

Our numerical models involve compression of the entire continental lithosphere, during the first 10 Ma of orogenic growth. While the models presented are nonunique, they are, nonetheless the culmination of a number of tests that provide a mechanism explaining crustal deformation data (surface structures, GPS, electrical conductivity, seismic velocities). They also demonstrate the important control of the strength of the lithospheric mantle on the geometry of the collision process. Results are compared to the Southern Alps of New Zealand and discussed in respect to previous studies:

1) Advection appears to be the main mechanism that accounts for the deflection/uplift of isotherms after 7 Ma of shortening (Figure 4b). This confirms the validity of the thermokinematical approach of [START_REF] Shi | Thermal modelling of the Southern Alps[END_REF].

2)

Crustal deformation is characterized by a major shear zone, that one can associate with the Alpine Fault, which from being initially vertical, inclines to about 45 o and curves at depth to merge with a ductile shear in the lower crust (Figure 4e). This illustrates the inclination through time of an initially vertical transcurrent shear zone, and is consistent with 3D modeling (Gerbault et al., in press).

3)

Flexure of the brittle upper crust above the crustal root is responsible for faulting being concentrated 6080 km east from the Alpine Fault (Figure 5) and is coincident with several active faults and with the boundary of the observed Ushape electrical conductor [START_REF] Wannamaker | Fluid generation and pathways beneath an active compressional orogen, the New Zealand Southern Alps, inferred from magnetotelluric data[END_REF]. This result, consistent with the argument that midcrustal fluids accumulate in a locally extended zone (less compressed) in the lower crust and then migrate to the surface, finds further support here in our dynamical modeling.

4)

Comparison with GPS data of convergence velocity shows that the change in slope at a distance of 80 to 100 km from the Alpine Fault can be related to the place where the upper crust detaches from the mantle lithosphere (Figures 5 and6). The increase in velocity at the Alpine Fault, in the model, is sharper than in GPS data, because geodetic measurements account for short timescale aseismic deformation.

5)

As compression increases, flow in the ductile lower crust propagates to the east, so that the upper crust detaches from the lithospheric mantle. The tip of this ductile shear zone, located about 100 to 130 km from the Alpine Fault, can connect to thrust faults in the upper crust, and generate surface uplift. This characteristic distance and the amount of surface uplift associated with it, are not only dependent on the (low) strength of the ductile lower crust, but also on the (strong) strength of the lithospheric mantle.

6)

The crustal root develops and resembles the image obtained from seismic velocity models in the central Southern Alps (Figure 7), with 'subduction of Pacific' mantle. If the strength of the lithospheric mantle is weaker than 200 MPa, the amplitude and wavelength of the crustal root are too small (Figure 8). If the rate of softening is too slow then the resulting surface deformation is too wide (Figure 9), a result that was also shown by [START_REF] Pysklywec | Modeling the behaviour of the continental mantle lithosphere during plate convergence[END_REF]. Satisfactory models with softening ranging from 200 MPa to 50 MPa are consistent with a number of geophysical studies that estimated maximum stress contrast that a lithosphere can support (e.g. [START_REF] Mcnutt | Implications of regional gravity for state of stress in the earth's crust and upper mantle[END_REF][START_REF] Molnar | Temperatures, heat flux, and frictional stress near major thrust faults[END_REF][START_REF] Gerbault | At what stress level is the Indian Ocean lithosphere buckling ?[END_REF]. The typical horizontal integrated force in our models is 26 10 12 N.m.

7)

A good fit is obtained to the measured gravity anomaly for our reference model, without needing to include localised denser mantle lithosphere based on measured Pwave delays [START_REF] Molnar | Continuous deformation versus faulting through the continental lithosphere of New Zealand[END_REF][START_REF] Stern | Teleseismic P wave delays and modes of shortening the mantle lithosphere beneath South Island, New Zealand[END_REF]. This difference warrants further investigation. What is the role of mechanical anisotropy in relating seismic velocities to density, is 100 km of shortening within 7 Ma sufficient to generate significant RayleighTaylor instability below the Southern Alps? Asymmetric and localised shear strain together with distributed symmetric mantle thickening may develop in continental collision zones, but would require additional detailed study. Our modelling indicates that for the Southern Alps, an asymmetric crustal root but symmetric lithosphere thickening would require more complex constitutive laws than those used here.

8)

The 20 km offset between the maximum topography and maximum depth of the crustal root were not reproduced in the models, despite a number of parametric tests. In addition, the Bouguer gravity low increases progressively to the southeast and away from the plate boundary, indicating that deep crustal thickening does not strike parallel to the plate boundary. This suggests that the dynamic equilibrium of the orogen is three dimensional. Although many of the observed deformational structures at an obliquely converging zone like the Southern Alps may be modeled using 2D orthogonal compression, structures that do not parallel strike require 3D modeling.

Finally, we modeled here the first stage of mountain growth, also called the linear growth stage (e.g. [START_REF] Shen | Largescale crustal deformation of the Tibetan Plateau[END_REF]. At times greater than 10 Ma, several studies emphasize the role of ductile, channel flow in the lower crust in monitoring the lateral expansion of the orogen [START_REF] Royden | Coupling and decoupling of crust and mantle in convergent orogens: implications for strain partitioning in the crust[END_REF][START_REF] Vanderhaeghe | Role of partial melting during orogenic evolution; insights from thermomechanical modeling[END_REF]. In our models, the change in shear direction within the ductile lower crust (see Figure 5) illustrates the complexity of ductile flow processes at the edge of a continental orogen, in its initial stages.

Table 1: Rheological properties for reference model: powerlaw creep parameters from [START_REF] Ranalli | Rheology of the Earth[END_REF]. Mohr Coulomb parameters are cohesion So, friction angle φ, and fluid pressure ratio λ (normal stress proportional to lithostatic pressure). Strain softening acts within both given values of the total shear strain, and modifies the value located on the same line in the previous column. Where specified, strain softening acts instead on viscosity µ. [START_REF] Poliakov | Diapirism and topography[END_REF], Poliakov et al., 1993) belongs to the Fast Lagrangian Analysis of Continua (FLAC) family [START_REF] Cundall | A microcomputer program for modelling largestrain plasticity problems[END_REF], Cundall, 1989), which is a large strain fully explicit timemarching numerical algorithm exploiting the Lagrangian moving grid method. The latter allows for solution of Newton's equation of motion in large strain mode, holding a locally symmetric small strain formulation commonly used in continuum mechanics. The method can reproduce initialization and evolution of nonpredefined faults (treated as large shear bands).

The timemarching scheme of the algorithm means that one loop represents one timestep, which must be small enough to prevent the physical information propagating from one element to its neighbour during this interval of time. Each loop contains the following procedure: velocities are calculated from Newton's law, with density ρ, time t, velocity vector v, and gravity acceleration g: ρ.δv i /δt =δσ ij /δx j + ρg. Components for the deformation rate are deduced from ? ij =1/2 (δv i /δx j + δv j /δx i ). Brittleelasticviscous constitutive laws then provide the stress distribution and equivalent forces ρ.δv i /δ t, which in turn provide input for the next timestep. The algorithm employs a dynamic relaxation technique based on the introduction of artificial masses in the dynamic system, permitting to handle strain localisation. The adaptive remeshing technique allows the models to be run up to about 25% of total shortening.

Brittleelasticviscous behaviour is modelled so that the minimum deviatoric stress produced by Mohr Coulomb elastoplasticity and Maxwell viscoelasticity is retained, at each timestep and for each element. Thus the brittleductile transition is selfdefined in the model, and is referred to as the depth at which the deviatoric shear stress becomes lower than about 20 MPa, regardless of which brittle or ductile behavior actually dominates (e.g. Rannali, 1995). The failure criterion is reached when the relationship between normal and tangential stresses, σ n and τ, along any given orientation equals τ=S o tanφ.σ n , S o being the cohesion and φ the friction angle for a given material property. Dilatancy angle ψ=0, so that nonassociated plasticity occurs (see [START_REF] Vermeer | The orientation of shear bands in biaxial tests[END_REF].

Viscous behaviour is prescribed by powerlaw dislocation creep, and Chen and Morgan's (1995) expression of the effective viscosity is used, with Q the activation energy, R the gas constant, A the material constant, and T the absolute temperature (see text for definition of σ s and ? s ):

µ* = ?.(4/3A) 1/n . ? s 1/n-1

.exp (Q/nRT) = σ s /2? s . When the stresses are evaluated, the increase in strain rate is greater than the decrease in effective viscosity, so that the shear stress actually increases in a shear zone. For this reason, in the Pacific crust, we introduce softening of the viscosity by a factor 10 within a certain range of the total shear strain ε s (see table 1).

Extrapolation of the powerlaw creep parameters for olivine in the upper mantle lithosphere may provide greater shear stress than the MohrCoulomb yield stress. Since dominant deformation processes at mantle depth are poorly known, a number of studies then assume a brittle behaviour (See [START_REF] Ranalli | Rheology of the Earth[END_REF][START_REF] Burov | The effective elastic thickness of continental lithosphere: what does it real ly mean?[END_REF]. For this reason, we chose to empirically assign zero friction and test the cohesion for the mantle lithosphere. Also, we test cohesion strain softening of the mantle lithosphere, as it appears necessary to allow for localised thickening (see Table 1 for parameters and section 3.3). Pysklywec et al. (in press) for example, apply friction strain softening. Cohesion or friction strain softening are achieved by reducing their value according to a linear equation that depends on the amount of total shear strain (see table 1). Softening of friction angle in the Australian crust is introduced to allow deformation to propagate to the right of the modelled plate boundary, and allows the maximum of topography to develop to the right of the thickened crust. Softening of the viscosity of the Pacific crust allows the deformation to remain localised close to the plate boundary.

Temperature is updated through time. The heat flux through an element of volume is equal to the heat produced by internal sources plus the variation of volume. If heat is advected with velocity V i , H is the heat production per unit mass, k is the thermal diffusivity, and ρ the density, the general equation of heat is [START_REF] Ranalli | Rheology of the Earth[END_REF]:

∂T/∂t + V i .∇T k ∇ 2 T = ρH.
The initial continental temperature field is calculated according to the age dependent procedure developped by [START_REF] Burov | The effective elastic thickness of continental lithosphere: what does it real ly mean?[END_REF], with a crustal heat production H = H s .exp(y/h r ), and H s =9.10 10 W/kg, h r =10 km, and crust and mantle thermal conductivities κ c = κ m = 3 W/K/m (κ = k.ρ.C p , C p = 1 kJ/kg/K the specific heat). For a lithosphere of thermal age of 100 Ma, a temperature T b = 1350 o C at the base of the lithosphere h l = 60 km, the resulting geotherm leads to the 600 o C isotherm lying at about 25 km depth. After 7 Ma of shortening, the far field surface heat flux is about 60 mW/m/K, increasing to about 120 mW/m/K at the topography maximum.

Erosion and sedimentation processes are taken into account with a classic equation of diffusion. The equation of mass conservation is written in terms of a linear relationship between the variation of surface elevation h with time, and the surface slope derivative. We adopt a constant value for the coefficient of diffusion k er (set to 1000 m 2 /yr), so that ∂h/∂t = k er .∂ 2 h/∂x 2 . If the area of an element increases by more than half of its initial value, then the element's 'phase' switches to a 'sediment' rheology (Table 1).

The Bouguer gravity anomaly is modelled from the progressive displacement of layers since time 0 when the model is assumed to be in equilibrium with no anomaly. The gravity anomaly relates to the vertical displacement dz of element (i) of width dx located at the boundary between layers of different density ∆ρ. Element (i) contributes to the gravity anomaly by dg i = 2G ∆ρ.dx.dz, where the gravitational constant G = 6.67 10 11 m 3 /kg/s 2 . Each element (k) at the surface of the model is thus subjected to anomalous gravity acceleration created by the displacement of elements (i) that compose the layer boundary (l) of contrasted density. Three layers (l) of contrasted density exist: the upper and lower crust boundary (∆ρ=250kg/m 3 ), the crust and the mantle boundary at the Moho (∆ρ=300kg/m 3 ), and the base of the model (∆ρ=50kg/m 3 ). For elements (i) at depth z i , located on (l), and with r i the distance between (k) and (i): gb (k) = Σ l Σ i (dg i . z i ) / r i 2 .

Figures captions

Figure 1: a) Geological setting in the Southern Alps of New Zealand [START_REF] Mortimer | Basement geology of South Island of New Zealand[END_REF]: Eastern province terranes on the Pacific plate deform against Western province terranes. Dark grey lines show approximate location of geophysical transects for the SIGHT experiment. b) Seismic velocity model across SIGHT line 2 from [START_REF] Scherwath | Crustal structure from seismic data in the vicinity of the Alpine Fault[END_REF]. Note that horizontal coordinates do not coincide with the following model coordinates. Vertical exaggeration is 4.

Figure 2: Initial conditions for the 2d model (see text). Figure 3: Surface deformation for reference model: topography (a, b, top) and erosion rates (c, d, bottom) plotted against horizontal distance (a,c, left) and time (b ,d, right). Symbols 'o' and 'x' correspond to the spatial maximum of topography after 3.2 and 7 Ma. After 7 Ma this maximum is located about 15 km east from the maximum erosion rate. Figure 5: Shear strain ellipses after 7 Ma. a) Reference model: numerical values of the shear strain are same as for Figure (4d). b) Higher resolution model (400x40elements). Superimposition of the high electrical conductivity signal in thick grey contour: note the 'flip' of shear strain ellipses in the ductile lower crust, at X~25km. c) Most compressive principal stress lines for the reference model. Grey circles are earthquakes distribution, dashed lines are presumed geometry of faults at depth, thick grey line is the 12 km deep seismogenic zone (all data from [START_REF] Leitner | A focused look at the Alpine Fault, New Zealand, focal mechanisms and stress Observations[END_REF].

Figure 6: Comparison of surface convergence velocity with GPS measurements [START_REF] Beavan | Contemporary horizontal velocity and strainrate fields of the PacificAustralian plate boundary zone through New Zealand[END_REF]. Thick light grey line is the reference model, with distance from the plate boundary to the eastern border of the model Xl = 400km, and applied convergence rate Vx = 15mm/yr. Thin dark grey line is an alternative model with smaller distance Xl=200km and slower convergence rate Vx = 8 mm/yr. Note the difference in velocity increase at the Alpine Fault, and the similar change in slope at a distance of about 80 km to the East.

Figure 7: Comparisons of reference model with observed a) Topography, b) Moho, and c) gravity anomaly. Black lines: observed topography, Moho from seismic velocity models along Transect 2 [START_REF] Scherwath | Crustal structure of the central South Island, New Zealand from an onshoreoffshore seismic data[END_REF], and measured Bouguer gravity anomaly. Grey lines: same horizons obtained with the reference model. In (a), vertical dashed lines are the observed (in black) and modelled (in grey) positions of the Alpine Fault (A.F.) and of the main divide (M.D.). 500 meters have been subtracted from the model topography (a), and 100 mGals have been added to the model gravity (c), (see text for justification).

Figure 8: Models with mantle lithosphere strength of 50 MPa, reducing to 20 MPa. Total shear strain (maximum in dark grey) and topography. a) Model with strainsoftening in the range ε =[0.01;1], and k er =1000 m 2 /yr: underthrusting of the 'Australian' mantle. b) Model with strainsoftening in the range ε =[[0.05;0.5] and k er =500 m 2 /yr: symmetric thickening of the mantle lithosphere. Because the lithospheric mantle has a small strength, amplitude and wavelength of Moho deflection are small.

Figure 9: Model with mantle lithosphere cohesion set at 200 MPa, with slower rate of strainsoftening (linear decrease from 200 MPa to 50 MPa within the range ε=0.12). Deformation propagates 'more' into the ductile lower crust, producing thrust faulting from X~90 to 130 km. As a result, the secondary maximum of topography reaches the same height as the primary one, after only 6 Ma of shortening.

Figure 4 :

 4 Figures captionsFigure1: a) Geological setting in the Southern Alps of New Zealand[START_REF] Mortimer | Basement geology of South Island of New Zealand[END_REF]: Eastern province terranes on the Pacific plate deform against Western province terranes. Dark grey lines show approximate location of geophysical transects for the SIGHT experiment. b) Seismic velocity model across SIGHT line 2 from[START_REF] Scherwath | Crustal structure from seismic data in the vicinity of the Alpine Fault[END_REF]. Note that horizontal coordinates do not coincide with the following model coordinates. Vertical exaggeration is 4.Figure2: Initial conditions for the 2d model (see text). Figure3: Surface deformation for reference model: topography (a, b, top) and erosion rates (c, d, bottom) plotted against horizontal distance (a,c, left) and time(b ,d, right). Symbols 'o' and 'x' correspond to the spatial maximum of topography after 3.2 and 7 Ma. After 7 Ma this maximum is located about 15 km east from the maximum erosion rate.Figure4: Reference model after 7 Ma of shortening (100 km) A.F. is Alpine Fault: a) surface topography, b) temperature field, c) deviatoric shear stress, d) total shear strain (elastic, viscous and brittle), e) shear strain rate and instantaneous velocity vectors, f) rheological phases (layers). Dark elements close to the surface are sediments. The vertical dashed line is positioned with reference to Figure (e) shear zone rising to the surface, namely, which we call the 'Alpine Fault zone'.Figure5: Shear strain ellipses after 7 Ma. a) Reference model: numerical values of the shear strain are same as for Figure(4d). b) Higher resolution model (400x40elements). Superimposition of the high electrical conductivity signal in thick grey contour: note the 'flip' of shear strain ellipses in the ductile lower crust, at X~25km. c) Most compressive principal stress lines for the reference model. Grey circles are earthquakes distribution, dashed lines are presumed geometry of faults at depth, thick grey line is the 12 km deep seismogenic zone (all data from[START_REF] Leitner | A focused look at the Alpine Fault, New Zealand, focal mechanisms and stress Observations[END_REF].Figure6: Comparison of surface convergence velocity with GPS measurements[START_REF] Beavan | Contemporary horizontal velocity and strainrate fields of the PacificAustralian plate boundary zone through New Zealand[END_REF]. Thick light grey line is the reference model, with distance from the plate boundary to the eastern border of the model Xl = 400km, and applied convergence rate Vx = 15mm/yr. Thin dark grey line is an alternative model with smaller distance Xl=200km and slower convergence rate Vx = 8 mm/yr. Note the difference in velocity increase at the Alpine Fault, and the similar change in slope at a distance of about 80 km to the East.Figure7: Comparisons of reference model with observed a) Topography, b) Moho, and c) gravity anomaly. Black lines: observed topography, Moho from seismic velocity models along Transect 2[START_REF] Scherwath | Crustal structure of the central South Island, New Zealand from an onshoreoffshore seismic data[END_REF], and measured Bouguer gravity anomaly. Grey lines: same horizons obtained with the reference model. In (a), vertical dashed lines are the observed (in black) and modelled (in grey) positions of the Alpine Fault (A.F.) and of the main divide (M.D.). 500 meters have been subtracted from the model topography (a), and 100 mGals have been added to the model gravity (c), (see text for justification).Figure8: Models with mantle lithosphere strength of 50 MPa, reducing to 20 MPa. Total shear strain (maximum in dark grey) and topography. a) Model with strainsoftening in the range ε =[0.01;1], and k er =1000 m 2 /yr: underthrusting of the 'Australian' mantle. b) Model with strainsoftening in the range ε =[[0.05;0.5] and k er =500 m 2 /yr: symmetric thickening of the mantle lithosphere. Because the lithospheric mantle has a small strength, amplitude and wavelength of Moho deflection are small.Figure9: Model with mantle lithosphere cohesion set at 200 MPa, with slower rate of strainsoftening (linear decrease from 200 MPa to 50 MPa within the range ε=0.12). Deformation propagates 'more' into the ductile lower crust, producing thrust faulting from X~90 to 130 km. As a result, the secondary maximum of topography reaches the same height as the primary one, after only 6 Ma of shortening.

  Figure1

Figure 2 Figure 3

 23 Figure 2

Figure 5 Figure 6 Figure

 56 Figure 5

Acknowledgements S. Cox, R. Sutherland, D. Woodward, D. EberhartPhillips, M. McSaveney, P.Koons, P. Upton, G. Caldwell, M. Scherwath provided extended information on the Southern Alps geology and geophysics. S. Ellis, D. Woodward, J. Beavan provided useful comments on the manuscript. We also thank R. Pysklywec, T. Stern and J. Braun for their thoughtful reviews. This study was funded by New Zealand Foundation for Research, Science and Technology (contract 05X0010) provided to the Institute of Geological & Nuclear Sciences (contribution 2222). Vrije Universiteit in Amsterdam has kindly provided facilities for M.G. to continue working on this project.