Crustal melting and continent uplift by mafic underplating at convergent boundaries

Domain: earth_science
Published: 2024-10-19
Editorial outcome: minor_revision

Abstract
The thick crust of the southern Tibetan and central Andean plateaus includes high-conductivity, low-velocity zones ascribed to partial melt. The melt origin and effect on plateau uplift remains speculative, in particular if plateau uplift happens before continental collision. The East Anatolian Plateau (EAP) has experienced similar, more recent uplift but its structure is largely unknown. Here we present an 80 km deep geophysical model across EAP, constrained by seismic receiver functions integrated with interpretation of gravity data and seismic tomographic, magnetotelluric, geothermal, and geochemical models. The results indicate a 20 km thick lower crustal layer and a 10 km thick mid-crustal layer, which both contain pockets of partial melt. We explain plateau uplift by isostatic equilibration following magmatism associated with roll-back and break-off of the Neo-Tethys slab. Our results suggest that crustal thickening by felsic melt and mafic underplate are important for plateau uplift in the EAP, Andes and Tibet.

Introduction
The presence of partial melt within the crust of the southern Tibetan and central Andean plateaus is inferred from low seismic velocities and low electrical resistivity 1–9. Despite similarity in geophysical observations, the origin and dynamics of the crustal partial melt are attributed to different processes in the two plateaus. Models include partial melting due to crustal thickening in the Tibetan Plateau 2 and by intrusion of asthenospheric material in the central Andes 10. The timing and cause of plateau uplift and crustal thickening remain controversial, and non-uniqueness of interpretations of crustal structure has contributed to competing geodynamic models. The rise of the Tibetan Plateau has been attributed to underthrusting of Indian lower crust 11, oblique subduction of lithospheric blocks within the plateau 12, and lateral crustal flow 13,14. In the central Andes, the removal of dense mantle lithosphere and the addition of magma to the lower and middle crust are proposed to promote plateau uplift 10,15.

The uplift of the Gangdese Arc in the Lhasa Block of the Tibetan Plateau began >60 Ma 16,17, predating continent-continent collision 18, and the central Andean Plateau has been uplifting for 23 Myr 19. The long uplift history of these plateaus introduces challenges in inferring their deep structure at the onset of uplift. Therefore, processes leading to the initiation of plateau uplift may better be studied at young collisional plateaus. The uplift of the East Anatolian Plateau (EAP) in Turkey with a topography of ~1.5–2.0 km began around 12 Ma in relation to the collision between the Arabian Plate and Eurasia 20 with a current convergence rate of 20–30 mm/yr 21.

The EAP has very low Vs anomaly (around −6 to −9% with respect to regional average) in the uppermost mantle 22,23, low Pn velocity (generally Vp < 8.0 km/s and some zones have Vp < 7.8 km/s) 24,25, and high lithosphere temperature (800–900 °C at Moho) 26. It is broadly covered by Neogene to Quaternary volcanic rocks 27,28. Dynamic support of the plateau by upwelling asthenosphere has been proposed 29 based on the assumption of a 38–46 km thick crust, as indicated by low-frequency (f ≈ 0.72 Hz) Receiver Function (RF) studies, which interpret a strong P-wave converter as the Moho under the assumption of constant average crustal Vp/Vs = 1.73 throughout the region 30. This assumption of constant Vp/Vs ratio contradicts seismic tomography models which include significant low-velocity bodies, especially in the active volcanic provinces 22,23, suggesting the presence of intra-crustal, partially molten zones. These zones would require high crustal Vp/Vs and should also produce strong P-RF converters at their base, which easily would be mistaken for the seismic Moho.

Here we apply a refined method 31 for imaging the crustal structure of the EAP by seismic RF. We combine data from a temporary deployment of seismic stations at ~50 km interval 32 and permanent stations from the Kandilli Observatory across the EAP (Fig. 1 and Supplementary Table 1). Our 450 km long RF profile extends from N to S across the EAP in its central part with abundant volcanic fields, the Bitlis Suture, and the northern part of the Arabian Plate. We image strong and weak P-RF converters, and we apply a new gravity-based interpolation method to identify the Moho (Methods) 31. We demonstrate that the base of the zone that includes partially molten bodies was earlier interpreted as the Moho 30. To reduce the bias in our interpretation of the RF results, we constrain the upper 80 km of the lithosphere by integrating the RF results with regional seismic tomography 22,23,33,34, magnetotelluric 35, geothermal 26 and geochemical models 27.

Results
Crustal receiver function profile
Our P-wave RF model (Fig. 2a) shows the following:
1. The Arabian Plate has a crustal thickness of 35–40 km, as identified by a strong converter in an otherwise non-converting environment.
2. Shallow positive converters close to the margin of the Arabian Plate mark a main detachment which we interpret as the Main Bitlis Thrust (MBT) at a depth of 10–15 km. This blind thrust 36 shares strong similarity to the Main Himalaya Thrust in Nepal 11,37,38.
3. The Bitlis Massif to the south of the Muş Suture has a crustal thickness of ~40 km, with a gradually N-dipping Moho converter (M2) to ~50 km depth in the EAP.
4. The EAP includes two positive converters: M1 at ~30 km and M2 at ~50 km depth (Fig. 2a and Supplementary Table 1). Typically, M1 has larger amplitude than M2 (Supplementary Fig. 1). Additionally, we observe a strong negative converter TL which is located ~10 km shallower than M1 within the EAP. These sub-horizontal converters are observed on all stations within the EAP and are not identified outside the EAP.

The crustal layer between TL and M1 has a very high Vp/Vs ratio of ~2.0 and a low Vs velocity ranging from 2.9 to 3.1 km/s based on calculations of synthetic RFs (Fig. 2b and Supplementary Fig. 2). Results from inverting the RFs for velocities also suggest a low Vs of 2.96-3.33 km/s except for station CMCY in northern EAP (Fig. 3a, Supplementary Figs. 3 and 4). The lower layer between M1 and M2 has a high P-wave velocity (Vp ~ 7.2 km/s) and a high Vp/Vs ratio (~2.0) in its lower part (Fig. 2b).

Other geophysical observations
We use the long-wavelength filtered Bouguer anomaly for resolving the density contrast across the Moho and for regional scale isostatic calculations, whereas the resolution of the gravity data (~9 km laterally 39) may be too low to model the dimensions of individual melt pockets. A body wave travel time tomographic velocity model 33 shows that the crustal layer between M1 and M2 (30–50 km depth) in the EAP has a Vp velocity between 6.6 and 7.2 km/s (Fig. 3f). Gravity modeling (Supplementary Fig. 5b) suggests that the layer between M1 and M2 has a high density of 3.05 g/cm3, assuming a density of 3.20 g/cm3 in the hot uppermost mantle at 50-70 km depth (Fig. 2c). This layer also includes several bodies with low electrical resistivity of 10–30 Ωm 35 (Fig. 3g).

The low-velocity anomaly between TL and M1 (20–30 km depth) cannot be clearly resolved by body wave tomography due to vertical resolution limitations 33. It is, however, observed in Rayleigh wave data at a period of 30 s 22. Strong low-density (2.60 g/cm3, Fig. 2 and Supplementary Fig. 5b) and low-resistivity (~1 Ωm) zones 35 (Fig. 3g) are unique for the EAP in the study region. The converters TL, M1, and M2 along with magnetotelluric anomalies are only observed in the volcanic province of the EAP, which has 0.5–1.0 km higher topography than neighboring blocks (Fig. 3).

Crustal converters and the presence of magma
The parameters of the layer between TL and M1 (~1 Ωm electric resistivity 35, Vs ~ 3.0 km/s 22, and Vp/Vs ratio ~ 2.0 at 20–30 km depth) are similar to a layer containing crustal silica melt observed in the central Andes, where joint inversion of surface wave dispersion and RF constrains Vs = 1.9-3.1 km/s 40, with average Vs ~ 3.1 km/s from RF synthetic modelling 3, and electric resistivity lower than 1 Ωm 7,41. The Vp/Vs value of ~2.0 between TL and M1 indicates a distinct Vs decrease in this layer which may be caused by the presence of either fluids or partial melt. Considering the young magmatic activity in the EAP (Fig. 4a) and the temperature close to the felsic solidus at 30 km depth (Fig. 4b), the presence of pockets with felsic partial melt is the most likely explanation for this depth interval at 20–30 km. Heat flow data are absent along our profile, but the 12–18 km deep Curie temperature depth indicates high heat flow 26. For a dihedral angle of 20–40° in the equilibrium geometry model, with an aspect ratio of 0.10–0.15 42, we estimate the melt volume fraction by a forward model for partial melting 43 to be between 3 and 13% for a Vs value between 2.96 and 3.33 km/s (Supplementary Fig. S6a).

The MT model (Fig. 3g) indicates the presence of three major, individual melt bodies in the EAP separated by ~70 km, but the resolution of the available data is insufficient to provide measures of their dimension due to leakage of the low-resistivity in the model. The depth to the top of these bodies (~20 km) corresponds to the top (TL) of the upper low-velocity zone (LVZ), which extends continuously across the EAP (Fig. 3c), and the LVZ thickness may indicate the vertical extent of the melt pockets. The receiver functions and the seismic tomographic models have insufficient resolution to outline the individual partially molten bodies, but the recurrent low-velocity zone at all stations in the EAP indicates that the material between TL and M1 generally is close to the solidus 44, also outside the bodies containing partial melt.

Previous RF interpretations selected M1 as the Moho 30,45 because it is the strongest converter in the RF section. By the new method 31 based on long-wavelength gravity-interpolation of converters from the Arabian Plate we, instead, find that the strong converter represents the strong contrast between the low-velocity partially molten layer and the high-velocity (Vp ~ 7.2 km/s) and high-density (~3.05 g/cm3) lower crust, most likely with mafic composition (Supplementary Fig. 5b). The lower part of this layer has a Vp/Vs of ~2.0 and an electric resistivity of 10–30 Ωm 35 which indicates the presence of a, not yet fully solidified, mafic magmatic material. High Vp of ~7.2 km/s is typical for mafic rocks at 1.4 GPa and 900 °C (~50 km depth in the EAP) 43. This Vp value is lower than in the underplated lower crust with Vp = 7.4-7.6 km/s at the Baikal Rift Zone 46, but it is similar to other underplated layers like in the Kuril Arc, the Norwegian-Danish Basin and the Donbas Basin 47. It is unlikely that this layer may represent the underthrusted Neo-Tethys oceanic slab, since tomographic images and geochemical results indicate that the oceanic slab has broken off 27,34.

Geodynamic evolution and magmatism
Most of the recent magmatism 28 along the profile is concentrated in the EAP with variable magma composition ranging from felsic to mafic and a clear pattern of north to south younging of the onset of magmatic activity since 12 Ma (Fig. 4a and Supplementary Fig. 7). The occurrence of very young magmatism, with age less than 1 Myr, suggests the presence of magma chambers and partially molten bodies within the crust that can still be detected using geophysical data. The predicted geotherms 26 intersect the experimental wet mafic rock solidus 48 at depths of ~30–50 km and cross the wet felsic solidus 49 at a depth of ~30 km (Fig. 4b). This indicates that the mafic lower crust between M1 and M2 at 30 to 50 km depth, especially its lower part with high temperature, may still be partially molten which is consistent with the Vp/Vs of ~2.0. Further, the low velocity zone between TL and M1 at 20–30 km depth probably includes partial melt of felsic composition, considering that the estimated temperature of ~600 °C at a depth of 30 km 26 is below the wet mafic melting solidus, the high Vp/Vs of ~2.0 and the low resistivity of ~1 Ωm.

We explain the lateral change of magmatism from 12 Ma to 3 Ma by spatial changes in subduction dynamics and geometry. Assuming that the subducting Neo-Tethys slab induces melting of the mantle wedge at a depth of ~70 km (Fig. 4c), and observing the magmatism of various ages at distances of 345 ± 15 km (~12 Ma), 295 ± 36 km (~5 Ma), 230 ± 21 km (~3 Ma), the subduction angle would have steepened from 20 ± 1.5° (~12 Ma) to 26 ± 6° (~5 Ma). Alternatively, the magmatism might have been caused by delamination of the lithospheric mantle and lower crust 28,50, but this mechanism does not explain the clear southward younging of magmatism initiation. We propose that underplating 47 (addition of mafic magma to the lower crust or below the original Moho) led to magma fractionation and that the accompanying release of fluid and heat initiated melting of preexisting continental crust, which may explain the presence of the mid-crustal partially molten bodies.

Slab roll-back and break-off in the EAP have earlier been proposed based on geochronology, geochemistry and seismic tomography 27,34. We extend this interpretation by our geophysical model (Fig. 4d), which suggests that the magmatism moved laterally from north to south due to slab roll-back and steepening followed by slab break-off, and that this evolution caused intensive addition of mafic magma to the lower crust and intra-crustal partial melting.

Methods
Receiver functions (RF)
The (radial) Ps RF is calculated by time domain iterative deconvolution by ITERDECON 58 with a Gaussian filter of 2.5, which corresponds to a cut off frequency of ~1.24 Hz. RFs were processed for epicentral distances of 30–90° and each RF from each station has been manually inspected (Supplementary Fig. 9). RFs with a large difference from the overall characteristics of other RFs were rejected in the final RF results. We manually picked the P-wave arrival time of each RF (Supplementary Fig. 9c). For some seismic phases with weak amplitude, we also used the multiple-taper correlation (MTC) 59 to estimate the RFs (Supplementary Fig. 9d); a comparison of the RF results from different deconvolution methods ensures recognition that the weak seismic phases are not noise artifacts. Synthetic RF tests were calculated by the IRFFM (Interactive Receiver Function Forward Modeller) software 60.

Assuming a single-layer medium, the crustal thickness H and the crustal average κ = Vp/Vs are determined by H-κ stacking 61 of the converted Ps wave and the reflected multiples from the Moho interface. We use weight parameters of 0.6, 0.3, 0.1 for the Ps, PpPs, PsPs + PpSs phases to stack the peak amplitudes (Supplementary Fig. 1b). The sensitivity of the results to the assumed constant Vp = 6.3 km/s is insignificant within uncertainty (Supplementary Fig. 10). We record all possible H-κ pairs, and use the gravity field data to determine the final H-κ values according to a new method 31 (see details in the section on Gravity interpolation of RF results).

RFs were migrated into the depth domain (Figs. 2a, 3b) based on the assumed crustal P-wave velocity of 6.3 km/s, and the Vp/Vs ratio (Supplementary Table 1) retrieved from the H-κ stacking with application of gravity interpolation of the RFs (Supplementary Fig. 11). The final RF and velocity profiles are obtained by applying a moving average filter with a radius of 30 × 2 km (horizontal × vertical) to the migrated RF rays and the inverted Vp and Vs profiles.

Velocity inversion of RFs
We estimate vertical P- and S-wave velocity profiles from the RFs by applying three different methods: linear RF inversion 62, non-linear RF inversion by the transdimensional Markov-chain Monte Carlo (MCMC) algorithm 63 and the Neighborhood Algorithm (NA) 64,65. The results of the non-linear algorithms are similar, including the identification of the low-velocity zones in the middle crust. The first algorithm assumes a constant Vp/Vs = 1.73. We first check that this method leads to similar results as an earlier application for the EAP 30, when constrained to the upper 40 km of the crust. Application of the same method to more realistic depths of 50 and 56 km for detecting the Moho, produces results that include the pronounced low-velocity zone between 20 and 30 km depth as well as a Moho at around 50–60 km depth (Supplementary Fig. 3a). By extending the covered depth interval, the visual fit between calculated and observed RF is substantially improved (Supplementary Fig. 3a).

The non-linear transdimensional MCMC inversion algorithm 63 estimates both Vp and Vs profiles (and thereby Vp/Vs). Given the strong variation in Vp/Vs along the EAP profile, this method provides better resolution of depth and velocity variation, and it leads to a better fit between calculated and observed RF (Supplementary Fig. 3b). The results identify a clear Moho transition in addition to a strong velocity contrast at the base of the mid-crustal layer with partially molten bodies. The determined vertical velocity profiles for individual stations correlate very well between each other along the profile as documented by the velocity profiles (Fig. 3c, e). The resulting vertical Vs profiles also correlate with the occurrences of high-conductivity bodies in the MT model (Fig. 3g). The NA inversion 64,65 identify the same low-velocity zones beneath stations within EAP although the results are less detailed as the NA is restricted to 6 layers in the models (Supplementary Fig. 12). Each inversion included 1000 iterations with 20 models generated at each iteration and 10 Voronoi cells re-sampled at each iteration. The results are consistent with results from the MCMC inversion which confirms that the inverted Vs profiles are robust.

As the inversion was carried out independently for each station, the strong lateral continuity of the resulting vertical velocity profiles testifies to the presence of anomalous bodies (melt pockets) in the same depth interval along the profile, and it justifies the horizontal interpolation of the results between the stations.

Gravity calculations
The horizontal resolution of ~9 km and the amplitude accuracy of ~9 mGal of the EIGEN-6C4 gravity data 39 maybe insufficient for modelling the individual molten bodies in the crust with unknown size. Instead, we use the long wavelength gravity signal for calculation of the isostatic topography and for resolving the density contrast across the Moho. The 2-D gravity forward calculations are made with the FastGrav software. The Bouguer gravity anomaly is calculated for the EIGEN-6C4 gravity data 39 and the ETOPO1 topography data 66, assuming a crustal density of 2.67 g/cm3. To model the gravity signal from a deep Moho, the short-wavelength Bouguer anomaly fluctuations, originating from the shallow crust, are removed with a Gaussian filter of 500 km 31.

Model 1 (Supplementary Fig. 5b) explains the observed gravity anomaly and is consistent with the RF profile. It includes the crustal layer with partially molten bodies (with bulk layer density of 2.60 g/cm3, Vs ~ 3.0 km/s, Vp ~ 6.0 km/s) and the mafic underplated layer (density of 3.05 g/cm3, Vs of 3.6-4.1 km/s, Vp ~ 7.2 km/s) for an asthenosphere density D1 of 3.2 g/cm3. Model 2 without the mid-crustal low-velocity layer and the underplated lower crust (Supplementary Fig. 5c) cannot fit the gravity data unless the density of the asthenosphere beneath the EAP 22 is unrealistically high (>3.3 g/cm3, cf. green line) or the density of the lower crust is high (≥2.95 g/cm3, cf. black line), corresponding to an underplated lower crust, which supports the presence of an underplated lower crustal layer.

Gravity interpolation of RF results
We apply a new method 31 for regional determination of Moho depth variation by minimizing the difference between the long period gravity signal and potential Moho depths along the profile, selected from the results of H-κ stacking 31. The method selects the optimum Moho depth at each seismic station from several possible candidates derived from the H-κ stacking by testing the crustal models against gravity. We adopt an empirical correlation between regional long-wavelength (>248 km) Bouguer gravity anomaly (ΔgB) and Moho depth (D), so that a lower value of the Bouguer gravity anomaly corresponds to a deeper Moho: ΔgB = a + bD, where a and b are empirical regional parameters. For a known ΔgB, we minimize the difference between the Moho depth obtained as (ΔgB − a)/b and multiple Moho candidates from the regional RFs by scanning the a, b parameter space. This procedure 31 determines the approximate Moho depth for any given ΔgB value, and allows us to select the closest RF Moho depth (Supplementary Fig. 11) 31.

Predicted topography
To calculate the predicted topography in the three plateaus, we assume regional isostasy with pressure equilibrium at a depth of 80 km below the different tectonic units. The Arabian Plate is close to isostatic equilibrium with a free air (FA) gravity anomaly of ~20 mGal, and we adopt it as reference model. The lithostatic pressure is Pref = t0 ρ0 + Σhi ρi, where t0 and ρ0 are thickness and density of the material above sea level, hi and ρi are thickness and density of different layers from sea level to 80 km depth. Isostatic topography for other tectonic units is therefore tiso = (Pref − Σhi ρi)/ρ0. The predicted free air gravity anomaly is FA = 2πG(t0 ρ0 + Σhi ρi). At full isostatic equilibrium (and neglecting lithosphere flexure) FA should be near-zero, but the central Andes and Tibet have FA anomalies of ~90 mGal and ~40 mGal, respectively. We assume that the extra ΔFA anomaly originates from a non-isostatic topography component of Δt = ΔFA/2πG/ρ0 in addition to the isostatically predicted topography, tiso, with the Arabian Plate adopted as reference. This non-isostatic correction allows for including the observed FA into the isostatic calculation.

Quantitative comparison with tomography and MT models
Age data of magmatic rocks in Eastern Anatolia 28 within a ±100 km wide corridor is projected onto the profile (Fig. 4; Supplementary Fig. 7). The absolute P-wave and S-wave velocity profiles 33 are located along the RF profile. The MT model 35 is smoothed by parameters of alpha = 1, tau = 335, to compare to the RF result, as these parameters provide the smallest misfit between the magnetotelluric data and the resistivity model 35.

We calculate the P-wave and S-wave velocities at different melt fraction, pressure and temperature by forward modeling 43 of seismic velocity and electrical conductivity and compare with the absolute P- and S-wave velocity profiles 33 and the smoothed resistivity model 35 along the RF profile (Fig. 3; Supplementary Fig. 6).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The present work aims to construct an 80 km deep geophysical model across the East Anatolian Plateau (EAP). To achieve this target, seismic receiver function was done and its results were integrated with interpretation of gravity data, and seismic tomographic, magnetotelluric, geothermal and geochemical models. I believe that this study will make a significant contribution to understanding the internal dynamics of the region. After these minor revisions, I believe it can be published in this journal.\nMinor comments.\n1. The presence of an intra-crustal low-velocity zone in Eastern Anatolia, or perhaps more significantly, an \"intracrustal asthenosphere,\" appears inevitable in my view. While the overall suggestion holds merit, the explanation for the existence of the low-velocity zone in the text is notably weak. Specifically, the assertion that the \"MBT main Bitlis thrust\" is associated with subduction lacks sufficient supporting data. The assumption that the segments observed in tomographic studies belong to the Arabian plate is entirely speculative. Recent papers, such as \"Portner et al 2018 P-wave tomography\" and \"Kounoudis et al 2020, seismic imaging of the Eastern Mediterranean mantle,\" present models with segments that cannot be reasonably inferred. As a result, the \"underplating\" assumption does not seem to have a well-supported foundation. The subduction and the subduction location of the Arabian plate should be discussed in detail. If necessary, alternative option(s) for the \"underplating\" assumption should be proposed and argued. The subduction zone from the front of the Bitlis massif should be considered in light of the results of this study.\n2. The content of abstract looks very weak, I think it would be better if it will rewrite again.\n3. Please prepare a list consisting the number of each station’s measurement.", "reviewer_2": "Based on the new results of Receiver functon(RF), in combination with the results of regional seismic tomography, magnetotelluric, geothermal and geochemical models, this paper reveals the architecture of the upper 60 km of the East Anatolian Plateau lithosphere and emphasizes the importance of mafic underplating in generating the high attitude of the EAP, Andes and Tibetan plateaus. The seismic images are of high quality and the interpretations are mostly reasonable and reliable. This paper is well organized in general, and the results would have broad impacts for geologic society. I would recommend its possible publication in Nature Communications after moderate revisions.\nMy comments are listed below in an order of line numbers:\nLines 13-14 in abstract: Pre-collisional break-off of the Neo-Tethyan slab ------ break-off generally occurs after collision; granitoid melting ---- it is strange, please considered to revise it, e.g., generation of felsic rocks …\nLine 4 in the Introduction section: This depends on the timing for generating the crustal partial melt.\nLines 5-7: These models have not been mentioned in later discussions.\nLine16: Reference 16 should be moved forward and expressed as 14-16 and new critical reference (e.g., Hu et al., 2016, ESR; Zhu et al., 2015, SR) should be added to the original site of 16.\nLines 105-106: This half sentence seems contradicts with later description. I would suggest delete it.\nLine 117: This pattern is not clear and too simplified. The original paper by Schleiffarth et al. (2018) interpreted it as lithospheric delamination as traditionally interpreted. Can you please add a new figure to show the spatial and temporal variations of rocks since 12 Ma?\nLines 129-144: This part is very weak as (1) lithospheric delamination can also result in intensive addition of mafic magma to the lower crust and intra-crustal partial melting and (2) you did not consider the effects of trench retreat and tectonic shortening since 12 Ma.\nLine 167: The seismic “bright spot” in southern Tibet is located in middle-upper crust, not middle-lower crust, given that the Tibetan crust is very think since 20 Ma (up to 70-80 km).\nLines 177-180: Please note that the high-density, high velocity layer in southern Tibet is not mafic, but eclogitized mafic and ultramafic rocks as summarized by Zhu et al. (2022, NC; 2023, AREPS).\nLines 180-183: I am not sure the logic here as granitoids are related to partial melting of the (ultra)mafic lower crust or fractional crystallization of parental magma, and its thickness is also related to exhumation. When did the granitoid molten form? Is it coeval with the onset of plateau uplift? Otherwise, this comparison is meaningless.\nLine 190: Felsic melting is strange. This sentence would be changed to: but the thickness of continental crust may not be insufficient to result in extensive crustal melting generating felsic rocks.\nLine 195: This sentence is unnecessary.\nLines 199-201: Essentially, pre-collision plateau uplift is related to isostasy of thickened crust due to mafic underplating. This has been addressed by Zhu et al. (2017).\nLines 201-203: the presence of continental crust above a subduction zone ----- do you mean a continental margin/arc? If yes, do you think the EAP is under an continental arc setting, not intraplate?\nI hope these comments are helpful and looking forward to seeing the publication of this important paper.", "reviewer_3": "The manuscript proposes that mafic lowermost crust and partial melt at shallower levels characterize convergent boundaries, such as in Tibet, in the Andes, and in the Eastern Anatolian Plateau, the focus of this study. The geoscientific analysis is based on seismic receiver functions, to highlight the main discontinuities in the lithosphere, and complemented with the analysis of gravity anomalies, magnetotelluric results, thermal field, and some geochemical data. The conclusions draw strong similarity between the three areas, based on which the authors generalize the prior conditions and formation processes in question for plateau uplift.\nThe main topic is very relevant and interesting to a broader community, and the multi-disciplinary approach of the authors is laudable. There are, however, also clear weaknesses, which do prevent claiming that the results are robust. These are (1) the scale and resolution of the investigations, (2) the lack of quantitative uncertainty assessment and the lack of quantitative link between most (not all) of the different observables, (3) some missing but key references. I develop these points in more detail below, and add further remarks.\n(1) Typical scales and resolution of the results.\nThe authors claim that the lower crustal mafic layer as well as partial melt prevails over several hundreds of kilometers, in a fairly homogeneous way. This statement cannot be supported by the data presented in this study, because the seismic stations (a total of 16, spread heterogeneously over 450 km distance, see Fig. 1) leave very large gaps along the profile (see ray coverage in Extended Data Figure 8d). Moreover, the gravity data, being derived mainly from satellite measurements, is not able to resolve any feature smaller than ca. 100 km, and it has been filtered with a Gaussian of 500 km (!). In the meantime, the MT results of prior studies (Extended Data Figure 4d) clearly shows local anomalies, and NOT a generalized layer of low resistivity. This is very similar to what is observed in more recent MT studies across the Tibetan plateau, and what seismological analysis of Hi-CLIMB stations show: local, isolated low-velocity zones, with ca. 10 km thickness and a horizontal extent no longer than 50 km. For these reasons, the results are not robust and generalization of local features seen at low spatial resolution and coverage is not convincing.\nTo highlight this problem with an example: in Extended Data Figure3b bottom panel, there are 8 receiver functions in the zone of interest, within the 200 km long zone with an upper crust, partial melt, D2 and D1 layers. We should see positive RF amplitudes at the bottom of the partial melt layer (bottom of low-velocity zone), and of the D2 layer (Moho). On two of the eight RFs the top converter is not visible (not visible on the 3rd from the left, and offset in depth on the 5th from the left). At the Moho level, the 4th from the left has no positive amplitude (there are two adjacent ones).\n(2) Uncertainties, simplifying assumptions.\nAlthough the study features a very nice diversity of datasets, there is limited direct link between them, and the simplifying assumptions of the RF analysis is a clear limitation. The H-kappa stacking method is very good for the first order, and I do agree with the distinction between deeper and shallower converter, however the depth estimate can be seriously off due to the need to assume a single average Vp value for the entire crust with this approach (here: Vp=6.3km/s). The uncertainty analysis does not take into account this factor. Other approaches to invert for the Vs structure with depth (e.g. with the Neighbourhood Algorithm method) would provide much better constraints, and would allow to come up with more credible quantitative depth and Vs estimates.\nThe connection between RFs and gravity using ref.27 seems to be good idea, however there is no other quantitative link between any other pair of datasets.\nA related point: it is not addressed how the distribution of various density layers would behave on the long term, on the timescale of convergent boundary evolution. If the partial melt layer is so thick, at a density of 2.6 g/cm3 and a melt fraction of ca. 7-14% (see Ext. Data Fig. 5a) which makes it mobile, how does it stay below a layer with 2.7 g/cm3 density?\nHeat flow data and any thermal model are not shown in the study, although that would likely be a very useful argument.\n(3) References.\nThe manuscript starts with several references on Tibet from the INDEPTH experiment, which are 20-25 years old. There are numerous, more recent references that show that partial melt is not generalized, neither in MT results, nor in seismological images. The INDEPTH results were also low spatial resolution (see comment 1 above) and were located in a major rift zone of the plateau. As a result, showing 40 km of partial melt in Figure 4c is not acceptable, even schematically.\nMinor comments\n-reference 6 relates to a volcanic study, not generalized to convergent settings over 100s of km\n-ms. Line 9-12: these are partly contradictory explanations\n-line 179-180: not clear what is meant here – the lowermost crust of Tibet is eclogitized\n-l.195: is this sentence attached to the previous, or needs to be completed?\n-extended data fig. 3: is the gray band around the gravity data represent some sort of uncertainty? How is this estimated?\n-extended data fig.4c,d: it would be interesting to show a profile of Vp/Vs as well\n-it could be explained why some figures/arguments use free-air and others Bouguer anomalies"}, "author_response": "Response to reviews:\nReviewer #1 (Remarks to the Author):\nThe present work aims to construct an 80 km deep geophysical model across the East Anatolian Plateau (EAP). To achieve this target, seismic receiver function was done and its results were integrated with interpretation of gravity data, and seismic tomographic, magnetotelluric, geothermal and geochemical models. I believe that this study will make a significant contribution to understanding the internal dynamics of the region. After these minor revisions, I believe it can be published in this journal.\nThank you for this clear recommendation and the succinct comment to the manuscript. The enclosed specific proposals for revision have helped us improve the manuscript.\nMinor comments.\n1. The presence of an intra-crustal low-velocity zone in Eastern Anatolia, or perhaps more significantly, an \"intracrustal asthenosphere,\" appears inevitable in my view. While the overall suggestion holds merit, the explanation for the existence of the low-velocity zone in the text is notably weak.\nWe unfortunately apparently have formulated ourselves unclear about the layers containing the partial melts. Our model includes the partial melts in discrete pockets within each of the layers interpreted from the geophysical data. We have modified the text throughout the manuscript to emphasise this aspect, e.g. in the abstract: “The results reveal a 20 km thick lower crustal layer and a 10 km thick mid-crustal layer, which both contain pockets of partial melt.”\nSpecifically, the assertion that the \"MBT main Bitlis thrust\" is associated with subduction lacks sufficient supporting data.\nWe do not associate the Bitlis Suture with subduction. Instead, we associate it with the subsequent continental collision as already in the manuscript and as implicitly given in the reviewer’s comment. We have deleted one sentence which might have caused the misunderstanding (at former lines 105-106): “which may be associated with the NE-ward subduction of the Neo-Tethys oceanic slab” in order to avoid possible misunderstanding.\nThe assumption that the segments observed in tomographic studies belong to the Arabian plate is entirely speculative. Recent papers, such as \"Portner et al 2018 P-wave tomography\" and \"Kounoudis et al 2020, seismic imaging of the Eastern Mediterranean mantle,\" present models with segments that cannot be reasonably inferred. As a result, the \"underplating\" assumption does not seem to have a well-supported foundation. The subduction and the subduction location of the Arabian plate should be discussed in detail. If necessary, alternative option(s) for the \"underplating\" assumption should be proposed and argued. The subduction zone from the front of the Bitlis massif should be considered in light of the results of this study.\nThe various tomographic models show differing features, which is not unusual. Our interpretation follows the interpretation by Lei and Zhao (EPSL, 2007) which indicates northward subduction from a shallow level around the Bitlis suture (Fig. R1). Based on this result, we interpret the newer tomographic models, referenced by the reviewer, with a shallower remnant of the Neotethyan slab as shown in the enclosed profile (Fig. R2), which roughly coincides with our profile.\nThe two cited references both interpret their images with the Neotethyan slab being only present within the transition zone in eastern Anatolia, which indicates that it broke off at least 30 My ago. However, the images actually show significant positive velocity anomalies above the transition zone which are not taken into account in their interpretations and, more important, their interpretation requires that the Neotethys slab (part of the Arabian plate here) would have subducted within the continental Arabian plate, which is unrealistic, whereas our interpretation matches the Bitlis suture.\nOn this basis we maintain our interpretation.\nFigure R1 from https://www.sciencedirect.com/science/article/abs/pii/S0012821X11002032 has been redacted.\nFigure R2 from https://pubs.geoscienceworld.org/gsa/geosphere/article/14/3/907/530928/Subduction‐termination‐through‐progressive‐slab has been redacted.\nal. instead interpret the slab solely located in the mantle transition zone, which would require that the trench of the subducting Arabian oceanic plate was located far inside the continental Arabian plate (illustrated by stippled yellow line); this is certainly unrealistic. Our interpretation indicates that the detached slab still is located at a depth of 150-200 km or shallower – the exact depth cannot be determined from this profile.\n2. The content of abstract looks very weak, I think it would be better if it will rewrite again.\nWe have rephrased the abstract with emphasis on the discrete appearance of the melt inclusions, see above.\n3. Please prepare a list consisting the number of each station’s measurement.\nThis list was already enclosed as Extended Data Table 1 in our original submission – now Supplementary Table 1. The reviewer may not have noticed it at the end of the Extended Data file.\nReviewer #2 (Remarks to the Author):\nBased on the new results of Receiver functon(RF), in combination with the results of regional seismic tomography, magnetotelluric, geothermal and geochemical models, this paper reveals the architecture of the upper 60 km of the East Anatolian Plateau lithosphere and emphasizes the importance of mafic underplating in generating the high attitude of the EAP, Andes and Tibetan plateaus. The seismic images are of high quality and the interpretations are mostly reasonable and reliable. This paper is well organized in general, and the results would have broad impacts for geologic society. I would recommend its possible publication in Nature Communications after moderate revisions.\nThank you for this clear recommendation and the succinct comment to the manuscript. The enclosed detailed suggestions for minor changes have helped us improve the manuscript.\nMy comments are listed below in an order of line numbers:\nLines 13-14 in abstract: Pre-collisional break-off of the Neo-Tethyan slab ------ break-off generally occurs after collision; granitoid melting ---- it is strange, please considered to revise it, e.g., generation of felsic rocks …\nWe have deleted “pre-collisional” at the break-off, which was included by mistake, and also clarified that break-off probably is later than the slab roll-back, which could be misunderstood in the earlier version. We have changed the term granitoid melt throughout the manuscript to felsic melt.\nLine 4 in the Introduction section: This depends on the timing for generating the crustal partial melt.\nWe agree as we also write it in the immediately following part of the manus. No change.\nLines 5-7: These models have not been mentioned in later discussions.\nThese references are given to introduce the study problem. Ref. 2 is already used later. We have now included ref. 8 into the later discussion where it naturally belongs.\nLine16: Reference 16 should be moved forward and expressed as 14-16 and new critical reference (e.g., Hu et al., 2016, ESR; Zhu et al., 2015, SR) should be added to the original site of 16.\nWe give separate reference to the new paper (ref. 16) because it provides the first data evidence for pre-collisional uplift in order to emphasise this aspect. The proposed references do not provide this evidence. No change.\nLines 105-106: This half sentence seems contradicts with later description. I would suggest delete it.\nThank you. The sentence could be misunderstood, now deleted, cf. comment to Reviewer 1.\nLine 117: This pattern is not clear and too simplified. The original paper by Schleiffarth et al. (2018) interpreted it as lithospheric delamination as traditionally interpreted. Can you please add a new figure to show the spatial and temporal variations of rocks since 12 Ma?\nFig. 4a (former 3a), lower panel already provides the requested information along the profile, and demonstrates the southward younging of magmatism onset while magmatism continues in the north. Schleiffarth et al. interpreted the recent pattern addressed by us by lithospheric delamination and slab tearing/rollback (my underlining). We have added a comment about this dual possibility for the interpretation. We also have introduced the new Supplementary Fig. 7 with a detailed map, as requested by the reviewer.\nLines 129-144: This part is very weak as (1) lithospheric delamination can also result in intensive addition of mafic magma to the lower crust and intra-crustal partial melting and (2) you did not consider the effects of trench retreat and tectonic shortening since 12 Ma.\nDelamination certainly can cause magmatism, but the southward younging of the magmatism is difficult to associate with delamination; therefore we favour the subduction model with roll-back and subsequent break-off, similar to earlier studies as referenced. We have added a comment on this aspect: “Alternatively, the magmatism could have been caused by delamination of the lithospheric mantle and lower crust 24, but this mechanism does not explain the clear southward younging of magmatism initiation”.\nLine 167: The seismic “bright spot” in southern Tibet is located in middle-upper crust, not middle-lower crust, given that the Tibetan crust is very think since 20 Ma (up to 70-80 km).\nWe agree, and we did write that the “bright spot” may extend downward to the top of the lower crust or be only 10 km thick; in both cases it is located in the upper to middle crust. We have added a specific comment (underlined by me): “two end-member models for the thickness of the upper to middle crustal layer with partially molten bodies have been proposed”.\nLines 177-180: Please note that the high-density, high velocity layer in southern Tibet is not mafic, but eclogitized mafic and ultramafic rocks as summarized by Zhu et al. (2022, NC; 2023, AREPS).\nWe are aware that the lower crust in Tibet may be partially metamorphosed into eclogite facies, as we also already include four earlier citations. The existence of the eclogitic lower crust has recently been disputed in part of the Lhasa block, where there is no mafic lower crust to 80 km depth (Wang et al., Nat.Comm. 2021, as already referenced). Anyway, this question is relatively peripheral to our manuscript.\nLines 180-183: I am not sure the logic here as granitoids are related to partial melting of the (ultra)mafic lower crust or fractional crystallization of parental magma, and its thickness is also related to exhumation. When did the granitoid molten form? Is it coeval with the onset of plateau uplift? Otherwise, this comparison is meaningless.\nThank you. As described earlier, we have modified the term from granitoid to felsic melt throughout the manuscript in order to avoid discussion of the exact formation process, which is not our topic.\nLine 190: Felsic melting is strange. This sentence would be changed to: but the thickness of continental crust may not be insufficient to result in extensive crustal melting generating felsic rocks.\nThank you. On reviewer’s suggestion, we have changed to “We speculate that underplating may occur, but the thickness of the continental crust may be insufficient to generate large volumes of felsic rocks by extensive crustal melting”.\nLine 195: This sentence is unnecessary.\nThank you for spotting. The punctuation was wrong. The logics has been fixed by changing the text to “The final high plateau is built by inflating the weakened, low viscosity continental crust containing pockets of partial melt, while the uplift induced regional climate change may affect the high topography and vice-versa 44,45”.\nLines 199-201: Essentially, pre-collision plateau uplift is related to isostasy of thickened crust due to mafic underplating. This has been addressed by Zhu et al. (2017).\nWe agree. We have inserted the reference.\nLines 201-203: the presence of continental crust above a subduction zone ----- do you mean a continental margin/arc? If yes, do you think the EAP is under an continental arc setting, not intraplate?\nWe mean that ocean-continent subduction is required at the onset of the plateau rise. The sentence has been modified to clarify: “Our results indicate that subduction below continental lithosphere is a necessary condition for collisional thickening by underplating, intra-crustal melting and related plateau uplift”.\nI hope these comments are helpful and looking forward to seeing the publication of this important paper.\nThank you very much; the comments certainly have helped us to improve the manuscript.\nReviewer #3 (Remarks to the Author):\nThe manuscript proposes that mafic lowermost crust and partial melt at shallower levels characterize convergent boundaries, such as in Tibet, in the Andes, and in the Eastern Anatolian Plateau, the focus of this study. The geoscientific analysis is based on seismic receiver functions, to highlight the main discontinuities in the lithosphere, and complemented with the analysis of gravity anomalies, magnetotelluric results, thermal field, and some geochemical data. The conclusions draw strong similarity between the three areas, based on which the authors generalize the prior conditions and formation processes in question for plateau uplift.\nThank you for this succinct comment to the manuscript. The enclosed specific proposals for revision have helped us improve the manuscript.\nThe main topic is very relevant and interesting to a broader community, and the multi-disciplinary approach of the authors is laudable. There are, however, also clear weaknesses, which do prevent claiming that the results are robust. These are (1) the scale and resolution of the investigations, (2) the lack of quantitative uncertainty assessment and the lack of quantitative link between most (not all) of the different observables, (3) some missing but key references. I develop these points in more detail below, and add further remarks.\nWe have made new calculations, resulting in new figures, to address the specific suggestions, as explained below.\n(1) Typical scales and resolution of the results.\nThe authors claim that the lower crustal mafic layer as well as partial melt prevails over several hundreds of kilometers, in a fairly homogeneous way. This statement cannot be supported by the data presented in this study, because the seismic stations (a total of 16, spread heterogeneously over 450 km distance, see Fig. 1) leave very large gaps along the profile (see ray coverage in Extended Data Figure 8d). Moreover, the gravity data, being derived mainly from satellite measurements, is not able to resolve any feature smaller than ca. 100 km, and it has been filtered with a Gaussian of 500 km (!). In the meantime, the MT results of prior studies (Extended Data Figure 4d) clearly shows local anomalies, and NOT a generalized layer of low resistivity. This is very similar to what is observed in more recent MT studies across the Tibetan plateau, and what seismological analysis of Hi-CLIMB stations show: local, isolated low-velocity zones, with ca. 10 km thickness and a horizontal extent no longer than 50 km. For these reasons, the results are not robust and generalization of local features seen at low spatial resolution and coverage is not convincing.\nThe reviewer nicely summarises our observations and the obtainable resolution by various methods, from which they derive the same conclusions as us: the presence of isolated partially molten bodies. Unfortunately, we must have formulated ourselves unclear about the layers containing the partial melts as both rev. 1 and 3 comment on this aspect. Our model does include the partial melts in discrete pockets within each of the layers interpreted from the seismic data, as demonstrated by the MT profile (former Suppl. Fig. 4e, now Fig. 3g). We have modified the text throughout the manuscript to emphasise this aspect, clarified it in the caption to Fig. 3 and 4, and we have moved Suppl. Fig. 4 to Fig. 3 in the main manuscript.\nTo highlight this problem with an example: in Extended Data Figure3b bottom panel, there are 8 receiver functions in the zone of interest, within the 200 km long zone with an upper crust, partial melt, D2 and D1 layers. We should see positive RF amplitudes at the bottom of the partial melt layer (bottom of low-velocity zone), and of the D2 layer (Moho). On two of the eight RFs the top converter is not visible (not visible on the 3rd from the left, and offset in depth on the 5th from the left). At the Moho level, the 4th from the left has no positive amplitude (there are two adjacent ones).\nWe agree that the resolution of receiver functions does not allow imaging of all details. However, the example discussed by the reviewer largely provide full resolution, e.g. the third RF actually shows a positive signal superimposed on the long blue (negative) signal. The two other examples show positive signals close to the interfaces, which reflect the short-wavelength variation in exact depth to the bottom of the low-velocity zones. In order to provide more solid evidence, we now include vertical velocity profiles at each seismic station derived by inversion of the receiver functions, as suggested later by reviewer 3. The results are shown as vertical profiles (new Fig. 3a, d, f, g) and as interpolated profiles (Fig. 3c, e). The inversion is illustrated in new Supplementary Fig. 3. The results of the inversion support our interpretation of pronounced low-velocity zones with some variation in the exact depth along the profile, as expected for localised partially molten bodies.\n(2) Uncertainties, simplifying assumptions.\nAlthough the study features a very nice diversity of datasets, there is limited direct link between them, and the simplifying assumptions of the RF analysis is a clear limitation. The H-kappa stacking method is very good for the first order, and I do agree with the distinction between deeper and shallower converter, however the depth estimate can be seriously off due to the need to assume a single average Vp value for the entire crust with this approach (here: Vp=6.3km/s). The uncertainty analysis does not take into account this factor. Other approaches to invert for the Vs structure with depth (e.g. with the Neighbourhood Algorithm method) would provide much better constraints, and would allow to come up with more credible quantitative depth and Vs estimates.\nWe agree about the uncertainty of depth estimates for RFs which, however, usually is most pronounced for large scale interpretation with huge variation in Vp/Vs; therefore, the exact thickness differences between the Arabian (without melt) and the Anatolian (with melt) side of the Bitlis Suture may be subject to some uncertainty due to Vp/Vs variation, which nevertheless is insignificant for our conclusions.\nReviewer 3 specifically mentions the possible importance of the choice of Vp= 6.3 km/s for the Hk stacking, which originally was addressed in the original H-k stacking paper by Zhu and Kanamori (JGR 2000). We have carried out a specific test of its importance by recalculating the Hk stacking for Vp in the interval between 6.0 and 6.6 km/s (new Suppl. Fig. 10). The results confirm the statements by Zhu and Kanamori; they are only marginally affected by the exact choice of Vp; and the result for Vp=6.3 km/s actually fits the inversion result. As also proposed by the reviewer, we have carried out an inversion for the velocity structure from the RFs (New Fig. 3 and Suppl. Fig. 3, see also above).\nThe connection between RFs and gravity using ref.27 seems to be good idea, however there is no other quantitative link between any other pair of datasets.\nA related point: it is not addressed how the distribution of various density layers would behave on the long term, on the timescale of convergent boundary evolution. If the partial melt layer is so thick, at a density of 2.6 g/cm3 and a melt fraction of ca. 7-14% (see Ext. Data Fig. 5a) which makes it mobile, how does it stay below a layer with 2.7 g/cm3 density?\nHeat flow data and any thermal model are not shown in the study, although that would likely be a very useful argument.\nThere cannot be a direct quantitative link between the geophysical and geochemical data! The link is through the processes discussed. The link between the various geophysical datasets is our geophysical model and follows from the distinct similar variation in the different physical parameters, which together provide evidence for the material content and physical state.\nRegarding the partially molten layers, we emphasise that the melt does not occur as one homogenous layer along the whole 200 km long zone, but instead it occurs in isolated bodies or pockets along the profile, which implies that this depth interval is not weak on long term. We also notice that the East Anatolian Plateau (EAP) is characterised by continuing volcanism, which shows that part of the melt actually migrates to the surface. We now have emphasised the aspect of isolated pockets throughout the manuscript, cf. comments to Reviewer 1. We do not have a thick continuous molten layer, which indeed would have affected long-term stability. These isolated melt pockets in our models cannot have such effect.\nWe agree that Heat Flow could provide useful data. We refer to Fig. 4b, where we discuss the geotherm and the formation of melts (was also in the original submission). The measured heat flow in EAP is highly heterogenous between 30 and >150 mW/m2, cf. Artemieva and Shulgin (2019), but data is absent in our study region. The published Curie Temperature Depth along our profile (Aydin et al., 2005) is generally 12-18 km in EAP, corresponding to a hot geotherm. We have added one sentence: “Heat flow data are absent along our profile, but the 12-18 km deep Curie temperature depth indicates high heat flow26”\n(3) References.\nThe manuscript starts with several references on Tibet from the INDEPTH experiment, which are 20-25 years old. There are numerous, more recent references that show that partial melt is not generalized, neither in MT results, nor in seismological images. The INDEPTH results were also low spatial resolution (see comment 1 above) and were located in a major rift zone of the plateau. As a result, showing 40 km of partial melt in Figure 4c is not acceptable, even schematically.\nWe refer to “classical” papers on Tibet of our choice, which still provide the current end-member models. We also referred to e.g. Ibarra et al. (Nat.Geo. 2023). We now have added two recent papers in the introduction part for the low-velocity zones (Huang et al., JGR, 2020) and low-resistivity zones (Jin et al., JGR, 2022) in Tibet, they are all regional research, not only along the graben zones. Journal limitations do not allow to include all relevant references.\nWe agree that a homogenous 40 km thick layer of partial melt is unrealistic and we do not propose that it exists (see our comments above), but a 40 km thick layer with patches of melt inclusion could perhaps be a possibility. However, we also find the idea of a 10 km thick zone of melt patches the most realistic, and only for completeness we include the thick layer option in our isostatic calculations because it has been widely discussed after it was published by Nelson et al. (1996). We have clarified in the caption to Fig. 5 (previous Fig. 4) and in the figure legend that the melt zone is not homogenous.\nMinor comments\n-reference 6 relates to a volcanic study, not generalized to convergent settings over 100s of km\nTrue. We have changed to Unsworth et al., Geosphere 2023.\n-ms. Line 9-12: these are partly contradictory explanations\nYes, the presented three models are different, but they have all been proposed in the cited literature. No change.\n-line 179-180: not clear what is meant here – the lowermost crust of Tibet is eclogitized\nEclogitization has been proposed for Tibet. To avoid confusion we have deleted the additional sentence on eclogitisation.\n-l.195: is this sentence attached to the previous, or needs to be completed?\nThank you. This sentence was incomplete. As it bears no immediate implication for our results, we have deleted: “although this mafic lower crust does not account for possible eclogitization in Tibet”.\n-extended data fig. 3: is the gray band around the gravity data represent some sort of uncertainty? How is this estimated?\nOne cannot provide direct estimate of the uncertainty of the gravity data. The thickness of the grey line has been chosen to make it visible below the calculated lines.\n-extended data fig.4c,d: it would be interesting to show a profile of Vp/Vs as well\nUnfortunately the Vp/Vs calculated from the Vp- and Vs- tomography models of the crust is not meaningful. This is often the case when calculating Vp/Vs for two models that have been inverted independently. Instead, we have added profiles of Vp and Vs (Fig. 3c, e) derived from our new velocity inversion of the RFs.\n-it could be explained why some figures/arguments use free-air and others Bouguer anomalies\nWe use the Free Air gravity anomaly for isostatic calculations and the Bouguer anomaly for structural interpretation according to common practise. We have added the following sentence to clarify: “Gravity modeling based on the Bouguer anomaly”. The importance of using the free air anomaly for isostasy was already described in the manuscript."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Reviewer #1 had no further comments", "reviewer_2": "Dear Editor,\nI have read the responses to Reviewers’ comments and the revised manuscript. I found that the logic, lines of evidence, interpretations, and conclusions are very clear and thus recommend its possible acceptance in its present form.\nCongratulations!", "reviewer_3": "The Authors have changed (or clarified) a key element of their earlier interpretation, namely: partial melt occurs in localized pockets in the crust. This is not a major surprise, and corresponds to a geologically more classical perception of how melt is distributed. It would be extremely useful if the Authors could provide their ideas on the typical scale (or range of scale, both horizontally and vertically) of these bodies. Unfortunately, the seismological data is still too sparse along the profile to do this in a convincing way – and this remains the main limitation of this study.\nMoreover, there is now an issue with Figure 2: the CCP migration, which is now very misleading. There is too much lateral smoothing, and one has the impression of continuous LVZ while the new interpretation of the authors is localized partial melt pocket. This figure does not correspond to the interpretation provided by the authors.\nOverall, the Authors come back from a 200-km long partially molten layer interpretation to localized, 10 km-thick pockets interpretation. This is a major change – but it seems from the seismological images that it is due to the spare spatial coverage. The discussion on the “effect on plateau uplift” is much less direct as it was with a continuous weak layer, namely: we don’t really learn new knowledge how melt and mafic underplating influences plateau uplift. I’m really doubtful about the potential impact of this study.\nRegarding my earlier (Reviewer3) comments, some of the issues remain or raise even more concern.\n(1) Typical scales\nThe answer is partially satisfactory.\nA main problem is that with localized partial melt, the sparse RFs are even less suitable to constrain the size and structural variations of the crust! This doubt is unfortunately confirmed by the new Figure 3, panel b, also panel a, which shows the actual data coverage – all the rest is interpolation!\nAlso, the resolution of the gravity anomalies has not been answered. With 500 km filtering wavelength, gravity is likely not bringing any constraint – I’d remove it completely from the study of the EAP.\n(2) Uncertainty\nThe issue with RF inversion is the linearized approach, and also the need to chose a _single_ Vp value for the _entire_ crust in the H-kappa grid search. For high-quality results with proper velocity-depth functions, I strongly encourage the authors to perform NA inversion, as mentioned in the initial review. The trans-dimensional MCMC inversion result the authors show seems too perfect, I suspect it has fallen in a local minimum, and it also misses to reduce the number of layers to something reasonable (there seem to be fixed-thickness layers). The example in Fig. S3 does not allow to interpret an LVZ as wide as the red band.\nAlso, many numbers changed for reasons that are not so clear (see lines 115-122 in the tracked changes version). How reliable is the model, then? There is no quantification.\nThe answer “One cannot provide direct estimate of the uncertainty of the gravity data.” is incorrect, but it does not affect this manuscript.\nRegarding the answers to the previous reviews, I’ve also read the communication between the authors and Reviewer 1. In my reading, the answer on the “subductions” shown in Figure R1 and R2 is not credible, as the tomographic velocity anomalies are hard to interpret, which leads to very subjective descriptions."}, "author_response": "Reviewer #1 had no further comments\n\nReviewer #2 (Remarks to the Author):\nDear Editor,\nI have read the responses to Reviewers’ comments and the revised manuscript. I found that the logic, lines of evidence, interpretations, and conclusions are very clear and thus recommend its possible acceptance in its present form.\nCongratulations!\n\nWe are pleased to notice that original Reviewer 2 finds that our manuscript has been improved to a high standard.\n\nReviewer #3 (Remarks to the Author):\nThe Authors have changed (or clarified) a key element of their earlier interpretation, namely: partial melt occurs in localized pockets in the crust. This is not a major surprise, and corresponds to a geologically more classical perception of how melt is distributed. It would be extremely useful if the Authors could provide their ideas on the typical scale (or range of scale, both horizontally and vertically) of these bodies. Unfortunately, the seismological data is still too sparse along the profile to do this in a convincing way – and this remains the main limitation of this study.\nWe have not changed our interpretation, but we have put more emphasis (clarified) on the aspect of discrete pockets of melts as we also wrote in the former rebuttal:\n“Unfortunately, we must have formulated ourselves unclear about the layers containing the partial melts as both rev. 1 and 3 comment on this aspect. Our model does include the partial melts in discrete pockets within each of the layers interpreted from the seismic data, as demonstrated by the MT profile (former Suppl. Fig. 4e, now Fig. 3g). We have modified the text throughout the manuscript to emphasise this aspect, clarified it in the caption to Fig. 3 and 4, and we have moved Suppl. Fig. 4 to Fig. 3 in the main manuscript” (from our rebuttal).\nThe MT model gives an idea about the lateral distance between the three major, individual melt bodies, which is around 70 km, but the resolution of the available data is insufficient to provide measures of their dimension (new text added to lines 79-82 and 118-129; all line numbers in version with tco). The MT model may have some leakage of the high-conductivity bodies into larger depth, and their depth extent is therefore unclear. This aspect has been added to the manuscript together with an explanation that the seismic data has insufficient resolution to identify the individual partially molten bodies but, nevertheless, the seismic data identifies two layers with relatively low velocity (new lines 309-321). New text in lines 125-129 estimates the vertical scale.\nMoreover, there is now an issue with Figure 2: the CCP migration, which is now very misleading. There is too much lateral smoothing, and one has the impression of continuous LVZ while the new interpretation of the authors is localized partial melt pocket. This figure does not correspond to the interpretation provided by the authors.\nWe did not change Fig. 2. It is identical to the figure of our original submission. Neither Reviewer 3 nor any other reviewer commented on Fig. 2 in the original submission. The new text in lines 69-70 and 124-129 explains the continuity of the low-velocity layer. Figure 2 provides key support for our model.\nOverall, the Authors come back from a 200-km long partially molten layer interpretation to localized, 10 km-thick pockets interpretation. This is a major change – but it seems from the seismological images that it is due to the spare spatial coverage. The discussion on the “effect on plateau uplift” is much less direct as it was with a continuous weak layer, namely: we don’t really learn new knowledge how melt and mafic underplating influences plateau uplift. I’m really doubtful about the potential impact of this study.\nWe did not change our interpretation, cf. our reply above. We emphasise that Reviewers 1 and 2 find the paper of interest and properly documented. The low-velocity layers may not contain partial melt everywhere, but the seismic results include a continuous layer with low velocity which indicates that the material generally is close to the solidus, cf. new lines 124-129. Together with our density estimates, these observations provide strong constraint on the isostatic uplift mechanism.\nReviewer 3 may have mistaken regional with local isostasy, when assuming that knowledge on the melt pocket geometry is needed to explain plateau uplift. Our key conclusion about the mechanism for plateau uplift is very significant and novel. We now include a new Fig. 6 to make our general conclusion clear. The reviewer’s comments do not relate to this conclusion.\nOn this background we cannot see how we possibly can modify our manuscript based on these comments, other than our new lines 69-70, 79-82, 118-129, 306-316.\nRegarding my earlier (Reviewer3) comments, some of the issues remain or raise even more concern.\n(1) Typical scales\nThe answer is partially satisfactory.\nA main problem is that with localized partial melt, the sparse RFs are even less suitable to constrain the size and structural variations of the crust! This doubt is unfortunately confirmed by the new Figure 3, panel b, also panel a, which shows the actual data coverage – all the rest is interpolation!\nWe already answered this comment in the former rebut, cf. above. We have no goal to constrain the size of local structural variations of the crust. The receiver functions cannot identify the geometry of the individual melt pockets, but they identify the overall layering, including a low-velocity depth interval where the individual melt bodies are located. The geothermal and petrological data indicate that this low-velocity layer is generally close to the solidus and therefore the material has low density, also outside the pockets with melt (lines 124-129), as further supported by our gravity modelling. The MT model identifies the partially molten pockets along the profile, as we already included in the revised manuscript and wrote in our rebut, as also appreciated by Reviewers 1 and 2. It is not the purpose of this paper to image the individual pockets of melt. We provide an integrated model of the physical parameters from which we can make a significant conclusion about the regional scale isostatic uplift mechanism.\nThe interpolation along the profile is justified because the models derived by single station inversion all are similar (lines 317-321), more detail in the comments to the following review paragraph. Therefore, our model is well constrained by consistent receiver functions and then interpolated, as in all other studies where profiles are illustrated based on broad-band seismic data acquired in 2D grids; our data sampling is denser than most other similar surveys because we combine data from two different sources, providing higher resolution. Further, our model is based on integration with tomography and MT models, which further justifies the interpolation.\nAlso, the resolution of the gravity anomalies has not been answered. With 500 km filtering wavelength, gravity is likely not bringing any constraint – I’d remove it completely from the study of the EAP.\nWe have added a comment about the resolution of the gravity data which is ~9 km horizontally and ~9 mGal for the anomalies (lines 323-325).\nThe 500 km wavelength filtered gravity data is not (and cannot be) used for the detailed interpretation of melt pockets, but is instead used for other purposes:\n1. Identification of the long-wavelength variation of Moho depth, i.e. for application of the new receiver function/gravity method by ref. 31, which documents numerically the best choice of wavelength filter. This paper is already referenced in the manuscript.\n2. Calculation of the isostatic topography in East Anatolia, Tibet and Andes, which is our main result!\n3. Resolving the density contrast across the Moho.\nFor these purposes, it is important to only consider the long-wavelength signal. We have added comments about this aspect in the main text (lines 79-82) and in Methods (lines 271-274; 323-327).\nThe reviewer apparently has misunderstood the purpose of the use of gravity here. We do not use the gravity, and it cannot be used, for estimation of the size of the individual melt pockets. The gravity is required and cannot be removed from the paper as it is essential for the interpretation and justification of the final geodynamic model (Fig. 5), as also appreciated by Reviewers 1 and 2.\n(2) Uncertainty\nThe issue with RF inversion is the linearized approach, and also the need to chose a _single_ Vp value for the _entire_ crust in the H-kappa grid search. For high-quality results with proper velocity-depth functions, I strongly encourage the authors to perform NA inversion, as mentioned in the initial review. The trans-dimensional MCMC inversion result the authors show seems too perfect, I suspect it has fallen in a local minimum, and it also misses to reduce the number of layers to something reasonable (there seem to be fixed-thickness layers). The example in Fig. S3 does not allow to interpret an LVZ as wide as the red band.\nWe now use forward modelling of the receiver functions as in the original submission, linear inversion where we show that previous inversion missed the low-velocity zone because the model did not extend sufficiently deep, and the non-linear Monte Carlo inversion algorithm (MCMC, lines 284-288 and 300-310, also present in the former resubmission). The reviewer indicates that we use linear inversion, but the MCMC inversion is non-linear, similar to the proposed NA algorithm which however has smaller vertical resolution than MCMC. There is no evidence for this reviewer statement.\nThe H-k stacking method has been applied according to state-of-the art. Reviewer 3 gave the same comment in the last review, after which we added supplementary figure 10, which demonstrates that the choice of constant Vp does not affect the results within uncertainty as already demonstrated in the original H-k paper by Zhu and Kanamori (JGR, 2000) as well as already stated in our previous rebut, but not acknowledged in the new review.\nWe already provided an uncertainty analysis of the parameter space in the revision (Suppl. Figs. 3 and 4), but the reviewer does not comment on this new addition and indicate that we neglected this aspect.\nReviewer 3 recommended us to include results from velocity inversion, which we did in the revision, but the present review does not comment on our new results except stating that the inversion is linear, despite it is not, and stating that the result is too perfect; but how can the result be too perfect? We used the standard, tested and widely applied MCMC inversion method by letting the program determine the best vertical velocity profiles for each station individually. The results show consistency between the individual stations which exactly demonstrates the robustness of our results, and which justifies the interpolation.\nThere cannot be local minimum as the reviewer claims: Consistent results were obtained for all stations along the profile, when inverted independently (lines 300-310, 317-321). The reviewer further states that the algorithm misses (fails?) to reduce the number of layers, but the results are consistent despite a large number of layers are allowed in the inversion, which emphasises the robustness of the result, as appreciated by reviewers 1 and 2. The algorithm actually first ran 300,000 inversions to define the optimum number of layers for each station, so the number of layers is not fixed but instead optimised by the algorithm. All inversion parameters are and were listed in Supplementary Table 3.\nReviewer 3 believes that the algorithm uses fixed-layer-thickness, but this is not the case. They recommend us to use another specific inversion method (NA), with the argument that our results seem to have fallen into a “local minimum”, despite Monte Carlo inversion cannot fall into local minima! We have now tested the NA method. Both the MCMC method used by us and NA are non-linearised and the results of the two methods are similar as we now show in New Supplementary Fig. 12, although the vertical resolution of the NA method, preferred by the reviewer, is much smaller than for the MCMC method.\nAlso, many numbers changed for reasons that are not so clear (see lines 115-122 in the tracked changes version). How reliable is the model, then? There is no quantification.\nThe values in lines 115-122 of the track-change version (now lines 112-117) did not change the Curie T depth, it was merely added based on a publication. Melt fraction changed from 7-14 to 3-13% because the velocity range changed from 2.90-3.10 to 2.96-3.33 km/s after the new velocity inversion of the receiver functions ended up with a slightly wider range of values than our previous forward modelling. The values are very consistent, contrary to the statement by the reviewer.\nThe answer “One cannot provide direct estimate of the uncertainty of the gravity data.” is incorrect, but it does not affect this manuscript.\nThis our answer, taken out of context, referred to the thickness of a line in one of the figures, and had no relation for uncertainty as we wrote in the previous rebut. It was made thick to be visible. We now give the uncertainties of the data (lines 325-326).\nRegarding the answers to the previous reviews, I’ve also read the communication between the authors and Reviewer 1. In my reading, the answer on the “subductions” shown in Figure R1 and R2 is not credible, as the tomographic velocity anomalies are hard to interpret, which leads to very subjective descriptions.\nReviewer 1 accepted our comment, which was basically that the interpretation of earlier tomography is very uncertain, and we agree with the reviewer that their interpretation is non-unique. However, the tomography certainly does not contradict our model!\nAs described above we have revised the manuscript substantially based on review 3. Some further revision was also introduced in order to improve clarity of our arguments, including the addition of new Fig. 6."}]

Decision letter