Catherine Dorbath

Muriel Gerbault

Gabriel Carlier

Michel Guiraud

The Double Seismic Zone of the Nazca Plate in Northern Chile: High Resolution Velocity Structure, Petrological Implications and ThermoMechanical Modelling

This paper presents an interdisciplinary study of the Northern Chile Double Seismic Zone. First, a high resolution velocity structure of the subducting Nazca Plate has been obtained by a new doubledifference tomography method. The double seismic zone (DSZ) is observed between 80 and 140 km depth and the two seismic planes are 20 km apart. Then, the chemical and petrologic characteristics of the oceanic lithosphere associated to this DSZ are deduced by using current thermalpetrologicalseismological models, and are compared to pressuretemperature conditions provided by a numerical thermo mechanical model. Our results agree with the common hypothesis that seismicity in both upper and lower planes are related to fluid releases associated to metamorphic dehydration reactions. In the seismic upper plane located within the upper crust, these reactions would affect material of basaltic (MORB) composition and document different metamorphic reactions occurring within highP (> 2.4 GPa) and low T (< 570°C) jadeitelawsonite blueschists and, at greater depth (> 130 km), lawsoniteamphibole eclogite conditions. The lower plane lying in the oceanic mantle can be associated to serpentinite dehydration reactions. The Vp and Vs characteristics of the region in between both planes are consistent with a partially (~2530 vol.% antigorite, ~010% vol. % brucite and ~410 vol. % chlorite) hydrated harzburgitic material. Discrepancies persist that we attribute to complexities inherent to heterogeneous structural compositions. While various geophysical indicators evidence particularly cold conditions in both the descending Nazca plate and the continental forearc, thermomechanical models indicate that both seismic planes delimit the innerslab compressional zone around the 400°C (±50°C) isotherm. Lower plane earthquakes are predicted to occur in the slabs flexural neutral plane, where fluids released from surrounding metamorphic reactions could accumulate and trigger seismicity. Fluids migrating upwards from the tensile zone below could be blocked in their ascension by the compressive zone above this plane, thus producing a sheeted layer of free fluids, or a serpentinized layer. Therefore earthquakes may present either downdip compression and downdip extension characteristics. Numerical tests indicate that innerslab compression is not only favored by the slab's thermal structure such as plate age. i) A weak ductile subduction channel, and ii) a cold mantle forearc both favor innerslab compression by facilitating transmission of compressional stresses from the continental lithosphere into the slab. iii) Decreasing the radius of curvature of the slab broadens the depth of innerslab compression, whereas iv) decreasing upper plate convergence diminishes its intensity. All these factors indicate that if indeed DSZs contour innerslab compression, they cannot only be linked to slab unbending, but also to the transmission of high compressional stresses from the upper plate into the slab.

Introduction

Since [START_REF] Hasegawa | Doubleplaned deep seismic zone and upper mantle structure in the northeastern Japan arc[END_REF] showed the existence of a doubleplaned deep seismic zone in the Northern Japan arc, more and more double seismic zones (DSZ) have been recognised in worldwide subduction zones like Tonga, New Britain, Kuriles, Kamchatka, Aleutians or Alaska [START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentinite dehydratation in subducting oceanic mantle?[END_REF]. Actually, a recent analysis of various subducting plates reveals the global prevalence of DSZ, which were found in segments of 16 subduction zones [START_REF] Brudzinski | Global Prevalence of Double Benioff Zones[END_REF]. For these zones, the intermediate depth earthquakes (80 to 180 km) define two dipping planes separated by a distance increasing with plate age, from ~8 km for a ~12Myold slab to ~30 km for a 160Myold slab. Seismicity in both upper and lower seismic zones occurs where stable hydrous minerals are predicted (Hacker et al., 2003b). Upper plane seismicity is linked to dehydration of metabasalts from the oceanic crust, and lower plane seismicity would be promoted by dehydration embrittlement of serpentinised mantle, e.g. antigorite (Hacker et al. 2003b;[START_REF] Yamasaki | Double seismic zone and dehydratation embrittlement of the subducting slab[END_REF][START_REF] Brudzinski | Global Prevalence of Double Benioff Zones[END_REF], with mantle hydration acquired by infiltration of seawater through outerrise faults prior to subduction [START_REF] Seno | Double seismic zones, compressional deep outerrise events, and superplumes, in Subduction: Top to Bottom[END_REF][START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentinite dehydratation in subducting oceanic mantle?[END_REF].

Despite this observation of a first order dependancy of DSZ separation on plate age by [START_REF] Brudzinski | Global Prevalence of Double Benioff Zones[END_REF], the mechanical cause for their occurence remains misunderstood. Numerous authors have suggested that DSZs are caused by unbending of the lithosphere after passing the position of maximum curvature (e.g. [START_REF] Engdahl | A double Benioff zone beneath the central Aleutians: An unbending of the lithosphere[END_REF]Isacks & Barazangi, 1977;[START_REF] Wang | Unbending combined with dehydration embrittlement as a cause for double and triple seismic zones[END_REF], inducing upper plane compressional earthquakes in the top part of the slab, and lower plane downdip tension earthquakes in its bottom part. But many slabs are nearly planar where the lower plane seismicity continues (e.g. [START_REF] Fujita | Double seismic zones and stresses of intermediate depth earthquakes[END_REF][START_REF] Abers | Relationship between shallowand intermediate depth seismicity in the Eastern Aleutian subduction zone[END_REF], suggesting that a continuous stress acts on the slab long after it has passed the point of unbending. Faultplane solutions are also found poorly consistent with plate bending and unbending about a neutral fiber (e.g. [START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentinite dehydratation in subducting oceanic mantle?[END_REF]. A while ago, [START_REF] Fujita | Double seismic zones and stresses of intermediate depth earthquakes[END_REF] had suggested that double seismic zones are observed in slabs for which the plate sinking rate (function of its age) differs from surface plate convergence velocity, sufficiently to generate stresssegmented seismic zones. When such slabs are on average, neither dominantly in tension or compression, seismogenic deviatoric stresses may be generated by local factors such as, for example, thermal stresses or compressional stresses transmitted from the upper plate into the slab.

It is known that a double seismic zone also exists in the subducting Nazca Plate beneath Northern Chile, just south of the Arica elbow. Yet, seismic data lead controversal DSZ characteristics. In this area, the 50 Ma old Nazca plate (e.g. [START_REF] Cande | Eocene propagating rifts in the southwest Pacific and their conjugate features on the Nazca plate[END_REF] is subducting beneath the South American plate at a rate of 84 mm/yr. From a microseismic field experiment between 20 and 24°S, [START_REF] Comte | An Inverted Double Seismic Zone in Chile: Evidence of Phase Transformation in the Subducted Slab[END_REF] observed, at depths of 80 to 150 km, two families of events with opposite focal mechanisms; the average separation between these two sets of events is about 15 km, suggesting a double layered slab. Nevertheless, the high quality results from the PISCO'94 experiment, deployed in the same region, gave no evidence for a double seismic zone [START_REF] Graeber | Threedimensional models of P wave velocity and PtoS velocity ratio in the southern central Andes by simultaneous inversion of local earthquake data[END_REF]. On the other hand, locally recorded events between 18.5 and 19.5°S display, at depths greater than 100 km, two distinct planes separated by 20 to 25 km, with an extreme variability of the focal mechanisms observed in both planes [START_REF] Comte | A DoubleLayered Seismic Zone in Arica, Northern Chile[END_REF]. About 300 km south, around 22°S, [START_REF] Rietbrock | A narrowky spaced doubleseismic zone in the subducting Nazca plate[END_REF] used the doubledifference algorithm (hypoDD, [START_REF] Waldhauser | hypoDD: A computer program to compute doubledifference hypocenter locations[END_REF] to evidence a 9 km thick DSZ at 80130 km depth, and with predominant extensional events in both upper and lower planes.

Recently, [START_REF] Zhang | Highresolution subductingslab structure beneath Northern Honshu, Japan, revealed by doubledifference tomography[END_REF] provided high resolution constraints on the petrological processes involved in double seismic zones by studying the subduction beneath Northern Honshu, using an improved doubledifference tomography method, tomoDD. Where previous seismic tomography studies lack information, particularly on the velocity structure of the deeper part of the slabs, this method is able to produce more accurate velocity images in the region near the sources [START_REF] Zhang | DoubleDifference Tomography: The Method and its Application to the Hayward Fault, California[END_REF]. In this paper, we first apply this new doubledifference tomography code to the Arica zone data (section 2.1), to obtain a detailed velocity model of the subducting Nazca plate where a double seismic zone has been observed (section 2.2). From this velocity model we then compare our results with those obtained for Honshu, and discuss the implications on petrological compositions when using [START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF] calculation scheme (section 2.2). Amongst these pecular implications, we argue for a cold thermal environment, and therefore develop in section 3 a thermomechanical numerical model of the Northern Chile subduction zone. We start this section 3 by reviewing the litterature, surface heatfluxes, thermal models and geophysical data of the area (section 3.1). We then expose our numerical approach, with initial conditions that are based on these data, and results that display the evolution of the subduction zone towards a transient thermomechanical state (section 3.2) that accounts for elasticviscous and brittle behaviour. We then reconsider our initial petrological results in comparison to the pressuretemperature conditions obtained with these thermomechanical models (section 3.3), and finally display alternate thermomechanical models, in which we show that not only temperatures but also stresses play a crucial rôle in the width and depth of DSZs (section 3.4). This paper is consequently rather dense as it contains a number of interdisciplinary considerations, but it is also this interdisciplinary approach that provided us with a global view of the processes in play, allowing to frame the pecular conclusions we propose.

2 High resolution velocity structure

Data and method

From June to August 1996, a dense temporary seismic network was deployed in the Arica region from the coast to the Bolivian border, between the Peruvian border and 19°20'S. This short period array of 24 vertical and 10 threecomponent stations complemented the Arica permanent network of 9 stations (Figure 1). P and S arrival times in these 43 stations were used to determine hypocentral locations. A description of the network and location process may be found in [START_REF] Comte | A DoubleLayered Seismic Zone in Arica, Northern Chile[END_REF]. In this paper, only the travel times of the very best located earthquakes recorded during this experiment, complemented by those recorded by the permanent network from January 1996 to December 1997, have been used.

We kept 2052 events, which met all the following restrictive criteria: number of readings higher than 10 including 2 S, distance to the nearest station lower than depth, rms lower than 0.25, azimuthal gap lower than 300 and conditioning factor lower than 100 (Figure 1). We prefered using these few very selected events instead of numerous redundant and poorly located events that will introduce noise in the data and degrade the results. We obtained about 36,000 absolute travel times (~55% of P waves), from which we constructed 171,000 differential travel times for events pairs. We used the same interevent maximal distance as [START_REF] Zhang | Highresolution subductingslab structure beneath Northern Honshu, Japan, revealed by doubledifference tomography[END_REF], e.g. 10 km at common stations, in order to reduce the velocity model dependancy on doubledifference earthquake relocation.

In addition to the hypocentral relocations, we simultaneously solve the P and S wave velocity structure at each grid node. The chosen inversion grid is presented on Figure 1 in map view. The origin of the grid is at X=69°W, Y=18°S, Z=0 km. Rotation of the reference frame to N77°E corresponds to the convergence direction of the Nazca plate underneath South America at that latitude. The distance between nodes is 10 km in the three directions, down to 160 km, below the maximal depth of the hypocenters in this area. We first inverted these traveltimes data for a 1D velocity model, starting from the very simple model previously used for the standard hypocentral location procedure in Northern Chile (in blue in Figure 2), and interpolating it in 17 regularly spaced layers defined between 0 and 160 km. The final model (in red in Figure 2) is very close to the model obtained in the same geographic area by [START_REF] Comte | Highresolution Subducting Slab Structure Beneath Northern Chile Using the DoubleDifference Method[END_REF] from a 1D inversion of 10,000 events recorded by 13 stations. In particular, the high velocity zone at ~50 km depth and the low velocity zone between 120 and 150 km are observed on both models. The obtained velocities were used as initial velocities for the 3D tomoDD inversion.

Numerous tests have been performed to check the robustness of the 3D solution. As we ran tomoDD in LSQR mode, we controlled the damping at each iteration by keeping the condition number value between 40 and 50 [START_REF] Waldhauser | hypoDD: A computer program to compute doubledifference hypocenter locations[END_REF]. Various interevent maximum distances have been tested, and show that inversions using different cutoffs give similar results. Several inversions were performed with various strategies of weighting absolute data relatively to differential data. Our final weighting scheme first a priori downweights the differential data, in order to obtain a largescale velocity model from the absolute data. In a second step, we favor the differential data, in order to improve the location of nearby events and refine the velocity model. In the last set of iterations, we then reweight the data, according to misfit and event separation, in order to remove or downweight outliers. Various distances between the grid nodes have been tested in order to get the best tradeoff between grid spacing and data resolution. The square root of the derivative weight sum (DWS) along profile P2 (Figure 1) is presented for Vp and Vs models in Appendix 1 (A) and (B) respectively. Thereafter, the isovalue 10, considered as the limit for resolved structures, is drawn on all figures. The influence of the chosen initial velocity model on the final velocity structures is presented in Appendix 1. Figures (C) and (D) present the results obtained for Vp and Vs starting from the final 1D model (red dots, Figure 2), whereas figures (E) and (F) present the results obtained starting from the initial model (blue dots, Figure 2). We see that the influence of the initial velocity model is minor.

Three resolution tests are presented in Appendix 2 for Vp, and Appendix 3 for Vs. In each test, the calculated synthetic absolute and differential traveltimes are included only for the same data distribution as for the real data. Following the standard procedure, we add uniformly distributed random noise to the synthetic data, in the interval +/ 0.025s for P waves and +/ 0.05s for S waves. The inverted models are obtained with the same inversion scheme as the real data set and starting from the same 1D final model. Similar resolution tests have been performed in DSZ contexts by [START_REF] Zhang | Highresolution subductingslab structure beneath Northern Honshu, Japan, revealed by doubledifference tomography[END_REF] or [START_REF] Shelly | Highresolution subduction zone seismicity and velocity structure beneath Ibaraki Prefecture, Japan[END_REF].

-Checkerboard test: the synthetic input P and Swave velocity models consist of 20 x 20 x 20 km patches (covering 2 nodes in horizontal and vertical directions) with velocities alternatively +5% and 5% compared to the normal starting model (Appendices 2 and 3, Figures A). Recovered P and S wave models are presented in figures B of Appendices 2 and 3. It is clear that P and S wave velocities features are well resolved in the subduction zone. The features are recovered in position and magnitude for the upper plane, but with a bit lower amplitudes for the lower plane. In the test presented here, the diagonal step of the grid (45°) differs from the dip of the slab (35°), and the colour alternance of diagonal high/low velocities, if similar to our results for Vp, is opposite to our results for Vs. The continental lithosphere area is poorly resolved down to 80 km to the east, a result that is evident in both following tests. -Restoring resolution test [START_REF] Zhao | Tomographic imaging of P and S wave velocity structure beneath northeastern Japan[END_REF] Except in the easternmost part of the recovered velocity models (F), no marked anomaly is observed. This third test shows that, if no anomaly exists, our data do not generate it, therefore if no low Vp/high Vs zone is present, the inversion does not generate it as an artifact.

Tomography results and discussion

We present on Figure 3 three cross sections through our final velocity models Vp and Vs and through the relocated hypocentres along profiles P1 (Y=20 km), P2 (Y=60 km) and P3 (Y=100 km) (Figure 1). These sections are representative of the slab in Northern Chile, as other cross sections show similar features. Figure 4 displays Vp and Vs zoomed on the DSZ along profile P2, which is in the center of the studied area and thus the best resolved. On both figures, the white lines correspond to the isovalue 10 of the square root of derivative weight sum, considered as the limit for resolved structures. The Vp/Vs model is obtained by simple division of Vp by Vs models in regions where both models are resolved. Appendix 1 (A) and (B) shows that the resolution of both models allows to calculate it. However, it is important to keep in mind that such direct division can carry instabilities.

At intermediate depth (80-140 km) the double seimic zone is situated under our seimic array, well defined and reasonably well resolved. Both seismically active planes are 5 to 10 km thick, with a clear separation of 20 km (Figure 3). In this part of the velocity model (Figure 4), the planes present outstanding distinct features which are summarized in Table 1.

The Forearc mantle:

In the forearc mantle below about 40 km and above ~110 km depths, Vp values are in the 7.58 km/s range, and Vs values are in between 4 and 4.5 km/s, very similar to the Honshu region. The Vp/Vs values do not present any large size anomaly, but oscillate from about 1.75 to 1.85. In contrast, in the Honshu region between 2090 km depth below the volcanic arc, Vp/Vs values are homogeneously high (~1.81.9), and are interpreted by [START_REF] Zhang | Highresolution subductingslab structure beneath Northern Honshu, Japan, revealed by doubledifference tomography[END_REF] to be associated with melt and water inclusions transported through hydraulic fractures from the slab where dehydration reactions occur. These dehydration reactions are attributed to the occurrence of serpentinized peridotite containing meltfilled pores, and are consistent with thermal studies of the Honshu region indicating temperatures above ~700°C. In Northern Chile, however, several studies indicate that the forearc mantle is significantly colder, such as thermal models by [START_REF] Honda | Thermal structure beneath Tohoku, Northeast Japan a case study for understanding the detailed thermal structure of the subduction zone[END_REF] or [START_REF] Springer | Interpretation of heatflow density in the Central Andes[END_REF], and flexural rigidity estimations by [START_REF] Tassara | Interaction between the Nazca and South American plates and formation of the Altiplano-Puna plateau: Review of a flexural analysis along the Andean margin (15°-34°S)[END_REF] and [START_REF] Tassara | Elastic thickness structure of South America estimated using wavelets and satellitederived gravity data[END_REF]. In addition, geophysical studies show an exceptionally low seismic attenuation (e.g. [START_REF] Haberland | Coincident anomalies of seismic attenuation and electrical resistivity beneath the southern Bolivian Altiplano plateau[END_REF], more details in the thermal structure section 3.2), and a specific seismic anisotropy pattern [START_REF] Kneller | Btype olivine fabric in the mantle wedge: Insights from highresolution nonNewtonian subduction zone models[END_REF][START_REF] Kneller | Stress, strain, and Btype olivine fabric in the forearc mantle: Sensitivity tests using highresolution steadystate subduction zone models[END_REF], details in section 3.2), also indicative of low temperature and high stress environment. Therefore, many lines of evidence show that the Northern Chile forearc mantle differs from other warm subduction zones that are characterised by extensively serpentinized, meltequilibrated forearc mantle (e.g. [START_REF] Kelemen | Thermal structure due to solidstate flow in the mantle wedge beneath arcs[END_REF][START_REF] Hyndman | Serpentinization of forearc mantle[END_REF], see also discussion in section 3.3).

The geological evolution of the south Peruvian and north Chilean forearc is complex as indicated by the occurrence of Precambrian blocks (i. e. Arequipa Massif) preserved within predominantly volcanic Phanerozoic material. Thus, the forearcmantle is probably also complex and predominantly composed of more or less metasomatized harzburgitic material, for which general compositions systematically lead to high Vp (> 8. km/s) and Vs (> 4.6 km/s) values. At depths of 40110 km, our recorded velocity values (Vp = 7.58.0 km/s ; Vs = 4.04.5 km/s) militate for a partially hydrated harzburgitic forearc mantle, when deduced from petrological simulations (1.22.9% H 2 O, 1024% modal antigorite; Hacker et al. 2003a;[START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF]. These petrological simulations are also consistent with the presence of a cold mantle in Northern Chile (<600°C).

The DSZ, petrological candidates deduced from tomographic results

Our results (Table 1) are very close to the slab velocity structure obtained by [START_REF] Zhang | Highresolution subductingslab structure beneath Northern Honshu, Japan, revealed by doubledifference tomography[END_REF] in the Honshu DSZ at 39°N, and slightly different from those obtained by [START_REF] Shelly | Highresolution subduction zone seismicity and velocity structure beneath Ibaraki Prefecture, Japan[END_REF] at 37°N. In Honshu, the old (130 Ma) Pacific plate subducts beneath the Eurasian plate at 91 mm/yr [START_REF] Demets | Current plate motions[END_REF]. The two seismic planes of the DSZ are separated by ~30 km at intermediate depths of 50-140 km.

In New Zealand, the observed DSZ are not as clear as in Japan, however the velocity structure obtained by [START_REF] Zhang | Double difference seismic tomography method ans its application[END_REF] with tomoDD in the Wellington region, may be compared to the Chilean results. Nevertheless in this region, the Pacific plate is about 120 Ma old, and the subduction process is much slower (42 mm/yr). Both seismic planes are detected at shallower depths, from 50 km to only 80 km due to landsea limits, and are separated by ~10 km. The results of these three tomoDD studies are summarized in Table 2. The structure of the Hikurangi subduction zone beneath the central North Island, New Zealand, to the north of Zhang's study [START_REF] Reyners | Imaging subduction from the trench to 300 km depth beneath the central North Island, New Zealand, with Vp and Vp/Vs[END_REF] was also obtained by simultaneous inversion for both hypocentres and the 3D Vp andVp/Vs structure [START_REF] Thurber | Earthquake locations and threedimensional crustal structure in the Coyote Lake area, central California[END_REF][START_REF] Thurber | Local earthquake tomography: velocities and Vp/Vstheory[END_REF]. Results define the DSZ down to 130-140 km depth. Earthquakes in the lower plane are generally associated with regionally high Vp values, and a specifically low Vp/Vs, which is similar to lower plane earthquakes in Northern Honshu.

Focal mechanisms would help identifying the physical processes and metamorphic reactions involved in Northern Chile DSZ. Focal mechanisms constructed for the events recorded during the 1996 campaign have been presented by [START_REF] Comte | A DoubleLayered Seismic Zone in Arica, Northern Chile[END_REF]. For the present study, we used our relocated events to construct focal mechanisms when possible. Only 45 solutions were constrained, among which 41 located in the upper plane, and 4 remaining in the lower plane. Our results do not differ from those obtained by [START_REF] Comte | A DoubleLayered Seismic Zone in Arica, Northern Chile[END_REF]: fault plane solutions for these intermediate depth events vary significantly, even between nearby events. The very few events associated to the lower subduction plane also present an extreme variability (with both tensional and compressional mechanisms). Finally, polarities read from the very few earthquakes located in between both planes do not enable to draw any reliable mechanism. Such results were also observed by [START_REF] Abers | Relationship between shallowand intermediate depth seismicity in the Eastern Aleutian subduction zone[END_REF] in the Aleutian DSZ.

In order to assess the petrological compositions consistent with our tomographic data, the velocities proposed and discussed hereafter have been computed by averaging the values obtained in the central part of our model for the same area of the DSZ (upper plane, lower plane, and region in between), from P1 to P3 (Figure 1). Mineral rock compositions involved in upper and lower planes, and in between, are obtained by using Hacker & Abers's procedure (2004). For suitable Vp velocity data, we have selected rock compositions displaying the Vs velocities that most closely approximate the observed data. Then, the results are compared with available, independent petrogenetic data.

2.2.3The upper plane

The upper plane of earthquakes in Northern Chile lies in a region of relatively low Vp (7.7 km/s) and normal Vs (4.6 km/s) values, thus relatively low Vp/Vs ratio (~1.67) (Table 1), from 80 to 130 km depth. Between 80 km (~2.4 GPa) and 140 km (~4.5 GPa) depth and, for rocks of basic composition (MORB), mineral velocities estimated by using the calculation scheme of [START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF] imply mineral assemblages of relatively low Vp and high Vs values and, relatively high Vp and low Vs values, respectively for each depth range. Common minerals at such pressure and temperature are coesite, jadeite, kyanite and garnet, that show high Vs (>4.6 km/s) but also high Vp (>7.8 km/s) values. Therefore, whatever their respective modal proportions, these minerals must necessarily be associated to amphibole and/or chlorite, which significantly decrease the bulkrock Vp values. Conversely, the observed low Vp velocity values preclude anhydrous metabasic rocks such as eclogite. Velocity simulations with [START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF] tables agree with this conclusion.

Petrogenetic grids modelling high pressure metabasic rocks (CNFMASH system + silica phase and H 2 O in excess) provide additional constraints (e.g., [START_REF] Clarke | Eclogites and Blueschists of the Pam Peninsula, NE New Caledonia: a reappraisal[END_REF][START_REF] Carson | Calculated mineral equilibria for eclogites in CaONa 2 OFeOMgOAl 2 O 3 SiO 2 H 2 O: application to the Pouébo Terrane, Pam Peninsula, New Caledonia[END_REF]. In theses systems, at temperatures above 550°C, garnet and clinozoisite occur at the expense of chlorite and lawsonite, respectively. In contrast at low temperatures (<570°C), chlorite coexists with lawsonite, amphiboles and omphacite over a large pressure domain (e.g., lawsonite is stable above 1.6 GPa at 550°C). Since the relatively low Vp values of Northern Chile upper seismic plane imply the presence of chlorite and amphibole rather than garnet and clinozoisite, we propose that the involved mafic metamorphic rocks are characterized by low temperature (<570°C) mineral assemblages in presence of water. Petrological simulations using low temperature (500°C) and lawsoniteomphaciteamphibole chlorite assemblages with 4.05.7 wt % H 2 O and, at greater depth, lawsonite amphibole eclogite with 3.0 wt % H 2 O [START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF], provide Vp values that are close to observations. Accordingly, we suggest that between 80 and 140 km depths, the seismic upper plane in Northern Chile is related to fluid release associated to metamorphic reactions occurring within jadeitelawsonite blueschists and, at greater depth, under lawsoniteamphibole eclogite conditions. Note that these conditions occur at relatively low temperature (<570°C), consistent with indications of the cold state of the North Chilean subduction zone (see section 3.1).

The intermediate zone

In between both seismic planes of the Northern Chile DSZ, mean Vp values are high (8.5 km/s), mean Vs values are slightly low to normal (4.5 km/s), leading a high Vp/Vs ratio (~1.8) (Table 1). Nevertheless, we observe a global decrease of the velocity values with depth, from [Vp ~8.8, Vs ~4.7] at ~80 km depth, to [Vp ~8.3,Vs ~4.4] at ~130 km depth. Vp values as high as 8.8 km/s are indicative of unmetamorphosed or metamorphosed harzburgite composition (Hacker et al., 2003a) under very high pressure (above 4.5 GPa for unmetamorphosed harzburgite, above 7 GPa for metamorphosed harzburgite), much greater than those inferred from the lithostatic conditions of the intermediate zone of Northern Chile(< 3.5 GPa). Simulations using [START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF] indicate that, at 80 km depth and temperature <500°C, observed mean Vp values are consistent with anhydrous or slightly hydrated harzburgitic material (<0.6 wt% H 2 O, ~10 vol. % antigorite). At 130 km depth, observed Vp values of 8.3 km/s are consistent with a similar material of water content up to 1.2 wt% (~10 vol. % antigorite) at low temperature (300°C) and 0.6 wt% (~5 vol. % antigorite) at higher temperature (600°C).

2.2.5The lower plane

In Northern Chile, the lower plane of seismicity is characterized at all depths down to 130 km, by low Vp values (7.4 km/s), high Vs values (4.7 km/s), and very low Vp/Vs ratio (~1.57) (Table 1). Similar data have only been reported in the lower plane of Wellington DSZ, New Zealand [START_REF] Zhang | Double difference seismic tomography method ans its application[END_REF]; Table 2). Such results have been generally attributed to serpentinites or chloritebearing harzburgites. However, at depths greater than 80 km (above ~2.5 GPa), occurrence of fully serpentinized rocks (15 wt% H 2 O) is ruled out because they display too low Vp values (<6.2 km/s) with respect to the observed data of Northern Chile. Better results are obtained when partially hydrated harzburgite compositions are considered (57 wt% H 2 O, with ~2530 vol.% antigorite, ~010% vol. % brucite and ~410 vol. % chlorite, from petrological simulations,). At considered depths and in order to maintain relatively high Vs, such compositions imply that temperature probably does not exceed 500 °C in brucitebearing or, depending on water activity, 500650°C in brucitefree assemblages (see, i. e., Hacker et al., 2003a;[START_REF] Perrillat | Kinetics of antigorite dehydration: a realtime Xray diffraction study[END_REF].

In conclusion, like for other worldwide DSZs, Northern Chile lower plane seismicity can be ascribed to dehydration reactions occuring within partially hydrated harzburgites. As noted above, the most probable involvement of brucite breakdown in these dehydration reactions and the fact that water activity is probably <1, implies that temperatures do not exceed 600°C.

Serious limitations to previous interpretations arise when considering the observed Vs values and Vp/Vs ratios (Table 1). Even worse, no mineral assemblage displaying both observed Vp and Vs values can be simulated using Hacker & Abers's calculation tables. The simulations systematically result in mineral assemblages displaying higher or lower Vs velocity for a suitable Vp velocity and vice versa. Such discrepancy may be ascribed to various parameters such as the presence of more complex oceanic crust or mantle compositions, occurrence of freefluid phases (i.e., [START_REF] Nakajima | Threedimensional structure of Vp, Vs, and Vp/Vs beneath northeastern Japan: Implications for arc mgmatism and fluids[END_REF]Hacker et al., 2003a), pore geometry (i. e., [START_REF] Takei | Effect of pore geometry on Vp/Vs: From equilibrium geometry to crack[END_REF], reaction kinetics [START_REF] Evans | The serpentinite multisystem revisited: chrysotile is metastable[END_REF][START_REF] Perrillat | Kinetics of antigorite dehydration: a realtime Xray diffraction study[END_REF], suitable mixing models used to convert velocity data in mineral composition (i. e., [START_REF] Wang | Elastic properties of forsteriteenstatite composites up to 3.0 GPa[END_REF], or medium to large scale rock heterogeneity (lattice preferred orientations, compositional layering; Ji & Wang, 1999).

Contribution of metamorphic sediments and/or greywackes within the subducting upper crust, as proposed for the Wellington DSZ, New Zealand [START_REF] Zhang | Double difference seismic tomography method ans its application[END_REF], could drastically reduce Vp/Vs ratios [START_REF] Christensen | Poisson's ratio and crustal seismology[END_REF]. However, such a hypothesis is unlikely in Northern Chile because the subducting Nazca plate is known to have experienced very few pelagic deposits (< 850 m; [START_REF] Prince | Crustal rupture and the initiation of imbricate thrusting in the PeruChile Trench[END_REF]. This feature is a direct consequence of the arid climate affecting the southwest intertropical coast of South America (e.g. Lamb & Davis, 2003, and references therein).

Petrological diversity in mantle material is rather limited. Nevertheless, addition of magnetite and clinopyroxene (diopside +jadeite) tend to increase the Vs values of the partially hydrated harzbugitic material used in our similations. Occurrence of magnetite is common in serpentinized rocks where it reaches up to 10 vol. % of the rock (i. e., [START_REF] Coleman | A chemical study of serpentinisation, Burro Mountain, California[END_REF]. Significant amounts of clinopyroxene (> 8 vol. %) suggests that lherzolitic or wehrlitic layers could be interstratified within the harzburgitic material, a feature observed in some ophiolitic complexes that are considered to be typical of oceanic lithosphere (i. e., Oman ophiolites; [START_REF] Ceuleneer | Mantle flow patterns at an oceanic spreading centre : the Oman peridotites record[END_REF][START_REF] Benn | Mantlecrust transition zone and origin of wehrlitic magmas : evidence from the Oman Ophiolite[END_REF][START_REF] Takazawa | Geochemistry and origin of the basal lherzolites from the northern Oman Ophiolite (northern Fizh Block)[END_REF]. This presence of clinopyroxenerich layers in harzburgitic materiel brings up an another problem generally not considered in largescale models (i. e., Hacker et al., 2003a, b;[START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF]) : rock heterogeneity. The main characteristic of natural outcrops of mafic or ultramafic materials is their textural heterogeneity (i. e., Ceulener et al., 1988;[START_REF] Benn | Mantlecrust transition zone and origin of wehrlitic magmas : evidence from the Oman Ophiolite[END_REF][START_REF] Clarke | Eclogites and Blueschists of the Pam Peninsula, NE New Caledonia: a reappraisal[END_REF][START_REF] Carson | Calculated mineral equilibria for eclogites in CaONa 2 OFeOMgOAl 2 O 3 SiO 2 H 2 O: application to the Pouébo Terrane, Pam Peninsula, New Caledonia[END_REF][START_REF] Takazawa | Geochemistry and origin of the basal lherzolites from the northern Oman Ophiolite (northern Fizh Block)[END_REF]. For example, jadeitelawsonite blueschists and eclogites commonly exhibit metric, decametric or kilometric eclogitic, gabbroic lenses or layers, and peridotites generally contain pyroxenitic layers or dykes. In addition, some layers are more intensively deformed whereas others are less deformed (i. e., gabbroic or eclogitic lenses in blueschists, serpentinite layers in peridotites). Such features introduce significant anisotropy that could have important consequences on measured seismic velocities. In particular, it is remarkable that the Northern Chile DSZ lower plane displays Vp velocity values indicative of partly serpentinized peridotite (7.4 km/s) but Vs velocity values typical of anhydrous peridotite (4.7 km/s). The possibility that these data represent a fully serpentinized matrix containing anhydrous harzburgite lenses or layers should be explored. Nevertheless, such a hypothetical concept needs to be experimentally tested.

Many authors (i.e., [START_REF] Nakajima | Threedimensional structure of Vp, Vs, and Vp/Vs beneath northeastern Japan: Implications for arc mgmatism and fluids[END_REF][START_REF] Takei | Effect of pore geometry on Vp/Vs: From equilibrium geometry to crack[END_REF]Hacker et al., 2003a,b) suggest that the presence of free H 2 Orich fluids phases can explain low Vp/Vs ratios. The occurrence of water in excess is certainly required by the modelling of amphibolerich eclogitic rocks in [START_REF] Carson | Calculated mineral equilibria for eclogites in CaONa 2 OFeOMgOAl 2 O 3 SiO 2 H 2 O: application to the Pouébo Terrane, Pam Peninsula, New Caledonia[END_REF]. Hydrous minerals, mainly chlorite, hornblende, epidote or lawsonite for mafic compositions, and antigorite, chrysotile, talclike phase, or chlorite for ultramafc compositions are the source of the water released during prograde metamorphism. Discontinuous dehydration reactions and/or entrapment of water by physical processes during continuous dehydration could lead to the observed low Vp/Vs, and hence favor seismicity. Kinetics is an important factor that can affect the modeling. Petrogenetic grids are based on equilibrium thermodynamics. Yet, as the slab moves downwards, rapid changes in P and T can lead to kinetics driving reactions. Experiments by [START_REF] Perrillat | Kinetics of antigorite dehydration: a realtime Xray diffraction study[END_REF] suggest that antigorite could breakdown at temperatures much below equilibrium petrogenetic grids (e.g. 520°C at 4GPa), thus consistent with release of water at low temperature, especially if fluid is taken out from the medium at the same rate as it is produced. [START_REF] Jung | Intermediatedepth earthquake faulting by dehydration embrittlement with negative volume change[END_REF] have also shown that antigorite breakdown under shear stress does not follow equilibrium thermodynamics, and that dehydration embrittlement could be the process for continuous dehydration.

Another parameter that can significantly affect the rock Vp/Vs ratios is the pore geometry [START_REF] Watanabe | Effects of water and melt on seismic velocities and their application to characterization of seismic reflectors[END_REF][START_REF] Takei | Acoustic properties of partially molten media studied on a simple binary system with a controllable dihedral angle[END_REF][START_REF] Takei | Effect of pore geometry on Vp/Vs: From equilibrium geometry to crack[END_REF]. For example, a specific pore configuration (i.e., tube geometry, equilibrium geometry, Watanabe 1993) corresponding to low hydration and waterundersaturated pores could be responsible for the observed velocities in the intermediate zone of the DSZ.

All our previous petrological interpretations of seismic velocities are based on data provided by Hacker et al. (2003a) and [START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF], for which oceanic lithosphere materials are assumed restricted to an average MORB composition for the crust and an average harzburgite composition for the mantle. In order to explain our measured seismic velocities, we had to assume very specific compositional characteristics, and especially very cold conditions. Obviously other parameters such as rock heterogeneity and fluid configurations still require additional investigations in order to accurately evaluate their relative contribution. Meanwhile, an attempt to provide additional constraints is made in the following section, first by presenting the thermal and geophysical characteristics of both the subducting and the overriding plates, and then by developing thermomechanical models of the Northern Chile subduction zone.

3Thermomechanical modelling

The idea is to setup a model with "reasonable" initial thermal and rheological structure, close to the Northern Chile conditions, and to record stress buildup and propagation of deformation when applying upper plate convergence and basal slabpull, over a relatively shorttime scale of 2 Ma. We present first the available data on the thermal structure of the region (section 3.1), referring to surface heatflow measurements, existing thermal studies, and geophysical data. We then describe our numerical thermomechanical approach (section 3.2), which differs from standard thermal models of subduction zones. A reference model is then presented, illustrating the general dynamics of deformation, resulting from specific initial conditions. With this model, we discuss the PressureTemperature specificities for the Northern Chile DSZ (section 3.3). Finally, we illustrate with alternate models how the oceanic and continental thermal structures, the subducting channel rheology, and relative plate convergence, all affect these PressureTemperature conditions (section 3.4).

Thermal structure in Northern Chile inferred from data

Considering existing thermal studies of the Northern Chile subduction zone, we first recall the heat flow data, studied by [START_REF] Springer | Heatflow density across the Central Andean subduction zone[END_REF], who combined existing values (e.g. [START_REF] Hamza | Heat flow map of South America[END_REF] with boreholes data. The mean heatflow density along a WE generalized lithospheric structure shows the following trends: (1) Within the oceanic Nazca Plate low values of about 30 mW/m 2 occur in the region of the PeruChile country boundary: these are low compared to the 70mW/m 2 deduced from a 50 Ma old plate cooling model. [START_REF] Springer | Heatflow density across the Central Andean subduction zone[END_REF] explained this value by a probable faster cooling resulting from heat convection compared to pure heat conduction models. (2) A minimum heat flow density of about 20 mW/m 2 is observed in the Coastal Cordillera, which would reflect cooling in the continental crust as a result of the subduction of oceanic lithosphere. Values increase toward the forearc region to 4060 mW/m 2 . (3) Heat flow varies from about 50 to 180 mW/m 2 in the area of the magmatic arc (sparse and equivocal) and the Altiplano. This large variability is interpreted by [START_REF] Springer | Heatflow density across the Central Andean subduction zone[END_REF] as heat sources resulting from isolated magma chambers at shallow depths (~10km), in agreement with geological and geophysical information (see review in [START_REF] Schilling | Partial Melting in the Central Andean Crust: a Review of Geophysical, Petrophysical, and Petrologic Evidence[END_REF]. ( 4) High values of about 80 mW/m 2 are typical for the backarc region of the Eastern Cordillera, attributed to either doubling of the crust or rise of the asthenosphere. ( 5) Heat flow density is about 40 mW/m 2 in the Subandean Ranges and Chaco Basin, somewhat lower than the 5060 mW/m 2 of the adjacent Brazilian Shield farther to the east [START_REF] Hamza | Heat flow map of South America[END_REF].

Geophysical investigations also give insight on the thermal structure of Northern Chile. We summarize here a review by [START_REF] Schilling | Partial Melting in the Central Andean Crust: a Review of Geophysical, Petrophysical, and Petrologic Evidence[END_REF]. From magnetotelluric investigations, the volcanic arc does not appear as a conductive feature along the three MT profiles carried north of 21°S, but show instead a low conductivity zone, and very low conductivity in the forearc beneath the Coastal Cordillera, extending into the upper mantle (e.g. [START_REF] Brasse | The Bolivian Altiplano conductivity anomaly[END_REF]. However, there is a broad and intense anomaly in the midcrust beneath the plateau at least 4050 km deep (note that this is more than 350 km east of the trench), with evidences for a conductive root reaching the upper mantle, possibly representing the mantle wedge [START_REF] Brasse | The mantle wedge in the Bolivian orocline in the view of deep electromagnetic soundings[END_REF]. Seismic waves through the forearc crust indicate a cold crust containing only minor free fluids, and Qp tomography shows no indication of a fluid "curtain" extending from the descending slab into the forearc region. Some regions of the forearc show extremely low attenuation, indicative of a strong, cold lithospheric block (e.g. [START_REF] Haberland | Coincident anomalies of seismic attenuation and electrical resistivity beneath the southern Bolivian Altiplano plateau[END_REF][START_REF] Schurr | Deep seismic structure of the Atacama basin, Northern Chile[END_REF]. This may impede magma generation, and explain the observed shift of the volcanic arc to a depth to the top of the slab reaching 110130 km (e.g. global study by Syracuse & Abers, 2007).

The thermal structure of our numerical model should thus be constructed to comply with these data, and will be mostly based on previous thermal models developped by [START_REF] Springer | Interpretation of heatflow density in the Central Andes[END_REF], and gravity models by [START_REF] Tassara | Threedimensional density model of the Nazca plate and the Andean continental margin[END_REF]. This setup will be described in detail further in section 3.2.2, after having first presented our general numerical approach.

The thermomechanical modelling approach

We do not wish to account for full thermodynamical modeling of subduction processes in Northern Chile, for various reasons: 1) this subduction zone has been ongoing for almost 80 millions of years with still debated alternances of extensional and compressional upper plate arc regimes. For example, the very characteristic Altiplano feature that has built on the upper plate obviously indicates specific deep structural properties, and its rise, linked with the formation of the Arica elbow, is still debated to have occurred since either within the last 10 Myrs or 35 Myrs ago (e.g. Roperch et al., 2006, and references therein). 2) full thermomechanical models have already been published by [START_REF] Sobolev | What drives orogeny in the Andes?[END_REF], 2006), over a time duration of about 20 Myrs. Since in fact, a number of geophysical data and conceptual models of the present day structure of the Northern Chile margin exist, instead in this study, we wish to reproduce on a relatively shorttime scale of about 2 Myrs, specific deformation and stress patterns surrounding the subduction zone interface. Therefore, we start from an hypothetical thermal structure, a geometrical structure and a rheological layering, that are built as close as we can to available published information, and which, implicitely, assume a relative thermal steady state and dynamical equilibrium of the forces in play.

We use a Lagrangian approach that simulates elasticviscousbrittle behavior, Parovoz [START_REF] Poliakov | Diapirism and topography[END_REF]), a 2D finitedifferences numerical code based on the FLAC technique [START_REF] Cundall | A microcomputer program for modeling largestrain plasticity problems[END_REF]. It resolves differential equations alternately, with the output for the solution of the equations of motion used as input to the constitutive equations for a progressive, timeexplicit, calculation. FLAC is based on a quadrilateral mesh formed by superimposing two constant stress/constant strain triangular finite elements. Each element behaves according to prescribed stress/strain constitutive laws in response to kinematical boundary conditions: MohrCoulomb nonassociative elastoplasticity is combined with temperaturedependent viscoelasticity (Maxwell) so that locally, the minimum of both stresses is chosen. Accounting for elasticbrittleductile properties permits to model formation of faults or shear zones and lithospheric extension and collision (e.g. [START_REF] Poliakov | Fractal plastic shear bands[END_REF][START_REF] Lavier | Selfconsistent rollinghinge model for the evolution of largeonset lowangle normal faults[END_REF][START_REF] Burov | Erosion and rheology controls on syn and postrift evolution: verifying old and new ideas using a fully coupled numerical model[END_REF][START_REF] Gerbault | Numerical models of lithospheric deformation forming the Southern Alps of New Zealand[END_REF]. The free upper surface boundary allows to adequately resolve for surface processes, which are accounted for with a diffusion equation with variable coefficient (Table 3). The basic equations of motion and heat transport in Lagrangian formulation for Parovoz can be found in a number of publications (e.g. [START_REF] Burov | Ascent and emplacement of magma reservoirs in brittleductile upper crust[END_REF][START_REF] Gerbault | Numerical models of lithospheric deformation forming the Southern Alps of New Zealand[END_REF][START_REF] Pourhiet | Rifting through a stack of inhomogeneous thrusts (the dipping pie concept)[END_REF][START_REF] Yamato | Burial and exhumation in a subduction wedge: Mutual constraints from thermomechanical modeling and natural PTt data (Schistes Lustrés, western Alps)[END_REF]. The heat transport equation accounts for heat conduction, radioactive and shear heating, while heat advection is accounted for by the displacements of the Lagrangian grid. This version of Parovoz does not account for phases reactions and subsequent density evolution, and employs a remeshing technique using passive markers similar to [START_REF] Yamato | Burial and exhumation in a subduction wedge: Mutual constraints from thermomechanical modeling and natural PTt data (Schistes Lustrés, western Alps)[END_REF]. See appendix 4 for more details.

Rheologies and boundary conditions

The model (Figure 5) simulates an area 300 km deep by 1800 km wide, comprising a flexuring oceanic plate dipping at 35°. With 900 x 150 grid elements, horizontal resolution is 2 km, and vertical resolution decreases downwards from 1 to 4 km. Geometrical, rheological and thermal properties are summarized in Tables 3, 4, 5. Elastic, brittle and powerlaw creep parameters are chosen according to various sources from rock mechanics extrapolated to the lithospheric scale (e.g. [START_REF] Kirby | Intermediatedepth intraslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs[END_REF]Rannalli, 1995, and references therein), and combined with the known special conditions for the Andes (e.g. [START_REF] Springer | Interpretation of heatflow density in the Central Andes[END_REF][START_REF] Tassara | Interaction between the Nazca and South American plates and formation of the Altiplano-Puna plateau: Review of a flexural analysis along the Andean margin (15°-34°S)[END_REF][START_REF] Sobolev | What drives orogeny in the Andes?[END_REF]. Composition and densities are defined mostly according to the recent threedimensional density model from [START_REF] Tassara | Threedimensional density model of the Nazca plate and the Andean continental margin[END_REF], who compiled appropriate characteristics of the Chilean Andes.

The following heterogeneities are inserted (see Tables 4,5). First, in the continental crust (various red to brown colours in Figure 5): the upper crust is supposed to be of weak composition in the Bolivian foreland (e.g. [START_REF] Sobolev | What drives orogeny in the Andes?[END_REF], as opposed to the forearc crust, which is of rather mafic composition, coherent with geophysical data and geological history (e.g. flexural rigidity study by [START_REF] Tassara | Interaction between the Nazca and South American plates and formation of the Altiplano-Puna plateau: Review of a flexural analysis along the Andean margin (15°-34°S)[END_REF][START_REF] Tassara | Factors controlling the crustal density structure underneath active continental margins with implications for their evolution[END_REF]; seismology and magnetotellurics reviewed for example by [START_REF] Schilling | Partial Melting in the Central Andean Crust: a Review of Geophysical, Petrophysical, and Petrologic Evidence[END_REF] geology indicating longlived arc magmatism, [START_REF] Lucassen | Composition and density model of the continental crust in an active continental margin -the Central Andes between 18° and 27[END_REF]. The Altiplano domain is characterised by a prethickened lower crust (65 km with respect to 35 km elsewhere), which has a more felsic behaviour (e.g. [START_REF] Lucassen | Composition and density model of the continental crust in an active continental margin -the Central Andes between 18° and 27[END_REF][START_REF] Beck | The nature of orogenic crust in the central Andes[END_REF].

Second, the subduction channel, defined 10 km thick, is composed of three portions: the upper 10 km mimick unconsolidated oceanic crust, with a plagioclase rheology and a moderate friction angle of 5°. Below, weaker wet granite rheology is employed (this choice is discussed later on) and friction is 0.025 (angle of 1.43°), corresponding to the mean value proposed by [START_REF] Lamb | Shear stresses on megathrusts: implications for mountain building behind subduction zones[END_REF]. At 120 km the density of the channel increases, mimicking eclogitization. The driving boundary conditions are 1) a westward velocity Vr= 3 cm/yr applied from the eastern edge of the model lithosphere, simulating South America plate motion, and 2) a velocity traction Vb = 5 cm/yr inclined at 35°downwards exactly where the oceanic plate intersects the base of the model at 300 km depth (similar to [START_REF] Sobolev | What drives orogeny in the Andes?[END_REF]. Lithostatic normal stresses are applied at the base and lateral borders of the model at asthenospheric depths, and along the extremity of the oceanic plate along the left border.

Mechanically, we start from a geometrical configuration of layers of different density, which are compensated only in the vertical direction with the appropriate initial topography. About 1 millions years are needed for the model to build up stresses and deform internally in response to applied assymetric slab pull and upper plate driving forces. This timescale is defined by elastic and yield parameters and the loading strainrate, and can be roughly estimated from 1D Hooke's law: with a maximum yield stress close to 1 GPa, a standard Young's modulous equal to 30 GPa, and ~7 cm/yr of horizontal compression applied over a model length of 1800 km.

-Thermal structure

Thermally, we proceed similarly as in the mechanical approach, in the sense that we presume of initial temperatures close to steady state, that will not evolve significantly during the 2 Myr of the model duration: heat conduction, radioactive heating and shear heating, as well as convective heat transport are all directly limited by the short duration of the run. This short 2 Myr time scale also allows us to neglect evolving phase transformations (eclogitization or serpentinization of the oceanic crust), partial melting, or thermochemical effects. The initial thermal structure is thus based on previous estimates by [START_REF] Springer | Interpretation of heatflow density in the Central Andes[END_REF] and [START_REF] Tassara | Threedimensional density model of the Nazca plate and the Andean continental margin[END_REF]. We first calculate oceanic and continental geotherms as a function of their thermal age [START_REF] Parsons | An analysis of the variation of ocean floor bathymetry and heat flow with age[END_REF][START_REF] Turcotte | Geodynamics: applications of continuum physics to geological problems[END_REF][START_REF] Burov | The effective elastic thickness of continental lithosphere: what does it really mean?[END_REF]. We choose an age of 60 Myrs for the oceanic lithosphere, in order to deal with a lithosphere, only slightly but reasonably colder than the unsatisfying "theoretical" 50 Mrsold age. We thus have surface values of 60 mW/m 2 , still greater than the measured 3040mW/m 2 . Implications will be discussed later on. For continental temperatures, we define the thermal depth of the continental plate at 180 km [START_REF] Tassara | Threedimensional density model of the Nazca plate and the Andean continental margin[END_REF], standard radiogenic crustal heating (see table 3), and insert a gaussian shape thermal anomaly initiating from the interplate contact, extending below the Altiplano thickened crust to the Subandean area, that superimposes on [START_REF] Tassara | Threedimensional density model of the Nazca plate and the Andean continental margin[END_REF] estimation of the lithosphereasthenosphere boundary. [START_REF] Tassara | Threedimensional density model of the Nazca plate and the Andean continental margin[END_REF] obtained a lithosphereastenosphere boundary (isotherm 1350°C) moving up to 60 km depth below the Western Cordillera and Altiplano, descending down to 160180 km depth further East (Figure 5). At the interplate contact, Springer (1999) used a heating shear stress of 15 Mpa, whereas [START_REF] Lamb | Shear stresses on megathrusts: implications for mountain building behind subduction zones[END_REF] employs 10 MPa, from seismic coupling arguments: we insert an interface shear heating stress of 10 MPa.

In fact, we will not follow here the standard approach for thermal modelling of subduction zones, in which a kinematically prescribed slab drives a dynamic cornerflow in the wedge under a stationary overriding plate (e. g. [START_REF] Peacock | Seismic Consequences of Warm Versus Cool Subduction Metamorphism: Examples from Southwest and Northeast Japan[END_REF][START_REF] Van Keken | Highresolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle[END_REF][START_REF] Kelemen | Thermal structure due to solidstate flow in the mantle wedge beneath arcs[END_REF][START_REF] Conder | A case for hot slab surface temperatures in numerical viscous flow models of subduction zones with improved fault zone parameterization[END_REF][START_REF] Kneller | Btype olivine fabric in the mantle wedge: Insights from highresolution nonNewtonian subduction zone models[END_REF][START_REF] Kneller | Stress, strain, and Btype olivine fabric in the forearc mantle: Sensitivity tests using highresolution steadystate subduction zone models[END_REF][START_REF] Abers | The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow[END_REF]. These thermal models use an Eulerian method that resolves the temperature equation for a steady state that is reached over several tens of millions of years (ideally 100 Myrs), and allows to study the thermal influence of viscous mantle corner flow on the surrounding plates. Our approach here differs since we only aim at reproducing 2 Ma of thermo mechanical dynamics, and therefore assume that the thermal state at the onset of the run is equivalent (e.g. does not change significantly) to that at the end of the run: this is a simplified hypothesis, that limits the development of significant mantle corner flow. We will discuss this limitation later on in the light of our results. On the other hand, the advantage of our approach compared to standard thermal models of mantle wedge corner flow, is that we do not fix the position or motion of plates in the subduction zone area during the run, thus we do not need to define special wedge boundary conditions (e.g. [START_REF] Conder | A case for hot slab surface temperatures in numerical viscous flow models of subduction zones with improved fault zone parameterization[END_REF][START_REF] Van Keken | Highresolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle[END_REF], such as a progressive transition between the fixed part and the subducting part of the oceanic plate. Here in our models, no condition other than the initial temperature is applied, since the 10 km thick channel interface, like all other layers, is selfconsistently defined with a temperature dependent brittleductile behaviour, located in between both deformable oceanic and continental plates. The approach is thus resolutely different.

3.2.3Mechanical behaviour of the Reference model

Elements of the model progressively load and deform as they are submitted to the gravity force and the kinematic boundary conditions. They reach a stress threshold that depends on their rheological properties: warm zones deform ductilely, while cold zones deform brittlely. We will here mainly describe deep stresses and temperatures in order to connect them with PressureTemperature (PT) petrological inferences from our tomography data. Stress and strain patterns in the continental crust will be only briefly described here, as it is not the aim of this paper to discuss this part of the modeling study, but will be done in another paper.

-In the oceanic lithosphere (green domain in Figures 5 and6c), flexural stresses develop due to its bended geometry. Before the trench at X~50 km, brittle shear zones develop at the surface, extending down to several tens of kms into the oceanic mantle (Figure 6b). These outerrise faults are compatible with previous observations of deep outerrise events along the Chilean margin [START_REF] Clouard | Outer rise stress changes related to the subduction of the Juan Fernandez Ridge, central Chile[END_REF]. As the oceanic lithosphere dips under the continent, tensile stresses develop in direction parallel to the applied slabpull velocity (Figure 6e), together with high compressional stress on the top portion of the slab down to about 130 km depth (commented in section 3.3). -In the subduction channel (yellow to light orange zones in figures 5, 6c), shear deformation develops, at relatively low shear stress (Figures 6d). At the surface a trench forms about 8 km deep. The shear stress reaches ~40 MPa down to ~60 km depth (Figure 7a), consistent with previous estimates (e.g. [START_REF] Lamb | Shear stresses on megathrusts: implications for mountain building behind subduction zones[END_REF], and references therein), and corresponds approximatively to a brittleductile transition. Above, deformation occurs dominantly by brittle failure according to the MohrCoulomb criterion at constant friction angle (1.43), and corresponding to the seismogenic zone. Below 80km depth, the subduction channel behaves ductile with a low effective viscosity (<10 20 Pa.s., Figure 7b). -The continental upper crust (brown, Figures 5,6c) deforms mostly brittle, and is partly ductile at ~10 km depth in the Bolivian area (X>~550 km, Figure 6b). The lower crust (red, Figure 6c) remains competent and deforms little in the forearc and in the cratonic easternmost domain, in contrast to the Altiplano, where weaker rheology and warmer temperatures enable to develop ductile shear deformation. This area below the plateau links deformation from the subducting slab, eastwards to the surface in the Eastern Bolivian weak upper crust (e.g. [START_REF] Sobolev | What drives orogeny in the Andes?[END_REF]. This mode of simple shear deformation is consistent with previous propositions such as by [START_REF] Tassara | Interaction between the Nazca and South American plates and formation of the Altiplano-Puna plateau: Review of a flexural analysis along the Andean margin (15°-34°S)[END_REF][START_REF] Tassara | Elastic thickness structure of South America estimated using wavelets and satellitederived gravity data[END_REF]. -As a whole, the forearc area behaves as a rigid lithosphericscale wedge, bounded to the east by the Western Cordillera, and extending down to ~110 km depth. This rigid block is cold and composed of strong crust and mantle rheologies, therefore deforms little (Figure 6b, e.g. [START_REF] Tassara | Elastic thickness structure of South America estimated using wavelets and satellitederived gravity data[END_REF]. Observed seismicity occurs in the range <400°C for crustal domains and <600°C for mantle domains (Figure 6c). However, stresses are not homogeneous (Figure 6de): seismicity occurs first at the high stress contact of the crust with the subduction channel, and then further east in slightly lower deviatoric stress areas. Principal stresses orientations indicate that this seismicity occurs in a downdip tensional zone at 3080 km depth (red crosses at X=140230 km), and where the most compressive stress (blue crosses) rotates from horizontal to vertical at 1040 km depth (X~220300 km). -The continental mantle (blue, Figures 5,6c) deforms little at temperatures lower than ~600°C, while it deforms ductile at greater temperatures, and with negligible shear stress (power law temperature dependant olivine rheology). East of and along the subduction channel below the forearc, cold temperatures induce significant shear stress (around 100 MPa, figure 6d) in the lithospheric mantle down to about 100 km depth. This cold "rigid mantle wedge" (e.g. [START_REF] Honda | Thermal structure beneath Tohoku, Northeast Japan a case study for understanding the detailed thermal structure of the subduction zone[END_REF], or more recently named the cold nose [START_REF] Abers | The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow[END_REF] is delimited on top by the base of the forearc crust, at its base by the warmer and ductile corner wedge, and to the west by the nearly vertical ~600°C olivine brittleductile transition. While recent thermal studies show the high capacity of mantle corner flow in heating the adjacent slab and the overriding plate (e.g.; [START_REF] Billen | A low viscosity wedge in subduction zones[END_REF][START_REF] Van Keken | Highresolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle[END_REF][START_REF] Currie | The thermal structure of subduction zone back arcs[END_REF][START_REF] Kelemen | Thermal structure due to solidstate flow in the mantle wedge beneath arcs[END_REF], coherent with heat flow and geophysical characteristics of several subduction zones, other thermal models also show that the development of a cold forearc mantle wedge is favored when the oceanic plate is old and subducts fast, and when a strong coupling is accounted for in between the fixed and descending parts of the oceanic plate (e.g. [START_REF] Kneller | Btype olivine fabric in the mantle wedge: Insights from highresolution nonNewtonian subduction zone models[END_REF][START_REF] Abers | The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow[END_REF]. The data from Northern Chile indicates such characteristics, with [START_REF] Springer | Interpretation of heatflow density in the Central Andes[END_REF] thermal model indicating a shallowest depth of the asthenospheric wedge located 70 km below the Western Cordillera, and vertically descending isotherms down to about 100 km depth before intersecting with the slab (such verticalisation is also assumed by Abers et al., 2006, for Cascadia). This cold forearc mantle remains the strongest layer of the continental lithosphere, and acts as a butress to subduction, thus controlling deformation of the upper plate. It transmits compressional stresses from the upper plate, that apply perpendicular to the subduction zone, into the slab. The stresses that we obtain are in the range of those cited by [START_REF] Kneller | Btype olivine fabric in the mantle wedge: Insights from highresolution nonNewtonian subduction zone models[END_REF][START_REF] Kneller | Stress, strain, and Btype olivine fabric in the forearc mantle: Sensitivity tests using highresolution steadystate subduction zone models[END_REF], who explain trenchparallel fast shearwave splitting directions by the presence of Btype olivine fabric in the forearc mantle. Heading East below the high plateau, high temperatures (1350°C at ~60km depth) provide a low viscosity behaviour that may be associated to active mantle corner flow, coherent with magnetotelluric and seismological data. Below the Eastern Cordillera (X~500 km) and eastward, progressive cooling recomposes a rigid, cratonic like behavior down to >150 km depth.

3.3DSZ PressureTemperature conditions inferred from the thermomechanical model

We display now the modeled Pressure-Temperature (PT) conditions in the double seismic zone, spread between 80 and 140 km depth along the slab. Figure 7 displays the thermal and stress conditions within the subducting channel crust, and Figure 8 shows PressureTemperature conditions within the slab. We superimpose the earthquake locations and the model (Fig. 6) by assuming that at least the top part of the upper plane of earthquakes should be located inside the subducting oceanic crust channel (yellow in figure 5). The uncertainty in this «visual» fit is about 5 km, due to the projection of seismic profiles on the curved in plane geometry of the trench at these latitudes, and to the thermomechanical grid size. Whereas it is common to assume that WadatiBenioff earthquakes delineate the subduction channel, an uncertainty remains on this assumption, see for example a depth dependent increasing decay in Abers et al. (2006), for the Cascades subduction zone. However, note that in this latter study, this correspondenc is preserved in the upper part of the WadatiBenioff zone.

Upper DSZ plane

Figure 7c plots minimal and maximal temperatures within the subduction channel: whereas maximum temperatures reach 750°C at the top of the subducting crust, they drastically decrease to about 300°C at the bottom base of the subducting crust, which in the model was defined 10 km thick. In between 90 and 130 km depth, minimum temperatures are ~100°C lower than those estimated by Abers et al. (2006) for the, yet shallower, seismic portion of central Cascadia subduction zone, where a younger 38 Ma oceanic lithosphere converges at a speed of 55mm/yr. Maximum temperatures, however are close to 700°C, which is about 100°C less than what [START_REF] Van Keken | Highresolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle[END_REF] obtained for the top of the Honshu slab below the volcanic arc. Considering the equivalent surface heatflow between both subducting plates in Honshu and Northern Chile, we evaluate that this 100°C difference is linked with the obviously warmer general temperatures of the Honshu upper plate with respect to Northern Chile (see previous sections).

We have here a possible explanation for the distance measured by [START_REF] Abers | The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow[END_REF] in Cascadia, between the earthquakes of the WadatiBenioff zone and the low seismic velocity zone: the latter could represent the warmest and thus ductile top of the crust, whereas the former could represent the coldest and thus seismically active base of the subducting crust.

In our prior estimates of petrological compositions using [START_REF] Hacker | Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and H 2 O contents of minerals and rocks at pressure and temperature[END_REF] procedure, Vp and Vs values of Northern Chile DSZ upper plane were best modeled with waterrich (up to 5.7 wt% H 2 O), and mineral assemblages comprising omphacite, lawsonite, amphibole and chlorite (section 2.2.3). Thermodynamically calculated mineral equilibria for metabasaltic compositions (i. e., [START_REF] Guiraud | Calculated mineral equilibria in the greenschistblueschisteclogite facies in Na 2 OFeOMgOAl 2 O 3 SiO 2 H 2 O. Methods, results and geological applications[END_REF][START_REF] Clarke | Eclogites and Blueschists of the Pam Peninsula, NE New Caledonia: a reappraisal[END_REF]Carlson et al., 1999) all indicate that such mineral assemblages only occur at temperatures below 570°C and are typical of blueschist conditions. In addition, the probable persistence of lawsonite in mineral assemblages suggests that temperatures, even at 3.5 GPa (~140 km) do not exceed ~700°C, the temperature at which lawsonite is no longer stable. Thus, independent mineralogical data are consistent with the low temperature conditions of our thermomechanically modelled crust. Then, a fluid release associated to progressive breakdown of H 2 Orich phases, such as chlorite and amphibole, appears as the most probable scenario to explain this upper plane seismicity.

Intermediate zone

Our previous petrological estimates for this zone was a partially hydrated harzburgitic oceanic mantle, with however, remaining high Vp values with respect to Vs and Vp/Vs values (section 2.2.4). We note from the present thermomechanical models that along the subduction zone, the DSZ seismicity contours and bounds the shape of the innerslab compressive zone (Figure 8a) from 80 to 130 km depth. The intermediate zone is thus in a compressional state of stress, at approximately at least 400 MPa above lithostatic stress. Such overlithostatic pressures of several hundreds of MPa tend to increase Vp (Hacker et al., 2003a), thus offering an explanation to the relative high numbers measured in Northern Chile. We may therefore evaluate that this intermediate zone is even less hydrated than previously estimated.

3.3.3Lower DSZ plane

If the upper plane's seismicity represents indeed the subducting oceanic crust, then the observed lower plane's seismicity falls in the model, on the neutral plane of the subducting mantle located about 20 km below the top of the slab (Figures 8a,8c). There, pressure becomes lithostatic as it switches from downdip compression to downdip tension (Figure 6e,8c). This neutral plane could form a secondary channel of «relaxed» pressure, into which fluids from lateral areas could migrate: tensile stresses below the neutral plane should help open pores and let light fluids migrate upwards, whereas compressional stresses above the neutral plane should close pores and create a permeability barrrier to fluid ascent. Consequently, free fluids could remain trapped in this neutral plane. Earthquakes may thus be triggered close to and on either side of this neutral plane, due to the presence of freefluids. Within a partially hydrated harzburgitic material (the most probable candidate deduced insection 2.2.5), such a freefluid phase pore configuration could explain the observed low Vp/Vs ratios, as suggested by many authors (i.e., [START_REF] Nakajima | Threedimensional structure of Vp, Vs, and Vp/Vs beneath northeastern Japan: Implications for arc mgmatism and fluids[END_REF][START_REF] Takei | Effect of pore geometry on Vp/Vs: From equilibrium geometry to crack[END_REF]Hacker et al., 2003a).

Seismic events occurring in this lower DSZ plane should therefore not present any dominant compressional or extensional focal mechanism, but will rather occur depending on their position either slightly above or slightly below this neutral fiber. In fact, previous studies worldwide cannot isolate a common trend for the location of compressional or extensional events in the upper and lower plane of DSZs respectively, nor the opposite (cf. review in [START_REF] Peacock | Are the lower planes of double seismic zones caused by serpentinite dehydratation in subducting oceanic mantle?[END_REF]. Our proposition here differs from previous studies that link upper and lower seismic planes to high deviatoric stresses in the downgoing crust and mantle (e.g. [START_REF] Wang | Unbending combined with dehydration embrittlement as a cause for double and triple seismic zones[END_REF], and for which the elastic core is located in between both seismic planes. Here instead, we propose that the plate's stress neutral fiber has not kept a significantly thick elastic core (e.g. [START_REF] Sleep | The double seismic zone in downgoing slabs and the viscosity of the mesosphere[END_REF], but consists of an already fractured area (such as by outerrise faulting) in which fluids have accumulated and trigger neighbouring seismogenic dehydration reactions.

Note that in Figure 8c, the minimal distance between the subducting channel and the neutral lithostatic stress is closer to 15 km than to 20 km. We think that this is due to our choice of a relatively warm subducting oceanic lithosphere (60 Myr old) with respect to real conditions infered from surface heatflow values (corresponding to classical 100 Myr old oceanic lithospheres). With a colder lithosphere, this distance increases, as indicated by Brudzinski et al.'s (2007) study and an alternate cold model (section 3.4). Furthermore, if instead we try to associate the upper and lower plane seismicity to the compressional and extensional zones of the subducting lithosphere respectively, then we see that their distance in the model is 25 km (Fig. 8c), which is not satisfying either, and will even increase if a colder oceanic lithosphere was assumed.

In general, our modeled temperatures are about 100°C lower than those predicted for the Honshu lower seismic plane. Petrological interpretations indicate conditions lower than about 500600°C, while our thermomechanical models indicate conditions around 400°C. There are several possible explanations:

1 Temperatures in the Northern Chile subduction zone are anomalously cold: first the Nazca plate has anomalously low surface heat flow (30mW/m 2 instead of the theoretical 70mW/m 2 , e.g. [START_REF] Springer | Heatflow density across the Central Andean subduction zone[END_REF], then the overriding plate presents geophysical characteristics of a rigid cold forearc as opposed to a very weak and warm Altiplano, indicative of an eastward shift of the warming mantle wedge corner flow (e.g. [START_REF] Haberland | Coincident anomalies of seismic attenuation and electrical resistivity beneath the southern Bolivian Altiplano plateau[END_REF][START_REF] Brasse | The mantle wedge in the Bolivian orocline in the view of deep electromagnetic soundings[END_REF]. [START_REF] Honda | Thermal structure beneath Tohoku, Northeast Japan a case study for understanding the detailed thermal structure of the subduction zone[END_REF] had already suggested that the mantle forearc of central South America is probably cooler than that of the Honshu region, possibly associated to the different nature of volcanic rocks in both regions. More recently, [START_REF] Kneller | Stress, strain, and Btype olivine fabric in the forearc mantle: Sensitivity tests using highresolution steadystate subduction zone models[END_REF] proposed to explain trenchparallel fast shearwave splitting directions in the forearc mantle of several subduction zones, including Northern Chile, by the presence of Btype olivine fabric, associated to a relatively cold and a high stress environment. Since we used in our models a Nazca plate age of 60 Myr giving surface heatflow values of 60mW/m 2 , we may still be overestimating innerslab real temperatures.

2 The numerical models carried here have assumed simple thermal properties that are basically diffusing and advecting during a short 2 Ma time: there is no account for thermal expansion, nor for latent heat linked to phase transformations. Certainly the longlived subduction zone of Chile incorporates greater complexity than the present models do. Several studies showed that powerlaw viscous flow in the mantle wedge is capable of heating the plate interface temperatures by about 200°C compared to isoviscous cases (e.g. [START_REF] Van Keken | Highresolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle[END_REF][START_REF] Conder | A case for hot slab surface temperatures in numerical viscous flow models of subduction zones with improved fault zone parameterization[END_REF]. Although here we use a similar powerlaw rheology, the short duration of the model (2 Ma) limits the efficiency of convective heat transport. We may thus be underestimating slab interface temperatures. However, [START_REF] Van Keken | Highresolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle[END_REF] also showed that thermal gradients towards the innercore of the slab may remain high, especially in cases of fast subduction of cold lithosphere, thus preserving its initial cold state to great depths. Consequently, we may be underestimating temperatures at the top of the slab, but not that much innerslab temperatures.

Testing various parameters

The reference model of Figure 6 is our preferred model with which we discussed the PT conditions for the DSZ of Northern Chile. But in order to obtain these appropriate PT conditions, we had to tune a number of initial parameters. Figure 9 shows a set of 5 alternate models, which illustrate the effect of varying some parameters with respect to the reference model (Figure 9a, Table 6):

-First, if we take a 100 Ma old subducting oceanic plate (Figure 9b), then the modeled surface heat flow better matches the data from Northern Chile. The 400°C isotherm and the compressional domain on top of the slab extend deeper down to about 180 km, in good correlation with the deeper DSZ of Honshu (plate age of 130 Ma). Our petrological interpretations, which strongly suggest a cold oceanic plate, still agree with this thermal setting. -Second, if we take warmer conditions in the overriding plate (Figure 9c), the intensity of the compressional zone vanishes. In this model, the continental thermal thickness (hl) was decreased to 120 km, and the thermal anomaly was increased. There, the warm forearc mantle flows ductilely and cannot transmit that much high stresses from the upper plate into the subducting plate. Other lines of evidence indicate that the forearc mantle in Northern Chile is competent, such as thermal studies and gravity modeling (e.g. [START_REF] Honda | Thermal structure beneath Tohoku, Northeast Japan a case study for understanding the detailed thermal structure of the subduction zone[END_REF][START_REF] Tassara | Elastic thickness structure of South America estimated using wavelets and satellitederived gravity data[END_REF]. Also, it must be competent enough in order to sustain the stresses generated by the extremely thick crust of the Central Andes. In the opposite case, the Altiplano would collapse and present day deformation would concentrate within the forearc and highplateau areas, and this is not observed. Deformation now mainly occurs in the Bolivian fore land.

-Third, the composition of the subduction channel (Figure 9d) was modified with ductile creep parameters corresponding to plagioclase, replacing wet granite in the reference model. Now, viscosities are greater in the subduction channel, and higher shear stresses develop there (comparison with the reference model in Figure 7a,b). Shear stresses in the subducting channel below about 80 km depth, absorb part of the normal (and thus compressional) stresses acting on one side and the other of the subduction zone, thus reducing the compressional component that is transmitted from the upper plate into the slab. -Fourth, the radius of curvature of the subduction zone was reduced to 250 km instead of 450 km (Figure 9e). The radius of curvature, r, estimated from WadatiBenioff zones, mostly ranges from 200 to 500 km (e.g. [START_REF] Wang | Unbending combined with dehydration embrittlement as a cause for double and triple seismic zones[END_REF]. While Honshu and Chile have a similar radius of 450500 km, New Zealand's dipping seismic zone approximates a 300 km arc radius [START_REF] Reyners | Imaging subduction from the trench to 300 km depth beneath the central North Island, New Zealand, with Vp and Vp/Vs[END_REF]. In this model, compression on top of the slab extend to a shallower depth, consistent with shallowing of the plates unbending depth, but is also extends to greater depth, indicating that plate unbending there, has no more role to play. The slab is seen here to act like a vertical wall that « butresses » horizontal compression.

-Fifth, no upper plate convergence was applied (Figure 9f), the right boundary of the model was maintained fixed in the horizontal direction. In addition, the overthickened crust of the Altiplano was removed, in order to inhibit the local lateral compression exerted by gravitational forces. Now, the intensity of compression on top of the slab is much lower, showing that upper plate convergence is indeed important in the development of innerslab compression.

These alternate models demonstrate that the innerslab flexural compressional zone does not only depend on the subducting plate age, but also on the regional stress environment. Transmission of compressional stresses from the upper plate into the subducting plate in turn, depends on lubrification of the subducting channel, thermal and rheological states of the overriding plate, and relative plate convergence.

Conclusion

The results of a high resolution tomography of the Nazca plate Double Seismic Zone beneath North Chile, combined with a thermomechanical numerical modelling of the subduction zone, provides new insights on the petrological and thermodynamical properties of the oceanic plate and overriding continental forearc. The velocity structure of the upper 140 km of the slab has been obtained, taking advantage of the two seismic planes observed at depths 80 to 140 km. The results highlight the specific characteristics of the zones where earthquakes occur. Moreover, the geodynamical setting leading to the observation of the DSZ, as well as its pressure and temperature conditions have been investigated by the thermomechanical modeling. In summary:

the DSZ in Northern Chile is observed between 80 to 140 km depth. Both seismic planes are separated by ~20 km. -The upper seismic plane, corresponding to oceanic crust, is characterised by intermediate Vp (~7.7 km/s) and low Vp/Vs (1.67). These values suggest that between 80 and 140 km depth, the seismic upper plane is related to fluids releases associated to metamorphic reactions occurring within jadeite lawsonite blueschists and perhaps, at higher depths (~130 km) under lawsoniteamphibole eclogite conditions. -Seismicity in the lower plane, lying in the oceanic mantle, yields low Vp (~7.4 km/s) and low Vp/Vs (~1.57). These values are consistent with the common hypothesis of seismic activity associated to water release resulting from brucite, chlorite and antigorite breakdown in dehydration reactions occurring within serpentinized harzburgite (antigorite ≤ 25 vol. %). -The region in between the two planes, with high Vp (8.5 km/s) and high Vp/Vs (1.9), presents characteristics corresponding to a moderately serpentinized harzburgite (antigorite ≤10 vol. %). -Our petrological estimates infer general thermal conditions lower than 500°C. This low temperature range in comparison to estimates from other subduction zones, is justified by the extremely anomalous structure of the Nazca plate, which has surface heatflow values of ~30mW/m 2 compared to its theoretical 70mW/m 2 (considering its 50 Myr age, [START_REF] Springer | Heatflow density across the Central Andean subduction zone[END_REF]. Furthermore, the forearc mantle into which the Nazca plate subducts is also probably cold and only partially serpentinised, as indicated by previously acquired geophysical data such as Qp attenuation, magnetotellurics or shear wave splitting data. This generally cold state of the subduction zone may have significant implications on the kinetics of metamorphic reactions and on the rock structures.

Although water saturated, the slab could display offsets between real parageneses and predicted ones estimated from equilibrium thermodynamics petrogenetic grids. This feature, together with the fact that North Chile DSZ most probably contains partially hydrated rock zones, suggests texturally heterogeneous environments in which rock fragments of metastable mineral assemblages (i. e., gabbro, peridotite) persist within metamorphic rocks of higher metamorphic grade (i. e., blueschist, eclogite) or more hydrated rocks (i. e., serpentinite). Global velocity properties of such environments could be strongly affected by this textural heterogeneity, such as Vp values being more sensitive to the water content of a weak matrix and Vs values being more sensitive to rigid nodules. -Thermomechanical models consistent with a cold environment reproduce the deformation and stress patterns of the subduction zone. In agreement with previous studies, deformation propagates from the subduction channel shearzone, into the Altiplano lower crust and the Subandean zone, whereas the forearc area remains rigid and sustains high stresses down into forearc mantle. When these thermo mechanical models are superimposed on the seismic data, the DSZ is predicted to lie within the 350 450°C range that contours an innerslab compressional zone from 80 to 140 km depth. This cold temperature range is about 100°C lower than predictions from thermal models for other, yet probably warmer subduction zones. Alternate tests show that warmer slab and forearc conditions reduce the extent and intensity of this innerslab compressional zone. -From the thermomechanical models, lower plane earthquakes are predicted to occur in the neutral fiber plane between downdip tension and downdip compression, where we suggest that fluids released from surrounding metamorphic reactions can migrate. We think that tensile stresses below this neutral plane open pores and favor upward migration of freefluids, whereas compressive stresses above this neutral plane closes pores and creates a permeability barrier able to trap fluids along a thinly sheeted serpentinized layer. This would be consistent with the generally admitted relation between dehydration reactions and seismicity [START_REF] Nakajima | Threedimensional structure of Vp, Vs, and Vp/Vs beneath northeastern Japan: Implications for arc mgmatism and fluids[END_REF][START_REF] Takei | Effect of pore geometry on Vp/Vs: From equilibrium geometry to crack[END_REF]Hacker et al., 2003a). Therefore these lower plane earthquakes are likely to present both downdip compression and downdip extension characteristics, precisely because of their location about the neutral plane. In between both seismic planes, the models produce compressional stresses of several hundreds MPa: this should affect seismic velocities, and in particular, reduce the amount of fluids required to explain observed values. -Numerical tests indicate that the innerslab compressional zone depends on the following parameters:

1) an older plate yields deeper isotherms and also deeper compressional zone. 2) a weak ductile subduction channel favors innerslab compression, by facilitating transmission of compressional stresses from the upper plate. 3) A cold structure of the continental forearc also favors innerslab compression, for the same reason. 4) Decreasing subducting plate radius of curvature broadens the innerslab compression upward and downwards. 5) decreasing upper plate convergence diminishes it. All these results indicate that if indeed DSZ's contours innerslab compression in general, they cannot only be linked to slab unbending (e.g. see discussion in [START_REF] Abers | Relationship between shallowand intermediate depth seismicity in the Eastern Aleutian subduction zone[END_REF], but also to the transmission of compressional stresses from the upper plate into the slab, as had suggested [START_REF] Fujita | Double seismic zones and stresses of intermediate depth earthquakes[END_REF] Aknowledgements: C. D. thanks Clifford Thurber for inviting her to attend the workshop on seismic waveform alignment, earthquake location, and seismic tomography, held at the Department of Geology and Geophysics, UW Madison, from May 1721, 2004. M.G. thanks Andres Tassara for helping constraining the models for another issue, Jaime Campos for providing office facilities at Univ. de Chile, system administrators Rodrigo Sanchez in Chile and Guillaume Pernot in Toulouse without who the several hundreds of runs couldn't have been done, Mathurin Lebris for his preliminary implementation of passive markers, and Marko Testi for his magical Matlab tricks.

Figure 1: Epicenters of the 2052 events used in the tomography, the darker, the deeper. Red and blue triangles are stations. In gray, the inversion grid used for the inversion. The black cross is the origin of the model (18°S, 69°W). Thin red lines locate the three profiles P1, P2 and P3 presented Figure 3, and the thick red line, the cross section location shown in Figure 4. Insert: Central Andes tectonic setting between 17° and 20.5°S, from [START_REF] Tassara | Interaction between the Nazca and South American plates and formation of the Altiplano-Puna plateau: Review of a flexural analysis along the Andean margin (15°-34°S)[END_REF]. cc: Coastal Cordillera, id: Intermediate Depression, dc: Domeyko Cordillera, wc: Western Cordillera, ap: Altiplano. Black vector: Nazca plate absolute motion. Figure 2: Initial layered model (blue) and final Vp model (red) obtained by 1D inversion. Figure 3: Cross section through the Vp and Vs values along the 3 profiles shown on Figure 1. Plotted earthquakes are within a 20 km wide section. At depths 80-140 km, the double seismic zone is well defined. Both seismic planes are 5 to 10 km thick with a separation of 20 km. Structures above the white line, where derivative weight sum for each node is more than 10, is meaningful. Figure 4: Zoom on the velocity structure of the DSZ: cross sections through the Vp and Vs values in the central part of profile P2 (Y=60 km). Figure 5: Initial conditions for the thermomecanical model: each colored layer has specific rheological density and thermal properties (see Tables 3,4,5). Isotherms 400, 600, 1000, 1300°C are contoured. Asthenospheric and continental mantle in blue, crustal layers in red to brown, oceanic plate in green (crust is in orange), subducting channel in yellow and light orange). Initial topography is evaluated from initial isostatic equilibrium. Boundary conditions are lithostatic compensation on the sides and at the base of the model (narrow black arrows), except at the base of the slab and along the right edge of the continental lithosphere, where velocities of respectively 5cm/yr (slabpull) and 3cm/yr (upper plate convergence) are applied (bold arrows). Vertical sections of pressure (tension > 0 here, and the lithostatic component is subtracted) throughout the forearc, locating from left to right the neutral plane at depths 80, 100, 120 and 140 km (green circles). Note the ~15 km distance of this neutral plane with the subducting channel crust above, whereas the distance between maximum compression and maximum tension is greater than 25 km. Note the variable stress state, including tension, in the forearc mantle and in the channel area (pink ellipses). Appendix 4: Thermomecanical markers procedure for tracking pressure and temperature.

Remeshing in Lagrangian methods is a common source of numerical diffusion, that the use of passive markers diminish, but cannot completely inhibit (e.g. [START_REF] Yamato | Burial and exhumation in a subduction wedge: Mutual constraints from thermomechanical modeling and natural PTt data (Schistes Lustrés, western Alps)[END_REF]. In Parovoz like in other FLAC methods [START_REF] Cundall | A microcomputer program for modeling largestrain plasticity problems[END_REF], stresses and strains are computed at each timestep on the triangle subelements of the mesh (each element is made of 2 superimposed pairs of triangles that allow accurate calculation of stress and strain components), and temperatures and velocities are nodal values). At the onset of the run, we inserted 9 regularly spaced markers in each mesh element, that passively move with the Lagrangian grid during normal computation, and carry linearly interpolated, elementary stresses and temperatures. All 9 markers of a same element are assigned identical element material properties, grouping density, conductivity, and rheological parameters, represented by a colour "phase" in the setup Figure 5. When the angle of any triangular element becomes smaller than a critical value, set to 13°here, remeshing occurs and a new regular mesh replaces the old one. Then, the markers are counted in each new element, in order to identify the dominant material properties ("phase") that composes it. During remeshing, markers are also used to evaluate subelement stress components (and therefore, pressure): new stress components are averaged over that of all markers located inside the new subelement. For temperature, new nodal and markers temperatures are linearly interpolated from the old nodal mesh.

In Figure 8b, we plot markers values in the pressure neutral plane of the slab (Plithos in green), and in the plane located 200 MPa above and below (Plithos+200MPa in blue and Plithos200MPa in red respectively), within a ±15 MPa range. Coloured dots are markers, coloured lines show averaged local values, and black lines indicate depth zone. We see that markers of a same colour and at similar depths follow inclined straight patterns, that correspond to markers belonging to a same subelement: since pressure is an elementary value, markers thus have constant pressure, but slightly different lithostatic pressure because they are located at different depths (within the ±15 MPa interval). On the other hand, markers temperatures are interpolated from nodal values, and therefore do not present such artificially grouped patterns. We see that temperatures vary, at given depth, by about 20°C (horizontal extent of dots): this value is constrained by our choice of the ±15 MPa pressure interval. If we took a smaller pressure interval, then "holes" would appear at certain depths zones, because there would be no markers of that pressure range, in that zone. These "holes" are typically 2 to 3 km wide, equal to the mesh resolution there. At a given depth, the three coloured lines show a relative pressure offset of 200 MPa (as prescribed), and a temperature difference of about 30°C, which is greater than the 20°C range that individuates each pressure range (width area of same colour dots on Figure 8b). Consequently our mesh resolution allows to differentiate these three domains.

  : this test estimates the ability of the real data to recover the model obtained from the inversion. The input velocity models, which are the inverted true models, are shown in Figures C of Appendices 2 and 3, and the corresponding recovered velocity models are shown in Figures D. The velocity features of the subduction zone are well recovered in position and magnitude within the upper plane. Within the lower plane, the recovered velocities are decreased. Nevertheless, this test shows that if the observed anomalies exist, the dataset enables to recover them. -the last test starts from an input model without any velocity anomaly (Appendices 2 and 3, Figures E).

Figure 6 :

 6 Figure1: Epicenters of the 2052 events used in the tomography, the darker, the deeper. Red and blue triangles are stations. In gray, the inversion grid used for the inversion. The black cross is the origin of the model (18°S, 69°W). Thin red lines locate the three profiles P1, P2 and P3 presented Figure3, and the thick red line, the cross section location shown in Figure4. Insert: Central Andes tectonic setting between 17° and 20.5°S, from[START_REF] Tassara | Interaction between the Nazca and South American plates and formation of the Altiplano-Puna plateau: Review of a flexural analysis along the Andean margin (15°-34°S)[END_REF]. cc: Coastal Cordillera, id: Intermediate Depression, dc: Domeyko Cordillera, wc: Western Cordillera, ap: Altiplano. Black vector: Nazca plate absolute motion. Figure2: Initial layered model (blue) and final Vp model (red) obtained by 1D inversion. Figure3: Cross section through the Vp and Vs values along the 3 profiles shown on Figure1. Plotted earthquakes are within a 20 km wide section. At depths 80-140 km, the double seismic zone is well defined. Both seismic planes are 5 to 10 km thick with a separation of 20 km. Structures above the white line, where derivative weight sum for each node is more than 10, is meaningful. Figure4: Zoom on the velocity structure of the DSZ: cross sections through the Vp and Vs values in the central part of profile P2 (Y=60 km). Figure5: Initial conditions for the thermomecanical model: each colored layer has specific rheological density and thermal properties (see Tables3,4,5). Isotherms 400, 600, 1000, 1300°C are contoured. Asthenospheric and continental mantle in blue, crustal layers in red to brown, oceanic plate in green (crust is in orange), subducting channel in yellow and light orange). Initial topography is evaluated from initial isostatic equilibrium. Boundary conditions are lithostatic compensation on the sides and at the base of the model (narrow black arrows), except at the base of the slab and along the right edge of the continental lithosphere, where velocities of respectively 5cm/yr (slabpull) and 3cm/yr (upper plate convergence) are applied (bold arrows). Figure 6: Reference model,with superimposed seismicity. From top to Bottom: a) modeled topography in blue, erosion rate in green, and real topography at 18.5°S in black. b) Shear strain rate and velocity vectors, indicating in red areas of greater shear deformation, either brittle or ductile (note the outerrises). c) Rheological layers and isotherms: note the 400°C isotherm. d) Deviatoric shear stress indicating in red areas of stress greater than 250 MPa, in the slab, in the crust. Note some stress in the deep portions of the forearc. e) Crosses of the principal stresses, zoomed at the subduction zone, showing downdip compression above downdip tension, and rotating and complexe stresses in the forearc. Figure 7: Thermomechanical properties of the subduction channel for reference model (red) and model with prescribed higher viscosity in the subduction channel (blue, model d figure 8). a) deviatoric shear stress, indicating for the reference model a maximum at ~30 km depth of 40 MPa, reducing to nearly 0 below ~80 km depth. b) effective viscosity, which becomes lower than 10^20 Pa.s. in the reference model, to the difference of model d.c) Minimum (dashed lines) and maximum (plain lines) temperatures. Note that the great increase around ~90km depth corresonds to the detachement and elevation of light oceanic crust material throughout the continental mantle (see yellow bubble in figure6c). Figure8: PressureTemperature conditions obtained for the suggested lower seismic plane, located at the depth of the neutral plane in the subducting slab: a) Background colors show pressure, with subtraction of the lithostatic component (compression is negative and blue). b) Pressure temperature paths: green markers (dots) and green line (mean value ) show the neutral plane within the slab (labeled Plithos), indicating temperatures of 320~380°C at lithostatic pressure of 2.55 GPa. Pressure and temperature paths are also shown in blue and red, for areas above and below this neutral plane, which are at an overpressure and underpressure of 200 Mpa respectively, corresponding to the compressive stress (labeled Plithos+200MPa) and tensile stress (Plithos200MPa) of Figure8a. See appendix 4 for numerical sampling details. c)

Figure 7 :

 7 Figure1: Epicenters of the 2052 events used in the tomography, the darker, the deeper. Red and blue triangles are stations. In gray, the inversion grid used for the inversion. The black cross is the origin of the model (18°S, 69°W). Thin red lines locate the three profiles P1, P2 and P3 presented Figure3, and the thick red line, the cross section location shown in Figure4. Insert: Central Andes tectonic setting between 17° and 20.5°S, from[START_REF] Tassara | Interaction between the Nazca and South American plates and formation of the Altiplano-Puna plateau: Review of a flexural analysis along the Andean margin (15°-34°S)[END_REF]. cc: Coastal Cordillera, id: Intermediate Depression, dc: Domeyko Cordillera, wc: Western Cordillera, ap: Altiplano. Black vector: Nazca plate absolute motion. Figure2: Initial layered model (blue) and final Vp model (red) obtained by 1D inversion. Figure3: Cross section through the Vp and Vs values along the 3 profiles shown on Figure1. Plotted earthquakes are within a 20 km wide section. At depths 80-140 km, the double seismic zone is well defined. Both seismic planes are 5 to 10 km thick with a separation of 20 km. Structures above the white line, where derivative weight sum for each node is more than 10, is meaningful. Figure4: Zoom on the velocity structure of the DSZ: cross sections through the Vp and Vs values in the central part of profile P2 (Y=60 km). Figure5: Initial conditions for the thermomecanical model: each colored layer has specific rheological density and thermal properties (see Tables3,4,5). Isotherms 400, 600, 1000, 1300°C are contoured. Asthenospheric and continental mantle in blue, crustal layers in red to brown, oceanic plate in green (crust is in orange), subducting channel in yellow and light orange). Initial topography is evaluated from initial isostatic equilibrium. Boundary conditions are lithostatic compensation on the sides and at the base of the model (narrow black arrows), except at the base of the slab and along the right edge of the continental lithosphere, where velocities of respectively 5cm/yr (slabpull) and 3cm/yr (upper plate convergence) are applied (bold arrows). Figure 6: Reference model,with superimposed seismicity. From top to Bottom: a) modeled topography in blue, erosion rate in green, and real topography at 18.5°S in black. b) Shear strain rate and velocity vectors, indicating in red areas of greater shear deformation, either brittle or ductile (note the outerrises). c) Rheological layers and isotherms: note the 400°C isotherm. d) Deviatoric shear stress indicating in red areas of stress greater than 250 MPa, in the slab, in the crust. Note some stress in the deep portions of the forearc. e) Crosses of the principal stresses, zoomed at the subduction zone, showing downdip compression above downdip tension, and rotating and complexe stresses in the forearc. Figure 7: Thermomechanical properties of the subduction channel for reference model (red) and model with prescribed higher viscosity in the subduction channel (blue, model d figure 8). a) deviatoric shear stress, indicating for the reference model a maximum at ~30 km depth of 40 MPa, reducing to nearly 0 below ~80 km depth. b) effective viscosity, which becomes lower than 10^20 Pa.s. in the reference model, to the difference of model d.c) Minimum (dashed lines) and maximum (plain lines) temperatures. Note that the great increase around ~90km depth corresonds to the detachement and elevation of light oceanic crust material throughout the continental mantle (see yellow bubble in figure6c). Figure8: PressureTemperature conditions obtained for the suggested lower seismic plane, located at the depth of the neutral plane in the subducting slab: a) Background colors show pressure, with subtraction of the lithostatic component (compression is negative and blue). b) Pressure temperature paths: green markers (dots) and green line (mean value ) show the neutral plane within the slab (labeled Plithos), indicating temperatures of 320~380°C at lithostatic pressure of 2.55 GPa. Pressure and temperature paths are also shown in blue and red, for areas above and below this neutral plane, which are at an overpressure and underpressure of 200 Mpa respectively, corresponding to the compressive stress (labeled Plithos+200MPa) and tensile stress (Plithos200MPa) of Figure8a. See appendix 4 for numerical sampling details. c)

Figure 8 :

 8 Figure1: Epicenters of the 2052 events used in the tomography, the darker, the deeper. Red and blue triangles are stations. In gray, the inversion grid used for the inversion. The black cross is the origin of the model (18°S, 69°W). Thin red lines locate the three profiles P1, P2 and P3 presented Figure3, and the thick red line, the cross section location shown in Figure4. Insert: Central Andes tectonic setting between 17° and 20.5°S, from[START_REF] Tassara | Interaction between the Nazca and South American plates and formation of the Altiplano-Puna plateau: Review of a flexural analysis along the Andean margin (15°-34°S)[END_REF]. cc: Coastal Cordillera, id: Intermediate Depression, dc: Domeyko Cordillera, wc: Western Cordillera, ap: Altiplano. Black vector: Nazca plate absolute motion. Figure2: Initial layered model (blue) and final Vp model (red) obtained by 1D inversion. Figure3: Cross section through the Vp and Vs values along the 3 profiles shown on Figure1. Plotted earthquakes are within a 20 km wide section. At depths 80-140 km, the double seismic zone is well defined. Both seismic planes are 5 to 10 km thick with a separation of 20 km. Structures above the white line, where derivative weight sum for each node is more than 10, is meaningful. Figure4: Zoom on the velocity structure of the DSZ: cross sections through the Vp and Vs values in the central part of profile P2 (Y=60 km). Figure5: Initial conditions for the thermomecanical model: each colored layer has specific rheological density and thermal properties (see Tables3,4,5). Isotherms 400, 600, 1000, 1300°C are contoured. Asthenospheric and continental mantle in blue, crustal layers in red to brown, oceanic plate in green (crust is in orange), subducting channel in yellow and light orange). Initial topography is evaluated from initial isostatic equilibrium. Boundary conditions are lithostatic compensation on the sides and at the base of the model (narrow black arrows), except at the base of the slab and along the right edge of the continental lithosphere, where velocities of respectively 5cm/yr (slabpull) and 3cm/yr (upper plate convergence) are applied (bold arrows). Figure 6: Reference model,with superimposed seismicity. From top to Bottom: a) modeled topography in blue, erosion rate in green, and real topography at 18.5°S in black. b) Shear strain rate and velocity vectors, indicating in red areas of greater shear deformation, either brittle or ductile (note the outerrises). c) Rheological layers and isotherms: note the 400°C isotherm. d) Deviatoric shear stress indicating in red areas of stress greater than 250 MPa, in the slab, in the crust. Note some stress in the deep portions of the forearc. e) Crosses of the principal stresses, zoomed at the subduction zone, showing downdip compression above downdip tension, and rotating and complexe stresses in the forearc. Figure 7: Thermomechanical properties of the subduction channel for reference model (red) and model with prescribed higher viscosity in the subduction channel (blue, model d figure 8). a) deviatoric shear stress, indicating for the reference model a maximum at ~30 km depth of 40 MPa, reducing to nearly 0 below ~80 km depth. b) effective viscosity, which becomes lower than 10^20 Pa.s. in the reference model, to the difference of model d.c) Minimum (dashed lines) and maximum (plain lines) temperatures. Note that the great increase around ~90km depth corresonds to the detachement and elevation of light oceanic crust material throughout the continental mantle (see yellow bubble in figure6c). Figure8: PressureTemperature conditions obtained for the suggested lower seismic plane, located at the depth of the neutral plane in the subducting slab: a) Background colors show pressure, with subtraction of the lithostatic component (compression is negative and blue). b) Pressure temperature paths: green markers (dots) and green line (mean value ) show the neutral plane within the slab (labeled Plithos), indicating temperatures of 320~380°C at lithostatic pressure of 2.55 GPa. Pressure and temperature paths are also shown in blue and red, for areas above and below this neutral plane, which are at an overpressure and underpressure of 200 Mpa respectively, corresponding to the compressive stress (labeled Plithos+200MPa) and tensile stress (Plithos200MPa) of Figure8a. See appendix 4 for numerical sampling details. c)

Figure 9 :

 9 Reference model and 5 alternate models, showing pressure (the opposite of compression which is negative here, in blue) minus the lithostatic component, and isotherms along the slab. a) Reference model. b) Cold 150 Ma old oceanic plate, 300°C isotherm exceeds 160 km depth and innerslab compression is high, broadens and deepens. c) Warm continent, 400°C isotherm narrows and compression shallows to 100 km depth. d) High viscosity of the deep subduction channel, tension propagates into the forearc and innerslab compression vanishes at 90 km depth. e) Small slab curvature of 250 km, isotherms deepen and innerslab compression elongates from ~30 to 160 km depth. f) no applied upper plate convergence and no pre thickened crust, blurry compression appears from 80120 km depth.

  Figure 1

Figure 9

 9 Figure 9

  . In the case of Honshu and New Zealand, both subducting plates are much older than Northern Chile, but display equivalent surface heat fluxes capable to create the appropriate cool environment for segmented compressional stresses to develop. Furthermore, if we applyBrudzinsky et al.'s (2007) rule of thumb of DSZ separations equal to 0.14 km/Myr, a plate age of 50 Myr should give a separation of 7 km. If we use an effective plate age of 100 Myr more consistent with the observed surface heatflow in Northern Chile, then this separation should be 14 km, still lower than what is actually obtained in this study. Consequently, there must be another factor than plate age alone, that controls the distance between both seismic planes. The previous models indicate that indeed, the compressional state of stress of the upper plate transmits into the subducting plate, and increases the thickness and magnitude of its own compressional state, thus increasing the DSZ separation. Further global studies should help confirm this effect.

Table 1 :

 1 Vp, Vs values and Vp/Vs ratios in the Nazca Plate DSZ obtained in this study.

		Vp (km/s)	Vs (km/s)	Vp/Vs
	upper plane	7.7	4.6	~1.67
	region in between	8.5	4.5	~1.89
	lower plane	7.4	4.7	~1.57

Table 2 :

 2 Vp, Vs values and Vp/Vs ratios obtained by tomoDD in other DSZ. Honshu 39°N: Zhang et al. (2004), Honshu 37°N: Shelly et al. (2006), New Zealand 41°S: Zhang (2003).

		Region	Vp (km/s)	Vs (km/s)	Vp/Vs
	upper plane	Honshu 39°N Honshu 37°N	7.8 -8.0 8.5	4.4 -4.5	1.72 -1.85 1.82
		NewZealand 41°S	7.5		1.65
	region in between Honshu 39°N Honshu 37°N	8.2 -8.4	4.6	1.8 -1.85
		NewZealand 41°S	8.5		1.9
	lower plane	Honshu 39°N Honshu 37°N	7.6 -7.8 8.0	4.7 -4.9	1.6 -1.7 1.7
		NewZealand 41°S	7.0 -7.5		1.55

Table 3 .

 3 Notations and physical values common for all experiments[START_REF] Turcotte | Geodynamics: applications of continuum physics to geological problems[END_REF]. Bold values are those that were tested and displayed Figure8.

	Parameter P	Values and units Pa, MPa	Definition Stress, pressure
	ε µ* φ So		s 1 10¹ 0 10² 0 Pa s 1.520°Friction angle (see Table 4) Strain Effective viscosity 10 MPa Cohesion
	, G Vx λ Vb r		30 GPa 3 cm/yr 0 5 cm/yr 250450 km	λ Lamé elastic constants ( = G) Horizontal upper plate convergence Basal slabpull velocity Radius of curvature of oceanic plate
	hc hl T, To, Tb	3555 km 120180 km 10°1350°Temperature, surface, base Normal and prethickened continental Moho depth Thickness of lithosphere
	Cp Δ Τ		10 J/kg /°C 912°C/km	Specific heat Initial arc anomaly, with maximum at X=330530 km
	α kc		3×10 5 /°C 2.5 2.93.3W/m/°C	Thermal expansivity Thermal conductivity, crust, mantle, channel
	H hr		9.×10 10 W kg1 Radiogenic heat production, constant 10 km Radiogenic heat production exponential decay depth
	AC AO lc		400 Ma 60150 Ma 10 km	Thermotectonic age of the continental lithosphere Age of the oceanic lithopshere Thickness of subduction channel
	ke		25, 2000 m 2 /yr	Erosional diffusion coef., left & right of the Eastern Cordillera

Table 4 :

 4 Model parameters: density (ρ), compared to those chosen by[START_REF] Tassara | Threedimensional density model of the Nazca plate and the Andean continental margin[END_REF], dominant composition (labels given in Table5), friction (φ ), conductivity (k), for models layers. Value in bold was modified for model 8d.

		Oceanic	Oceanic	Asthenospher	Continental	Continental	Continental	Subducting	Subductio
		mantle	Crust	e	upper crust	lower crust	mantle	sediments	n channel
	ρ	3350	3050	3350	2750	3050	3340	2750	32503370
	ρ Tassara 32603500	2980	33103450	2700	3100	32403320	2980	29803550
	compos.	oldry	mfgr	oldry	qz	mfgr	olwet	plg	wetgrplg
	k	3.3	2.9	3.3	2.5	2.5	3.3	2.9	2.9
	φ	20°15°5°15°15°5°5°1.5°T
	able 5: Dislocation creep parameters (n power exponent, A material constant, Q activation energy), after references provided in Ranalli, 1995.
		Dry Olivine	Wet Olivine	Mafic granulite	Plagioclase	Quartz	Wet Granite
			oldry	olwet	mfgr	Plg	qz	wetgr
	n		3	2.5		4.2	3.2	2	1.9
	A (MPa n/s )	7.e4	3.e4	1.4e4	3.3e4	1.e3	2.e4
	Q (Jmol 1 )	5.2e5	4.44e5	4.45e5	2.38e5	1.67e5	1.37e5

Table 6 :

 6 Models tested parameters. Specific variation for each model is shown in bold (See Figure8, and Tables3, 4, 5 for more details).

		Oceanic plate age	Continental temperature	Creep parameters for channel	Plate curvature	Upper plate
	a) Reference model	60 Ma	hl=180 km + 9°/km	wet granite	450 km 3cm/yr + prethickened crust
	b) Old oceanic plate	160 Ma	hl=180 km + 9°/km	wet granite	450 km	3cm/yr + prethickened crust
	c) Warm continent	60 Ma	hl=120 km + 10°/km	wet granite	450 km	3cm/yr + prethickened crust
	d) High viscosity channel	60 Ma	hl=180 km + 9°/km	plagioclase	450 km	3cm/yr + prethickened crust
	e) Small plate curvature	60 Ma	hl=180 km + 9°/km	wet granite	250 km	3cm/yr + prethickened crust
	f) No upper plate compression	60 Ma	hl=180 km + 9°/km	wet granite	450 km	0cm/yr

(A) The square root of the derivative weight sum (DWS) for P waves and (B) for S waves along profile P2. Vp (C) and Vs (D) models obtained starting from the final 1D model (red dots, Figure 2). Vp (E) and Vs (F) models obtained starting from the initial model (blue dots, Figure 2). We see that the influence of the initial velocity model is minor.

Appendix 2: Three resolution tests performed to study the resolution of our inversion for Vp.

Checkerboard