New Constraints for the On‐Shore Makran Subduction Zone Crustal Structure

The Makran Subduction Zone is the primary seismic/tsunami hazard of the northwestern Indian Ocean, but little is known of its on‐shore seismic structure. We derived a shear wave velocity model extending to > $ > $ 100 km depth beneath a ∼400 km‐long seismic profile oriented parallel to the convergence vector of the Arabian Sea Plate. Receiver function/surface wave analysis shows that the average structure in the coastal region comprises a ∼22–28 km‐thick low wavespeed sedimentary cover and a 6–8 km‐thick gradient zone overlying > $ > $ 100 km‐thick high wavespeed upper mantle. The ocean‐basement interface dips gently northward, remaining a positive impedance contrast to ∼50 km depth at ∼250 km north of the coast where it disappears as the basaltic/gabbroic oceanic crust has probably transformed to eclogite. Further north, a weak arrival at ∼5 s in the receiver functions appears, grading northward into the Moho arrival of the continental Iranian Plateau. This disruption in the seismic signature of the Moho occurs in the forearc region where the dip of the subducting oceanic plate steepens. The southern Iranian Plateau's continental crust has an average Vs of 3.55 ± 0.05 km s−1, an almost flat Moho 40–45 km deep, and a sub‐Moho mantle Vs of 3.75 ± 0.05 km s−1 in the 50–80 km depth range. Weak Moho conversions probably result from ∼20% serpentinization of peridotite in the mantle wedge. Receiver functions indicate a flat continental Moho – no crustal root beneath the high topography region of the volcanic belt, which therefore must be compensated by low upper mantle densities. The high Vp/Vs ratio observed for the mantle wedge suggests ∼1%–2% partial melt.


Tectonic Setting and Results of Previous Work
The Arabian Sea Plate (Figure 1), the last remnant of the Neo-Tethys, is presently subducting roughly northward beneath the coast of southern Iran and Pakistan. Plate reconstructions imply that the convergence rate increases from 3.7 cm yr −1 in the western Makran to 4.2 cm yr −1 in the eastern Makran (DeMets et al., 1990), but GPS measurements estimate the present convergence rate of the western Makran to be only 2.7 cm yr −1 (Khorrami et al., 2019;Vernant, Nilforoushan, Hatzfeld, et al., 2004). In the west, the Makran is separated from the Arabian-Eurasian continental collision zone in the Zagros by the right-lateral Zendan-Minab-Palami Fault system (Figure 1) and in the east from the Indo-Eurasian continental collision zone by the left-lateral Ornach-Nal-Chaman Fault system (to the east of the area shown in Figure 1). The Zendan-Minab-Palami Fault system experiences a low level of seismicity (Yamini-Fard et al., 2007) and does not appear to be a major through-going lithospheric fault system. Instead, it is a localized transition zone accommodating the differential motion between the Zagros and Makran (Gholamzadeh et al., 2009). The western Makran experiences a lower level of seismicity than does the eastern Makran (Byrne et al., 1992) and a short-duration deployment of ocean bottom seismometers to the south of the western Makran recorded no earthquakes (Niazi et al., 1980). Subduction beneath the Makran probably started during Late Cretaceous time (Arthurton et al., 1982;Dykstra & Birnie, 1979). The age of the presently-subducting oceanic crust of the Arabian Sea Plate is not well determined due to a lack of seafloor magnetic anomalies, but it is probably of Jurassic or Cretaceous age (e.g., Hutchison et al., 1981;McKenzie & Sclater, 1971). The accretionary prism (Figure 1) is more than 350 km wide with only the frontal 100-150 km currently submarine. The topography of the accretionary wedge varies from ∼3,000 m below sea level to ∼1,500 m above sea level; however, there is no bathymetric trench. Most of the north-south convergence is absorbed in the subduction of the Arabian seafloor; however, a small amount of shortening has been observed between the Makran coast and the volcanic arc ∼250 km to the north (Khorrami et al., 2019;Vernant, Nilforoushan, Chery, et al., 2004). The shortening within the on-shore portion of the accretionary prism may contribute to the uplift of ∼0.2 mm yr −1 (Reyss et al., 1999). The prism contains extensive zones of Mesozoic age melange and large intact ophiolite ( Figure 1) (Ghazi et al., 2004). To the north of the Makran accretionary prism lies the subduction-related Jazmurian and Mashkel Depressions (Figure 1), and further north lies a sub-parallel andesitic-to-rhyolitic volcanic arc ( Figure 1) with the prominent centers of Bazman, Taftan and Koh-i-Sultan (Farhoudi & Karig, 1977).
The off-shore seismic structure of the Makran accretionary prism has been studied extensively (e.g., Fowler et al., 1985;Grando & McClay, 2007;Kopp et al., 2000;Kukowski et al., 2001;Niazi et al., 1980;Smith et al., 2012). However, there have been only a few seismological studies of the on-shore structure and they reach conflicting conclusions. Al-Lazki et al. (2014) found that a very low P n velocity ( 7.7 km s −1 ) beneath the western Makran and low P n wavespeed (∼7.9 km s −1 ) extends to the north of the Jazmurian Depression and the volcanic arc. Several surface wave studies covered portions of the Makran, but all have poor path coverage. The partitioned surface waveform tomography of Shad Manaman et al. (2011) depicts a shallow, high-wavespeed feature extending from beneath the Arabian seafloor to below the coastal region of the Makran and then dipping northward to about 250 km below central Iran. The ambient noise Rayleigh wave (16-40 s period) group velocity tomography of Abdetedal et al. (2015) found that the Mohorovicȋć discontinuity (Moho) beneath the Makran deepened from less than 25 km in the coastal region to more than 50 km beneath the volcanic arc. Abdollahi et al. (2018) inverted the Abdetedal et al. (2015) dispersion data jointly with free-air gravity data, finding that the crust varied from an average V s ∼3.95 km s −1 and ∼33 km thick near the coast to an average V s ∼3.75 km s −1 and ∼48 km thick beneath the volcanic arc. There are three published receiver function results for CHBR ( Figure 1). Taghizadeh-Farahmand et al. (2015) used HK-stacking and forward modeling and found the crust to be ∼33 ± 2 km thick with a V p /V s of ∼1.77. Penney et al. (2017) inverted the CHBR receiver function with a Rayleigh wave group velocity and found a low V s ∼3 km s −1 , 26 ± 2 km-thick layer, which they suggested were sediments of the accretionary wedge, overlying a ∼4.1 km s −1 oceanic basement. Penney et al. (2017) found a similar structure below the TURB (Figure 1) in the Makran of southwestern Pakistan. Motaghi et al. (2020) used HK-stacking and forward modeling of receiver functions from recordings at CHBR and nearby stations, concluding that the sediments were as much as 10.1029/2021JB022942 4 of 23 20 km thick beneath the coastal plain and overlay a 4.5-7 km-thick oceanic crust. Haberland et al. (2021) present a P-wave model for the on-shore structure of the accretionary prism from wide-angle seismic recordings. We discuss these published results in more detail in Section 5.

Data
The crust and uppermost mantle structure of the western Makran presented in this paper is derived from the analysis of teleseismic body-wave and local and regional surface wave recordings. These seismic recordings come primarily from a network of 23 temporary seismographs we operated between mid-2016 and late 2020 over the western portion of the Makran Subduction Zone in southern Iran (IASBS/CAM -red-outlined circles in Figure 1). In addition to using data from the temporary network, we include data from two seismographs operated by the International Institute of Earthquake Engineering and Seismology (INSN -blue circles, Figure 1) and one seismograph operated by the University of Tehran (IRSC -purple circles, Figure 1). Details of the profile stations are given in Table 1  from the coast more or less parallel to the plate convergence direction and with an average of ∼20-25 km station spacing. There is some deviation from a strict linear profile due to access and security constraints.
The fundamental mode group velocity dispersion observations for the western Makran are taken from the work of Ahmadzadeh- Irandoust et al. (2022). They used multiple-filter analysis (Dziewonski et al., 1969) to extract longer period (⪆15 s) dispersion information from vertical component seismograms of magnitude 4 and greater earthquakes recorded on Iranian seismographs, including the Makran stations. Their short period dispersion was determined from ambient noise analysis (Bensen et al., 2007) from the Makran network recordings. Ahmadzadeh- Irandoust et al. (2022) merged their dispersion measurements with those of Gilligan and Priestley (2018). They then employed the procedure of Debayle and Sambridge (2004) to invert path-average surface wave measurements on a 1° × 1° grid at sixteen periods between 5 and 70 s to obtain the geographic variation of the group velocity. We interpolated these dispersion maps to build 5-70 s pseudo-dispersion curves for each Makran network recording site.

Methodology
Teleseismic receiver function analysis (Burdick & Langston, 1977;Phinney, 1964;Langston, 1977) is now an established seismological tool for determining the V s and V p /V s structure beneath a seismograph site, and we apply this technique to determine the V s structure beneath the on-shore portion of profile AA' (Figure 1). The waveforms from events greater than m b 5.5 in the distance 30°-90° were band-pass filtered between 0.05 and 2.0 Hz, rotated from NS and EW to radial and transverse with respect to the event-station great-circle path, and receiver functions were computed using the iterative deconvolution method of Ligorria and Ammon (1999). To stabilize the time domain deconvolution, we applied a 2.5 Gaussian width filter (corner frequency ∼1.0 Hz) in the deconvolution. In addition, we computed receiver functions for each event using several different Gaussian filters (1.0, 1.6, 2.5, 3.5, 4, 5) to examine the variations in crustal response with frequency, but the Gaussian 2.5 receiver functions are the main focus of our analysis in this paper. To examine the backazimuth variation we form backazimuthal gathers for radial and transverse receiver functions after correcting for moveout (e.g., Figure 2a). The moveout-corrected receiver functions are then stacked to determine a mean site receiver function with error bounds (e.g., Figure 2b).
We then invert the mean site receiver function jointly with the pseudo-dispersion curve for each site using the linearized least-squares inversion algorithm of Herrmann (2013). Receiver functions are sensitive to high wave-number features beneath the recording location and do not provide strong constraints on the absolute wavespeed. Surface wave dispersion constrains the absolute wavespeed in a broad region surrounding the recording site but has a low sensitivity to discontinuities. The starting inversion model consisted of a finely parameterized structure defined by 2 km-thick flat-lying isotropic layers with S-velocity and V p /V s for each layer and extending to 280 km depth. The wavespeed of the model used in the inversion (e.g., Figure 2d) is based on the velocity model AK135 (Kennett et al., 1995) but with the crustal wavespeeds replaced by the uppermost mantle wavespeed of AK135 (4.48 km s −1 to 110 km depth). Since there is no a priori crustal model imposed, the inversion procedure will introduce layering as required by the data. We repeated the inversion by replacing the shallow structure with velocities of 4.28 and 4.68 km s −1 to examine the effect of the starting model on the outcome and for both of these, the resulting inversion model converged to the model found from the 4.48 km s −1 starting model. The receiver functions are given a higher weight (87.5%) compared to the surface waves (12.5%) in the inversion.
For the four representative recording sites whose results are shown in Figures 2-5), we first simplified the inversion model by combining similar adjacent 2-km-thick layers into thicker structures with either constant V s or constant V s -gradient to form a minimum parameter model which explained the data. We then tested the main features of the simplified model to evaluate their importance in fitting the data.
Thus, our underlying assumptions are that the structure beneath each seismograph site is close to one-dimensional (1-D), can be approximated by plane isotropic layers, and that the effect of side scattering from local crustal heterogeneity is minimal. Receiver function analysis normally assumes simple horizontal layering with the velocity changing only in the vertical dimension. Subduction zone environments can be significantly more complex due to dipping structures; however, a number of studies (e.g., Ferris et al., 2003;Owens et al., 1988;Yuan et al., 2000) have shown that the dipping structures can be successfully modeled in subduction zone environments using receiver function analysis. In this study, we do not attempt to reproduce all of the signal in each receiver function, 6 of 23 but we create a simple model for the on-shore portion of the Makran Subduction Zone which explains the stable features in the observed receiver functions and is consistent with other geophysical and geological observations.
The difference in spatial resolution of the receiver functions and surface wave dispersion suggests that we interpolate the receiver function wavefield to equalize the lateral sensitivity of the two data sets and simplify the complexity in the receiver functions. Following Chai et al. (2015), we interpolated the receiver functions recorded along the profile to acquire the interpolated receiver function wavefield. The moveout-corrected receiver functions were weighted and stacked with the receiver functions within ±20 km of the interpolated point weighted unity and with the weight linearly decreasing to zero between 20-80 km from the interpolation point. (2) simplified model; and (3) effect of removing the velocity increase beginning at ∼100 km depth. The long-dashed black line is the inversion starting model. The V s (z) structure to 150 km depth is shown in the inset of (d). We identify the Moho as a significant break in slope of V s (z) where V s first is very close to or exceeds 4.0 km s −1 . This is a clear transition for all sites except those marked with a "?" in Table 1. where, in addition to the mean stacked 2.5 Gaussian radial receiver function, stacked radial receiver functions for eastern (upper trace) and northwestern (lower trace) backazimuths are included. In (d), V s (z) models for the crust and uppermost mantle beneath AJNT: (1) V s model derived from inverting the mean receiver function and pseudo-dispersion data; (2) simplified model; (3) effect of removing the LVZ centered at ∼52 km depth; (4) effect of removing the LVZ centered at ∼68 km depth; (5) effect of removing the LVZ below ∼80 km; and (6) model derived from inversion of the northwestern backazimuth receiver function and pseudodispersion data. The dashed black line is the inversion starting model. 8 of 23

The On-Shore Makran Accretionary Prism
The on-shore coastal plain of the Makran is separated from the off-shore portion of the accretionary prism by a narrow coastal strip where normal faults and mud volcanoes are prominent (e.g., Back & Morley, 2016;Von Rad et al., 2000). In the coastal plain, the surface geology consists primarily of thick sequences of relatively homogeneous slope and shelf facies rocks which have experienced ∼20%-30% shortening with broad synclines up to 20 km across, separated by tighter thrusted anticlines. Further inland, in the central Coast Ranges, only turbidites are exposed. These have been tightly folded and imbricated, with shortening strains of 60% or more (Platt et al., 1985).
The recordings from the permanent station CHBR (Figure 1) are representative of the data from the temporary coastal sites, but the CHBR data are much more numerous. The receiver function vs. backazimuth observed at CHBR (Figure 2a) shows that over most backazimuths, there is a strong arrival at ∼4 s. In the 70-260° backazimuth range (all backazimuth ranges are clockwise from north), there is a significant phase at ∼2 s, and a broad arrival at ∼13 s. The first motion of the transverse receiver function changes from negative to positive at a backazimuth of ∼170° and back to negative at a backazimuth of ∼350°. The CHBR stacked, 2.5 Gaussian receiver function ( Figure 2b) contains a pronounced peak at ∼4 s, a minor peak at ∼2 s and a broad peak at ∼13 s. The higher frequency 5.0 Gaussian receiver function shows that the broad ∼13 s peak in the 2.5 Gaussian receiver function is composed of two closely-spaced but distinct peaks. The CHBR fundamental mode Rayleigh group velocity curve shows low wavespeeds at periods less than ∼20 s with an Airy phase at ∼18 s, a sharp rise in wavespeeds between 20-45 s period, and high wavespeeds ( 3.5 km s −1 ) at periods greater than ∼50 s.
The inversion result for the CHBR receiver function and dispersion ( Figure 2d) shows the average V s model for this site has a very thin (∼2 km, ∼2.1 km s −1 ) surface layer, a ∼10 km-thick shallow crustal layer with V s ∼2.9 km s −1 , below which there is a ∼10 km thick, ∼3.1 km s −1 layer. Beginning at ∼22 km depth, there exists a strong gradient (∼0.10 s −1 ) extending to ∼32 km depth where the wavespeed is ∼4.25 km s −1 . Beneath ∼38 km depth, there is a low velocity zone (LVZ) extending to ∼70 km depth. However, forward testing indicates that this  Figure 3. In (d), V s (z) models for the crust and uppermost mantle beneath TAFT: (1) model derived from inverting the mean receiver function and pseudo-dispersion data; (2) simplified model; (3) effect of removing the LVZ centered at ∼28 km depth; (4) effect of removing the LVZ centered at ∼76 km depth; (5) effect of removing the LVZ centered below ∼110 km depth; and (6) model derived from the inversion of the western receiver function stack and pseudo-dispersion data. The long-dashed black line is the inversion starting model. LVZ is not required by the receiver function and is only weakly constrained by the surface wave dispersion. Forward modeling tests also show that the double peak at ∼13 s consists of two peaks, one from a positive velocity increase starting at ∼100 km depth, and a second one which is a crustal reverberation. Thus, below ∼36 km, V s is nearly constant to ∼100 km depth below which there is a second strong positive V s gradient.

Jazmurian Depression
A significant change in the nature of the receiver functions and surface wave dispersion occurs near the latitude of the Jazmurian Depression, the sedimentary basin north of the Makran Range (e.g., Burg, 2018;Shahabpour, 2010). The basin, which may have only begun subsiding in the Late Pliocene (McCall & Kidd, 1982), has been the site of shallow marine sedimentation and volcanism, but the thickness of the basin has remained a question (Burg, 2018). The contact between the Jazmurian Depression and the Makran accretionary wedge has steeply-dipping normal faults related to the ongoing subsidence (Dolati & Burg, 2013).
Recordings from AJNT (Figures 1 and 3), which is located at the eastern edge of the Jazmurian Depression, are representative of the data from sites across this region. The receiver functions are complicated because they vary with azimuth, suggesting considerable structural complexity. The radial receiver functions have a positive arrival at ∼2 s over the 30-270° backazimuth range and a strong negative arrival at ∼5-6 s delay in the ∼260-50° backazimuth range. The first arrival on the tangential receiver function is primarily negative for eastern backazimuth and primarily positive for northwestern backazimuths but, because of gaps in the recordings, it is not possible to constrain the azimuth where the change in polarity takes place. There is a prominent coherent positive tangential arrival at 5-6 s spanning the ∼40-125° backazimuth range and a strong negative conversion at ∼5 s in the 215-270° backazimuth range. The mean AJNT receiver function ( Figure 3b) has a positive arrival at ∼2 s, broad bounds between ∼3.5-6.5 s that reflect the variability of the radial receiver functions at these delay times, and a negative arrival at ∼8 s. Figure 3b shows stacked receiver functions computed from events with east and northwest backazimuths. Because a majority of the events recorded at AJNT are from eastern backazimuths, the mean and eastern receiver function stacks are similar. The northwestern receiver function stack, which consists of a smaller number (11) of individual receiver functions, differs from the mean receiver function with a secondary arrival interfering with the direct arrival, a positive arrival at ∼3 s, and a strong, negative arrival at ∼6 s. The AJNT short period ( 20 s period) group velocity is somewhat higher (∼2.7 km s −1 ) than that observed in this period range at CHBR (∼2.50 km s −1 ), but the long period dispersion ( 40 s) is slower (∼3.25 km s −1 ) than that which is observed in this period range at CHBR (∼3.45-3.50 km s −1 ).
AJNT has a complicated inversion V s (z) model with multiple LVZs. Teleseismic arrivals from the east partially sample the crust of the Sistan Suture Zone whereas arrivals from the northwest sample the crust of the Jazmurian Depression. Inversion of the mean AJNT receiver function results in a crustal model consisting of a low V s (∼2.6 km s −1 ) surface layer, a ∼22 km-thick layer with a 0.015 s −1 gradient where V s increases from ∼3.0-3.3 km s −1 at ∼28 km depth, and a ∼14 km thick layer with a 0.05 s −1 gradient where V s increases from ∼3.3-4.0 km s −1 at 42 km depth. Below this, there is a weak LVZ which forward modeling shows is not a robust feature of the model. Beneath this, V s increases to ∼4.1 km s −1 before again decreasing to below 4.0 km s −1 . Forward modeling shows that this decrease is also not a robust feature of the V s (z) model. V s then increases to ∼4.5 km s −1 at ∼82 km depth but then decreases to less than ∼3.8 km s −1 at ∼105 km depth. This thick, deep LVZ is not constrained by the receiver function but is required by the dispersion data. The V s (z) model derived from the northwestern stack is similar to the model derived from the mean receiver function to a depth of ∼30 km, but between 40 and 50 km depth there is a significant LVZ with V s dropping to ∼3.2 km s −1 . This feature is the source of the strong negative arrival at ∼6 s. Between ∼50 and 65 km depth there is a 0.08 s −1 gradient with a V s 4.5 km s −1 at 65 km depth. Below this V s drops to ∼3.7 km s −1 and remains low to at least 150 km depth.

Sistan Suture Zone, Iranian Plateau
North of the Jazmurian Depression, the profile crosses the Sistan Suture Zone and the volcanic arc. The Sistan Suture is an abandoned, deformed accretionary prism (Bröcker et al., 2013;Tirrul et al., 1983) extending northward from the central Makran for ∼700 km along the Iran-Afghanistan border. This feature formed when a branch of the Neo-Tethys Ocean closed between ∼65 and 90 Ma, uniting the continental Lut and Helmand Blocks.
The data from ZHSF (Figures 1 and 4) are representative of the recordings from stations located on the Sistan Suture Zone away from the volcanic arc. There is a positive arrival at ∼5 s spanning the entire backazimuth range, a significant negative arrival ∼7.5-8 s spanning most of the ∼25-130° backazimuth range and a second significant negative arrival at ∼335-125° at ∼2 s. The tangential receiver function is smaller in amplitude compared to the tangential receiver function observed at sites to the south. The tangential component first arrival is negative in the ∼75-215° backazimuth range and positive in the ∼220-50° backazimuth range. The stacked ZHSF receiver function (Figure 4b) is relatively simple with a prominent positive arrival at ∼5 s, a significant negative arrival at ∼2 s, a minor negative arrival at ∼8 s, and minor positive arrivals at ∼9.5 and 17 s. The fundamental mode Rayleigh group velocity curve is remarkably flat with short-period ( 25 s) group wavespeeds of 2.9-3.0 km s −1 and long-period ( 30 s) group wavespeeds of ∼3.25 km s −1 .
The ZHSF inversion model is simple compared to the V s models for the more southern sites. The shallow portion of the ZHSF model consists of a thin, ∼2.85 km s −1 surface layer, a ∼8 km thick ∼3.25 km s −1 layer, a ∼28 km thick ∼3.4 km s −1 layer, a strong gradient (0.14 s −1 ) between 38 and 42 km, a weak gradient (∼0.01 s −1 ) between 42 and 64 km depth, and a deep V s ( 3.65 km s −1 ) LVZ. The sensitivity test shows that the upper mantle LVZ is weakly constrained by the ZHSF receiver functions but is required to fit the ZHSF surface wave dispersion.

Makran Volcanic Arc
TAFT is located on the slopes of the ∼4,000 m high andesite strato-volcano Taftan (Figure 1) which is built on the ∼2 km-high Sistan Suture Zone of the Iranian Plateau. Taftan consists of two main summits, Narkuh and Matherkuh. Narkuh is highly eroded while Matherkuh is partly covered by fresh-appearing, thick andesite lava flows. The Taftan edifice is primarily built of lava flows, but ignimbrites and pyroclastic flows also occur (Biabangard & Moradian, 2008). Zircon U-Pb data yield a cluster of ages at ∼3.2 Ma but with some samples giving dates as young as 0.8 Ma (Pang et al., 2014), whereas 40 Ar/ 39 Ar dates range from ∼7 to 0.7 Ma (Biabangard & Moradian, 2008). There are no credible reports of historical eruptions, although there is some speculation of minor activity in 1902 (Global Volcanism Program, 2020). The geochemical signatures of the Taftan magmas suggest that they were derived from the lithospheric-enriched melts by subduction metasomatic processes (Biabangard & Moradian, 2008). Saadat and Stern (2011) propose that, based on Sr and Nd isotropic ratios, the magmas that erupted at the Makran volcanic arc have no significant crustal contamination, implying a relatively thin ( 40 kmthick) crust. They postulate that enrichment of large ion-lithophile elements, relative to light rare-earth-elements and depletion in Nb relative to large ion-lithophile elements are similar to most other convergent plate boundary arc basalts, suggesting that Makran basalts are formed by the melting of subcontinental mantle modified by dehydration of subducted Arabian Sea Plate. Furthermore, Saadat and Stern (2011) find that the Pb isotropic ratios of the basalts are consistent with a contamination of subducted sediments.
There is a coherent positive arrival at ∼2 s in the radial receiver function in the ∼10-120° backazimuth range followed by a coherent negative arrival at ∼3.5 s and a significant negative arrival at ∼2 s between ∼170-190° which possibly extends to ∼340° backazimuth. The first arrival of the tangential receiver function is positive in the ∼100-190° range and negative across the ∼300-75° backazimuth range. There is a strong, negative arrival at ∼2 s for southeastern backazimuths. The mean and eastern backazimuth receiver functions have a prominent positive arrival at ∼2 s and a minor positive arrivals at ∼5 and 7 s. The western backazimuth receiver function has a strong negative arrival at ∼2 s and a minor positive arrival at ∼5 s. The TAFT group velocity curve is very flat with the mean group wavespeed never exceeding ∼3 km s −1 .
The V s structure derived from inverting the TAFT mean receiver function and the surface wave group velocity (Figures 5b-5d) is unusual in that at no point shallower than 150 km depth does V s exceed ∼4.1 km s −1 . The inversion model consists of a ∼4 km-thick surface layer with V s less than 3 km s −1 , a ∼20 km-thick, ∼3.2 km s −1 layer, a small LVZ centered at ∼28 km depth, below which there is a positive gradient (0.050 s −1 ) to ∼58 km depth where V s reaches ∼3.85 km s −1 . Below this, there is a second negative gradient with V s dropping to ∼3.25 km s −1 at ∼78 km depth. This structure overlies a third LVZ with V s dropping to ∼3.25 km s −1 at ∼140 km depth. Sensitivity tests show that the shallower of the two deep LVZs is only weakly constrained by the surface wave dispersion and the deep LVZ is not constrained by the data. Inversion of the western backazimuth receiver function results in a structure above ∼55 km depth, which is similar to that of the inversion result of the mean receiver function but reaching a slightly higher V s of ∼4.0 km s −1 at ∼58 km depth. Below this, V s remains approximately constant to ∼72 km depth but beneath this, there is an extensive LVZ extending to at least 130 km depth.

Interpolated Receiver Function Profile
In the previous subsections, we discussed inversion results for seismic data recorded at sites located on four representative tectonic structures across the Makran Subduction Zone of southern Iran. The inversion results for the remaining recording sites are summarized in Table 1. To equalize the lateral sampling of the receiver functions and surface wave dispersion, we have interpolated the receiver functions as described in Section 3.2, and the interpolated wavefield and its inversion model are shown in Figure 6. As pointed out by Chai et al. (2015), interpolation of the receiver function wavefield accentuates the first-order features in the receiver function response and suppresses the large-amplitude fluctuations at individual stations that complicate receiver function analysis. In addition, the interpolation of the receiver function wavefield compensates for the difference in spatial sampling of receiver functions and surface waves when the two data are jointly inverted.  (Kopp et al., 2000); the three red lines denote the top of the sedimentary layer, the top of the oceanic basement and oceanic Moho found in the controlled-source experiment located somewhat to the southwest of the profile (Niazi et al., 1980). The short-dashed red line beneath the southern portion of the receiver function profile indicates our interpretation of the oceanic Moho with the dotted line where the Moho is less clear. The long-dashed red line starting at ∼400 km indicates a possible Moho beneath the volcanic arc and Sistan Suture Zone. Over-plotted on (c) are the earthquake mechanisms ±200 km from the profile taken from Penney et al. (2017). The focal spheres have been rotated by 90° and displayed from the side to correspond to the cross-section view of the wavespeed structure. A wavespeed scale for the V s structure is shown at the left in (c). (d) The variation of GOCE free-air gravity along the profile. All distances refer to the southern end of the plot. Figure 6b shows a smooth variation in a strong, positive arrival at ∼5 s delay near the coast, to ∼4 s delay immediately south of the Makran Range. This arrival disappears at ∼7 s delay at ∼375 km distance along the profile beneath the eastern Jazmurian Depression. There is a second positive arrival at ∼2 s delay which has significant amplitude at the coast but then decreases in amplitude and disappears ∼100 km north of the coast. North of ∼410 km, there is a low amplitude positive arrival at a delay time of ∼5 s which increases in amplitude north of ∼425 km. Beneath the region of the volcanic arc, there is a strong positive phase at ∼2 s immediately followed by a strong negative arrival at ∼3.5 s with another significant negative phase at 7-8 s delay which has its peak amplitude beneath the volcanic arc. Figure 6c displays the variation in the V s structure beneath the profile derived from the inversions of the interpolated receiver functions and surface wave dispersion. The shallow structure is composed of material with < 3.5 km s −1 across the whole profile. In the southern portion of the profile there is a thin band of material with a wavespeed of 3.6-3.9 km s −1 at a depth ∼25-30 km which broadens to the north of ∼350 km. North of ∼425 km the very slow wavespeed, near-surface material thins and approaches zero thickness at the northern end of the profile. However, very slow wavespeed (∼3.6 km s −1 ) occurs from ∼40 km depth and V s 4.1 km s −1 wavespeed material persists to at least 100 km depth beneath the volcanic arc.

Discussion
We have derived a shear wave velocity model that extends to a ∼100 km depth beneath a ∼400 km-long profile of seismographs oriented approximately parallel to the convergence vector and stretching from the Makran coast to a point ∼400 km north on the southern Iranian Plateau (Figure 1). The model was obtained by jointly inverting receiver functions from seismic data recorded at 26 sites along the profile with fundamental mode Rayleigh group velocity dispersion determined from ambient noise and regional earthquake analysis. Figure 7 summarizes the Makran structure we discussed in Section 4. This model confirms the conjectures of the large-scale features proposed in early papers on the region (e.g., Farhoudi & Karig, 1977;Jackson & MKenzie, 1984). The incoming Arabian Sea Plate begins subducting at a very shallow angle beneath the accretionary complex. Sediments as thick as 25 km accumulate over the incoming plate (brown scalloped), and the trench/forearc is clogged with sediment. Stresses due to the thickening sediment and tectonic compression cause shortening and uplift (e.g., Berberian & King, 1981;Haghipour & Burg, 2014;Jackson & MKenzie, 1984;Penney et al., 2017;Page et al., 1979;Snead, 2002) of the sediments that form the coastal plane and the Makran Range.
As the oceanic crust is buried deeper beneath the sediment, the shear wavespeed of the basaltic/gabbroic crust (black to brown transition) likely alters to blueschist metamorphism. When the oceanic crust reaches a depth of ∼50 km, the basaltic/gabbroic crust begins to transform to eclogite (green) and the increased density causes this portion of the plate to become gravitationally unstable and the angle of subduction becomes steeper. This takes place in the vicinity of the Jazmurian Depression. The entrained sediments and the metamorphic reactions transforming the basaltic/gabbroic oceanic crust to eclogite release large volumes of H 2 O-rich fluids (red arrows) into the overlying mantle wedge (pink). The rising fluids lead to the serpentinization of the mantle wedge peridotite, lowering the density and velocity. An effect of the metamorphic reaction taking place in the oceanic crust and the overlying mantle wedge is to render the Moho invisible in the receiver functions from the nose of the mantle wedge (Moho gap). Farther north, the seismic profile crosses the continental crust of the Sistan Suture Zone and the volcanic arc. The crust is ∼40 km thick beneath the Sistan Suture Zone, similar to the crustal thickness beneath central Iran (Ahmadzadeh- Irandoust et al., 2022) and there is no discernible continental root beneath the elevated region of the volcanic arc.

Coastal Sedimentary Structure
The high sedimentation rate and shallow angle of the incoming plate has led to the formation of an accretionary prism ( Figure 1) ∼350 km wide, of which only the frontal 100-150 km is currently submarine. During the initial period of the Arabian Sea Plate subduction, the accretionary wedge grew by accreting Himalayan turbidites deposited in the paleo-Indus deep-sea fan (Grando & McClay, 2007;Kukowski et al., 2001). The uplift of the Murray Ridge System during the Early Miocene (Gaedicke et al., 2002) shifted the Indus sediments to the south and, since that time, the Makran accretionary complex has been maintained largely by the recycling of sediments coming from the erosion of the older prism and conveyed to the front of the modern accretionary complex through a series of river and canyon systems (Haghipour & Burg, 2014;Haghipour et al., 2015), and by incorporating abyssal sediments scraped off the subducting Arabian Sea Plate (Burg, 2018;Grando & McClay, 2007;Platt et al., 1985). The present-day Himalaya-derived sediments contribute little to the Makran accretionary wedge . Mud volcanoes are prevalent along the Makran coastal region (Negaresh, 2008;Wiedicke et al., 2001), suggesting over-pressured sediments at deeper depths in the accretionary wedge. Tectonic compression and stresses due to the thickening sediment cause shortening and deformation of the sediments as seen in off-shore folding and faulting (e.g., Grando & McClay, 2007;Kopp et al., 2000), the uplifted coastal marine terraces (e.g., Page et al., 1979;Snead, 2002), the on-shore folding (e.g., Farhoudi & Karig, 1977), thrust faulting (e.g., Berberian & King, 1981), uplift (e.g., Haghipour & Burg, 2014) of the Makran Range, and the high angle (40-50°) reverse-faulting mechanisms for shallow earthquakes (e.g., Jackson & MKenzie, 1984;Penney et al., 2017).
Modeling of receiver function and surface wave data from coastal sites shows the presence of a ∼22-28 km-thick low wavespeed sedimentary cover. Near the coast, the sedimentary section V s increases from ∼2.5 km s −1 at the surface to ∼3.7 km s −1 at ∼22 km depth. The on-shore sedimentary structure of our model is consistent with the off-shore structures determined from controlled-source seismic recordings (e.g., Fruehn et al., 1997;Grando & McClay, 2007;Kopp et al., 2000;Niazi et al., 1980;;Smith et al., 2012). Similar thick sedimentary structures are observed in other present-day accretionary structures and in accretionary structures in the geologic record (Cawood et al., 2009).
The sparse wide-angle recordings of Niazi et al. (1980) (Figures 1 and 6c) ∼200 km southwest of our profile, found the oceanic crust to be ∼20 km-thick and consisting of a 13-14 km-thick low wavespeed (V p ∼4 km s −1 ) layer over a 6-7 km thick oceanic crust. Near-vertical controlled-source seismic recordings (Grando & Mc-Clay, 2007) along several profiles ∼100 km southwest of our seismic profile found that at the deformation front, abyssal sediments are being scraped from the underthrusting Arabian Sea Plate, forming an imbricate fan of active thrust faults. Further north, they found that a previously deposited 4 km-thick Himalayan turbidite sequence is underplating and uplifting the wedge and the mud volcanoes are spatially related to the onset of the underplating of sediments (Grando & McClay, 2007). Kopp et al. (2000) found, from active-source data collected off the Pakistan coast ∼200 km east of our profile, that the upper part of the accretionary prism consists of ∼10 km of unconsolidated sediments overlying a thick layer of more consolidated underplated sediment (Figure 6c). They determined that the shallow sediment layer is characterized by a low velocity gradient (∼0.16 s −1 ). From their wide-angle measurement, Kopp et al. (2000) conclude that the décollement steps down toward the north, thus allowing for substantial underplating, and they suggest that ∼3 km of the incoming sedimentary section is thrust beneath the deformation front. Mass balanced calculations suggest that ∼50% of the down-going sediments are incorporated into the wedge by subcretion (Platt et al., 1985(Platt et al., , 1988. The projection of the off-shore structure derived from the controlled-source seismic experiments (e.g., Ellouz-Zimmermann, Lallemant, et al., 2007;Fowler et al., 1985;Fruehn et al., 1997;Kopp et al., 2000;Smith et al., 2012) is consistent with the on-shore structure we find in our study (Figure 6c). Haberland et al. (2021) determined a compressional velocity model for the on-shore accretionary wedge structure from active seismic source recordings along three ∼200 km-long lines perpendicular to the coast. Their easternmost profile is coincident with the receiver function profile in this study. Figure 8 compares the average V p from the active seismic experiment with the average V s from the receiver function/surface wave analysis at two points, one at the coast and a second ∼50 km north of the coast (Figure 1). The V p averages at the two points are formed from the tomographic model for the easternmost seismic profile of Haberland et al. (2021). The average V p structure near the coast (point a, Figure 1) consists of a strong, ∼2-3 km thick near-surface gradient zone where V p increases from ∼3.2-5 km s −1 , overlying a ∼10 km-thick, V p ∼5 km s −1 upper crustal layer, below which there is positive gradient where V p increases to ∼7 km s −1 at ∼22 km depth (Figure 8a). Haberland et al. (2021) do not see the basement near the coast. Further inland (point b, Figure 1), the average V p structure consists of a 12 km thick gradational section with V p increasing from ∼4.5 km s −1 at the surface to ∼5.9 km s −1 . Below this, there is a ∼10 km thick section with a near uniform V p of ∼5.9 km s −1 above a positive gradient (0.10 s −1 ) extending to the base of their model (Figure 8b). Between 60-90 km inland from the coast, there is a band of reflecting elements at 30-35 km depth dipping northward at ∼8° ± 2° which Haberland et al. (2021) associate with the base of the subducting oceanic crust. Also included in Figure 8 are V p /V s profiles for the two points. In the upper ∼15 km near the coast the V p /V s -ratio is close to 1.73 but deviates below ∼15 km depth with V p /V s reaching ∼2.1 at ∼23 km depth. This high V p /V s layer could indicate over-pressure in the subducting sediments and the source of the mud volcanoes near the coast. At the point ∼50 km from the coast the V p /V s is close to 1.73 at all depths.

The Subducting Arabian Sea Plate
Beneath the sedimentary cover near the coast at depths greater than ∼23 km, there is a 6-8 km-thick strong gradient layer in which V s increases from ∼3.5-4.2 km s −1 . This feature is compatible with the receiver function signature expected from layers 2 and 3 of the oceanic crust. The base of the gradient has a shallow 4 ± 2° northward dip (Figure 6c). This is consistent in depth and dip with the band of reflectors observed in the controlled-source seismic data (Haberland et al., 2021). The receiver function arrival at 5-6 s is a composite of energy from several conversions, but is primarily from the oceanic Moho (Figure 2b). This phase disappears ∼250 km north of the coast at the northern fringe of the Jazmurian Depression at a delay time corresponding to a depth of ∼50 km (Figure 6b). Between the coast and the Jazmurian Depression where seismology shows the flat slab, the free-air gravity anomaly is 25 ± 5 mgals (Figure 6d).
The only information on the structure of the basaltic/gabbroic oceanic crust off the western Makran is the low-resolution model from the short-duration OBS recording of Niazi et al. (1980) who found a 6-7 km-thick, V p ∼6.7 km s −1 layer below the sediments which they hypothesized was the oceanic crust. This V p wavespeed and oceanic crustal thickness is consistent with the V s wavespeed and thickness observed near the coast in our model. In a more detailed active-source seismic study conducted off the eastern Makran coast, Kopp et al. (2000) found an oceanic crust ∼9 km thick dipping about 3° northward and consisting of a shallow low velocity gradient and a strong velocity contrast (V p 5.7-7.0 km s −1 ) between oceanic layers 2 and 3. Figure 6b shows a decrease in the Ps delay time of approximately a second at ∼220 km distance beneath the near-coastal sites, suggesting a significant shallowing of the Moho. This is not a prominent feature in the inverted velocity model (Figure 6c). Haberland et al. (2021) found a similar uplift in the Moho structure from their controlled-source recordings along a refraction profile that coincides approximately with our receiver function profile. Unpublished work (Irandoust, 2021;Mokhtarzadeh, 2022) indicates that this feature extends east-west along the coast almost perpendicular to the profile. The position of this apparent uplift is approximately where one might expect the outer rise uplift (Melosh & Raefsky, 1980), but the amplitude is much too great for this to be an outer rise structure similar to that observed in other subduction zones. The amplitude of this uplift decays to the west, but it has a similar amplitude eastward into Pakistan (Irandoust, 2021). The receiver function profile location coincides approximately with the boundary between the western and eastern Makran. Subduction in eastern Makran is hindered by the stable block to its north, whereas subduction in the west has less opposition from the overriding plate. It may possibly be that this change in subduction has caused the uplift of the Moho.  2015) recognize a converting interface 33 ± 2 km below CHBR which they identified as the Moho and concluded from this that there was no crustal thickening beneath CHBR nor evidence of collisional processes taking place in this region. Penney et al. (2017) surmised that the strong gradient at 26 ± 2 km depth in their CHBR receiver function inversion model was the top of the oceanic basement; they found a similar crustal structure to the east beneath the Makran of Pakistan (TURB, Figure 1). Motaghi et al. (2020) used HK-stacking of receiver functions from recordings from CHBR and the nearby station NGCH and arrived at a similar conclusion to that of Penney et al. (2017). Motaghi et al concluded that the double peak at ∼13 s in the receiver functions consisted of the arrival of reverberations from the top and bottom of a 4.5 km-and 7 km-thick hidden-layer oceanic crust beneath NGCH and CHBR, respectively. However, their conclusion implies that the oceanic crust would have an anomalously low V p /V s . Our sensitivity test (Figure 2) shows that this double peak at ∼13 s in the CHBR receiver functions does arise from interfering arrivals, one being the crustal multiple and the second a conversion from a V s increase at ∼100 km depth.
The nature of the Moho conversion below the receiver function profile is similar to the Moho pattern observed beneath Cascadia (e.g., Audet et al., 2009;Bostock et al., 2002;Nicholson et al., 2005). In the south across the Makran accretionary prism, there is a strong Moho conversion. Below the Jazmurian Depression the Moho conversion in the receiver functions becomes weak and disappears ( Figure 6b) at a delay time corresponding to a depth of ∼50 km at a distance about ∼250 km north of the coast. This implies a weak or non-existent velocity contrast which probably signifies the transformation of the basaltic/gabbroic oceanic crust to eclogite. V s for basalt and eclogite at this depth are ∼4.0 km s −1 and ∼4.7 km s −1 respectively, compared to peridotite ∼4.6 km s −1 (Christensen, 1996;Christensen & Mooney, 1995), rendering the oceanic Moho nearly invisible. The accompanying increase in density, possibly by as much as 15% (e.g., Brocher, 2005) would cause the incoming Arabian Sea Plate to become gravitationally unstable and descend more steeply into the mantle. However, it cannot be ruled out that the disappearance of the oceanic Moho arrival is due to a sudden steepening in the dip of the subducting Arabian Sea Plate.
Our profile recording allows us to track the subducting oceanic lithosphere from beneath the accretionary prism to ∼50 km depth. Earthquake mechanisms with well-constrained hypocentral depths (e.g., Penney et al., 2017) ( Figure 6c) allow the subducting oceanic plate to be traced to somewhat deeper depths. Near the coast, there are several shallow thrust-faulting events. Farther north, there is a group of normal-faulting earthquakes at depths ∼50 km (at 400-500 km on Figure 6c) in the vicinity of where the Moho arrival disappears. The M w 7.7 Khash earthquake (Barnhart et al., 2014) resulted from normal faulting at ∼75-80 km deep somewhat to the east of the receiver function profile. Penney et al. (2017) pointed out that the normal-faulting earthquakes have an ∼ENE-WSW linear distribution in map view parallel to the strikes of their nodal-planes and are observed across the longitudinal extent of the subduction zone. Penney et al. (2015) took the fault planes of the shallow thrust earthquakes ( Figure 6c) to be the northward-dipping nodal planes and then, if these events occur on the subduction interface, the dips of 8-10 ± 5° of their nodal planes provide an estimate of the dip of that interface near the coastline. This is close to, but somewhat steeper than, the dip of the subducting interface inferred from the receiver-function/surface-wave and active-source models. Penney et al. (2017) further suggested that the dip of the seismogenic zone inferred from the earthquake distribution appears to steepen to the north of the normal-faulting earthquakes and proposed that the normal-faulting events probably represent extension at a hinge in the subducting Arabian Sea Plate.
One of the normal-faulting earthquakes occurred almost directly beneath the Bazman volcano at a depth of 74 km (Jacob & Quittmeyer, 1979). This implies that the subduction interface is no deeper than ∼90 km under this part of the volcanic arc, allowing for a ±15 km possible error in the earthquake depth. The location of the volcanic arc is affected by various factors besides the depth of the slab, but globally, the depth of the subducting slab has been found to be in the 100-120 km depth range below the volcanic arc (Syracuse & Abers, 2006). Global body-wave tomographic studies (e.g., Van der Meer et al., 2018) show a steeply (∼45°) northward-dipping high V p -wavespeed structure in the upper mantle (∼280-600 km depth) beneath central Iran. Models of the shallow upper mantle structure (e.g., Shad Manaman et al., 2011) show the presence of a northward-dipping high wavespeed structure. When projected back to a shallow depth, both features intersect approximately with the region where the receiver function signature of the oceanic plate disappears.
One feature of the oceanic slab structure is perplexing -the south-to-north V s variation. High V s wavespeeds in our model occur near the coast (100-125 km) and further inland (200-325 km) (Figure 6c), but there are zones with lower wavespeed (125-200 km and 325 km) (Figure 6c). This portion of the structure is primarily constrained by the surface wave dispersion and there can be a reduction in the amplitude of the deeper structure due to the decreased sensitivity of the surface waves with depth. However, the lateral variation in V s along the profile suggests that this is not the cause since the sensitivity of the surface waves is very close to uniform along the profile. The observed variation may arise from faulting in the crust and uppermost mantle of the subducting oceanic lithosphere. Bending-related faulting in subduction zones has been recognized world-wide (e.g., Contreras-Reyes et al., 2011Grevemeyer et al., 2007;Ivandic et al., 2008;Ranero & Sallarès, 2004) and these faults can provide pathways that allow water to infiltrate the crust and uppermost mantle, altering the basalt/gabbro of the oceanic crust and peridotite of the uppermost mantle and affecting the seismic wavespeeds. V p and V s in the uppermost mantle can decrease from 8.2 to 7.7 km s −1 and 4.7 to 4.2, respectively (e.g., Christensen, 1972Christensen, , 2004Grevemeyer et al., 2018). If this is the cause of the lower wavespeed of the subducting Arabian Sea Plate, it is puzzling as to why this would affect portions of the Plate differently.

Sistan Suture Zone, Volcanic Arc and Mantle Wedge
North of the Jazmurian Depression, the receiver function/surface wave inversion shows that the continental crust of the southern Iranian Plateau has an average V s of 3.55 ± 0.05 km s −1 , an almost flat Moho at a depth of 40-45 km based on the nearly uniform Ps delay time (Fig. 6b), and a sub-Moho mantle V s of 3.75 ± 0.05 km s −1 in the 50-80 km depth range. However, the location of the Moho beneath the volcanic arc is not clear in our model ( Figure 5d). In our inversion model (Figure 6c), the Moho appears as a weak velocity increase. In a global subduction zone survey, Bostock (2013) suggested that the large volume change associated with the transformation of the basaltic/gabbroic oceanic crust to eclogite may rupture the plate boundary, allowing fluids to penetrate the forearc mantle wedge leading to the serpentinization and melting of the mantle wedge rocks. The serpentinization and the free fluids can diminish, erase or even invert the velocity contrast of the Moho of the overriding plate. Below the Sistan Suture Zone, there is a thin upper mantle lid below which is a substantial mantle LVZ where V s drops to ∼3.5 km s −1 at ∼130 km depth. Sensitivity tests (Figures 4 and 5) show that the very low V s in the uppermost mantle is required.
Our seismic profile cuts across the volcanic belt between the active volcanoes of Bazman and Taftan; one of our recording sites is on Taftan (Figure 1). There is a pronounced low velocity layer (LVL) in the crust seen in the receiver function as the significant negative arrival that starts at ∼410 km along the profile and terminates at ∼580 km ( Figure 6b). This appears as a LVL in the lower crust in the inversion structures of the individual receiver functions ( Figure 5).
In subduction zones, the crust of the incoming plate comes in contact with the crust of the overriding plate in the forearc region and then descends into the mantle. Wilson (1954) was the first to suggest that H 2 O-rich fluids are generated from the Wadati-Benioff zones. In addition to the release of H 2 O from the entrained sediments as the oceanic crust is carried to deeper depths in the mantle, the metamorphic reactions transforming the basalt/ gabbroic oceanic crust to eclogite release large amounts of H 2 O-rich fluids into the overlying mantle wedge (e.g., Christensen, 2004;Peacock, 1993Peacock, , 1996. The rising fluids cause extensive hydrothermal alteration and lower the solidus of the material in the overlying mantle wedge, causing serpentinization and melting of the mantle peridotite. Serpentinization of the mantle wedge has now been documented in most subduction zones where detailed seismic studies have taken place (e.g., Eberhart-Phillips et al., 2005;Grevemeyer et al., 2018;Husen et al., 2003;Kamiya & Kobayashi, 2000;Kuwatani et al., 2011;Reynard, 2013;Seno et al., 2001;Sodoudi et al., 2006;Xia et al., 2008). Serpentinization causes a substantial lowering of the wavespeed and density of the mantle wedge (Peacock, 1996). One effect of the serpentinized peridotite is to reduce the seismic signature of the overlying continental Moho, and, in extreme instances, can result in an inverted Moho (Bostock et al., 2002). The conversion from the Arabian Sea Plate Moho is clear beneath the accretionary prism but deepens and disappears at ∼400 km (Figure 6b). In this same region, earlier weak arrivals exist which are probably the first indication in the receiver functions of the conversion from the continental Moho. This disruption in the seismic signature of the Moho occurs in the Makran Subduction Zone forearc region where the dip of the subducting Arabian Sea Plate steepens. Bostock et al. (2002) proposes that the weak conversions from the continental Moho may result from serpentinization of the peridotite in the mantle wedge of the Cascadia Subduction Zone.
Between the continental Moho and the inferred depth of the Arabian Sea Plate, the V p (Al-Lazki et al., 2014) and V s (this study) are low and give a V p /V s of ∼2.05, significantly higher than possible for 100% antigorite (Ji et al., 2013). The receiver functions indicate a nearly flat continental Moho across the region (Figure 6b), and there is no crustal root beneath the elevated region of the volcanic belt. The free-air gravity anomaly over the volcanic belt is ∼75 mGal (Figure 6d) indicating partial compensation. However, the nearly flat Moho implies the partial compensation of the volcanic zone topography occurs from a low density sub-Moho mantle. The average elevation of the volcanic arc is ∼1.75 km. Assuming a mean density of 2850 kg m −3 for the crust, the free-air gravity anomaly would be ∼200 mGal. Therefore, partial compensation of the elevation must occur in the mantle. If this compensation occurs in a 60 km-thick layer between the base of the crust at ∼40 km depth and the inferred depth of the subducting Arabian Sea Plate at ∼100 beneath the volcanic arc, it implies a reduction in density of 50 kg m −3 to account for the observed gravity anomaly or an upper mantle density of 3,150 kg m −3 , similar to that inferred for other subduction zones (Hyndman & Peacock, 2003).
Although the mantle wedge V s beneath the volcanic arc is consistent with 100% serpentinization of the peridotite to antigorite (Ji et al., 2013), this degree of serpentinization is incompatible with the observed free-air gravity anomaly over the volcanic belt and the inferred V p /V s of the mantle wedge from the V p (Al-Lazki et al., 2014) and V s found in this study. The high V p /V s ratio compared to that for 100% antigorite suggests that the low V s we observe for the mantle wedge results from a fluid-rich or partially molten mantle wedge. Hammond and Humphreys (2000) find ∼8% reduction for upper mantle V s for realistic partial melt geometries per percent of partial melt. The V s we observe in the mantle wedge beneath the volcanic arc therefore implies the presence of ∼1%-2% partial melt. This is consistent with the findings of Saadat and Stern (2011) that Makran arc basalts formed as a result of a low degree of partial melting of the mantle wedge which is only moderately modified by subducted components.

Summary
We have derived a shear wave velocity model extending to ∼100 km depth beneath a ∼400 km-long profile of seismographs oriented approximately parallel to the convergence vector and stretching from the Makran coast to a point ∼400 km north on the southern Iranian Plateau (Figure 1). The model was obtained by jointly inverting receiver functions and fundamental mode Rayleigh group velocity dispersion from seismic data recorded at 26 sites along the profile. Modeling of receiver function and surface wave data from coastal sites shows the presence of a ∼22-28 km-thick low wavespeed sedimentary cover. Near the coast, the sedimentary section V s increases from ∼2.5 km s −1 at the surface to ∼3.7 km s −1 at ∼22 km depth. A comparison of the V s structure of the accretionary prism sediment cover with the V p results of Haberland et al. (2021) shows that in the upper ∼15 km near the coast the V p /V s -ratio is close to 1.73 but deviates below ∼15 km depth with V p /V s reaching ∼2.1 at ∼23 km depth. This high V p /V s layer could indicate over-pressure in the subducting sediments and the source of the mud volcanoes near the coast. At the point ∼50 km from the coast the, V p /V s is close to 1.73 at all depths.
Beneath the sedimentary cover near the coast there is a 6-8 km-thick strong gradient layer in which V s increases from ∼3.5-4.2 km s −1 , a feature that is compatible with the receiver function signature expected from layers 2 and 3 of the oceanic crust. The base of the gradient has a shallow 4 ± 2° northward dip (Figure 6c), consistent in depth and dip with the band of reflectors observed in the controlled-source seismic data (Haberland et al., 2021). The nature of the oceanic Moho conversion below the Makran profile is similar to the oceanic Moho pattern observed beneath a number of other subduction zones (e.g., Fukao et al., 1983;Maekawa et al., 1993;Ohmi & Hurukawa, 1996;Peacock & Wang, 1999;Rondenay et al., 2008). In the south across the Makran accretionary prism, there is a strong Moho conversion. Below the Jazmurian Depression, the Moho conversion in the receiver functions becomes weak and disappears ( Figure 6b) at a delay time corresponding to a depth of ∼50 km. This implies a weak or non-existent velocity contrast which probably signifies the transformation of the basaltic/gabbroic oceanic crust to eclogite.
North of the Jazmurian Depression, the receiver function/surface wave inversion demonstrates that the continental crust of the southern Iranian Plateau has an average V s of 3.55 ± 0.05 km s −1 , an almost flat Moho at a depth of 40-45 km, and a sub-Moho mantle V s of 3.75 ± 0.05 km s −1 in the 50-80 km depth range. The conversion from the Arabian Sea Plate Moho is clear beneath the accretionary prism but deepens and disappears at ∼400 km ( Figure 6b). In this same region, the weak arrival at ∼5 s in the receiver functions appears and grades to the north into the continental Moho arrival. This disruption in the seismic signature of the Moho occurs in the Makran Subduction Zone forearc region where the dip of the subducting Arabian Sea Plate steepens.
The fact that weak continental Moho conversions may result from serpentinization of the peridotite of the mantle wedge has been proposed for Cascadia (Audet et al., 2009;Bostock et al., 2002;Nicholson et al., 2005). The V p ≈ 7.75 km s −1 observed in the Makran mantle wedge (Al-Lazki et al., 2014) suggests that the degree of serpentinization is ∼20%, in the range inferred for other subduction zones (e.g., Chou et al., 2009;Kamiya & Kobayashi, 2000;Ranero & Sallarès, 2004;Van Avendonk et al., 2011;Wada et al., 2008). Serpentinization of 15%-30% would reduce the density to ∼3,000 kg m −3 (Hyndman & Peacock, 2003). The receiver functions indicate a nearly flat continental Moho across the region ( Figure 6b); that is, there is no crustal root and the high topography region of the volcanic belt must be partially compensated by low densities (∼3,150 kg m −3 ) in the mantle, similar to that inferred for other subduction zones (Hyndman & Peacock, 2003). The high V p /V s ratio compared to that for 100% antigorite suggests that the low V s we observe for the mantle wedge likely results from ∼1%-2% partial melt of the mantle wedge, consistent with the petrochemistry of the Makran arc basalts (Saadat & Stern, 2011).