Ultra-low-velocity anomaly inside the Pacific Slab near the 410-km discontinuity

doi:10.21203/rs.3.rs-1128948/v1

Download PDF

Article

Ultra-low-velocity anomaly inside the Pacific Slab near the 410-km discontinuity

https://doi.org/10.21203/rs.3.rs-1128948/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 03 May, 2023

Read the published version in Communications Earth & Environment →

Version 1

posted

You are reading this latest preprint version

Under equilibrium conditions, the mineral olivine (α) transforms to wadsleyite (β) near 410 km depth^1,2. The phase transition is a global boundary known as the “410-km discontinuity”, where seismic wave speed and density increase by 3–10 percent^3-6. We present new observations of triplicated P-waves from dense seismic arrays in northeastern China that constrain the structure of the Pacific slab near the 410-km discontinuity beneath the northern Sea of Japan. Our analysis of P-wave travel times and waveforms at periods as short as T = 2 s indicates the presence of an ultra-low-velocity layer within the cold slab, with a P-wave velocity that is at least 20% lower than in the ambient mantle and an apparent thickness of ~ 16 km along the wave path. This ultra-low-velocity layer could contain unstable material (e.g., poirierite) with reduced grain size where diffusionless transformations are favored^7-12

Near cold subducting slab, the 410-km can be elevated due to the positive Clapeyron slope of the α-β phase transition¹³. Inside the slab where temperatures are lower than 1000°C, the nucleation and growth mechanism of the α-β phase transformation¹⁴ involving temperature-dependent atomic diffusion may be inhibited¹⁵. New pathways in the form of shear transformation mechanisms have been supported by the recent discovery of a new shear-induced high-pressure olivine polymorph, the ω-olivine phase (also known as ε*-phase), in heavily-shocked meteorites^7-12. If this new diffusionless process occurs in Earth’s mantle, it will cause the seismic wave speed in the cold slab to decrease, associated with a grain size reduction.

Such a scenario could delay the arrival times of seismic waves that propagate a long distance through the slab in the vicinity of the α-β phase transition. In this paper, we analyze P waves generated by October 10, 2009, Mw 5.9 Kuril earthquake, which occurred at a depth of 114 km and was recorded by seismic stations in northeastern China. The position of the Kuril earthquake with respect to the seismic stations is ideal for studying the structure of the Pacific slab near the 410-km because P waves propagate hundreds of kilometers above or through the slab as illustrated by the Wadati-Benioff zone¹⁶ (Fig. 1a) and by the high-resolution FWEA18¹⁷ image of the mantle beneath the northwestern Pacific (Fig. 1b and 1c). In addition, the seismic stations are at epicentral distances between 1700 and 2400 km and record the triplication of P waves. The waveforms have up to three distinct pulses related to the direct wave that propagates above the 410-km and the upper interface of the slab and refracted waves that cross them. The 700-km long array aperture and the 50-km average station spacing allow us to link the traveltime and waveform characteristics of the P wave triplication to anomalous layering at a ~ 15-km scale corresponding to the width of the Fresnel zone.

We separate the seismic data into two subsets. The waveforms from stations along the source-receiver path RR’ are the reference data. Along RR’ the upper interface of the Pacific slab is located at about 450 km depth well below the 410-km. Since P waves along RR’ turn above the slab interface (Fig. 1b), we expect that the behavior of the P wave triplication can be explained by standard 1-D seismic velocity profiles of the mantle (Extended Data Fig. 1). In contrast, P waves recorded by stations along SS’ to the south of RR’ propagate several hundreds of kilometers through the high-velocity slab and are complicated by layering in the slab in the vicinity of the 410-km (Fig. 1c).

The onset times of the refracting P waves (Fig. 2a) provide the first clue that P waves along SS’ have decelerated by a low-velocity zone within the slab. First, the traveltime gradients indicate that the apparent velocities of the refracted waves are smaller for SS’ (i.e., 9.89 km/s) than for RR’ (i.e., 10.77 km/s) even though they have propagated within the high-velocity slab along SS’. Second, the travel times of the refracted waves along SS’ are similar or up to ~ 1 s longer than refracted waves along RR’ (Fig. 2a). For RR’, the refracted waves cross the 410-km discontinuity in the “normal” mantle (above the slab interface). However, for SS’, the refracted waves cross the 410-km discontinuity and propagate horizontally within the slab. Therefore, the similar travel times of the refracted waves point to a low-velocity zone inside the slab.

An alternative explanation for the delayed arrival time of the refracted waves along SS’ can be the lower velocities near the stations (e.g., at an epicentral distance of about 1500 km in Fig. 1c). Since such velocity anomaly will affect the direct waves and refracted waves to the same extent due to their almost identical ray paths at shallower depths (see the delayed arrival times for the direct waves along SS’ in Fig. 2a), the differential travel time between the refracting and direct P waves can cancel out such influence. The almost identical differential travel times between SS’ and RR’ (Fig. 2b) confirm our speculation of a low-velocity zone inside the slab.

By comparing recorded and computed waveforms of the P waves, we constrain the location, thickness, and velocity reduction of this low-velocity anomaly. We compute waveforms up to 0.5 Hz using the 2-D finite-difference method¹⁹ to avoid the high computational cost of 3-D solvers at 0.5 Hz. However, we acknowledge that complex 3-D wave propagation along the slab interface may cause some of the complexity in the P waveforms. We use the FWEA18 tomographic model of the mantle as a reference structure and perturb the layered wave-speed structure near the P-wave turning points (the region indicated by the box in Fig. 1c) to better explain the waveform. We argue that our quasi 1-D modeling of P-wave perturbations is justified given the 1-D character of the wave speed anomalies in FWEA18 near the P wave turning points, although more complex structures away from the turning points (e.g., near the source region) could contaminate the results (see Methods). We compare the synthetic and recorded waveforms after aligning each trace on the first arrival to isolate the influence of heterogeneity near the P wave turning points on the travel time differences and amplitudes of the direct and refracted waves. Waveforms are band-pass filtered between 0.02 and 1 Hz to minimize side-lobes. Since the waveforms with ray paths sub-parallel to the strike direction are very sensitive to the depths of the subducting slab, we further divided the stations along SS’ into sub-linear arrays S₁S₁’ and S₂S₂’ for source-station azimuth ranges of 262^o-264^o and 264^o -266^o, respectively.

There are two significant differences between the recorded and the synthetic P wave triplication computed for FWEA18 for the S₁S₁’ and S₂S₂’ profiles (Fig. 3a and 3b). For epicentral distances smaller than 1850 km, when the P waves propagate near the upper surface of the slab, the recorded P wave is weaker than the computed P wave because its energy is divided into two separate signals (see the waveforms for stations DNI, NE5E, NE4A along S₁S₁’ and station NE5D along S₂S₂’). Between 2000 and 2200 km, when the P waves propagate through the slab, the recorded P waveforms have three distinct peaks near about 46 s, 49 s, and 52 s. The FWEA18 synthetics match the first and the third peak (see for example stations HST, JCT, and LHTJ along S₁S₁’ and NE3A and PST along S₂S₂’) but do not predict the presence of the second peak due to the P-wave refraction through the slab’s upper surface. The waveform mismatch discrepancy at 0.5 Hz between the FWEA18 synthetics and the recorded waveforms is not surprising because the FWEA18 image is developed from waveforms with frequencies lower than 0.125 Hz (Extended Data Fig. 3) and has been parameterized to accommodate velocity variations over spatial scales larger than 50 km¹⁷. In contrast, along RR’ where seismic waves propagate in the ambient mantle (i.e., above the slab), the recorded waveforms are simpler with at most two peaks (Fig. 3c). The mismatches between the recorded data and the FWEA18 synthetics are mainly due to arrival time differences between the two peaks.

Our preferred seismic model, called “FWEA18-LVZ”, is a modification of FWEA18 in mainly two ways (see Fig. 4). First, the strength of the waveform anomaly in FWEA18 is tripled. This reduces the arrival time misfit and the stronger wave speed contrast across the slab’s upper interface predicts the double-peaked P waveform for the nearest stations (< 1850 km). The slab’s P-wave velocity anomaly of +6% in FWEA18-LVZ (increased from 2% in FWEA18) with respect to the ambient mantle is consistent with previous waveform modeling studies of the northeastern Pacific region^20,21.

Second, a low-velocity layer centered at around 400 km depth (390 - 395 km along S₁S₁’ and 402 - 412 km along S₂S₂’) predicts the second of the three P-wave signals recorded between 2000 and 2200 km that is missing in the FWEA18 synthetics. The P-wave speed reduction must be at least 20%, with an apparent thickness of up to 16 km along the wave path (Fig. 4). The length of the low-velocity layer is fixed at 600 km, based on the length of the P wave segments through the slab along the S₁S₁’and S₂S₂’ profiles (see Fig. 1). The low-velocity layer decelerates the direct wave traveling through the slab and it extends the branch (cusp) of the triplication, generating a second pulse at ~ 49 s in the waveform (see Extended Data Fig. 4). Only a low-velocity layer (rather than a high-velocity discontinuity) can generate the missing second phase and predict the travel time difference of ~ 6 s between the first (at ~ 46 s) and the third (at ~ 52 s) pulses since the average wave speed in the slab is higher in FWEA18-LVZ than in FWEA18.

We explore the thickness and the velocity reduction of the low-velocity layer by a hybrid modeling approach. For each value of the velocity reduction from ~ 5% to 90% (with intervals of 5%), we determine the optimal thickness of the low-velocity zone and the vertical velocity structure (inside the box in Fig. 1c), via the Niching Genetic Algorithm^22,23 (see Methods), to minimize the differences between recorded and synthetic waveforms. Fig 4 shows all the acceptable models that fit the recorded waveforms equally well. For instance, the match between the recorded and computed waveforms is similar for a low-velocity layer with either a thickness of 16 km and a velocity reduction of 20 % or a thickness of 2 km and a velocity reduction of 90%.

The nonuniqueness of the velocity reduction and the layer thickness is due to the minimum spatial resolution since the Fresnel zone of a P-wave (~ 0.5 Hz) is 15–20 km wide near the 410-km. This trade-off is evident in the acceptable models in the last iteration (Fig. 4), and also in the posterior distribution of all the 1,280 models (see Extended Data Fig. 5).

By experimentation, we find that the estimates of the thickness of the low-velocity layer and its velocity reduction depend slightly on the width of the low-velocity layer. However, minimum width of 200 km is necessary to produce the second P wave signal (Extended Data Fig. 6). The existence of such a low-velocity layer inside the slab is not only identified by the non-gradient-based modeling approach (with the Niching Genetic Algorithm^22,23) but also confirmed by the sensitivity kernels derived from the adjoint method²⁴ (Extended Data Fig. 7). Finally, our results do not have changed if we had chosen a different reference seismic structure (i.e., the 1-D IASP91 model instead of the FWEA18 model). Although the waveform fits are worse, a low-velocity layer centered at 370 km depth with a velocity reduction of 10–20% and a thickness of 20 km also explained the essential waveform features.

The velocity reduction is much stronger than the P-wave velocity contrast of 5% between the slab and the ambient mantle. Despite the trade-off between the layer thickness and the velocity reduction, the product of these two parameters is a constant (Extended Data Fig. 5). If we assume a velocity reduction of -5%, the low-velocity zone would be at least 50 km thick, implying that the lower half of the slab would be detached. Our waveform modeling also indicates that a -5% velocity reduction is not sufficiently strong to produce a second P wave signal. Therefore, we rule out the possibility that slab detachment could explain the low-velocity layer. Moreover, there is neither evidence for a gap in the Wadati-Benioff zone nor tomographic evidence that the Pacific slab would be weak and detached near the P-wave turning points, near 139^oE longitude and 46^o N latitude (see Fig. 1a). The presence of metastable olivine, whose P wave speed is 5% smaller than that of wadsleyite, is not a viable explanation either for the same reason.

Partial melting is inferred from the transition-zone water filter hypothesis²⁵ and it is reasonable to consider it as a possible explanation for the wave speed anomaly observed above the 410-km, as usually done to interpret low-velocity layers in equilibrated mantle regions ^26,27. However, there is no evidence for such a low velocity between 383 and 415 km depth for the reference source-receiver path RR’ where the slab is below the 410-km (Extended Data Fig. 1). It appears that the resolved low-velocity zone is confined to the Pacific slab. Thus, partial melting unlikely explains the velocity anomaly, as further confirmed by the fact that the location of the ultra-low-velocity layer along S₁S₁’ (383 to 402 km) is shallower than that along S₂S₂’ (395 to 415 km), correlating well with the relative depths of the Wadati-Benioff zone in these two regions (Fig. 1a).

Phase transformation, with volume reductions, is a mechanism for the reduction of the effective bulk modulus^28,29. In zones (e.g, thickness of 5 – 30 km^30-32 for the 410-km) where olivine (low-pressure phase) and wadsleyite (high-pressure phase) coexist, pressure perturbation by a passing seismic wave (e.g., on the order of ~ 10^-7 GPa)³³ will disrupt the equilibrium and induce the phase transformation. For a pyrolytic Earth model, the theoretical P-wave velocity reduction caused by the softening of the effective bulk modulus can be at most ~ 30%^33,34, which falls well into our estimated value between 20% and 90%. However, such softening of the elastic modulus is unexpected in the cold slab, due to the slower kinetics³⁵. We observe the ultra-low-velocity zone only inside the slab (along SS’) but not in the ambient mantle (along RR’) where softening is more likely.

Our preferred explanation of this ultra-low-velocity zone inside the slab is the existence of a layer of destabilized olivine which is associated with a significant grain size reduction: spineloid phase (β- and/or ω-olivine) would form as layers crosscutting the host olivine crystal, disorganizing the structure homogeneity and yielding a seismic wave speed reduction. In an experiment with ultrasonic interferometry and in situ X-ray diffraction, reduction of shear-wave velocity has been observed at the onset of the olivine-wadsleyite transformation and is interpreted as resulting from the existence of an intermediate spineloid phase (ω-olivine)³⁶. We propose that such destabilized layer could consist of an intermediate phase transiently existing during the olivine–wadsleyite phase transformation within the cold subducting slab, where diffusion-assisted processes are inhibited. This intermediate phase, known either as ε*-phase or ω-olivine^7-12, was at first theoretically predicted⁷, then observed in meteorites^11,12. It was recently named poirierite¹² and its transient (meta)stability could impact the rheology of the mantle^36,37. Our seismological observations could highlight the transient (meta)stability of poirierite, under substantial shear stress, within the cold sinking slab. Nonetheless, any phase transformation at relatively low temperatures necessarily induces long-lived grain size reduction^38,39, which is also known to reduce P-wave velocity⁴⁰ and is expected to contribute to the observed reduction.

It is the first time an ultra-low-velocity zone (centered around 400 km, with a wave speed reduction of at least -20%) is identified inside a subducting slab. Viscosity reduction is usually associated with low-velocity zones⁴¹, which for instance would favor deformation and fragmentation of stagnant slabs in the mantle transition zone^42-44. Further theoretical and experimental studies on grain size reduction caused by poirierite are needed to cross-check the value range we have provided in this study. And global-scale seismic studies are needed to provide more insights into olivine phase transitions in the deep Earth and planetary materials.

Anderson, D.L. 1967. Phase changes in the upper mantle. Science 157(3793), 1165–1173.
Ringwood, A.E. 1975. Composition and Petrology of the Earth's Mantle. MacGraw-Hill 618.
Niazi, M. & Anderson, D.L. 1965. Upper mantle structure of western North America from apparent velocities of P waves. Journal of Geophysical Research 70(18), 4633–4640.
Johnson, L.R. 1967. Array measurements of P velocities in the upper mantle. Journal of Geophysical Research 72(24), 6309–6325.
Flanagan, M.P. & Shearer, P.M. 1998. Global mapping of topography on transition zone velocity discontinuities by stacking SS precursors. JGR Solid Earth 103(B2), 2673–2692.
Kennett, B.L.N. & Engdahl, E.R. 1991. Traveltimes for global earthquake location and phase identification. Geophysical Journal International 105(2), 429–465.
Poirier, J. P. 1981. Martensitic olivine-spinel transformation and plasticity of the mantle transition zone. Anelasticity in the Earth 4, 113–117.
Hyde, B. G., White, T. J. & Johnson, A. W. S. 1982. Structures related to those of spinel and the β-phase, and a possible mechanism for the transformation olivine↔spinel. Zeitschrift für Kristallographie-Crystalline Materials 160(1–4), 53–62.
Madon, M. & Poirier, J.P. 1983. Transmission electron microscope observation of α, β and γ (Mg,Fe)₂SiO₄ in shocked meteorites: planar defects and polymorphic transitions. Physics of the Earth & Planetary Interiors 33(1), 31–44.
Reynard, B., Petit, P. E., Guyot, F. & Gillet, P. 1994. Pressure-induced structural modifications in Mg₂GeO₄-olivine: A Raman spectroscopic study. Physics & Chemistry of Minerals 20(8), 556–562.
Tomioka, N. & Okuchi, T. 2017. A new high-pressure form of Mg₂SiO₄ highlighting diffusionless phase transitions of olivine. Scientific Reports 7(1), 1–9.
Tomioka, N., Bindi, L., Okuchi, T., Miyahara, M., Iitaka, T., Li, Z., Kawatsu, T., Xie, X., Purevjav, N., Tani, R. & Kodama, Y. 2021. Poirierite, a dense metastable polymorph of magnesium iron silicate in shocked meteorites. Communications Earth & Environment 2(1),1–8.
Bina, C.R. & Helffrich, G. 1994. Phase transition Clapeyron slopes and transition zone seismic discontinuity topography. JGR Solid Earth 99(B8), 15853–15860.
Sung, C.M. & Burns, R.G. 1976. Kinetics of the olivine→spinel transition: implications to deep-focus earthquake genesis. Earth & Planetary Science Letters 32(2),165–170.
Kirby, S.H., Stein, S., Okal, E.A. & Rubie, D.C. 1996. Metastable mantle phase transformations and deep earthquakes in subducting oceanic lithosphere. Reviews of Geophysics 34(2),261–306.
Hayes, G.P., Moore, G.L., Portner, D.E., Hearne, M., Flamme, H., Furtney, M. & Smoczyk, G.M., 2018. Slab2, a comprehensive subduction zone geometry model. Science 362(6410),58–61.
Tao, K., Grand, S.P. & Niu, F., 2018. Seismic structure of the upper mantle beneath eastern Asia from full waveform seismic tomography. Geochemistry, Geophysics, Geosystems 19(8),2732–2763.
White, M.C., Fang, H., Nakata, N. & Ben-Zion, Y. 2020. PyKonal: a Python package for solving the eikonal equation in spherical and Cartesian coordinates using the fast-marching method. Seismological Research Letters 91(4),2378–2389.
Li, D., Helmberger, D., Clayton, R.W. & Sun, D., 2014. Global synthetic seismograms using a 2-D finite-difference method. Geophysical Journal International, 197(2),1166–1183.
Ding, X.Y. and Grand, S.P., 1994. Seismic structure of the deep Kurile subduction zone. JGR Solid Earth, 99(B12),23767–23786.
Zhan, Z., Helmberger, D.V. and Li, D., 2014. Imaging subducted slab structure beneath the Sea of Okhotsk with teleseismic waveforms. Physics of the Earth & Planetary Interiors 232,30–35.
Koper, K.D., Wysession, M.E. & Wiens, D.A. 1999. Multimodal function optimization with a niching genetic algorithm: A seismological example. Bulletin of the Seismological Society of America 89(4),978–988.
Li, J., Chen, M., Koper, K.D., Zhou, T., Xi, Z., Li, S. & Li, G., 2021. FastTrip: A Fast MPI-Accelerated 1D Triplication Waveform Inversion Package for Constraining Mantle Transition Zone Discontinuities. Seismological Research Letters 92(4), 2647–2656.
Tromp, J., Komatitsch, D. & Liu, Q. 2008. Spectral-element and adjoint methods in seismology. Communications in Computational Physics 3(1),1–32.
Bercovici, D. & Karato, S.I. 2003. Whole-mantle convection and the transition-zone water filter. Nature 425(6953),39–44.
Revenaugh, J. & Sipkin, S.A. 1994. Seismic evidence for silicate melt atop the 410-km mantle discontinuity. Nature, 369(6480),474–476.
Wei, S.S. & Shearer, P.M. 2017. A sporadic low-velocity layer atop the 410 km discontinuity beneath the Pacific Ocean. JGR Solid Earth, 122(7),5144–5159.
Anderson, D.L. 1998. The scales of mantle convection. Tectonophysics 284(1–2),1–17.
Jackson, I. 2007. In Treatise on Geophysics Vol. 2 (ed. Schubert, G.) 493–525.
Shearer, P.M. 2000. Upper mantle seismic discontinuities. Geophysical Monograph – AGU 117,115–132.
Van der Meijde, M., Marone, F., Giardini, D. & Van der Lee, S. 2003. Seismic evidence for water deep in Earth's upper mantle. Science 300(5625),1556–1558.
Ricard, Y., Matas, J. & Chambat, F. 2009. Seismic attenuation in a phase change coexistence loop. Physics of the Earth & Planetary Interiors 176(1–2),124–131.
Li, L. & Weidner, D.J. 2008. Effect of phase transitions on compressional-wave velocities in the Earth’s mantle. Nature 454(7207), 984–986.
Ricard, Y., Matas, J. & Chambat, F. 2009. Seismic attenuation in a phase change coexistence loop. Physics of the Earth & Planetary Interiors 176(1–2),124–131.
Durand, S., Chambat, F., Matas, J. & Ricard, Y., 2012. Constraining the kinetics of mantle phase changes with seismic data. Geophysical Journal International 189(3),1557–1564.
Perrillat, J.P., Chantel, J., Tauzin, B., Jonfal, J., Daniel, I., Jing, Z. & Wang, Y. 2017, December. Acoustic Velocities Across the Olivine-Wadsleyite-Ringwoodite Transitions and the Seismic Signature of the 410 km Mantle Discontinuity. In AGU Fall Meeting Abstracts, DI13A-0279.
Ferrand, T. P. & Deldicque, D. 2021. Reduced Viscosity of Mg₂GeO₄ with Minor MgGeO₃ between 1000 and 1150°C Suggests Solid-State Lubrication at the Lithosphere–Asthenosphere Boundary. Minerals 11(6), 600.
Green II, H. W., Shi, F., Bozhilov, K., Xia, G. & Reches, A. Z. 2015. Phase transformation and nanometric flow cause extreme weakening during fault slip. Nature Geoscience 8(6), 484–489.
Incel, S., Hilairet, N., Labrousse, L., John, T., Deldicque, D., Ferrand, T., Wang, Y., Renner, J., Morales, L. & Schubnel, A. 2017. Laboratory earthquakes triggered during eclogitization of lawsonite-bearing blueschist. Earth & Planetary Science Letters 459,320–331.
Faul, U. H. & Jackson, I. 2005. The seismological signature of temperature and grain size variations in the upper mantle. Earth & Planetary Science Letters 234(1–2), 119–134.
Stern, T.A., Henrys, S.A., Okaya, D., Louie, J.N., Savage, M.K., Lamb, S., Sato, H., Sutherland, R. & Iwasaki, T. 2015. A seismic reflection image for the base of a tectonic plate. Nature 518(7537),85–88.
Lee, C. & King, S.D. 2011. Dynamic buckling of subducting slabs reconciles geological and geophysical observations. Earth & Planetary Science Letters 312(3–4),360–370.
Zhou, Y. 2018. Anomalous mantle transition zone beneath the Yellowstone hotspot track. Nature Geoscience 11(6), pp.449–453.
Faccenna, C., Funiciello, F., Civetta, L., D Antonio, M., Moroni, M. & Piromallo, C. 2007. Slab disruption, mantle circulation, and the opening of the Tyrrhenian basins. GSA special papers 418,153.

Acknowledgments

We thank M. Chen (passed away on 2021 July 18), who built the connection between J. Li and T. P. Ferrand through an online seminar, during the COVID-19 pandemic. Seismic records used in this study came from the CEArray and the NECESSArray, and we thank the team members for their deployments. We thank T. Bao, J. Ning, Z. Xi, S. van der Lee for discussions. We thank the Institute for Cyber-Enabled Research (ICER) at Michigan State University for providing high-performance computing resources. This research was supported by the NSF grant 1802247 and startup fund of Min Chen at Michigan State University and by the research fellowship of T. P. Ferrand (Alexander von Humboldt Foundation).

Author contributions

J. Li processed the seismic data and developed the inversion code. T. P. Ferrand, J. Li, M. Chen, T. Zhou, and J. Ritsema analyzed the results. T. P. Ferrand proposed the preferred interpretation and discussed with J. Li and L. Stixrude on the relevance of alternative explanations. J. Li. and J. Ritsema took the lead in writing the manuscript, and all authors discussed the results and edited the manuscript.

Competing interests

The authors declare no competing interests.

Data processing

The waveform data used in this paper are recorded by the broad-band stations from the permanent CEArray⁴⁵ and the temporary NECESSArray (the NorthEast China Extended Seismic Array) in northeast China. The intermediate-depth (114 km) earthquake used in this paper occurred on 10/12/2009, with a moment magnitude of 5.9⁴⁶.

We applied a Butterworth filter (first-order, zero-phase shift) with a frequency band from 0.05 to 1 Hz to the displacement data, after removing the instrument response. For a better illustration of the waveform (Fig. 3), we applied a reduced slowness of 10.1 s/km to the time axis to allow stations at different epicentral distances to appear in the same time window.

Waveform inversion

Since triplications have the advantage of minimizing the influence from the shallower part and emphasizing structures near the turning points of the ray paths, researchers^{47, 48} used a 1-D model throughout the research region. In this study, to increase the accuracy, we chose the 2-D tomographic model FWEA18¹⁷ as the background P-wave model in a certain cross-section.

During the inversion, we kept the model outside of the inversion box (dashed box in Fig. 1c) unchanged, and only modified the region inside the box by adding 1-D model perturbations. Since we only modify the velocity structure of FWEA18 near the P-wave turning points, unconstrained parts of FWEA18 (with probably underestimated wave speed) near the source region might lead to some overestimation of the inverted values near the turning points. Nevertheless, the current P-wave velocity anomaly in the FWEA18 model is not large enough to predict the double-peaked P wave.

For the model setup, for anchor points away from the interface, we set the interval to be 20 km and allowed the P-wave speed at each anchor point to change within 0.5 km/s. In addition, there are also anchor points immediately on the discontinuity to capture the velocity jump. The location of the interface and thickness are also free parameters, with a relatively large variable range of 35 km and 50 km, respectively. The variable range of the velocity in the ultra-low-velocity zone will be discussed in the following ‘tradeoff’ section.

The inversion code is based on the non-gradient-based Niching Genetic Algorithm^{22, 23}, which searches the model space via massive forward modeling. As for the misfit window, we choose a continuous window from 39 to 57 seconds (reduced time). We note that before calculating the L2 norm of the differences between the data and synthetics, we aligned them by cross-correlation to emphasize the relative information between triplicated waveforms and to mitigate baseline shifts.

We selected a GPU-based 2-D finite-difference method¹⁹ for our forward modeling tool to reduce the computational costs for the non-gradient-based inversion approach. For one 2-D simulation (up to 1 Hz), the computation takes 6 seconds on one NVIDIA V100 GPU card. We note that several corrections, e.g., out-of-plane, point source excitation, and Earth-flattening, have been incorporated to better account for the 3-D wavefield spreading²².

The tradeoff between model parameters

To avoid falling into local minima, in the inversion code²³, we divided the model population into subpopulations. And a penalty term for the model similarities between subpopulations is applied. In each inversion, we performed 100 iterations. And in each iteration, there are 80 models in 10 subpopulations, to ensure model diversity. In total, we searched 80,000 models in each inversion.

To explore all possible candidates in the model space, for each region, we performed 16 inversions, instead of one. For each inversion, the location and thickness of the low-velocity layer are still free parameters. The only difference is that we forced the velocity reduction to be fixed (i.e., from -0.5 km/s to -8.0 km/s with intervals of 0.5 km/s) to reduce the total number of free parameters at the expense of performing more inversions. With randomly distributed model parameters as inputs, most of the models in the last iteration are clustered and follow a curve (see Extended Data Fig. 5), indicating the trade-off between velocity reduction and layer thickness is observed. Nevertheless, the product of these two parameters is constant (Extended Data Fig. 5).

Finally, we further examined the waveform fitting of all the models in the last iteration and selected the final acceptable models (Fig. 4).

Robustness of the ultra-low-velocity layer.

For the inversion region (dashed box in Fig. 1c), the vertical height is about 250 km, which is large enough (about twice the thickness of the slab). In the horizontal direction, we set the length to be 600 km, which contains most of the ray paths near the turning points. To test how the inversion region affects the results, we conducted extra inversions with different horizontal lengths of 200 km, 300 km, 400 km, and 800 km, respectively. Results show that when the length is too short (e.g., 200 km), no model can fit the data. For other lengths, the inverted models are consistently pointing to the existence of the ultra-low-velocity zone, with slightly different locations and amplitudes (Extended Data Fig. 6).

To test the dependency of our results on the initial model of FWEA18, we performed another 1-D inversion. In this case, we chose the QSEIS⁴⁹ code as the forward modeling tool. In the inverted 1-D model sets, a significant velocity reduction (-10% to -20% when the layer thickness is about 20 km) is also observed (Extended Data Fig. 8).

Besides the non-gradient-based inversion approach, we also used SPCFEM2D²⁴ to calculate the sensitivity kernels. In practice, we calculate the adjoint source for certain stations (NE3A, PST, and LYT for the northern region. HST, JCT, LHTJ, and QYU for the southern region). With the adjoint sources, we run the adjoint simulation and then get the misfit function gradient, also known as the sensitivity kernel (Extended Data Fig. 7).

Data availability

Raw seismic data for NECESSArray (the NorthEast China Extended Seismic Array) are available at the Data Management Center of the Incorporated Research Institutions for Seismology (http://www.iris.edu/dms/nodes/dmc) under network YP. CEArray data from the Data Management Center of China National Seismic Network at Institute of Geophysics, China Earthquake Administration cat be requested at [email protected].

Extended references

45 Zheng, X.F., Yao, Z.X., Liang, J.H. & Zheng, J. 2010. The role played and opportunities provided by IGP DMC of China National Seismic Network in Wenchuan earthquake disaster relief and researches. Bulletin of the Seismological Society of America100(5B), 2866-2872.

46 Ekström, G., Nettles, M. & Dziewoński, A.M. 2012. The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. Physics of the Earth & Planetary Interiors200, 1-9.

47 Grand, S.P. and Helmberger, D.V. 1984. Upper mantle shear structure of North America. Geophysical Journal International76(2), 399-438.

48 Brudzinski, M.R. & Chen, W.P. 2000. Variations in P wave speeds and outboard earthquakes: evidence for a petrologic anomaly in the mantle transition zone. Journal of Geophysical Research: Solid Earth105(B9), 21661-21682.

49 Wang, R. 1999. A simple orthonormalization method for stable and efficient computation of Green's functions. Bulletin of the Seismological Society of America89(3), 733-741.

There is NO Competing Interest.

ExtendedData.docx

Download PDF

Journal Publication

published 03 May, 2023

Read the published version in Communications Earth & Environment →

Version 1

posted

You are reading this latest preprint version

Ultra-low-velocity anomaly inside the Pacific Slab near the 410-km discontinuity

Status:

Journal Publication

Version 1

Abstract

Figures

Full Text

References

Declarations

Acknowledgments

Author contributions

Competing interests

Methods

Data processing

Waveform inversion

The tradeoff between model parameters

Robustness of the ultra-low-velocity layer.

Data availability

Extended references

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1