The Southern Zagros Collisional Orogen: New Insights From Transdimensional Trees Inversion of Seismic Noise

Imaging and resolving the lateral continuity of 3-D crustal structures enhance our ability to interpret seismicity and to understand how orogens are created. We apply a Bayesian, hierarchical inversion approach based on a transdimensional trees-structured wavelet parameterization to recover phase velocity maps at 2–40 s periods. We then invert phase velocity dispersion to constrain a 3-D shear velocity model of the crust beneath south-central Iran. Together with accurate earthquake centroid depths and focal mechanisms, the pattern of 3-D velocity variations supports recent suggestions that most large earthquakes in the Zagros occur within the lower sedimentary cover or close to the sediment-basement interface. Furthermore, we find evidence for Arabian basement underthrusting beneath central Iran, although only in one location does it appear to generate earthquakes. Our new 3-D tomographic model clarifies and throws new light on the crustal structure of the SE Zagros and its relation to seismicity and active faulting.


Introduction
Seismically active orogenic belts, such as the Zagros fold and thrust belt (ZFTB, Figure 1), result from shortening continental crust. The ZFTB is well known for its high earthquake activity and involves thrusting and folding of a thick Cambrian-to-Quaternary sedimentary sequence of the former Arabian continental margin, separated from central Iran by the suture of the Main Zagros Thrust Fault (MZTF) (e.g., Berberian, 1995). In the eastern Fars arc, the Infracambrian Hormuz salt, a layer with maximum thickness likely to exceed 1 km in places (Edgell, 1996), decouples the sedimentary cover from its underlying basement, allowing the style of deformation between the two to be different or spatially separated. Beneath the sedimentary cover of the Zagros lies the Precambrian eastern platform of the Arabian plate, which is far less understood than its counterpart to the west (in the Arabian shield), due to lack of exposure.
Earthquake focal mechanisms show that the dips of seismically active thrust faults in the Zagros are typically between 30 • and 65 • , perhaps inherited from the reactivation of normal faults formed during the Mesozoic rifting of the Arabian continental margin (Jackson, 1980). From detailed InSAR and seismological studies it has become clear that much of the seismicity occurs in the upper Arabian basement or in the lower part of the overlying sedimentary cover (e.g., Nissen et al., 2011Nissen et al., , 2014. Coseismic surface faulting is extremely rare in the Zagros, presumably because of the numerous decoupling evaporite horizons within the sedimentary section. Much of our understanding of the current deformation has come from accurate focal mechanisms and centroid depths for M w > 5 earthquakes determined from long-period teleseismic P and SH body waves ( Figure 1 and supporting information Figure S1). There are now more than 80 earthquakes analyzed in this way, particularly in the Fars area (as published in; Maggi et al., 2000;Nissen et al., 2011;Talebian & Jackson, 2004). Most of these earthquakes in the ZFTB are confined to the upper crust, rupturing at depths less than 20 km, and there is no evidence of a seismically active subducting slab beneath central Iran. However, structural imaging and insights in the Zagros are limited and localized geographically to a few 1-D (from receiver function and dense microearthquake networks; e.g., Gholamzadeh et al., 2009;Hatzfeld et al., 2003;Nissen et al., 2011;2014) and 2-D models (from balanced cross sections or seismic reflection lines; e.g., Jahani et al., 2009;McQuarrie, 2004;Molinaro et al., 2005;Vergés et al., 2011) or depth-to-magnetic basement maps Figure 1. Simplified tectonic map of the northeastern corner of the Arabian plate as modified from the National Geoscience Database of Iran. Black lines are major faults; in the Zagros, these are usually inferred from abrupt changes in the stratigraphic level within the Mesozoic sedimentary section. Red line is the suture boundary (MZTF: Main Zagros Thrust Fault) between the Arabian platform (to the south) and central Iran (to the north). Colored dots are teleseismically determined earthquake centroid depths, in the Zagros region only. Seismic stations used in this study are depicted as white squares. Green squares denote the stations used to measure the EGFs shown in Figure S2. White lines show the location of vertical tomographic sections presented in Figure 3. Thick black line shows the location of the receiver function transect from Paul et al. (2006). Black arrows are GPS velocities, for some representative locations, relative to Eurasia-fixed frame (Khorrami et al., 2019). The inset in the upper right corner shows the location of the study area. Background map represents elevation, black is high, and white is low. MKS: Makran Subduction front; MZPF: Minab-Zendan-Palami Fault; MP: Musandam Peninsula. (Teknik & Ghods, 2017). In this study, we contribute 3-D images and discuss their compatibility with the earthquake, structural, and seismic reflection data.
We exploit recently acquired data from a passive seismometer network and recent developments in transdimensional tomography to produce a 3-D shear velocity model of the crust beneath the ZFTB. Our tomographic model is an advance on previous geophysical modeling efforts in southeastern Iran and introduces new information on the primary structure of the crust. In particular, our 3-D depth-to-basement interpretations are used to examine the distribution of large earthquakes in the sedimentary layer and the underlying Precambrian basement, with the 3-D pattern able to show the geographical continuity of structure and seismicity that complements 1-D and 2-D investigations. Moreover, we use the velocity variations in our tomographic model to test the hypothesis of Arabian basement underthrusting beneath Iran, as has been claimed from the presence of low-angle (dip <15 • ) thrusting earthquakes north of the MZTF in the far southeastern Zagros (Nissen et al., 2011;Talebian & Jackson, 2004).

Data and Methods
Continuous seismic records used in this study were recorded between 2014 and 2016 from multiple networks. These include 23 permanent stations maintained by the National Center of Meteorology of the United Arab Emirates, 31 temporary stations run by the Petroleum Institute (Abu Dhabi), and 23 broadband stations of the Iranian Seismic Network (Figure 1). Empirical Green's functions were extracted from the long-term cross correlation and stacking of vertical component data. To increase the signal-to-noise ratio, we stacked the daily cross correlations for every station pair using a nonlinear phase-weighted procedure (Pilia et al., 2016; Schimmel et al., 2011;Ventosa et al., 2017), which is insensitive to the amplitude content of the seismic signal and reduces sensitivity to incoherent noise (see Text S1 and Figure S2). Rayleigh wave phase velocity dispersion measurements were carried out using the image transformation technique of Yao et al. (2006) as modified by Young et al. (2011). Fundamental mode path-averaged phase velocities were obtained for periods between 2 and 40 s (see Text S1 and Figure S3).
Lateral variations of period-dependent phase velocity maps were obtained using a novel inversion technique based on transdimensional sampling over tree structures of variable complexity, developed by Hawkins and Sambridge (2015). The method initially predicts ray path trajectories via a robust grid-based Eikonal solver known as the fast-marching method (Rawlinson & Sambridge, 2004;Sethian, 1996). Subsequently, it employs a wavelet parameterization of the velocity field with a hierarchy of wavelet coefficients constructing models that range from coarse to fine scale, in which seven different levels of perturbation are allowed ( Figure S4). A key feature of the algorithm is that trees are dynamic, since the number and distribution of model parameters are implicitly controlled by the data. Hence, the inversion is stabilized through information in the data without the need for explicit regularization. The ensemble of velocity models is generated by sampling the model space using the reversible-jump Markov chain Monte Carlo (rj-McMC) sampler, following the formalism of Green (1995). At each iteration, we can either add one or more nodes to a tree or remove one or more nodes from a tree or perturb one or more node values within an existing tree, making the inverse problem a transdimensional one as the number of wavelet coefficients is allowed to vary throughout the inversion. The posterior probability density function of any model created within a Markov chain is subsequently evaluated, and the model can be either accepted or rejected. Parallel tempering (Earl & Deem, 2005;Sambridge, 2014), a parallel interacting chain approach, was used to aid convergence of the Markov chains by more effectively sampling the model space and overcoming the potential issue of local modes and multimodal posterior distributions. Both Laplacian prior and noise in the data are treated as hierarchical parameters, hence considered unknown during the inversion (Hawkins & Sambridge, 2015;Hawkins et al., 2018). The recovery of a phase velocity map for a single period takes approximately one million iterations (500,000 discarded as "burn-in" period) computed in 1.5 hr using 320 CPUs, distributed in eight parallel Markov chains with five temperatures and eight cores per chain. Although the inversion scheme is computationally expensive, the improved quality of the results and provision of uncertainty estimates fully justifies the additional cost when tomographic inversion results are compared to those obtained with a conventional, regularized inversion method that implements an arbitrary parameterization (see Text S2 and Figures S5 to S10). The final phase velocity maps are obtained by computing the mean of the ensemble of models from varying parameterization across all chains.  Figure 1 for location). White circles represent teleseismic centroid depths within ±35 km distance from their respective tomographic cross sections. Errors on the centroid depths are typically in the order of ±3 km (Nissen et al., 2011), roughly comparable to the size of the circles. The black dashed line in profile BB ′ delineates the Moho discontinuity as derived by Paul et al. (2006) from a parallel receiver function profile about 200 km to the NW (see Figure 1 for location of the profile); the black dash-dotted line in profile BB ′ outlines a velocity feature that may be associated with the Arabia basement (see Section 4). Vertical exaggeration is 1:2 for all profiles but BB ′ , which has a 1:1 aspect ratio. DF: Deformation Front; MZTF: Main Zagros Thrust Fault; SSZ: Sanandaj-Sirjan metamorphic zone; MZPF: Minab-Zendan-Palami Fault.
The phase velocity field is subsequently discretized into 1,186 evenly spread grid points from which Rayleigh wave phase velocities are extracted from each period-dependent map to generate a pseudo 1-D phase velocity dispersion curve. To invert for 1-D shear wave speed profiles, we adopt the rj-McMC and assume that the number, position, and velocity of the layers are all unknown during the inversion; the posterior is then a transdimensional function and can be sampled with the rj-McMC algorithm (Bodin et al., 2012;Crowder et al., 2019;Pilia et al., 2014Pilia et al., , 2015. More details of the methods are given in the supporting information.

Tomographic Model
Convergence of the algorithms was monitored by qualitatively comparing the trends of independent chains for a set of indicators (i.e., likelihood, hierarchical, prior, and number of wavelet coefficients) to ensure that  (Nissen et al., 2011). (b) Schematic cross section along profile DD ′ illustrating the relationship between earthquake distribution and main tectonic structure as determined by this study. The boundary between the Arabian basement and central Iran is unknown and it does not appear as a velocity contrast in the tomographic model. (c) Depth distribution over the ZFTB of teleseismic centroid depths compared to the base of sediments depth as inferred from this study. Our results confirm that there is a tendency for medium-to-large earthquakes to occur in the lower part of the sedimentary cover. (d) Zoom-in of the eastern ZFTB where earthquakes occur beyond the MZTF with a low-angle north dipping thrusting mechanism, as suggested by the fault plane solutions. Background is topography. Numbers are centroid depths in km. they sample about the same region of the model space (further explanation is provided in Text S2, together with an example for a 2-D map at 20 s in Figure S10). To assess the quality of our inversion results, we separately perform resolution tests based on synthetic data for 2-D and 1-D velocity structure by using identical source-receiver combinations to the observational data set (see Text S3 and Figures S11 and S12). The spatial resolution on the 2-D maps was evaluated through a predetermined checkerboard input, which involves perturbations of different size and strong velocity gradients. An intrinsic property of the transdimensional parameterization applied here is that spatial resolution can vary throughout the model, as it adapts to the information content of the data. Along with mean velocity information that can be extracted from the ensemble of models, one can also calculate uncertainty associated with structure at the same spatial resolution (i.e., an error map) by calculating the standard deviation of the solution ensemble, as shown in Figure 2a for a map at 10 s period. The lowest standard deviation values correspond to the better resolved areas, which ultimately match areas where ray path coverage is more dense. Although it is difficult to precisely determine the exact spatial resolution in each area of the model, we estimate that we can recover structure as small as 50 × 50 km 2 in areas of low standard deviation (less than 0.1 km/s) and up to 150 × 150 km 2 in others.
In order to investigate the reliability of the 1-D inversion, we assume a known shear velocity model (as that found along profile DD ′ near the Strait of Hormuz) and calculate the corresponding phase velocity 10.1029/2019GL086258 curve (Figures 2b and 2c). The latter is inverted for a shear velocity model and the results compared to the initial profile. The exact depth of each discontinuity is well retrieved, but we anticipate that layers thinner than 5 km cannot be accurately resolved by our data set. However, it is important to note that the posterior probability of having a discontinuity is greatest at approximately the same depths as the discontinuities contained in the synthetic model (Figure 2c). The results of our synthetic tests show that we can recover both pattern and amplitude of velocity anomalies with confidence, except perhaps in the peripheral areas of the network where smearing in the 2-D maps can hinder the construction of reliable dispersion curves that are subsequently inverted for shear velocity structures. Figure 3 illustrates shear velocities via a set of cross sections representing the mean of the ensemble of models produced by the transdimensional, Bayesian, hierarchical inversion procedures (see Section 2). The Arabian domain (including the ZFTB), particularly in the foreland area of the continental collision, appears to be dominated by a broadly stratified velocity structure. The model shows a thick layer of relatively low velocities in the uppermost crust, in regions where Figure 2a suggests the model is well resolved. Between 8 and 13 km depth, the tendency is to move to faster velocities (increase of up to 1 km/s), with the boundary, or steep gradient, presumably indicating a depth limit on the thickness of Cambrian-Quaternary sedimentary formations of the Zagros. Elevated crustal velocities extend for about 8 km deeper, after which they revert to a pattern of slightly lower velocities that persists for nearly 15 km. Of particular relevance to further discussion is the observation of a localized area, centered on 28 • N/57 • E (see section DD ′ on Figure 3), where earthquakes occur NE of the MZTF and extend to depths of about 30 km, with focal mechanisms consistent with low-angle thrusting of the Arabian platform beneath central Iran (Figure 4). A distinct change in the general pattern of velocities occurs as the profiles transition from the Arabian platform to central Iran through the MZTF.

Discussion
Previous 3-D tomographic models of the Zagros have been coarse (resolution ∼250 km horizontally and ∼20 km vertically) and principally focused on the lithospheric mantle, illuminating a thick lithospheric root beneath the ZFTB, compatible with accommodation of continental shortening by thickening (Motaghi et al., 2015;Priestley et al., 2012). Similarly, thickening has been documented in the crust (Paul et al., 2006;2010) through a number of 2-D receiver function transects across the Zagros mountain belt, showing that the Moho discontinuity deepens by about 20 km beneath the Sanandaj-Sirjan metamorphic zone to 65 km, from ∼45 km on either side. These studies offer few constraints on 3-D or upper-midcrustal structure. The only previous 3-D crustal study of the area is that of Mottaghi et al. (2013), who applied ambient noise surface wave tomography. However, due to the sparse station coverage, lack of shorter-period data (shortest period used was 8 s), and the use of an inversion scheme based on regular parameterization, the resulting tomographic images are not of comparable resolution to ours. Our 3-D shear wave model, obtained by exploiting a new seismic network in the Arabian peninsula and recent developments in robust transdimensional, Bayesian-driven tomography, is an improvement, especially in estimating variations in depth-to-basement, and reveals some significant tectonic features of this young active continental collision. We examine our results below, with the aid of a robust earthquake catalog that features well-constrained depths (uncertainties are typically ±3 km) and focal mechanisms. By contrast, it is worthwhile noting that the EHB catalog earthquake depths, until recently the most accurate teleseismic compilation for the larger events, are generally characterized by 10-15 km errors (Engdahl et al., 2006). Newly calibrated earthquake relocations by Karasözen et al. (2019) are a significant improvement, though uncertainties on focal depths are larger than those obtained with body waveform modeling.
Using data from a dense temporary seismic network, Hatzfeld et al. (2003) investigated 1-D crustal and upper-mantle velocity variations beneath the Ghir region (centered at about 28 • N/53 • E and close to profile BB ′ in Figure 3) through a microearthquake survey and receiver functions. They inferred that the ∼11 km thick sediments are underlain by an 8 km thick upper crystalline crust and that the lower crust, which shows unusually low P wave speeds on their receiver function analysis, extends from 19 km depth to the bottom of the crust, found to be at 46 ± 2 km. Despite the inherent inability of surface waves (i.e., phase velocity in our study) to localize sharp discontinuities (as shown in Figures 2 and S12), the results produced here fundamentally concur with the 1-D profile derived by Hatzfeld et al. (2003).
The uppermost part of the crust in the ZFTB is clearly dominated by low S-wave velocities, ranging ∼2.5-3.3 km/s. The shallowest part of the white line in Figure 3 (3.3 km/s contour) is in most places congruent with a rapid velocity increase at depth, which we interpret as the base of the Phanerozoic sedimentary cover at a depth that largely agrees with previous 1-D or 2-D local studies (e.g., Hatzfeld et al., 2003;Jahani et al., 2009) and the magnetic-basement model of Teknik and Ghods (2017). Note that our tomographic model cannot resolve the Hormuz salt layer, which is too thin (<5 km) to be recovered by our data set and may anyway have a velocity indistinguishable from either the upper basement or the bottom of the sediments. As inferred from the NW-SE profile A-A ′ in Figure 3, our velocity model reveals a variable thickness of Phanerozoic cover from ∼8 to 13 km depth, along a trend parallel to the strike of the orogen. Indeed, a prominent feature is that the greatest thickness of the low-velocity layer corresponds to the SE end of the MZTF (north of Musandam Peninsula). This is unsurprising, since the present-day GPS convergence rate is rather uniform across the ZFTB (Walpersdorf et al., 2006), but the width of the seismically active zone, and the fold belt itself, is roughly half the width further northwest. The shortening of sediments is therefore confined to a narrower width with a corresponding greater thickening. This concurs with the observations of Jahani et al. (2009), who used several seismic reflection profiles across the Zagros. Although they use two-way traveltime transects (processing of the seismic data is not mentioned, but we assume they show time-migrated profiles), with conversion to depth simply approximated by considering a fixed velocity value, they estimate a general increase in thickness of Cambrian-to-Quaternary cover ranging from 9 to 15 km, along strike from NW toward the SE and the Hormuz Strait. The pattern of teleseismically derived centroid depths appears to follow a similar trend (Figure 3, profile AA ′ ). The new velocity model presented in this work is consistent with the suggestion by Nissen et al. (2011Nissen et al. ( , 2014) that medium-to-large earthquakes are concentrated in the middle-lower part of the sedimentary section throughout the ZFTB (see histogram in Figure 4c). Similarly, nearly half of the revised focal depths from Karasözen et al. (2019) are consistent with nucleation within the sedimentary cover. Thus, for instance, earthquakes in Qeshm island (around 27 • N/55.5 • E) and in Iran north of the Musandam peninsula (see profiles AA ′ , CC ′ , and DD ′ from Figure 3) are distinctly deeper than those found in the central Fars arc (see, e.g., line BB ′ from Figure 3), following the pattern of increasing sediment thickness. In CC ′ and DD ′ , where the velocity model is well constrained, the profiles confirm the suture of the MZTF as the NE limit to the thick ZFTB sediments. This same feature is suggested in BB ′ , but here the model could be less well resolved in the upper crust (see Section 3 and Figure 2a).
At depths beneath the Zagros sediments, profile DD ′ (Figures 4a and 4b) shows that earthquakes continue NE beyond the MZTF and the limit of the low-velocity sedimentary cover, where low-angle thrust mechanisms indicate that the underlying Arabian Precambrian basement is transported beyond the surface exposure of the suture (Figure 4d), as suggested by Talebian and Jackson (2004), Gholamzadeh et al. (2009), andNissen et al. (2011). It is probable that here the Arabian basement undergoes greenschist facies metamorphism and may have wave speeds similar to those found further north in central Iran. Thus, it may not be surprising that we are unable to identify a boundary between the Arabian and Eurasian domains at depth. Yamini-Fard et al. (2007) used a temporary seismic network to study the transition between the southern Zagros collision and Makran subduction zone, claiming to image underthrusting of the Arabian lower crust beneath central Iran using a P-wave velocity model. However, their results come from a spot measurement centered on 27.8 • N/57.8 • E (covering an area less than 30 km by 30 km at 20 km depth) and their network is between 50 and 100 km outside the cluster of seismicity we discuss in section DD ′ . Our results suggest that this underthrusting is widespread on a broader region.
Profile BB ′ might also provide evidence for Arabian basement underthrusting beneath Iran, as the high-velocity region at around 10 km depth appears to be persistent for ∼80 km beyond the surface expression of the MZTF (see dash-dotted black line in Figure 3). This relatively subdued feature apparently dips at an angle that follows that of the Moho (as shown by Paul et al., 2006) as well as that of the low-angle thrusting earthquakes shown in Figures 4a and 4b. Nonetheless, this observation leads to the obvious question as to why seismicity beyond the MZTF is localized solely in the southeastern end of the belt (along profile DD ′ ). Is this the only site of the ZFTB where the basement of Arabia slid beneath and across the MZTF? Or has underthrusting occurred in other areas that are currently aseismic? As previously mentioned, profile BB ′ is located in an area where the Zagros orogen is wider than further east (i.e., along profile DD ′ ) and the seismicity is concentrated in the lower, SW part of the orogenic belt, as does the surface shortening inferred from GPS measurements (Walpersdorf et al., 2006). It is also apparent from GPS-based maximum shear strain rate maps that the highest deformation rates are presently localized in the southeastern Zagros, as recently shown by Khorrami et al. (2019). A possible conjecture is that underthrusting along the MZTF may have been a diachronous process, developing earlier in the central part of the belt. Subsequent internal heating would allow the Arabian Precambrian basement to become hot and weak, therefore rendering it unable to generate earthquakes. This is the same effect discussed by Craig et al. (2012) for southern Tibet, where temperature is the principal factor controlling how far seismicity extends to the north beneath Tibet within the underthrusting Indian shield. Eventually, the low-angle north dipping thrusting earthquakes of the eastern ZFTB will cease when strength of the Arabian basement will decrease due to heat production.

Conclusions
We have applied a fully probabilistic tomographic scheme that involves a tree-structured wavelet parameterization of the velocity field to a recent seismic network in the eastern Middle East. Our new 3-D shear velocity model is able to illuminate variations in the thickness of the sedimentary cover of the southern Zagros mountains to a greater level of detail than more conventional tomographic methods. When combined with accurate teleseismic centroid depths, the velocity pattern confirms that larger earthquakes in the Fars arc tend to occur in the lower sedimentary cover of the former Arabian passive margin. We also demonstrate that both the base of the sediments and maximum earthquake depths tend to increase along strike to the SE. Although we do not observe a clear velocity contrast between the basements of Arabia and Iran in the southernmost part of the ZFTB, integration of our 3-D tomographic model with known seismic focal mechanisms and depths provides confirmation that the lower crust (Proterozoic basement) of Arabia underthrusts beneath central Iran. In the central Fars arc, present-day shortening and seismicity within the sedimentary section are concentrated in the southern part of the orogenic belt. Further north, the underlying Arabian basement is sliding beneath Iran, although it appears to be unable to generate earthquakes.