A Displaced Lower Mantle Source of the Hainan Plume in South China Revealed by Receiver Function Imaging of the CEArray

We analyzed 49,592 teleseismic receiver functions (RFs) recorded by 278 CEArray stations to image the mantle transition zone (MTZ) beneath the South China Block to understand the origins of deep velocity anomalies and their potential links to subduction and intraplate volcanism. We employed a fast‐marching method and a high‐resolution 3‐D velocity model (FWEA18) derived from full waveform inversion in computing P‐to‐S conversion times to better image the 410‐ and 660‐km discontinuities. Our results indicate that the common‐conversion‐point stacking of RFs using 3‐D conversion times yielded better migration images of the two discontinuities. The images revealed a slightly depressed 410‐km with a few small uplifted patches, and showed that the 660‐km beneath the western Yangtze Craton is depressed by 10–25 km, which is likely caused by the stagnant Paleo‐Pacific slab. The 660‐km beneath the southern Cathaysia Block has a 5–15 km high plateau with a topographic low at its central part. The lateral dimension of the topographic low is ∼150 km and is located beneath the central Pearl River Mount Basin near Hong Kong. We speculate that the topographic low occurs within the Hainan plume with a temperature excess of ∼300–400 K and is caused by the garnet phase transition. The displaced deep plume enters the MTZ and spreads nearly horizontally at the base. The plume evolves into two channels with a minor one toward the northeast and a major one toward the southwest, which keep moving upward to the 410‐km. The southwest channel is likely the source that feeds the Hainan volcanoes.


Introduction
The South China Block (SCB) was formed during the Neoproterozoic when the Yangtze Craton (YC) in the northwest and the Cathaysia Block (CB) in the southwest amalgamated along the Jiangnan orogenic belt (e.g., Zhao, 2015).SCB is situated on the southeastern part of the Eurasian plate and has experienced a long history of deformation and magmatism associated with subduction of oceanic plates, collision of continental plates and upwelling of mantle plume (Briais et al., 1993;Y. Zheng & Zhang, 2007).Tomographic investigations have unveiled significant high-and low-velocity anomalies inside the mantle transition zone (MTZ) under the SCB (e.g., J. Huang & Zhao, 2006;C. Li & van der Hilst, 2010;L. Zhou et al., 2012).While the high-velocity anomalies are usually referred to the subducted slabs of oceanic plates or delaminated continental lithospheres, the low-velocity anomalies have been attributed to the Hainan plume arising from the lower mantle (e.g., J. Huang, 2014;Tao et al., 2018;Xia et al., 2016).The Hainan plume was first proposed by Lebedev et al. (2000) and Lebedev and Nolet (2003) based on their S-wave tomography model.However, the exact location of the Hainan plume in the lower mantle is not well resolved partially because traveltime tomography intrinsically has poor resolution on deep low-velocity structures.In addition, the sensitivity of seismic velocity to temperature becomes weaker in the lower mantle (e.g., Trampert et al., 2001), making it harder to image temperature anomalies in the lower mantle.
The MTZ, situated between the upper and lower mantle, is characterized by sharp increases in seismic velocity and density in the top and bottom of the zone, which are known as the 410-and 660-km seismic discontinuities (hereafter referred to as 410-and 660-km).The two seismic discontinuities occurred at depths where mantle minerals undergo phase changes in crystal structure observed by experimental studies (e.g., Katsura & Ito, 1989).The 410-km is associated with a phase transition of olivine to wadsleyite at a pressure of ∼14 GPa, while the 660km corresponds to the transition from ringwoodite to perovskite (also known as bridgmanite) and ferropericlase at a pressure of ∼24 GPa (e.g., Ita & Stixrude, 1992).Experimental studies (e.g., Ita & Stixrude, 1992;Katsura & Ito, 1989) showed that the 410-km phase transition is an exothermic reaction with a positive Clapeyron slope and the 660-km one is an endothermic reaction with a negative Clapeyron slope; a decrease in temperature (i.e., cold slab) results in a shallower 410-km and a deeper 660-km, leading to a thicker MTZ, while an excess of mantle temperature (i.e., hot plume) results in a thinner MTZ owing to a depressed 410-km and an uplifted 660-km.In summary, the depths of the 410-and 660-km in different tectonic settings are expected to vary with mantle temperature, and thus their depths can be used as thermometers to make in situ measurements of temperature inside the MTZ.Therefore, mapping lateral variations of the 410-and 660-km can provide information complementary to velocity structures in constraining temperature anomalies around the MTZ.
Receiver functions (RFs) have been widely used in imaging the 410-and 660-km beneath a seismic array.In the SCB, many studies using RF data recorded by arrays with various sizes and densities have been conducted (e.g., Ai et al., 2007;Han et al., 2020;He & Santosh, 2021;H. Huang et al., 2015;R. Huang et al., 2014;Wei & Chen, 2016) to investigate the MTZ beneath different parts of the SCB.These studies have found significant lateral variations in the depth of the two discontinuities and some correlations between the observed depth structures and velocity anomalies in tomographic models.However, we noticed that several of these studies used 1-D reference models (e.g., Ai et al., 2007) in computing the traveltime moveouts of the P-to-S conversions at the two discontinuities; therefore, the estimated depths of the 410-and 660-km could be significantly biased by unmodeled 3-D velocity structures.For example, if the upper mantle has a vertically extended anomaly with its Pand S-wave velocity being 2% and 3% higher than that of the reference model, then the imaged depths of the 410and 660-km could be, respectively, ∼15 and ∼24 km shallower than their actual depths.Due to the lack of P-and S-wave velocity models with comparable accuracy and resolution, the other studies (e.g., H. Huang et al., 2015) used either a 3-D P-or S-wave velocity model or a constant scaling relationship of δlnVs/δlnVp to compute the Pto-S conversion times, which could also introduce some artifacts to RF images.Therefore, there are large discrepancies in the estimated depths of the two discontinuities as well as the MTZ thickness among these studies, which could be caused by differences in data coverage and reference velocity models.Consequently, questions such as where the Hainan plume emerges from the lower mantle and what trajectory it takes in the MTZ are yet to be answered by RF images.
In this study, we utilized a high-resolution 3-D P-and S-wave velocity model, FWEA18 (Tao et al., 2018), which was derived from full waveform inversion of both body wave and surface wave data, as the reference velocity model and the fast marching method (FMM) in numerically solving the 3-D eikonal equation (de Kool et al., 2006;Guan & Niu, 2018;Rawlinson & Sambridge, 2004a, 2004b) to compute accurate relative traveltimes of P-to-S converted waves with respect to the direct P arrivals.We compiled 49,592 high-quality RFs recorded by 278 broadband seismic stations that provide a good coverage of the SCB (Figure 1b).We also employed a common-conversion-point (CCP) stacking technique to perform the depth migration of the RFs.The P-to-S conversion events at the 410-and 660-km are clearly shown in the CCP images.We compared lateral variations of the discontinuity depths with velocity anomalies of the FWEA18 model, which provide valuable insights into the deep geodynamic processes in the SCB region, especially regarding the deep origin of the Hainan plume as well as its kinematics in the MTZ and upper mantle.

CEArray Data
The waveform data used in this study are recorded by the permanent broadband seismic network operated by the China Earthquake Administration (CEA), which is hereafter referred to as the CEArray.The CEArray consists of a backbone national seismograph network (CNDSN), 31 regional networks, and several small aperture arrays with more than 1,000 stations (X.Zheng et al., 2009).We selected 261 CEArray stations located in a rectangular area between 100°and 122°E and 18°-30°N, which includes the SCB and its surrounding areas (Figure 1b).We also included 17 temporary stations in the southeastern margin of the Tibetan Plateau (SE Tibet) to ensure a relatively even distribution of the broadband stations.We examined waveforms from 623 earthquakes, each with a magnitude greater than 5.0 and an epicentral distance ranging from 30°to 90°(Figure 1a inset), which were recorded by the CEArray stations between 08/2007 and 07/ 2010 and the temporary stations between 09/2003 and 09/2004.These earthquakes offer good coverage in terms of both distance and azimuth (Figures 2a and 2b).

RFs
We followed Z. Liu et al. (2015) to compute RFs from teleseismic data and performed CCP stacking to image the 410-and 660-km.We first used the method of Niu and Li (2011) to estimate true orientations of the two horizontal components for each station using the collected teleseismic data, and then applied a correction if a station showed a misalignment of greater than 5°.The north (N) and east (E) components were then rotated to the radial (R) and transverse (T) directions.We applied a bandpass filter to the R-T-Z with a corner frequency of 0.01-5 Hz and down-sampled the data to 20 samples per second (sps).We applied the "water-level" deconvolution technique (Clayton & Wiggins, 1976) to the R and Z records to generate RFs.The water level was set to 0.01 and the Gaussian lowpass filter parameter was set to 1.5, equivalent to a corner frequency of ∼0.5 Hz.We manually examined all the RFs and selected a total of 49,592 RFs with a high signal-to-noise ratio (SNR) for CCP stacking.The P-to-S converted waves at the 410-km (P410s) and the 660-km (P660s) are clearly observed around the expected arrival times (Figure 2c), indicating sharp velocity contrasts at the two depths.

3-D Pds Traveltimes
We applied CCP stacking based migration technique to enhance P-to-S conversions that share a common geographic location.Here we first computed the raypaths of a P-to-S conversion at a depth, d, below the seismic station, Pds, by raytracing the 1-D IASP91 model (Kennett & Engdahl, 1991).Figures 3a and 3b show the geographic distributions of the P410s and P660s, which suggest that the study area is well sampled by the RF data set.Seismic migration is the process of back projecting signals recorded on seismograms in the time domain to their subsurface source structures in the depth domain.To perform the time-to-depth conversion, the relative traveltime of Pds with respect to the direct P arrival needs to be computed through a reference model.As the Pds traveltime is affected by conversion depth and Vp and Vs structure above the conversion, an accurate 3-D reference model is crucial for the depth migration.Conventionally, the relative Pds traveltime is first computed through a 1-D average model, such as the IASP91 model, and then a 3-D traveltime correction is computed by integrating anomalies along the 1-D raypath of Pds using the 3-D reference model (Z.Liu et al., 2015).This raytracing approach is very time-consuming when tens of thousands of RFs are used.In addition, the 1-D raypath can also introduce significant errors in Pds traveltime when prominent 3-D velocity anomalies are present in the upper mantle.We adopted the technique developed by Guan and Niu (2018) to calculate the 3-D Pds traveltimes with the 3-D FWEA18 model using the numerical eikonal solver, FMM (de Kool et al., 2006).Here we briefly review the three steps in computing Pds traveltime for each RF: (a) for each earthquake, we define a 3-D volume (hereafter referred to as EP-volume) that covers the study areas plus a margin and a depth range of 0-800 km.The margin width is determined by the epicentral distance of the earthquake, which varies from 8°to 3°in the distance range of 30°-90°.We employed the IASP91 model to compute the traveltimes of the grids at the boundaries of the volume and applied the FMM to the FWEA18 model to compute the 3-D P-wave traveltimes within the volume, TP(l,m,n), which were stored as the traveltime table for interpolation.(b) For each station, we applied the FMM to FWEA18 to calculate a 3-D S-wave traveltime table (T S (i,j,k)) for another 3-D volume (RS-volume), which is centered at the station with a lateral dimension of 6°× 6°and extends from surface to 800 km deep.When we computed the two traveltime tables, we employed a cell size of 0.05°and 5 km in the horizontal and vertical directions, respectively.(c) For every RF, we computed the absolute traveltime of Pds (t Pds ) and the direct P (t P 0 ), as well as the Pds relative traveltime with respect to P (Δt Pds = t Pds − t P 0 ) through a linear interpolation of the two 3-D traveltime tables (Y.Zhang et al., 2022).We employed the 1-D Pds conversion point location of the IASP91 model as an approximation in computing Δt Pds .To confirm this assumption, we randomly selected 2000 RFs to track their 3-D Pds conversion points and to calculate the associated Pds traveltimes, we found that the 3-D Pds conversion traveltimes are very similar to those computed with 1-D Pds conversion locations.More than 90% of the P660s show a difference less than 0.1 s in traveltime.The location discrepancy is smaller than the Fresnel zone for a P660s wave with a dominant period of 1.5 s, and the traveltime difference is comparable to the precision of the FMM.Thus, our assumption on the Pds raypath should not cause systematic biases on the estimated depths of the two discontinuities.

CCP Stacking
To perform CCP stacking, we first digitized the rectangular study area of 100°-124°E and 18°-30°N with a set of 0.1°× 0.1°grids, resulting in a total of 29,161 (241 × 121) grids.For each grid, we set up a circular bin centered at the grid with a radius of 1°and gathered all the RFs within the bin.Hereafter we refer to the number of RFs in a bin as the hit count of the bin.Figures 3c and 3d show the hit counts at depths of 410 and 660 km, respectively.If a bin has a hit count ≥150, we stacked all the RFs within the bin using a 0.2 s long time window centered at Pds traveltime (Δt Pds ).We varied the conversion depth, d, from 300 to 800 km with an increment of 1 km.It should be noted here that the 3-D image volume (100°-124°E, 18°-30°N, 300-800 km) is different from the EPvolumes and RS-volumes used in computing the P-and S-wave traveltime tables, respectively.We used two sets of Pds traveltimes (Δt IASP91 Pds and Δt FWEA18 Pds ), one computed from the 1-D IASP91 model and another from the 3-D FWEA18 model, to migrate the RF data.

Calibration of FMM With Raytracing the IASP91 Model
To verify the accuracy of the FMM in computing Pds traveltimes, we employed the 1-D IASP91 model to conduct calibration tests using the raytracing method.For each RF, we calculated the 1-D traveltimes of P410s and P660s using the raytracing (Δt IASP91 RT P410s , Δt IASP91 RT

P660s
) methods, and further computed their absolute differences (δt The P410s traveltime differences (δt IASP91 P410s ) of all the 49,592 RFs lie within 0.1 s, with an average and standard deviation of 0.027 and 0.019 s (Figure S1a in Supporting Information S1).Similarly, all the P660s traveltime differences (δt IASP91 P660s ) fall within a 0.2 s range, with an average and standard deviation of 0.067 and 0.035 s (Figure S1b in Supporting Information S1).We also confirmed the similarity of all the Pds times in the depth range of 300-800 km and concluded that the FMM provides an accurate way in computing Pds traveltimes.

CCP Stacking With the 1-D IASP91 Model
We organized the 3-D image volumes into depth profiles at each of the 0.1°× 0.1°grids.Due to the large number of RFs used in the CCP stacking, almost all the profiles showed clear P410s and P660s.We manually picked the depths of the two conversion signals from the depth profiles.More specifically we tried to pick the P410s and P660s within the depth ranges of 380-440 km and 640-690 km, respectively.Figures 4a and 4b show the histograms of the depths to the 410-and 660-km.The 410-km varies from ∼400 to 440 km, with an average depth of 418.8 km and a standard deviation of 8.30 km, while the 660-km is located in the depth range of ~650-692 km, with an average of 672.0 km and a standard deviation of 7.87 km.Figures 5a and 5b display the lateral variations of depths to the 410-and 660-km derived from Pds traveltimes calculated with the 1-D IASP91 model, and Figure 5c shows the lateral variations of MTZ thickness, which is the depth difference of the two discontinuities.The 410-km exhibits an apparent depth of ≥420 km in a large area covering the entire CB and the SE Tibet but is Geochemistry, Geophysics, Geosystems 10.1029/2023GC011292  shallower than ∼400 km beneath the YC (Figure 5a).The southeastward dipping trend in the apparent depth of the 410-km is also seen in the apparent depth to the 660-km discontinuity (Figure 5b).The apparent depth of the 660km beneath the southeastern part of the study area is greater than 670 km but is less than 660 km beneath the YC.The similar trend in the depth of the two discontinuities may also suggest that the upper mantle velocity beneath the CB is lower than that beneath the YC.The apparent MTZ thickness is relatively thin beneath the CB and the southeastern edge of the YC.The MTZ beneath the Hainan volcanoes and the Qiongzhou Strait is ∼240 km thick, ∼10 km thinner than the global average.The MTZ beneath the SE Tibet and western YC is, however, ∼10 km thicker than the global average.

CCP Stacking With the 3-D FWEA18 Model
In the time-to-depth conversion of Pds, there is a tradeoff between the conversion depth and the upper mantle velocity.Thus, the measured depths of the two discontinuities can be biased by unmodeled 3-D velocity structures in the upper mantle and MTZ that are absent in the 1-D IASP91 model.This is clearly shown in the Pds traveltimes calculated with the IASP91 (Δt IASP91

Pds
) models.Figures S1c and S1d in Supporting Information S1 show the absolute traveltime differences of P410s (δt and P660s (δt , respectively.δt P410s has a maximum of ∼3 s with an average and standard deviation of 0.850 and 0.487 s, respectively.Meanwhile, δt P660s shows a maximum of ∼4 s with an average and standard deviation of 1.278 and 0.741 s. Figures 4c and 4d show the histograms of the depths to the 410-and 660-km measured using the FWEA18 model.The 410-km depth varies from 393 to 427 km, with an average of 413.9 km and a standard deviation of 4.25 km (Figure 4c), and the 660-km depth ranges from 646 to 685 km, with an average of 661.1 km and a standard deviation of 7.00 km (Figure 4d).Both depths are significantly different from those measured using the IASP91 model (Figures 4a and 4b).In particular, we noticed that the standard deviation of the depth to the 410-km is almost half as 1-D estimates.
Figures 5d and 5e show mapviews of the depths to the 410-and 660-km obtained by the 3-D model, which are also significantly different from those estimated with the 1-D model.The maximum depth change of the 660-km is approximately 40 km, corresponding to the ∼4 s maximum difference in the P660s traveltime calculated using the two models.The 3-D results revealed that the 410-km is predominantly characterized by depressions ranging from 5 to 15 km with a few small uplifted patches (Figure 5d).Notably, there are two nearly parallel depression zones along the NW-SE direction on the topographic map of the 410-km (P1, P2 in Figure 5d).In contrast, the 660-km in the northwestern part of the study region is significantly deeper compared to the southeastern part (area C in Figure 5e).Additionally, the northern area of the YC exhibited a large-scale depression of the 660-km, with a magnitude ranging from 10 to 25 km.There is a broad region covering the southwestern half of the CB and its coastal area that shows an uplifted 660-km with an amplitude of 5-15 km (area A in Figure 5e).In the middle of the uplifted area A, the 660-km has a low topography with an almost normal depth, which is noted by B in Figure 5e.The depression seems to peak at the central part of the Pearl River Mouth Basin (PRMB) near Hong Kong (area B in Figure 5e).We also noticed that several high-amplitude anomalies are located at the edges of the study region, for example, E1 (Figure 5d) and E2 (Figure 5e) in the northeast corner and E3 and E4 (Figure 5e) at the southeast edge.Since they are relatively less well sampled by our RF data set, we decided not to interpret them in this study as the results require further verification.
To further show the lateral depth variations of the two discontinuities, we further showed seven depth profiles of the CCP stacking in Figure 6, including 3 E-W profiles along latitudes 22°N, 24°N, and 26°N, 3 N-S profiles along longitudes of 102°E, 110°E, and 114°E, and one NE-SW profile across the Hainan volcanic zone and PRMB.The locations of the 7 profiles are indicated in Figures 5a and 6h.In each profile, we also showed the Swave velocity perturbations of the FWEA18 model.All the profiles clearly showed P410s and P660s, allowing accurate picking of the depths of the two discontinuities.These profiles also clearly showed the major features about the two discontinuities revealed by  is clearly visible along the 114°E profile (Figure 6f) and the SW-NE profile (Figure 6g) and the uplift is surrounded by low-velocity anomalies.
Figure 5f shows the lateral distribution of the MTZ thickness estimated from the FWEA18 model.Compared to the depths of the two discontinuities, the MTZ maps estimated from IASP91 and FWEA18 models displayed more similarities, albeit with some disparities in degree and extent.The MTZ thickness ranges from 221 to 270 km, with an average value of approximately 247.3 km (Figure 5f).As indicated in Figure 4c, the 410-km exhibits a smaller lateral variation than the 660-km.Therefore, the MTZ thickness map is primarily determined by fluctuations in the depth of the 660-km.The MTZ is much thicker in the western and northern parts of the study area and becomes thinner in the central and southeastern parts of the study region.The variations in the MTZ thickness beneath the YC and CB align coherently with velocity anomalies in the MTZ of the FWEA18 model (Figure S2 in Supporting Information S1).A significant proportion of the high-velocity anomalies are situated in the northwestern part of the study area, which corresponds to the thickening of the MTZ.
Conversely, in the southeastern region, large-scale low-velocity anomalies are observed, correlating with a thinning of the MTZ.

Discussion
It is known that an accurate reference velocity model is critical in correctly positioning events in migration of seismic data.In the case of imaging MTZ with RFs, P-and S-wave velocity models with similar accuracy and resolution are required to obtain the reliable depths of the 410-and 660-km.Early RF studies of the SCB used 1-D reference models (e.g., Ai et al., 2007).H. Huang et al. (2015) employed a 3-D S-wave velocity model and a scaled P-wave velocity model by assuming a fixed δlnVs/δlnVp.We thus believe that there is a need to employ a high-resolution velocity model to better image the 410-and 660-km beneath the SCB.We further argue that FWEA18 fits the need as a reference model for RF imaging, since the P-and S-wave velocity models possess the similar model uncertainty and spatial resolution as Vp and Vs were inverted simultaneously from full waveform data including both body and surface waves.
To verify that the large discrepancy of the estimated depths of the two discontinuities is attributed to the reference models, we used the P-and S-wave velocity as well as the migrated depth of the 660-km of the 3-D FWEA18 model to generate synthetic waveforms of P660s using the similar source-receiver geometry of the data.We then employed the IASP91 and FWEA18 as the reference models to migrate the synthetic RFs.The migrated depths of the 660-km agree almost perfectly with the two sets of depths derived from the data.Figure S3 in Supporting Information S1 shows an example of the migrated depths of the 660-km of the synthetic (dotted lines) and data (solid lines) using the IASP91 (top) and the FWEA18 (bottom) models, respectively.
A highly accurate velocity model can lead to constructive stacking of the P410s and P660s waves, resulting in high amplitudes of the two phases.Thus, the amplitude ratio of the two phases between RF images migrated with different reference models is an index to evaluate the quality of a reference model.To do so, we first measured the amplitudes of P410s and P660s from CCP stacking profiles based on the IASP91 (A IASP91 P410s , A IASP91 P660s ) and FWEA18 (A FWEA18 P410s , A FWEA18 P660s ) models, and then computed their ratios (AR P660s / A IASP91 P660s ).Figures S4a and S4b in Supporting Information S1 show the histograms of the P410s and P660s amplitude ratios; Over 60% and 90% of the profiles showed an amplitude ratio >1 and >0.9, respectively.A map view of the two amplitude ratios (Figures S4c and S4d in Supporting Information S1) also suggested that most of the study area showed a similar or an increased amplitude of P410s and P660s with the FWEA18 model.For these reasons, our discussion and interpretation of the MTZ structure in the following sections will be based on the CCP images derived from the FWEA18 model.

Reconcile Depth Variations of the Two Discontinuities and Velocity Anomalies to MTZ Thermal Structures
To better constrain the origin of the topography of the 410-and 660-km, we conducted a quantitative comparison with the P-and S-wave velocity anomalies in the MTZ shown in the FWEA18 model.More specifically, we would like to investigate whether temperature anomalies are the primary cause of the observed depth and velocity anomalies.Although SS precursor studies (Flanagan & Shearer, 1998;Gu & Dziewonski, 2002) suggested that the global averaged depths of the two discontinuities are 418 and 660 km, respectively, we still used 410 and 660 km as the reference depths to compute the depth anomalies of the two discontinuities.The average depths of the 410-and 660-km are 413.9 and 661.1 km, resulting in an average MTZ thickness of 247.3 km beneath the study region.
As shown in the previous section, the MTZ exhibits different structures beneath the YC and the CB.Thus, we chose two rectangular areas that have high hit counts of RF data (100°-110°E, 24°-30°N) from the YC and (110°-118°E, 20°-26°N) from the CB, for more detailed comparison.The two areas are outlined in Figure 3.The average depths of the 410-km beneath the two areas are identical: 414.8 km for the YC and 414.8 km for the CB.These values closely resemble the overall average depth of the 410-km in the whole SCB, indicating a widespread depression throughout the study area.The average 660-km depth beneath the YC is 666.4 km and 654.1 km beneath the CB, resulting in a difference of approximately 12 km.
If the observed depth anomalies are caused by temperature-sensitive phase transitions of the upper mantle minerals, and if there are thermal anomalies extending the entire MTZ, then we would expect to observe correlations/anticorrelations among the depths of the two discontinuities as well as the MTZ thickness.If phase transitions of olivine are the major causes of the two seismic discontinuities, then there is an anticorrelation between the depths of the 410-and 660-km.In Figure 7, we plot the deviations of MTZ thickness versus the 410km depth, MTZ thickness versus the 660-km depth, and 410-km depth versus the 660-km depth beneath the whole region, the YC, and the CB.The 410-km depth beneath the CB (Figure 7c) displays a stronger negative correlation with the MTZ thickness than that beneath the YC (Figure 7b).Meanwhile, there is a strong positive correlation between the MTZ thickness and the 660-km depth across the entire study area (Figure 7d), including the YC (Figure 7e) and CB (Figure 7f).The correlation between the depths of the 410-and the 660-km is relatively low (Figures 7g-7i), with the highest of 0.294 from the CB (Figure 7i).Based on these observations, we speculated that the thicker MTZ beneath the YC is primarily caused by a depressed 660-km, which is likely originated from low temperature.This thermal anomaly also does not seem to extend throughout the entire MTZ.In comparison, the high temperature anomalies within the MTZ below the CB are widespread across the transition zone, affecting both the 410-and 660-km topographies.
Using the above scaling relationship, we converted the maximum 15 km depression of the 410-km beneath western YC (area C in Figure 5d) to an elevated mantle temperature of 182-227 K. Similarly, the deepened 660km by 20 km observed beneath the same region (area C in Figure 5e) corresponds to a cold surrounding mantle with a relative temperature of − 270 to − 394 K. On the other hand, the temperature in the P1 and P2 areas (Figure 5d) of the CB is uplifted by ∼5 km, corresponding to a mantle warmer by 61-76 K in temperature.The broad uplift of 660-km beneath the southern CB with a peak amplitude of 15 km (area A in Figure 5e) can be resulted from a hotter-than-normal mantle with a temperature anomaly of 203-296 K.We further converted these temperature anomalies to P-and S-wave velocity perturbations based on their temperature sensitivities ((∂lnV P / ∂T) 410 , (∂lnV S / ∂T) 410 ; (∂lnV P / ∂T) 660 , (∂lnV S /∂T) 660 ) estimated by Cammarano et al. (2003).The calculated velocity perturbations are listed in Table S1 of Supporting Information S1.These values are roughly comparable to those of the FWEA18 model (Figure S2 in Supporting Information S1).However, the temperatureconverted S-wave perturbations at 410 km beneath the YC and at 660 km beneath the CB are much larger than those in the FWEA18 model.We further computed the correlation between depth variations of the two discontinuities and S-wave velocity perturbations (dlnVs) of FWEA18 at each grid (Figure S5 in Supporting Information S1).In general, the correlation between the two is not so obvious (Figure S5 in Supporting Information S1), which may suggest that temperature anomalies might not be the sole cause of the observed velocity and depth perturbations.The low correlation could also be caused by errors in the depth and velocity estimates.In particular, the area northwest of Taiwan (hatched area in Figure S2b of Supporting Information S1) shows a strong low-velocity anomaly, implying that there is an elevated temperature.In contrast, the region is featured by a depressed 660-km discontinuity (Figure 5e), which usually means a delayed post-spinel phase transition under a low temperature condition.We noticed that the P660s amplitude migrated with the FWEA18 model in this area is lower than that of the IASP91; therefore, the migrated depths of the 660-km could be incorrect in this area, which might be the cause of the above inconsistency.Another area is north of Hainan Island, outlined by the black dashed line (Figures S2b and S2d in Supporting Information S1), where a high S-wave velocity and shallow 660km are observed.We confirmed that the CCP stacked amplitude of P660s based on the FWEA18 model is higher than that of the IASP91 model, indicating that the depth estimate is likely robust.Moreover, the correction of the mentioned high velocity is expected to lead to a deeper 660-km; thus, the elevated 660-km in this area is not the artifact of the high-velocity anomaly.We notice that the high-velocity anomaly in Figure S2b of Supporting Information S1 corresponds to an area with less elevated temperature as compared to its surroundings (Figure S2d in Supporting Information S1).

Slab Segments Within the MTZ
The most prominent area with a depressed 660-km is in the YC region (region C in Figure 5e), except for the northeastern corner and southeastern edge of the study area (E2, E3, and E4 in Figure 5e), which need further verification with more data.At the C area, the maximum depression of the 660-km is ∼25 km compared to the global average (660 km).The most plausible explanation for this depression is the presence of low temperature structures, which are imaged as high velocity bodies around the 660-km by many tomographic studies (e.g., Sun et al., 2016;Tao et al., 2018;R. Zhang et al., 2017).Various interpretations have been offered regarding the possible origin of the high-velocity anomalies, such as detached lithosphere (e.g., R. Zhang et al., 2017), the eastward subducting Indo-Burma slab (e.g., Xu et al., 2018), and the stagnant slab of the Paleo-Pacific plate (e.g., Sun et al., 2016).Our CCP stacking images showed that the 660-km in the C area has a maximum depression of 25 km, suggesting that the mantle at the base of the MTZ in the C area is ∼337-492 K colder than its adjacent areas.This temperature difference could result in 5%-7% perturbations in S-wave velocity, which is consistent with the FWEA18 model (Figure S2 in Supporting Information S1).
Our CCP stacking images, however, do not offer additional constraints on the origin of the high-velocity anomalies, as the depth of the 660-km is affected by not only temperature but also water content (e.g., Higo et al., 2001;Litasov et al., 2005).The 3-D perspective view of high-velocity anomalies of the FWEA18 model beneath the study area (Figure 8a) indicates that the high-velocity anomalies across the 660-km in the C area are isolated structure without obvious connection with the subducting high-velocity slabs of the Philippine Sea plate from the east and the India plate from the west.Since the northern part of the depression is located beneath the southern Sichuan Basin where the cratonic lithosphere is expected to be intact as a high-velocity anomaly down to ∼200 km deep is visible from many tomographic models (e.g., Tao et al., 2018), geometrically it is difficult to attribute the high velocity structures to the delaminated continental lithosphere.The depression is elongated along the NE-SW direction, which may imply a piece of slab being subducted from either the southeast or the northwest.Therefore, we favor more on an origin related to the Paleo-Pacific subduction.

Upwelling of the Hainan Plume
Compared to the YC, the CB region has experienced more active magmatic activity (e.g., X. M. Lhou & Li, 2000).Several studies suggest that the magmatic activity in this region may be associated with either lateral mantle flow resulting from plate collision (e.g., M. Liu et al., 2004) or mantle wedge circulation related to plate subduction (e.g., Maruyama et al., 2009).Lebedev et al. (2000) employed tomography to identify a low velocity zone surrounding the Hainan Island, which was subsequently designated as the Hainan plume.The mantle plume hypothesis, which involves the transport of hot materials from the deep mantle, has been effectively applied to explain intraplate magmatism in various locations with hotspots.Nevertheless, the precise location where the Hainan plume enters the MTZ, along with its connection with the Hainan volcanoes in the upper mantle, remains unknown.
Both global and local seismic traveltime tomography have been utilized to investigate the deep structure beneath this region.Some findings indicate that mantle upwelling may be rooted within the MTZ, while others propose that it may originate from the core-mantle boundary (e.g., He & Santosh, 2021; J. Huang & Zhao, 2006).J. Regarding where the Hainan plume enters the MTZ from the lower mantle, studies utilizing RFs to image the thinning position of the MTZ have produced multiple results.Some researchers proposed that the intersection of the Hainan plume and the MTZ occurs beneath the Qiongzhou Strait (H.Huang et al., 2015), while others suggested it is to the northeast of the Hainan Island, somewhere at the coastal area of the Guangdong Province (Wei & Chen, 2016).
In this study, we found a significant uplift of the 660-km beneath the southern CB (Figure 5e), which also displays notable thinning of the MTZ thickness (Figure 5f).The most prominent uplift area sits in the coastal region to northeast of the Hainan volcanoes (A in Figure 5e).One possible explanation of the uplift is the presence of a mantle plume in the area.A close look of the broad uplift also revealed its annular shape with a low area that is centered at 114.2°E and 22.6°N (near Hong Kong in Figure 1b) with a diameter of approximately 150 km.The 660-km in the area displays a depression relative to the surrounding uplifted areas, which is labeled as B in Figure 5e.The depth profile shown in Figure 6g extending from the Hainan Island reveals a 660-km with an undulation of uplift-depression-uplift from SW to NE. Figure 8b shows 3-D distribution of low-velocity anomalies of the FWEA18 model beneath the study area.We found that the B area is at the center of a large low-velocity anomaly in the MTZ, which seems to extend to the lower mantle, suggesting a scenario that hot material arising from the lower mantle enters the transition zone at the B area and then spreads horizontally within the MTZ.The hot materials continue to rise and evolve into two branches corresponding to the two nearly parallel depression zones of the 410-km (P1, P2 in Figure 5d).The major branch keeps moving upward to the southwest and feeds the Hainan volcanoes (yellow arrows in Figure 6g).
In general, the depth of the 660-km is mainly controlled by the phase transformation of olivine mineral, the major constituent of the upper mantle.However, in a substantially high temperature condition, phase transitions of the non-olivine minerals, such as the majorite-garnet to perovskite transition, can play a major role in determining seismic properties (Hirose, 2002).The majorite-garnet to perovskite transition has a positive Clapeyron slope, meaning that a deeper, rather than a shallower 660-km is expected when mantle temperature is significantly higher than the global mean.This indicates that in the presence of substantial high temperature anomalies, the garnet phase transformation in mantle minerals may play a more important role, resulting in a deeper, rather than shallower, and seismologically observed 660-km.In fact, several previous studies have reported a deepened 660km within hot mantle beneath Kenya, Tanzania, and Hawaii that may be associated with mantle plumes (Cao et al., 2011;Deuss, 2007;Huerta et al., 2009).Thus, we postulate that the B area has a mantle temperature higher than its surroundings and the depression in the area is caused by the positive Clapeyron slope of the majoritegarnet to perovskite transition.
To investigate this hypothesis, we generated 2-D synthetic models of mantle with a pyrolite composition (Ringwood, 1962) and laterally varying geotherms.Since the garnet-perovskite phase change only occurs at elevated temperatures above 200-300 K (Jenkins et al., 2016), we designed two 2-D models based on the geometry of the NE-SW profile shown in Figure 6g: the first model with a lateral potential temperature (Tp) variation of 1600-1800-1900-1800-1600 K (Figure 9a; maximum Tp = 1900 K) and the second model with a varying Tp of 1600-1800-2000-1800-1600 K (Figure 9b; maximum Tp = 2000 K).The lateral extensions of the 1800 and 1900 K/2000 K segments are 200 and 150 km, respectively, which are roughly similar to the dimensions of A and B in the NE-SW profile (Figure 6g).We first employed Perple_X (Connolly, 2005) to compute seismic velocities and density of the two 2-D models.The calculated S-wave velocity of the two models are shown in Figures 9a and 9b, while the P-wave velocity and density of the two models are shown in Figure S6 of Supporting Information S1.We also showed depth profiles of Vp, Vs, and density at 5 grid points that represent the background mantle, moderately hot mantle and substantially hot mantle in Figure S7 of Supporting Information S1.
We then employed the 2-D finite-difference method (D.Li et al., 2014) to compute synthetic waveforms of the two models for all the event-receiver pairs used in the CCP stacking of the NE-SW profiles shown in Figure 6g.
The resulting synthetic seismograms were processed in a similar manner to the real data to obtain RFs, which were subsequently employed in the CCP stacking.The results of the CCP stacking of the synthetic RFs of the two models are shown in Figures 9c and 9d.The 660-km of the two CCP images exhibits undulations comparable to the real data image.In the two edges away from the hot center, the depth to the 660-km is close to normal (Figures 9c and 9d).Both images show that the 660-km rises and then falls when it approaches to the center with increasing temperature.We noticed that the depressions of the 660-km at the central area in the two synthetic images are much larger than that observed in the B area, suggesting either the Clapeyron slope used in the mineral physical models is overestimated or the ringwoodite to perovskite phase change of the olivine and water content might have played a role in the observed 660-km in the B area.
Our results thus may suggest that the Hainan volcanoes are likely fed by a displaced hot plume in the lower mantle.The CCP images indicate the following scenario: a conduit of hot plume materials with a diameter of ∼150 km rises from the lower mantle and hits the 660-km at central PRMB near Hong Kong, which shows a lowvelocity anomaly in the uppermost lower mantle and a depression of 660-km in the B area.Temperature within the plume is ∼300-400 K higher than that of a normal mantle at the same depths.Hot plume materials accumulate at the base of the MTZ and spread horizontally before progressing further toward the 410-km as two rising channels of hot materials (Figure 10).While the minor branch moves upward to the northeast, the major channel flows to the southwest, rising to the shallow depths and feeding the Hainan volcanoes at the surface.

Conclusions
We conducted CCP migration of a large RF data set with a 3-D high-resolution velocity model to image the 410and 660-km.By comparing velocity structures around the two discontinuities and numerical analyses with synthetic models of mantle minerals and different geotherms, we reached the following conclusions: (a) The 660km beneath the western YC shows a 10-25 km depression and is surrounded by high-velocity anomalies, suggesting that the depression is likely caused by a cold stagnant slab with higher-than-normal seismic velocity.
The NE-SW elongated shape of the depression may imply that the stagnant slab is associated with the northwestward subducted Paleo-Pacific plate.The cold stagnant slab does not seem to reach the upper part of the MTZ and affect the 410-km.(b) The 660-km beneath the southern part of the CB is uplifted by 5-15 km in a broad area with a lateral extension of ∼500-600 km.There is a low topographic value at the central part of the plateau with a diameter of ∼150 km.We attributed the topographic low to the positive slope of the garnet-perovskite phase transition in the hot mantle with an excess temperature of 300-400 K.By integrating seismic observations and mineral physical modeling, we propose that the Hainan plume has a displaced origin in the lower mantle, which is located at the central PRMB near Hong Kong.Once the plume enters the MTZ, it spreads horizontally and diverges into two channels, one minor branch toward the northeast and one major branch toward the southeast.Both

Figure 1 .
Figure 1.(a) Topographic map showing the Southeast Asia and the South China Block.The white solid and dashed lines denote the boundaries of plates and active tectonic blocks, respectively.The rectangle outlined by the red dashed lines represents the study area that includes the Cathaysia Block, the Yangtze Block and part of the southeastern margin of the Tibetan plateau.The inset in the upper left corner is an azimuthal equidistant map showing the 623 earthquakes (red stars) used in computing receiver functions.The blue triangle in the center indicates the center of the study area, and three circles denote the epicentral distances of 30°, 60°, and 90°, respectively.(b) Distribution of the 278 stations used in the study (blue triangles).The red stars denote the Hainan volcanoes.PRMB: Pearl River Mouth Basin.

Figure 2 .
Figure 2. (a) Histograms showing the distribution of receiver functions (RFs) along epicentral distance and (b) back azimuth.(c) Stacked RFs in each 1°-bin plotted against the epicentral distances.Red and blue colors denote the positive and negative pulses, respectively.Dashed lines denote the theoretical arrival times of P410s and P660s predicted by the IASP91 model.

Figure 3 .
Figure 3. (a) Lateral distribution of the P410s and (b) P660s piercing points are shown in red dots.(c) Distribution of the P410s hit count and (d) P660s hit count.The black rectangles indicate the representative areas of the Yangtze Craton and Cathaysia Block, respectively, for statistical analyses of the depths of the 410-and 660-km.

Figure 4 .
Figure 4. (a) Histograms of the depth to the 410-km and (b) the 660-km estimated from each circular bin with a hit count ≥150 using the 1-D IASP91 model as the reference model based on fast marching method.Std: standard deviation.Panels (c) and (d) are the same as (a) and (b) except for the reference model, which is the 3-D FWEA18 model.

Figure 5 .
Figure 5. (a) Color contour maps showing lateral variations in the depth to the 410-km and (b) to the 660-km as well as in the thickness of the mantle transition zone (c).The solid purple lines show the geographic locations of the seven profiles in Figure 6.The color scales are shown at the bottom of each panel.Depths and thickness are estimated with a time-to-depth conversion computed from the 1-D IASP91 model.Panels (d)-(f) are the same as (a)-(c) except for the reference model, which is the 3-D FWEA18 model, used in computing the time-to-depth conversion.P1, P2, A, B, and C denote anomalies mentioned in the main text.E1-E4 are the anomalies located at the edge of the study area and are not discussed in the main text.YC, Yangtze Craton; CB, Cathaysia Block; PRMB, Pearl River Mouth Basin.
Figure 5: (a) most of the profiles exhibit a slightly depressed 410-km; (b) the 660-km at the western part of the 26°N profile and the northern part of the 102°E profile (C in Figures6c and 6d) located beneath the western YC are depressed by a maximum of ∼25 km, and the area with a depressed 660km shows a high S-wave velocity; (c) the annular shaped uplifting of 660-km beneath the southern part of the CB Figure 6.

Figure 6 .
Figure 6.Panels (a-c) show three E-W depth profiles of the CCP stackings (solid black lines) along latitudes 22°N, 24°N, and 26°N.S-wave velocity perturbations of the FWEA18 model are also shown on the profiles for comparison.Red and blue colors denote the positive and negative anomalies, respectively.The two horizontal white dashed lines denote the reference depths of 410 and 660 km.Amplitude variances (Amp_Var: Amp bin /Amp mean ) of P410s and P660s along the profiles are shown above the profile with blue dashed and red solid lines, respectively.Panels (d-f) are the N-S depth profiles along longitudes 102°E, 110°E, and 114°E.Panel (g) is a depth profile along the NE-SW direction with an origin aligning with the volcano in Hainan Island (red star on the top left of the profile).Yellow arrows show the plausible flow direction of the southwest branch.The areas marked with A, B, and C correspond to those in Figure 5. (h) Topographic map showing the geographic locations of the seven profiles (solid purple lines) and the seismic stations used in the study.YC, Yangtze Craton; CB, Cathaysia Block; PRMB, Pearl River Mouth Basin.

Figure 7 .
Figure 7. (a) The mantle transition zone thickness perturbations are shown as a function of the depth variations of the 410-km measured in the entire study region, (b) the Yangtze Craton, and (c) the Cathaysia Block.The dotted line with a negative slope represents a perfect anticorrelation, and CC indicates the calculated cross-correlation coefficient between the two data sets.Panels (d-f) are the same as (a-c) except for the depth variations of the 660-km.The dotted line with a positive slope signifies an ideal correlation.Panels (g-i) are the depth variations of the 660-km are shown as a function of the depth variations of the 410-km measured in the three regions.
Huang (2014) observed a low-velocity anomaly beneath the Hainan volcanoes, which tilts toward the northeast and extends to 700 km deep.The shape of the low-velocity anomaly does not fall into a classic image of a mantle plume which ascends vertically from the deep mantle.Xia et al. (2016) proposed a mushroom-like continuous low-velocity anomaly, which is characterized by a narrow columnar tail extending from the lower mantle beneath the northeastern part of the Hainan hotspot and a broad head spreading laterally within and around the MTZ.

Figure 8 .
Figure 8. 3-D illustration of the (a) high (dlnVs ≥ 2.0%) and (b) low (dlnVs ≤ − 2.0%) velocity bodies of the FWEA18 in the depth range of 0-800 km beneath the South China Block region.The colored map shown at the bottom of each 3-D illustration represents the 660-km topography, which is the same as Figure 5e.

Figure 9 .
Figure 9. (a) Synthetic 2-D velocity profile showing S-wave velocity computed from a mineralogical model with a maximum Tp of 1900 K at the center.Tp decreases toward the edge with increments of 100 K and reaches 1600 K at the two edges of the profile.S-wave velocity is shown in color and the color scale is shown at the bottom of the profile.G1-G5 show the location of five grids with their 1-D velocity profile being shown in Figure S7 of Supporting Information S1.Panel (b) is the same as (a) except for a maximum Tp of 2000 K at the center of the profile.Synthetic tests for the 660-km undulations under extremely high temperatures.Panel (c) shows the CCP stacking of synthetic receiver functions computed with the 2-D velocity model shown in (a).Red color denotes the positive pulses and the light red rectangle represents the maximum Tp of 1900 K region.Panel (d) is the same as (c) except for a maximum Tp of 2000 K.

Figure 10 .
Figure 10.Schematic illustration showing the deep origin of the Hainan plume and its pathway in the mantle transition zone and upper mantle toward the Hainan volcanoes on the surface.CCP stackings along the NE-SW profile shown in Figure 6g are also plotted for comparison.