https://orcid.org/0000-0001-9636-3268
SUMMARY
In seismology, wavefield injection refers to the propagation of seismic waves generated by remote sources into local domains bounded by enclosed surfaces. The simulations of wavefield injection, primarily focused on the interaction between incoming seismic waves and local structures, are key to earthquake hazard modelling and full-waveform seismic tomography using tele-seismic waves. In this paper, we show that simulating wavefield injection is equivalent to solving the wave equation subject to interface discontinuity conditions. To provide a general framework to study wavefield injection, we formally define the interface discontinuity problem, and discuss its representation theorem and uniqueness. We also develop an efficient interface-discontinuity-based numerical algorithm to solve the wavefield injection problem through implementations of spectral-element methods, and show with numerical examples that wavefield injection can be accurately simulated at different scales with this algorithm. Under this framework, we draw connections with previously proposed wavefield injection algorithms/hybrid methods, and clarify several theoretical questions on wavefield injection from previous research. We demonstrate the efficiency and accuracy of our approach through wavefield injection examples at local and continental scales. Furthermore, we illustrate the applicability of the interface discontinuity approach to performing kinematic fault simulations through a numerical example.
1 INTRODUCTION
Simulations of seismic wave propagation are essential for various fields in seismology, such as full-waveform seismic tomography (e.g. Fichtner et al. 2010; Tape et al. 2010; French et al. 2013; Zhu et al. 2015) and earthquake hazard modelling (e.g. Graves 1998; Pitarka 1999; Zhao et al. 2007). Simulation tools that are based on analytical solutions of the seismic wave equation in simple media (e.g. Takeuchi & Saito 1972; Kennett & Kerry 1979; Wang 1999; Masters et al. 2011; Wang et al. 2017), although highly efficient, cannot adapt to real-world complexities of interest, such as surface topography and heterogeneity. With the rapid advances in numerical methods and high-performance computing facilities, many solvers of seismic wave propagations are developed based on finite-difference method (Graves 1996), spectral-element method (Komatitsch & Tromp 1999b; Gokhberg & Fichtner 2016; Cupillard & Capdeville 2018), discontinuous Galerkin method (Dumbser & Käser 2006) and distributional finite difference method (Masson 2022; Lyu et al. 2024; Masson et al. 2024), achieving accurate and efficient seismic wave simulations with complex geometry and 3-D heterogeneity at all scales.
Since it is still computationally too expensive to simulate seismic waves at high frequency within the whole Earth, it is useful and feasible to simulate the interaction between seismic waves and local structures. For example, local crustal variations (Wang et al. 2016, 2022; Xu et al. 2023; He et al. 2024), lithospheric structures (Clouzet et al. 2018; Wang et al. 2021) and lowermost mantle heterogeneities (Lin et al. 2019) such as large low-shear-velocity provinces and ultra-low velocity zones interact with tele-seismic incoming waves and the scattered waves can be used for full-waveform seismic tomography (often referred to as box tomography; Masson & Romanowicz 2017). Local sedimentary basins can strongly affect ground motions and are of interest in earthquake hazard modellings (e.g. Yoshimura et al. 2003). In those studies, the numerical computations are confined to a localized study region in which local structures exist, whereas the incoming or outgoing waves outside the study region are approximately computed using some efficient solvers, often ones based on analytical solutions, assuming some simplified (e.g. homogeneous or 1-D) background structural model. Such strategy is often called hybrid methods. Depending on the problem setup, the hybrid simulation may involve either wavefield injection, which propagates a background wavefield into a local domain, or wavefield extrapolation, which propagates a perturbed wavefield away from the local domain. Since wavefield extrapolation can be viewed as injecting the wavefield from the local domain into the external domain, its derivations and applications can be carried out similar to those for wavefield injection. Therefore in this study, we focus on wavefield injection.
Early schemes for wavefield injection, such as Shtivelman (1985), Wen & Helmberger (1998), Robertsson & Chapman (1999), Opršal et al. (2002), Zhao et al. (2008), etc., modify the finite-difference scheme to allow the coupling of a background wavefield and the finite-difference solver. Capdeville et al. (2003) constructed a Dirichlet-to-Neumann operator to couple in real time the spectral-element solver in the outer shell of the Earth with the normal-mode solution in the inner sphere to accurately calculate the full wavefield including high-order scatterings. Bielak & Christiano (1984) proposed to solve for the total and scattered wavefield inside and outside the local domain, respectively, as a general strategy to simulate wavefield injection, which was further elaborated by Opršal et al. (2009) using the concept of interface discontinuity. This strategy was further developed by subsequent works (e.g. Bielak et al. 2003; Yoshimura et al. 2003) in the discretized finite-element form of the seismic wave equation and was named the domain reduction method. Leng et al. (2020) followed a similar procedure, but derived the algorithm in the weak form of the seismic wave equation. An independent approach was developed by Masson et al. (2014) based on the concept of force mirrors, which are forces that can reconstruct the background wavefield in the local domain. Force mirrors are equivalent to the boundary-value representation theorem (chap. 2.5; Aki & Richards 2002) in the continuous form, but they can accommodate different numerical methods and geometries in discretized forms. Lin et al. (2019) derives the injection algorithm directly using the boundary-value representation theorem. Lyu et al. (2022) proposed an injection algorithm by combining the physical and numerical representation theorem. Another category of methods derive an equivalent boundary condition on the injection boundary by imposing the Stacey absorbing boundary condition (Clayton & Engquist 1977) on the scattered wavefield (e.g. Monteiller et al. 2013, 2015; Tong et al. 2014, 2015). Although the methods based on equivalent boundary conditions are in general easy to derive and implement, they cannot accommodate more accurate absorbing boundaries such as the perfectly matched layer (PML; Bérenger 1999; Komatitsch & Martin 2007) and therefore may suffer from artifical reflections of the scattered waves. Hybrid simulations can be performed by combining these methods with background wavefield generated by 1-D solvers, such as generalized ray theory (Wen & Helmberger 1998), normal-mode solution (Capdeville et al. 2003), direct solution method (Monteiller et al. 2013), frequency-wavenumber method (FK; Tong et al. 2014, 2015) and AxiSEM (Monteiller et al. 2021; Pienkowska et al. 2021), or 3-D solvers, such as AxiSEM3D (Leng et al. 2020) and SPECFEM3D_GLOBE (Clouzet et al. 2018).
Although the techniques to simulate wavefield injection are already relatively mature, several theoretical questions remain to be clarified. First, how are the different methods for wavefield injection mentioned above related? Second, what information on the background wavefield is needed to solve the wavefield injection problem? Previously proposed methods commonly require the displacement and traction of the background wavefield on the injection interface to simulate wavefield injection. Is it possible to simulate wavefield injection with only displacement or traction on the injection interface? Finally, the methods based on the boundary-value representation theorem (Lin et al. 2019; Pienkowska et al. 2021; Lyu et al. 2022) or force mirrors (Masson et al. 2014) are designed to reconstruct the wavefield in the local domain using its boundary values, while, in practice, the boundary value of the background wavefield is used because the true wavefield is unknown. Can the true wavefield be accurately reconstructed using the boundary values of the background wavefield? To clarify these questions, we formally define the interface discontinuity problem of the wave equation, and recognize that the wavefield injection problem is a special case of the interface discontinuity problem. We discuss the representation theorem and the uniqueness of the interface discontinuity problem, and derive a numerical scheme to solve the interface discontinuity problem using the spectral-element method. Using the framework of the interface discontinuity problem, we reveal that previously proposed methods based on domain reduction, force mirrors, representation theorem and equivalent boundary conditions are theoretically related as they all solve the same interface discontinuity problem. Moreover, through the uniqueness of the interface discontinuity problem, we show that the displacement and traction of the background wavefield are both required in order to solve the wavefield injection problem. Finally, based on the similarity between the representation theorems of the interface discontinuity problem and the boundary value problem, we state that replacing the boundary values of the true wavefield with those of the background wavefield, as in Masson et al. (2014), Lin et al. (2019), Pienkowska et al. (2021) and Lyu et al. (2022), can indeed reconstruct the true wavefield in the local domain.
The rest of the manuscript is organized as follows: In Section 2, we formally define the wavefield injection problem and the interface discontinuity problem, and recognize that the former is a special case of the latter. In Section 3, the representation theorem of the interface discontinuity problem is studied and compared with the boundary-value representation theorem. In Section 4, the uniqueness of the interface discontinuity problem is discussed. In Sections 5 and 6, we derive a numerical scheme to solve the interface discontinuity problem using the spectral-element method, and show its effectiveness in simulating wavefield injection with numerical examples. In the discussion Section 7, we discuss the relationships between our method and previously proposed methods to simulate wavefield injection, and illustrate the potential usefulness of our numerical scheme to kinematic fault simulations.
Note that in the discussions throughout this paper, we assume that the injection interface is honoured by the mesh in the local domain, that is, the interface is defined only by points on element boundaries and does not run through any element. It is worth mentioning that the force mirror method (Masson et al. 2014) is not subject to such restriction, and therefore, with additional storage cost, offers more flexibility than other methods when the local and background solvers do not share a common interface. However, in most cases, we can expect the background solver to be able to accurately compute displacement and traction at any point in its simulation domain via interpolation, and therefore, as long as the injection interface is in the simulation domain of the background solver, the displacement and traction of the background wavefield can be easily obtained.
2 WAVE EQUATION, WAVEFIELD INJECTION AND THE INTERFACE DISCONTINUITY PROBLEM
We discuss the wavefield injection in the context of linearly elastic wave propagation. The wave equation that governs the evolution of displacement |$\boldsymbol {u}$| in space |$\boldsymbol {x}$| and time t due to a body force |$\boldsymbol {f}$| is given by
where |$\rho$| is density, and |$\boldsymbol {C}$| is the fourth-order elastic tensor. The wave equation eq. (1) is solved in some global domain|$\oplus$|, which could be the whole Earth. Certain boundary condition needs to be assigned at the boundary |$\partial \oplus$|, most commonly a traction-free boundary condition at the Earth’s surface
where |$\hat{\boldsymbol {n}}$| is the normal vector on |$\partial \oplus$| pointing outward.
Suppose that we are interested in obtaining the solution |$\boldsymbol {u}$| of eq. (1) subject to eq. (2) only in a source-free local domain|$\Omega \subset \oplus$|. With classical numerical methods, numerical computations need to be performed in a computational domain that includes at least both the local domain |$\Omega$| and the source |$\boldsymbol {f}$|, with appropriate absorbing boundary conditions if the computational domain is not the global domain |$\oplus$|. If the source |$\boldsymbol {f}$| is far away from the local domain |$\Omega$|, then the computational domain must be much larger than |$\Omega$|, leading to large computational costs. Instead of directly solving eq. (1), we seek to reconstruct |$\boldsymbol {u}$| in |$\Omega$| by carrying out numerical computations only in |$\Omega$| (and possibly also a small neighbourhood surrounding |$\Omega$|), with the help of a known background wavefield|$\boldsymbol {u}^0$|. This background wavefield satisfies the wave equation with background parameters|$\rho ^0$| and |$\boldsymbol {C}^0$| as
and a traction-free boundary condition
Here and throughout this paper, we use |$\overline{A}:= A\bigcup \partial A$| to denote the closure of set A, and |$A\setminus B:=\lbrace \boldsymbol {x}|\boldsymbol {x}\in A, \boldsymbol {x}\not\in B\rbrace$| to denote the set difference between A and B. Note that the background wavefield |$\boldsymbol {u}^0$| does not necessarily need to be well-defined in the local domain |$\Omega$|. For example, local topography exists in |$\Omega$| may not be in the background model that generates wavefield |$\boldsymbol {u}^0$|. Eqs (3) and (4) together cannot uniquely determine |$\boldsymbol {u}^0$| because the boundary in eq. (4) does not entirely enclose the domain |$\oplus \setminus {\bar{\Omega }}$|. The background wavefield can be arbitrarily chosen from the infinite solution space to eqs (3) and (4), and can even be the total wavefield, eqs (1) and (2), as will be discussed in Section 3.
If the following assumptions hold:
- the true and background parameters are identical outside the local domain(5)$$\begin{eqnarray} \rho = \rho ^0, \qquad \boldsymbol {C} = \boldsymbol {C}^0 \qquad \text{in}\ \overline{\oplus \setminus \Omega } \end{eqnarray}$$
- the local domain is source-free(6)$$\begin{eqnarray} \boldsymbol {f} = 0 \qquad \text{in}\ \overline{\Omega }, \end{eqnarray}$$
then reconstructing |$\boldsymbol {u}$| in |$\Omega$| using |$\boldsymbol {u}^0$| is called the wavefield injection problem. Conversely, in the wavefield extrapolation problem, the wavefield propagates away from the local domain |$\Omega$|. In this paper, we focus on the wavefield injection problem, but all discussions and derivations can be easily adapted to the wavefield extrapolation problem by swapping the local and external domains |$\Omega$| and |$\oplus \setminus \Omega$|.
As pointed out earlier, a source that is far away from the local domain is the cause of the large computational burden. Therefore, to restrict the numerical computation in |$\Omega$|, it is key to eliminate the source term f from eq. (1). To achieve this, an intuitive manipulation is to subtract eq. (1) by eq. (3). Due to eq. (5), the subtraction can be carried out outside |$\Omega$|:
and the boundary condition can be also subtracted:
|$\boldsymbol {u}-\boldsymbol {u}^0$| arises from local structures in |$\Omega$|, and therefore is commonly referred to as the scattered wavefield. The true wavefield |$\boldsymbol {u}$| can be viewed as the sum of the background and scattered wavefield, and is therefore also called the total wavefield.
Eq. (7) implies that the scattered wavefield |$\boldsymbol {u}-\boldsymbol {u}^0$| satisfies the source-free wave equation outside the local domain |$\Omega$|. Due to eq. (6), we know that the total wavefield |$\boldsymbol {u}$| satisfies the source-free wave equation inside the local domain |$\Omega$|. Using the same treatment as Bielak & Christiano (1984), we write a unified source-free equation in the entire global domain |$\oplus$| by defining a new wavefield
|$\boldsymbol {w}$| satisfies the source-free equation in the entire global domain |$\oplus$|
with boundary condition
However, the displacement and traction of wavefield |$\boldsymbol {w}$| is discontinuous across the injection interface|$\Gamma$| between |$\Omega$| and |$\oplus \setminus \Omega$|
where |$\Big |_{\Gamma +}^{\Gamma -}$| denotes the difference of displacement or traction across |$\Gamma$|, the values inside |$\Omega$| minus those outside, and |$\hat{\boldsymbol {n}}$| is the normal vector on |$\Gamma$| pointing away from |$\Omega$|. If |$\Omega$| is strictly inside |$\oplus$|, then |$\Gamma = \partial \Omega$| (Fig. 1b). However, if |$\oplus$| and |$\Omega$| have a common boundary, then |$\Gamma = \partial \Omega \setminus \partial \oplus$| (Fig. 1a).
Geometry of the local domain |$\Omega$| (blue areas) and the injection interface |$\Gamma$| (red lines). The thick black lines represent the free surface of the Earth |$\partial \oplus$|. In case (a), the boundary of the local domain and the free surface have an intersection. The injection interface |$\Gamma$| is defined as the part of |$\partial \Omega$| that is not on the free surface. In case (b), the boundary of the local domain does not intersect with the free surface, and the injection interface is the same as the boundary of the local domain |$\partial \Omega$|.
Inspired by eqs (10)–(13), we define the interface discontinuity problem as the problem of solving a source-free wave equation eq. (10) subject to certain boundary condition eq. (11) and discontinuity on interface |$\Gamma$|
where |$\hat{\boldsymbol {n}}$| points from |$\Gamma -$| towards |$\Gamma +$|. Note that eqs (14) and (15) mean that the displacement and traction is continuous everywhere other than on the interface |$\Gamma$|. Clearly, the wavefield injection problem is a special case of the interface discontinuity problem with |$\Delta \boldsymbol {w} = \boldsymbol {u}^0\Big |_{\Gamma }$| and |$\Delta \boldsymbol {T} = \boldsymbol {n}\cdot (\boldsymbol {C}:\nabla \boldsymbol {u}^0)\Big |_{\Gamma }$| in eqs (14) and (15). Interface |$\Gamma$| always encloses some local domain |$\Omega$| for wavefield injection problems, but it can be an arbitrary surface in |$\oplus$| for general interface discontinuity problems.
The kinematic fault simulation problem, which computes the seismic wavefield due to a prescribed slip on the fault (e.g. Magnoni et al. 2013), can also be viewed as a special case of the interface discontinuity problem with displacement discontinuity |$\Delta \boldsymbol {w}$| equal to the slip history and traction discontinuity |$\Delta \boldsymbol {T}=0$|. In the discussion Section 7.2, we illustrate that the numerical scheme we propose to solve the interface discontinuity problem can directly be used for kinematic fault simulation. The general interface discontinuity problem does not correspond to any physically meaningful scenario in seismology, but the interface discontinuity of electromagnetic waves, which describes the surface charges and currents (e.g. Law et al. 2020), is widely studied and used.
3 THE REPRESENTATION THEOREM OF THE INTERFACE DISCONTINUITY PROBLEM
The representation theorem of the boundary value problem (chap. 2.5; Aki & Richards 2002), which states that the wavefield can be represented by boundary values of displacement and traction, is well-known in seismology. In this section, we discuss the representation theorem of the interface discontinuity problem, and compare it with that of the boundary value problem.
Due to the close relationship with the fault kinematic problem, the derivation of the representation theorem for the interface discontinuity condition problem is almost identical to that of the well-known source representation theorem (eq. 3.2; Aki & Richards 2002). Therefore, we directly write the equation here and refer the readers to chap. 3.1 of Aki & Richards (2002) for the detailed derivation. Suppose that wavefield |$\boldsymbol {w}$| satisfies the source-free wave equation eq. (10), traction-free boundary condition eq. (11) and interface discontinuity eqs (14) and (15), |$\boldsymbol {G}(\boldsymbol {x}, t; \boldsymbol {x}^\prime )$| is the Green’s function in |$\oplus$| with the true parameters |$\rho$| and |$\boldsymbol {C}$| and satisfying the traction-free boundary condition on |$\partial \oplus$|, then
Considering the relationship between the wavefield injection problem and the interface discontinuity problem discussed in Section 2, we have
Note that this equation states that the total wavefield in the local domain can be represented by the displacement and traction of the background wavefield on the injection interface. This is different from the boundary-value representation theorem, which represents the total wavefield by the boundary values of the total wavefield itself:
Here |$\int _{\partial \Omega }$| can be replaced by |$\int _\Gamma$|, because either |$\Gamma =\partial \Omega$| or |$\hat{\boldsymbol {n}}\cdot (\boldsymbol {C}:\nabla \boldsymbol {u}) = \hat{\boldsymbol {n}}\cdot (\boldsymbol {C}:\nabla \boldsymbol {G}^T) = 0$| on |$\partial \Omega \setminus \Gamma \subset \partial \oplus$|. In fact, eq. (18) is a special case of eq. (17), because the total wavefield itself can be viewed as a special choice of the background wavefield, due to its satisfaction to eqs (3) and (4).
Strictly speaking, the wavefield injection problem should be solved using eq. (17) instead of eq. (18), because the boundary values of the true wavefield |$\boldsymbol {u}(\boldsymbol {x}^\prime ,t^\prime )$| and |$\boldsymbol {n}\cdot (\boldsymbol {C}:\nabla \boldsymbol {u}(\boldsymbol {x}^\prime ,t^\prime ))$| are unknown. Several previous studies (e.g. Masson et al. 2014; Lin et al. 2019; Pienkowska et al. 2021; Lyu et al. 2022) derived their equations based on the boundary-value problem of the total wavefield, but replaced these boundary values with those of the background wavefield in their implementations, which cannot assume to be valid without proof. By comparing eqs (17) and (18), it is clear that such replacement can still reconstruct the total wavefield in the local domain. However, theoretically, the wavefield injection problem should be viewed as an interface discontinuity problem rather than a boundary value problem.
4 THE UNIQUENESS OF THE INTERFACE DISCONTINUITY PROBLEM
We show in this section that the interface discontinuity eqs (14) and (15) can uniquely determine the wavefield. Similar to the proof of the uniqueness of the boundary value problem, we assume |$\boldsymbol {w}_1$| and |$\boldsymbol {w}_2$| both satisfy the source-free wave equation eq. (10), the traction-free boundary condition eq. (11), zero initial condition, and the interface discontinuity eqs (14) and (15). To show the uniqueness, it suffices to show |$\delta \boldsymbol {w}\equiv \boldsymbol {w}_1-\boldsymbol {w}_2$| must be zero.
We define the energy of wavefield |$\delta \boldsymbol {w}$| as
in which the first and second term on the right-hand side correspond to the kinetic and elastic energy, respectively. Obviously, the energy is always non-negative, and vanishes if and only if the wavefield is zero everywhere.
Since |$\boldsymbol {w}_1$| and |$\boldsymbol {w}_2$| both satisfy the source-free wave equation eq. (10), their difference |$\delta \boldsymbol {w}$| satisfy the same equation
Multiplying both sides of eq. (20) by |$\dfrac{\partial \delta \boldsymbol {w}}{\partial t}$| using the dot product, integrating over the global volume and applying integration by parts, we obtain
If |$\boldsymbol {w}_1$| and |$\boldsymbol {w}_2$| both satisfy the zero-traction boundary condition eq. (11) and the interface discontinuity eqs (14) and (15), then their difference |$\delta \boldsymbol {w}$| satisfies the zero-traction boundary condition
and has no interface discontinuity
Eq. (22) guarantees that the integral on |$\partial \oplus$| in eq. (21) is zero, and eqs (23) and (24) together guarantee the integral on |$\Gamma$| is zero. This leads to |$\dfrac{\mathrm{d}E[\delta \boldsymbol {w}]}{\mathrm{d}t} = 0$|. With zero initial condition, we have |$E[\delta \boldsymbol {w}](t) = E[\delta \boldsymbol {w}](0) = 0$|, indicating |$\delta \boldsymbol {w}$| must be zero. Therefore, |$\boldsymbol {w}_1$| must be the same as |$\boldsymbol {w}_2$|, meaning that the interface discontinuity on displacement and traction, together with boundary and initial conditions, determines a unique solution of the wave equation eq. (10).
Eq. (21) also suggests that displacement and traction discontinuities must be simultaneously specified to guarantee uniqueness. To see this, suppose that the displacement discontinuity or traction discontinuity is unspecified on |$\Gamma$|, then one of eqs (23) and (24) may not be satisfied, meaning the second integral in eq. (21) may be non-zero, hence two different solutions may exist.
By considering the relationship between the wavefield injection problem and the interface discontinuity problem, the uniqueness statement above implies that the displacement and traction of the background wavefield on the injection interface are both needed to simulate wavefield injection. This conclusion also suggests that the wavefield injection problem and the boundary value problem are fundamentally different: the former requires both displacement and traction information on the injection interface, whereas the latter, as is well-known, can be uniquely solved with displacement or traction information on the boundary.
5 NUMERICAL SCHEME OF THE INTERFACE DISCONTINUITY PROBLEM BASED ON THE SPECTRAL-ELEMENT METHOD (SEM)
In this section, we derive an SEM scheme for the interface discontinuity problem, as defined by eq. (10) subject to eqs (11), (14) and (15). To obtain the weak form, we multiply both sides of eq. (10) by a test function |$\phi$|, and integrate over |$\oplus$|. Using integration by parts, and considering the traction-free boundary condition eq. (11) and the traction discontinuity eq. (15), we obtain
Note that here we require |$\phi$| to be continuous across |$\Gamma$|. To solve the interface discontinuity problem, we need to find |$\boldsymbol {w}$| that satisfies eq. (25) for any admissible test function |$\phi$|, subject to the displacement discontinuity (14) on |$\Gamma$|. Based on SEM, the computational domain can be first discretized into elements with Gauss–Legendre–Lobatto (GLL) points, and the wavefield |$\boldsymbol {w}$| is approximated on these GLL points.
We suppose that the space is discretized in a way that the interface is honoured, that is, the interface is defined only by points on element boundaries and does not run through any element. We use |$\mathsf {w}_{\Gamma }^-$| and |$\mathsf {w}_{\Gamma }^+$| to denote the displacements for points on the negative and positive side of |$\Gamma$|, and |$\mathsf {w}_i$| for points not on |$\Gamma$|, then the discretized version of eq. (25) can be written as:
in which |$\mathsf {M}_{\Gamma \Gamma }^-$|, |$\mathsf {M}_{\Gamma \Gamma }^+$| and |$\mathsf {M}_{ii}$| are positive-definite diagonal matrices in the SEM formulation. The displacement discontinuity eq. (14) can be written as
Eqs (26) and (27) together form a solvable system of ordinary differential equations. Directly solving this system requires treating each point on the interface as two separated GLL points, representing the positive and negative sides of the interface, respectively. This leads to a split-node scheme (e.g. Kaneko et al. 2008), which needs some special treatment to be implemented in SEM. To avoid this, we use eq. (27) to eliminate |$\mathsf {w}_{\Gamma }^-$| from eq. (26), resulting in
This scheme does not require the split-node implementation and only keeps the wavefield on the positive side of the interface |$\mathsf {w}_{\Gamma }^+$|, which is easy to implement in SEM. Comparing with a typical wave simulation, only two additional steps are needed to incorporate interface discontinuity: (1) before applying the stiffness matrix in the elements that are on the negative side of the interface, an additional step is needed to add the displacement discontinuity |$\Delta \mathsf {w}_{\Gamma }$| for the points on the interface (Fig. 2), and (2) an additional force term |$\Delta \mathsf {T}_{\Gamma } - \mathsf {M}_{\Gamma \Gamma }^-\Delta \ddot{\mathsf {w}}_{\Gamma }$| needs to be included for points on the interface. The acceleration term |$\Delta \ddot{\mathsf {w}}_{\Gamma }$| can either be computed by the background solver as in the examples illustrated in Section 6, or by taking numerical derivatives of |$\Delta \mathsf {w}_{\Gamma }$|. The implementation of this scheme is outlined in Algorithm 1. Note that eq. (28) is equivalent to eq. (16) of Lyu et al. (2022), which is obtained based on the representation theorem.
Elements on the negative side of the interface (blue areas). Before applying the stiffness matrix in these elements, an additional step is needed to add the displacement discontinuity for the points on the interface. The injection interface |$\Gamma$| is marked by the thick black line, ‘+’ and ‘−’ mark the positive and negative sides of the interface.
An outline of proposed numerical scheme to solve the interface discontinuity problem using SEM. Steps required to incorporate the interface discontinuity in addition to the standard steps for SEM scheme are shown in blue color.
 |
|
An outline of proposed numerical scheme to solve the interface discontinuity problem using SEM. Steps required to incorporate the interface discontinuity in addition to the standard steps for SEM scheme are shown in blue color.
|
|
In principle, we can eliminate either |$\mathsf {w}_\Gamma ^-$| or |$\mathsf {w}_\Gamma ^+$| from eq. (26) using eq. (27). However, when using this algorithm to simulate wavefield injection, it is often desirable to implement an absorbing boundary, most commonly either the Stacey boundary condition Clayton & Engquist (1977) or PML Bérenger (1999); Komatitsch & Martin (2007), immediately outside the local domain to absorb the scattered wavefield. Eliminating |$\mathsf {w}_\Gamma ^-$| and keeping |$\mathsf {w}_\Gamma ^+$| allows fully independent implementation of absorbing boundary from wavefield injection. This is also noticed by Lyu et al. (2022; section 2.6.3 therein).
6 NUMERICAL EXAMPLES
We verify the correctness of our algorithm for the wave injection problem and illustrate the effect of absorbing boundary conditions using three numerical examples, at regional scale with and without local structural perturbations and at continental scales without local structural perturbations. For all examples, the waves in the local domain are simulated using the SPECFEM3D_Cartesian package (Komatitsch & Tromp 1999b). All numerical computations in this paper are performed on the Niagara cluster at the SciNet HPC Consortium (Ponce et al. 2019).
6.1 Injecting a plane wave into a regional domain without structural perturbation
In the first example, we simulate the injection of a plane wave into a Cartesian domain of |$x\in [0, 200\, \mathrm{km}], y\in [0, 200\, \mathrm{km}], z\in [-200\, \mathrm{km}, 0]$|. The background model is a two-layer 1-D model, including an infinitely-thick homogeneous mantle and a |$35\, \mathrm{km}$|-thick crust on top (Fig. 3a), with the properties of each layer given in Table 1. The model inside the local domain (Fig. 1a) is identical to the background model. The background wavefield, generated by the FK programme (Tong et al. 2014), includes an incoming P wave in the mantle from the west (azimuth = |$270^\circ$|) at an incident angle of |$15^\circ$| (Fig. 1a), and its reflections, transmissions and conversions at the Moho and the free surface. The source time function of the incoming P wave is a Gaussian function with dominant period of |$1\, \mathrm{s}$|.
Crustal and mantle properties of the two-layer 1-D model used for the numerical experiments in Section 6.1.
| Density (Kg m−3) | |$V_P$| (m s−1) | |$V_S$| (m s−1) | |
|---|---|---|---|
| Crust | 2600 | 5800 | 3198 |
| Mantle | 3380 | 8080 | 4485 |
| Density (Kg m−3) | |$V_P$| (m s−1) | |$V_S$| (m s−1) | |
|---|---|---|---|
| Crust | 2600 | 5800 | 3198 |
| Mantle | 3380 | 8080 | 4485 |
Crustal and mantle properties of the two-layer 1-D model used for the numerical experiments in Section 6.1.
| Density (Kg m−3) | |$V_P$| (m s−1) | |$V_S$| (m s−1) | |
|---|---|---|---|
| Crust | 2600 | 5800 | 3198 |
| Mantle | 3380 | 8080 | 4485 |
| Density (Kg m−3) | |$V_P$| (m s−1) | |$V_S$| (m s−1) | |
|---|---|---|---|
| Crust | 2600 | 5800 | 3198 |
| Mantle | 3380 | 8080 | 4485 |
The local domain is discretized into elements with size of |$2.5\, \mathrm{km}$| in each dimension, enabling accurate simulations of waves at |$\ge 1\, \mathrm{s}$| period. A four-element-thick PML is placed outside the injection interface to absorb scattered waves. 19 receivers are placed on the free surface at |$(x=i\times 10\, \mathrm{km}, y=100\, \mathrm{km}, z=0)$| for |$i=1,2,\ldots ,19$| (Fig. 1).
The Z and R components of the seismograms recorded at the receivers are shown in Fig. 4. All waveforms are lowpass filtered at |$1.5\, \mathrm{s}$|. Since there is no local structural perturbation, the total wavefield and the background wavefield should be identical. Fig. 4 shows that the difference between the waveforms computed by the injection algorithm (Algorithm 1) and the FK method are negligible, which proves that Algorithm 1 can accurately inject the background wavefield into the local domain.
The Vp model used in Sections (a) 6.1, without local structural perturbation and (b) 6.2, with local structural perturbations of |$\pm 30$| per cent at |$120\, \mathrm{km}$| depth, sliced at |$y=100\, \mathrm{km}$|. The blue triangles at |$z=0$| mark the positions of receivers. The black arrow in (a) shows the direction in which the incoming plane wave propagates.
(a) Vertical (Z) and (b) radial (R) component of the seismograms recorded at the receivers, with the direct P arrivals aligned at |$0\, \mathrm{s}$|. All the waveforms are lowpass filtered at |$1.5\, \mathrm{s}$|. The black lines show the seismograms computed using the injection algorithm, and the red dashed lines mark the difference between the waveforms computed using the injection algorithm (Algorithm 1 ) and the FK method, amplified by factor of 10. Since there is no local structural perturbation, the total wavefield and the background wavefield are identical. Therefore, the difference between the waveforms computed by injection and FK are negligible.
6.2 Injecting a plane wave into a regional domain with structural perturbations
We then add structural perturbations of |$\pm 30$| per cent in |$V_P$| and |$V_S$| in the local domain, as shown in Fig. 3(b). The setup of the numerical experiment is identical to Section 6.2. To verify the correctness of the injection result, we enlarge the local domain in x and z direction by factor of 2, that is, we inject the incoming P wave into an enlarged local domain of |$x\in [-100\,\mathrm{km}, 300\, \mathrm{km}], y\in [0, 200\, \mathrm{km}], z\in [-400\, \mathrm{km}, 0]$|. Since the structural perturbations are far away from the injection interface, we can safely assume that the results on the enlarged domain is accurate and can be treated as the ground truth. We also compare with the results using the method of Tong et al. (2014), which absorbs the scattered waves using the Stacey boundary condition, on the original local domain.
The waveforms at the receivers simulated by Algorithm 1 with PML, as well as a comparison between the accuracy of Algorithm 1 with PML and that of Tong et al. (2014), in terms of the difference between the waveforms simulated on the original and enlarged local domain, are displayed in Fig. 5. We can observe that using Algorithm 1 with PML, the discrepancy between the waveforms simulated on the original and enlarged local domain are negligible, indicating that Algorithm 1 with PML is very accurate. In comparison, using Tong et al. (2014)’s method on the original local domain, the discrepancy is apparent when amplified by factor of 10, especially at the receivers near the injection interface.
(a) Z and (b) R components of the seismograms recorded at the receivers when local structural perturbations have been added.
6.3 Injecting waves into a continental domain without structural perturbation
In the final numerical example, we test the capability of Algorithm 1 to accommodate a larger scale and a realistic wavefield excited by a tele-seismic earthquake. The local domain is a |$10^\circ \times 20^\circ$| chunk (Fig. 6a), truncated at |$1017.6\, \mathrm{km}$|, with the Earth’s spherical curvature incorporated. Both the structural model in the local domain and the background model are a smoothed version of AK135 (Fig. 6b), without attenuation, anisotropy or ellipticity. With |$104\times 200$| elements in the topmost layer, and a doubling layer immediately beneath the Moho (|$35\, \mathrm{km}$|) to accommodate the vertical velocity variation, the mesh enables accurate simulations of waves with period of |$\ge 6.5\, \mathrm{s}$|. To absorb the potential scattered waves, we implement a PML with thickness of four elements on the sides and two elements at the bottom. The mesh generation and the wave simulation are performed using the Cube2sph toolkit Liu et al. (2024), which is based on SPECFEM3D_Cartesian but generates meshes that accounts for the Earth’s spherical curvature and implements PML on curvilinear grids. The background wavefield is excited by an earthquake that is |$\sim 40^\circ$| away from the centre of the mesh (Fig. 6a), and is simulated by the AxiSEM package Nissen-Meyer et al. (2014). The source time function is a Gaussian function with dominant period of |$6.0\, \mathrm{s}$|. Wavefield injection is simulated using Algorithm 1 , with the injection interface co-located with the inner boundary of the PML. Fig. 7 shows the Z component of the seismograms recorded at the receivers as indicated on Fig. 6(a), computed by the AxiSEM package and the SPECFEM3D_Cartesian package with the proposed wavefield injection Algorithm 1 . The P, PP, S and Rayleigh phases match perfectly between the seismograms computed by AxiSEM and wavefield injection, suggesting that Algorithm 1 can be used for the injection of tele-seismic waves at a continental scale.
(a) The domain of simulation (thick black contour), source (red beachball) and receivers (blue triangles). (b) The mesh for computation. The elements are coloured by shear velocity, with warm and cold colours representing low and high shear velocities, respectively. A doubling layer is placed immediately beneath the Moho.
(a) Z component of seismograms computed using the AxiSEM package (red solid lines), the injection method (blue dashed lines), and their difference (green dashed lines), bandpass filtered at periods of |$8-50\, \mathrm{s}$|, recorded at the receivers shown in Fig. 6(a). The black dots mark the arrival times of P, PP and S phases. Panel (b) is a zoom-in of (a) on the body-wave part (black dotted contour). Note that the difference is very small in amplitude but is amplified by factor of 50 so that it can be seen visually.
7 DISCUSSION
7.1 Relationships with previously proposed injection methods
We discuss the connection between the method proposed in this paper based on interface discontinuity and previously proposed methods in three main categories: (1) domain reduction methods, which eliminate the terms that are related to the background wavefield outside the local domain from the global system of partial differential equations (PDE), (2) equivalent force methods, which replace the background wavefield with equivalent body forces at the injection interface, and (3) equivalent boundary condition methods, which transform the background wavefield into an equivalent boundary condition on the injection interface. Although these existing methods are derived from different approaches, we show in the following that all of them can be casted as solving the same interface discontinuity problem.
Moreover, we show that Algorithm 1 is numerically equivalent (i.e. having equivalent discretized forms) to the domain reduction method (Bielak et al. 2003) and the force mirror method (Masson et al. 2014). To clearly illustrate such equivalence, we rewrite eq. (28) in the context of wavefield injection. As discussed in Section 2, unlike the general interface discontinuity problem, for wavefield injection, interface |$\Gamma$| always encloses some local domain |$\Omega$|. Therefore, we can explicitly write |$\mathsf {w}_i:=[\mathsf {w}_l, \mathsf {w}_e]^T$|, with subscripts l and e representing points in and outside the local domain |$\Omega$|, respectively. With this, eq. (28) can be rewritten as
in which |$\mathsf {u}^0_{\Gamma }$| and |$\mathsf {T}^0_{\Gamma }$| represent the displacements and tractions of the background wavefield for points on the injection interface |$\Gamma$|.
7.1.1 Domain reduction methods
The original domain reduction method was formally proposed by Bielak et al. (2003), where the finite-element forms of the PDE system are written for the total wavefield in the global domain and the background wavefield outside the local domain, and the latter is eliminated from the former to obtain the algorithm to simulate wavefield injection. Leng et al. (2020) followed a similar approach, but performed the elimination upon the undiscretized weak form of the PDE system. As a result of the elimination, both methods introduce the scattered wavefield as the difference of the total and background wavefield outside the local domain. We perform the elimination and define the scattered wavefield (eq. 9) upon the original PDE, which leads to the introduction of the more general interface discontinuity problem. Nevertheless, no matter which form the elimination is performed upon, the original PDE, the undiscretized weak form or the discretized finite-element form, they all solve the same problem and therefore lead to equivalent algorithms. In fact, the weak form we derived (eq. 25) is identical to that in Leng et al. (2020, eq. 10 therein).
One minor difference is that the algorithms in Leng et al. (2020, eq. A3–A4 therein) and this paper (eq. 28) keep the scattered wavefield on the injection boundary, whereas Bielak et al. (2003) keeps the total wavefield. As discussed in Section 5, both treatments are equivalent, but keeping the scattered wavefield on the injection boundary leads to a more convenient implementation to incorporate the absorbing boundary.
Using the notations in this paper, and considering that the mass matrices are diagonal for SEM, the original domain reduction method (eq. 8 of Bielak et al. 2003) can be written as
in which |$\mathsf {u}^0_e$| represents displacements of the background wavefield for points outside the local domain. After a reorganization using |$\mathsf {w}_{\Gamma }^--\mathsf {w}_{\Gamma }^+=\mathsf {u}^0_{\Gamma }$|, eq. (30) becomes
The background wavefield satisfies the wave equation outside the local domain, resulting in a discretized form
Substituting eq. (32) into eq. (31) leads to eq. (29), which proves that the domain reduction method (Bielak et al. 2003) is numerically equivalent to Algorithm 1. Note that eq. (32) has a similar form with eq. (10) of Lyu et al. (2022).
7.1.2 Equivalent force methods
The equivalent force methods (Masson et al. 2014; Lin et al. 2019; Pienkowska et al. 2021; Lyu et al. 2022) seek to find body forces at the injection interface that are equivalent to the incoming wavefield. This is achieved by invoking the representation theorem (Lin et al. 2019; Pienkowska et al. 2021; Lyu et al. 2022), or force mirrors (Masson et al. 2014), which, as pointed out therein, are equivalent to the representation theorem when the injection interface is honoured in the local and background solvers. However, the representation theorem used in those studies are for boundary value problems (eq. 18), where the total wavefield in the local domain is represented by its own boundary values. Similarly, the force mirrors to reconstruct the total wavefield in Masson et al. (2014) are expressed in terms of the boundary values of the total wavefield itself. As discussed in Section 3, although it is not obvious theoretically that the total wavefield on the boundary can be replaced by the background field, it leads to correct algorithms because of the similarity between the representation theorems of the boundary value problem and the interface discontinuity problem. A similar statement is also made in Pienkowska et al. (2021) that using equivalent forces in terms of the background wavefield and the Green’s functions with local structural perturbations can reconstruct the total wavefield. Since the equivalent force methods can be seen as solving the interface discontinuity problem eqs (10)–(13) using its representation theorem eq. (17), methods in this category and our method are equivalent and must lead to the same results if the scattered waves are properly absorbed.
Selecting the window function to be 1 for points in the local domain |$\Omega$| and on the injection interface |$\Gamma$| and 0 for points outside the local domain |$\Omega$|, the force mirror method (eqs 27 and 28 of Masson et al. 2014) has the same form as eq. (30), and therefore, is numerically equivalent to Algorithm 1 and the domain reduction method based on the discussion in Section 7.1.1. However, as pointed out in Lyu et al. (2022), such selection of the window function requires a buffer layer outside the injection interface. To avoid this, Lyu et al. (2022) suggested to choose the window function to be 1 only for points in the local domain |$\Omega$| and 0 for points on the injection interface |$\Gamma$| and outside the local domain |$\Omega$| (following Lyu et al. (2022), we refer to such choice as the MYM window function). With this window function, the force mirror method (eqs 27 and 28 of Masson et al. 2014) can be written as
which can be reorganized into
Eq. (34) is identical to eq. (29) if and only if the background wavefield satisfies the discretized wave equation in the layer of elements immediately inside the injection interface
Unlike eq. (32), eq. (35) does not always hold, because the background wavefield does not necessarily satisfy the wave equation in the local domain. Nevertheless, in most practical cases, eq. (35) is true, because we typically leave a buffer zone between the injection interface and the local structures. In those cases, the force mirror method with the MYM window function is numerically equivalent with Algorithm 1.
7.1.3 Equivalent boundary condition methods
Enforcing the Stacey boundary condition (Clayton & Engquist 1977), which is a first-order paraxial approximation, on the scattered wavefield, an equivalent traction boundary condition can be derived (Monteiller et al. 2013, 2015; Tong et al. 2014, 2015). This can be seen as solving the interface discontinuity problem eqs (10)–(13) using the Stacey boundary condition at the injection interface. Since the traction boundary condition can be very easily incorporated in SEM, the equivalent boundary condition methods have simpler implementations than other methods, but suffer from the limitation that they cannot accommodate more accurate boundary conditions, such as PML, and may lead to inaccurate simulation results, as shown in Section 6.2.
7.2 Kinematic fault simulations
As mentioned in Section 2, the kinematic fault simulation problem is another special case of the interface discontinuity problem. Comparing with the widely used split-node scheme, the non-split-node scheme Algorithm 1 has advantages in its easy mesh generation and implementation. Here, we provide a numerical experiment to show that Algorithm 1 can indeed be used to perform kinematic fault simulations.
The mesh and structural model used in this numerical example are identical to those in Section 6.1. We define a |$80\, \mathrm{km}$|-long, |$15\, \mathrm{km}$|-deep vertical strike-slip fault along the Y-axis in the |$200\, \mathrm{km}\times 200\, \mathrm{km}\times 200\, \mathrm{km}$| simulation domain (Fig. 8a). Using Algorithm 1, we simulate the seismic wavefield excited by three prescribed Gaussian slip distributions
with a constant peak slip of |$s_0 = 3\, \mathrm{m}$|, different horizontal and vertical dimensions |$\sigma _y$| and |$\sigma _z$|, and a step-like source time function with a half-duration of |$\sigma _t = 3.76\, \mathrm{s}$|. The spatial slip distributions on the fault and the corresponding Y-component seismograms on the receivers aligning perpendicular to the fault plane with epicentral distances of |$0-90\, \mathrm{km}$| (Fig. 8a) are shown in Figs 8(b)–(d), and are compared with the results of point-source moment tensors with the same scalar moment computed by |$M_w = \int \mu s(y,z,\infty )\mathrm{dy}\mathrm{dz}$| (eq. 3.16 of Aki & Richards 2002), where |$\mu$| is the shear modulus, and the integral is numerically evaluated on the fault plane using Gauss–Lobatto–Legendre quadrature (eq. 30 of Komatitsch & Tromp 1999a). The waveforms are low-pass filtered at |$1.5\, \mathrm{s}$|.
The setup and results of the numerical experiment of performing kinematic fault simulations based on interface discontinuity. (a) The wave simulation is performed in a |$200\, \mathrm{km}\times 200\, \mathrm{km}\times 200\, \mathrm{km}$| domain (grey box), and an |$80\, \mathrm{km}$|-long, |$15\, \mathrm{km}$|-deep vertical strike-slip fault (red plane) is defined along the Y-axis, with the slip direction marked by the green arrows. The top panels of (b)–(d) show the spatial slip distributions, with the equivalent moment magnitude labelled on the upper-right corner, and the bottom panels show the corresponding Y-component seismograms at the receivers marked by the blue triangles in panel (a). We consider three Gaussian slip distributions (36), with spatial dimensions much smaller (b, |$\sigma _y = \sigma _z = 0.2\, \mathrm{km}$|), comparable (c, |$\sigma _y = \sigma _z = 2.5\, \mathrm{km}$|) and larger (d, |$\sigma _y = 15\, \mathrm{km}, \sigma _z = 2.5\, \mathrm{km}$|) than the smallest wavelength, and we compare the results based on interface discontinuity (red dashed lines) and point-source moment tensor (black lines).
In the first case (Fig. 8b), the dimensions of the slip distribution |$\sigma _y = \sigma _z = 0.2\, \mathrm{km}$| are much smaller than the minimum wavelength |$\sim 4.8\, \mathrm{km}$|, resulting in a moment magnitude of |$M_\mathrm{ w}=4.8$|. The waveforms computed based on the interface discontinuity problem using Algorithm 1 are identical to those using a point-source moment tensor. In the second case (Fig. 8c), |$\sigma _y = \sigma _z = 2.5\, \mathrm{km}$|, which are comparable with the minimum wavelength, resulting in |$M_\mathrm{ w}=6.1$|. Only the waveform at the receiver with epicentral distance of |$10\, \mathrm{km}$| is different using Algorithm 1 and the point-source moment tensor, due to the finite-fault effect. In the third case (Fig. 8d), |$\sigma _y = 15\, \mathrm{km}$|, which is larger than the minimum wavelength, and |$\sigma _z = 2.5\, \mathrm{km}$|, resulting in |$M_\mathrm{ w}=6.8$|. The difference between the results using Algorithm 1 and the point-source moment tensor can be clearly observed at all the receivers, as far as |$90\, \mathrm{km}$| from the fault plane.
This numerical experiment shows that when using Algorithm 1 to perform kinematic fault simulations, the results match well with the point-source moment tensor when the dimensions of the slip distribution is much smaller than the minimum wavelength, and can capture the finite-fault effect when the slip patch is comparable or larger than the minimum wavelength. This indicates that the non-split-node numerical scheme for the interface discontinuity problem Algorithm 1 can effectively perform kinematic fault simulations. It remains to be further investigated whether a similar non-split-node scheme can be implemented to perform dynamic fault simulations.
8 CONCLUSION
We formally define the interface discontinuity problem of the seismic wave equation, discuss its representation theorem and uniqueness, and design an efficient SEM scheme to solve it. We study the wavefield injection problem under the framework of the interface discontinuity problem, as the former can be seen as a special case of the latter. We draw connections with the previously proposed methods for wavefield injection based on domain reduction, equivalent forces and equivalent boundary conditions by recognizing that they all solve the same interface discontinuity problem. We also show that the SEM scheme proposed in this paper is numerically equivalent to the domain reduction method and the force mirror method.
Furthermore, we clarify that wavefield injection should be seen as an interface discontinuity problem, rather than a boundary value problem. As required to uniquely solve the interface discontinuity problem, we state that the displacement and traction of the background wavefield must be simultaneously supplied to simulate wavefield injection. However, due to the similarity between the representation theorems of the boundary value problem and the interface discontinuity problem, previously proposed methods citing the representation theorem or the force mirrors of the boundary value problem are still correct.
Based on numerical experiments, we verify that the SEM scheme for the interface discontinuity problem can accurately simulate wavefield injection at both regional and global scales. With background wavefield generated by some global wave simulator, the proposed SEM scheme can be used in box tomography based on tele-seismic scattered waves targeted at various structures. The proposed scheme can be also applied to investigating the effect of local structures on ground motions and modelling earthquake hazards. Finally, we show that kinematic fault simulations can be easily performed based on SEM using this scheme if the fault plane can be honoured by the mesh.
ACKNOWLEDGMENTS
The authors acknowledge that open-source visualization software Paraview and GMT, and Python packages Numpy, SciPy, Matplotlib and Obspy are heavily used to produce the figures in this paper. TL was funded by the Ontario Trillium Scholarship. We thank Barbara Romanowicz, Chao Lyu and an anonymous reviewer for their valuable comments and suggestions. This study was part of the TREMOR project supported by Air Force Research Laboratory contract FA9453-24-9-0001 and funded through the Geophysical Detection of Nuclear Proliferation University Affiliated Research Center at University of Alaska Fairbanks. TL, ND, TL and QL are also partially supported by NSERC Discovery grant #487237. The high-performance computing performed for this research was enabled by support provided by the Digital Research Alliance of Canada (alliance.can.ca).
DATA AVAILABILITY
The computer programs for the implementations of Algorithm 1 in SPECFEM3D_Cartesian has been merged into the SPECFEM3D Github repository (https://github.com/SPECFEM/specfem3d/pull/1706). The generation of background wavefield database from AxiSEM is handled by the AxiSEMLib available at https://github.com/nqdu/AxiSEMLib.
