AccScience Publishing / JSE / Article
Submit to JSE
Cite this article
3
Download
47
Views
Journal Browser
Volume | Year
Issue
Search
Forthcoming Issue
Current Issue
View All
News and Announcements
View All
ARTICLE

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

NAM YUN,  CHOL SIM
Show Less
Institute of Mathematics, State Academy of Sciences, Pyongyang, D.P.R. Korea.,
JSE 2020, 29(2), 99–124;
Submitted: 11 July 2018 | Accepted: 12 December 2019 | Published: 1 April 2020
© 2020 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
  1. Thus, combining of geometric integrations with NAD may lead to decreasesin numerical dispersion and numerical error over other classical methods.
  2. Motivated by this fact, Zhang et al. (2014) developed a nearly-analyticexponential time difference (NETD) method. NETD uses an exponentialtime difference (ETD) scheme to approximate the time derivatives and thenearly-analytic discrete (NAD) operator to approximate the spacederivatives. The ETD scheme used in this method can be also derived fromthe idea of the Lie group method. Therefore, the NETD method bycombining the Lie group method with the NAD method to solve seismicwave equations effectively suppresses the numerical dispersion andpreserves the qualitative properties of PDEs. Through theoretical andnumerical analysis, they showed that NETD method provides more accurateresults, has high computational efficiency, and a low memory requirementcompared to the Lax-Wendroff correction (LWC) and the staggered-grid(SG) methods. However, the NETD method has some disadvantages; forexample, its stability criterion is relatively compact and has onlysecond-order time accuracy. It is important to improve the stability of the
  3. NETD method and increase its time accuracy for practical calculations.
  4. From this view point, we focus on this problem. We present animproved NETD (INETD), for solving acoustic- and elastic-wave equations.
  5. The INETD scheme is composed of the following two steps: (1) Derivation anew implicit ETD scheme, (ii) Combination of the implicit ETD for thetemporal discretization with the NAD method for the spatial discretization.
  6. This produces an improvement in the stability and the time accuracy overthe NETD, because the explicit ETD scheme is used to approximate the timederivatives in the NETD. In general, the implicit scheme has a good stabilityproperty, but it requires more expensive computations than the explicitscheme, because it solves the system of equations at every temporal step.
  7. We show that the INETD method needs the same computational costs andstorage space as the NETD method. On the other hand, we find anotherpossibility for improving the stability of the INETD method. It is to use theterms up to the third order in the truncated power series expansion used toexplicitly solve the INETD system of linear equations. Our analysis shows itbrings about a remarkable improvement in the stability of the method. Onthe whole, the our INETD method are shown to be superior to the FDM, the
  8. PS and the NETD methods in some aspects such as suppressing numericaldispersions, high accuracy, computational efficiency, a low memory and theimprovement of stability.
  9. This article is organized as follows. First, we review the basic NETDscheme for solving the acoustic wave equation in a 2D homogeneousmedium. Then, we derive the INETD scheme. Next, we investigate in detailthe properties of the INETD method. We derive the stability condition forthe 1D and 2D cases. We estimate the theoretical and numerical errors of the
  10. INETD method, and compare against the NETD method. And then, wederive the dispersion relation and give the numerical dispersion results.
  11. Finally, for the 2D case we use the INETD, the original NETD and thefourth order LWC methods to simulate acoustic wave fields, and presentsome comparisons of numerical results for different methods. And then ahomogeneous elastic model and a two-layer elastic model are furtherselected to investigate the computational validity of elastic wave-fieldsimulations.
  12. Our theoretical and numerical results show that the INETD is animprovement over the original NETD in stability, numerical error andnumerical dispersion.BASIC NEARLY-ANALYTIC EXPONENTIAL TIME DIFFERENCE(NETD) METHOD
  13. In this section, we review and analyze the derivation of the original
  14. NETD method (Zhang et al., 2014) for solving the acoustic wave equation.Exponential time difference methodConsider the stiff systemdu—=Lu+N(u,talt (ut), (1)where L is a higher-order linear term and N is a lower-order nonlinearterm. 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 ¢,,. Toapproximate the integration in eq. (3), using the approximation asN(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
  15. The acoustic wave equation in 2D homogeneous medium is written asduwe co Au, (6)Cewhere z and c, are the displacement and the acoustic velocity, respectively,2 2oe, . .and A= Pye + ra is the Laplace operator. The same notations as that inx azthe basic NETD (Zhang et al., 2014) are used in our present study, i.e.,7 7ou du 07/ dw dwarya 0 人 “有 二 | > we (me .tEq. (6) can be rewritten asou -WwWataW _ aU aU)’ (7)ot 0 ax? az?
  16. The high-order spatial derivatives in eq. (7) are determined by thedisplacement U and its gradients using the fourth-order NAD operator asshown in Appendix A of Zhang et al. (2014). Then eq. (7) can betransformed into an ODE system as follows:d(U) fo 1(U) <| ae @
  17. The first part of the right side of eq. (8) is a linear combination of u,wwith its gradients at the grid point (x;, z,), and the second part is a linearcombination of u,w with its gradients at its neighboring grid points. I isa3x3 identity matrix and 2 is a 3x3 matrix that depends on the spatialdiscrete operator. Using the notations7 - (0 IY,,=(U,W)',, and 人and regarding N(Y,, ) as a nonlinear function of Y,,, eq. (8) can berewritten asqY,, - ~的 + N(Y, 外 ， (9)Eq. (9) is rewritten as =LY+N(Y) ， (10)where Y= (Y,, .): Using the exponential time difference scheme (5) to solveeq. (10), Zhang et al. (2014) obtained a numerical scheme
  18. Y™'' =exp(LAr)Y' +L' (exp(LAr)-1)N(¥') ， (11)where At is the temporal step size. To improve the stability of the NETDmethod, by replacing Y' with stv Y'), scheme (11) is reformed intoan implicit schemeel (exp(LLar)-1)N y' =1 (12)[er + 5h (exp(LAr) “1)N}¥'
  19. To explicitly solve the system of linear equations in scheme (12), it isrewritten as
  20. 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/) -UNIMPROVED NEARLY-ANALYTIC EXPONENTIAL TIMEDIFFERENCE (INETD) METHOD
  21. Now take a closer look at the approximation (4) used in the original
  22. NETD method. The exponential time difference scheme (5) obtained byusing the approximation (4) is explicit and has only second order accuracy,because of the use of low-order approximation (4). Therefore, the numericalscheme (11) also becomes explicit, which results in instability, and hassecond order accuracy in time.
  23. In view of the above-mentioned points, we try to increase the stabilityand the time accuracy of the original NETD, using the followingapproximation instead of (4),
  24. N(u(t, +5),¢, +s)=N(u't,)+2(N(u'stuu:)-N(u'st,)) , (15)to obtain the exponential time difference schemeu'™! =exp(Lr)u' + day (exp(Lz) -Lr-I)N'' +T1 (16)=(L ) (LexptL7 -exp(LZ+DN'，where N” =N(u',t,) , and N”™ =N(u't, +T =N(u™' th.) Usingthe exponential time difference scheme (16) to solve eq. (10), gives thenumerical scheme
  25. Y'*' = exp(LAt)Y” + “ (L')’(exp(LAt) — LAt -I)NY'' +t1 (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”.
  26. Now, following the original NETD method, we can solve the system oflinear equations in scheme (18) explicitly. LetC:= 元 (L ) (exp(LAr) - LAt -I)Nand use the truncated power series expansion(I-C)'=1+C+C . (19)Then, we can rewrite scheme (18) explicitly asY =(1+C+C’)[ed 十 ~ (L') (LAr exp(LAt) - exp(LAt) + DN) 六 (20)
  27. Comparing, we find that the INETD method (20) needs the same storagespace and computational costs as the NETD method (13). For example, thetotal number of arrays required for computation in both is 12 for 1D and18 for 2D in the acoustic content. There are differences in computationalcosts between the coefficient matrices of the NETD and the INETD schemes,but compared with the total computational costs this difference can beneglected, because there is no need to compute the coefficient matrix at eachiteration, so it is sufficient to estimate it once before the iterations.
  28. Therefore, computational costs of two methods are the same. Also, the
  29. INETD method is an improvement over the original NETD method in boththe stability and the computational accuracy in time, as verified in ourtheoretical analyses and numerical experiments below. The numericaldispersion error of the INETD method is less than that of the NETD method.
  30. The exponential time difference scheme (16) used in the INETDmethod can be also obtained from the Lie group method. The Lie groupmethod is the geometric integration method and preserves qualitativeproperties of differential equations.
  31. STABILITY CRITERIA AND NUMERICAL DISPERSION ANALYSISStability condition
  32. To keep numerical calculation stable, the temporal increment Atmust satisfy the stability condition of the INETD method. Following theanalysis process presented by Yang et al. (2006) and Lioyd (1996), weobtain the stability condition of the INETD method for solving the acousticwave 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 Cowhere fh denotes the space increment. For the 2D homogeneous case, thestability condition of the INETD method (20) under the conditionAx =Az=h is given byAtsa,,. § ~03136. (22)Co Co
  33. Table 1. The approximate maximum Courant numbers of the original and improvedmethods.Method NETD INETD@yax(1-D) 0.4132 0.5178Oyyqy,(2-D) 0.2583 0.3136
  34. For comparison, Table 1 contains the approximate maximum Courantnumbers 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 thanthat of the NETD method.Numerical dispersion analysis
  35. In wavefield simulations, numerical dispersion (or grid dispersion) isthe most significant numerical problem limiting the usefulness of point-wisediscretization of wave equations. Following the analysis of previous work(Dablain, 1986; Lioyd, 1996; Yang et al., 2006), we obtain the dispersionrelationship of the INETD method for 1D and 2D acoustic wave equations.(a) (b)4 41一 一 INETD 一 一 INETDNETD NETD, | —— 加. 「 ~ 性9| 0.98 0.80 014 02 03 04 085 0 01 02 03 04 O85Sp Sp(c) (d)44 1.1 -一 一 INETD 一 一 INETD05), -NETD 一 NETD885 0.1 0.2 0.3 0.4 0.5 Oo 0.1 0.2 0.3 0.4 0.5Sp Sp
  36. Fig. 1. The ratio R of the numerical velocity to the exact phase velocity versus thesampling rate Sp for the INETD and NETD methods for the 1D homogeneous case. The
  37. Courant numbers are a 0.1, b 0.2, ec 0.3, and d 0.4, respectively .
  38. Table 2. The average dispersive errors of the INETD and the original NETD methods fordifferent Courant numbers for the 1D Case.Courant number (C )Method 0.1 0.2 03 04INETD 1.924 2.051 2.402 3.028NETD 1.949 2.321 3.207 4.665
  39. We define the spatial sampling ratio first proposed by Moszo et al.(2000) as S, =h/A, in which A, denotes the wavelength. The dispersionratio R is defined as R=c,,,/c,, in which c,,, and c, denote thenumerical and true velocities, respectively. R=1 implies the method has nonumerical dispersion. The dispersion relation as a function of the samplingrate (grid spacing per wavelength) is shown for the 1D case in Fig. 1. Table2 shows the average numerical dispersive errors of the INETD and theoriginal NETD methods for different Courant numbers (Fig. 1). The averagenumerical dispersive error is given by5f |I-R| ds,Error =x100(%).
  40. From Fig. 1 and Table 2, we can see the average dispersion error andthe maximal dispersion error of the INETD method for the 1D case are lessthan that of the NETD method. Fig. 2 shows the average numericaldispersive errors of the INETD and the original NETD methods for different
  41. Courant numbers and various values of propagation angles in the 2D case. In
  42. Fig. 2 we observe that the average dispersion error of the INETD method isless than that of the NETD.57、 NETD(a=0.25) .~” 7人 一Average Numerical Dispersive Error. Ss ,,SNETD(a=0.1) 4— CeeINETD(o=0.1)0 mW 2 308 48 5 6 78Propagation Angle
  43. Fig. 2. The average numerical dispersive errors of the INETD (solid line) and the original
  44. NETD (dashed line) methods for different Courant numbers and propagation angles forthe 2D Case.ERROR ANALYSIS
  45. To investigate the accuracy of the INETD method, we estimate thetheoretical and relative numerical errors of the INETD method, and comparewith the NETD method.Theoretical analysis
  46. Using the Taylor series expansion, we conclude that the INETDmethod is accurate third-order in time and fourth-order in space. Thus, thetemporal accuracy of the INETD method is increased from second order inthe original NETD method to third order, and the spatial accuracy remainsthe same as that of the original.Numerical errors
  47. To further illustrate the accuracy of the INETD method, in thefollowing we compare the numerical errors caused by the INETD methodand the original NETD method for the 2D acoustic wave equation.
  48. 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 Cowhere cu is the velocity of the plane wave, / is the frequency, and 6, isan incident angle at time +=0. Obviously, the analytical solution for theinitial value problem (23) isMinx2) wes|2a/|1-Z0s0h sind) (24)Co Co
  49. For comparison, we use the INETD and the original NETD methodsto solve the initial value problem (23). In our numerical experiments, theparameters are chosen as follows: the computational domain 0
  50. Fig. 3 shows the computational results of the relative error E, atdifferent times for different spatial and temporal increments, where the twolines of £, for the INETD and the original NETD methods are shown in asemi-log scale. From Fig. 3, we conclude that the numerical errorsintroduced by the INETD method measured by 五; are consistently less thanthose of the original NETD method.(a) Ax=Az=40m, At=0.002s b) Ax=Az=40m, At=0.0015s10?= pag treet ON)S er_°=山由soOo —— INETD —— INETD---- NETD ---- NETD10° 10°0 0.4 0.8 1.2 0 0.4 0.8 1.2c) Ax=Az=30m, At=0.001s4102 spree WET=_=山o&ao& 一 一 INETD 一 一 INETD---- NETD ---- NETD107 10'0 0.4 0.8 1.2 0 0.4 0.8 1.2Time(s) Time(s)
  51. 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). Thespatial 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
  52. As shown above, the INETD method is superior in the numericalerrors 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
  53. Fig.3b, 465 s to generate Fig. 3c, and 620 s to generate Fig. 3b. All theexperiments 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 INETDmethod to reduce numerical dispersion which is another important aspectcharacterizing numerical schemes, wecompare the numerical solutionsgenerated by the INETD and the NETD methods, and by the fourth-order
  54. LWC method (Dablain, 1986) with the results from the analytic method (Akiand Richards, 1980) for the acoustic wave equation in a homogeneousmedium model.
  55. In this example, we choose a homogeneous medium with an acousticvelocity of co = 4 km/s. The model size is 0sx<7.2km and 0
  56. Fig. 4. Numerical errors of the INETD, the NETD and the LWC methods at the receiver.
  57. Note that these figures are plotted on different scales.
  58. 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.
  59. Fig. 4 shows a comparison among the INETD, the NETD and the
  60. LWC methods made of their numerical error on the different spatial andtemporal increments at the receiver.
  61. From Fig. 4, we see that the INETD and the NETD methods havelower numerical dispersion than that of the fourth-order LWC method. Andwe observe that the errors and the numerical dispersion of the INETD are thesmallest of the three methods. Meanwhile, the errors and the numericaldispersion of the INETD are gradually smaller as the spatial and temporalgrid sizes decrease, indicating both the convergence and the validity of theINETD.NUMERICAL EXAMPLES
  62. In this section we present four numerical examples to investigate thevalidity of the INETD method. In all the experiments, the model parametersare the same as those used in Zhang et al. (2014), and the explosive source[eq. (26)] is used.Example 1: Three-layer model
  63. First we choose a computational domain of 0sx<6 km and 0
  64. Ax=Az=20m and At=5x10~“s, resulting in 4 grid points per minimumwavelength.
  65. Fig. 5 shows the wavefield snapshots generated by INETD at fourtimes; these show the wave propagation phenomena including reflection andtransmission at the interfaces without visible numerical dispersion. Thisexperiment demonstrates that INETD can deal well with models that containstrong velocity contrasts between the adjacent layers on coarse spatial grids.Distance (km) Distance (km)©Depth (km)(a) T=0.3s ° (b) T=0.5sDepth (km)(d) T=1.0s
  66. Fig. 5. Snapshots of acoustic wavefields at time a 0.3 s, b 0.5 s, ¢ 0.8 s, and d 1.0 s forthe three-layer model generated by the INETD method.Example 2: Marmousi model
  67. In the second experiment, we choose the Marmousi model (Versteegand Grau, 1991) to demonstrate the performance of the INETD method inrealistically 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 spatialincrements Ax =Az=24m and the temporal increment Ar =0.000436s. Thenumber of grid points is 384x122 so that the computational domain is0
  68. Fig. 7 shows the wavefield snapshots generated by the INETD attimes T=0.6s, T=0.9s, T=1.2s, and 7 =1.5s, respectively; these snapshotsare clean and have no visible numerical dispersion for the complicatedgeological model even though the spatial increment is large. In thisexperiment, we use the perfectly matched layer absorbing boundarycondition (Komatitsch and Tromp, 2003; Ma et al., 2015) to treat the gridedges.Example 3: a homogeneous elastic model
  69. We consider the following elastic wave equations in a 2Dhomogeneous TI medium:au ou a°w ouPa = Cl (cb +Cy4) axdz Cy athaw 2 2 aw ， (27)w 0 1Pe = Cu aa tC + Cu) tO oa thDistance (km)0 9.192WDepth (km)904(a) T=0.6s局192Depth (km)904全192Depth (km)9040 9.192Depth (km)和(d) T=1.5s904
  70. 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 forthe Marmousi model generated, by the INETD method.where uw and w are the displacement components in the x- andz-directions, respectively. cll, C13, C33 and c44 are the elastic constants; p isthe medium density; f and f, are the source force components in thex-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 401x401points, the spatial grid increment is Av =Az=30 m, and the time incrementis At=0.00067 s.
  71. The source located at the center of the computational domain is a
  72. Ricker wavelet with a dominant frequency of 1 = 30 Hz. Fig. 8 shows thex- and z-component snapshots at T=0.6s generated by the INETD method.
  73. The snapshots show clear wavefield information of the displacement withoutvisible numerical dispersion. The cusps and the anisotropy of the velocity ofthe wave propagation can be clearly seen. In short, the numerical resultsdemonstrate that the INETD method can simulate the elastic wavepropagation effectively.Distance (km) Distance (km)Depth (km)Depth (km)
  74. 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
  75. For this example, we choose a two-layer elastic model with largevelocity contrasts across the interface. The size of the computational domainis O
  76. The source, which has a frequency of /) = 20 Hz, is a symmetric Rickerwavelet located at (x=2km, z=1.85km).
  77. Fig. 9 shows the x-direction and z-direction displacement snapshots at
  78. T =0.6s, generated by the INETD. The numerical results clearly show theelastic wave propagation phenomena in the two-layer model, whichindicates that the INETD can also be used in a multilayer medium withstrong interfaces for elastic wave propagation.Distance (km) Distance (km)Depth (km)Depth (km)
  79. Fig. 9. Snapshots of elastic wavefields in the two-layer medium at time 0.6 s, generatedby the INETD method for (left) the x-direction and (right) the z-direction displacementcomponents.ANOTHER CONSIDERATION OF THE INETD SCHEME
  80. Now take a closer look at the INETD scheme (20). We can replaceapproximation (19) by the approximation(I-C)'=14+C+C+C , (28)resulting in another INETD scheme asY = (I+C+C +C’)1 (29)[expan + rei (L ) (LAt exp(LAt) — exp(LAt) + DN ]Y'.
  81. Table 3. The approximate maximum Courant numbers of the original and improvedmethods.Method NETD INETD1 INETD2w (1-D) 0.4132 0.5178 0.8263Ogg, (2-D) 0.2583 0.3136 0.7811
  82. Table 3 shows the approximate maximum Courant numbers of the
  83. NETD scheme (13), the INETD scheme (20) (INETD1) and the INETDscheme (29) (INETD2) for the 1D and 2D cases. From Table 3, we see thatthe INETD scheme (29) exhibits higher stability limits than the INETDscheme (20) and the NETD scheme (13). The INETD scheme (29) givesmore accurate numerical results than the INETD scheme (20) and the NETDscheme (13) when the same grids are used. However, the INETD scheme(29) increases the computational cost by about 30% compared with theINETD scheme (20).CONCLUSIONS
  84. We have improved the original NETD for solving acoustic and elasticwave equations, and investigated in detail the properties of the INETDmethod, including its stability, numerical dispersion, and accuracy. The
  85. INETD method, compared with the NETD method, has more relaxedstability. In 1D, the approximate maximum Courant number of the INETDmethod 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 (forthe INETD scheme (29) 0.7811). The INETD method increases the accuracyin time from the second order of the original NETD method to the thirdorder, while keeping the spatial accuracy the same as the original NETDmethod. This conclusion is also verified by our numerical experiments ofcomputing 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 lessthan that of the NETD method.
  86. When using the same temporal and spatial increments, thecomputational costs and storage space of the INETD (20) and the NETDmethods are the same, while the INETD scheme (29) increases thecomputational cost by about 30% compared with the NETD scheme.
  87. However, because the INETD scheme (29) has a much higher stability limit,higher order of temporal accuracy and lower dispersion than the NETD, itcan lead to high computational efficiency.
  88. Meanwhile, wavefield modelling illustrates that the improved methodcan effectively reduce numerical dispersion when too-coarse computationgrids are used. The INETD method is related to the Lie group method.ACKNOWLEDGEMENTS
  89. We thank the editor and an anonymous reviewer for their detailedcomments and helpful revisions that greatly improved this paper. We alsothank Prof. Dinghui Yang for his fruitful suggestions that enabled us tohighly improve the paper.REFERENCES
  90. Aki, K. and Richards, P.G., 1980. Quantitative Seismology: Theory and Methods. W.H.Freeman and Co., San Francisco.
  91. Chen, S., Yang, D.H. and Deng, X.Y., 2010. A weighted Runge-Kutta method with weaknumerical dispersion for solving wave equations. Commun. Comput. Phys., 7:1027-1048.
  92. Dablain, M.A., 1986. The application of high-order differencing to the scalar waveequation. Geophysics, 51: 54-66.
  93. Fei, T. and Larner, K., 1995. Elimination of numerical dispersion in finite differencemodelling and migration by flux-corrected transport. Geophysics, 60: 1830-1842.
  94. Feng, L.L., Yang, D.H. and Xie, W., 2015. An efficient symplectic reverse timemigration method using a local nearly-analytic discrete operator in acoustictransversely isotropic media with a vertical symmetry axis. Geophysics, 80(2):3-S112.
  95. Fornberg, B., 1998. A Practical Guide to Psuedospectral Methods. Cambridge UniversityPress, Cambridge.
  96. Hairer, E., Lubich, C.H. and Wanner, G., 2006. Geometric Numerical Integration. 2'4 ed.,Springer Verlag, Berlin.
  97. Huang, B.S., 1992. A program for two-dimensional seismic wave propagation by thepseudospectrum method. Comput. Geosci., 18 : 289-307.
  98. Iserles, A., Munthe-Kaas, H., Norsett, S.P. and Zanna, A., 2000. Lie group methods. ActaNumer., 9: 215-365.
  99. Kelly, K., Ward, R., Treitel, S. and Alford, R., 1976. Synthetic seismograms: afinite-difference approach. Geophysics, 41: 2-27.
  100. Komatitsch and Tromp, J., 2003. A perfectly matched layer absorbing boundarycondition for the second-order seismic wave equation. Geophys. J. Internat., 154:146-153.
  101. Kosloff, D. and Baysal, E., 1982. Forward modeling by a Fourier method. Geophysics,47: 1402-1412.
  102. Kosloff, D., Reshef, M. and Loewenthal, D., 1984. Elastic wave calculations by the
  103. Fourier method. Bull. Seismol. Soc. Am., 74 : 875-891.
  104. Li, J.S., Yang, D.H., Wu, H. and Ma, X., 2017. A low-dispersive method using thehigh-order stereo-modeling operator for solving 2D wave equations. Geophys. J.Internat., 210: 1938-1964.
  105. Li, J.S., Yang, D.H. and Liu, F.Q., 2013. An efficient reverse-time migration methodusing local nearly analytic discrete operator. Geophysics, 78(1): -S23.
  106. Liu, S.L., Yang, D.H., Lang, C., Wang, W.S. and Pan, Z.D., 2017. Modified symplecticschemes with nearly-analytic discrete operators for acoustic wave simulations.Comput. Phys. Communicat., 213: 52-63.
  107. Lioyd, N., 1996. Finite Difference and Spectral Methods for Ordinary and Partial
  108. Differential Equations. Ph.D. thesis, Cornell University, Ithaca.
  109. Ma, X., Yang, D.H. and Liu, F., 2011. A nearly analytic symplectically partitioned
  110. Runge-Kutta method for 2-D seismic wave equations. Geophys. J. Internat., 187:480-496.
  111. Ma, X., Yang, D.H., Song, G.J. and Wang, M.X., 2014. A low-dispersive symplecticpartitioned Runge-Kutta method for solving seismic wave equations - I. Schemeand theoretical analysis. Bull. Seismol. Soc. Am., 104: 2206-2225.
  112. Ma, X., Yang, D.H. and Song, G.J., 2015. A low-dispersive symplectic partitioned
  113. Runge-Kutta method for solving seismic wave equations - II. Wave-fieldsimulations. Bull. Seismol. Soc. Am., 105: 657-675.
  114. Mizutani, H., Geller, R.J. and Takeuchi, N., 2000. Comparison of accuracy and efficiencyof time domain schemes for calculating synthetic seismograms. Phys. Earth Planet.Inter., 119 : 75-97.
  115. Moszo, P., Kristek, J. and Halada, L., 2000. 3-D fourth-order staggered-gridfinite-difference schemes: stability and grid dispersion. Bull. Seismol. Soc. Am., 90:587-603.
  116. Munthe-Kaas, H., 1999. High order Runge-Kutta methods on manifolds. Appl. Numer.Math., 29: 115-227.
  117. Robert, L.B., 1991. An Introduction to Lie Groups and Symplectic Geometry. Ph.D.thesis,Duke University, Park City.
  118. Tong, P., Yang D.H., Hua, B.L. and Wang, M.X., 2013. A high-order stereo-modelingmethod for solving wave equations. Bull. Seismol. Soc. Am., 103: 811-833.
  119. Versteeg, R.J. and Grau, G., 1991. The Marmousi experience. Proc. EAGE workshop
  120. Practical Aspects of Seismic Data Inversion, Copenhagen, 1990.
  121. Wang, S.Q., Yang, D.H. and Yang, K.D., 2002. Compact finite difference scheme forelastic equations. J. Tsinghua Univ., 42: 1128-1131 (in Chinese).
  122. Xie, W., Yang, D.H., Liu, F.Q. and Li, J.S., 2014. Reverse-time migration in acoustic
  123. VTI media using a high-order stereo operator. Geophysics, 79(3): WA3-WA11.
  124. Yang, D.H., Liu, E.R., Zhang, Z.J. and Teng, J.W., 2002. Finite-difference modelling in2D anisotropic media using a flux corrected transport technique. Geophys. J.Internat., 148: 320-328.
  125. Yang, D.H., Lu, M., Wu, R.S. and Peng, J. M., 2004. An optimal nearly analytic discretemethod for 2D acoustic and elastic wave equations. Bull. Seismol. Soc. Am., 94:1982-1991.
  126. Yang, D.H., Tong, P. and Deng, X.Y., 2012. A central difference method with lownumerical dispersion for solving the scalar wave equation. Geophys. Prosp., 60:885-905.
  127. Yang, D.H., Teng, J.W., Zhang, Z.J. and Liu, E., 2003. A nearly-analytic discrete methodfor acoustic and elastic wave equations in anisotropic media. Bull. Seismol. Soc.Am., 93: 882-890.
  128. Yang, D.H., Peng, J.M., Lu, M. and Terlaky, T., 2006. Optimal nearly analytic discreteapproximation to the scalar wave equation. Bull. Seismol. Soc. Am., 96:1114-1130.
  129. Yang, D.H., Song, G.J. and Lu, M., 2007. Optimally accurate nearly-analytic discretescheme for wave-field simulation in 3D anisotropic media. Bull. Seismol. Soc. Am.,97: 1557-1569.
  130. Yang, D.H., Wang, N., Chen, S. and Song, G.J., 2009. An explicit method based on theimplicit Runge- -Kutta algorithm for solving the wave equations. Bull. Seismol. Soc.Am., 99: 3340-3354.
  131. Yang, DH. and Wang, L., 2010a. A split-step algorithm with effectively suppressing thenumerical dispersion for 3D seismic propagation modeling. Bull. Seismol. Soc.Am., 100: 1470-1484.
  132. Yang, D.H., Song, G.J., Hua, B.L. and Henri, C., 2010b. Simulation of acousticwave-fields in heterogeneous media: A robust method with automaticallysuppressing numerical dispersion for large grid steps. Geophysics, 75(3):T99-T1
Share
Back to top
Journal of Seismic Exploration, Electronic ISSN: 0963-0651 Print ISSN: 0963-0651, Published by AccScience Publishing
We use cookies on our website to ensure you get the best experience.
Read more about our cookies here
Accept