Element-by-element parallel spectral-element methods for 3-D 1 acoustic-wave-equation-based teleseismic wave modeling 2

11 The increasing demand for the high-resolution imaging of deep lithosphere structures 12 requires the utilization of a teleseismic dataset for waveform inversion. The construction of an 13 efficient algorithm for the teleseismic wavefield modeling is valuable for the calculation of 14 misfit kernels or Fréchet derivatives when the teleseismic waveform is used for adjoint 15 tomography. Here, we introduce an element-by-element parallel spectral-element method 16 (EBE-SEM) for the efficient modeling of teleseismic wavefield propagation in a localized 17 geology model. Under the assumption of the plane wave, the frequency-wavenumber (FK) 18 technique is implemented to compute the boundary wavefield used for constructing the 19 boundary condition of the teleseismic wave incidence. To reduce the memory required for the 20 storage of the boundary wavefield for the incidence boundary condition, an economical 21 strategy is introduced to store the boundary wavefield on the model boundary. The perfectly 22 matched layers absorbing boundary condition (PML ABC) is formulated by the EBE-SEM to 23 Solid Earth Discuss., doi:10.5194/se-2017-38, 2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c © Author(s) 2017. CC-BY 3.0 License.

absorb the scattered wavefield from the model interior.The misfit kernel (derivatives of the waveform misfit with respect to model parameters) can be easily constructed without extra computational effort for the calculation of the element stiffness matrix per time step during the calculation of the adjoint wavefield.Three synthetic examples demonstrate the validity of EBE-SEM for use in teleseismic wavefield modeling and the misfit kernel calculation.

Introduction
Teleseismic waves provide tremendous amount of information for the detection of crust and upper mantle structures (Rondenay, 2009).In the past fourth years, many techniques have been established to analyze teleseismic wave datasets, including receiver function analysis (Langston, 1977, Kind et al., 2012), teleseismic wave travel-time tomography based on ray theory (Zhang et al., 2011), teleseismic migration (Shragge et al., 2006), and teleseismic scattering tomography (Roecker et al., 2010;Tong et al., 2014a).For the requirements of high-resolution teleseismic wave imaging, the adjoint-state method has become the state-of-the-art technique for teleseismic wave imaging (Tong et al., 2014a;Monteiller et al., 2015).
The technique of adjoint tomography constructs the Fréchet derivatives of the objective function with respect to model parameters by two times of numerical solving full seismic wave equation (Tromp et al., 2005;Liu and Tromp, 2006).Adjoint tomography has been successfully implemented to investigate crust (Tape et al., 2009) and continental subsurface heterogeneity (Chen et al., 2015).Indeed, the seismic wave-equation-based adjoint waveform tomography, which has a greater resolution than the ray-based traveltime tomography for the same dense seismic ray coverage (Liu and Gu, 2012), is able to image the small-scale (half of the minimum wavelength) heterogeneity (Virieux and Operato, 2009).The main drawback of Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License. the adjoint tomography is its huge computational burden.The computational requirement is linearly related to the earthquake events used for tomography inversion and the iterations required by the optimization technique.For a typical local scale model, such as southern California region, to perform adjoint tomography inversion, several thousand of 3-D full-wavefield simulations are required (Tape et al., 2007(Tape et al., , 2009)).
Because most earthquakes occur in the crust and seismic rays cannot easily illuminate the deep lithosphere in local region seismic tomography, it can be difficult to image the deep lithosphere structures (Tong et al., 2014b).Increasing the model size enables more deep reflections and refractions to be included in the inversion dataset; as a result, deep structures can be inversed by fitting these reflected and refracted waveforms (Chen et al., 2015).
However, for continental-scale models, it is difficult to invert short-period seismic data on a standard computing cluster, such as 1-2 s for P waves and 3-6 s for S waves (Tong et al., 2015).
To reduce the amount of computation involved in solving the full seismic wave equation, many hybrid methods have been developed to localize the 2-D/3-D numerical solvers.The boundary conditions for the localized simulation model are provided by fast-computing 1-D analytical solutions for the 1-D background Earth model (Capdeville et al., 2003;Monteiller et al., 2013Monteiller et al., , 2015;;Tong et al., 2014aTong et al., , 2015)).The 2-D/3-D responses to the heterogeneity inside the localized model contribute to the coda waves of the teleseismic phases, and the 2-D/3-D effects outside the model are neglected.This assumption of a 1-D background layered Earth model is similar to that in receiver function analysis and is often effective for a station that is sufficiently far from the source (Langston, 1977;Rondenay, 2009).
Although computation efforts can be efficiently reduced by these hybrid methods, the Solid Earth Discuss., doi:10.5194/se-2017Discuss., doi:10.5194/se- -38, 2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.computational costs are still excessive for a small workstation when we are faced with the several thousand of forward simulations required in 3-D teleseismic adjoint tomography.To further reduce the computational costs, Roecker et al., (2010) constructed a frequency domain 2.5-D hybrid method for teleseismic wave simulations.To simplify the teleseismic wavefield invariant along a particular axis, the 2.5-D formulation can significantly reduce the computational demands.However, the 2.5-D formulation may restrict the application of the method in an arbitrarily heterogeneous media (Tong et al., 2015).
In addition to the great computational demand (CPU time) associated with the 3-D hybrid methods, satisfying the memory requirements for storing the boundary wavefields to construct teleseismic incident boundary conditions is another important issue that should be carefully considered.Tong et al., (2015) adapted the Clayton and Engquist-type (CE-type) boundary condition (1977) to interface the 1-D background analytical solution with a numerical solver on the boundary of a localized model.This treatment can not only assign the teleseismic incident condition for the computational domain, but also absorb the scatter wavefield from the interior of the heterogeneous model.The implementation of the CE-type boundary condition is extremely simple and dese not substantially increase the required CPU time.
However, the CE-type boundary condition can efficiently absorb the waves only at approximately normal incident angles, and the incident waves at grazing angles may be reflected back to the model (Yang et al., 2003), which may reduce the accuracy of the forward and adjoint wavefields in teleseismic adjoint tomography and thus decrease the accuracy of the constructed Fréchet derivatives.Note that the CE-type boundary condition requires 9 wavefield components (3 displacement components and 6 stress components) to be stored on the model boundary; this requirement may be a great burden to the computer memory for the Here, we introduce EBE-SEM for the efficient modeling of teleseismic plane-wave propagation in localized models.One great advantage of EBE-SEM is the easy parallelization of the spectral-element algorithm, which does not require assembly of the global stiffness matrix.The spectral elements are equally allocated to every CPU core, and the product of the stiffness matrix with the solution vector is calculated element by element; these aspects are quite useful to ensure load balance among the CPU cores.In addition, the misfit kernel can be efficiently constructed, because each element stiffness matrix is stored in CPU memory and it is simultaneously available when calculating the misfit kernel.The perfectly matched layers (Collino and Tsogka, 2001;Komatitsch and Tromp, 2003;Liu et al., 2014)

EBE-SEM
A schematic of a teleseismic plane wave entering a localized model is depicted in Figure 1, where the localized mode is delineated by the blue lines.We first introduce EBE-SEM for acoustic wave propagation in an infinite half space, which includes the localized model.We denote the total wavefield t u , which is the summation of the background wavefield FK u in layered media and the scattered wavefield s u .FK u can be efficiently calculated with the wavenumber-frequency domain method (FK) (Haskell, 1953;Zhu and Rivera, 2002;Tong et al., 2014a

EBE-SEM for total wave simulation
We assume that the total wavefield (P or S wave) obeys the following acoustic wave equation (Tong et al., 2014c): where the two dots above t u denote the second-order time derivative, c is either the P or S wave velocity model, and ( ) is the gradient operator.Following the classical spectral-element method (SEM) (Seriani, 1997;Komatitsch and Vilotte, 1998;Komatitsch and Tromp, 2002), using an arbitrary test function v to multiply both sides of Eq.
(1) and integrating over the domain Ω , we obtain the following weak form wave equation: where n is the normal vector of the boundary ∂Ω .Under the natural boundary condition, the right side of Eq. ( 2) is zero.If the infinite space is decomposed into N non-overlapping hexahedral elements, then Eq. ( 2) can be written as: Each hexahedral element is mapped to the master cube [ ] 3 1, 1 − , and we choose the interpolation points to be consistent with the Gauss-Legendre-Lobatto (GLL) quadrature points (Cohen, 2002).The interpolation basis functions in each element is as follows: where the subscript of φ denotes the i-th basis function; ξ , γ , and η are the three coordinates in the isoparametric coordinate system; and ϕ is the Lagrange polynomial.The i-th interpolation point in the physical element is mapped to the ( ) where n is the polynomial order of the interpolation basis.Substituting Eq. ( 5) into Eq.(2), and using the GLL quadrature rule, we obtain the ordinary differential equation in terms of time: ( ) Eq. ( 6), we have .
N e e e t t e T T An explicit temporal scheme, such as the Newmark scheme (Liu et al., 2017a), and structure-preserving schemes (Liu et al., 2017b) can be adopted for time integration.Here, we simply use the second-order central difference method (Dablain, 1986) for the time evolution.
The product of the stiffness matrix and the solution vector is performed element by element at the element level, rather than at the global level.The computational burden of SEM mainly stems from the matrix-vector product, and the significant advantage of the element-by-element scheme of Eq. ( 9) is that the great computational balance among CPU cores when the parallel algorithm is performed on a workstation.Our parallel code is based on Message Passing Interface (MPI) library (www.mpich.org/downloads),and we equally distribute spectral-elements to each of the CPU cores.Because e K is symmetric, only the upper triangle of the e K must be stored in computer memory.

EBE-SEM for PML formula
In our previous work, a second-order PML absorbing boundary condition (PML ABC) are formulated by the mixed-grid finite-element method (Liu et al., 2014).Here, we use this type of PML ABC to absorb the scattered wavefield s u .The formula of the PML ABC can be written as Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.Based on the natural boundary condition, the corresponding wake form of Eq. ( 10) The discretized version of Eq. ( 11) is .
From Eq. ( 14), we observe that element damping matrices do not need to be stored, and only the mass matrix and the element stiffness matrix must be kept in computer memory.Note that because e x K , e y K , and e z K are not symmetric, and all the elements of each element matrix should be stored.

A discussion of EBE-SEM
Although EBE-SEM is specially designed for teleseismic wave modeling, i.e., Eq. ( 9) is for teleseismic total wavefield propagation if the proper teleseismic incident boundary condition is added and Eq. ( 14) is for absorbing the scatted wavefield, EBE-SEM can be directly used for wavefield simulation for an earthquake that occurred interior of the model (computation domain) if a source term is added in Eq. ( 9).
The seminal work for introducing EBE-SEM was that of Seriani (1997) large spare linear system of equations, which may be not efficient for larger-scale seismic waveform modeling.Another aspect of the BEB algorithm designed by Seriani (1997) is that this algorithm is for use with only rectangular elements; the restriction may be problematic for curved elements.In contrast, our EBE-SEM can be used for general cases.
If appropriate boundary conditions are added for the computation domain and the PML domain, then the computation will be isolated in the counterpart domain, i.e., Eq. ( 9) for the computation domain and Eq. ( 14) for the PML domain.The discussion of the boundary condition is presented below.

Teleseismic wave incident boundary conditions
For the completeness of the paper, we first simply introduce the FK method for determining the 3-D acoustic-wave-equation-based plane-wave propagation in layered media.
For more details, the reader is referred to Haskell (1953) and Tong et al. (2014a).Our focus is to construct the teleseismic wave incident boundary conditions and develop a highly efficient method for storage of the boundary wavefields.

Plane wave propagation in 1-D layered media
Transforming the acoustic wave equation into the ( ) .
Eq. ( 16) is changed to the following ordinary differential equations: where 1 . By calculating the eigenvalue and eigenvector of F , the general solution of Eq. ( 17) can be written as where 1 C and 2 C correspondto the amplitudes of the down-going and up-going acoustic waves, respectively, and is the ray parameter.As shown in Figure 2, the velocity and thickness of the m-th layer are denoted as m c and m h , respectively.We define the following matrices for the m-th layer: The relationship between the wavefields at .
n n a a a a 2 y at the free surface can be written as Based on the free surface boundary condition 2 0 y = , we have Substituting Eq. ( 25) and ( 23) into Eq.( 18), we obtain the wavefield at m z z = .Based on Eq.
(21), we can calculate the wavefield in each layer.

Incident boundary conditions
We assume the infinite space is composed of the hexahedron set , , , , , , .
B e e e e The first 1 N elements of B compose the computation domain (the domain bounded by the blue lines in Figure 1); the elements from 1 1 N + to 1 2 N N + compose the PML domain.We define the set of elements in the computational domain as A total of 3 N elements in C B that contact the boundary of PML domain are collected in the set: The set of 2 N elements that compose the PML domain is given by 1 , , .We define an operator that maps the elements of vector A to a target vector t A according to the numbering set n A .The index number of each elements of set t A in set A is gathered From Eq. ( 9), we have .
where C A is the numbering set of the nodes in the computational domain.Eq. ( 32) shows only matrix-vector products in the elements of the computational domain and the elements of the PML domain in contact with boundary of the computational domain contribute to the acceleration wavefield on the nodes in the computational domain.From Eq. ( 32), we further have is the numbering set of nodes of the computational domain located on the interface between the computational and PML domains.The second term of the right-side of Eq. ( 33) can be written as Eq. ( 33) is the plane-wave incident condition of the computational domain.From Eq. ( 34), we need to store only the boundary wavefield extremely important for large-scale models.This storage technique is similar to the method of linear combination of a boundary wavefield, which was proposed in our previous work (Liu et al., 2015).
To simplify the discussion, the boundary condition of the first equation of Eq. ( 14) is discussed.For Eq. ( 14), we have

U
does not need to be calculated by Eq. ( 14).From Eqs. ( 35) and ( 36 (37) In addition to the boundary condition Eq. ( 36), the Dirichlet boundary condition is added on the outer boundaries of the PML domain, and the natural boundary condition is added on the planes that connecting the free surface of the computation domain.Because of the boundary conditions, the element-level matrix-vector product is restricted to only the element in the computational domain and the PML domain.Because boundary conditions Eq. ( 33) and ( 36) involve the plane wavefields, and we call Eqs. ( 33) and ( 36) teleseismic wave incident conditions.

Analysis of the computational costs
We use the model in Figure 1 to quantitatively discuss the computational cost of EBE-SEM for teleseismic wave modeling.Because the thickness of the PML domain is only three elements wide in our numerical examples, the number of float point operatorations is trivial compared with the computational domain, and the main computational cost of the PML domain is mainly from the storage requirement of the boundary condition (Eq.( 36)).
The model is decomposed into a total of 75000 cubic elements with a size of 2 2 2 km km km × × .The fifth-order interpolation polynomial is used in the space.The time interval is 0.01 s, and the total time length for numerical modeling is 120 s.
The float point operation in the computation domain is mainly from the element matrix-vector product, and a total of This increased compuational amount may be compensated for the great load balance of EBE-SEM in parallel computing, as will be discussed in the numerical examples.
The memory requirements of EBE-SEM and the conventional SEM for teleseismic wave modeling are presented in Table 1.A total of 13.64 GB is required to store the boundary wavefield to construct the teleseismic incident conditions.If the classical compressed sparse row (CSR) storage format (Greathouse and Daga, 2014) is used for the storage of spare stiffness matrix of the conventional SEM, then the storage requirement is 24.04 GB, which is nearly four times as large as the storage requirement of EBE-SEM to store the element stiffness matrices.Two factors contribute to this storage difference.First, CSR must store the row and column information of the global stiffness matrix in addition to the storage of non-zero elements.Second, the element stiffness matrix of EBE-SEM is symmetric, i.e., only the upper triangle of element stiffness matrix must be stored.

Numerical examples
Three numerical examples are provided to validate the efficiency of EBE-SEM for teleseismic wave modeling.In all the examples, the Gaussian source-time function with a cut-off frequency of 1 Hz is used (Tong et al., 2014a).

Benchmark for the 1-D crust-upper mantle model
Except for parameters presented in Section 4, the upper velocity and the lower layers are 3000 m/s and 4500 m/s, respectively; the incidence angle is 15 o , and the azimuth angle is 170 o .
The cross section of the plane-wave front with Moho is (-300 km, 0, 30 km) at the initial time t=0.A seismic station is located at (70 km, 50 km, 0 km).The computed results are shown in  To qualitatively evaluate the accuracy of EBE-SEM for teleseismic wave modeling, the synthetic seismology at the station is compared with the reference solution, which is generated by FK.The result is shown in Figure 4.As depicted in Figure 4(a), excellent agreement is achieved between the synthetic and reference waveforms.The direct wave with largest amplitude is followed by crust multiples with relatively small amplitude.Although the amplitude of the crust multiples is small, the phases, amplitudes, and travel times of these multiples are correctly modeled.The error curve is shown in Figure 4(b).The maximum difference between the numerical solution and the reference solution is approximately 5 7 10 − × , which is extremely small compared with the amplitudes of the direct wave and crust multiples.
The computation was performed on a workstation with 2 Intel Xeon CPUs (E5-2680 v3) and 128 GB of RAM, and the total number of CPU cores is 24.The master-slave communication pattern is used in our parallel algorithm, which is similar to that of Komatitsch and Tromp (2002).The parallel efficiency is shown in

Teleseismic waveform misfit kernel
One key advantage of the EBE-SEM is the great convenience of constructing the misfit ( ) where C Ω denotes the computational domain.The i-th element of the misfit kernel can be calculated by the following equation: where NT is the total time step.Eq. ( 41) shows that element stiffness matrices are used for the construction of the misfit kernel.
We consider the model in the second numerical example to be the real model, and the 1-D layered model in Figure 1 is the initial domain.The observed and synthetic waveforms at the station (red triangle in Figure 1) are shown in Figure 7   wave.The incident angle is θ , which is the angle between the normal vector and -z.The azimuth angle is ϕ , which is the angle between the projection of normal vector into x-y plane and -y.The azimuth angle is not explicitly shown in this figure.
Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.case of a relatively large scale model decomposed by a dense numerical mesh.
are formulated by EBE-SEM to effectively absorb scattered waves.A detailed discussion is presented to incorporate the teleseismic incident boundary condition for EBE-SEM, and only two layers of the boundary wavefield must be stored in computer memory.The high efficiency of EBE-SEM for teleseismic wave modeling and misfit kernel construction is demonstrated by three numerical examples.

KM
are presented in Appendix A. Because the interpolation points are consistent with the GLL quadrature points, e M is diagonal; this is one great advantageous of Legendre or polynomial-based SEM because explicit inverse of the mass matrix can be easily obtained.The assembled global diagonal mass is as follows:  is used to distinguish the difference between the back-projected mass matrix and the element mass matrix in Eq. (6).Because the elements of the element mass matrix corresponding to the common points shared by the adjacent elements are Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.summed together to form the global mass matrix in Eq. (7), some elements of e M  are greater than those of e M .If we denote the global solution vector as t U , then the projected back element solution z d are damping coefficients along the three coordinate axes; the 180 subscript of d denotes the first spatial derivative of the damping coefficient;

D
are the element damping matrices.The detailed expressions of these 189 element matrices are presented in Appendix A. The damping matrices are also diagonal 190 matrices, whose element values are scaled by a constant compared with the element values of 191 e M .If we denote the damping matrices as x D , y D , and z D , then we obtain the 192 element-by-element scheme for Eq.(12): 193
. In the 2-D case, the element-by-element scheme is combined with the Chebyshev orthogonal polynomial-based SEM.Because the Chebyshev orthogonal polynomial is orthogonal associated with the weight 2 1 / 1 ξ − , the Chebyshev orthogonal polynomial-based SEM cannot lead to a diagonal mass matrix and an iterative algorithm is generally used to solve a Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.
we consider that the incident plane wave has only an up-going component and has unit amplitude, then we have

N
elements in P B that contact the boundary of computation domain are gathered in the set Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.
of nodes located on the interface between the computational and PML domains; and it is not necessary to store FK U on the nodes of the elements in , are required in each time step.If a global stiffness matrix is assembled, then the product of the global stiffness matrix with the global solution vector requires a total of 9 6.435220951 10 × float point operations.The computational burden of EBE-SEM is larger than that of the conventional SEM because the acceleration on common nodes shared by the adjacent elements is calculated Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.morethan once.However, the number of operations EBE-SEM is increased by only 8.5%.

Figure 3
Figure 3 shows snapshots at three time instants.When t=10 s, the plane waves are still Figure 5.Because the master-slave communication pattern requires at least two CPU cores to perform the parallel algorithm, the computation time of single CPU core is not shown in Figure 5.When the number of CPU cores is 2, the communication amount between the CPU cores is negligible Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.compared with the computation amount.If we denote the CPU time is as 2 T when 2 CPU cores are used in parallel computation, the CPU time of n CPU cores can be estimated by line in Figure 5 is the estimated CPU times (theoretical times), and the blue line with circles is the practical CPU times.It can be clearly observed that the actual CPU times are extremely close to the theoretical times, even if the number of CPU cores is 24.The red arrow in Figure 5 indicates the abnormal CPU time.This phenomenon is attributed to the excessive communication amount when the number of CPU cores is 9. However the anomaly is not large compared with the neighboring CPU times.5.2Plane-wave incidence to a 3-D modelTo demonstrate that EBE-SEM works well in 3-D heterogeneous models, an abnormal structure with cube shape of an additional 15% plus wave speed is added at the center of model in Figure1.The size of the abnormal structure is velocity structure, the other computational parameters are the same as those of the first numerical example.The simulation results are shown in Figure6.The red arrows in the upper plots of Figure 6 indicate the distortion of the wave front because of the velocity anomaly of the media.Because of the plus velocity anomaly, the distorted wave front travels faster than the undistorted wave.From the lower plots of Figure 6, we can observe strong scatted waves.As the yellow arrows indicates, the strong scattered waves do not reflected into the interior of the model because of the efficiency of the PML ABC constructed in the paper.
Figure 8.Because of the singularity of the source, large amplitudes are distributed in the average value of the damping coefficient on the element e.

Figure 1 .
Figure 1.Schematic of a teleseismic plane wave entering a localized model, which is a

Solid
Figure 2. 1-D layered model.The thickness, velocity, and ray parameter of each layer are 579

Figure 3 .
Figure 3. Plane wave incident on a 1-D crust-upper mantle model.The snapshots at the upper,

Figure 4 .Figure 5 .
Figure 4. (a) Comparison of the waveforms generated by EBE-SEM (solid line) and FK

Figure 6 .
Figure 6.Snapshots in a heterogeneous model.Snapshots in the upper and lower panels are

Figure 7 .Figure 8 .
Figure 7.The construction adjoint source.(a) The synthetic waveform from the second ). s u is excited by the incident plane wave because of the heterogeneity of the node in the Solid Earth Discuss., doi:10.5194/se-2017-38,2017 Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.
kernel, because the element stiffness matrix is kept in computer memory and no extra effort is Solid Earth Discuss., doi:10.5194/se-2017-38,2017Manuscript under review for journal Solid Earth Discussion started: 8 May 2017 c Author(s) 2017.CC-BY 3.0 License.require to recalculate the element stiffness matrix.To illustrate this advantage, we first define x is the teleseismic incident parameters, d is the observed dataset, and s is the synthetic dataset.The cross section between the plane wave and the Moho at y=0, the incident angle, and the azimuth angle constitute the teleseismic incident parameters.To simplify the discussion, only one source and one station are considered in this example.Based on the continuous adjoint method discussed by Fichtner s