Data processing
The waveform data used in this paper are recorded by the broad-band stations from the permanent CEArray45 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.946.
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, researchers47, 48 used a 1-D model throughout the research region. In this study, to increase the accuracy, we chose the 2-D tomographic model FWEA1817 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 Algorithm22, 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 method19 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 spreading22.
The tradeoff between model parameters
To avoid falling into local minima, in the inversion code23, 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 QSEIS49 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 SPCFEM2D24 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.