Imaging mantle transition zone discontinuities in southwest China from dense array ambient noise interferometry

1.   Introduction

Mantle transition zone (MTZ), bounded by the 410-km and 660-km discontinuities, includes sharp gradients of velocity and density generated by polymorphic phase transitions (Deuss, 2009; Shearer, 2000). MTZ plays a pivotal role in mantle dynamics and the seismic properties of the MTZ discontinuities provide an important constraint on the mantle composition and temperature together with mineral physics data (e.g. Hirose, 2002).

The MTZ discontinuities can be detected by triplications in traveltime and waveform modeling as well as some secondary seismic phases reflected or converted at the interface (Shearer, 2000). However, most of these methods are limited by the earthquake-station geometry. The body waves reflected at the MTZ discontinuities from ambient noise interferometry can provide high-resolution seismic constraints without the earthquake-station geometric limitation (Feng et al., 2017; Poli et al., 2012).

Surface waves extracted from ambient noise cross-correlation functions (NCFs) based on continuous seismic array data have been widely used to investigate the crustal and upper mantle structure in the past decade (Lin et al., 2008; Shapiro et al., 2005; Yang et al., 2007; Yao et al., 2006). Besides the dominant surface waves recovered from ambient noise, more and more body wave signals in various frequency bands have been identified on the NCFs although with much weaker amplitudes (Boue et al., 2013; Lin et al., 2013; Nishida, 2013; Poli et al., 2015; Zhan et al., 2010). And part of these body waves has been successfully utilized to investigate structure of the Earth’s interior, such as the MTZ discontinuities (Feng et al., 2017; Poli et al., 2012). Body waves reflected at the MTZ discontinuities can provide high-resolution constraints, comparable with receiver functions and higher than PP and SS precursor observations, on the MTZ structure without the earthquake-station geometric limitation.

Clear P₄₁₀P and P₆₆₀P phases can be observed on the NCFs section within the secondary microseism frequency band (0.1–0.2 Hz), but invisible within the primary microseism band (0.05–0.1 Hz) (Feng et al., 2017). The secondary microseism, which doubles the frequency of an ocean swell, is proposed to be generated through nonlinear wave-wave interaction (Longuethiggins, 1950). The pressure, which is exerted on the ocean bottom associated with nonlinear wave-wave interaction, can be estimated by single vertical forces that dominantly generate P waves (Nishida and Takagi, 2016). Clear teleseismic P waves can be observed within the secondary microseism band on the beamforming image (Wang et al., 2018a). The generated teleseismic P waves within the secondary microseism band will be a powerful tool to investigate the internal structure of the Earth.

In this study, we observed clear P₄₁₀P and P₆₆₀P signals recovered from the NCFs from a dense seismic array in southwest China. These reflected waveforms were first converted into the depth-domain signals based on a global 1-D layered velocity model. Then the depth-domain signals were stacked with the common middle point (CMP) stacking technique within different bins to reveal the lateral variations of the MTZ discontinuities.

2.   Raw waveform data and NCF calculation

In southwest China, the continuous waveform data of totally 438 available broadband seismic stations, including 350 portable stations from the Himalaya Project (ChinArray, 2006) and 88 permanent provincial stations from the China National Seismic Network (Data Management Center of China National Seismic Network, 2007; Zheng et al., 2010), were collected to calculate NCFs (Figure 1). Only the continuous records of the vertical component from March 2011 to November 2013 were utilized to calculate NCFs because vertical component data have higher signal-to-noise ratio than radial and transverse components.

Figure  1.  Distribution of seismic stations in southwest China, including portable stations (purple triangles) and permanent stations (blue triangles)

DownLoad: Full-Size Img PowerPoint

The standard procedure described by Bensen et al. (2007) was adopted to calculate the NCFs. The continuous records were first cut into daily segments followed by single-station data pre-processing and station pair cross-correlation. Single-station data processing includes the remove of the mean and trend of each daily segment, the deconvolution of the instrument responses, down-sampling to 1 Hz, time domain normalization and spectral whitening. Pre-processed daily segments were then correlated for different station pairs and all the daily NCFs for the same pair were linearly stacked as the final NCF to increase the signal-to-noise ratio and the final NCF was filtered to 0.1–0.2 Hz. Finally, 95,633 vertical-vertical component NCFs were obtained in total. All the final NCFs were arranged by distance as shown in Figure 2a. The number of NCFs within each 19-km width distance bin is shown in Figure 2b and most interstation distances of NCFs are between 60–800 km.

Figure  2.  (a) The stacked NCFs arranged by distance, (b) The number of NCFs within each 19-km width bin

DownLoad: Full-Size Img PowerPoint

Dominant Rayleigh waves and clear body wave signals at ~100 s and ~150 s can be identified on the NCF section (Figure 2a). Feng et al. (2017) have confirmed that the body waves emerging at ~100 s and ~150 s are P₄₁₀P and P₆₆₀P phases, respectively, by calculating the 4th-root vespagrams and comparing with the theoretical arrival times of P₄₁₀P and P₆₆₀P. As there are more NCFs in our study region, only the linear time-slowness slant stacking (Rost and Thomas, 2002) was performed and the stacking results of four distance bins (0–200 km, 100–300 km, 200–400 km, and 300–500 km) are shown in Figure 3. The local maxima of the slant stacking amplitude results have a good agreement with the theoretical predictions (the white circles) of P₄₁₀P and P₆₆₀P phases at corresponding distances. The slant stacking results have high signal-to-noise ratios and the retrieved body wave signals are not severely biased by the predominant Rayleigh waves. The early arrival body waves (0–50 s) are teleseismic P waves with high apparent velocity (~23 km/s) and detailed analysis of these early body waves can be found in Wang et al. (2018a, b).

Figure  3.  Envelope of the time-slowness slant stacking of NCFs within four distance bins: 0–200 km, 100–300 km, 200–400 km, and 300–500 km. The white circles indicate the theoretical arrival times and ray parameters of P₄₁₀P and P₆₆₀P phases at 100 km, 200 km, 300 km, and 400 km epicenter distances calculated with the Taup toolkit (Crotwell et al., 1999) from the ak135-f model (Kennett et al., 1995)

DownLoad: Full-Size Img PowerPoint

3.   Depth-domain conversion of NCFs and CMP stacking

Unlike the time-domain stacking in our previous study (Feng et al., 2017), all the NCFs were first converted to the depth-domain from the time-domain before they were stacked within each bin. And the depth-domain conversion is partly equivalent to the continuous time-domain traveltime correction (normal moveout correction, NMO). As the body wave signals, instead of the Rayleigh waves, are our targets, the whole NCF trace is regarded as the reflected waveforms from discontinuities at different depths. For a horizontally layered model, the relationship between the reflected traveltime t and the offset H (interstation distance here) can be expressed by a hyperbola (Sheriff and Geldart, 1995):

where t₀ is the traveltime at the zero offset, indicating that the event and seismic station are located in the same place; v is the velocity of the medium above the reflecting interface.

For a horizontally layered flat Earth model, there are n layers above the reflecting interface and the thickness and velocity of each layer are Z_i and v_i, respectively, with i=1, 2, …, n. The traveltime at the zero offset is

When the reflecting depth is larger than 2 times the offset, the velocity v in equation (1) can be approximately replaced by the root mean square velocity:

For a given offset, the ak135-f model was first interpolated to a layered model with all one-kilometer thick layers in the depth-domain and the traveltime of body wave reflected at each interface can be calculated with equation (1). The relationship between the reflected traveltime and the reflecting depth is shown in Figure 4c for four given epicenter distances. Then the corresponding time-domain NCF for a station pair was interpolated to obtain the depth-domain NCF according to the traveltime-depth curve for the given offset H (Figure 4c). Feng et al. (2017) have tested that the influence of a 3-D velocity model on the traveltimes of P₄₁₀P and P₆₆₀P phases is usually secondary and insignificant. Here all the NCFs were converted to the depth-domain NCFs simply based on the 1-D ak135-f Earth model.

Figure  4.  (a) Schematic plot of P₄₁₀P and P₆₆₀P ray paths. The black star and triangle denote the virtual source and receiver, respectively. CMB: core mantle boundary; ICB: inner core boundary. (b) Travel time curves of P₄₁₀P and P₆₆₀P calculated with the Taup toolkit (Crotwell et al., 1999) (red and green stars) and equation (1) (red and green lines) from the ak135-f model (Kennett et al., 1995), respectively. (c) The relationship between the reflecting depth and traveltime for given offsets: 50 km, 100 km, 150 km, and 200 km. The traveltime was calculated with equation (1)

DownLoad: Full-Size Img PowerPoint

As shown in Figure 2a, the interference of the dominant Rayleigh waves and P₄₁₀P phase occur within the interstation distance range of 200–400 km. In order to avoid the interference of Rayleigh waves and keep the relative amplitude information of P₄₁₀P and P₆₆₀P phases, only the NCFs with interstation distances less than 200 km were selected for final stacking. As our target depth range is 400–700 km (larger than 2 times the maximum offset), equation (1) and the root mean square velocity given by equation (3) can be utilized to calculate the traveltime for a given reflecting depth and offset.

As discussed above, the reflection time-distance equation can be expressed as a hyperbolic function for a horizontally layered velocity model. But the true Earth is spherical thus equation (1) is not accurate. In addition, the root mean square velocity is an approximation only when the reflecting depth is much larger than the offset. All these approximations will cause errors in the traveltime calculation and these errors need to be assessed before further analysis. To evaluate the accuracy of traveltime calculation, the traveltimes of P₄₁₀P and P₆₆₀P (noted as t_{_P410P} and t_{_P660P} respectively in Figure 4) calculated with equation (1) and the spherical Earth model ray-tracing from the ak135 1-D Earth model should be compared. Figure 4 shows the ray paths of P₄₁₀P and P₆₆₀P phases and their traveltimes at different offsets calculated with two different methods. The maximum discrepancy of t_{_P410P} is less than 0.2 s for the offset range 0–200 km and the difference is even smaller for t_{_P660P}. Converted to the depth-domain, 0.2 s traveltime error will cause about 1 km bias in reflecting interface depth, which is negligible in our study.

Following Feng et al. (2017), we stacked the depth-domain NCFs along the two profiles marked as the purple arrowed lines in Figure 5. Figure 5 shows an example bin as the black circle and the bin centers are indicated as the yellow stars. All the depth-domain traces whose reflecting points are located within the example bin (green dots in Figure 5) are shown in Figure 6a and the stacked traces are shown in Figure 6b. As the body waves reflected at shallower interfaces are severely interfered by the Rayleigh waves, the shallow signals were tapered off using the black taper shown in Figure 6b. The waveform distortion is insignificant when the depth is larger than two times the offset. Therefore, the phase weighted stacking (PWS) technique (Schimmel and Paulssen, 1997) can still be applied to the depth-domain NCFs. The linear stacking result and the PWS result have high similarity. However, the PWS result has a higher signal-to-noise ratio.

Figure  6.  (a) All the depth-domain NCFs whose middle reflecting points are located within the example bin shown as the black circle in Figure 5; (b) The linear stacking (green line) and PWS (red line) results of all the NCFs shown in (a) and the taper for extracting the main MTZ reflection phases used in the depth domain

DownLoad: Full-Size Img PowerPoint

Figure  5.  Seismic station distribution and the locations of two profiles (purple arrowed lines) for NCF stacking. The yellow stars indicate the stacking bin centers and the black circle denotes the stacking bin circle with the radius of 1 degree. The green dots are all the reflection points (the middle points of station pairs) within the example bin

DownLoad: Full-Size Img PowerPoint

Repeating the above procedure for each bin along the two profiles, we thus can obtain the waveforms of the MTZ reflected body waves along the two profiles (Figure 7). Both P₄₁₀P and P₆₆₀P phases can be clearly identified along the two profiles with obvious lateral variations in both amplitudes and depths implying the lateral heterogeneity of the MTZ discontinuity architectures within our study area.

Figure  7.  The tapered PWS traces along the two profiles P1 (a) and P2 (b) as shown in Figure 5. The number of stacked NCFs along each profile is given above. The two red dashed lines indicate the 410-km and 660-km depths

DownLoad: Full-Size Img PowerPoint

4.   Results and discussion

On the stacked profiles, we observe clear reflected body wave signals and significant lateral variations on the stacked depth-domain NCF profiles (Figure 7). The P₄₁₀P waveforms are relatively simple with high similarity except at the edge of two profiles. However, the P₆₆₀P waveforms are much more complicated with significant variations in both waveforms and amplitudes. For example, the P₆₆₀P waveforms change dramatically near 101°E along the east-west P1 profile and exhibit apparent north-south contrast in amplitudes along the north-south P2 profile. No clear P₅₂₀P signal can be observed on our stacked profiles though it has been recognized in SS precursor observations in many studies (see the review by Deuss, 2009). Then we will focus on discussing the accuracy and reliability of our method.

To test the stability of the final average depth-domain NCFs and the influence of ambient noise source seasonality, all the monthly stacked NCFs within the example bin shown in Figure 5 were obtained and compared with the final average trace of all available NCFs (Figure 8a). Most of the monthly stacked traces show clear reflected body wave signals and are stable and similar to the average trace. It should be noted that P₄₁₀P phase is more coherent and stable than P₆₆₀P phase in the monthly-stacked NCFs. As shown in Figure 8b, the seasonality of reflected body waves was also checked by comparing the quarterly stacked NCFs and the final average trace. Only the NCF recovered from the northern hemisphere winter data is very different from NCFs of other seasons, which shows much weaker and unstable P₆₆₀P phase but quite stable P₄₁₀P phase. The average traces including and excluding the winter trace in Figure 8b were also compared in Figure 8c and both traces are very similar. Though the winter NCF is very different from other three quarterly stacked results, it does not have significant influence on the final average NCF. As long as the time span is long enough, all the daily NCFs can be stacked to obtain stable reflected body waves from MTZ interfaces. However, if the time span is not long enough, the winter NCFs should be removed to obtain more stable reflected body wave signals in our study region. In this study, all available daily NCFs were just stacked for simplicity.

Figure  8.  The monthly and quarterly stacked depth-domain NCFs within the example bin shown in Figure 5 and their comparison with the final average NCF. (a) The monthly stacked NCFs (thin black lines) and the final average trace of all the available NCFs (the red line); (b) The quarterly stacked NCFs and their final average trace; (c) The average traces of all the available NCFs including (the red line) and excluding (the blue line) the northern hemisphere winter NCFs. Here the NCFs were stacked linearly instead of with the PWS method

DownLoad: Full-Size Img PowerPoint

The relationship between the NCF and the Green’s function is critical to investigate the true discontinuity architecture. As the major teleseismic P wave microseism sources are far from our study region (Wang et al., 2018b), the source influence can be regarded as the same for different bins in our study region. Therefore, the lateral variations of the stacked waveforms are assumed mainly from the complicated mantle structure instead of the influence of the microseism source distribution. The originally stacked NCFs, instead of the Green’s function estimated by the differentiation of NCFs (Feng et al., 2017; Nakahara, 2006), can still be utilized to investigate the lateral variations of the MTZ discontinuities. Therefore, we focus on the lateral variations of P₄₁₀P and P₆₆₀P waveforms and related discontinuity architectures instead of the absolute depths of the MTZ discontinuities.

We also test the influence of bin size on the final stacked correlations. We found that the final reflected body wave phases tend to be stable when the radius of the circle bin is larger than 0.8 degree. However, the circle bin radius is finally set to be one degree considering the balance between stability and lateral resolution. For a horizontally layered velocity model, the error caused by the time-to-depth conversion is negligible (see Figure 4). In addition, 3-D tomographic models can be easily adopted in the time-to-depth conversion to get more accurate results. However, most of the 3D tomographic models are very smooth and the amplitude of velocity anomalies is about 1% (e.g. Li et al., 2008). From our previous analysis, the 3-D velocity correction will not change the interface depths significantly (Feng et al., 2017). Also, this paper mainly focuses on the methodology instead of the detailed structure of the MTZ discontinuities. For the sake of simplification, all the time-to-depth conversion was based on the ak135 1-D Earth model.

The major MTZ discontinuities are interpreted as mineralogical phase transitions (e.g. Hirose, 2002; Ringwood, 1975; Tschauner et al., 2014). Olivine transforms into wadsleyite around 410 km depth and further transforms into ringwoodite around 520 km depth, which finally transforms into bridgmanite around 660 km depth. Olivine likely accounts for 40%–60% of the upper mantle rock composition and the remaining pyroxene and garnet still have additional phase transitions interacting with the olivine phase transitions.

For the 410-km discontinuity, only the phase transition from olivine to wadsleyite is involved. And the 410-km discontinuity can be identified in most data types as a simple single peak (Deuss, 2009). Thus, the 410-km discontinuity has been proved to be simple by both seismology and mineral physics studies, which is also consistent with our observations from ambient noise interferometry. Most of the final stacked P₄₁₀P waveforms show great consistency along the two profiles except the significant variations at the edge of the two profiles. As the number of stacked NCFs within the bins at the edge of each profile is much less, the waveform variations may be partly due to the lower signal-noise ratios (e.g. the trace at the north end of the south-north P2 profile). As shown in Figure 8a, the P₄₁₀P waveforms appear more coherent and stable than the P₆₆₀P waveforms. However, some traces show both varied P₄₁₀P waveforms and clear coherent P₆₆₀P waveforms (e.g. the three traces at the west end of the east-west P1 profile), which is unlikely just caused by low signal-to-noise ratios. Hence, significant lateral variations of the 410-km discontinuity architecture may be caused by composition and/or temperature variations.

Though the 520-km discontinuity have been previously observed in seismology, especially in SS precursors (Deuss, 2009; Shearer, 1990; 1991), no clear coherent P₅₂₀P phases can be identified on the stacked NCF section. The 520-km discontinuity seems highly variable, which is mostly observed by longer period precursor data and absent in high-frequency observations (Deuss, 2009; Shearer, 2000). The 520-km discontinuity may be broad in its sharpness and thus likely transparent to our relative high-frequency data.

The 660-km discontinuity is much more complicated compared with the 410-km discontinuity due to the coexistence of the phase transitions from ringwodite to bridgmanite and from garnet to Ca-perovskite within the similar depth range (Hirose, 2002; Tschauner et al., 2014). The 660-km discontinuity was also indicated absent from the global stacked PP precursor data but becomes visible on regional stacked results due to its highly varying character (Deuss, 2009; Deuss et al., 2006). And our observed P₆₆₀P waveforms show significant lateral variations in consistency with mineralogical predictions and other seismological observations. The amplitudes of P₆₆₀P waveforms exhibit obvious changes from 98°E to 105°E along the east-west P1 profile, whereas the waveforms of P₄₁₀P remain high consistency. Especially the small amplitudes near 101°E, where high velocity anomaly appears at about 660 km depth, may be related to the stagnant slab (Huang and Zhao, 2006; Li et al., 2008). The waveforms of P₆₆₀P along the north-south P2 profile show apparent north-south contrast: small amplitudes in the north while large amplitudes in the south bounded by 27°N. Compared with the 410-km discontinuity, much more complicated 660-km discontinuity features were also reported in our study region in two recent receiver function studies (Xu et al., 2018; Zhang et al., 2017). Both these two studies and our result indicate significant changes of the 660-km discontinuity around 101°E, that is, depression and higher amplitudes from receiver functions but smaller amplitudes of the reflected body waves from NCFs. As the upper mantle discontinuities are more likely caused by velocity gradients rather than sharp interfaces, the reflected and converted body waves may have different relative amplitude characteristics in different frequency bands. Receiver functions usually have higher frequency contents than the reflected body waves in this study. The amplitude contrast of the converted and reflected body waves may be due to the difference in data frequency content and thus their different sensitivity to interface structures. And the significant lateral variations of the 660-km discontinuity are probably attributed to the eastward subduction and the stagnant Indian oceanic plate. The lateral variations in the waveforms of P₆₆₀P may be caused by complicated discontinuity features due to lateral changes of temperature and/or composition.

5.   Conclusions

This study proposed a new technique to study the lateral variations of the MTZ discontinuities in the depth domain from dense seismic array ambient noise interferometry. First, all available continuous waveforms from dense array stations in southwest China were correlated for each station pair to generate the final time-domain NCFs, which were then converted into the depth-domain NCFs based on a 1-D layered Earth model. Then these depth-domain NCFs were stacked with the CMP stacking technique within each bin to enhance the signal to noise ratio of the retrieved signals. On the stacked profiles, clear reflected body wave signals from the MTZ discontinuities (around the 410-km and 660-km depths) can be observed with significant lateral variations. Simple and laterally coherent P₄₁₀P and much more complicated P₆₆₀P waveforms were revealed along these two profiles, which generally agree with mineralogical predictions and recent receiver function studies. As broadband observations are critical for studying the detailed structure of these discontinuities and this method can provide reliable complementary observations within the secondary microseism band (0.1 to 0.2 Hz), which is in a higher frequency band compared to PP or SS precursors. Portable dense seismic array observations in China and the world can provide new opportunities to study the interior interface structures of the Earth with ambient microseisms.

Acknowledgements

We appreciate the comments from two anonymous reviewers, which help to improve the original manuscript. Waveform data of permanent stations in SW China for this study are provided by Data Management Centre of China National Seismic Network at Institute of Geophysics, China Earthquake Administration (doi:10.11998/SeisDmc/SN, http://www.seisdmc.ac.cn). Waveform data of portable were provided by China Seismic Array Data Management Center at Institute of Geophysics, China Earthquake Administration (doi:10.12001/ChinArray.Data). This work is supported by China Earthquake Science Experiment Project, China Earthquake Administration (Nos. 2017CESE0101 and 2016CESE0201) and the National Natural Science Foundation of China (No. 41574034).