Introduction

The mantle transition zone (MTZ) bounded by the 410 and 660 km seismic discontinuities (hereafter, the 410 and 660) connects the upper and lower mantles and plays an important role in slab subduction and mantle convection. Experimental and theoretical mineral studies indicate that a partially molten low-velocity layer (LVL) atop the MTZ can be expected, with water and/or other volatiles (i.e. CO2)1,2,3,4. Such LVLs have been extensively observed against different tectonic backgrounds around the world5,6,7,8,9. In oceanic subduction zones, e.g. the present circum-Pacific subduction zone6,7,10,11, the ancient Asian oceanic subduction zone5,12,13 and the Tethyan oceanic subduction zone14,15, these LVLs are widely explained as partial melting due to hydrous and/or other volatile materials upwelling from the MTZ triggered by slab subduction on a regional or global scale5,7,10,11,12,13,14,15. Although numerical modelling shows that this process can be expected on a regional4 or global1 scale, direct observational evidence to construct the relationship between such an LVL and a subducted slab is still absent. In continental subduction zones, LVL observations are relatively few. We also know little about how the LVL affects slab subduction processes.

Generally, in the Wilson cycle, after the closure of an ocean, oceanic subduction transforms into continental collision, followed by continental subduction. The previously subducted oceanic slab that is recycled into the mantle is expected to exert a strong influence on the physical and chemical properties of the deep mantle, in turn influencing the subsequent continental subduction. However, ongoing continental subduction in the deep mantle is rare and thus hinders our understanding of the related processes.

After the subduction and closure of the former Tethyan oceans, the Indian and Eurasian plates collided and formed the Tibetan Plateau and Himalayan mountains during the Cenozoic era. These landforms comprise the Earth’s highest plateau and active continental collision zone. The Hindu Kush and Pamir are located west of the Himalayas and are geologically named the western Himalayan syntaxis (Fig. 1). Interestingly, the Hindu Kush–Pamir mantle exhibits vigorous intermediate-depth (75–300 km) seismicity16,17,18, characterised by ~30 recurring large earthquakes of magnitudes (M) ≥ 7 with focal depths larger than 180 km since 190019. Ongoing subducted slab stretching, tearing and break-off have been suggested to create sufficiently high strain rates to trigger seismicity20,21,22,23,24, which likely results from the subducted buoyant Indian continental crust being pulled to a depth of ~160 km together with the slab21. Seismic observations show that strong high-velocity anomalies possibly representing the subducted slab have detached into the MTZ in the Hindu Kush16,20,21,22,25,26. This rare active continental subduction provides a good opportunity to detect the interaction between the deep upper mantle and slab subduction in continental collision zones.

Fig. 1: Map of tectonic setting and seismic data used.
figure 1

a Tectonic setting of the study area. The white line marks the subducted front of the Indian slab16,65. The purple lines represent the Tethyan sutures65,66. The yellow solid line indicates the boundary between the Indian and Eurasian plates67. The grey arrow shows the direction of Indian plate movement relative to stable Eurasia. The black arrow marks the location of the cross-section shown in (b). b P-wave velocity anomalies from the seismic tomography21. The black circles are the projections of intermediate–deep earthquakes21 onto the cross-section within a distance of 0.5°. c Earthquakes and seismic stations used. The beach balls represent focal mechanisms from the global centroid moment tensor (gCMT)68 for the two earthquakes. The black and blue triangles mark the locations of the seismic stations of the Xinjiang network and the Incorporated Research Institutions for Seismology (IRIS), respectively. The coloured dots and grey lines represent the surface projections of the turning and reflection points, and ray paths of triplicated waves for events 20110118 (red and blue) and 20140925 (brown). The solid and dashed black rings outline the main sampling areas around 410 and 660, respectively. The two black arrows mark the locations of the cross-sections shown in Fig. 2a. The thick purple lines represent the Tethyan sutures67.

Full size image

The sparse seismic station distribution in the Hindu Kush–Pamir limits our understanding of the detailed structure around the MTZ. Seismic triplications produced by the 410 and the 660 of the MTZ are at epicentral distances of 10–30°. They include direct seismic waves that propagate above the 410 (branch AB), in the MTZ (branch CD) and below the 660 (branch EF), and reflected waves from the 410 (branch BC) and the 660 (branch DE)5 (Supplementary Fig. S1). These branches are sensitive to structure perturbations around the MTZ (see the Methods section). Perturbed velocities can be obtained at the turning and reflected points of seismic triplications because of the highest sensitivity at these points7,27. Thus, seismic triplications are ideal for constraining the MTZ structure in the region with sparse stations.

Triplications from the earthquake configuration in southwestern Pakistan and the seismic stations of the Xinjiang network28 around the Tianshan mountains (Fig. 1c) are in the range of epicentral distances of 10–30° and can be used to sample the MTZ of the Hindu Kush–Pamir region. Here, two intermediate-depth earthquakes of Mw ≥ 5.5 (Supplementary Table S1) with simple source time functions are selected. The focal depths for the earthquakes are redetermined by modelling separate teleseismic P and sP waveforms29 (Supplementary Fig. S2). The observed triplicated waveforms are preprocessed to improve the signal-to-noise ratios (SNRs) (see the Methods section). Based on the variation in the differential times of the observed triplications, we divide the study region into three subregions, namely, A, B and C, from NW to SE (Fig. 1c).

We construct an initial model from the IASP9130 model for the mantle and the CRUST1.031 model for the crust (see the Methods section). One-dimensional (1-D) synthetic waveforms and travel time curves are calculated based on the orthonormal propagator algorithm32 and the dynamic raytracing33, respectively. Relative to the initial model, there should be an LVL atop the MTZ and a depressed 660 and/or high-velocity anomaly (HVA) in the MTZ in the Hindu Kush–Pamir region (see the Methods section). To further quantitatively constrain the velocity perturbations, we define the LVL using dup_LVL for the depth of its upper interface, d410 for its lower interface, and dVp_LVL for the P-wave velocity drop. Similarly, we define the parameters for the HVA as dup_HVA, d660 and dVp_HVA. We do not consider a gradual interface because the triplicated seismic waveforms in the frequency band of 0.02–0.5 Hz are only weakly sensitive to the sharpness of the upper and lower interfaces11. We adopt the iterative grid search method we developed15 to match the synthetic waveforms with observed triplicated P and sP waveforms to find the optimal LVL and HVA models (see the Methods section). The observations provide strong seismic evidence documenting upper mantle melting triggered by slab subduction and offer new insights into slab subduction-related processes in an intracontinental setting.

Results

Figures S5 and S6 show the grid searches for the iterative processes to find the optimal LVL and HVA models by simultaneously fitting the triplicated P and sP waveforms. The optimal models for each iteration (Supplementary Tables S2 and S3) are obtained from the convergent points with the maximum cross-correlation coefficients (CCs) between the observed and synthetic waveforms. Relative to those in the initial iteration, the CCs for the optimal models obtained from the final iteration are markedly increased. The optimal LVL models (dup_LVL, dVp_LVL and d410) from the final iteration for sampling subregions A, B and C are (357 km, −4.2%, 406 km), (347 km, −3.3%, 415 km) and (338 km, −3.0%, 426 km), respectively. No HVAs are present in the MTZ, i.e. the MTZ velocities are normal, and the corresponding depths of the 660 are 687, 682 and 684 km for subregions A, B and C, respectively.

The upper interface depth and velocity drop of the LVL can produce similar differential travel times between the AB and BC/CD branches12 but different amplitude ratios, which contribute to resolving the trade-off between them. Thus, the focal mechanism is important because the radiation pattern affects the amplitudes of the AB, BC and CD branches. Moreover, errors in the focal depth influence the differential travel times between the OA and OC branches11. The two earthquakes considered here have relocated focal depths of 71 km and 51 km and different focal mechanisms of reverse and strike faults. To further estimate the accuracies of the focal depths and focal mechanism solutions (strike, dip, and slip) for these two earthquakes, we perform tests to analyse how the results are affected by reasonable errors in the source parameters (Supplementary Figs. S7, S8 and S9). The resulting optimal models change little with the perturbations of different focal parameters (Supplementary Tables S4, S5 and S6). Except for the additional 10° dip in subregion B for event 20110118, the CCs for the optimal models with the changed focal parameters are smaller than those with the original focal parameters for the two earthquakes. These results suggest that the relocated focal depths and focal mechanism solutions and their associated optimal models are acceptable.

Structural heterogeneity along the ray paths of triplicated waves may introduce errors into the optimal models in the 1-D waveform modelling method. To test the possible effects of structural heterogeneity near sources, we perform a grid-search for the optimal models by fitting the triplicated P and sP waveforms separately in subregions A, B and C based on the previously obtained results (Supplementary Figs. S10, S11 and S12). Although the P and sP paths near the sources are distinctly different, the resulting optimal models are similar (Supplementary Tables S7 and S8). This observation indicates that the optimal models are independent of the possible existing structural heterogeneity near sources.

Previous seismic tomographic models16,20,21,22,25 consistently show an HVA in the upper mantle and MTZ. Recent finite-frequency seismic tomography21 using more regional data provides higher-resolution HVAs, possibly representing a subducted slab in the deep mantle beneath subregions A, B and C (Supplementary Fig. S13). In accordance with the tomographic result, we design a simple two-dimensional (2-D) slab model at depths of 250 to 600 km with a width of 180 km and a velocity increment of 3.0% modified based on the obtained 1-D optimal LVL and HVA models for subregion A (Supplementary Fig. S14a). For subregions B and C, the slab models are at depths of 200 to 650 km with a width of 180 km and a velocity increment of 1.0%, and at depths of 300 to 650 km with a width of 100 km and a velocity increment of 1.5% (Supplementary Fig. S14b, c). The distance between the seismic source and the left boundary of the slab model is 1000 km for subregion A and 900 km for subregions B and C. Using the 2-D synthetic waveform-calculating software package SPECFEM2D34 based on the spectral element method, we calculate synthetic waveforms based on the simplified 2-D slab models and the 1-D LVL and HVA models (Supplementary Fig. S15). The resulting triplicated waveforms are similar between the 2-D and 1-D models, with high CCs of 0.9786, 0.9989 and 0.9966 for subregions A, B and C, respectively. This could indicate that the slabs have limited effects on the triplicated waveforms.

To quantitatively estimate the effects of the slabs, we perform a grid-search for the optimal LVL and HVA model parameters by fitting the 2-D synthetic triplicated waveforms calculated with the slab models using the 1-D synthetic waveforms without the slab models (Supplementary Fig. S16). The resulting optimal models (dup_LVL, dvp_LVL, d410, d660) for subregions A, B and C are (362 km, −3.8%, 404 km, 694 km), (349 km, −3.2%, 413 km, 685 km) and (342 km, −3.0%, 423 km, 685 km), respectively. Relative to the forward models (Supplementary Tables S3), the slabs produce errors of (+5 km, +0.4%, −2 km, +7 km), (+2 km, +0.1%, −2 km, +3 km) and (+4 km, +0.0%, −3 km, +1 km) for subregions A, B and C, respectively.

To test the effects of shallow structural heterogeneity near seismic stations, we construct a low-velocity anomaly at depths of 0 to 150 km with a velocity drop of 3.0% and a width of 150 km modified based on a 1-D LVL and HVA model (350 km, −4.0%, 400 km, 680 km) (Supplementary Fig. S17). Synthetic waveforms are calculated with and without this low-velocity anomaly (Supplementary Fig. S18). The resulting triplicated waveforms are very similar, with high CCs of 0.9984, which indicates that the low-velocity anomaly at shallow depths has limited effects on the triplicated waveforms. The resulting optimal model we found through the 1-D grid search (Supplementary Fig. S19) is (348 km, −4.1%, 400 km, 680 km), which indicates that the effect of shallow structure near seismic stations on the final results is negligible.

To estimate the possible effects of the waveform data, we randomly sample seismic stations to perform a grid search for the optimal model parameters. The number of stations for each round of sampling is set to be more than 80% of all stations. The total number of sampling rounds is 1000. Thus, we obtain 1000 optimal models (Supplementary Fig. S20). The resulting possible optimal model parameters are listed in Supplementary Table S9. After correcting the effects of the slabs analysed above, the final optimal LVL and HVA models are (351 ± 5 km, −4.7 ± 0.3%, 407 ± 3 km, 681 ± 4 km), (342 ± 3 km, −3.5 ± 0.2%, 416 ± 2 km, 679 ± 4 km) and (334 ± 6 km, −2.8 ± 0.2%, 428 ± 3 km, 682 ± 3 km) for subregions A, B and C, respectively (Fig. 2a).

Fig. 2: Comparison of observed and synthetic triplicated P and sP waveforms calculated with the optimal structure models.
figure 2

a The P-wave velocity models. The red solid and white dashed lines correspond to the optimal and IASP91 models. The background colour represents P-wave velocity anomalies as determined from the seismic tomography21. The velocities outside and within the yellow dashed lines are obtained from cross sections I-I’ and II-II’ (shown in Fig. 1c), respectively. b–d Observed (black lines) and synthetic (red lines) triplicated P and sP waveforms for subregions A, B and C, respectively. The blue dashed and green solid lines are the travel time curves calculated using the IASP91 and optimal models, respectively. The cyan lines mark the time windows for waveform cross-correlation.

Full size image

Figure 2b, c and d shows the observed and corresponding synthetic waveforms calculated using the optimal LVL and HVA models. Compared to the initial model, the optimal LVL and HVA models show improved similarity to the observations. The delayed and extended BC branches for both the P and sP waveforms are well-matched by the optimal LVL models. The differential times between the OA and OC/BC branches are more sensitive to the depth of the 410 than the upper interface depth and velocity drop of the LVL14, which markedly increase from subregion A to subregion C for both P and sP triplications and are well simulated by the optimal depths of the 410. Similarly, the observed EF and DE branches can be well matched by using the optimal depressed 660.

Discussion

Observational evidence for upper mantle melt related to slab subduction

We detected three LVLs in subregions A, B and C with obvious lateral structural variation. The LVLs are present underneath the subducted Indian slab, which extends into the MTZ in the Hindu Kush region16,20,21,22,25. The subregion A is very close to the subducted Indian slab, whereas the subregion C is approximately 200–300 km away from the subducted Indian slab (Fig. 1c or Fig. 2a). The LVLs are characterised by thicknesses of 56–94 km and P-wave velocity drops of −2.8 to −4.7%. This cannot be explained by the phase transition of serpentine into dense hydrous magnesium silicate35, nor can it be explained by the phase transformation of olivine36. Similar LVLs are found worldwide5,6,7,8,9,10,11,12,13,14,15, and these results are usually explained as being caused by partial melting induced by water and/or other volatiles (i.e. CO2). The LVLs in the three subregions show an obvious negative correlation between the velocity drops and the lower-interface depths (Supplementary Fig. S21), supporting the conclusion that the LVLs are mainly the result of melt, rather than the temperature. The observed LVLs can be explained by a small amount of water-rich melt completely wetting the grain boundaries37.

The melts formed at the 410 are denser than the surrounding solid minerals in the upper mantle but lighter than the transition zone minerals1; consequently, they will accumulate and become trapped above the 410. Therefore, we infer that the lower interface of the LVL is the 410, and the topography is sensitive to temperature variation because of the mineral phase transition of olivine and wadsleyite with a positive Clapeyron slope38. A low temperature uplifts the 410, and a high temperature depresses the 410. Compared with the IASP91 model, the depth of the 410 in our observations deepens by 21 ± 6 km away from the downwelling cold slab over ~300 km. Considering the Clapeyron slope of 4.0 MPa/K39 and the mantle pressure gradient of 35 MPa/km40, a pure temperature variation of 184 (−53/+53) K is required. This large gradient in the temperature is mainly the result of the transition from low to high temperature induced by regional mantle convection triggered by the subducted slab. The low-temperature region, closest to the subducted slab, is dominated by the cold downwelling flow controlled by the subducted slab. The high-temperature region is the result of upwelling hot hydrous MTZ materials triggered by the subducted slab, as expected from numerical modelling4,41. This kind of fluctuation in the 410 provides direct observational evidence that an LVL atop the 410 can be generated by a subducted slab on a regional scale.

Numerical modelling4,42 indicates that there should be an obvious gap between the subducted slab and its resulting LVL. However, the LVLs that we observed are very close to the subducted slab and lack an obvious gap. The hydrous partially molten LVL is likely to have formed before the collision of the Indian and Eurasian plates in the pre-Cenozoic era due to regional convection induced by long-term subduction of the Palaeo-, Meso- and Neo-Tethyan oceanic plates5 and/or global diffuse upwelling1. The subsequently subducted Indian slab penetrated the LVL, leading to the absence of a gap between the subducted slab and the LVL.

Hydrous MTZ in the western Himalayan syntaxis

The existence of the partially melted LVL can be an indicator that the MTZ is rich in water in the Hindu Kush7. In addition, the depressed 660 and normal velocity MTZ that we observed in the Pamir also suggest a hydrous MTZ. With respect to the IASP91 model, the depth of the 660 descended by 19–22 km with little lateral variation. The depth of the 660 is more sensitive to temperature than chemical composition, as it mainly originates from phase transitions of ringwoodite to bridgmanite and ferropericlase with a negative Clapeyron slope43. Possible low temperatures from the subducted cold slab to the southwest and/or the disappearance of the Tethyan slabs could depress the 660, but they would also produce high seismic velocities. However, we observed a normal velocity in the MTZ, which is consistent with previous seismic tomographic images16,20,21,22. The subducted Indian and Tethyan slabs could also release rich water and store in the MTZ because of the high water solubility of its dominant minerals, wadsleyite and ringwoodite. In the central part of the front of the Indian continental subduction zone, i.e. the Qiangtang terrane, water in the MTZ is also indicated by various seismic observations15,44. Water can depress the 66045,46 and reduce the seismic velocity47. Thus, the opposing action of the low temperature and the existence of water may be the cause of the normal MTZ velocity and the greatly depressed 660.

The effects of melt and water in the deep mantle on slab subduction

Our observations show that the upper interface of the LVL becomes deeper (from 334 km to 351 km) and the velocity drop becomes larger (from −2.8% to −4.7%) with decreasing distance from the subducted slab. The fluctuation of the upper interface is consistent with numerical simulations41,42. In this case, the LVL is assumed to have a considerably lower viscosity than the ambient mantle. The hydrous melt in the LVL can wet the solid olivine grain boundaries to produce a low-viscosity layer48. Moreover, the existence of water in the MTZ can also reduce the viscosity of nominally anhydrous minerals, leading to a lower viscosity around the MTZ than in the upper and lower mantles49,50.

The evolution of a subducted slab is strongly influenced by the viscosity of the mantle51. When the present subducted slab penetrates through the ancient subduction-induced low-viscosity LVL and into the MTZ, the low-viscosity LVL is likely to produce relatively weak resistance against slab movement. In contrast, the asthenosphere, with a relatively high viscosity above the LVL, can be more easily entrained into the deeper mantle through the subducted slab. This process may lead to compression of the LVL, which is consistent with our observations of the deepened upper interface and the increased velocity drop of the LVL close to the subducted slab.

It is usually believed that the difference in the subduction rate of the Indian slab between the shallow and deep domains induces slab tearing and break-off in the middle deep domain, as revealed by intermediate-depth earthquakes and seismic tomographic images16,20,21,22. At shallow depths, the subducted buoyant continental crust, which extends to a depth of 160 km, is thought to resist the downwards movement of the underlying lithospheric mantle, causing a slow displacement rate. However, in the deep domain, the influence of the viscosity change in the mantle is less considered in slab subduction. Although numerical simulations indicate that slab break-off with brittle failure can occur, this phenomenon depends on many factors, i.e. the plate convergence rate, slab rheology and three-dimensional geometry, and thus may result in different evolutionary patterns of slab subduction52,53. The decreased viscosity of the LVL and MTZ may make the subducted slab move readily into the MTZ, thereby more easily leading to a large displacement rate than a normal viscosity would. This behaviour, combined with the resistance of the shallow buoyant crust, seems to be more likely to produce sufficiently high strain rates to induce the numerous recurring large earthquakes observed in the past one hundred years. Therefore, the rheology of the LVL and MTZ is likely a crucial parameter constraining the speed of slab detachment and thus should be considered in numerical modelling.

Geodynamic implication for continental slab subduction

After the closure of the Tethyan oceanic plate, the Indian subcontinent collided with and subducted beneath the Eurasian continent during the Cenozoic era54. Seismic tomography studies show high-velocity anomalies under the present Indian craton, which are separated from the shallower velocity anomaly layers and inferred to as subducted Tethyan oceanic slabs55. Geological studies suggest that there are no oceanic slabs in the Hindu Kush since the Late Cretaceous56,57. The shortening of the collision system at the western syntaxis is at least 670 km58. Although there are necking or thinning structures for the subducted lithospheric mantle at the depth range of 200–250 km16,20,21,22, the intermediate-depth seismicity is widely spread until 300 km depth in the Hindu Kush region21,59. All these observations support the idea that the continental lithospheric mantle of the greater India can subduct into the depth at least up to the MTZ.

In accordance with the above analysis, we describe a slab subduction-related process in a continental subduction zone (Fig. 3). A hydrous MTZ and a low-viscosity partially molten LVL generated by ancient oceanic slab subduction likely have formed before the subduction of the continental slab. When the subducted continental lithospheric mantle penetrates into the existing LVL, it may readily pass through the LVL and enter the MTZ because of the lubricating action of the melt. In contrast, the shallow slab may have difficulty moving downwards due to the resistance of the subducted buoyant continental crust coupled with its lithospheric mantle. This differential slab movement between the shallow and deep domains is likely to induce strong slab deformation with stretching, tearing and eventual break-off in the middle domain, accounting for frequently recurring large intermediate-depth earthquakes in an intracontinental setting. Part of the subducted slab in the MTZ disturbs hydrous (and/or other volatile) MTZ materials and causes them to upwell into the upper mantle in its vicinity, inducing neo-melting in the LVL. Our work cannot quantify how much the low-viscosity layer materials has contributed to slab subduction, but our seismic observations indicate that successful models must consider the reduced viscosity in the deep mantle.

Fig. 3: Conceptual geodynamic diagram illustrating the slab subduction-related process in a continental subduction zone.
figure 3

The stars indicate intermediate-depth earthquakes resulting from slab stretching, tearing or break-off. The orange regions delineated by solid and dashed black lines represent observed and inferred partially molten low-velocity layers, respectively.

Full size image

Methods

Seismic triplication sampling of the MTZ

Because the ray paths for the triplicated waves near the foci and receivers are similar, the differential travel times and amplitude ratios between different branches of the triplicated waves are mainly related to the radial velocity variations near the 410 or the 6605,11,44,60,61. Direct P waves and the depth phase, sP, can both generate clear triplicated waveforms5. The ray paths near the source for sP triplication are distinctly different from those for P triplication (Supplementary Fig. S1). Joint use of the P and sP triplications has recently been applied to reduce the influence of possible heterogeneity near the source5.

Preprocessing of triplicated seismic waveforms

We apply second-order Butterworth filters with corner frequencies of 0.02–0.5 Hz (i.e. periods of 2−50 s) to the observed triplicated waveforms. This filter can detect the sharpest triplicated P and sP pulses with the highest SNRs (Supplementary Fig. S3), making it possible to capture the main frequency of ~3 s at an acceptable computing cost for the waveform calculations. We select the P and sP waveforms after the removal of instrumental responses with SNRs larger than 5.0.

Selection of the reference Earth model

We calculate synthetic waveforms for event 20110118 using three 1-D reference Earth models with different upper mantle and MTZ structures: PREM62, IASP9130 and STW10563 (Supplementary Fig. S4). Their crustal structures are replaced by the CRUST1.031 model because the Moho in our study region is much deeper than the global average value. The CCs between the synthetic and observed waveforms are calculated to evaluate the waveform-fitting quality. To reduce the effects from complex coda waves, the time window for cross-correlation is varied with the differential travel times between the AB and BC/CD branches, i.e. 1 s before the synthetic travel times of the AB phase and 3 s after the observed peak times of the BC phase. The average CC for IASP91 is 0.78, larger than the values of 0.74 and 0.75 for PREM and STW105, respectively. Similar to related former studies in the central part15,60,61, northern margin5, and southeastern margin14 of the Tibetan Plateau, we thus select the modified IASP91 model as our initial model, with the Q parameter from PREM model.

Seismic signal for the velocity anomaly

Compared with the IASP91 model predictions, the observed differential times between the AB and BC/CD branches are longer, and the epicentral distance of branch BC is also longer for both the P and sP waveforms (Supplementary Fig. S4b). These observations indicate that there is an LVL atop the MTZ5,14. In addition, the differential times between the O’F and O’D branches decrease with increasing epicentral distances from 23° to 28°. Such decreases can be produced by a depressed 660 and/or HVAs in the MTZ14,15. These observed features related to the 410 and 660 are similar to those in the Qiangtang terrane in the central Tibetan Plateau15. The optimal model in the Qiangtang terrane is characterised by an LVL atop the MTZ and an HVA in the lower MTZ, with a depressed 410 and 660. The CC between the observed and synthetic waveforms (Supplementary Fig. S4e) calculated with the optimal Qiangtang model is improved to 0.81, larger than those of IASP91, PREM and STW105 models.

Iterative grid search

As the triplicated seismic branches AB/BC for the LVL and DE/EF for the HVA arrive at almost the same time for epicentral distances of 19–22° (Supplementary Fig. S1b), the uncertainties in the LVL and HVA models may affect each other in terms of the time and amplitude of the theoretical triplicated waveforms. We developed an iterative grid search method in which the models are updated in accordance with the CC between the observed and synthetic waveforms15. In the first iteration, we search for the optimal parameters for the LVL with respect to the IASP91 model. In the second iteration, the obtained LVL is included to further calculate the HVA. In the following iterations, we search for the updated optimal parameters for both the LVL and HVA based on the earlier iterations until the maximum CC becomes stable. Given the expense of numerical simulations, this kind of iterative grid search can eliminate much of the work for waveform modelling. In waveform modelling, the density perturbations are proportional to the velocity variations.

The optimal structural parameters for the LVL are sought in the range of dup_LVL = 300–400 km, dlow_LVL = 390–430 km and dVp_LVL = −0.5 to −6.0%, with intervals of 5 km and 0.5% based on global variations in this structure as obtained in previous studies5,6,7,8,9,10,11,12,13,14,15. Similarly, the parameter ranges for the HVA are set to dup_HVA = 440–620 km, dlow_HVA = 650–700 km, and dVp_HVA = 0.5–2.5%14,15.