Correlated crustal and mantle melting documents proto-Tibetan Plateau growth

ABSTRACT The mechanism that causes the rapid uplift and active magmatism of the Hoh-Xil Basin in the northern Tibetan Plateau and hence the outward growth of the proto-plateau is highly debated, more specifically, over the relationship between deep dynamics and surface uplift. Until recently the Hoh-Xil Basin remained uncovered by seismic networks due to inaccessibility. Here, based on linear seismic arrays across the Hoh-Xil Basin, we present a three-dimensional S-wave velocity (VS) model of the crust and uppermost mantle structure beneath the Tibetan Plateau from ambient noise tomography. This model exhibits a widespread partially molten crust in the northern Tibetan Plateau but only isolated pockets in the south manifested as low-VS anomalies in the middle crust. The spatial correlation of the widespread low-VS anomalies with strong uppermost mantle low-VS anomalies and young exposed magmatic rocks in the Hoh-Xil Basin suggests that the plateau grew through lithospheric mantle removal and its driven magmatism.

~50 km (Fig. 1a and Supplementary Fig. 1a).As deployment periods of the temporary seismic networks used in this study did not substantially overlap in time, surface wave signals at long periods extracted from ambient noise interferometry would be limited by the short interstation distances within each sub-network, and consequent loss of resolution in the deep crust and mantle.
Thus, to improve the data coverage in the northern Tibetan Plateau, we also processed 38 permanent broadband seismic stations from the China National Seismic Network (CDSN) which overlap in time with almost all temporary networks (Supplementary Figs.1a and 2).
We extract Rayleigh wave signals from ambient noise similarly to previous studies [2,3].
After removing the instrument response, resampling (1 Hz), low-pass filtering (corner 0.01 Hz), applying time-domain normalization and spectral whitening, we calculate cross-correlations for all simultaneously recorded station pairs using 1-hour length segments.The hourly crosscorrelations in a day are saved and linearly stacked to a daily one.Then, we average the positive and negative lags and linearly stack all daily cross-correlations over the entire operating period.
In total, we obtain stacked cross-correlations for more than 22000 station pairs.Significant surface wave signals with velocities of 2-4 km/s can be recognized in the stacked crosscorrelations at different periods (Supplementary Fig. 1b).
We pick group and phase velocity dispersions of Rayleigh wave by using the Automatic Frequency-Time ANalysis (AFTAN) package.After this automatic procedure, we manually checked all dispersion curves and retained dispersions that are continuous over a wide range of periods.We also impose a minimum interstation distance control following previous studies [2,3], and the threshold is set to be greater than two wavelengths.For both group velocity and phase velocity, we finally obtained more than 13000 dispersion curves (Supplementary Fig. 3c).More than 10000 dispersion measurements are available at each period between 8-30 s, decrease gradually in the longer (35-65) periods.Less than 1500 valid measurements are available at long periods of 70-100 s (Supplementary Fig. 3c).Thus, we only use measurements in the periods of 6-65 s in the inversion of 3-D Vs model.These group and phase velocity measurements cover the study region well (Supplementary Fig. 4).At short periods, the measurements along ray paths across basins, i.e., Qaidam Basin in the northeastern Tibetan Plateau and Hoh-Xil Basin in North Tibet, show significantly low velocities.At long periods, the velocity is high across the Qaidam Basin but still low across the Hoh-Xil Basin.

Direct inversion of 3-D V S model
We apply the direct inversion method [4] to construct the three-dimensional (3-D) V S model.
The key to the inversion problem is to seek a model that minimizes the travel time residuals, that is, the deviation of the observed travel times ( ) from the ones calculated from the model ( ).

/ 22
Accordingly, this can be expressed in the classical matrix form and can be solved by several iterations to minimize the objective function as following: where is a model smoothing operator, which is set to the second order of the spatial derivative operator, and is the weighting parameter balancing data residual term ‖ ‖ and model regularization ‖ ‖ .The sensitivity matrix can be calculated from the initial model, which is updated after each iteration.
Here, we mesh the study region as 24-layer grids with the depth interval of 5 km at 0-80-km depth, 10 km at 90-100-km depth, 15 km at 115-145-km depth, and 20 km at 160-180-km depth (Supplementary Fig. 5c).We set the V S in the initial model as the 1-D Vs model inverted from the average group and phase velocity dispersion curves (Supplementary Figs.3a and 3b) using the Computer Programs in Seismology package [5] (Supplementary Fig. 5c).The meshed grids have horizontal intervals of 0.4° in latitude and longitude, which was chosen based on available ray coverage (Supplementary Fig. 5d) and tests with different grid sizes (Supplementary Fig. 6a).
We can find a decrease of the residual when the grid size decreases from 0.5° to 0.4°, but the residual decreases only slightly when the grid size further decreases to 0.3° and even increases when using large weighting parameters (Supplementary Fig. 6a).Different weighting parameter ( ) were also tested and 70 is chosen for the preferred model, as it provides an optimal trade-off between model norm and data residual term (Supplementary Fig. 6a).With our meshed space and preferred weighting parameter, the group and phase travel time residuals show a Gaussian type distribution centered near zero after 10 iterations (Supplementary Fig. 6b).We compare the inversion with the optimal weighting parameter with those with smaller or larger weighting parameters (Supplementary Figs.6c-e).Anomaly patterns are very similar for all these cases, demonstrating that the inversion is stable with regard to reasonable variations of the weighting parameter.

Resolution tests
Bootstrap tests are adopted to estimate the uncertainty from data errors and initial model [6].
We perform the bootstrap procedure by resampling with replacement the dispersion curve set for 50 times and then inverting for velocity models according the aforementioned procedure.The initial model for each inversion is set as the V S randomly perturbing each layer of the average 1-D model with a standard deviation of 0.15 km/s.The average values of the output Vs models are similar to the final model of this study with small uncertainties (standard deviation < 0.05 km/s) in regions with good data constraints (Supplementary Fig. 7).
We then perform a series of checkerboard tests to estimate the resolution ability of used data and method (Supplementary Fig. 8).The input models, constructed by 5% maximum perturbation relative to the average 1-D V S model, have various scales in horizontal and separated patterns in depths.We synthesize group and phase travel times for the paths same as in the real dataset.
According to the measurement error of Rayleigh wave dispersion estimated in previous studies [2], we add 1% random noise to the synthetic dataset and invert them following the aforementioned procedure.For small anomalies (0.8 °), the patterns can be recovered at shallow depths (< 80 km) but with some smearing at larger depths (> 80 km).For larger anomalies (1.2 ° and 1.6 °), the entire depth range of the model can be reasonably resolved, except some underrecovery of amplitudes at larger depths (> 80 km).Considering the coverage of ray paths (Supplementary Figs. 4 and 5d) and the uncertainty estimated in bootstrap tests (Supplementary Fig. 7), we exclude the edge region that exists smearing in checkerboard tests (Supplementary Fig. 8) in the following discussion, which is delineated by black frames in Figs. 2 and 4.
To evaluate the resolution capability of the dataset and method with respect to the connectivity of crustal LVZs in the Tibetan Plateau, we construct the synthetic tests for connected and segmented low-V S anomalies at depths of 25-45 km (Supplementary Fig. 9).The input model is meshed onto the same grid as that in the direct inversion procedure.Then, we synthesize group and phase travel times for the paths same as in the real dataset with added 1% random noise and invert them following the same direct inversion procedure as for the real data.
Generally, the absolute velocities of the output models are smoother than the inputs, and the strength of crustal low-velocity anomaly is under-recovered.However, the horizontal and indepth distribution of LVZs can be acceptably reconstructed in the output models.

Estimating the melt fractions of crust and uppermost mantle
The direct effects of partial melting on seismic velocity have been studied over a wide range of pore shapes, and considered to be dependent on the contiguity of grain boundaries, which is a function of dihedral angle in an equilibrium geometry [7].The contiguity of grain boundaries is thus not solely a function of the melt fraction but also strongly dependent on the composition of the melt-or fluid-bearing rock system, and was experimentally derived for various compositions [8].Here, combining the theoretical seismic velocity-contiguity relationships [7] and the experimental contiguity-melt fraction relationships for variable systems [8], we can obtain the relationship between V S and melt fraction to estimate the melt fractions of crust and uppermost mantle in the Tibetan Plateau (Fig. 4).
To estimate the melt fractions from the V S -melt fraction relationship, a reference V S at a specific depth is required.Considering the composition inferred from the xenoliths and the radial anisotropy induced by the mica lattice, Hacker et al. [9] concluded that the lowest vertically polarized V S is ~3.35 km/s at a depth of 25-35 km.Thus, for V S lower than ~3.35 km/s, partial melting should be considered.Using this value as the reference V S , we can estimate the melt fractions at the depth of 30 km which is almost the center of the crustal LVZs imaged in this study (Fig. 2).Following the V S -melt fraction relationship of the crustal mineral system, our results imply the presence of widespread partially melted crust in North Tibet with the highest inferred melt fraction of ~8% at 30-km depth beneath the western Hoh-Xil Basin (Fig. 4b).
Referring to a typical vertically polarized V S of ~4.5 km/s in the uppermost mantle [10], we can infer that the uppermost mantle beneath the partially melted crust in North Tibet contains ~4-6% melt from the low-velocity anomalies at depth of 100 km (Figs.4c and 4d).In contrast, in South Tibet, the partially melted crust is only sporadically distributed and with lower melt fractions.In the uppermost mantle, evidence for partial melting is only evident near the YGR (Figs. 4b-d).It should be noted that the reference V S should be variable due to the heterogeneity in the crust and upper mantle, hence there are some uncertainties in the estimation of melt fraction when using the same reference V S for the whole Tibetan Plateau.Here, we use different reference V S to calculate the melt fraction at 30-km depth (Supplementary Fig. 12).The amplitude of melt fractions does vary with the reference V S , but the patterns of partially molten crust are almost similar by using different reference V S .These tests indicate that the partially molten crust widespread in the northern Tibetan Plateau but sparsely isolated in the south revealed in this study is a reliable pattern.