Anne Paul

Ayoub Kaviani

Denis Hatzfeld

Jérôme Vergne

Mohammad Mokhtari

Seismological evidence for crustal-scale thrusting in the Zagros mountain belt (Iran)

Crustal receiver functions computed from the records of 45 temporary seismological stations installed on a 620-km long profile across central Zagros provide the first direct evidence for crustal thickening in this mountain belt. Due to a rather short 14-km average station spacing, the migrated section computed from radial receiver functions displays the Moho depth variations across the belt with good spatial resolution. From the coast of the Persian Gulf to 25 km southwest of the Main Zagros Thrust (MZT), the Moho is almost horizontal with slight depth variations around 45 km. Crustal thickness then increases abruptly to a maximum of ~70 km beneath the Sanandaj-Sirjan metamorphic zone, between 50 and 90 km northeast of the surface exposure of the MZT. Further northeast, the Moho depth decreases to ~42 km beneath the Urumieh-Dokhtar magmatic assemblage and the southern part of the Central Iranian micro-continent. The region of thickest crust is located ~75 km to the northeast of the Bouguer anomaly low at -220 mgals. Gravity modelling shows that the measured Moho depth variations can be reconciled with gravity observations by assuming that the crust of Zagros underthrusts the crust of central Iran along the MZT considered as a crustal-scale structure. This hypothesis is compatible with shortening estimates by balanced cross-sections of the Zagros folded belt, as well as with structural and petrological studies of the metamorphic Sanandaj-Sirjan zone.

Introduction

The Zagros mountain belt results from the closure of the Neotethys oceanic domain and the collision of the northern margin of the Arabian platform with the micro-plates of central Iran, accreted to the southern margin of Eurasia during the Mesozoic [e.g. [START_REF] Besse | Late Permian to Late Triassic palaeomagnetic data from Iran: constraints on the migration of the Iranian block through the Tethyan Ocean and initial destruction of Pangaea[END_REF].

Estimates for the age of the initial collision between Arabia and Eurasia along the Zagros suture vary between late Cretaceous [e. g. [START_REF] Berberian | Towards a paleogeography and tectonic evolution of Iran[END_REF] and Pliocene [see review in [START_REF] Allen | Late Cenozoic reorganization of the Arabia-Eurasia collision and the comparison of short-term and long-term deformation rates[END_REF]. The latest plate motion reconstruction by [START_REF] Mcquarrie | Cenozoic evolution of Neotethys and implications for the causes of plate motions[END_REF] based on updated maps of sea floor magnetic anomalies in the North and Central Atlantic and reconstructions across the Red Sea shows that ocean closure occurred no later than 10 Ma. Such recent ages lead some authors to consider the Zagros as the archetype of a continent-continent collision belt at an initial stage of its evolution. For example, the pioneering thermo-mechanical models of a continental collision zone by [START_REF] Bird | Thermal and mechanical models of continentcontinent convergence zones[END_REF] and [START_REF] Bird | Finite element modeling of lithosphere deformation: the Zagros collision orogeny[END_REF] relied on the Zagros case.

Studying the Zagros belt could thus be of great help to better understand the dynamics of mountain ranges. Nevertheless, published data on its crustal structure are very scarce, including Moho depth estimates. The only available profiles of crustal thickness variations have been computed from Bouguer anomaly modelling by [START_REF] Dehghani | The gravity field and crustal structure of Iran[END_REF] for whole Iran, and [START_REF] Snyder | Deep crustal structure and flexure of the Arabian plate beneath the Zagros collisional mountain belt as inferred from gravity observations[END_REF] for Zagros. Only two reports of direct Moho depth measurements exist in the Zagros. [START_REF] Giese | The crustal structure in Southern Iran derived from seismic explosion data[END_REF] provided questionable estimates of crustal thickness based on a non-reversed poor-quality refraction profile from quarry blast sources. [START_REF] Hatzfeld | Seismological constraints on the crustal structure beneath the Zagros mountain belt (Iran)[END_REF] estimated a crustal thickness of 46±2 km from receiver functions computed at a single station close to the town of Ghir in central Zagros. The numerous seismic reflection profiles recorded for oil exploration do not penetrate the crust beyond the thick Hormuz salt layer at 9-to-12 km depth [see e.g. [START_REF] Blanc | Structural styles in the Zagros simple folded zone, Iran[END_REF][START_REF] Mcquarrie | Crustal-scale geometry of the Zagros fold-thrust belt, Iran[END_REF][START_REF] Sherkati | Variation of structural style and basin evolution in the central Zagros (Izeh zone and Dezful Embayment), Iran[END_REF].

To improve our knowledge on crustal thickness variations across the Zagros, we deployed a temporary seismological network along a 620-km long profile between Busher on the coast of the Persian gulf and Posht-e-Badam, 160 km northeast of Yazd (Figure 1). The profile crosses all the morpho-tectonic units of the Zagros collision zone. In the southwest, it first cuts the Zagros fold-and-thrust belt (hereinafter referred to as ZFTB) well-known for its spectacular folds. The lower Cambrian to Pliocene continuous sedimentary sequence of the Arabian platform is strongly folded and thrusted with a basal decollement horizon in the lower Cambrian Hormuz salt at 9-to-12 km depth [START_REF] Stöcklin | Structural history and tectonics of Iran: A review[END_REF]. Its deformation is thus believed to be decoupled from the crustal shortening of the basement on distributed reverse faults [START_REF] Jackson | Errors in focal depth determination and the depth of seismicity in Iran and Turkey[END_REF][START_REF] Jackson | Basement faulting and the focal depths of the larger earthquakes in the Zagros mountains (Iran)[END_REF][START_REF] Berberian | Master blind thrust faults hidden under the Zagros folds: active basement tectonics and surface morphotectonics[END_REF]. The ZFTB is affected by frequent earthquakes of magnitude generally smaller than 7 concentrated at a depth of 8-15 km in the upper crystalline crust beneath the sedimentary sequence [START_REF] Talebian | A reappraisal of earthquake focal mechanisms and active shortening in the Zagros mountains of Iran[END_REF][START_REF] Tatar | Tectonis of the central Zagros (Iran) deduced from microearthquake seismicity[END_REF]. Northeast of the High Zagros, the profile cuts the Main Zagros Thrust (MZT) which is believed by most authors to be the suture zone between Arabia and Central Iran [e.g. Dewey and Grantz, 1973]. It separates the ZFTB from the Sanandaj-Sirjan zone (SSZ), a highly-deformed and moderately-metamorphosed remnant of the southern active margin of the Iranian continental block [e.g. [START_REF] Stöcklin | Structural history and tectonics of Iran: A review[END_REF][START_REF] Agard | Convergence history across Zagros (Iran): constraints from collisional and earlier deformation[END_REF]. The SSZ is bounded to the northeast by the Urumieh-Dokhtar magmatic assemblage (UDMA) characterized by almost continuous volcanic activity from Eocene to present [e.g. [START_REF] Berberian | Towards a paleogeography and tectonic evolution of Iran[END_REF]Berberian and Berberian, 1981], and believed to be the Andean-type arc related to the subduction of the Neotethys [START_REF] Berberian | Late Cretaceous and Early Miocene Andean-type plutonic activity in northern Makran and Central Iran[END_REF]. The 200-km-long northeasternmost segment of the profile crosses the southwestern part of the central Iranian microcontinent (CIMC) block, northeast of the town of Yazd.

Many important questions on the Zagros collision belt do not have a final answer yet.

The precise timing of the orogenic process and of the kinematics of the Arabia-Eurasia plate motion is still questioned [e.g. [START_REF] Mcquarrie | Cenozoic evolution of Neotethys and implications for the causes of plate motions[END_REF]Molinaro et al., 2005b]. Even though most authors agree on the location of the suture at the MZT, [START_REF] Alavi | Tectonics of the Zagros orogenic belt of Iran: new data and interpretations[END_REF] still proposes that it coincides with the northeastern boundary of the SSZ. Balanced and restored cross-sections through the fold-thrust belt suggest that cumulative northeast-southwest shortening lies between 25 km [Sherkati and Letouzey, 2005] and 49 km [START_REF] Blanc | Structural styles in the Zagros simple folded zone, Iran[END_REF] in the Dezful Embayment, and between 45 km [Molinaro et al., 2005b] and 69 km [START_REF] Mcquarrie | Crustal-scale geometry of the Zagros fold-thrust belt, Iran[END_REF] in the Fars arc. In the basement, shortening is accommodated at least partially by distributed thrusts as suggested by fault plane solutions [e.g. [START_REF] Tatar | Tectonis of the central Zagros (Iran) deduced from microearthquake seismicity[END_REF]. But the deformation processes that enable the lithosphere to accommodate the shortening at greater depths remain unknown. Do the basement thrusts connect at depth on a single decollement level? Is deformation in the lower crust continuous or localized? Does the Arabian lithosphere underthrust the Iranian microcontinent? A better understanding of the dynamics of the belt requires constraints on the lithospheric structure, including mapping of Moho depth which is the main goal of this work.

We use the receiver function technique which is an efficient tool to image structural discontinuities at crustal or mantle depths beneath seismological stations [START_REF] Vinnik | Detection of waves converted from P to SV in the mantle[END_REF]Langston, 1979]. The processing aims at enhancing P-to-S converted phases at velocity discontinuities in earthquake records at teleseismic distances.

Data selection and processing

From November 2000 to April 2001, we operated a network of 66 seismological stations along a profile trending N42 almost perpendicularly to the tectonic strike of Zagros. Our choice of this location relied both on scientific (most two-dimensional part of the range suitable for a transect) and logistical (easy and continuous road) considerations. The 45 stations whose locations are plotted in Figure 1 were equipped with broadband (Guralp CMG40-T and CMG3-ESP, Streckheisen STS-2) or Lennartz Le-3D-5s sensors, and continuously-recording Agecodagis Minititan digitizers. The average inter-station spacing was 14 km. We recorded 111 earthquakes of magnitude ≥ 5.5 at teleseismic distances between 25° and 98°. From their records, we computed radial and transverse receiver functions using the time domain iterative deconvolution method of [START_REF] Ligorria | Iterative deconvolution and receiver-function estimation[END_REF]. Radial receiver functions (RF) were selected on quality criteria including the arrival time of the most energetic pulse (that corresponds to the direct P at 0s in the absence of a thick low-velocity sedimentary cover), the amplitude ratio of the direct P to secondary arrivals, and the overall RMS of the RF. A time-section plot of radial RF stacked by station is shown in Figure 2. A move-out correction has been applied before stack to a constant slowness of 7 s.deg -1 using the same velocity model as in the depth migration (see details in the following section).

Traces are plotted at the abscissae of the stations, after projection onto the line trending N42 centred at location [30.6ºN; 53ºE] where it crosses the MZT (see Figure 1). Stacked radial RF exhibit waveforms of variable complexity, but a more-or-less energetic phase is visible all along the profile between 5 and 8s after the direct P wave (marked by arrows for selected traces in Figure 2). Since this phase has all the attributes of a P-to-S conversion at an interface with a rather strong velocity (and density) increase with depth (rather strong amplitudes, positive polarity) located at a depth between 40 and 70 km (time after P onset between 4.5 and 7.5s), with a good lateral continuity on many tens of km, we infer that it corresponds to the converted phase at the Moho (hereinafter named Ps). The spatial coverage is particularly weak under the SSZ at abscissae between 50 and 100 km due to a lack of stations in a basin with soft sand site conditions. However, the section documents clear lateral variations in the Ps-P time difference with the latest Ps being observed beneath the SSZ. The High Zagros region, between the Kazerun fault (KZF in Figures 1 and2) and the MZT, exhibits very strong and laterally coherent Ps converted phases. The radial RF at station B2 (km -66) plotted in Figure 3a (left panel) display this converted pulse of very strong amplitude. Figure 3b shows the radial and transverse RF at station M5 located at km 28 inside the SSZ.

Although the distribution in back-azimuth (center panel) is roughly the same at stations B2 and M5, the waveforms at M5 are less coherent than at B2. Nevertheless, the largest amplitudes pulses at M5 clearly arrive at later times than at B2. At both profile ends, southwest of the KZF and northeast of the SSZ, the Ps conversion is less clear and less laterally continuous, as documented by Figures 2 and3c (station H1, km 285). As a whole, the waveforms of the radial RF are much more complicated in the Zagros than in other regions where we have experience of similar recording conditions (e.g. in Tibet: [START_REF] Vergne | Seismic evidence for stepwise thickening of the crust across the NE Tibetan plateau[END_REF] or [START_REF] Wittlinger | Teleseismic imaging of subducting lithosphere and Moho offsets beneath Western Tibet[END_REF]). Moreover, signals of non-negligible amplitude remain after deconvolution on the transverse RF as documented by Figure 3 (right panels). The coherency of the transverse RF for different back-azimuths is an indication for departures from the basic hypothesis that crustal layers are horizontal and isotropic. Indeed, a close analysis of the transverse RF shows that there is no clear coherency from one station to its neighbours, indicating that possible interface dips and/or anisotropy vary laterally at inter-station scale.

However, locally dipping interfaces or anisotropy mostly change the amplitudes of converted signals on the radial RF. They have very small influences on the arrival times, which is the only information we are using in forthcoming analyses.

Migrated depth section

To image the crustal thickness variations, we performed a common-conversion-point (CCP) depth migration of the radial RF using the method described by [START_REF] Zhu | Crustal structure across the San Andreas Fault, southern California from teleseismic converted waves[END_REF]. It is based on the projection of the amplitude vector of each radial RF along the corresponding ray path computed in a particular velocity model. The volume beneath the network is binned in the three directions of space (with very wide bins in the direction perpendicular to the N42 profile) and the amplitude of each bin is obtained by averaging amplitudes of radial RF with ray paths crossing the bin.

The migrated image is highly dependent on the velocity model. Here, we used a modified version of the IASP91 standard earth model with a Moho at 75-km depth. This depth is chosen immediately above the maximum Moho depth observed after migration to ensure that Ps-P travel times are converted to Moho depths using crustal and not mantle velocities. For stations located southwest of the MZT, we used the 3-layer velocity model computed by [START_REF] Hatzfeld | Seismological constraints on the crustal structure beneath the Zagros mountain belt (Iran)[END_REF] from the inversion of arrival times of micro-earthquakes and radial RF in the Ghir region. It includes a 11-km-thick sedimentary layer with Vp=4.7 km.s -1 and Vp/Vs=1.77, a 9-km-thick crystalline crust with Vp=5.8 km.s -1 and a lower crust with Vp=6.5 km.s -1 , both with a normal Vp/Vs of 1.73. For stations northeast of the MZT where no information on crustal velocities is available, we assumed that there was no sedimentary layer and we increased the thickness of the upper crystalline crust to 20 km. The existence of a thick sedimentary sequence beneath the Zagros results in a Moho pull up of ~4 km after migration.

The Vp/Vs ratio has a critical influence on the time-to-depth conversion of Ps converted phases. It can be evaluated together with the crustal thickness from a grid-search stacking of the Ps and its first multiple, the PPs phase, at each station where the epicentral distance coverage is sufficient and the PPs is clear enough [START_REF] Zhu | Moho depth variation in southern California from teleseismic receiver functions[END_REF]. Since none of these two conditions is fulfilled in our dataset, it is impossible to estimate the Vp/Vs ratio in conjunction with the crustal thickness by this method. Another way of checking that the assumption of a laterally constant crustal Vp/Vs is not too crude is to stack the CCP migrated sections of the PPs and PSs multiples which are much more sensitive to Vp/Vs than the primary Ps [e.g. [START_REF] Wittlinger | Teleseismic imaging of subducting lithosphere and Moho offsets beneath Western Tibet[END_REF]. If the image obtained from the migrated multiples displays a converted phase at the same depth as the Moho in the Ps migrated section, the assumption of a laterally constant Vp/Vs of 1.73 is true. If the converted phases do not coincide, their depth difference can be used to measure the actual value of the crustal Vp/Vs ratio. We observe a converted phase in the section computed from the multiple phases only between km -170 and 0, that is in the part of the profile with strong primary Ps phases. In this region, the Moho depth obtained by migrating the multiples is close enough to the one obtained with the Ps to ensure that the assumption of a laterally constant Vp/Vs of 1.73 is reasonable. Unfortunately, no inference on the Vp/Vs ratio in other parts of the profile can be drawn from the multiples because they have too weak amplitudes.

The raw and smoothed CCP migrated depth sections are shown in Figure 4. The CCP migration significantly improves the signal-to-noise ratio and the spatial continuity of conversion interfaces. It is particularly useful for this dataset since the Ps pulse is hardly visible on a large number of individual radial RF. In the Moho depth range, the rather dense station spacing guarantees a sufficient multifold coverage (number of RF stacked in each bin) except beneath the SSZ (km 50 to 150) and at the south-western end of the profile. This is documented by the middle panel of Figure 4 where empty bins (crossed by no ray) filled in light grey are numerous beneath the SSZ. Furthermore, this plot shows that rays do not crisscross enough above ~30km depth with this station spacing to image the structure of the upper crust. The crust-mantle boundary appears as the only strong-amplitude laterally-continuous conversion interface. Its lateral continuity is enhanced by horizontal smoothing of the raw migrated section (bottom panel in Figure 4) using a Gaussian operator of half-width 10 km.

In the ZFTB, the image of the Moho changes at the Kazerun fault (KZF). Beneath the lower-elevation south-western Zagros, the radial RF display strong waveform changes between neighbour stations, resulting in a laterally discontinuous conversion interface in the depth section. This effect is reinforced by the rather poor coverage due to instrumental problems. However, the Moho can be picked with a slight north-eastward dip from a depth of 43±2 km beneath the coast of the Persian gulf (km -370, station F13) to 47±2 km beneath the region of Kazerun (km -155, station F5). A secondary pulse gives the red patch at ~20 km depth between km -240 and -180 in Figure 4. [START_REF] Hatzfeld | Seismological constraints on the crustal structure beneath the Zagros mountain belt (Iran)[END_REF] made a similar observation on radial RF computed at a single station in central Zagros and interpreted this pulse as a conversion on the top of the lower crust. Unfortunately, the station spacing is too large and the ray coverage is too poor in this depth range to precise the lateral extension of this interface. The High Zagros region is characterized by a more laterally-continuous Ps phase of stronger amplitude. The Moho depth seems to decrease slightly between 47±2 km and 43±2 km, 30 km southwest of the surface trace of the MZT. Unfortunately, good-quality RF are rare at stations further northeast due to bad site conditions. The maximum Moho depth of 69±2 km is mainly documented by a few RF recorded at stations S9 and S10. It is reached between 50 and 90 km northeast of the surface exposure of the MZT. As mentioned above, we checked by migrating the multiples that the hypothesis of a laterally constant Vp/Vs ratio of 1.73 is acceptable beneath the High Zagros.

However, due to the weak ray coverage beneath the SSZ where even observations of the primary converted pulse are scarce, this hypothesis cannot be checked in the region of thickened crust. Yet an underestimate of the crustal average Vp/Vs and an overestimate of the crustal average P-wave velocity would induce an artificial pull down of the Moho on the migrated depth section. Explaining the whole of the observed 25-km pull down would require either a very high Vp/Vs ratio larger than 2.1, which is impossible, or a P wave velocity of 5.2 km.s -1 combined with a Vp/Vs ratio of 2.0. Such a low value for the crustal average Vp is contradictory with petrological studies that document exposures of obducted ophiolite remnants in the MZT crush zone and metamorphosed phyllites and metavolcanics together with large-scale calc-alkaline plutons in the SSZ [e.g. [START_REF] Agard | Convergence history across Zagros (Iran): constraints from collisional and earlier deformation[END_REF]. Moreover, it is incompatible with the high Vp/Vs ratio of 2.0 which is even larger than the typical ratio of an oceanic crust (1.89 by [START_REF] Christensen | Poissons's ratio and crustal seismology[END_REF]). Therefore, the 25-km crustal thickening cannot be entirely an artefact of the migration. Still, part of it could be an artefact since an underestimate of 0.16 for Vp/Vs would result in an artificial pull-down of ~14 km. With an average Vp/Vs of 1.89, the crust beneath the SSZ would have the composition of an oceanic crust, quite far from the average continental crust of [START_REF] Christensen | Seismic velocity structure and composition of the continental crust: A global view[END_REF] [START_REF] Christensen | Poissons's ratio and crustal seismology[END_REF].

Although this composition is unrealistic, we consider that 0.16 is the upper bound for a possible underestimate of the average crustal Vp/Vs ratio beneath the SSZ in the migration, and that the observed maximum Moho depth could be overestimated by ~14 km at the maximum.

Most RF computed at stations in the UDMA (Urumieh-Dokhtar magmatic assemblage) and central Iran block display very complex waveforms with strong azimuthal variations and non-negligible energy on the transverse component (e.g. Figure 3c). Unfortunately, the azimuthal coverage is not sufficient to allow any analysis of the RF in terms of crustal anisotropy or dipping layers. We can only conclude on a complex crustal structure. However, good-quality records at stations S12 (km 140) display an unambiguous Ps pulse giving a Moho depth at 45±2 km. The precise geometry of the crustal thinning beneath the northeastern margin of the SSZ cannot be imaged due to the lack of ray coverage. At km 220, the Moho picked from the migrated depth section displays a 8-km Moho step which is particularly clear on the unsmoothed section (Figure 4). Further north-eastward, between km 270 and the end of the profile, the crustal thickness decreases gently from 50±2 km to 42±2 km. [START_REF] Dehghani | The gravity field and crustal structure of Iran[END_REF] computed a coarse Moho depth map of Iran from the conversion of the Bouguer anomaly map, and they determined a maximum Moho depth of 52 to 53 km right beneath the MZT along the profile considered in the present paper. A more precise modelling of the Bouguer anomaly data in the Zagros region has been conducted by [START_REF] Snyder | Deep crustal structure and flexure of the Arabian plate beneath the Zagros collisional mountain belt as inferred from gravity observations[END_REF]. Along a profile transverse to the range located ~50 km to the southeast of our profile, they proposed a crustal thickening from 40 km beneath the Persian Gulf to 65 km beneath the MZT. The deepest Moho is located in their model right beneath the surface trace of the MZT and the Bouguer anomaly minimum of -200 mgals. This disagrees with our results where the thickest crust is located beneath the SSZ, 50 to 90 km northeast of the exposure of the MZT. Our direct observation of the largest lapse times between the P and the Ps phases at stations in the SSZ is a strong indication for an thickening of the crust beneath the SSZ and not beneath the MZT. To reconcile seismological and gravity observations, we computed the Bouguer anomalies produced by a series of two-dimensional crustal models of the Zagros belt assuming the Moho depth profile picked from the CCP migrated depth section of Figure 4, and compared them to observations. To obtain the observed Bouguer anomaly along the N42 average profile used in the migration of the receiver functions, we first calculated a gravity map of Zagros by smoothing and interpolating the gravity data retrieved from the world database of the BGI (Bureau Gravimétrique International). The Bouguer anomaly map and the locations of the measured points are displayed in Figure 5. It shows that gravity measurements are dense enough in the study region to ensure the quality of the interpolation. The observed anomaly curve (thin continuous line in the top panels of Figures 6 and7) and its error bars (grey-shaded zone within the dotted lines) result from a sliding-window average of interpolated data points located in a 30-km-wide strip along the N42 profile (see location in Figure 5). We checked that the measured data points located within the same strip do fall within the error bars. The average topographic profile was computed in the same way from the GTOPO30 world digital elevation model of the U.S. Geological survey.

Gravity modelling

We used the method derived by [START_REF] Talwani | Rapid gravity computations for twodimensional bodies with application to the Mendocino submarine fracture zone[END_REF] In order to correct for the offset between the minima of the computed and the observed anomalies, we need shallower and/or thicker high-density material northeast of the MZT than southwest of it. A possible crustal model fulfilling this condition is shown in Figure 6b. It assumes that high-density rocks from the crystalline upper crust and the lower crust of the SSZ are raised by crustal-scale underthrusting of the southern margin of the central Iranian microcontinent by the northern margin of the Arabian platform which holds the Zagros belt. A similar hypothesis has been tested by Molinaro et al. [2005a] in a combined modelling of gravity, geoid and topography data along a profile across the Fars arc. Note however that the crustal thickening of their model was totally hypothetical. Strong geological evidences suggest that the Sanandaj-Sirjan zone overlaps the highly-deformed rocks of the so-called Crush zone in the MZT region, which in turn overthrust the Zagros folded belt [e.g. [START_REF] Stöcklin | Structural history and tectonics of Iran: A review[END_REF][START_REF] Ricou | Le Zagros[END_REF]. Moreover, [START_REF] Agard | Convergence history across Zagros (Iran): constraints from collisional and earlier deformation[END_REF] suggest that the amount of shortening, the deformation style and the lack of high-pressure metamorphic rocks along the past oceanic suture zone in the Crush zone indicate that the thrust branched onto the MZT at the surface is a major structure of crustal scale, possibly rooted to Moho depths. Very few other geological or geophysical data are available to constrain the density model. Most authors agree on a crystalline basement depth between 8 km and 12 km in the Zagros according to the fact that earthquakes interpreted as being related to reverse basement faults are concentrated in depth between 8 and 12 or 15 km [e.g. [START_REF] Berberian | Master blind thrust faults hidden under the Zagros folds: active basement tectonics and surface morphotectonics[END_REF][START_REF] Talebian | A reappraisal of earthquake focal mechanisms and active shortening in the Zagros mountains of Iran[END_REF][START_REF] Tatar | Tectonis of the central Zagros (Iran) deduced from microearthquake seismicity[END_REF]. We used the estimate by [START_REF] Hatzfeld | Seismological constraints on the crustal structure beneath the Zagros mountain belt (Iran)[END_REF] of an average lower crustal thickness of 25 km beneath Zagros. As no data on the crustal structure is available for the region northeast of the MZT, we assumed the depth to the basement to be 8 km and the thickness of the lower crust 25 km, exactly as beneath the Zagros. We know that this assumption is probably wrong and that the sediment layer is not that thick in the CIMC, if it exists. However, since we aim at comparing different hypotheses for the origin of crustal thickening to gravity observations, this assumption can be made for the sake of simplicity. All other geometric parameters of the final model were set by trial-and-error modelling of the observed Bouguer anomaly curve. We used the density values discussed by [START_REF] Snyder | Deep crustal structure and flexure of the Arabian plate beneath the Zagros collisional mountain belt as inferred from gravity observations[END_REF]: 2610 kg.m -3 for the sediments, 2800 kg.m -3 for the upper crystalline crust, 2950 kg.m -3 for the lower crust, and 3200 kg.m -3 for the upper mantle.

The top panel of Figure 6b shows that the crustal density model plotted in the bottom panel induces a Bouguer anomaly (thick continuous line) which satisfactorily fits the observations. We consider the fit as correct when the overall shapes of the observed and computed anomaly correspond, that is when the computed anomaly fits within the error bars of the observations. The crustal-scale overthrusting onto the MZT compensates for the greater Moho depth beneath the SSZ and shifts the anomaly minimum toward the MZT by increasing the lower crustal thickness and raising higher-density rocks. To improve the fit in the ZFTB, we assumed the existence of 3 basement faults in the vicinity of the MZT and beneath the KZF. However, as detailed modelling of the Bouguer anomaly is beyond the scope of this paper, these faults will not be discussed later. For the same reason, we did not try to improve the fit to the slightly more negative anomaly observed in the UDMA. This crustal model reconciles the Moho depth profile estimated from the RF with the gravity observations. Note however that the continuation of the thrust to the Moho is speculative since we have no direct observation of this structure at depth.

Discussion

On the origin of crustal thickening

The average crustal thickness beneath the ZFTB (measured from -300 to -25 km on the Moho depth profile of Figure 4) is 45±2 km, very close to the crustal thickness measured for the undeformed part of the Arabian platform, 45 km, from receiver functions at Ryad [START_REF] Sandvol | Lithospheric seismic velocity discontinuities beneath the Arabian Shield[END_REF][START_REF] Julia | Joint inversion of receiver function and surface wave dispersion observations[END_REF]. As discussed by [START_REF] Hatzfeld | Seismological constraints on the crustal structure beneath the Zagros mountain belt (Iran)[END_REF], the thickness of the crystalline crust beneath the Zagros (i.e. 35 ±2 km assuming 10 km of sediments) is about the same as the thickness of the thinned crust of the pre-collisional Arabian platform, assuming a stretching factor of 1.2 as estimated by [START_REF] Trowell | Lithospheric stretching, subsidence and magmatism in Oman and the Middle East[END_REF] from stratigraphic data along the Zagros margin. Therefore, we confirm on a much more complete dataset the conclusion discussed in [START_REF] Hatzfeld | Seismological constraints on the crustal structure beneath the Zagros mountain belt (Iran)[END_REF] that the crystalline basement of the ZFTB has not been thickened much by the collision yet.

The most striking feature of the measured Moho depth profile is the 25-km thickening of the crust in a rather narrow region (160 km) beneath the north-easternmost Zagros, the MZT and the Sanandaj-Sirjan metamorphic belt. The thickened crust is not located beneath the high elevations of the so-called High Zagros (between the MZT and the KZF, see the topographic profile in Figure 4), but in the back side of the Main Zagros thrust beneath the lower elevations of the SSZ. [START_REF] Alavi | Tectonics of the Zagros orogenic belt of Iran: new data and interpretations[END_REF] postulated that the suture is located between the SSZ and the UDMA and redefined the SSZ as a zone of thrust faults that have transported Phanerozoic units of the Arabian margin and obducted ophiolites to the southwest. We compared this hypothesis to our Moho depth profile and the Bouguer anomaly data in Figure 7a. Fitting the gravity observations requires that the basement layer reaches the surface in the SSZ, which sounds realistic, and that the lower crust is twice thicker beneath the SSZ than beneath the ZFTB. This requirement of a very thick lower crust is hardly compatible with Alavi's model which assumes that the crust beneath the SSZ is made of a stack of thrusted slices (see cross-section in Alavi [2004]). Moreover, such a lower crustal thickening requires a shortening of 50 to 100% in the SSZ, whereas it would be almost null beneath the ZFTB. Such a strong strain difference between two contiguous parts of the same continental margin sounds unrealistic.

Therefore, we think that our Moho depth model is not compatible with the hypothesis of [START_REF] Alavi | Tectonics of the Zagros orogenic belt of Iran: new data and interpretations[END_REF] and that the suture is located at the MZT.

In the following, the SSZ is considered as the south-westernmost part of the Iranian margin. Still, two hypotheses can be proposed to explain the thickening. In the first one, it is a remnant of an Andean-type crustal thickening of the Iranian margin related to the subduction of the Neo-Tethys oceanic domain. The crustal thickening of the SSZ would thus be due to volcanic intrusions, magmatic or metasediments underplating, and would predate the collision. Since the MZT is the southern border of the SSZ, this hypothesis implies that the suture is more or less vertical, as shown by Figure 7b (bottom). Here again, the lower crust beneath the SSZ has to be twice thicker than beneath the UDMA or CIMC to reconcile the Moho depth model and the Bouguer anomaly data (Figure 7b). This is quite different from the crustal structure of a typical Andean margin inferred from gravity modelling since in the Andes, the Bouguer anomaly is inversely correlated with the crustal thickness [e. g. [START_REF] Götze | Interpretation of gravity and geoid in the central Andes between 20° and 29°S[END_REF]. We think that this discrepancy argues against the hypothesis that the SSZ crust was thickened during the phase of oceanic subduction prior to the continental collision, as an Andean-type passive margin.

The second hypothesis explains the thickening by crustal-scale underthrusting of the active margin of the Iranian micro-block by the passive margin of the Arabian continent. The MZT considered as the thrust surface is rooted at Moho depth, as shown in the gravity model of Figure 6b. We have no direct evidence of this thrust cutting the whole crust in the migrated depth section of Figure 4. However, gravity modelling shows that the assumption of crustalscale thrusting of the margin of the Iranian microcontinent over the ZFTB is compatible with Bouguer anomaly data. Two geological arguments in favour of this hypothesis have been presented in the previous section: 1. there are clear indications that the SSZ overlaps the Crush zone which in turn overthrusts the ZFTB; 2. no high pressure metamorphic rocks are found along the suture [START_REF] Agard | Convergence history across Zagros (Iran): constraints from collisional and earlier deformation[END_REF]. Another positive element is the consistency between the total shortening measured in the Zagros fold-thrust belt by balanced and restored cross-sections (67 km by [START_REF] Mcquarrie | Crustal-scale geometry of the Zagros fold-thrust belt, Iran[END_REF] in the Fars, 49 km by [START_REF] Blanc | Structural styles in the Zagros simple folded zone, Iran[END_REF] and 25 km by Sherkati and Letouzey [2005] in the Dezful Embayment, and 45 km by Molinaro et al. [2005b] in the eastern arc of Fars north of Bandar Abbas) and the one associated to the crustal-scale thrust in the cross-section of Figure 6b (~30 km).

Present dynamics of the mountain belt

Although the Zagros belt as a whole is characterized by an important seismic activity, the MZT region in central Zagros as well as the SSZ are almost devoid of earthquakes [START_REF] Jackson | Active tectonics of the Alpine-Himalayan belt between western Turkey and Pakistan[END_REF]. Furthermore, no events at depth greater than 30 km have been reliably located in the ZFTB or northeast of the MZT [START_REF] Maggi | A re-assessment of focal depth distributions in southern Iran, the Tien Shan and northern India: do earthquakes really occur in the continental mantle?[END_REF][START_REF] Talebian | A reappraisal of earthquake focal mechanisms and active shortening in the Zagros mountains of Iran[END_REF]. These are clear indications that the crustal-scale thrust coincident with the MZT, if it exists, is seismically quiet which could mean that it is locked. Measurements of surface displacements by satellite geodesy suggest that the present shortening rate across the MZT in the Fars is almost negligible and that deformation at a rate of 8±2 mm.yr -1 is concentrated between the Persian gulf and High Zagros [START_REF] Tatar | The present-day deformation of the central Zagros from GPS measurements[END_REF]Walpersdorf et al, submitted]. However, since deformation of the Zagros sediment cover is decoupled from deformation in the basement by decollement levels including the Hormuz salt layer, the strain pattern at the surface is not a reliable indication of the strain pattern in the basement.

Anyhow, the hypothesized locking of the crustal thrust can be explained in the evolution models of a convergent plate-boundary proposed by [START_REF] Regard | From subduction to collision: Control of deep processes on the evolution of convergent plate boundary[END_REF] and [START_REF] Bird | Finite element modeling of lithosphere deformation: the Zagros collision orogeny[END_REF]. In their laboratory experiments of the closure of an oceanic basin and the transition to continental collision, [START_REF] Regard | From subduction to collision: Control of deep processes on the evolution of convergent plate boundary[END_REF] showed that oceanic subduction is followed by an episode of continental subduction when the passive margin of the "southern" continent (here: Arabia) underthrusts the active margin of the "northern" one (here: Iran). During that phase, surface shortening of the subducting lithosphere remains negligible and strain localizes at the trench. When the negative buoyancy of the subducted continental material is not compensated any more by downward driving forces such as slap pull, continental subduction stops and it is followed by collision. The onset of collision coincides with the onset of surface shortening of the "southern" continent and the creation of a fold-thrust belt, close to the margin first, and progressively extending further south. In the Zagros, the present tectonics is characterized by 1. a well-developed fold-thrust belt, 2. distributed shortening in the basement between the Persian Gulf and the High Zagros on blind thrust faults, and 3. the lack of deformation in the MZT region documented by seismic quiescence. These observations match the results of the experiments by [START_REF] Regard | From subduction to collision: Control of deep processes on the evolution of convergent plate boundary[END_REF] and suggest that the phase when the mega-thrust coincident with the suture absorbs all the convergence has ceased and that Zagros is currently undergoing distributed compressional deformation. If we could demonstrate that the megathrust also implies the lithospheric mantle, we would write that the episode of continental subduction has been replaced by continental collision. [START_REF] Bird | Finite element modeling of lithosphere deformation: the Zagros collision orogeny[END_REF] reached similar conclusions using two-dimensional finite element modelling of lithospheric deformation in the Zagros assuming pure anelastic rheology. He showed that shear stress concentrates on the suture zone for all models with slap-pull boundary conditions, whatever the densities and flow laws considered. To account for distributed deformation in the Zagros crustal wedge, the Neo-tethys oceanic slab must be detached there, the Arabia-Iran convergence now occurring in a compressive regime.

The question of slab break-off

Kadinsky-Cade and [START_REF] Kadinsky-Cade | Seismotectonics of southern Iran: the Oman line[END_REF] considered that the occurrence of a single intermediate-depth event at 107 km beneath the Urumieh-Dokhtar volcanic belt might indicate that the slab is still attached to the Arabian margin. However, the much bigger hypocenter database of [START_REF] Maggi | A re-assessment of focal depth distributions in southern Iran, the Tien Shan and northern India: do earthquakes really occur in the continental mantle?[END_REF] clearly shows that the few deep-focus events are related with the western edge of the Makran slab and have nothing to do with anything that was once under the Zagros. Hypocenters in the Zagros region are concentrated at depths shallower than 20 km, precluding both the presence of an oceanic slab subducting beneath central Zagros and active continental subduction at the MZT. Seismic wave tomographies of the mantle from body wave travel times [START_REF] Kaviani | La chaîne de collision continentale du Zagros (Iran): structure lithosphérique par analyse de données sismologiques[END_REF] and surface wave studies at regional scale [START_REF] Maggi | Surface waveform tomography of the Turkish-Iranian plateau[END_REF][START_REF] Bourova | Etude de la structure lithosphérique par l'analyse d'ondes de surface dans deux zones de convergence: la mer Egée et l'Iran[END_REF] have failed to image any high-velocity anomaly beneath UDMA and Central Iran. On the contrary, the upper mantle north of the suture is characterized by lower velocities than beneath Zagros, both for P and S waves.

However, the rather poor lateral resolution of surface wave studies [START_REF] Maggi | Surface waveform tomography of the Turkish-Iranian plateau[END_REF] and the inhomogeneous azimuthal coverage of teleseismic residual inversion [START_REF] Kaviani | La chaîne de collision continentale du Zagros (Iran): structure lithosphérique par analyse de données sismologiques[END_REF] preclude any final conclusion on oceanic slab break-off. [START_REF] Snyder | Deep crustal structure and flexure of the Arabian plate beneath the Zagros collisional mountain belt as inferred from gravity observations[END_REF] concluded their work on Bouguer anomaly modelling by a thorough discussion on the mechanisms responsible for the strong Moho bend they had to consider beneath the High Zagros to fit the gravity data. The Moho pull-down that we find beneath the MZT surface exposure is even stronger, with a maximum dip of 11 to 17º (corresponding to a minimum crustal thickening of 15 km and a maximum of 25 km within 80 km distance), while it is 5º in their model. As noted by [START_REF] Snyder | Deep crustal structure and flexure of the Arabian plate beneath the Zagros collisional mountain belt as inferred from gravity observations[END_REF], the topographic load in the Zagros is by far insufficient to explain such a Moho deepening within such a short distance by simple elastic flexure of the lithosphere. This is even more critical with our Moho depth model where the thickest crust does not correspond to the highest elevations as documented by the topographic profile of Figure 4. [START_REF] Snyder | Deep crustal structure and flexure of the Arabian plate beneath the Zagros collisional mountain belt as inferred from gravity observations[END_REF] proposed that isostatic and elastic flexure forces are acting together with hydraulic thickening of the plastic lower crust due to horizontal compression to produce the observed localized crustal thickening. The lower crust deforms plastically due to high temperature and pressure conditions. Since it is confined between the rigid upper crust and upper mantle, it responds to horizontal compression like an incompressible hydraulic fluid by thickening and bending the Moho down. Their hypothesis is compatible with our observations. In particular, it could explain why we have to consider a thickened lower crust beneath the MZT region in our density model of Figure 6b. However, this hypothesis still has to be checked by thermomechanical modelling of lithospheric deformation.

On flexural characteristics of the Arabian plate

In their elastic flexure modelling of the Bouguer anomaly in the Himalayas of Nepal, Lyon-Caen andMolnar [1983, 1985] showed that the observed increase in gravity gradient requires a steeper Moho (10-15º) beneath the High Himalayas. They interpreted this steepening as due to a lateral weakening of flexural rigidity of the equivalent elastic plate.

Later, [START_REF] Cattin | Gravity anomalies, crustal structure and thermo-mechanical support of the Himalaya of Central Nepal[END_REF] modelled denser gravity data along the same profile with a more realistic 2-D thermo-mechanical model taking petrological changes into account. They showed that as the topographic load increases, the lithosphere is flexed down and progressively weakened by strain and high temperatures of deep rocks brought up along the MHT (Main Himalayan Thrust). This flexural and temperature weakening explains why the rigidity of the Indian plate decreases locally beneath the High Himalayas, in agreement with the observed Moho steepening. A comparison between Zagros and Himalayas leads to the paradoxical conclusion that, although the topographic load in the Zagros is much weaker, the Arabian plate is as strongly bended beneath the MZT than the Indian plate beneath the High Himalayas. As in Nepal, the local steepening of the Moho could be partly due to a local weakening of the Arabian lithosphere, either predating the collision and/or related to strain and thermal weakening [START_REF] Cattin | Gravity anomalies, crustal structure and thermo-mechanical support of the Himalaya of Central Nepal[END_REF]. However, the question of the force driving the flexure remains.

The Moho step northeast of Yazd

The last point of this discussion concerns the Moho step imaged in the northeastern part of the profile, 30-35 km north of the town of Yazd. As already noted in a previous section, such vertical Moho steps have been imaged beneath the Tibetan plateau, often right beneath surface traces of major faults or sutures between accreted terranes. The Moho step discovered here does not correlate with any of the major faults or sutures mapped on classical geological maps of Iran. The Nain-Baft ophiolitic alignment which follows the boundary between the SSZ and the UDMA (see Figure 1) suggests that the Sanandaj-Sirjan microcontinent was separated from the Central Iran block by a narrow ocean basin in Late Cretaceous time [START_REF] Arvin | The petrogenesis and tectonic setting of lavas from the Baft ophiolitic melange, southwest of Kerman, Iran, Can[END_REF]], but our Moho step is located more than 100 km northeast of this suture.

Nonetheless, a possible hypothesis is that the step might correspond to an unmapped ancient suture between 2 micro-blocks of the CIMC.

Conclusion

From the records of a seismological network installed for 4.5 months across central Zagros, we have imaged Moho depth variations with a lateral resolution of a few km along a 620-km long profile transverse to the range. The average crustal thickness is 45 km beneath most of the Zagros fold-and-thrust belt and 42 km beneath the Urumieh-Dokhtar volcanic zone and the southern part of the Iranian microcontinent. The region of the suture from the High Zagros to the Sanandaj-Sirjan metamorphic zone is characterized by a marked crustal thickening from ~45 to ~70 km in a narrow region ~160-km wide centred 70 km to the northeast of the surface trace of the MZT. The very steep northeastward dip of the Moho beneath the High Zagros, the MZT Crush zone, and the southern part of the SSZ is defined unambiguously by strong Ps converted phases. However, due to the complex waveforms of the receiver functions, the crustal multiples are hardly visible even after common conversion point migration of the PPs and PSs multiples and stack of the migrated sections. This prevents from any estimate of the crustal average Vp/Vs ratio, in particular beneath the Sanandaj-Sirjan zone where the crustal thickness could be overestimated by ~14 km at the maximum if the actual value of the average crustal Vp/Vs ratio is 1.89 (typical of an oceanic crust) and not 1.73.

To reconcile the gravity and seismological data, we propose that the Zagros wedge underthrusts the SSZ along a crustal-scale fault linked at the surface to the MZT. The horizontal shortening of ~30 km implied by this thrust agrees with the shortening measured from balanced and restored cross-sections of the ZFTB. The existence of this crustal thrust would also explain petro-structural observations made in the Crush zone such as the lack of high-pressure metamorphic rocks along the suture.

Relying on the distribution of earthquake epicenters in central Zagros and on published experimental and numerical models of the transition from oceanic subduction to continental collision, we propose that the MZT crustal thrust is presently locked and that continental subduction has ceased. The shortening of 8±2 mm.yr -1 measured in central Zagros by satellite geodesy [Walpersdorf et al, submitted] would thus be entirely accommodated by distributed shortening of the crust. A possible explanation for this change in tectonic regime would be slab break-off. Indeed, none of the presently available tomographies of the upper mantle beneath Central Iran display the high velocity anomaly that should be associated with the subducting oceanic lithosphere. Their resolution can be questioned. Nevertheless, this lack of high-velocity anomaly in the mantle north of the MZT is particularly puzzling when we consider the unambiguous image of the same Neotethysian slab in the mantle beneath the active Makran subduction zone, ~700 east of our profile [Fig. 9a in [START_REF] Bijwaard | Closing the gap between regional and global travel time tomography[END_REF]].

High-resolution tomographies of the mantle beneath the SSZ and central Iran are needed to confirm that continental collision has followed the initial stage of continental subduction and slab-pull conditions in the Zagros belt. Figure 7. Results of 2-D gravity modelling for: (a) a crustal model with the suture located at the boundary between the SSZ and the UDMA [START_REF] Alavi | Tectonics of the Zagros orogenic belt of Iran: new data and interpretations[END_REF], (b) a crustal model with an Andean-type thickened margin beneath the SSZ. Same legend as Figure 6.

  From km -35, the Moho dip increases abruptly. As documented by Figure2, rather clear Ps pulses at stations P6, P8, P9, M4 and M5 are convincing evidences of the rapidly increasing Moho depth beneath the southern part of the Sanandaj-Sirjan metamorphic zone.

  Figure 6a, bottom). It displays the gravity signature of the crustal thickness variations, with a broad and strong gravity low in the SSZ and a much weaker anomaly (both in amplitude and scale) over the Moho step beneath central Iran. The computed anomaly does not fit the observations at all.

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. Time-section of radial receiver functions stacked by station after move-out correction to a ray parameter of 7s.deg -1 . Each trace is plotted at the abscissa of the projection of the corresponding station onto the N42 profile. Distances are measured relative to point [30.6N, 53.0E] where the N42 profile intersects the MZT. Arrows point out pulses interpreted as the Ps phase converted at the Moho for selected traces with best signal-to-noise ratio. Phases are picked on different criteria including amplitude (positive pulse of strong amplitude), time (last pulse of strong amplitude), and lateral continuity between neighbouring stations. The top panel shows the number of receiver functions stacked at each station. The bar in the middle panel shows the locations of the main faults and of the boundaries of the structural units intersected by the profile. ZFTB: Zagros foldand-thrust belt; SSZ: Sanandaj-Sirjan zone; UDMA: Urumieh-Dokhtar magmatic assemblage; CIMC: Central Iranian microcontinent; KZF: Kazerun fault; MZT: Main Zagros thrust.

Figure 3 .

 3 Figure 3. Examples of receiver functions computed at 3 stations located in different parts of the profile. Radial (left) and transverse (right) receiver functions are sorted by backazimuth. The back-azimuth values are shown in the middle panels on a linear scale. Note that most RF correspond to back-azimuths between 10 and 120º. The stack trace for both radial and transverse RF is plotted in the top panel. a) RF at station B2 located in the ZFTB; b) RF at station M5 located in the SSZ; c) RF at station H1 located in the CIMC.

Figure 4 .

 4 Figure 4. Migrated depth section computed from radial receiver functions along the N42profile. The blue-to-red color map displays the average amplitude ratio of the P-to-S converted phase to the primary P for all rays crossing the bin. Top: Average elevations along the N42 profile; Inverted triangles show heights of seismological stations. Middle: Raw migrated depth section. Empty bins (without a single ray) are plotted in light grey. The dotted line is the Moho depth profile picked from the smoothed depth section. Bottom: Migrated depth section after filtering and smoothing.

Figure 5 .

 5 Figure 5. Bouguer anomaly map showing locations of seismological stations (black triangles) and gravity stations (small black dost) of the BGI database. The Bouguer anomaly map shown as gray-shaded contours was obtained by interpolation and smoothing of the gravity measurements. The dashed lines delimit the 30-km wide strip used to compute the observed Bouguer anomaly curve modelled in Figures 6 and 7.

Figure 6 .

 6 Figure 6. Results of 2-D gravity modelling for: (a) a homogeneous crust model with thickness variations measured from the migrated depth section of Figure 4, (b) our preferred model with crustal-scale overthrusting of the SSZ on the ZFTB. Top: observed and computed Bouguer anomaly curves. The observed anomaly is plotted as a thin black line. The grey zone bounded by dotted lines has a width of 2 standard deviations (see text). The thick

  Figure 1

Acknowledgments:

We are grateful to all the people who made the seismological experiment successful. In particular, we thank O. Coutant, G. Herquel, G. Poupinet and the IIEES field teams for their assistance in the field work, and to the people of Zagros for their welcome. We acknowledge C. Péquegnat's work on the BDsis software for data preparation and management, which was used for the first time on this experiment. The interpretations presented here emerged after fruitful discussions inside the team of French geoscientists working in Iran, in particular P. Agard, J. Lavé, J. Martinod, M. Molinaro, and V. Regard. AK received partial support from the French Embassy in Tehran for his Ph-D thesis in Grenoble. The experiment was funded by IIEES, Iran, and INSU-CNRS "Intérieur de la Terre" programme, France. Seismological stations used in the experiment belong to the Lithoscope and RLBM CNRS-INSU national pools.


S Correyero Plaza

Julien Jarrige

Denis Packan

Renacimiento Barcelo

SEVILLE, SPAIN / Hotel

S Correyero
email: scorreye@uc3m.es.

E Ahedo

Ion acceleration in the magnetic nozzle of an ECR thruster: Comparison of experimental measurements with a quasi 1D kinetic model

Keywords: Magnetic nozzle, acceleration, Laser induced Fluorescence, Langmuir probe, kinetic modelling, ambipolar potential, magnetized expansion, electron cyclotron resonance

teaching and research institutions in France or abroad, or from public or private research centers.

Ion acceleration in the magnetic nozzle of an ECR thruster: Comparison of experimental measurements

with a quasi 1D kinetic model

INTRODUCTION

During the past years, several authors have contributed to the physical understanding of different mechanisms involved in the magnetic nozzle expansion, such as plasma detachment, formation of ambipolar electric field, magnetic thrust, etc [START_REF] Arefiev | Ambipolar acceleration of ions in a magnetic nozzle[END_REF][START_REF] Ahedo | Two-dimensional supersonic plasma acceleration in a magnetic nozzle[END_REF][START_REF] Longmier | Ambipolar ion acceleration in an expanding magnetic nozzle[END_REF][START_REF] Takahashi | Approaching the theoretical limit of diamagnetic-induced momentum in a rapidly diverging magnetic nozzle[END_REF]. However, hypotheses on electron dynamics have been found to be crucial for characterizing the ion velocity profile. For instance, isothermal electrons lead to unbounded ion acceleration, and polytropic laws have demonstrated a limited valid range. In this context, a kinetic steady-state quasi 1D model was developed in Ref. [START_REF] Martinez-Sanchez | Electron cooling and finite potential drop in a magnetized plasma expansion[END_REF] in order to determine self-consistently the axial evolution of the ambipolar electric potential and the electron and ion distribution functions of a fully magnetized plasma expansion.

Contrary to the commonly used polytropic laws, the model in [START_REF] Martinez-Sanchez | Electron cooling and finite potential drop in a magnetized plasma expansion[END_REF] does not assume an expansion coefficient, but instead accounts for magnetic mirror effects and potential barriers to determine the electron velocity distribution function. The model determines the finite potential drop in the magnetized expansion, which for cold ions, only depends on the ion to electron mass ratio. On the experimental side, ONERA has been developing an Electron Cyclotron Resonance (ECR) thruster during the past years composed of an ECR coaxial cavity followed by a divergent magnetic nozzle [START_REF] Jarrige | Characterization of a coaxial ECR plasma thruster[END_REF]. In order to investigate the ion acceleration and electron dynamics, an experimental set-up has been carried out at ONERA, where Langmuir probes have been used to characterize the evolution of the plasma potential, plasma density and electron temperature along the axis of the magnetic nozzle, and a Laser Induced Fluorescence (LIF) diagnostic has been employed to determine the mean ion velocity in the thruster exit plane and along the center line of the expansion [START_REF] Jarrige | Investigation on the ion velocity distribution in the magnetic nozzle of an ECR plasma thruster using LIF measurements[END_REF]. In the work presented here, the kinetic model in [START_REF] Martinez-Sanchez | Electron cooling and finite potential drop in a magnetized plasma expansion[END_REF] has been adapted to compare the solution with the experimental data along the magnetic nozzle obtained at ONERA. This paper is organized as follows: Section 2-explains the set-up developed at ONERA and presents the main experimental results. Section 3-describes the model used to predict the plasma properties along the expansion. Finally, section 4-compares the experimental data with the theoretical model and section 5-concludes the main results of the research.

EXPERIMENTAL SET-UP

The ECR thruster of ONERA

The experiments carried out in this work were performed along the magnetic nozzle of the permanent magnet version of ONERA's ECR thruster [START_REF] Jarrige | Characterization of a coaxial ECR plasma thruster[END_REF]. It has an axisymmetric geometry and it is composed of a coaxial plasma source cavity of dimensions L s = 15 mm and R s = 13.5 mm followed by a divergent magnetic nozzle. The inner conductor has a diameter of 1.7 mm, and The thruster is electrically isolated by an in-line microwave component called Direct Current Blocker (DC block), which allows the thruster to be electrically floating. The propellant is injected at the back plate of the source through two holes of 1 mm diameter located symmetrically at two mid-points of the source radius. An external magnetic field is applied by means of a permanent magnet. The microwave power is injected by means of a solid state amplifier from Kuhne Electronic. It is constituted by a signal generator (between 2.3 and 2.6 GHz) and a linear amplifier, and the maximum power it can deliver is around 100 W. The power is transmitted by a 50 Ohm coaxial cable, whose losses have been measured using a vector network analyzer. The output of the generator is connected directly to a circulator with a 50 Ohm load, which allows the reflected power to be absorbed and dissipated in the load. Subsequently and before entering the feedthrough, the power line contains a bidirectional coupler connected to two different diodes that are in charge of measuring the forward and reflected power at this point of the line. Figure 1 shows an schematic of the the microwave power line setup.

Diagnostics

In order to obtain longitudinal profiles of the ion acceleration and plasma properties along the magnetic nozzle expansion, two different diagnostics were installed: Laser Induced Fluorescence (LIF) optical set-up and cylindrical Langmuir probes. The LIF calibration and specific details of the optical setup are explained in [START_REF] Jarrige | Investigation on the ion velocity distribution in the magnetic nozzle of an ECR plasma thruster using LIF measurements[END_REF], as well as the main results from the LIF measurements presented in this work. The laser consists of a SDL-TC10 tunable laser diode with 20 mW of maximum output power and the frequency span for this measurements was 25 GHz, which allowed measuring axial ion velocities up to 20 km/s.

The laser was aligned with the thruster axis and injected through a lens of 25 mm diameter. The detection was composed of a 50 mm lens and a 200 µm optic fiber. The probed volume was around 1 mm diameter. The optics were installed 450 mm far from the measurement point, to avoid perturbation and intrusion in the plasma beam. The detection was installed with an angle of 35 o and a camera was fixed to the thruster to verify alignment between the laser and the detection. All optics were fixed to the vacuum chamber, and the thruster was moved using a 3D translation stage. With this setup, any point in the thruster source and along the plume could be measured. Figure 2 schematizes the set-up. A cylindrical Langmuir probe of 0.1 mm diameter and 6 mm length was placed in a complementary 2D translation stage, so it could operate at different longitudinal positions and be placed outside the optical line while performing the LIF acquisitions. The I-V characteristic curve was successfully obtained at different locations without perturbing significantly the thruster performance (from 75 to 250 mm from the thruster exit plane). The plasma potential was obtained from the inflection point of the I-V characteristic curve, and the electron energy distribution function (EEDF) was computed with the Druveysteyn formula (Equation 1 [START_REF] Druyvesteyn | The mechanism of electrical discharges in gases of low pressure[END_REF], where F ( ) represents the EEDF, g( ) the electron energy probability function (EEPF), m e and q the electron mass and charge respectively and A the probe surface area). An effective electron temperature and density can be computed by taking the corresponding integrals of the EEDF, as it is shown in Equations 2 and 3. These formulas give reliable results when a magentic field is present if the condition r L >> r s is satisfied, where r L represents the electron Larmor radius or gyroradius and r s the collecting area of the probe, taking into account the sheath expansion.

F ( ) = g( ) √ = 2 √ 2m e q 3 A d 2 I dV 2 (Eq. 1
)

n e = 2 √ 2m e |q|A -∞ 0 d 2 I dV 2 V |q| dV (Eq. 2
)

T e = 4 √ 2m e 3n e A -∞ 0 d 2 I dV 2 V |q| 3/2 dV (Eq. 3)

Vacuum facility

The tests have been done at ONERA facilities, located in Palaiseau, Paris. The vacuum chamber, known as B09, consists of a cylindrical vessel with dimensions of 0.8 m diameter and 2 m long. The pumping is done by one primary pump, three Pfeiffer Hipace turbomolecular pumps, and one cryogenic pump, allowing to achieve a background pressure of 10 -7 mbar and a Xenon pumping speed of 13000 l/s. Table 1 shows the measured chamber pressure when the thruster was operating at different mass flow rates.

Experimental results

This section presents the experimental results obtained at 30 W of absorbed power and 1, 1.5 and 2 sccm of Xenon. The coordinate z represents the distance to the thruster exit plane. The magnetic field intensity as a function of this distance is plotted in Figure 3. Measurements with the Langmuir probe were taken from 75 mm to 250 mm downwards. It was not possible to place the probe closer to the thruster exit, since it is a very intrusive diagnostic, and it was disturbing the thruster performance. The estimated plasma potential, electron density and electron temperature along the expansion are shown in Figure 4 for the three different mass flow rates. The mean axial velocity is shown in Figure 5 for each case. As it can be seen, measurements were taken from the thruster exit plane to 115 mm downwards. Finally, the plasma potential measured with the Langmuir probe can be overlapped with the estimated potential from the LIF measurements assuming total ion energy is conserved. Figure 6 shows the plasma potential profile obtained combining both techniques. Since there is an effective electron cooling along the expansion, an isothermal model would not be able to represent accurately the results. For this reason, and in order to reduce the number of assumptions and avoid fitting experimental polytropic coefficients, a paraxial model is proposed which takes into account the depletion of the electron energy distribution function.

QUASI 1D EXPANSION MODEL

The model presented here is a variation of Martinez-Sanchez, Navarro and Ahedo's kinetic model detailed in [START_REF] Martinez-Sanchez | Electron cooling and finite potential drop in a magnetized plasma expansion[END_REF]. Here, an externally applied magnetic field B(z) vanishing at B = +∞ and with a single maximum located at z = 0, generates a fully divergent magnetic nozzle that expands a plasma produced at z = 0. The model attempts to determine the quasi-1D, steady-state response of the plasma beam which is assumed collisionless.

In this approach, the dispersion in the ion velocity distribution function is neglected and ions are treated as a single cold species. Since the magnetic field is only divergent, all ions that leave the thruster would be accelerated downstream, so their motion is much simpler than for electrons. As well, ion total energy and ion flux are conserved along the expansion. This is shown in Equations 4 and5, where m i is the ion mass, u i is the ion velocity, e is the elementary charge, and φ is the ambipolar electric potential, which is monotonically decreasing.

E i = m i u 2 i 2 + eφ = m i u 2 i0 2 + eφ 0 (Eq. 4) Γ i B = n i u i B = n i0 u i0 B 0 (Eq. 5)
The electron population is modeled taking into account the effective potential barriers due to the magnetic mirror effect. The forward electron distribution f e+ is assumed Maxwellian upstream, at B 0 , where plasma starts expanding, following the normalization f e+ d 3 w = n e0 . f e+ is defined in Equation 1, where m e is the electron mass, E e is the electron energy, and n e * and T e * are reference parameters (the actual source density and electron temperature also involve the backward distribution). For each energy, there is a maximum magnetic moment µ M where w e = 0 at each location. Within this assumption, electrons from the plasma reservoir could belong to two different groups (free or reflected), while along the expansion there could be a population of electrons double-trapped in a region of the divergent nozzle. This population would be unconnected to the upstream source and it would be bouncing back and forward between the two spatial locations where their parallel energy vanishes. The double-trapped electrons could appear in the transient formation of the stationary flow [START_REF] S Ánchez-Arriaga | Kinetic features and non-stationary electron trapping in paraxial magnetic nozzles[END_REF] or due to infrequent electron collisions. We assume that the distribution function of these electrons has the same functional dependence than the one for electrons. (Eq. 9) The self-consistent solution for the electric potential comes from solving the quasineutrality equation (Eq.10) with the current-free condition (Eq.11).

f e+ (E e ) = n e * m e 2πT e *
n i (B(z), φ(z)) = n e (B(z), φ(z))
(Eq. 10)

n i u i = n e u e (Eq. 11)
After normalizing Eq.10 and 11, the solution of the electric potential φ(B) depends on two parameters i and M , defined in equations 12 and 13 , respectively. In this work, solutions have been obtained for initial ion velocities equal to T e * /m i , and for Xenon ions (M = 2.39 • 10 5 ). Figure 7 shows the solution φ(B) and the normalized plasma properties n/n 0 , v i / (T e * /m i ), and T e /T e * for the mentioned parameters. The total potential drop obtained is -e(φφ 0 )/T e * = 7.4.

COMPARISON OF EXPERIMENTAL RE-SULTS WITH THE MODEL

This section presents the first comparison between the proposed analytical model and the experimental results. The longitudinal ion velocity and the plasma potential at the thruster axis can both be compared with the model, since the experimental data starts from the thruster exit plane. Since the theoretical model proposed in this work is valid only for a collisionless supersonic expansion, the starting point corresponds to an ion velocity equals to T e * /m i , so it is therefore necessary to identify the sonic transition in the experimental data, which is referred as z = z s0 as follows. The re-gion of minimum area (or maximum magnetic field) is not necessarily the "sonic" transition, since this point could be shifted due to ionization and/or diffusion processes [START_REF] Burm | Mach numbers for gases and plasmas in a convergent-divergent cascaded arc[END_REF]. Indeed, if the electron distribution function is non-isotropic already at the thruster exit, the temperature associated to the sonic transition is only the parallel component, so the results presented here should be corrected.

The reference temperature T e * has been defined by fitting the experimental data with the quasi-1D model. This temperature can be seen in Table 3, with its corresponding sonic point in the expansion. Figure 8 shows the comparison between the analytical model and the experimental data in terms of the mean ion velocity. The magnetic field has been normalized with its corresponding value at the sonic point B 0 = B s0 . As it can be seen, the three data sets normalized with their corresponding sound velocity T e * /m i follow the same tendency downwards of the "sonic" point, which is well represented by the model, but they follow different curves in the subsonic region. Figure 9 shows the measured plasma potential normalized with T e * where the red line represents the model results. Again, if the plasma potential is normalized from the sonic point with a reference temperature T e * , they follow the same curve independently of the mass flow rate. This is a very important result, since it reinforces the practicality of dimensionless 1D models. However, it is clear that the model should be improved to include ionization and diffusion processes, as well as anisotropy at the plasma reservoir, in order to properly identify the sonic transition in the expansion and avoid the calculation of a reference electron temperature T e * .

Xenon (mg

CONCLUSIONS

This paper presents experimental measurements of ion acceleration and plasma properties along the magnetic nozzle of the ECR thruster developed by ONERA. Langmuir probes and Laser Induced Fluorescence have been installed at ONERA facilities to determine the mean ion velocity, the plasma potential drop, electron density and electron temperature along the axis of the magnetic nozzle. The operating conditions were 30 W of absorbed power and 1, 1.5 and 2 sccm of Xenon. The ion velocity profile and plasma potential have been normalized with a reference temperature T e * and compared with a steady-state quasi 1D supersonic kinetic model. It has been found that the ion velocity and plasma potential profiles for different mass flow rates, follow the same curve when they are normalized and shifted to the sonic transition point. A quasi 1D steady state model has been developed for a fully divergent magnetic field, assuming quasinuetrality and current-free expansion. The ion population is treated as a cold fluid species, while the electron distribution function is solved self-consistently accounting for magnetic mirror effects and potential barriers. Electrons are defined Maxwellian in the forward direction at the starting point of the expansion, while the complete population accounts for the backward and forward distribution, and allows regions of double-trapped electrons. A unique solution for the ambipolar potential profile has been found, which depends only on the ion to electron mass ratio and a parameter relating the ion energy with a reference temperature T e * . A total potential drop of 7.4 T e * has been calculated for Xenon. Finally, the ion velocity and plasma potential profiles have been compared with the quasi-1D model. When normalizing the results, a reference tempera-ture T e * of 45 eV, 27 eV and 19 eV has been found for 1, 1.5 and 2 sccm, respectively. Despite the simple polytropic laws, the model presented here does not need to assume an expansion coefficient, and only depends on the ion to electron mass ratio. As a limitation when comparing with experimental results, the need of determining a reference electron temperature to define the sonic transition highlights the importance of obtaining accurate experimental data. As a next step in this research, the model must be improved to determine consistently the sonic transition of the expansion, and justify the normalization presented in this work. A first modification would be to incorporate electron temperature anisotropy in the plasma source, this is, modifying the definition of the electron distribution at the plasma reservoir. A more complete simulation should include collisions and diffusion processes in the near plume region.

Figure 1 :

 1 Figure 1: Schematic of the thruster and subsystems in the vacuum setup

Figure 2 :

 2 Figure 2: Schematic of optical LIF setup inside the vacuum chamber

Figure 3 :

 3 Figure 3: Magnetic field intensity along the thruster center line. z = 0 corresponds to the thruster exit plane.

Figure 4 :

 4 Figure 4: Longitudinal plasma properties: Plasma potential, electron density and electron temperature obtained with the Langmuir probe.

Figure 5 :

 5 Figure 5: Mean axial ion velocity in the thruster axis obtained with the LIF optical set-up

Figure 6 :

 6 Figure 6: Plasma potential obtained with the Langmuir probe and LIF measurements assuming total ion energy is conserved.

Figure 7 :

 7 Figure 7: Model results for xenon ions and initial ion velocity equal to T e * /m i

Figure 8 :

 8 Figure 8: Normalized ion velocity from LIF measurements with T e * /m i and validity with the analytical model (red).

Figure 9 :

 9 Figure 9: Normalized electric potential with the characteristic estimated T e * and comparison with analytical model (red).

Xenon mass flow rate (mg/s) Chamber pressure (mbar)

	0.1 0.15 0.2	2.8 • 10 -6 3.5 • 10 -6 3.7 • 10 -6

Table 1 :

 1 Measured chamber pressure at three different xenon mass flow rate. The absorbed power was 30 W for all cases.

Table 2 :

 2 Table 2 summarizes the main assumptions and characteristics of the model. Main assumptions and characteristics of the model

		Quasi 1D model
		Fully divergent
		Collisionless
		Quasineutral
		Current-free
		Electrons fully magnetized
		Cold fluid ions
		f e+ Maxwellian at z = 0
		Monotonic plasma potential
	Magnetic mirror effect and potential barriers
	Populated regions of double-trapped electrons
	Macroscopic variables are computed taking integral
	moments of the distribution function according to
	Equation 9, where the boundaries on the integra-
	tion domain are set according to the type of electron
	(free, reflected, double-trapped).
	χ e =	2πB (z) m 2

e χ e (E e , µ e ) f e+ (E e , µ e ) |w e | dE e dµ e

Table 3 :

 3 Estimated reference electron temperature at the sonic point and its location from the thruster exit plane

		/s) T e * (eV) z s0 (mm)
	0.1	45	7.4
	0.15	27	11.8
	0.2	19	12.5

ACKNOWLEDGMENTS

This work was made in the framework of project MINOTOR that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 730028. This work was additionally supported by the Spanish R D National Plan (Grant No. ESP2013-41052-P).