An improved nearly-analytic exponential time difference method for solving seismic wave equations

© 2025 by the Authors. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution -Noncommercial 4.0 International License (CC-by the license) ( https://creativecommons.org/licenses/by-nc/4.0/ )

Download PDF

Cite

XML

Abstract

Yun, N. and Sim, C., 2020. An improved nearly-analytic exponential time difference method for solving seismic wave equations. Journal of Seismic Exploration, 29: 99-124. We present an improved nearly-analytic exponential time difference (INETD) method, which is a version of the nearly-analytic exponential time difference (NETD) method proposed by Zhang. The NETD method is based on a combination of the exponential time difference (ETD) scheme with the nearly-analytic discrete (NAD) operator. In solving the acoustic and elastic wave equations, NETD first uses the NAD operator to approximate the high-order spatial differential operators, and then the resulting semi-discrete ODEs are solved by the ETD method. Because NETD is based on the structure of the Lie group method, it can achieve more accurate results than the classical methods. The main purpose of the present article is to improve the stability and the time accuracy in the original NETD method. To this end, we first derive an implicit ETD scheme and apply it to approximate the time derivatives instead of the explicit scheme in the original NETD method. For this new method, we analyze in detail its stability condition, the dispersion relation and the theoretical and numerical errors, and compare it with the NETD method. In addition, for the 2D case, using INETD, NETD and the fourth-order Lax-Wendroff correction (LWC) methods we simulate acoustic and elastic wave-fields, and compare numerical results for these different methods. Our theoretical analyses show that the stability condition of the INETD method is more relaxed than that of the original NETD method. The temporal accuracy of the INETD method is increased from second order in the original NETD method to third order, and spatial accuracy remains the same as that of the original. Moreover, the INETD method has the same computational costs and storage space as the NETD method and the numerical dispersion is less than the NETD method. The results of numerical experiments demonstrate the high efficiency of this method for acoustic and elastic wave modeling.

Keywords

nearly-analytic discrete method

exponential time difference method

Lie group method

numerical dispersion

numerical approximations and analysis

References

2000; Robert, 1991) and so on. Geometric integrations are shown to be
superior to the general-purpose methods in some aspects, especially in long
time computation, owing to their preservation of the qualitative properties.
Thus, combining of geometric integrations with NAD may lead to decreases
in numerical dispersion and numerical error over other classical methods.
Motivated by this fact, Zhang et al. (2014) developed a nearly-analytic
exponential time difference (NETD) method. NETD uses an exponential
time difference (ETD) scheme to approximate the time derivatives and the
nearly-analytic discrete (NAD) operator to approximate the space
derivatives. The ETD scheme used in this method can be also derived from
the idea of the Lie group method. Therefore, the NETD method by
combining the Lie group method with the NAD method to solve seismic
wave equations effectively suppresses the numerical dispersion and
preserves the qualitative properties of PDEs. Through theoretical and
numerical analysis, they showed that NETD method provides more accurate
results, has high computational efficiency, and a low memory requirement
compared to the Lax-Wendroff correction (LWC) and the staggered-grid
(SG) methods. However, the NETD method has some disadvantages; for
example, its stability criterion is relatively compact and has only
second-order time accuracy. It is important to improve the stability of the
NETD method and increase its time accuracy for practical calculations.
From this view point, we focus on this problem. We present an
improved NETD (INETD), for solving acoustic- and elastic-wave equations.
The INETD scheme is composed of the following two steps: (1) Derivation a
new implicit ETD scheme, (ii) Combination of the implicit ETD for the
temporal discretization with the NAD method for the spatial discretization.
This produces an improvement in the stability and the time accuracy over
the NETD, because the explicit ETD scheme is used to approximate the time
derivatives in the NETD. In general, the implicit scheme has a good stability
property, but it requires more expensive computations than the explicit
scheme, because it solves the system of equations at every temporal step.
We show that the INETD method needs the same computational costs and
storage space as the NETD method. On the other hand, we find another
possibility for improving the stability of the INETD method. It is to use the
terms up to the third order in the truncated power series expansion used to
explicitly solve the INETD system of linear equations. Our analysis shows it
brings about a remarkable improvement in the stability of the method. On
the whole, the our INETD method are shown to be superior to the FDM, the
PS and the NETD methods in some aspects such as suppressing numerical
dispersions, high accuracy, computational efficiency, a low memory and the
improvement of stability.
This article is organized as follows. First, we review the basic NETD
scheme for solving the acoustic wave equation in a 2D homogeneous
medium. Then, we derive the INETD scheme. Next, we investigate in detail
the properties of the INETD method. We derive the stability condition for
the 1D and 2D cases. We estimate the theoretical and numerical errors of the
INETD method, and compare against the NETD method. And then, we
derive the dispersion relation and give the numerical dispersion results.
Finally, for the 2D case we use the INETD, the original NETD and the
fourth order LWC methods to simulate acoustic wave fields, and present
some comparisons of numerical results for different methods. And then a
homogeneous elastic model and a two-layer elastic model are further
selected to investigate the computational validity of elastic wave-field
simulations.
Our theoretical and numerical results show that the INETD is an
improvement over the original NETD in stability, numerical error and
numerical dispersion.
BASIC NEARLY-ANALYTIC EXPONENTIAL TIME DIFFERENCE
(NETD) METHOD
In this section, we review and analyze the derivation of the original
NETD method (Zhang et al., 2014) for solving the acoustic wave equation.
Exponential time difference method
Consider the stiff system
du
—=Lu+N(u,t
alt (ut), (1)
where L is a higher-order linear term and N is a lower-order nonlinear
term. To solve eq. (1), the equation is multiplied by an integrating factor
< (exp(-Ltu) = exp(-La)N(u, ¢) , (2)
and then eq. (2) is integrated over a single temporal step of length 1,
u'*! =exp(Lr)u' + exp(L7) fexp(-Ls)N(u(, +5),t, +s)ds, (3)
where 4” denote the numerical solution on the n-th time layer ¢,,. To
approximate the integration in eq. (3), using the approximation as
N(u(t,+5),t,+5)=N(u'Qt,) , (4)
the numerical scheme for solving eq. (1) is obtained as:
u' = exp(Lr)u' +L! (exp(Lr)-I)N(u', t, ) . (5)
Basic computational formulae
The acoustic wave equation in 2D homogeneous medium is written as
du
we co Au, (6)
Ce
where z and c, are the displacement and the acoustic velocity, respectively,
2 2
oe, . .
and A= Pye + ra is the Laplace operator. The same notations as that in
x az
the basic NETD (Zhang et al., 2014) are used in our present study, i.e.,
7 7
ou du 07/ dw dw
arya 0 人 “有二 | > we (me .
t
Eq. (6) can be rewritten as
ou -WwW
at
aW _ aU aU)’ (7)
ot 0 ax? az?
The high-order spatial derivatives in eq. (7) are determined by the
displacement U and its gradients using the fourth-order NAD operator as
shown in Appendix A of Zhang et al. (2014). Then eq. (7) can be
transformed into an ODE system as follows:
d(U) fo 1(U) <
| ae @
The first part of the right side of eq. (8) is a linear combination of u,w
with its gradients at the grid point (x;, z,), and the second part is a linear
combination of u,w with its gradients at its neighboring grid points. I isa
3x3 identity matrix and 2 is a 3x3 matrix that depends on the spatial
discrete operator. Using the notations
7 - (0 I
Y,,=(U,W)',, and 人
and regarding N(Y,, ) as a nonlinear function of Y,,, eq. (8) can be
rewritten as
qY,, - ~
的 + N(Y, 外， (9)
Eq. (9) is rewritten as
=LY+N(Y) ， (10)
where Y= (Y,, .): Using the exponential time difference scheme (5) to solve
eq. (10), Zhang et al. (2014) obtained a numerical scheme
Y™'' =exp(LAr)Y' +L' (exp(LAr)-1)N(¥') ， (11)
where At is the temporal step size. To improve the stability of the NETD
method, by replacing Y' with stv Y'), scheme (11) is reformed into
an implicit scheme
el (exp(LLar)-1)N y' =
1 (12)
[er + 5h (exp(LAr) “1)N}¥'
To explicitly solve the system of linear equations in scheme (12), it is
rewritten as
Y' =(1+B+B? I[exp(Lay) + sh (exp(LAr) “1)N}¥' ， (13)
by using the truncated power series expansion
(I-B)'=I+B+B’ ， (14)
where B= sh (exp(LA/) -UN
IMPROVED NEARLY-ANALYTIC EXPONENTIAL TIME
DIFFERENCE (INETD) METHOD
Now take a closer look at the approximation (4) used in the original
NETD method. The exponential time difference scheme (5) obtained by
using the approximation (4) is explicit and has only second order accuracy,
because of the use of low-order approximation (4). Therefore, the numerical
scheme (11) also becomes explicit, which results in instability, and has
second order accuracy in time.
In view of the above-mentioned points, we try to increase the stability
and the time accuracy of the original NETD, using the following
approximation instead of (4),
N(u(t, +5),¢, +s)=N(u't,)+2(N(u'stuu:)-N(u'st,)) , (15)
to obtain the exponential time difference scheme
u'™! =exp(Lr)u' + day (exp(Lz) -Lr-I)N'' +
T
1 (16)
=(L ) (LexptL7 -exp(LZ+DN'，
where N” =N(u',t,) , and N”™ =N(u't, +T =N(u™' th.) Using
the exponential time difference scheme (16) to solve eq. (10), gives the
numerical scheme
Y'*' = exp(LAt)Y” + “ (L')’(exp(LAt) — LAt -I)NY'' +
t
1 (17)
vi (L')’ (LAt exp(LAt) - exp(LAt) + NY”.
Eq. (17) can be rewritten as follows:
( ~ 元 (LTD2z(exp(LAD -LAI-DN]Y =
t
(18)
[sd 十元 (LY (LAtexp(LAD - exp(LAr) + DN | Y”.
Now, following the original NETD method, we can solve the system of
linear equations in scheme (18) explicitly. Let
C:= 元 (L ) (exp(LAr) - LAt -I)N
and use the truncated power series expansion
(I-C)'=1+C+C . (19)
Then, we can rewrite scheme (18) explicitly as
Y =(1+C+C’)
[ed 十 ~ (L') (LAr exp(LAt) - exp(LAt) + DN) 六 (20)
Comparing, we find that the INETD method (20) needs the same storage
space and computational costs as the NETD method (13). For example, the
total number of arrays required for computation in both is 12 for 1D and
18 for 2D in the acoustic content. There are differences in computational
costs between the coefficient matrices of the NETD and the INETD schemes,
but compared with the total computational costs this difference can be
neglected, because there is no need to compute the coefficient matrix at each
iteration, so it is sufficient to estimate it once before the iterations.
Therefore, computational costs of two methods are the same. Also, the
INETD method is an improvement over the original NETD method in both
the stability and the computational accuracy in time, as verified in our
theoretical analyses and numerical experiments below. The numerical
dispersion error of the INETD method is less than that of the NETD method.
The exponential time difference scheme (16) used in the INETD
method can be also obtained from the Lie group method. The Lie group
method is the geometric integration method and preserves qualitative
properties of differential equations.
STABILITY CRITERIA AND NUMERICAL DISPERSION ANALYSIS
Stability condition
To keep numerical calculation stable, the temporal increment At
must satisfy the stability condition of the INETD method. Following the
analysis process presented by Yang et al. (2006) and Lioyd (1996), we
obtain the stability condition of the INETD method for solving the acoustic
wave equations in 1D and 2D. Through a series of mathematical operations,
we obtain the following stability condition for the 1D homogeneous case,
Ara -~0.5178, (21)
Co Co
where fh denotes the space increment. For the 2D homogeneous case, the
stability condition of the INETD method (20) under the condition
Ax =Az=h is given by
Atsa,,. § ~03136. (22)
Co Co
Table 1. The approximate maximum Courant numbers of the original and improved
methods.
Method NETD INETD
@yax(1-D) 0.4132 0.5178
Oyyqy,(2-D) 0.2583 0.3136
For comparison, Table 1 contains the approximate maximum Courant
numbers of the NETD and INETD methods for 1D and 2D. From Table 1,
we can see the stability condition of the INETD method is more relaxed than
that of the NETD method.
Numerical dispersion analysis
In wavefield simulations, numerical dispersion (or grid dispersion) is
the most significant numerical problem limiting the usefulness of point-wise
discretization of wave equations. Following the analysis of previous work
(Dablain, 1986; Lioyd, 1996; Yang et al., 2006), we obtain the dispersion
relationship of the INETD method for 1D and 2D acoustic wave equations.
(a) (b)
4 41
一一 INETD 一一 INETD
NETD NETD
, | —— 加
. 「 ~ 性
9| 0.9
8 0.8
0 014 02 03 04 085 0 01 02 03 04 O85
Sp Sp
(c) (d)
44 1.1 -
一一 INETD 一一 INETD
05), -NETD 一 NETD
885 0.1 0.2 0.3 0.4 0.5 Oo 0.1 0.2 0.3 0.4 0.5
Sp Sp
Fig. 1. The ratio R of the numerical velocity to the exact phase velocity versus the
sampling rate Sp for the INETD and NETD methods for the 1D homogeneous case. The
Courant numbers are a 0.1, b 0.2, ec 0.3, and d 0.4, respectively .
Table 2. The average dispersive errors of the INETD and the original NETD methods for
different Courant numbers for the 1D Case.
Courant number (C )
Method 0.1 0.2 03 04
INETD 1.924 2.051 2.402 3.028
NETD 1.949 2.321 3.207 4.665
We define the spatial sampling ratio first proposed by Moszo et al.
(2000) as S, =h/A, in which A, denotes the wavelength. The dispersion
ratio R is defined as R=c,,,/c,, in which c,,, and c, denote the
numerical and true velocities, respectively. R=1 implies the method has no
numerical dispersion. The dispersion relation as a function of the sampling
rate (grid spacing per wavelength) is shown for the 1D case in Fig. 1. Table
2 shows the average numerical dispersive errors of the INETD and the
original NETD methods for different Courant numbers (Fig. 1). The average
numerical dispersive error is given by
5
f |I-R| ds,
Error =
x100(%).
From Fig. 1 and Table 2, we can see the average dispersion error and
the maximal dispersion error of the INETD method for the 1D case are less
than that of the NETD method. Fig. 2 shows the average numerical
dispersive errors of the INETD and the original NETD methods for different
Courant numbers and various values of propagation angles in the 2D case. In
Fig. 2 we observe that the average dispersion error of the INETD method is
less than that of the NETD.
57
、 NETD(a=0.25) .~” 7
人一
Average Numerical Dispersive Error
. Ss ,,
SNETD(a=0.1) 4
— Cee
INETD(o=0.1)
0 mW 2 308 48 5 6 78
Propagation Angle
Fig. 2. The average numerical dispersive errors of the INETD (solid line) and the original
NETD (dashed line) methods for different Courant numbers and propagation angles for
the 2D Case.
ERROR ANALYSIS
To investigate the accuracy of the INETD method, we estimate the
theoretical and relative numerical errors of the INETD method, and compare
with the NETD method.
Theoretical analysis
Using the Taylor series expansion, we conclude that the INETD
method is accurate third-order in time and fourth-order in space. Thus, the
temporal accuracy of the INETD method is increased from second order in
the original NETD method to third order, and the spatial accuracy remains
the same as that of the original.
Numerical errors
To further illustrate the accuracy of the INETD method, in the
following we compare the numerical errors caused by the INETD method
and the original NETD method for the 2D acoustic wave equation.
We consider the following 2D initial value problem:
au du 1 du
+ =
ax? dz? co Ot?
u(0, x, Z) = cos(- 2m fu “cos Ox - anf. ‘sin 0 *Z)
Co co » (23)
8u(0,%52) - —20 fi sin( 77 cos 4, :一 am}. ‘sin 0, 2)
ot Co Co
where cu is the velocity of the plane wave, / is the frequency, and 6, is
an incident angle at time +=0. Obviously, the analytical solution for the
initial value problem (23) is
Minx2) wes|2a/|1-Z0s0h sind) (24)
Co Co
For comparison, we use the INETD and the original NETD methods
to solve the initial value problem (23). In our numerical experiments, the
parameters are chosen as follows: the computational domain 0
0
the incident angle of the plane wave 6, =7/4. The relative error for the 2D
case is defined by
W 7 100 | (25)
where u;, is the numerical solution and u(t,,%,.2,) is the analytical
solution. To avoid the effects of artificial reflections from the edges of the
computational domain, in this numerical experiment we use the exact
solution (24) of the initial value problem (23) at the artificial boundaries in
each of the two methods.
Fig. 3 shows the computational results of the relative error E, at
different times for different spatial and temporal increments, where the two
lines of £, for the INETD and the original NETD methods are shown in a
semi-log scale. From Fig. 3, we conclude that the numerical errors
introduced by the INETD method measured by 五; are consistently less than
those of the original NETD method.
(a) Ax=Az=40m, At=0.002s b) Ax=Az=40m, At=0.0015s
10?
= pag treet ON)
S er
_
°
=
山
由
s
oO
o —— INETD —— INETD
---- NETD ---- NETD
10° 10°
0 0.4 0.8 1.2 0 0.4 0.8 1.2
c) Ax=Az=30m, At=0.001s
4102 spree WET
=
_
=
山
o
&
ao
& 一一 INETD 一一 INETD
---- NETD ---- NETD
107 10'
0 0.4 0.8 1.2 0 0.4 0.8 1.2
Time(s) Time(s)
Fig. 3. The relative errors of the INETD and the original NETD methods measured by 五/
[formula (25)] are shown in a semi-log scale for the 2D initial value problem (23). The
spatial and temporal increments are (a) Ax =Az=40m and Ar=0.002s,
(b) Ax = Az = 40m and Ar=0.0015ms, (¢) Ar=Az=30m and Ar=0.001s,
(d) Ax = Az =30m and Ar =0.00075s, respectively.
NUMERICAL DISPERSION AND EFFICIENCY
As shown above, the INETD method is superior in the numerical
errors to NETD. The computational costs of the two methods are the same:
for example, both took about 130 s to generate Fig. 3a, 173 s to generate
Fig.3b, 465 s to generate Fig. 3c, and 620 s to generate Fig. 3b. All the
experiments in this work are done on a Core(TM)2 Duo CPU(2.40 GHz,
39 GHz) with 1.89 GB RAM. To investigate the ability of the INETD
method to reduce numerical dispersion which is another important aspect
characterizing numerical schemes, wecompare the numerical solutions
generated by the INETD and the NETD methods, and by the fourth-order
LWC method (Dablain, 1986) with the results from the analytic method (Aki
and Richards, 1980) for the acoustic wave equation in a homogeneous
medium model.
In this example, we choose a homogeneous medium with an acoustic
velocity of co = 4 km/s. The model size is 0sx<7.2km and 0
and the source has a dominant frequency of fo = 20 Hz located at the center
of the model, and the receiver is located at (x=4.2km, z=3.6km). The
source time-function is a Ricker wavelet
f(t) =-5.76 f; (1-16(0.6 ft) exp(-8(0.6ft-1)) (26)
h=60m, At=0.0025s h=60m, At=0.002s
2
L 0.1 iA
2 | 一一人 [Lins
山， -041|[ 一一 INETD V
一 -一 NETD | 六
2 LWC
0 0.1 0.2 0.3 0.4 0.5
05 KR
5 f
-可 0 f
山 . 一一 INETD|V
08 || 一 -一 NETD | Y
ELWC
-0.1
0 01 02 03 04 05 0 01 02 03 04 05
Error
“oO of 8602 03 04 05 oO 01 02 03 O04 05
Time (s) Time (s)
Fig. 4. Numerical errors of the INETD, the NETD and the LWC methods at the receiver.
Note that these figures are plotted on different scales.
We choose the spatial increments h=Ax=Az= 60 m, 50 m and 40 m,
resulting in 3.3, 4 and 5 grid points per minimal wavelength, respectively.
Fig. 4 shows a comparison among the INETD, the NETD and the
LWC methods made of their numerical error on the different spatial and
temporal increments at the receiver.
From Fig. 4, we see that the INETD and the NETD methods have
lower numerical dispersion than that of the fourth-order LWC method. And
we observe that the errors and the numerical dispersion of the INETD are the
smallest of the three methods. Meanwhile, the errors and the numerical
dispersion of the INETD are gradually smaller as the spatial and temporal
grid sizes decrease, indicating both the convergence and the validity of the
INETD.
NUMERICAL EXAMPLES
In this section we present four numerical examples to investigate the
validity of the INETD method. In all the experiments, the model parameters
are the same as those used in Zhang et al. (2014), and the explosive source
[eq. (26)] is used.
Example 1: Three-layer model
First we choose a computational domain of 0sx<6 km and 0
km and a three-layer model with the two horizontal interfaces at 2 km and
33 km depth, respectively. The velocity is 3.0 km/s in the top layer, 2.0
km/s in the middle layer, and 4.0 km/s in the bottom layer. The explosive
source [eq. (26)] with a dominant frequency of fo = 25 Hz is located at the
center of the domain and spatial and temporal increments are chosen as
Ax=Az=20m and At=5x10~“s, resulting in 4 grid points per minimum
wavelength.
Fig. 5 shows the wavefield snapshots generated by INETD at four
times; these show the wave propagation phenomena including reflection and
transmission at the interfaces without visible numerical dispersion. This
experiment demonstrates that INETD can deal well with models that contain
strong velocity contrasts between the adjacent layers on coarse spatial grids.
Distance (km) Distance (km)
©
Depth (km)
(a) T=0.3s ° (b) T=0.5s
Depth (km)
(d) T=1.0s
Fig. 5. Snapshots of acoustic wavefields at time a 0.3 s, b 0.5 s, ¢ 0.8 s, and d 1.0 s for
the three-layer model generated by the INETD method.
Example 2: Marmousi model
In the second experiment, we choose the Marmousi model (Versteeg
and Grau, 1991) to demonstrate the performance of the INETD method in
realistically complicated media. Fig. 6 shows the velocity structure c(x,z),
where the velocity varies from 1.5—5.5 km/s. We choose the spatial
increments Ax =Az=24m and the temporal increment Ar =0.000436s. The
number of grid points is 384x122 so that the computational domain is
0
frequency of fo = 15 Hz, is located in the middle on the surface, resulting in
17 grid points per minimal wavelength.
Distance (km)
Fig. 6. The Marmousi model.
Fig. 7 shows the wavefield snapshots generated by the INETD at
times T=0.6s, T=0.9s, T=1.2s, and 7 =1.5s, respectively; these snapshots
are clean and have no visible numerical dispersion for the complicated
geological model even though the spatial increment is large. In this
experiment, we use the perfectly matched layer absorbing boundary
condition (Komatitsch and Tromp, 2003; Ma et al., 2015) to treat the grid
edges.
Example 3: a homogeneous elastic model
We consider the following elastic wave equations in a 2D
homogeneous TI medium:
au ou a°w ou
Pa = Cl (cb +Cy4) axdz Cy ath
aw 2 2 aw ， (27)
w 0 1
Pe = Cu aa tC + Cu) tO oa th
Distance (km)
0 9.192
W
Depth (km)
904
(a) T=0.6s
局
192
Depth (km)
904
全
192
Depth (km)
904
0 9.192
Depth (km)
和
(d) T=1.5s
904
Fig. 7. Snapshots of acoustic wave fields at time a 0.6 s, b 0.9 s, ¢ 1.2 s, and d 1.5 s for
the Marmousi model generated, by the INETD method.
where uw and w are the displacement components in the x- and
z-directions, respectively. cll, C13, C33 and c44 are the elastic constants; p is
the medium density; f and f, are the source force components in the
x-and z-directions. The computational parameters are cll = 45, cl3 = 9.6,
033 = 37.5, and c44 = 12 GPa, and p=1.0 g/cm’. The mesh size is 401x401
points, the spatial grid increment is Av =Az=30 m, and the time increment
is At=0.00067 s.
The source located at the center of the computational domain is a
Ricker wavelet with a dominant frequency of 1 = 30 Hz. Fig. 8 shows the
x- and z-component snapshots at T=0.6s generated by the INETD method.
The snapshots show clear wavefield information of the displacement without
visible numerical dispersion. The cusps and the anisotropy of the velocity of
the wave propagation can be clearly seen. In short, the numerical results
demonstrate that the INETD method can simulate the elastic wave
propagation effectively.
Distance (km) Distance (km)
Depth (km)
Depth (km)
Fig. 8. Snapshots of elastic wave fields in the 2D homogeneous TI medium at time 0.6 s,
generated by the INETD method for (left) the x-direction displacement component and
(right) the z-direction displacement component.
Example 4: a two-layer elastic model
For this example, we choose a two-layer elastic model with large
velocity contrasts across the interface. The size of the computational domain
is O
layers is at 2.4 km depth. The Lamé constants and densities of the two-layer
model are p,=1.5g/em*, 4,=1.5GPa, y,=2.5GPa in the top layer, and
p,=2.0g/em*, A,=11.0 GPa, w,=15.0GPa in the bottom layer. These
correspond to the P-and SV -wave velocities of 1.63 and 1.29 km/s in the
top layer, and 3.61 and 2.74 km/s in the bottom layer. The chosen spatial
increment is Av =Az=15m, and the chosen temporal increment is Ar=0.0005 s.
The source, which has a frequency of /) = 20 Hz, is a symmetric Ricker
wavelet located at (x=2km, z=1.85km).
Fig. 9 shows the x-direction and z-direction displacement snapshots at
T =0.6s, generated by the INETD. The numerical results clearly show the
elastic wave propagation phenomena in the two-layer model, which
indicates that the INETD can also be used in a multilayer medium with
strong interfaces for elastic wave propagation.
Distance (km) Distance (km)
Depth (km)
Depth (km)
Fig. 9. Snapshots of elastic wavefields in the two-layer medium at time 0.6 s, generated
by the INETD method for (left) the x-direction and (right) the z-direction displacement
components.
ANOTHER CONSIDERATION OF THE INETD SCHEME
Now take a closer look at the INETD scheme (20). We can replace
approximation (19) by the approximation
(I-C)'=14+C+C+C , (28)
resulting in another INETD scheme as
Y = (I+C+C +C’)
1 (29)
[expan + rei (L ) (LAt exp(LAt) — exp(LAt) + DN ]Y'.
Table 3. The approximate maximum Courant numbers of the original and improved
methods.
Method NETD INETD1 INETD2
w (1-D) 0.4132 0.5178 0.8263
Ogg, (2-D) 0.2583 0.3136 0.7811
Table 3 shows the approximate maximum Courant numbers of the
NETD scheme (13), the INETD scheme (20) (INETD1) and the INETD
scheme (29) (INETD2) for the 1D and 2D cases. From Table 3, we see that
the INETD scheme (29) exhibits higher stability limits than the INETD
scheme (20) and the NETD scheme (13). The INETD scheme (29) gives
more accurate numerical results than the INETD scheme (20) and the NETD
scheme (13) when the same grids are used. However, the INETD scheme
(29) increases the computational cost by about 30% compared with the
INETD scheme (20).
CONCLUSIONS
We have improved the original NETD for solving acoustic and elastic
wave equations, and investigated in detail the properties of the INETD
method, including its stability, numerical dispersion, and accuracy. The
INETD method, compared with the NETD method, has more relaxed
stability. In 1D, the approximate maximum Courant number of the INETD
method is increased from 0.4132 in the original NETD method to 0.5178
(for the INETD scheme (29) 0.8263), and in 2D, from 0.2583 to 0.3136 (for
the INETD scheme (29) 0.7811). The INETD method increases the accuracy
in time from the second order of the original NETD method to the third
order, while keeping the spatial accuracy the same as the original NETD
method. This conclusion is also verified by our numerical experiments of
computing the relative error E,, via formula (25) for the 2D case. In addition,
the dispersion error of the INETD method for the 1D and 2D cases is less
than that of the NETD method.
When using the same temporal and spatial increments, the
computational costs and storage space of the INETD (20) and the NETD
methods are the same, while the INETD scheme (29) increases the
computational cost by about 30% compared with the NETD scheme.
However, because the INETD scheme (29) has a much higher stability limit,
higher order of temporal accuracy and lower dispersion than the NETD, it
can lead to high computational efficiency.
Meanwhile, wavefield modelling illustrates that the improved method
can effectively reduce numerical dispersion when too-coarse computation
grids are used. The INETD method is related to the Lie group method.
ACKNOWLEDGEMENTS
We thank the editor and an anonymous reviewer for their detailed
comments and helpful revisions that greatly improved this paper. We also
thank Prof. Dinghui Yang for his fruitful suggestions that enabled us to
highly improve the paper.
REFERENCES
Aki, K. and Richards, P.G., 1980. Quantitative Seismology: Theory and Methods. W.H.
Freeman and Co., San Francisco.
Chen, S., Yang, D.H. and Deng, X.Y., 2010. A weighted Runge-Kutta method with weak
numerical dispersion for solving wave equations. Commun. Comput. Phys., 7:
1027-1048.
Dablain, M.A., 1986. The application of high-order differencing to the scalar wave
equation. Geophysics, 51: 54-66.
Fei, T. and Larner, K., 1995. Elimination of numerical dispersion in finite difference
modelling and migration by flux-corrected transport. Geophysics, 60: 1830-1842.
Feng, L.L., Yang, D.H. and Xie, W., 2015. An efficient symplectic reverse time
migration method using a local nearly-analytic discrete operator in acoustic
transversely isotropic media with a vertical symmetry axis. Geophysics, 80(2):
$103-S112.
Fornberg, B., 1998. A Practical Guide to Psuedospectral Methods. Cambridge University
Press, Cambridge.
Hairer, E., Lubich, C.H. and Wanner, G., 2006. Geometric Numerical Integration. 2'4 ed.,
Springer Verlag, Berlin.
Huang, B.S., 1992. A program for two-dimensional seismic wave propagation by the
pseudospectrum method. Comput. Geosci., 18 : 289-307.
Iserles, A., Munthe-Kaas, H., Norsett, S.P. and Zanna, A., 2000. Lie group methods. Acta
Numer., 9: 215-365.
Kelly, K., Ward, R., Treitel, S. and Alford, R., 1976. Synthetic seismograms: a
finite-difference approach. Geophysics, 41: 2-27.
Komatitsch and Tromp, J., 2003. A perfectly matched layer absorbing boundary
condition for the second-order seismic wave equation. Geophys. J. Internat., 154:
146-153.
Kosloff, D. and Baysal, E., 1982. Forward modeling by a Fourier method. Geophysics,
47: 1402-1412.
Kosloff, D., Reshef, M. and Loewenthal, D., 1984. Elastic wave calculations by the
Fourier method. Bull. Seismol. Soc. Am., 74 : 875-891.
Li, J.S., Yang, D.H., Wu, H. and Ma, X., 2017. A low-dispersive method using the
high-order stereo-modeling operator for solving 2D wave equations. Geophys. J.
Internat., 210: 1938-1964.
Li, J.S., Yang, D.H. and Liu, F.Q., 2013. An efficient reverse-time migration method
using local nearly analytic discrete operator. Geophysics, 78(1): $15-S23.
Liu, S.L., Yang, D.H., Lang, C., Wang, W.S. and Pan, Z.D., 2017. Modified symplectic
schemes with nearly-analytic discrete operators for acoustic wave simulations.
Comput. Phys. Communicat., 213: 52-63.
Lioyd, N., 1996. Finite Difference and Spectral Methods for Ordinary and Partial
Differential Equations. Ph.D. thesis, Cornell University, Ithaca.
Ma, X., Yang, D.H. and Liu, F., 2011. A nearly analytic symplectically partitioned
Runge-Kutta method for 2-D seismic wave equations. Geophys. J. Internat., 187:
480-496.
Ma, X., Yang, D.H., Song, G.J. and Wang, M.X., 2014. A low-dispersive symplectic
partitioned Runge-Kutta method for solving seismic wave equations - I. Scheme
and theoretical analysis. Bull. Seismol. Soc. Am., 104: 2206-2225.
Ma, X., Yang, D.H. and Song, G.J., 2015. A low-dispersive symplectic partitioned
Runge-Kutta method for solving seismic wave equations - II. Wave-field
simulations. Bull. Seismol. Soc. Am., 105: 657-675.
Mizutani, H., Geller, R.J. and Takeuchi, N., 2000. Comparison of accuracy and efficiency
of time domain schemes for calculating synthetic seismograms. Phys. Earth Planet.
Inter., 119 : 75-97.
Moszo, P., Kristek, J. and Halada, L., 2000. 3-D fourth-order staggered-grid
finite-difference schemes: stability and grid dispersion. Bull. Seismol. Soc. Am., 90:
587-603.
Munthe-Kaas, H., 1999. High order Runge-Kutta methods on manifolds. Appl. Numer.
Math., 29: 115-227.
Robert, L.B., 1991. An Introduction to Lie Groups and Symplectic Geometry. Ph.D.thesis,
Duke University, Park City.
Tong, P., Yang D.H., Hua, B.L. and Wang, M.X., 2013. A high-order stereo-modeling
method for solving wave equations. Bull. Seismol. Soc. Am., 103: 811-833.
Versteeg, R.J. and Grau, G., 1991. The Marmousi experience. Proc. EAGE workshop
Practical Aspects of Seismic Data Inversion, Copenhagen, 1990.
Wang, S.Q., Yang, D.H. and Yang, K.D., 2002. Compact finite difference scheme for
elastic equations. J. Tsinghua Univ., 42: 1128-1131 (in Chinese).
Xie, W., Yang, D.H., Liu, F.Q. and Li, J.S., 2014. Reverse-time migration in acoustic
VTI media using a high-order stereo operator. Geophysics, 79(3): WA3-WA11.
Yang, D.H., Liu, E.R., Zhang, Z.J. and Teng, J.W., 2002. Finite-difference modelling in
2D anisotropic media using a flux corrected transport technique. Geophys. J.
Internat., 148: 320-328.
Yang, D.H., Lu, M., Wu, R.S. and Peng, J. M., 2004. An optimal nearly analytic discrete
method for 2D acoustic and elastic wave equations. Bull. Seismol. Soc. Am., 94:
1982-1991.
Yang, D.H., Tong, P. and Deng, X.Y., 2012. A central difference method with low
numerical dispersion for solving the scalar wave equation. Geophys. Prosp., 60:
885-905.
Yang, D.H., Teng, J.W., Zhang, Z.J. and Liu, E., 2003. A nearly-analytic discrete method
for acoustic and elastic wave equations in anisotropic media. Bull. Seismol. Soc.
Am., 93: 882-890.
Yang, D.H., Peng, J.M., Lu, M. and Terlaky, T., 2006. Optimal nearly analytic discrete
approximation to the scalar wave equation. Bull. Seismol. Soc. Am., 96:
1114-1130.
Yang, D.H., Song, G.J. and Lu, M., 2007. Optimally accurate nearly-analytic discrete
scheme for wave-field simulation in 3D anisotropic media. Bull. Seismol. Soc. Am.,
97: 1557-1569.
Yang, D.H., Wang, N., Chen, S. and Song, G.J., 2009. An explicit method based on the
implicit Runge- -Kutta algorithm for solving the wave equations. Bull. Seismol. Soc.
Am., 99: 3340-3354.
Yang, DH. and Wang, L., 2010a. A split-step algorithm with effectively suppressing the
numerical dispersion for 3D seismic propagation modeling. Bull. Seismol. Soc.
Am., 100: 1470-1484.
Yang, D.H., Song, G.J., Hua, B.L. and Henri, C., 2010b. Simulation of acoustic
wave-fields in heterogeneous media: A robust method with automatically
suppressing numerical dispersion for large grid steps. Geophysics, 75(3):
T99-T1

Previous article in this issue

Next article in this issue

Journal of Seismic Exploration, Electronic ISSN: 0963-0651 Print ISSN: 0963-0651, Published by AccScience Publishing

Publisher's Core Philosophy

We are committed to support the scientific community by publishing impactful research and enhancing communication among scientists. At AccScience Publishing, we are continuously looking for ways to accelerate scientific progress and to strive for transparency and open communication, making knowledge freely accessible without barrier.

9 Raffles Place, Republic Plaza 1 #06-00 Singapore 048619

+65 8182 1586

editorial@accscience.com