# High-order dynamic lattice method for seismic simulation in anisotropic media

High-order dynamic lattice method for seismic simulation in anisotropic media Summary The discrete particle-based dynamic lattice method (DLM) offers an approach to simulate elastic wave propagation in anisotropic media by calculating the anisotropic micromechanical interactions between these particles based on the directions of the bonds that connect them in the lattice. To build such a lattice, the media are discretized into particles. This discretization inevitably leads to numerical dispersion. The basic lattice unit used in the original DLM only includes interactions between the central particle and its nearest neighbours; therefore, it represents the first-order form of a particle lattice. The first-order lattice suffers from numerical dispersion compared with other numerical methods, such as high-order finite-difference methods, in terms of seismic wave simulation. Due to its unique way of discretizing the media, the particle-based DLM no longer solves elastic wave equations; this means that one cannot build a high-order DLM by simply creating a high-order discrete operator to better approximate a partial derivative operator. To build a high-order DLM, we carry out a thorough dispersion analysis of the method and discover that by adding more neighbouring particles into the lattice unit, the DLM will yield different spatial accuracy. According to the dispersion analysis, the high-order DLM presented here can adapt the requirement of spatial accuracy for seismic wave simulations. For any given spatial accuracy, we can design a corresponding high-order lattice unit to satisfy the accuracy requirement. Numerical tests show that the high-order DLM improves the accuracy of elastic wave simulation in anisotropic media. Numerical modelling, Computational seismology, Seismic anisotropy, Wave propagation INTRODUCTION In most cases, the seismic wave phenomenon is interpreted using the elastic wave equation, which is based on the continuum theory. Analytical solutions to the equation only apply to simple media (Aki & Richards 2002). Therefore, a variety of numerical methods, such as the finite-difference method (FDM, e.g. Kelly et al. 1976; Madariaga 1976; Virieux 1984, 1986; Saenger et al. 2000), have been developed to yield approximate solutions to the wave equation. In regard to the numerical simulation of seismic wave propagation, due to the limitations of digital computers, certain discretization in either the space domain or the wavenumber domain is required. However, wave equation-based numerical methods are governed by a spatially continuous differential equation. Special calculations of particular boundary conditions at discontinuities are required. As a result, although the model and wavefields are already discretized into discrete gridpoints and are no longer continuous, wave equation-based methods still need to smooth the model to satisfy the continuity requirement of the wave equation; additionally, it is difficult for these methods to incorporate the naturally discontinuous features of geological materials, such as fractures and faults. Meanwhile, another category of numerical methods has been quite successful in simulating dynamic fractures and fault propagation. Cundall & Strack (1979) developed one of the first discrete particle methods, which was called the ‘distinct element method’, by extending the solid-state physics models of crystalline materials (Hoover et al. 1974) to study the static deformation and fracturing of granular materials. This methodology has since been used to investigate earthquake dynamics (Mora & Place 1998) and tectonic processes (Saltzer & Pollar 1992). Instead of directly solving the partial differential equation, these discrete methods treat geological materials as assemblages of numerous interconnected discrete particles and focus on the ‘microscopic’ interactions of these particles. Therefore, these methods eliminate the continuity restriction of the wave equation, offering a possible way to incorporate geological discontinuity. For this reason, particle-based methods are of great flexibility and can be used to solve different geological problems. Toomey & Bean (2000) discussed a 2-D ‘discrete particle scheme’ for seismic wave simulation in isotropic media. Their method assumes the medium to comprise discrete particles that are closely packed in a hexagonal lattice. In this lattice, particles interact with their bonded neighbours through Hooke's law. Toomey et al. (2002) also used this discrete particle scheme to investigate fracture properties from seismic data. Further studies on this issue were carried out by Lubbe & Worthington (2006) and Möllhoff & Bean (2008) using the same scheme. This method only takes central forces into account; thus, the Poisson ratio of the media is fixed at 0.25. Del Valle-García & Sánchez-Sesma (2003) improved Toomey's method by including a bond-bending force term to discard this Poisson ratio restriction and implemented it to simulate Rayleigh waves on the free surface. O'Brien & Bean (2004) extended this elastic lattice method to a 3-D case by using a 3-D lattice model. Their method has since been employed to study the structural topography of volcanoes (O'Brien & Bean 2009). Hu & Jia (2016) discussed a particle-based dynamic lattice method (DLM) that can be used to simulate seismic wave propagation in transversely isotropic (TI) media. On the other hand, due to their use of discretization, numerical methods face spatial dispersion that causes error in seismic wave simulation. To achieve certain accuracy, these methods must increase the density of the gridpoints, which will significantly increase memory consumption and slow down the simulation. To improve the performance of wave simulation, numerical methods such as the FDM have developed corresponding high-order schemes (e.g. Levander 1988; Hixon 2000; Moczo et al. 2002; Saenger & Bohlen 2004) that can reduce this spatial dispersion. Originally, the DLM only considers the interactions between the nearest particles, which results in great spatial dispersion compared with other high-order numerical methods. Therefore, reducing the spatial dispersion of the DLM has become an important issue. Meanwhile, for the wave equation-based methods, such as the FDM, one can easily use high-order finite-difference operators to better approximate the partial derivative operator. For the particle-based DLM, however, this is not the case because the method does not solve the wave equation directly. Therefore, the challenge lies in building high-order particle-based methods. By retaining the nearest-particle restriction and including an additional local interaction force acting on individual particles, O'Brien et al. (2009) managed to decrease the numerical dispersion of their elastic lattice method for isotropic media. However, in the case of anisotropy, which is represented by the anisotropic interactions according to dynamic lattice theory, it is difficult to suppress numerical dispersion by simply adding such an interaction force without considering the micromechanical interactions between particles other than those between the closest neighbours. The particle-based irregular lattice method (O'Brien & Bean 2011), which takes more than the nearest neighbours into account, is brought up to simulate elastic wave propagation. Although the implementation of the irregular lattice removes the anisotropic nature of fracture propagation, it also results in more dispersion than using the regular lattice. To fully analyse the spatial dispersion of the DLM, we take the different interaction patterns of the particle lattice into account and establish the high-order DLM. According to our analysis, by including more particles into the lattice unit, this method will produce less spatial dispersion, as the magnitude of the dispersion is related to the number of the added particles in the lattice unit. Using the high-order DLM, we can achieve higher numerical accuracy with relatively larger lattice spacing, thereby reducing the memory consumption and improving the CPU time efficiency of elastic wave simulation in anisotropic media. BASIC LATTICE UNITS IN THE DLM The most prominent feature of the DLM is its unique method of discretizing the media into a particle lattice. Thus, the layout of the lattice unit plays an important role in this method (Monette & Anderson 1994). Here, the three different layouts of the basic lattice unit shown in Fig. 1 are investigated to obtain insight into this method and lay the foundation of the high-order DLM. In such a basic lattice unit, the central particle 0 is bonded with its neighbours. As Fig. 2 shows, particle interactions in the basic lattice unit are achieved through the stretching and bending of these bonds. The stiffness of the bonds is calibrated by the coefficients ki and bij when the bonds are stretched and bent, respectively. The bigger these coefficients are, the harder it is to stretch or bend these bonds. In the case of anisotropy, because these stiffness coefficients determine the elastic interactions between these particles, thereby defining the elastic properties of the particle lattice, we create an anisotropic lattice by setting the stiffness coefficients in different directions. According to Monette & Anderson (1994), the elastic energy in the basic square lattice unit is given by   $$E = \frac{1}{2}{\sum\limits_{i = 1}^{{n \mathord{\left/ {\vphantom {n 2}} \right.} 2}} {{k_i}\left( {\left| {{{{\bf r}}_{0i}}} \right| - \left| {{{\bf r}}_{0i}^0} \right|} \right)} ^2} + \frac{1}{2}{\sum\limits_{i = 1,j ={\bmod}\!(i,n) + 1}^n {{b_{ij}}\left( {\cos {\theta _{ij}} - \cos {\theta _0}} \right)} ^2},$$ (1)in which E is the elastic energy, and r0i and $${{\bf r}}_{0i}^0$$ are the vectors of the bond that connects particle 0 and i (bond 0–i for short) in the distorted and undistorted lattices, respectively. Similarly, θij and θ0 are the angles between bonds 0–i and 0–j in the distorted and undistorted lattices, respectively. mod(i, n) denotes the remainder function after i is divided by the amount of lattice bonds in the basic lattice unit, that is, n. Figure 1. View largeDownload slide Basic lattice units. (a) Crisscross-shaped unit, (b) X-shaped unit and (c) basic square unit. The central particle 0 is bonded with its neighbours. The coefficients ki and bij represent the stiffness of the bonds in terms of stretching and bending, respectively. a is the lattice spacing that denotes the spatial interval between adjacent particles. Figure 1. View largeDownload slide Basic lattice units. (a) Crisscross-shaped unit, (b) X-shaped unit and (c) basic square unit. The central particle 0 is bonded with its neighbours. The coefficients ki and bij represent the stiffness of the bonds in terms of stretching and bending, respectively. a is the lattice spacing that denotes the spatial interval between adjacent particles. Figure 2. View largeDownload slide Particle interactions through the stretching and bending of bonds. (a) It shows the stretching of the bond that connects particles 0 and i. The coefficient ki represents the stiffness of the bond in terms of stretching. (b) It demonstrates the bending of the bonds that connect particles i, j and 0 with particle 0 as the pivot. The bond-bending coefficient bij denotes the stiffness of these bonds in terms of bending. Figure 2. View largeDownload slide Particle interactions through the stretching and bending of bonds. (a) It shows the stretching of the bond that connects particles 0 and i. The coefficient ki represents the stiffness of the bond in terms of stretching. (b) It demonstrates the bending of the bonds that connect particles i, j and 0 with particle 0 as the pivot. The bond-bending coefficient bij denotes the stiffness of these bonds in terms of bending. Once the lattice unit is chosen, it is necessary to find the relationship between the micromechanical stiffness coefficients and the macroscopic parameters. According to Hu & Jia (2016), we can find the transformation by comparing the elastic energy in a lattice unit given by eq. (1) using the strain energy density φ of a 2-D TI medium with a vertical axis of symmetry (VTI medium) given by   $$\varphi = \frac{1}{2}{c_{11}}{\varepsilon _{xx}}^2 + \frac{1}{2}{c_{33}}{\varepsilon _{zz}}^2 + 2{c_{44}}{\varepsilon _{xz}}^2 + {c_{13}}{\varepsilon _{xx}}{\varepsilon _{zz}},$$ (2)where εxx, εzz and εxz denote the corresponding components of the strain tensor. c11, c33,c13 and c44 represent the corresponding elements of the elasticity tensor of the VTI medium. According to eqs (1) and (2), for the crisscross-shaped lattice unit in Fig. 1(a), where the horizontal interval between neighbouring particles is equal to the vertical interval, we obtain the relationship between the stiffness coefficients and the elasticity tensor for VTI media given by (see Appendix  A)   $${c_{13}} = 0,$$ (3a)  $${k_1} = {c_{11}},$$ (3b)  $${k_2} = {c_{33}},$$ (3c)and   $${b_{12}} = {b_{23}} = {b_{34}} = {b_{41}} = \frac{1}{4}{a^2}{c_{44}},$$ (3d)where a is the lattice spacing, indicating the interval between adjacent particles. For the x-shaped lattice unit in Fig. 1(b), this relationship has the form of   $${c_{11}} = {c_{33}},$$ (4a)  $${k_1} = {k_2} = 2{c_{44}},$$ (4b)  $${c_{13}} = 2{c_{44}} - {c_{11}}$$ (4c)and   $${b_{12}} = {b_{23}} = {b_{34}} = {b_{41}} = \frac{{{a^2}}}{4}\left( {{c_{11}} - {c_{44}}} \right).$$ (4d) According to eqs (3) and (4), both the crisscross- and the x-shaped lattice units are strongly restricted in terms of their elastic properties, that is, both units only represent some particular VTI media. The reason for this restriction is that the crisscross- and x-shaped lattice units only take four of the neighbouring particles into account. The interactions between these particles are thus insufficient to simulate all kinds of VTI media. The basic square unit plotted in Fig. 1(c) is fully discussed by Hu & Jia (2016); when b12 = b34 = b56 = b78 and b23 = b45 = b67 = b81, the relation between the stiffness coefficients of the square lattice and the elasticity tensor is given by   $${{\bf W}}\left( \theta \right)\left( {\begin{array}{c} {{k_1}}\\ {{k_2}}\\ {{k_3}}\\ {{k_4}}\\ {{b_{12}}}\\ {{b_{23}}} \end{array}} \right) = \left( {\begin{array}{c} {{c_{11}}}\\ {{c_{33}}}\\ {{c_{44}}}\\ {{c_{13}}}\\ 0\\ 0 \end{array}} \right),$$ (5a)where c11, c33, c13 and c44 are the elements in the elasticity tensor of the corresponding TI medium when the symmetry axis is in the vertical direction, θ is the tilted angle for the tilted TI medium (TTI medium) and W(θ) is a 6 × 6 matrix that can be written as   $${{\bf W}}\left( \theta \right) = \left( {\begin{array}{c{@}{\quad}c{@}{\quad}c{@}{\quad}c{@}{\quad}c{@}{\quad}c} {{{\cos }^4}\theta }&{2{{\cos }^4}\left( {\frac{\pi }{4} + \theta } \right)}&{{{\sin }^4}\theta }&{2{{\cos }^4}\left( {\frac{\pi }{4} - \theta } \right)}&{\frac{{{{\left( {\cos 2\theta - \sin 2\theta } \right)}^2}}}{{2{a^2}}}}&{\frac{{{{\left( {\cos 2\theta + \sin 2\theta } \right)}^2}}}{{2{a^2}}}}\\ {{{\sin }^4}\theta }&{2{{\sin }^4}\left( {\frac{\pi }{4} + \theta } \right)}&{{{\cos }^4}\theta }&{2{{\sin }^4}\left( {\frac{\pi }{4} - \theta } \right)}&{\frac{{{{\left( {\cos 2\theta - \sin 2\theta } \right)}^2}}}{{2{a^2}}}}&{\frac{{{{\left( {\cos 2\theta + \sin 2\theta } \right)}^2}}}{{2{a^2}}}}\\ {\frac{{{{\sin }^2}2\theta }}{4}}&{\frac{{{{\cos }^2}2\theta }}{2}}&{\frac{{{{\sin }^2}2\theta }}{4}}&{\frac{{{{\cos }^2}2\theta }}{2}}&{\frac{{{{\left( {\cos 2\theta + \sin 2\theta } \right)}^2}}}{{2{a^2}}}}&{\frac{{{{\left( {\cos 2\theta - \sin 2\theta } \right)}^2}}}{{2{a^2}}}}\\ {\frac{{{{\sin }^2}2\theta }}{4}}&{\frac{{{{\cos }^2}2\theta }}{2}}&{\frac{{{{\sin }^2}2\theta }}{4}}&{\frac{{{{\cos }^2}2\theta }}{2}}&{ - \frac{{{{\left( {\cos 2\theta - \sin 2\theta } \right)}^2}}}{{2{a^2}}}}&{ - \frac{{{{\left( {\cos 2\theta + \sin 2\theta } \right)}^2}}}{{2{a^2}}}}\\ {{{\cos }^2}\theta \sin 2\theta }&{2{{\cos }^2}\left( {\frac{\pi }{4} + \theta } \right)\cos 2\theta }&{ - \sin 2\theta {{\sin }^2}\theta }&{ - 2{{\cos }^2}\left( {\frac{\pi }{4} - \theta } \right)\cos 2\theta }&{\frac{{\cos 4\theta }}{{{a^2}}}}&{ - \frac{{\cos 4\theta }}{{{a^2}}}}\\ {\sin 2\theta \sin^2\theta }&{2{{\sin }^2}\left( {\frac{\pi }{4} + \theta } \right)\cos 2\theta }&{ - \sin 2\theta {{\cos }^2}\theta }&{ - 2{{\sin }^2}\left( {\frac{\pi }{4} - \theta } \right)\cos 2\theta }&{ - \frac{{\cos 4\theta }}{{{a^2}}}}&{\frac{{\cos 4\theta }}{{{a^2}}}} \end{array}} \right).$$ (5b) Since the determinant of this matrix is non-zero, we can always find the corresponding stiffness coefficients for any TI medium as long as the medium can be described by the elasticity tensor. Therefore, considering the limitations of the crisscross- and the x-shaped lattice units in representing realistic media, we build the high-order DLM based on the more applicable basic square lattice unit. HIGH-ORDER LATTICE UNITS Due to spatial discretization, the original DLM suffers from spatial dispersion in terms of seismic wave simulation. For the basic square lattice unit, micromechanical interactions only exist between the nearest pairs of particles. As a result, the dispersion of the DLM cannot be ignored as it can be in other high-order wave equation-based numerical methods. To suppress this spatial dispersion, we develop a high-order dynamic lattice that can decrease the dispersion based on the accuracy requirement. According to the theory of the DLM (Hu & Jia 2016), the interaction force acting upon an individual particle is calculated by Lagrange's equation. Hu & Jia (2016) simplified the calculation for this interaction force by linearizing Lagrange's equation. In the case of homogenous media, the interaction force on the central particle in the lattice unit can be written as a linear combination of the differences in the displacements between the central particle and its neighbours given by   $$\left( {\begin{array}{c} {f_0^x}\\ {f_0^z} \end{array}} \right) = \left( {\begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{{\bf A}}_{11}}}&{{{{\bf A}}_{12}}}&{{{{\bf A}}_{13}}}&{{{{\bf A}}_{14}}}\\ {{{{\bf A}}_{21}}}&{{{{\bf A}}_{22}}}&{{{{\bf A}}_{23}}}&{{{{\bf A}}_{24}}} \end{array}} \right) \times {\Big( {\begin{array}{*{20}{c}} {u_{01}^x}&\quad \cdots &\quad{u_{08}^x}&\quad{u_{01}^z}&\quad \cdots &\quad{u_{08}^z} \end{array}} \Big)^{\rm{T}}},$$ (6)where $$f_0^x$$ and $$f_0^z$$ are the x- and z-components of the interaction force acting upon the central particle that is, particle 0 in the lattice unit (Fig. 1c), respectively; and $$u_{0i}^x$$ and $$u_{0i}^z$$ are the x- and z-components of the difference in the displacement between particles 0 and i in the lattice unit, respectively, that is,   $$\left( {\begin{array}{c} {u_{0i}^x}\\ {u_{0i}^z} \end{array}} \right) = \left( {\begin{array}{c} {u_i^x - u_0^x}\\ {u_i^z - u_0^z} \end{array}} \right),$$ (7)where $$u_i^x$$ and $$u_i^z$$ are the x- and z-components of the displacement for particle i, respectively, whereas $$u_0^x$$ and $$u_0^z$$ are the x- and z-components of the displacement for particle 0, respectively. The coefficients for this linear combination are given by   $${{{\bf A}}_{11}} = {{{\bf A}}_{12}} = \left( {\begin{array}{c} {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{2{a^2}}} + {k_1}}\quad{ - \displaystyle\frac{{{b_{23}} - {b_{12}}}}{{4{a^2}}} + \frac{1}{2}{k_2}}\quad{\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{2{a^2}}}}\quad{ - \displaystyle\frac{{{b_{12}} - {b_{23}}}}{{4{a^2}}} + \displaystyle\frac{1}{2}{k_4}} \end{array}} \right),$$ (8a)  $${{{\bf A}}_{13}} = {{{\bf A}}_{14}} = {{{\bf A}}_{21}} = {{{\bf A}}_{22}} = \left( {\begin{array}{*{20}{c}} { - \displaystyle\frac{{{b_{23}} - {b_{12}}}}{{2{a^2}}}}&\quad{\displaystyle\frac{1}{2}{k_2}}&\quad{ - \displaystyle\frac{{{b_{12}} - {b_{23}}}}{{2{a^2}}}}&\quad { - \displaystyle\frac{1}{2}{k_4}} \end{array}} \right),$$ (8b)and   $${{{\bf A}}_{23}} = {{{\bf A}}_{24}} = \left( {\begin{array}{c} {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{2{a^2}}}}\quad { - \displaystyle\frac{{{b_{12}} - {b_{23}}}}{{4{a^2}}} + \displaystyle\frac{1}{2}{k_2}}\quad{\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{2{a^2}}} + {k_3}}\quad{ - \displaystyle\frac{{{b_{23}} - {b_{12}}}}{{4{a^2}}} + \displaystyle\frac{1}{2}{k_4}} \end{array}} \right).$$ (8c) Using eq. (6), one can obtain the interaction forces acting upon the particles as long as the displacements of these particles are provided. Using the initial displacements and the interaction forces of the particles, the displacements for the next time step are obtained through the Verlet accumulation, which is given by   $$\left( {\begin{array}{c} {u_i^x\left( {{t_0} + \Delta t} \right)}\\ {u_i^z\left( {{t_0} + \Delta t} \right)} \end{array}} \right) = \left( {\begin{array}{c}{\ddot{u}_i^x\left( {{t_0}} \right)}\\ {\ddot{u}_i^z\left( {{t_0}} \right)} \end{array}} \right)\Delta {t^2} + 2\left( {\begin{array}{c} {u_i^x\left( {{t_0}} \right)}\\ {u_i^z\left( {{t_0}} \right)} \end{array}} \right) - \left( {\begin{array}{c} {u_i^x\left( {{t_0} - \Delta t} \right)}\\ {u_i^z\left( {{t_0} - \Delta t} \right)} \end{array}} \right)$$ (9)where Δt is the time interval, $$u_i^x( t )$$ and $$u_i^z( t )$$ are the x- and z-components of the displacement for particle i at time t, respectively and $$\ddot{u}_i^x( t )$$ and $$\ddot{u}_i^z( t )$$ denote the x- and z-components of the acceleration for particle i, respectively, which can be written as   $$\left( {\begin{array}{c} {\ddot{u}_i^x\left( t \right)}\\ {\ddot{u}_i^z\left( t \right)} \end{array}} \right) = \frac{1}{{\rho {a^2}}}\left( {\begin{array}{c}{f_i^x}\\ {f_i^z} \end{array}} \right),$$ (10)where ρ is the density of the media and a is the lattice spacing. Based on this linearized interaction force and the Verlet accumulation, we carry out a plane wave analysis of the method according to Saenger & Bohlen (2004) to investigate the spatial dispersion of the basic square lattice unit for wave propagation. The plane wave that propagates in the lattice unit takes the form of   $$\left( {\begin{array}{c} {u_i^x\left( t \right)}\\ {u_i^z\left( t \right)} \end{array}} \right) = \left( {\begin{array}{c} {{U_x}\exp \left( {j{{\bf k}} \cdot {{\bf r}}_{0i}^0 - j\omega t} \right)}\\ {{U_z}\exp \left( {j{{\bf k}} \cdot {{\bf r}}_{0i}^0 - j\omega t} \right)} \end{array}} \right),$$ (11)where Ux and Uz are the x- and z-components of the amplitude of the plane wave, respectively; $${{\bf r}}_{0i}^0$$ is the vector for the bond that connects particles 0 and i in the undistorted lattice unit; k is the vector wavenumber; ω represents the angular frequency of the plane wave and $$j = \sqrt { - 1}$$. Substituting eq. (11) into eqs (10) and (6) yields   $$\left( {\begin{array}{c{@}{\quad }c} {{g_{11}}}&{{g_{12}}}\\ {{g_{21}}}&{{g_{22}}} \end{array}} \right)\left( {\begin{array}{c} {{U_x}}\\ {{U_z}} \end{array}} \right) = - {\omega ^2}\rho \left( {\begin{array}{c} {{U_x}}\\ {{U_z}} \end{array}} \right)$$ (12)in which   $${g_{11}} = {{{\bf A}}_{11}} \times {{\bf L}},$$ (13a)  $${g_{12}} = {g_{21}} = {{{\bf A}}_{12}} \times {{\bf L}},$$ (13b)and   $${g_{22}} = {{{\bf A}}_{23}} \times {{\bf L}}.$$ (13c) For the basic square lattice unit, L is a function of ak given by   $${{\bf L}}\left( {a{{\bf k}}} \right) = \left( {\begin{array}{c} {\displaystyle\frac{{\cos \left( {a{{\bf k}} \cdot {{{\bf e}}_{01}}} \right) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right),$$ (14)in which $${{{\bf e}}_{0i}} = \frac{{{{\bf r}}_{0i}^0}}{{| {{{\bf r}}_{0i}^0} |}}$$. Solving eq. (12) yields   $$v\left( {a{{\bf k}}} \right) = \sqrt { - \frac{1}{\rho }\left( {{g_{11}} + {g_{22}} \pm \sqrt {{{\left( {{g_{11}} - {g_{22}}} \right)}^2} + 4{g_{12}}{g_{21}}} } \right)} ,$$ (15)where v(ak) represents the phase velocity for the plane wave that propagates in the basic square lattice as a function of lattice spacing a and the vector wavenumber k, and ‘ ± ’ denotes two different wave modes, that is, P and S waves. When |ak| approximates zero, which implies that the lattice spacing is much less than the wavelength, the phase velocity given by eq. (15) tends to represent the theoretical phase velocity vtheo of the corresponding continua, that is,   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} v\left( {a{{\bf k}}} \right) = {v_{{\rm{theo}}}}.$$ (16) The dispersion comes from the difference between v(ak) and vtheo. According to eqs (13) and 15, v(ak) can also be regarded as a function of L; thus, eq. (16) can be expressed as   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} v\left( {a{{\bf k}}} \right) = \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} f\left( {{\bf L}} \right) = f ( {\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}}\left( {a{{\bf k}}} \right)} ) = {v_{{\rm{theo}}}},$$ (17)where f(L) represents the phase velocity as a function of L. Therefore, the dispersion is in fact caused by the difference between L(ak) and $$\mathop {\lim }\limits_{| {a{{\bf k}}} | \to 0} {{\bf L}}( {a{{\bf k}}} )$$. When |ak| approximates 0, we have   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}} = \left( {\begin{array}{c} { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _1}}\\ { - {{\cos }^2}{\theta _2}}\\ { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _3}}\\ { - {{\cos }^2}{\theta _4}} \end{array}} \right),$$ (18)where θi denotes the angle between k and e0i. The difference between L(ak) and $$\mathop {\lim }\limits_{| {a{{\bf k}}} | \to 0} {{\bf L}}( {a{{\bf k}}} )$$ can be written as   $$d\left( {a{{\bf k}}} \right) = {{\bf L}}\left( {a{{\bf k}}} \right) - \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}}\left( {a{{\bf k}}} \right) = \left( {\begin{array}{c} {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _1}}\\ { - {{\cos }^2}{\theta _2}}\\ { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _3}}\\ { - {{\cos }^2}{\theta _4}} \end{array}} \right).$$ (19) Carrying out Taylor's expansion on L(ak) at ak = 0 yields   $$d\left( {a{{\bf k}}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}} + \displaystyle\frac{1}{2}{{\cos }^2}{\theta _1}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}} + {{\cos }^2}{\theta _2}}\\ {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}} + \displaystyle\frac{1}{2}{{\cos }^2}{\theta _3}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}} + {{\cos }^2}{\theta _4}} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{1}{{24}}{{\cos }^4}{\theta _1}{a^2}{{{\bf k}}^2}}\\ {\displaystyle\frac{1}{6}{{\cos }^4}{\theta _2}{a^2}{{{\bf k}}^2}}\\ {\displaystyle\frac{1}{{24}}{{\cos }^4}{\theta _3}{a^2}{{{\bf k}}^2}}\\ {\displaystyle\frac{1}{6}{{\cos }^4}{\theta _4}{a^2}{{{\bf k}}^2}} \end{array}} \right) + {\rm{Order}}( {{{\left| {a{{\bf k}}} \right|}^4}} ).$$ (20) This indicates that the deviations in the phase velocities in the basic square lattice unit from the theoretical phase velocities are at Order(|ak|2). Fig. 3 shows the relative deviations of the phase velocities of both P and S waves in the basic square lattice unit from their theoretical phase velocities, which are given by   $${d_P} = \displaystyle\frac{{{v_P}\left( {a{{\bf k}}} \right) - v_P^{{\rm{theo}}}}}{{v_P^{{\rm{theo}}}}},$$ (21a)and   $${d_S} = \displaystyle\frac{{{v_S}( {a{{\bf k}}} ) - v_S^{{\rm{theo}}}}}{{v_S^{{\rm{theo}}}}},$$ (21b)where dP and dS are the relative deviations of the phase velocities of P and S waves, respectively; vP(ak) and vS(ak) are the P- and S-wave phase velocities in the basic lattice unit, respectively and $$v_P^{{\rm{theo}}}$$ and $$v_S^{{\rm{theo}}}$$ are the theoretical phase velocities of P and S waves, respectively, in the corresponding continua. To improve the accuracy of the DLM in seismic wave simulation, we need to reduce the difference between v(ak) and vtheo. According to eq. (19), in such a basic square lattice unit, we can only reduce the dispersion by using smaller lattice spacing; however, this will significantly increase the memory consumption as well as the CPU time. Through research on modifying the lattice unit, we find that the spatial accuracy benefits from adding more neighbour particles into the lattice unit. In an effort to simultaneously reduce the spatial dispersion in terms of all lattice spacing and speed up the convergence of the particle lattice model to the corresponding continua as the lattice spacing approximates zero, we extend the basic square lattice unit in all existing directions to form a second-order lattice unit, which is shown in Fig. 4(a). Due to this modification, the linear combination given by eq. (6) is reconstructed to adapt to the new lattice unit. The coefficients for the linear combination in eq. (6) are reassigned. The modified version of eq. (6) with respect to the extended second-order unit can thus be written as   \begin{eqnarray} \left( { \begin{array}{c} {f_0^x}\\ {f_0^z} \end{array}} \right)& = &{\left( { \begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{{\bf A}}_{11}}}&{{{{\bf A}}_{12}}}&{{{{\bf A}}_{13}}}&{{{{\bf A}}_{14}}}\\ {{{{\bf A}}_{21}}}&{{{{\bf A}}_{22}}}&{{{{\bf A}}_{23}}}&{{{{\bf A}}_{24}}} \end{array}} \right) \times {{{\bf H}}_1} \times \left( { \begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{\left( {u_{01}^x} \right)}_1}}& \cdots &{{{\left( {u_{08}^x} \right)}_1}}&{{{\left( {u_{01}^z} \right)}_1}}& \cdots &{{{\left( {u_{08}^z} \right)}_1}} \end{array}} \right)}\nonumber\\ &+& \left( {\begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{{\bf A}}_{11}}} &{{{{\bf A}}_{12}}} &{{{{\bf A}}_{13}}} &{{{{\bf A}}_{14}}}\\ {{{{\bf A}}_{21}}} &{{{{\bf A}}_{22}}} &{{{{\bf A}}_{23}}} &{{{{\bf A}}_{24}}} \end{array}} \right) \times {{{\bf H}}_2} \times {{\left( { \begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{\left( {u_{01}^x} \right)}_2}} &\cdots &{{{\left( {u_{08}^x} \right)}_2}}& {{{\left( {u_{01}^z} \right)}_2}} &\cdots &{{{\left( {u_{08}^z} \right)}_2}} \end{array}} \right)}^{\rm{T}}}, \end{eqnarray} (22)where the two factor matrices H1 and H2 that take the forms of H1 = h1I and H2 = h2I are assigned to the relative displacements between the central particle and its neighbouring particles in the inner layer, that is, $${( {u_{0i}^x} )_1}$$ and those between the central particle and the newly added particles, that is, $${( {u_{0i}^x} )_2}$$, respectively. To find the ideal factors h1 and h2, we analyse the spatial dispersion of the extended lattice unit. This numerical dispersion can be written as   $$v'\left( {a{{\bf k}}} \right) = \sqrt { - \displaystyle\frac{1}{\rho }\left( {{{g'}_{11}} + {{g'}_{22}} \pm \sqrt {{{\left( {{{g'}_{11}} - {{g'}_{22}}} \right)}^2} + 4{{g'}_{12}}{{g'}_{21}}} } \right)} ,$$ (23)where   $${g'_{11}} = {{{\bf A}}_{11}} \times {{\bf L'}},$$ (24a)  $${g'_{12}} = {g'_{21}} = {{{\bf A}}_{12}} \times {{\bf L'}},$$ (24b)and   $${g'_{22}} = {{{\bf A}}_{23}} \times {{\bf L'}}.$$ (24c) Figure 3. View largeDownload slide Relative errors of P- and S-wave phase velocities for the basic square lattice unit in terms of lattice spacing a and the wavenumber components kx and kz. Figure 3. View largeDownload slide Relative errors of P- and S-wave phase velocities for the basic square lattice unit in terms of lattice spacing a and the wavenumber components kx and kz. Figure 4. View largeDownload slide (a) A second-order lattice unit. (b) and (c) represent the relative errors of P- and S-wave phase velocities, respectively, in terms of lattice spacing a and the wavenumber components kx and kz. Figure 4. View largeDownload slide (a) A second-order lattice unit. (b) and (c) represent the relative errors of P- and S-wave phase velocities, respectively, in terms of lattice spacing a and the wavenumber components kx and kz. For this second-order lattice unit, we have   $${{\bf L'}} = {h_1}\left( {\begin{array}{c} {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right) + {h_2}\left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\cos ( {2a{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {2\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {2a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {2\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right).$$ (25) To maintain the elastic properties of the lattice unit, when |ak| approximates zero, the phase velocities of both P and S waves for the second-order lattice unit should be the same as those for the basic square unit. Therefore,   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} f( {{{\bf L'}}} ) = f\left( {\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L'}}( {a{{\bf k}}} )} \right) = \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} v'( {a{{\bf k}}} ) = \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} f( {{\bf L}} ) = f\left( {\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}}( {a{{\bf k}}} )} \right);$$ (26)and the restriction on L' is given by   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L'}} = \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}} = \left( {\begin{array}{@{}*{1}{c}@{}} { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _1}}\\ { - {{\cos }^2}{\theta _2}}\\ { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _3}}\\ { - {{\cos }^2}{\theta _4}} \end{array}} \right).$$ (27) Carrying out the Taylor expansion on L' at ak = 0 yields   $${{\bf L'}} = \left( {\begin{array}{@{}*{1}{c}@{}} { - \displaystyle\frac{{\left( {{h_1} + 4{h_2}} \right)}}{2}{{\cos }^2}{\theta _1}}\\ { - \left( {{h_1} + 4{h_2}} \right){{\cos }^2}{\theta _2}}\\ { - \displaystyle\frac{{\left( {{h_1} + 4{h_2}} \right)}}{2}{{\cos }^2}{\theta _3}}\\ { - \left( {{h_1} + 4{h_2}} \right){{\cos }^2}{\theta _4}} \end{array}} \right) + \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{{24}}{{\cos }^4}{\theta _1} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{6}{{\cos }^4}{\theta _2} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{{24}}{{\cos }^4}{\theta _3} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{6}{{\cos }^4}{\theta _4} \cdot {{\left| {a{{\bf k}}} \right|}^2}} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)} \end{array}} \right).$$ (28) Substituting eq. (28) into eq. (27) defines the restriction on h1 and h2 as   \begin{eqnarray}{h_1} + 4{h_2} &=& 1.\end{eqnarray} (29) The extended lattice unit will converge to the corresponding homogeneous continua under this restriction. However, it still leaves these two factors with one degree of freedom. On the other hand, similar to eq. (19), the dispersion for the extended lattice is also obtained from the difference between L' and $$\mathop {\lim }\limits_{| {a{{\bf k}}} | \to 0} {{\bf L}}$$, that is,   $$d' = {{\bf L'}} - \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}}.$$ (30) Substituting eqs (28) and (29) into eq. (30) yields   $$d' = \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{{24}}{{\cos }^4}{\theta _1} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{6}{{\cos }^4}{\theta _2} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{{24}}{{\cos }^4}{\theta _3} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{6}{{\cos }^4}{\theta _4} \cdot {{\left| {a{{\bf k}}} \right|}^2}} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)} \end{array}} \right).$$ (31) To reduce dispersion and speed up the convergence as |ak| approximates zero, the first term on the right-hand side of eq. (31) should be zero. This yields another restriction on h1 and h2 as   \begin{eqnarray}{h_1} + 16{h_2} &=& 0.\end{eqnarray} (32) Then, the magnitude of the dispersion term in eq. (31) is reduced to Order(|ak|4). Eqs (29) and (32) give the ideal relation between h1 and h2 that maintains the elastic properties of the lattice and reduces the numerical dispersion to Order(|ak|4). Fig. 4(b) shows the errors in phase velocities caused by the numerical dispersion in the extended second-order lattice unit. Applying this analysis to even higher order lattice units yields the general representation of the high-order DLM. As Fig. 5 illustrates, in the Nth-order lattice unit where the central particle is surrounded by N layers of particles, the interaction force acting on the central particles takes the form of   $$\left( {\begin{array}{@{}*{1}{c}@{}} {f_0^x}\\ {f_0^z} \end{array}} \right) = \sum\limits_{l = 1}^N {\left( {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} {{{{\bf A}}_{11}}}&{{{{\bf A}}_{12}}}&{{{{\bf A}}_{13}}}&{{{{\bf A}}_{14}}}\\ {{{{\bf A}}_{21}}}&{{{{\bf A}}_{22}}}&{{{{\bf A}}_{23}}}&{{{{\bf A}}_{24}}} \end{array}} \right) \times {{{\bf H}}_{\,l}} \times \left( {\begin{array}{*{20}{c}} {{{\left( {u_{01}^x} \right)}_l}}&\quad \cdots &\quad{{{\left( {u_{08}^x} \right)}_l}}&\quad{{{\left( {u_{01}^z} \right)}_l}}&\quad \cdots &\quad{{{\left( {u_{08}^z} \right)}_l}} \end{array}} \right)} ,$$ (33)where   $${{{\bf H}}_l} = {h_l}{{\bf I}},$$ (34) in which I is a 16 × 16 identical matrix; hl is the factor assigned to the lth layer and $${( {u_{0i}^x} )_l}$$ and $${( {u_{0i}^z} )_l}$$ represent the x- and z-components, respectively, of the difference in displacement between particle 0 and the corresponding particle in the lth layer that lies in the direction of r0i relative to particle 0. For this Nth-order lattice unit, we have   $${{{\bf L}}^{( N )}} = \sum\limits_{l = 1}^N {{h_l}\left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\cos ( {al{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 al{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {al{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 al{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right)} .$$ (35) Figure 5. View largeDownload slide The configuration of the Nth-order lattice unit. Figure 5. View largeDownload slide The configuration of the Nth-order lattice unit. Similar to the second-order lattice unit, for the Nth-order lattice unit, the dispersion also comes from the difference between L(N) and $$\mathop {\lim }\limits_{| {a{{\bf k}}} | \to 0} {{\bf L}}$$, that is, d(N). Substituting the Taylor expansion of L(N) at |ak| = 0 into this difference yields   $$\begin{array}{@{}c@{}} {d^{\left( N \right)}} = \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{1}{2}{{\cos }^2}{\theta _1} - \displaystyle\frac{{\sum\limits_{l = 1}^N {{l^2}{h_l}} }}{{2!}}{{\cos }^2}{\theta _1}}\\ {{{\cos }^2}{\theta _2} - \displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^2}{h_l}} }}{{2!}}{{\cos }^2}{\theta _2}}\\ {\displaystyle\frac{1}{2}{{\cos }^2}{\theta _3} - \displaystyle\frac{{\sum\limits_{l = 1}^N {{l^2}{h_l}} }}{{2!}}{{\cos }^2}{\theta _3}}\\ {{{\cos }^2}{\theta _4} - \displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^2}{h_l}} }}{{2!}}{{\cos }^2}{\theta _4}} \end{array}} \right) + \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\sum\limits_{l = 1}^N {{l^4}{h_l}} }}{{4!}}{{\cos }^4}{\theta _1} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^4}{h_l}} }}{{4!}}{{\cos }^4}{\theta _2} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\sum\limits_{l = 1}^N {{l^4}{h_l}} }}{{4!}}{{\cos }^4}{\theta _3} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^4}{h_l}} }}{{4!}}{{\cos }^4}{\theta _4} \cdot {{\left| {a{{\bf k}}} \right|}^2}} \end{array}} \right) \cdots \cdots - \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} }}{{\left( {2N} \right)!}}{{\cos }^{2N}}{\theta _1} \cdot {{\left| {a{{\bf k}}} \right|}^{2N - 2}}}\\ {\displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} }}{{\left( {2N} \right)!}}{{\cos }^{2N}}{\theta _2} \cdot {{\left| {a{{\bf k}}} \right|}^{2N - 2}}}\\ {\displaystyle\frac{{\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} }}{{\left( {2N} \right)!}}{{\cos }^{2N}}{\theta _3} \cdot {{\left| {a{{\bf k}}} \right|}^{2N - 2}}}\\ {\displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} }}{{\left( {2N} \right)!}}{{\cos }^{2N}}{\theta _4} \cdot {{\left| {a{{\bf k}}} \right|}^{2N - 2}}} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^{2N}}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^{2N}}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^{2N}}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^{2N}}} \right)} \end{array}} \right). \end{array}$$ (36) To reduce the dispersion and maintain the elastic properties of the Nth-order lattice unit, hl is given by   $$\left\{ {\begin{array}{@{}*{1}{c}@{}} {\sum\limits_{l = 1}^N {{l^2}{h_l}} = 1}\\ {\sum\limits_{l = 1}^N {{l^4}{h_l}} = 0}\\ {\sum\limits_{l = 1}^N {{l^6}{h_l}} = 0}\\ \vdots \\ {\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} = 0} \end{array}} \right..$$ (37) Under this relation, the dispersion is reduced to Order((ak)2N). This means that the high-order DLM using an Nth-order lattice unit will have 2Nth-order accuracy in terms of its lattice spacing. Based on the dispersion analysis, we can precisely choose the high-order lattice unit based on the accuracy requirement in terms of seismic wave simulation. Since the high-order DLM still focuses on the particle interactions that depend on the bonds that connect these particles, discontinuities in geological materials, such as fractures and faults, can also be interpreted by the breaking or weakening of these bonds. Therefore, the high-order DLM retains the advantage of using the discrete particle scheme in simulating complex geological materials, and it can be used to study both the static and dynamic deformations of the material. On the other hand, due to the particle feature of high-order DLMs, we form an absorbing pad by simply adding a friction force in the opposite direction of the particle velocity on each particle in the absorbing pad to gradually dampen their movements (see Appendix  B). In the aspect of 3-D seismic wave simulation, one can build a 3-D high-order lattice unit by combining three 2-D high-order lattice units in three orthogonal planes; then, the mathematical derivation of the 2-D high-order DLM can be easily implemented for the 3-D case. STABILITY CONDITION In regard to wave simulation in the time domain, due to the time iteration, the stability of the numerical simulation is also an important issue. We analyse the stability condition based on Von Neumann's method, which is described by Press et al. (1986). The numerical error is assumed as   $$\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _i^x\left( t \right)}\\ {\delta _i^z\left( t \right)} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {{U_x}\left( t \right)}\\ {{U_z}\left( t \right)} \end{array}} \right)\exp \left( {j{{\bf k}} \cdot {{\bf r}}_{0i}^0} \right),$$ (38)in which $$\delta _i^x( t )$$ and $$\delta _i^z( t )$$ are the x- and z-components of the numerical error for the displacement with respect to particle i, respectively. For this numerical error, the Verlet accumulation takes the form of   $$\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _i^x\left( {t + \Delta t} \right)}\\ {\delta _i^z\left( {t + \Delta t} \right)} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {\ddot{\delta }_i^x\left( t \right)}\\ {\ddot{\delta }_i^z\left( t \right)} \end{array}} \right)\Delta {t^2} + 2\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _i^x\left( t \right)}\\ {\delta _i^z\left( t \right)} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} {\delta _i^x\left( {t - \Delta t} \right)}\\ {\delta _i^z\left( {t - \Delta t} \right)} \end{array}} \right).$$ (39) According to eq. (33), we have   $$\left( {\begin{array}{c} {\ddot{\delta }_0^x\left( t \right)}\\ {\ddot{\delta }_0^z\left( t \right)} \end{array}} \right) = \displaystyle\frac{{{\bf G}}}{\rho }\left( {\begin{array}{c} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right),$$ (40a)where   $${{\bf G}} = \left( {\begin{array}{c{@}{\quad }c} {g_{11}^{\left( N \right)}}&{g_{12}^{\left( N \right)}}\\ {g_{21}^{\left( N \right)}}&{g_{22}^{\left( N \right)}} \end{array}} \right),$$ (40b)in which   $$g_{11}^{\left( N \right)} = {{{\bf A}}_{11}} \times {{{\bf L}}^{\left( N \right)}},$$ (40c)  $$g_{12}^{\left( N \right)} = g_{21}^{\left( N \right)} = {{{\bf A}}_{12}} \times {{{\bf L}}^{\left( N \right)}},$$ (40d)and   $$g_{22}^{\left( N \right)} = {{{\bf A}}_{23}} \times {{{\bf L}}^{\left( N \right)}}.$$ (40e) For the Nth-order lattice unit,   $${{{\bf L}}^{( N )}} = \sum\limits_{l = 1}^N {{h_l} \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\cos ( {la{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 la{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {la{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 la{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right)} .$$ (40f) Substituting eq. (40) into eq. (39), we have   $$\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( {t + \Delta t} \right)}\\ {\delta _0^z\left( {t + \Delta t} \right)} \end{array}} \right) - 2\left( {{{\bf I}} - \displaystyle\frac{{\Delta {t^2}}}{{2\rho }}{{\bf G}}} \right)\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right) + \left( {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( {t - \Delta t} \right)}\\ {\delta _0^z\left( {t - \Delta t} \right)} \end{array}} \right) = 0,$$ (41) We define the growth rate of the numerical error by   $$\gamma = \displaystyle\frac{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( {t + \Delta t} \right)}\\ {\delta _0^z\left( {t + \Delta t} \right)} \end{array}} \right\|}}{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right\|}} = \displaystyle\frac{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right\|}}{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( {t - \Delta t} \right)}\\ {\delta _0^z\left( {t - \Delta t} \right)} \end{array}} \right\|}}.$$ (42) Therefore, eq. (41) can be written as   $${\gamma ^2} - 2\eta \gamma + 1 = 0,$$ (43a)where   $$\eta = \displaystyle\frac{{\left\| {\left( {{{\bf I}}{\rm{ + }}\displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}{{\bf G}}} \right)\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right)} \right\|}}{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right\|}}.$$ (43b) The stability condition requires   $$\left| \gamma \right| \le 1.$$ (44)  This leads to the restriction on η given by   $$\left| \eta \right| \le 1.$$ (45)  Since   $$\left| \eta \right| \le \left| {1 + \displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}\mu } \right|,$$ (46)in which μ is the eigenvalue of matrix G, given by   $$\mu = \displaystyle\frac{1}{2}\left( {g_{11}^{\left( N \right)} + g_{22}^{\left( N \right)} \pm \sqrt {{{\left( {g_{11}^{\left( N \right)} - g_{22}^{\left( N \right)}} \right)}^2} + 4g_{12}^{\left( N \right)}g_{21}^{\left( N \right)}} } \right),$$ (47)eq. (45) can be rewritten as   $$- 2 \le \displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}\mu \le 0.$$ (48) The maximum values of μ are reached when   $${{\bf k}} = \left( {\begin{array}{@{}*{1}{c}@{}} {{k_x}}\\ {{k_z}} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {{\pi \mathord{\left/ {\vphantom {\pi a}} \right.} a}}\\ 0 \end{array}} \right),$$ (49a)  $${{\bf k}} = \left( {\begin{array}{@{}*{1}{c}@{}} {{k_x}}\\ {{k_z}} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} 0\\ {{\pi \mathord{\left/ {\vphantom {\pi a}} \right.} a}} \end{array}} \right),$$ (49b)or   $${{\bf k}} = \left( {\begin{array}{@{}*{1}{c}@{}} {{k_x}}\\ {{k_z}} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {{\pi \mathord{\left/ {\vphantom {\pi a}} \right.} a}}\\ {{\pi \mathord{\left/ {\vphantom {\pi a}} \right.} a}} \end{array}} \right).$$ (49c) Substituting eq. (49) into eq. (48) yields the general stability condition for the Nth-order lattice given by   $$0 \le \displaystyle\frac{{\Delta {t^2}}}{{\rho {a^2}}}\left( {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{{a^2}}} + {k_1}} \right)\sum\limits_{m = 1}^{{N \mathord{\left/ {\vphantom {N 2}} \right.} 2}} {{h_{2m - 1}}} \le 1,$$ (50a)  $$0 \le \displaystyle\frac{{\Delta {t^2}}}{{\rho {a^2}}}\left( {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{{a^2}}} + {k_3}} \right)\sum\limits_{m = 1}^{{N \mathord{\left/ {\vphantom {N 2}} \right.} 2}} {{h_{2m - 1}}} \le 1,$$ (50b)  $$0 \le \displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}\left( {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{{a^2}}} + {k_1} + {k_2} + {k_4} \pm \sqrt {k_1^2 + {{\left( {\displaystyle\frac{{{b_{23}} - {b_{12}}}}{{{a^2}}} + {k_4} - {k_2}} \right)}^2}} } \right)\sum\limits_{m = 1}^{{N \mathord{\left/ {\vphantom {N 2}} \right.} 2}} {{h_{2m - 1}}} \le 1,$$ (50c)and   $$0 \le \displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}\left( {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{{a^2}}} + {k_3} + {k_2} + {k_4} \pm \sqrt {k_3^2 + {{\left( {\displaystyle\frac{{{b_{23}} - {b_{12}}}}{{{a^2}}} + {k_4} - {k_2}} \right)}^2}} } \right)\sum\limits_{m = 1}^{{N \mathord{\left/ {\vphantom {N 2}} \right.} 2}} {{h_{2m - 1}}} \le 1,$$ (50d)where $${N \mathord{/ {\vphantom {N 2}}} 2}$$ represents the division of the integer N by 2. The stability analysis shows that the high-order dynamic lattice not only reduces the numerical dispersion but also changes the stability condition in seismic wave simulation. NUMERICAL EXAMPLES Numerical tests are implemented to validate the high-order DLM. We apply the first-, second-, third- and fourth-order DLMs in elastic wave simulation on the TI media shown in Fig. 6, where we load a single force source in the vertical direction at x = 2.25 km and z = 2.25 km to generate elastic waves, and use the absorbing boundary condition (see Appendix  B) at all boundaries. The locations of the receivers are listed in Table 1. To compare the spatial dispersions of different DLMs, we set the dominant frequency of the Ricker source wavelet as 50 Hz and the lattice spacing as 5 m.The lattice network for this simulation includes 900 × 900 particles. The time interval is 0.5 ms, and the simulation goes through 1300 time steps. Table 2 shows the CPU times for the different DLMs when simulating the elastic wave propagation. The snapshots of the x-component of the displacement fields from the DLMs are plotted in Fig. 7. These snapshots reveal that the spatial dispersion decreases as the order of the lattice unit increases. We also compare the seismograms from the high-order DLMs with those from the DRP (dispersion relation preserving)/opt MacCormack FDM (Zhu et al. 2009; Zhang et al. 2012) in Fig. 8. This particular FDM solves the first-order partial differential velocity–stress equations with low dispersion or dissipation; it is non-staggered and fourth-order accurate in space. The seismogram comparison is also conducted on the model given in Fig. 6. The lattice spacing for the DLMs is 1 m, whereas the grid spacing for the FDM is 0.5 m. The dominant frequency of the Ricker wavelet is 20 Hz, and the time interval, Δt, is 0.1 ms. Since the maximum frequency of the Ricker wavelet is estimated to be 2.5 times the dominant frequency, fc, the relation between the grid spacing and point-per-wavelength (PPW) is given by   $$\Delta h \approx \displaystyle\frac{{{V_S}}}{{2.5{f_c}{\rm{PPW}}}},$$ (51)where Δh is the grid spacing for the FDM or the lattice spacing for the DLM, fc is the dominant frequency and VS is the shear wave velocity. For the TTI media, we substitute VS with quasi-S-wave velocity, VS0, which, according to Thomsen (1986), is given by   $${V_{S{\rm{0}}}} = \sqrt {\displaystyle\frac{{{c_{44}}}}{\rho }} ,$$ (52)where c44 is the component of the elasticity tensor, and ρ represents the density. Therefore, the PPW for the DLM simulations is approximately 33, whereas the PPW for the FDM simulation is approximately 66. Seismograms are recorded at receivers 1–4. The comparison indicates that when the lattice spacing is much smaller than the wavelength, the seismograms from different DLMs and those from the FDM coincide with each other. On the other hand, it can also be witnessed that the differences between the high-order DLMs and the FDM are smaller than those between the first-order DLM and the FDM. Figure 6. View largeDownload slide Transversely isotropic (TI) model used for the seismic wave simulations. Figure 6. View largeDownload slide Transversely isotropic (TI) model used for the seismic wave simulations. Figure 7. View largeDownload slide Snapshots of the x-components of the displacement fields generated by different DLMs. Parts of the displacement fields in the red frames are enlarged for detailed comparison. Figure 7. View largeDownload slide Snapshots of the x-components of the displacement fields generated by different DLMs. Parts of the displacement fields in the red frames are enlarged for detailed comparison. Figure 8. View largeDownload slide Seismograms recorded at receivers 1–4. The differences (black lines) between the seismograms from the DLMs (red solid lines) and those from the DRP/opt MacCormack FDM (blue dashed lines) are amplified to be significant. Figure 8. View largeDownload slide Seismograms recorded at receivers 1–4. The differences (black lines) between the seismograms from the DLMs (red solid lines) and those from the DRP/opt MacCormack FDM (blue dashed lines) are amplified to be significant. Table 1. Locations of the receivers. Receiver  1  2  3  4  5  6  7  x (km)  2.70  2.60  2.45  2.20  2.25  2.49  2.49  z (km)  2.25  2.60  2.65  2.65  2.49  2.49  2.25  Receiver  1  2  3  4  5  6  7  x (km)  2.70  2.60  2.45  2.20  2.25  2.49  2.49  z (km)  2.25  2.60  2.65  2.65  2.49  2.49  2.25  View Large Table 2. CPU times for different DLMs Order  1  2  3  4  t (s)  248.0  593.9  857.6  1214.2  Order  1  2  3  4  t (s)  248.0  593.9  857.6  1214.2  View Large To investigate the spatial dispersion for different methods, we carry out elastic wave simulations on the same TTI model shown in Fig. 6 with different spacings using different DLMs The dominant frequency for the source wavelet is 30 Hz, and the time interval is 0.1 ms. Fig. 9 shows the vertical displacements recorded at receiver 6 from different numerical methods with lattice spacing varying from 2 to 8 m. The seismogram generated by the DRP/opt MacCormack FDM with a grid spacing of 0.5 m is used as a reference. The seismograms in Fig. 9 show the effect of changing the spacing on different methods. The first-order DLM appears to be most affected by spatial dispersion when the spacing varies. When the spacing gets bigger, the seismograms from high-order DLMs, especially those from the third- and fourth-order DLMs, still fit the referencing seismogram well. We also compare the convergence for different DLMs as the spacing decreases in Fig. 10, where the corresponding PPW is used as the second horizontal axis and the misfit error is given by   $${\rm{error}} = \displaystyle\frac{{\sum\limits_i^N {{{\left[ {S\left( {{t_i}} \right) - {S^0}\left( {{t_i}} \right)} \right]}^2}} }}{{\sum\limits_i^N {{{\left[ {{S^0}\left( {{t_i}} \right)} \right]}^2}} }},$$ (53)in which S(ti) is the seismogram generated by a certain numerical scheme with a given spacing from a certain receiver, and S0(ti) is the referencing seismogram generated by the DRP/opt FDM with a grid spacing of 0.5 m from the same receiver. Through this comparison, we can clearly see that the seismograms from the high-order DLMs quickly converge to the referencing seismogram and have less misfit error at a large lattice spacing than the first-order DLM. Figure 9. View largeDownload slide Seismograms from DLMs when the spacing is 2, 5 and 8 m (from top to bottom). The reference seismogram is obtained by the DRP/opt MacCormack FDM with a grid spacing of 0.5 m. Figure 9. View largeDownload slide Seismograms from DLMs when the spacing is 2, 5 and 8 m (from top to bottom). The reference seismogram is obtained by the DRP/opt MacCormack FDM with a grid spacing of 0.5 m. Figure 10. View largeDownload slide Convergence curves for the DLMs in terms of the misfit error with respect to the seismograms recorded at different receivers with decreasing spacing. Figure 10. View largeDownload slide Convergence curves for the DLMs in terms of the misfit error with respect to the seismograms recorded at different receivers with decreasing spacing. We test different DLMs on a more realistic model shown in Fig. 11. A single force source in the vertical direction is located at x = 2 km and z = 0.6 km. The dominant frequency of the Ricker wavelet is 20 Hz, and the time interval is 0.5625 ms. The spacing is set to be 2.5 m. Receivers are arranged in three seismic profiles. The locations of the receivers are listed in Tables 3 –5. Snapshots of the horizontal displacements obtained by different methods are shown in Fig. 12. Detailed seismogram comparisons are shown in Figs 13 and 14. The results obtained from different methods coincide, demonstrating the ability of the high-order DLM to handle realistic media. Figure 11. View largeDownload slide The TTI overthrust model. The elastic properties of the TTI overthrust are represented by the Thomsen moduli. Each of the four thrust blocks has its own tilted angle. The location of the single force source is marked by the red star. Receivers are arranged in three profiles to record the seismograms. Figure 11. View largeDownload slide The TTI overthrust model. The elastic properties of the TTI overthrust are represented by the Thomsen moduli. Each of the four thrust blocks has its own tilted angle. The location of the single force source is marked by the red star. Receivers are arranged in three profiles to record the seismograms. Figure 12. View largeDownload slide Snapshots of the horizontal displacement wavefields simulated by the DRP/opt MacCormack FDM and different DLMs. Figure 12. View largeDownload slide Snapshots of the horizontal displacement wavefields simulated by the DRP/opt MacCormack FDM and different DLMs. Figure 13. View largeDownload slide Seismogram comparison of horizontal displacements between different DLMs (red solid lines) and the DRP/opt MacCormack FDM (blue dashed lines). Records from the receivers in profiles 1, 2 and 3 are plotted in the first, second and third columns, respectively. Amplified differences between these two kinds of methods are represented by the black dotted lines. Figure 13. View largeDownload slide Seismogram comparison of horizontal displacements between different DLMs (red solid lines) and the DRP/opt MacCormack FDM (blue dashed lines). Records from the receivers in profiles 1, 2 and 3 are plotted in the first, second and third columns, respectively. Amplified differences between these two kinds of methods are represented by the black dotted lines. Figure 14. View largeDownload slide Seismogram comparison of vertical displacements between different DLMs (red solid lines) and the DRP/opt MacCormack FDM (blue dashed lines). Records from the receivers in profiles 1, 2 and 3 are plotted in the first, second and third columns, respectively. Amplified differences between the DLMs and the FDM are represented by the black dotted lines. Figure 14. View largeDownload slide Seismogram comparison of vertical displacements between different DLMs (red solid lines) and the DRP/opt MacCormack FDM (blue dashed lines). Records from the receivers in profiles 1, 2 and 3 are plotted in the first, second and third columns, respectively. Amplified differences between the DLMs and the FDM are represented by the black dotted lines. Table 3. Locations of the receivers in seismic profile 1. Receiver  1  2  3  4  5  6  x (km)  1.80  1.80  1.80  1.80  1.80  1.80  z (km)  0.20  0.35  0.50  0.65  0.80  0.95  Receiver  1  2  3  4  5  6  x (km)  1.80  1.80  1.80  1.80  1.80  1.80  z (km)  0.20  0.35  0.50  0.65  0.80  0.95  View Large Table 4. Locations of the receivers in seismic profile 2. Receiver  1  2  3  4  5  6  x (km)  2.10  2.25  2.40  2.55  2.70  2.85  z (km)  0.80  0.80  0.80  0.80  0.80  0.80  Receiver  1  2  3  4  5  6  x (km)  2.10  2.25  2.40  2.55  2.70  2.85  z (km)  0.80  0.80  0.80  0.80  0.80  0.80  View Large Table 5. Locations of the receivers in seismic profile 3. Receiver  1  2  3  4  5  6  x (km)  3.40  3.40  3.40  3.40  3.40  3.40  z (km)  0.20  0.35  0.50  0.65  0.80  0.95  Receiver  1  2  3  4  5  6  x (km)  3.40  3.40  3.40  3.40  3.40  3.40  z (km)  0.20  0.35  0.50  0.65  0.80  0.95  View Large On the other hand, it is important to discuss the computational cost of this high-order numerical scheme. However, the exact computational cost depends on many factors, such as the efficiency of the code, the architecture of the hardware and the compiler. Therefore, discussions of the computational cost of the high-order DLM depend on the amount of required variables and multiplications because they can be regarded as the theoretical measurements of the computational cost that eliminate other influences, such as coding and hardware. For a n × n lattice network that contains n2 particles, the Nth-order DLM needs 32N × n2 variables to store the coefficients of the linear combination that are used in eq. (33) to calculate the interaction force on each particle, and another 6n2 variables are needed to store the horizontal and vertical components of the particle displacements in three consecutive time steps. In terms of the CPU time consumption, the high-order DLM spends most of the CPU time calculating the interaction force acting upon the particles through eq. (33). For the Nth-order DLM, it takes 32N × n2 multiplications to obtain the interaction forces on every particle at each time step. Therefore, the CPU time spent calculating the interaction force for the Nth-order DLM is N times that of the first-order DLM. This trend is witnessed in the CPU times for different DLMs listed in Table 2. Because the high-order DLM significantly reduces the dispersion, it allows us to use a larger lattice spacing. This will considerably reduce the particle density as well as the computational cost. CONCLUSIONS In an effort to improve the accuracy and suppress the spatial dispersion of the original DLM, we analyse the properties obtained using the different layouts of lattice units. Due to the special arrangements of the crisscross- and x-shaped lattice units, these two types of basic lattice units are limited to certain elastic media. As a result, the basic square lattice is used to develop a high-order DLM. Our research indicates that the spatial dispersion decreases as we add more particles to the lattice unit. Based on the numerical dispersion analysis, we build the theory for a high-order DLM. According to our theory, the accuracy of the high-order DLM depends on the order of the lattice unit. The higher the order of the lattice unit is, the more layers of particles are included in the lattice unit. The high-order DLM using the Nth-order lattice unit is of 2Nth-order accuracy in terms of lattice spacing. On the other hand, to better elaborate the high-order DLM, we also derive the stability condition for the method. Numerical tests show that our high-order DLM significantly reduces the spatial dispersion and improves the accuracy as well as the convergence rate for seismic wave simulation in TI media. The detailed comparison of the accuracy between the referencing FDM and different DLMs validates the theory of the high-order DLM. Based on our theory, one can precisely choose the order of the DLM based on the accuracy requirement. The high-order DLM allows us to use a larger lattice spacing, which can greatly reduce the computational cost. This improvement will greatly benefit the application of the DLM in seismic wave simulations. ACKNOWLEDGEMENTS The authors thank W. Zhang for kindly providing the FD program and gratefully acknowledge X. Qiang, Q.-L. Li, Y. Jiang, L. Yang, C. Hu and Q.-H. Li for their fruitful discussions. This study received support from National Key Research and Development Program of China (no. 2017YFB0202903) and the National Natural Science Foundation of China (41774121 and 41374006). REFERENCES Aki K., Richards P.G., 2002. Quantitative Seismology , University Science Books. Cundall P.A., Strack O.D., 1979. A discrete numerical model for granular assemblies, Géotechnique , 29, 47– 65. https://doi.org/10.1680/geot.1979.29.1.47 Google Scholar CrossRef Search ADS   Del Valle-García R., Sánchez-Sesma F.J., 2003. Rayleigh waves modeling using an elastic lattice model, Geophys. Res. Lett. , 30, 1866, doi:10.1029/2003GL017600. https://doi.org/10.1029/2003GL017600 Google Scholar CrossRef Search ADS   Hixon R., 2000. Prefactored small-stencil compact schemes, J. Comput. Phys. , 165, 522– 541. https://doi.org/10.1006/jcph.2000.6631 Google Scholar CrossRef Search ADS   Hoover W.G., Ashurst W.T., Olness R.J., 1974. Two-dimensional computer studies of crystal stability and fluid viscosity, J. Chem. Phys. , 60, 4043– 4047. https://doi.org/10.1063/1.1680855 Google Scholar CrossRef Search ADS   Hu X., Jia X., 2016. A dynamic lattice method for elastic seismic modeling in anisotropic media, Geophysics , 81( 3), T131– T143. https://doi.org/10.1190/geo2015-0511.1 Google Scholar CrossRef Search ADS   Kelly K.R., Ward R.W., Treitel S., Alford R.M., 1976. Synthetic seismograms—a finite difference approach, Geophysics , 41, 2– 27. https://doi.org/10.1190/1.1440605 Google Scholar CrossRef Search ADS   Levander A.R., 1988. Fourth-order finite-difference P-SV seismograms, Geophysics , 53( 11), 1425– 1436. https://doi.org/10.1190/1.1442422 Google Scholar CrossRef Search ADS   Lubbe R., Worthington M.H., 2006. A field investigation of fracture compliance, Geophys. Prospect. , 54( 3), 319– 331. https://doi.org/10.1111/j.1365-2478.2006.00530.x Google Scholar CrossRef Search ADS   Madariaga R., 1976. Dynamics of an expanding circular fault, Bull. seism. Soc. Am. , 66, 163– 182. Moczo P., Kristek J., Vavrycuk V., Archuleta R., Halda L., 2002. 3D heterogeneous staggered-grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities, Bull. seism. Soc. Am. , 92, 3042– 3066. https://doi.org/10.1785/0120010167 Google Scholar CrossRef Search ADS   Monette L., Anderson M.P., 1994. Elastic and fracture properties of the two-dimensional triangular and square lattices, Model. Simul. Mater. Sci. Eng. , 2, 53– 66. https://doi.org/10.1088/0965-0393/2/1/004 Google Scholar CrossRef Search ADS   Mora P., Place D., 1998. Numerical simulation of earthquake faults with gouge: toward a comprehensive explanation for the heat flow paradox, J. geophys. Res. , 103, 21 067– 21 089. https://doi.org/10.1029/98JB01490 Google Scholar CrossRef Search ADS   Möllhoff M., Bean C.J., 2008. Validation of elastic wave measurements of rock fracture compliance using numerical discrete particle simulations, Geophys. Prospect. , 7( 5), 883– 895. O'Brien G.S., Bean C.J., 2004. A 3D discrete numerical elastic lattice method for seismic wave propagation in heterogeneous media with topography, Geophys. Res. Lett. , 31, L14608, doi:10.1029/2004GL020069. https://doi.org/10.1029/2004GL020069 Google Scholar CrossRef Search ADS   O'Brien G.S., Bean C.J., 2009. Volcano topography, structure and intrinsic attenuation: their relative influences on a simulated 3D visco-elastic wavefield, J. Volc. Geotherm. Res. , 183, 122– 136. https://doi.org/10.1016/j.jvolgeores.2009.03.004 Google Scholar CrossRef Search ADS   O'Brien G.S, Bean C.J., 2011. An irregular lattice method for elastic wave propagation, Geophys. J. Int. , 187( 3), 1699– 1707. https://doi.org/10.1111/j.1365-246X.2011.05229.x Google Scholar CrossRef Search ADS   O'Brien G.S., Bean C.J., Tapamo H., 2009. Dispersion analysis and computational efficiency of elastic lattice methods for seismic wave propagation, Comput. Geosci. , 35, 1768– 1775. https://doi.org/10.1016/j.cageo.2008.12.004 Google Scholar CrossRef Search ADS   Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1986. Numerical Recipes in FORTRAN: The Art of Scientific Computing , Cambridge Univ. Press. Saenger E.H., Bohlen T., 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid, Geophysics , 69, 583– 591. https://doi.org/10.1190/1.1707078 Google Scholar CrossRef Search ADS   Saenger E.H., Gold N., Shapiro S.A., 2000. Modeling the propagation of elastic waves using a modified finite-difference grid, Wave Motion , 31, 77– 92. https://doi.org/10.1016/S0165-2125(99)00023-2 Google Scholar CrossRef Search ADS   Saltzer S.D., Pollard D.D., 1992. Distinct element modeling of structures formed in sedimentary overburden by extensional reactivation of basement normal faults, Tectonics , 11, 165– 174. https://doi.org/10.1029/91TC02462 Google Scholar CrossRef Search ADS   Thomsen L., 1986. Weak elastic anisotropy, Geophysics , 51, 1954– 1966. Google Scholar CrossRef Search ADS   Toomey A., Bean C.J., 2000. Numerical simulation of seismic waves using a discrete particle scheme, Geophys. J. Int. , 141, 595– 604. https://doi.org/10.1046/j.1365-246x.2000.00094.x Google Scholar CrossRef Search ADS   Toomey A., Bean C.J., Scotti O., 2002. Fracture properties from seismic data—a numerical investigation, Geophys. Res. Lett. , 29 4, 1050, doi:10.1029/2001GL013867. https://doi.org/10.1029/2001GL013867 Google Scholar CrossRef Search ADS   Virieux J., 1984. SH-wave propagation in heterogeneous media—velocity-stress finite-difference method, Geophysics , 49, 1933– 1942. https://doi.org/10.1190/1.1441605 Google Scholar CrossRef Search ADS   Virieux J., 1986. P-SV wave propagation in heterogeneous media—velocity-stress finite-difference method, Geophysics , 51, 889– 901. https://doi.org/10.1190/1.1442147 Google Scholar CrossRef Search ADS   Zhang W., Shen Y., Zhao L., 2012. Three-dimensional anisotropic seismic wave modelling in spherical coordinates by a collocated-grid finite-difference method, Geophys. J. Int. , 188, 1359– 1381. https://doi.org/10.1111/j.1365-246X.2011.05331.x Google Scholar CrossRef Search ADS   Zhu H.J., Zhang W., Chen X.F., 2009. Two-dimensional seismic wave simulation in anisotropic media by non-staggered finite difference method (in Chinese), Chin. J. Geophys. , 52, 1536– 1546. Google Scholar CrossRef Search ADS   APPENDIX A: DERIVATION OF THE RELATION BETWEEN STIFFNESS COEFFICIENTS AND THE ELASTICITY TENSOR The elastic energy in the distorted lattice unit is defined by eq. (1). Assuming θij ≈ θ0, that is, the distortion is relatively much smaller than the lattice spacing, yields   $$\cos {\theta _{ij}} - \cos {\theta _0} \approx - \sin {\theta _0}\left( {{\theta _i} - {\theta _j}} \right),$$ (A1)where θi and θj are the angles by which the bonds 0–i and 0–j rotate from their original positions, respectively. θi can be expressed as   $${\theta _i} = \displaystyle\frac{{{{{\bf u}}_{0i}} \cdot {{\bf s}}_{0i}^0}}{{{{\left| {{{\bf r}}_{0i}^0} \right|}^2}}},$$ (A2)in which u0i is the difference in displacement between particle i and 0, $${{\bf s}}_{0i}^0$$ is the vector perpendicular to $${{\bf r}}_{0i}^0$$ in a right-handed system and $$| {{{\bf s}}_{0i}^0} | = | {{{\bf r}}_{0i}^0} |$$. On the other hand, when u0i is much smaller than $${{\bf r}}_{0i}^0$$, we also have   $$\left| {{{{\bf r}}_{0i}}} \right| - \left| {{{\bf r}}_{0i}^0} \right| = \displaystyle\frac{{{{{\bf u}}_{0i}} \cdot {{\bf r}}_{0i}^0}}{{\left| {{{\bf r}}_{0i}^0} \right|}}.$$ (A3) Thus, for the elastic energy given by eq. (1), we have   $$E \approx \displaystyle\frac{1}{2}{\sum\limits_{i = 1}^{n/2} {{k_i}\left( {\displaystyle\frac{{{{{\bf u}}_{0i}} \cdot {{\bf r}}_{0i}^0}}{{\left| {{{\bf r}}_{0i}^0} \right|}}} \right)} ^2} + \displaystyle\frac{1}{2}{\sum\limits_{i = 1,j = {\bmod}\!({i,n}){\rm{ + }}1}^n {{b_{ij}}{{\sin }^2}{\theta _0}\left( {\displaystyle\frac{{{{{\bf u}}_{0i}} \cdot {{\bf s}}_{0i}^0}}{{{{\left| {{{\bf r}}_{0i}^0} \right|}^2}}} - \displaystyle\frac{{{{{\bf u}}_{0j}} \cdot {{\bf s}}_{0j}^0}}{{{{\left| {{{\bf r}}_{0j}^0} \right|}^2}}}} \right)} ^2}.$$ (A4) To find the connection between the elastic energy in the lattice unit and the elastic energy in the continuum theory given by eq. (2), we substitute u0i with the strain tensor. According to the continuum theory, u0i can be expressed as   $${{{\bf u}}_{0i}} = \left( {\displaystyle\frac{{\partial {u_x}}}{{\partial x}}{{\left( {r_{0i}^0} \right)}_x} + \displaystyle\frac{{\partial {u_x}}}{{\partial z}}{{\left( {r_{0i}^0} \right)}_z}} \right){{{\bf e}}_x} + \left( {\displaystyle\frac{{\partial {u_z}}}{{\partial x}}{{\left( {r_{0i}^0} \right)}_x} + \displaystyle\frac{{\partial {u_z}}}{{\partial z}}{{\left( {r_{0i}^0} \right)}_z}} \right){{{\bf e}}_z},$$ (A5)where $${( {r_{0i}^0} )_x}$$ and $${( {r_{0i}^0} )_z}$$ are the x- and z-components of $${{\bf r}}_{0i}^0$$, respectively; meanwhile, the strain tensor is defined by   $${\varepsilon _{xx}} = \displaystyle\frac{{\partial {u_x}}}{{\partial x}},$$ (A6a)  $${\varepsilon _{zz}} = \displaystyle\frac{{\partial {u_z}}}{{\partial z}},$$ (A6b)and   $${\varepsilon _{xz}} = {\varepsilon _{zx}} = \displaystyle\frac{1}{2}\left( {\displaystyle\frac{{\partial {u_x}}}{{\partial z}} + \displaystyle\frac{{\partial {u_z}}}{{\partial x}}} \right),$$ (A6c)where $$\displaystyle\frac{{\partial {u_x}}}{{\partial x}}$$, $$\displaystyle\frac{{\partial {u_z}}}{{\partial z}}$$, $$\displaystyle\frac{{\partial {u_x}}}{{\partial z}}$$ and $$\displaystyle\frac{{\partial {u_z}}}{{\partial x}}$$ are the first-order derivatives of the displacement field. For the crisscrossed lattice unit (Fig. 1a), the x- and z-components of $${{\bf r}}_{0i}^0$$ are given by   $$\left( {\begin{array}{c@{\quad}c} {{{\left( {r_{01}^0} \right)}_x}}&{{{\left( {r_{01}^0} \right)}_z}}\\ {{{\left( {r_{02}^0} \right)}_x}}&{{{\left( {r_{02}^0} \right)}_z}}\\ {{{\left( {r_{03}^0} \right)}_x}}&{{{\left( {r_{03}^0} \right)}_z}}\\ {{{\left( {r_{04}^0} \right)}_x}}&{{{\left( {r_{04}^0} \right)}_z}} \end{array}} \right) = \left( {\begin{array}{c@{\quad}c} a&0\\ 0&a\\ { - a}&0\\ 0&{ - a} \end{array}} \right),$$ (A7)where a is the lattice spacing. Substituting eqs (A7), (A6) and (A5) into eq. (A4) yields   $$E = \displaystyle\frac{1}{2}{k_1}{a^2}{\varepsilon _{xx}}^2 + \displaystyle\frac{1}{2}{k_2}{a^2}{\varepsilon _{zz}}^2 + 2\left( {{b_{12}} + {b_{23}} + {b_{34}} + {b_{41}}} \right){\varepsilon _{xz}}^2.$$ (A8) Eq. (3) is obtained by comparing eq. (A8) with eq. (2) through   $$E = {a^2}\varphi ,$$ (A9)where φ is the elastic energy density given by eq. (2). In the case of the x-shaped lattice unit (Fig. 1b), the x- and z-components of $${{\bf r}}_{0i}^0$$ are given by   $$\left( {\begin{array}{c@{\quad}c} {{{\left( {r_{01}^0} \right)}_x}}&{{{\left( {r_{01}^0} \right)}_z}}\\ {{{\left( {r_{02}^0} \right)}_x}}&{{{\left( {r_{02}^0} \right)}_z}}\\ {{{\left( {r_{03}^0} \right)}_x}}&{{{\left( {r_{03}^0} \right)}_z}}\\ {{{\left( {r_{04}^0} \right)}_x}}&{{{\left( {r_{04}^0} \right)}_z}} \end{array}} \right) = \left( {\begin{array}{c@{\quad}c} {\displaystyle\frac{{\sqrt 2 }}{2}a}&{\displaystyle\frac{{\sqrt 2 }}{2}a}\\ { - \displaystyle\frac{{\sqrt 2 }}{2}a}&{\displaystyle\frac{{\sqrt 2 }}{2}a}\\ { - \displaystyle\frac{{\sqrt 2 }}{2}a}&{ - \displaystyle\frac{{\sqrt 2 }}{2}a}\\ {\displaystyle\frac{{\sqrt 2 }}{2}a}&{ - \displaystyle\frac{{\sqrt 2 }}{2}a} \end{array}} \right).$$ (A10) Similar to the crisscross-shaped lattice unit, the elastic energy for the x-shaped lattice unit is expressed as   $$\begin{array}{@{}*{3}{l}@{}} E& = &{\left( {{b_{12}} + {b_{23}} + {b_{34}} + {b_{41}} + \displaystyle\frac{{{a^2}{k_1}}}{4} + \displaystyle\frac{{{a^2}{k_2}}}{4}} \right)\left( {\displaystyle\frac{1}{2}{\varepsilon _{xx}}^2 + \displaystyle\frac{1}{2}{\varepsilon _{zz}}^2} \right)} + {\left( {\displaystyle\frac{{{a^2}{k_1}}}{4} + \displaystyle\frac{{{a^2}{k_2}}}{4} - {b_{12}} - {b_{23}} - {b_{34}} - {b_{41}}} \right){\varepsilon _{xx}}{\varepsilon _{zz}}}\\ && +\,{\displaystyle\frac{{{a^2}}}{2}\left( {{k_1} + {k_2}} \right){\varepsilon _{xz}}^2 + \displaystyle\frac{{{a^2}}}{2}\left( {{k_1} - {k_2}} \right)\left( {{\varepsilon _{xx}} + {\varepsilon _{zz}}} \right){\varepsilon _{xz}}} \end{array}.$$ (A11) Substituting eqs (2) and (A11) into eq. (A9) yields eq. (4). APPENDIX B: ABSORBING BOUNDARY CONDITION We use an absorbing pad to realize the absorbing boundary condition. Since the DLM is particle-based, we add a friction force in the opposite direction of the particle velocity on each particle in the absorbing pad. Therefore, the total force acting upon the particle in the absorbing pad is given by   $${{{\bf f}}_{{\rm{total}}}} = {{{\bf f}}^0} + {{{\bf f}}_{{\rm{friction}}}},$$ (B1)where ftotal is the total force, f0 is the interaction force given by eq. (33) and ffriction is the friction force, which takes the form of   $${{{\bf f}}_{{\rm{friction}}}} = - \mu {{\bf \dot{u}}},$$ (B2)where μ is the friction factor and $${{\bf \dot{u}}}$$ is the particle velocity. This friction force dampens the kinetic energy of the particle. In practice, we increase the friction factor as the absorbing pad expands to better attenuate the motions of the particles. Fig. B1 shows the effect of this absorbing pad. Figure B1. View largeDownload slide Snapshots of the wavefield at different time steps. Absorbing layers are added on the left and on the top of the area. The interfaces between the absorbing layers and the full elastic area are marked by the black lines. The boundary conditions on the right and bottom represent the fixed boundary conditions. Figure B1. View largeDownload slide Snapshots of the wavefield at different time steps. Absorbing layers are added on the left and on the top of the area. The interfaces between the absorbing layers and the full elastic area are marked by the black lines. The boundary conditions on the right and bottom represent the fixed boundary conditions. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# High-order dynamic lattice method for seismic simulation in anisotropic media

, Volume 212 (3) – Mar 1, 2018
22 pages

/lp/ou_press/high-order-dynamic-lattice-method-for-seismic-simulation-in-gLF7SJzrVF
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx508
Publisher site
See Article on Publisher Site

### Abstract

Summary The discrete particle-based dynamic lattice method (DLM) offers an approach to simulate elastic wave propagation in anisotropic media by calculating the anisotropic micromechanical interactions between these particles based on the directions of the bonds that connect them in the lattice. To build such a lattice, the media are discretized into particles. This discretization inevitably leads to numerical dispersion. The basic lattice unit used in the original DLM only includes interactions between the central particle and its nearest neighbours; therefore, it represents the first-order form of a particle lattice. The first-order lattice suffers from numerical dispersion compared with other numerical methods, such as high-order finite-difference methods, in terms of seismic wave simulation. Due to its unique way of discretizing the media, the particle-based DLM no longer solves elastic wave equations; this means that one cannot build a high-order DLM by simply creating a high-order discrete operator to better approximate a partial derivative operator. To build a high-order DLM, we carry out a thorough dispersion analysis of the method and discover that by adding more neighbouring particles into the lattice unit, the DLM will yield different spatial accuracy. According to the dispersion analysis, the high-order DLM presented here can adapt the requirement of spatial accuracy for seismic wave simulations. For any given spatial accuracy, we can design a corresponding high-order lattice unit to satisfy the accuracy requirement. Numerical tests show that the high-order DLM improves the accuracy of elastic wave simulation in anisotropic media. Numerical modelling, Computational seismology, Seismic anisotropy, Wave propagation INTRODUCTION In most cases, the seismic wave phenomenon is interpreted using the elastic wave equation, which is based on the continuum theory. Analytical solutions to the equation only apply to simple media (Aki & Richards 2002). Therefore, a variety of numerical methods, such as the finite-difference method (FDM, e.g. Kelly et al. 1976; Madariaga 1976; Virieux 1984, 1986; Saenger et al. 2000), have been developed to yield approximate solutions to the wave equation. In regard to the numerical simulation of seismic wave propagation, due to the limitations of digital computers, certain discretization in either the space domain or the wavenumber domain is required. However, wave equation-based numerical methods are governed by a spatially continuous differential equation. Special calculations of particular boundary conditions at discontinuities are required. As a result, although the model and wavefields are already discretized into discrete gridpoints and are no longer continuous, wave equation-based methods still need to smooth the model to satisfy the continuity requirement of the wave equation; additionally, it is difficult for these methods to incorporate the naturally discontinuous features of geological materials, such as fractures and faults. Meanwhile, another category of numerical methods has been quite successful in simulating dynamic fractures and fault propagation. Cundall & Strack (1979) developed one of the first discrete particle methods, which was called the ‘distinct element method’, by extending the solid-state physics models of crystalline materials (Hoover et al. 1974) to study the static deformation and fracturing of granular materials. This methodology has since been used to investigate earthquake dynamics (Mora & Place 1998) and tectonic processes (Saltzer & Pollar 1992). Instead of directly solving the partial differential equation, these discrete methods treat geological materials as assemblages of numerous interconnected discrete particles and focus on the ‘microscopic’ interactions of these particles. Therefore, these methods eliminate the continuity restriction of the wave equation, offering a possible way to incorporate geological discontinuity. For this reason, particle-based methods are of great flexibility and can be used to solve different geological problems. Toomey & Bean (2000) discussed a 2-D ‘discrete particle scheme’ for seismic wave simulation in isotropic media. Their method assumes the medium to comprise discrete particles that are closely packed in a hexagonal lattice. In this lattice, particles interact with their bonded neighbours through Hooke's law. Toomey et al. (2002) also used this discrete particle scheme to investigate fracture properties from seismic data. Further studies on this issue were carried out by Lubbe & Worthington (2006) and Möllhoff & Bean (2008) using the same scheme. This method only takes central forces into account; thus, the Poisson ratio of the media is fixed at 0.25. Del Valle-García & Sánchez-Sesma (2003) improved Toomey's method by including a bond-bending force term to discard this Poisson ratio restriction and implemented it to simulate Rayleigh waves on the free surface. O'Brien & Bean (2004) extended this elastic lattice method to a 3-D case by using a 3-D lattice model. Their method has since been employed to study the structural topography of volcanoes (O'Brien & Bean 2009). Hu & Jia (2016) discussed a particle-based dynamic lattice method (DLM) that can be used to simulate seismic wave propagation in transversely isotropic (TI) media. On the other hand, due to their use of discretization, numerical methods face spatial dispersion that causes error in seismic wave simulation. To achieve certain accuracy, these methods must increase the density of the gridpoints, which will significantly increase memory consumption and slow down the simulation. To improve the performance of wave simulation, numerical methods such as the FDM have developed corresponding high-order schemes (e.g. Levander 1988; Hixon 2000; Moczo et al. 2002; Saenger & Bohlen 2004) that can reduce this spatial dispersion. Originally, the DLM only considers the interactions between the nearest particles, which results in great spatial dispersion compared with other high-order numerical methods. Therefore, reducing the spatial dispersion of the DLM has become an important issue. Meanwhile, for the wave equation-based methods, such as the FDM, one can easily use high-order finite-difference operators to better approximate the partial derivative operator. For the particle-based DLM, however, this is not the case because the method does not solve the wave equation directly. Therefore, the challenge lies in building high-order particle-based methods. By retaining the nearest-particle restriction and including an additional local interaction force acting on individual particles, O'Brien et al. (2009) managed to decrease the numerical dispersion of their elastic lattice method for isotropic media. However, in the case of anisotropy, which is represented by the anisotropic interactions according to dynamic lattice theory, it is difficult to suppress numerical dispersion by simply adding such an interaction force without considering the micromechanical interactions between particles other than those between the closest neighbours. The particle-based irregular lattice method (O'Brien & Bean 2011), which takes more than the nearest neighbours into account, is brought up to simulate elastic wave propagation. Although the implementation of the irregular lattice removes the anisotropic nature of fracture propagation, it also results in more dispersion than using the regular lattice. To fully analyse the spatial dispersion of the DLM, we take the different interaction patterns of the particle lattice into account and establish the high-order DLM. According to our analysis, by including more particles into the lattice unit, this method will produce less spatial dispersion, as the magnitude of the dispersion is related to the number of the added particles in the lattice unit. Using the high-order DLM, we can achieve higher numerical accuracy with relatively larger lattice spacing, thereby reducing the memory consumption and improving the CPU time efficiency of elastic wave simulation in anisotropic media. BASIC LATTICE UNITS IN THE DLM The most prominent feature of the DLM is its unique method of discretizing the media into a particle lattice. Thus, the layout of the lattice unit plays an important role in this method (Monette & Anderson 1994). Here, the three different layouts of the basic lattice unit shown in Fig. 1 are investigated to obtain insight into this method and lay the foundation of the high-order DLM. In such a basic lattice unit, the central particle 0 is bonded with its neighbours. As Fig. 2 shows, particle interactions in the basic lattice unit are achieved through the stretching and bending of these bonds. The stiffness of the bonds is calibrated by the coefficients ki and bij when the bonds are stretched and bent, respectively. The bigger these coefficients are, the harder it is to stretch or bend these bonds. In the case of anisotropy, because these stiffness coefficients determine the elastic interactions between these particles, thereby defining the elastic properties of the particle lattice, we create an anisotropic lattice by setting the stiffness coefficients in different directions. According to Monette & Anderson (1994), the elastic energy in the basic square lattice unit is given by   $$E = \frac{1}{2}{\sum\limits_{i = 1}^{{n \mathord{\left/ {\vphantom {n 2}} \right.} 2}} {{k_i}\left( {\left| {{{{\bf r}}_{0i}}} \right| - \left| {{{\bf r}}_{0i}^0} \right|} \right)} ^2} + \frac{1}{2}{\sum\limits_{i = 1,j ={\bmod}\!(i,n) + 1}^n {{b_{ij}}\left( {\cos {\theta _{ij}} - \cos {\theta _0}} \right)} ^2},$$ (1)in which E is the elastic energy, and r0i and $${{\bf r}}_{0i}^0$$ are the vectors of the bond that connects particle 0 and i (bond 0–i for short) in the distorted and undistorted lattices, respectively. Similarly, θij and θ0 are the angles between bonds 0–i and 0–j in the distorted and undistorted lattices, respectively. mod(i, n) denotes the remainder function after i is divided by the amount of lattice bonds in the basic lattice unit, that is, n. Figure 1. View largeDownload slide Basic lattice units. (a) Crisscross-shaped unit, (b) X-shaped unit and (c) basic square unit. The central particle 0 is bonded with its neighbours. The coefficients ki and bij represent the stiffness of the bonds in terms of stretching and bending, respectively. a is the lattice spacing that denotes the spatial interval between adjacent particles. Figure 1. View largeDownload slide Basic lattice units. (a) Crisscross-shaped unit, (b) X-shaped unit and (c) basic square unit. The central particle 0 is bonded with its neighbours. The coefficients ki and bij represent the stiffness of the bonds in terms of stretching and bending, respectively. a is the lattice spacing that denotes the spatial interval between adjacent particles. Figure 2. View largeDownload slide Particle interactions through the stretching and bending of bonds. (a) It shows the stretching of the bond that connects particles 0 and i. The coefficient ki represents the stiffness of the bond in terms of stretching. (b) It demonstrates the bending of the bonds that connect particles i, j and 0 with particle 0 as the pivot. The bond-bending coefficient bij denotes the stiffness of these bonds in terms of bending. Figure 2. View largeDownload slide Particle interactions through the stretching and bending of bonds. (a) It shows the stretching of the bond that connects particles 0 and i. The coefficient ki represents the stiffness of the bond in terms of stretching. (b) It demonstrates the bending of the bonds that connect particles i, j and 0 with particle 0 as the pivot. The bond-bending coefficient bij denotes the stiffness of these bonds in terms of bending. Once the lattice unit is chosen, it is necessary to find the relationship between the micromechanical stiffness coefficients and the macroscopic parameters. According to Hu & Jia (2016), we can find the transformation by comparing the elastic energy in a lattice unit given by eq. (1) using the strain energy density φ of a 2-D TI medium with a vertical axis of symmetry (VTI medium) given by   $$\varphi = \frac{1}{2}{c_{11}}{\varepsilon _{xx}}^2 + \frac{1}{2}{c_{33}}{\varepsilon _{zz}}^2 + 2{c_{44}}{\varepsilon _{xz}}^2 + {c_{13}}{\varepsilon _{xx}}{\varepsilon _{zz}},$$ (2)where εxx, εzz and εxz denote the corresponding components of the strain tensor. c11, c33,c13 and c44 represent the corresponding elements of the elasticity tensor of the VTI medium. According to eqs (1) and (2), for the crisscross-shaped lattice unit in Fig. 1(a), where the horizontal interval between neighbouring particles is equal to the vertical interval, we obtain the relationship between the stiffness coefficients and the elasticity tensor for VTI media given by (see Appendix  A)   $${c_{13}} = 0,$$ (3a)  $${k_1} = {c_{11}},$$ (3b)  $${k_2} = {c_{33}},$$ (3c)and   $${b_{12}} = {b_{23}} = {b_{34}} = {b_{41}} = \frac{1}{4}{a^2}{c_{44}},$$ (3d)where a is the lattice spacing, indicating the interval between adjacent particles. For the x-shaped lattice unit in Fig. 1(b), this relationship has the form of   $${c_{11}} = {c_{33}},$$ (4a)  $${k_1} = {k_2} = 2{c_{44}},$$ (4b)  $${c_{13}} = 2{c_{44}} - {c_{11}}$$ (4c)and   $${b_{12}} = {b_{23}} = {b_{34}} = {b_{41}} = \frac{{{a^2}}}{4}\left( {{c_{11}} - {c_{44}}} \right).$$ (4d) According to eqs (3) and (4), both the crisscross- and the x-shaped lattice units are strongly restricted in terms of their elastic properties, that is, both units only represent some particular VTI media. The reason for this restriction is that the crisscross- and x-shaped lattice units only take four of the neighbouring particles into account. The interactions between these particles are thus insufficient to simulate all kinds of VTI media. The basic square unit plotted in Fig. 1(c) is fully discussed by Hu & Jia (2016); when b12 = b34 = b56 = b78 and b23 = b45 = b67 = b81, the relation between the stiffness coefficients of the square lattice and the elasticity tensor is given by   $${{\bf W}}\left( \theta \right)\left( {\begin{array}{c} {{k_1}}\\ {{k_2}}\\ {{k_3}}\\ {{k_4}}\\ {{b_{12}}}\\ {{b_{23}}} \end{array}} \right) = \left( {\begin{array}{c} {{c_{11}}}\\ {{c_{33}}}\\ {{c_{44}}}\\ {{c_{13}}}\\ 0\\ 0 \end{array}} \right),$$ (5a)where c11, c33, c13 and c44 are the elements in the elasticity tensor of the corresponding TI medium when the symmetry axis is in the vertical direction, θ is the tilted angle for the tilted TI medium (TTI medium) and W(θ) is a 6 × 6 matrix that can be written as   $${{\bf W}}\left( \theta \right) = \left( {\begin{array}{c{@}{\quad}c{@}{\quad}c{@}{\quad}c{@}{\quad}c{@}{\quad}c} {{{\cos }^4}\theta }&{2{{\cos }^4}\left( {\frac{\pi }{4} + \theta } \right)}&{{{\sin }^4}\theta }&{2{{\cos }^4}\left( {\frac{\pi }{4} - \theta } \right)}&{\frac{{{{\left( {\cos 2\theta - \sin 2\theta } \right)}^2}}}{{2{a^2}}}}&{\frac{{{{\left( {\cos 2\theta + \sin 2\theta } \right)}^2}}}{{2{a^2}}}}\\ {{{\sin }^4}\theta }&{2{{\sin }^4}\left( {\frac{\pi }{4} + \theta } \right)}&{{{\cos }^4}\theta }&{2{{\sin }^4}\left( {\frac{\pi }{4} - \theta } \right)}&{\frac{{{{\left( {\cos 2\theta - \sin 2\theta } \right)}^2}}}{{2{a^2}}}}&{\frac{{{{\left( {\cos 2\theta + \sin 2\theta } \right)}^2}}}{{2{a^2}}}}\\ {\frac{{{{\sin }^2}2\theta }}{4}}&{\frac{{{{\cos }^2}2\theta }}{2}}&{\frac{{{{\sin }^2}2\theta }}{4}}&{\frac{{{{\cos }^2}2\theta }}{2}}&{\frac{{{{\left( {\cos 2\theta + \sin 2\theta } \right)}^2}}}{{2{a^2}}}}&{\frac{{{{\left( {\cos 2\theta - \sin 2\theta } \right)}^2}}}{{2{a^2}}}}\\ {\frac{{{{\sin }^2}2\theta }}{4}}&{\frac{{{{\cos }^2}2\theta }}{2}}&{\frac{{{{\sin }^2}2\theta }}{4}}&{\frac{{{{\cos }^2}2\theta }}{2}}&{ - \frac{{{{\left( {\cos 2\theta - \sin 2\theta } \right)}^2}}}{{2{a^2}}}}&{ - \frac{{{{\left( {\cos 2\theta + \sin 2\theta } \right)}^2}}}{{2{a^2}}}}\\ {{{\cos }^2}\theta \sin 2\theta }&{2{{\cos }^2}\left( {\frac{\pi }{4} + \theta } \right)\cos 2\theta }&{ - \sin 2\theta {{\sin }^2}\theta }&{ - 2{{\cos }^2}\left( {\frac{\pi }{4} - \theta } \right)\cos 2\theta }&{\frac{{\cos 4\theta }}{{{a^2}}}}&{ - \frac{{\cos 4\theta }}{{{a^2}}}}\\ {\sin 2\theta \sin^2\theta }&{2{{\sin }^2}\left( {\frac{\pi }{4} + \theta } \right)\cos 2\theta }&{ - \sin 2\theta {{\cos }^2}\theta }&{ - 2{{\sin }^2}\left( {\frac{\pi }{4} - \theta } \right)\cos 2\theta }&{ - \frac{{\cos 4\theta }}{{{a^2}}}}&{\frac{{\cos 4\theta }}{{{a^2}}}} \end{array}} \right).$$ (5b) Since the determinant of this matrix is non-zero, we can always find the corresponding stiffness coefficients for any TI medium as long as the medium can be described by the elasticity tensor. Therefore, considering the limitations of the crisscross- and the x-shaped lattice units in representing realistic media, we build the high-order DLM based on the more applicable basic square lattice unit. HIGH-ORDER LATTICE UNITS Due to spatial discretization, the original DLM suffers from spatial dispersion in terms of seismic wave simulation. For the basic square lattice unit, micromechanical interactions only exist between the nearest pairs of particles. As a result, the dispersion of the DLM cannot be ignored as it can be in other high-order wave equation-based numerical methods. To suppress this spatial dispersion, we develop a high-order dynamic lattice that can decrease the dispersion based on the accuracy requirement. According to the theory of the DLM (Hu & Jia 2016), the interaction force acting upon an individual particle is calculated by Lagrange's equation. Hu & Jia (2016) simplified the calculation for this interaction force by linearizing Lagrange's equation. In the case of homogenous media, the interaction force on the central particle in the lattice unit can be written as a linear combination of the differences in the displacements between the central particle and its neighbours given by   $$\left( {\begin{array}{c} {f_0^x}\\ {f_0^z} \end{array}} \right) = \left( {\begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{{\bf A}}_{11}}}&{{{{\bf A}}_{12}}}&{{{{\bf A}}_{13}}}&{{{{\bf A}}_{14}}}\\ {{{{\bf A}}_{21}}}&{{{{\bf A}}_{22}}}&{{{{\bf A}}_{23}}}&{{{{\bf A}}_{24}}} \end{array}} \right) \times {\Big( {\begin{array}{*{20}{c}} {u_{01}^x}&\quad \cdots &\quad{u_{08}^x}&\quad{u_{01}^z}&\quad \cdots &\quad{u_{08}^z} \end{array}} \Big)^{\rm{T}}},$$ (6)where $$f_0^x$$ and $$f_0^z$$ are the x- and z-components of the interaction force acting upon the central particle that is, particle 0 in the lattice unit (Fig. 1c), respectively; and $$u_{0i}^x$$ and $$u_{0i}^z$$ are the x- and z-components of the difference in the displacement between particles 0 and i in the lattice unit, respectively, that is,   $$\left( {\begin{array}{c} {u_{0i}^x}\\ {u_{0i}^z} \end{array}} \right) = \left( {\begin{array}{c} {u_i^x - u_0^x}\\ {u_i^z - u_0^z} \end{array}} \right),$$ (7)where $$u_i^x$$ and $$u_i^z$$ are the x- and z-components of the displacement for particle i, respectively, whereas $$u_0^x$$ and $$u_0^z$$ are the x- and z-components of the displacement for particle 0, respectively. The coefficients for this linear combination are given by   $${{{\bf A}}_{11}} = {{{\bf A}}_{12}} = \left( {\begin{array}{c} {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{2{a^2}}} + {k_1}}\quad{ - \displaystyle\frac{{{b_{23}} - {b_{12}}}}{{4{a^2}}} + \frac{1}{2}{k_2}}\quad{\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{2{a^2}}}}\quad{ - \displaystyle\frac{{{b_{12}} - {b_{23}}}}{{4{a^2}}} + \displaystyle\frac{1}{2}{k_4}} \end{array}} \right),$$ (8a)  $${{{\bf A}}_{13}} = {{{\bf A}}_{14}} = {{{\bf A}}_{21}} = {{{\bf A}}_{22}} = \left( {\begin{array}{*{20}{c}} { - \displaystyle\frac{{{b_{23}} - {b_{12}}}}{{2{a^2}}}}&\quad{\displaystyle\frac{1}{2}{k_2}}&\quad{ - \displaystyle\frac{{{b_{12}} - {b_{23}}}}{{2{a^2}}}}&\quad { - \displaystyle\frac{1}{2}{k_4}} \end{array}} \right),$$ (8b)and   $${{{\bf A}}_{23}} = {{{\bf A}}_{24}} = \left( {\begin{array}{c} {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{2{a^2}}}}\quad { - \displaystyle\frac{{{b_{12}} - {b_{23}}}}{{4{a^2}}} + \displaystyle\frac{1}{2}{k_2}}\quad{\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{2{a^2}}} + {k_3}}\quad{ - \displaystyle\frac{{{b_{23}} - {b_{12}}}}{{4{a^2}}} + \displaystyle\frac{1}{2}{k_4}} \end{array}} \right).$$ (8c) Using eq. (6), one can obtain the interaction forces acting upon the particles as long as the displacements of these particles are provided. Using the initial displacements and the interaction forces of the particles, the displacements for the next time step are obtained through the Verlet accumulation, which is given by   $$\left( {\begin{array}{c} {u_i^x\left( {{t_0} + \Delta t} \right)}\\ {u_i^z\left( {{t_0} + \Delta t} \right)} \end{array}} \right) = \left( {\begin{array}{c}{\ddot{u}_i^x\left( {{t_0}} \right)}\\ {\ddot{u}_i^z\left( {{t_0}} \right)} \end{array}} \right)\Delta {t^2} + 2\left( {\begin{array}{c} {u_i^x\left( {{t_0}} \right)}\\ {u_i^z\left( {{t_0}} \right)} \end{array}} \right) - \left( {\begin{array}{c} {u_i^x\left( {{t_0} - \Delta t} \right)}\\ {u_i^z\left( {{t_0} - \Delta t} \right)} \end{array}} \right)$$ (9)where Δt is the time interval, $$u_i^x( t )$$ and $$u_i^z( t )$$ are the x- and z-components of the displacement for particle i at time t, respectively and $$\ddot{u}_i^x( t )$$ and $$\ddot{u}_i^z( t )$$ denote the x- and z-components of the acceleration for particle i, respectively, which can be written as   $$\left( {\begin{array}{c} {\ddot{u}_i^x\left( t \right)}\\ {\ddot{u}_i^z\left( t \right)} \end{array}} \right) = \frac{1}{{\rho {a^2}}}\left( {\begin{array}{c}{f_i^x}\\ {f_i^z} \end{array}} \right),$$ (10)where ρ is the density of the media and a is the lattice spacing. Based on this linearized interaction force and the Verlet accumulation, we carry out a plane wave analysis of the method according to Saenger & Bohlen (2004) to investigate the spatial dispersion of the basic square lattice unit for wave propagation. The plane wave that propagates in the lattice unit takes the form of   $$\left( {\begin{array}{c} {u_i^x\left( t \right)}\\ {u_i^z\left( t \right)} \end{array}} \right) = \left( {\begin{array}{c} {{U_x}\exp \left( {j{{\bf k}} \cdot {{\bf r}}_{0i}^0 - j\omega t} \right)}\\ {{U_z}\exp \left( {j{{\bf k}} \cdot {{\bf r}}_{0i}^0 - j\omega t} \right)} \end{array}} \right),$$ (11)where Ux and Uz are the x- and z-components of the amplitude of the plane wave, respectively; $${{\bf r}}_{0i}^0$$ is the vector for the bond that connects particles 0 and i in the undistorted lattice unit; k is the vector wavenumber; ω represents the angular frequency of the plane wave and $$j = \sqrt { - 1}$$. Substituting eq. (11) into eqs (10) and (6) yields   $$\left( {\begin{array}{c{@}{\quad }c} {{g_{11}}}&{{g_{12}}}\\ {{g_{21}}}&{{g_{22}}} \end{array}} \right)\left( {\begin{array}{c} {{U_x}}\\ {{U_z}} \end{array}} \right) = - {\omega ^2}\rho \left( {\begin{array}{c} {{U_x}}\\ {{U_z}} \end{array}} \right)$$ (12)in which   $${g_{11}} = {{{\bf A}}_{11}} \times {{\bf L}},$$ (13a)  $${g_{12}} = {g_{21}} = {{{\bf A}}_{12}} \times {{\bf L}},$$ (13b)and   $${g_{22}} = {{{\bf A}}_{23}} \times {{\bf L}}.$$ (13c) For the basic square lattice unit, L is a function of ak given by   $${{\bf L}}\left( {a{{\bf k}}} \right) = \left( {\begin{array}{c} {\displaystyle\frac{{\cos \left( {a{{\bf k}} \cdot {{{\bf e}}_{01}}} \right) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right),$$ (14)in which $${{{\bf e}}_{0i}} = \frac{{{{\bf r}}_{0i}^0}}{{| {{{\bf r}}_{0i}^0} |}}$$. Solving eq. (12) yields   $$v\left( {a{{\bf k}}} \right) = \sqrt { - \frac{1}{\rho }\left( {{g_{11}} + {g_{22}} \pm \sqrt {{{\left( {{g_{11}} - {g_{22}}} \right)}^2} + 4{g_{12}}{g_{21}}} } \right)} ,$$ (15)where v(ak) represents the phase velocity for the plane wave that propagates in the basic square lattice as a function of lattice spacing a and the vector wavenumber k, and ‘ ± ’ denotes two different wave modes, that is, P and S waves. When |ak| approximates zero, which implies that the lattice spacing is much less than the wavelength, the phase velocity given by eq. (15) tends to represent the theoretical phase velocity vtheo of the corresponding continua, that is,   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} v\left( {a{{\bf k}}} \right) = {v_{{\rm{theo}}}}.$$ (16) The dispersion comes from the difference between v(ak) and vtheo. According to eqs (13) and 15, v(ak) can also be regarded as a function of L; thus, eq. (16) can be expressed as   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} v\left( {a{{\bf k}}} \right) = \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} f\left( {{\bf L}} \right) = f ( {\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}}\left( {a{{\bf k}}} \right)} ) = {v_{{\rm{theo}}}},$$ (17)where f(L) represents the phase velocity as a function of L. Therefore, the dispersion is in fact caused by the difference between L(ak) and $$\mathop {\lim }\limits_{| {a{{\bf k}}} | \to 0} {{\bf L}}( {a{{\bf k}}} )$$. When |ak| approximates 0, we have   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}} = \left( {\begin{array}{c} { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _1}}\\ { - {{\cos }^2}{\theta _2}}\\ { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _3}}\\ { - {{\cos }^2}{\theta _4}} \end{array}} \right),$$ (18)where θi denotes the angle between k and e0i. The difference between L(ak) and $$\mathop {\lim }\limits_{| {a{{\bf k}}} | \to 0} {{\bf L}}( {a{{\bf k}}} )$$ can be written as   $$d\left( {a{{\bf k}}} \right) = {{\bf L}}\left( {a{{\bf k}}} \right) - \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}}\left( {a{{\bf k}}} \right) = \left( {\begin{array}{c} {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _1}}\\ { - {{\cos }^2}{\theta _2}}\\ { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _3}}\\ { - {{\cos }^2}{\theta _4}} \end{array}} \right).$$ (19) Carrying out Taylor's expansion on L(ak) at ak = 0 yields   $$d\left( {a{{\bf k}}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}} + \displaystyle\frac{1}{2}{{\cos }^2}{\theta _1}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}} + {{\cos }^2}{\theta _2}}\\ {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}} + \displaystyle\frac{1}{2}{{\cos }^2}{\theta _3}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}} + {{\cos }^2}{\theta _4}} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{1}{{24}}{{\cos }^4}{\theta _1}{a^2}{{{\bf k}}^2}}\\ {\displaystyle\frac{1}{6}{{\cos }^4}{\theta _2}{a^2}{{{\bf k}}^2}}\\ {\displaystyle\frac{1}{{24}}{{\cos }^4}{\theta _3}{a^2}{{{\bf k}}^2}}\\ {\displaystyle\frac{1}{6}{{\cos }^4}{\theta _4}{a^2}{{{\bf k}}^2}} \end{array}} \right) + {\rm{Order}}( {{{\left| {a{{\bf k}}} \right|}^4}} ).$$ (20) This indicates that the deviations in the phase velocities in the basic square lattice unit from the theoretical phase velocities are at Order(|ak|2). Fig. 3 shows the relative deviations of the phase velocities of both P and S waves in the basic square lattice unit from their theoretical phase velocities, which are given by   $${d_P} = \displaystyle\frac{{{v_P}\left( {a{{\bf k}}} \right) - v_P^{{\rm{theo}}}}}{{v_P^{{\rm{theo}}}}},$$ (21a)and   $${d_S} = \displaystyle\frac{{{v_S}( {a{{\bf k}}} ) - v_S^{{\rm{theo}}}}}{{v_S^{{\rm{theo}}}}},$$ (21b)where dP and dS are the relative deviations of the phase velocities of P and S waves, respectively; vP(ak) and vS(ak) are the P- and S-wave phase velocities in the basic lattice unit, respectively and $$v_P^{{\rm{theo}}}$$ and $$v_S^{{\rm{theo}}}$$ are the theoretical phase velocities of P and S waves, respectively, in the corresponding continua. To improve the accuracy of the DLM in seismic wave simulation, we need to reduce the difference between v(ak) and vtheo. According to eq. (19), in such a basic square lattice unit, we can only reduce the dispersion by using smaller lattice spacing; however, this will significantly increase the memory consumption as well as the CPU time. Through research on modifying the lattice unit, we find that the spatial accuracy benefits from adding more neighbour particles into the lattice unit. In an effort to simultaneously reduce the spatial dispersion in terms of all lattice spacing and speed up the convergence of the particle lattice model to the corresponding continua as the lattice spacing approximates zero, we extend the basic square lattice unit in all existing directions to form a second-order lattice unit, which is shown in Fig. 4(a). Due to this modification, the linear combination given by eq. (6) is reconstructed to adapt to the new lattice unit. The coefficients for the linear combination in eq. (6) are reassigned. The modified version of eq. (6) with respect to the extended second-order unit can thus be written as   \begin{eqnarray} \left( { \begin{array}{c} {f_0^x}\\ {f_0^z} \end{array}} \right)& = &{\left( { \begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{{\bf A}}_{11}}}&{{{{\bf A}}_{12}}}&{{{{\bf A}}_{13}}}&{{{{\bf A}}_{14}}}\\ {{{{\bf A}}_{21}}}&{{{{\bf A}}_{22}}}&{{{{\bf A}}_{23}}}&{{{{\bf A}}_{24}}} \end{array}} \right) \times {{{\bf H}}_1} \times \left( { \begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{\left( {u_{01}^x} \right)}_1}}& \cdots &{{{\left( {u_{08}^x} \right)}_1}}&{{{\left( {u_{01}^z} \right)}_1}}& \cdots &{{{\left( {u_{08}^z} \right)}_1}} \end{array}} \right)}\nonumber\\ &+& \left( {\begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{{\bf A}}_{11}}} &{{{{\bf A}}_{12}}} &{{{{\bf A}}_{13}}} &{{{{\bf A}}_{14}}}\\ {{{{\bf A}}_{21}}} &{{{{\bf A}}_{22}}} &{{{{\bf A}}_{23}}} &{{{{\bf A}}_{24}}} \end{array}} \right) \times {{{\bf H}}_2} \times {{\left( { \begin{array}{c{@}{\quad }c{@}{\quad }c{@}{\quad }c{@}{\quad }c{@}{\quad }c} {{{\left( {u_{01}^x} \right)}_2}} &\cdots &{{{\left( {u_{08}^x} \right)}_2}}& {{{\left( {u_{01}^z} \right)}_2}} &\cdots &{{{\left( {u_{08}^z} \right)}_2}} \end{array}} \right)}^{\rm{T}}}, \end{eqnarray} (22)where the two factor matrices H1 and H2 that take the forms of H1 = h1I and H2 = h2I are assigned to the relative displacements between the central particle and its neighbouring particles in the inner layer, that is, $${( {u_{0i}^x} )_1}$$ and those between the central particle and the newly added particles, that is, $${( {u_{0i}^x} )_2}$$, respectively. To find the ideal factors h1 and h2, we analyse the spatial dispersion of the extended lattice unit. This numerical dispersion can be written as   $$v'\left( {a{{\bf k}}} \right) = \sqrt { - \displaystyle\frac{1}{\rho }\left( {{{g'}_{11}} + {{g'}_{22}} \pm \sqrt {{{\left( {{{g'}_{11}} - {{g'}_{22}}} \right)}^2} + 4{{g'}_{12}}{{g'}_{21}}} } \right)} ,$$ (23)where   $${g'_{11}} = {{{\bf A}}_{11}} \times {{\bf L'}},$$ (24a)  $${g'_{12}} = {g'_{21}} = {{{\bf A}}_{12}} \times {{\bf L'}},$$ (24b)and   $${g'_{22}} = {{{\bf A}}_{23}} \times {{\bf L'}}.$$ (24c) Figure 3. View largeDownload slide Relative errors of P- and S-wave phase velocities for the basic square lattice unit in terms of lattice spacing a and the wavenumber components kx and kz. Figure 3. View largeDownload slide Relative errors of P- and S-wave phase velocities for the basic square lattice unit in terms of lattice spacing a and the wavenumber components kx and kz. Figure 4. View largeDownload slide (a) A second-order lattice unit. (b) and (c) represent the relative errors of P- and S-wave phase velocities, respectively, in terms of lattice spacing a and the wavenumber components kx and kz. Figure 4. View largeDownload slide (a) A second-order lattice unit. (b) and (c) represent the relative errors of P- and S-wave phase velocities, respectively, in terms of lattice spacing a and the wavenumber components kx and kz. For this second-order lattice unit, we have   $${{\bf L'}} = {h_1}\left( {\begin{array}{c} {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right) + {h_2}\left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\cos ( {2a{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {2\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {2a{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {2\sqrt 2 a{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right).$$ (25) To maintain the elastic properties of the lattice unit, when |ak| approximates zero, the phase velocities of both P and S waves for the second-order lattice unit should be the same as those for the basic square unit. Therefore,   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} f( {{{\bf L'}}} ) = f\left( {\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L'}}( {a{{\bf k}}} )} \right) = \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} v'( {a{{\bf k}}} ) = \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} f( {{\bf L}} ) = f\left( {\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}}( {a{{\bf k}}} )} \right);$$ (26)and the restriction on L' is given by   $$\mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L'}} = \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}} = \left( {\begin{array}{@{}*{1}{c}@{}} { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _1}}\\ { - {{\cos }^2}{\theta _2}}\\ { - \displaystyle\frac{1}{2}{{\cos }^2}{\theta _3}}\\ { - {{\cos }^2}{\theta _4}} \end{array}} \right).$$ (27) Carrying out the Taylor expansion on L' at ak = 0 yields   $${{\bf L'}} = \left( {\begin{array}{@{}*{1}{c}@{}} { - \displaystyle\frac{{\left( {{h_1} + 4{h_2}} \right)}}{2}{{\cos }^2}{\theta _1}}\\ { - \left( {{h_1} + 4{h_2}} \right){{\cos }^2}{\theta _2}}\\ { - \displaystyle\frac{{\left( {{h_1} + 4{h_2}} \right)}}{2}{{\cos }^2}{\theta _3}}\\ { - \left( {{h_1} + 4{h_2}} \right){{\cos }^2}{\theta _4}} \end{array}} \right) + \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{{24}}{{\cos }^4}{\theta _1} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{6}{{\cos }^4}{\theta _2} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{{24}}{{\cos }^4}{\theta _3} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{6}{{\cos }^4}{\theta _4} \cdot {{\left| {a{{\bf k}}} \right|}^2}} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)} \end{array}} \right).$$ (28) Substituting eq. (28) into eq. (27) defines the restriction on h1 and h2 as   \begin{eqnarray}{h_1} + 4{h_2} &=& 1.\end{eqnarray} (29) The extended lattice unit will converge to the corresponding homogeneous continua under this restriction. However, it still leaves these two factors with one degree of freedom. On the other hand, similar to eq. (19), the dispersion for the extended lattice is also obtained from the difference between L' and $$\mathop {\lim }\limits_{| {a{{\bf k}}} | \to 0} {{\bf L}}$$, that is,   $$d' = {{\bf L'}} - \mathop {\lim }\limits_{\left| {a{{\bf k}}} \right| \to 0} {{\bf L}}.$$ (30) Substituting eqs (28) and (29) into eq. (30) yields   $$d' = \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{{24}}{{\cos }^4}{\theta _1} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{6}{{\cos }^4}{\theta _2} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{{24}}{{\cos }^4}{\theta _3} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\left( {{h_1} + 16{h_2}} \right)}}{6}{{\cos }^4}{\theta _4} \cdot {{\left| {a{{\bf k}}} \right|}^2}} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^4}} \right)} \end{array}} \right).$$ (31) To reduce dispersion and speed up the convergence as |ak| approximates zero, the first term on the right-hand side of eq. (31) should be zero. This yields another restriction on h1 and h2 as   \begin{eqnarray}{h_1} + 16{h_2} &=& 0.\end{eqnarray} (32) Then, the magnitude of the dispersion term in eq. (31) is reduced to Order(|ak|4). Eqs (29) and (32) give the ideal relation between h1 and h2 that maintains the elastic properties of the lattice and reduces the numerical dispersion to Order(|ak|4). Fig. 4(b) shows the errors in phase velocities caused by the numerical dispersion in the extended second-order lattice unit. Applying this analysis to even higher order lattice units yields the general representation of the high-order DLM. As Fig. 5 illustrates, in the Nth-order lattice unit where the central particle is surrounded by N layers of particles, the interaction force acting on the central particles takes the form of   $$\left( {\begin{array}{@{}*{1}{c}@{}} {f_0^x}\\ {f_0^z} \end{array}} \right) = \sum\limits_{l = 1}^N {\left( {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} {{{{\bf A}}_{11}}}&{{{{\bf A}}_{12}}}&{{{{\bf A}}_{13}}}&{{{{\bf A}}_{14}}}\\ {{{{\bf A}}_{21}}}&{{{{\bf A}}_{22}}}&{{{{\bf A}}_{23}}}&{{{{\bf A}}_{24}}} \end{array}} \right) \times {{{\bf H}}_{\,l}} \times \left( {\begin{array}{*{20}{c}} {{{\left( {u_{01}^x} \right)}_l}}&\quad \cdots &\quad{{{\left( {u_{08}^x} \right)}_l}}&\quad{{{\left( {u_{01}^z} \right)}_l}}&\quad \cdots &\quad{{{\left( {u_{08}^z} \right)}_l}} \end{array}} \right)} ,$$ (33)where   $${{{\bf H}}_l} = {h_l}{{\bf I}},$$ (34) in which I is a 16 × 16 identical matrix; hl is the factor assigned to the lth layer and $${( {u_{0i}^x} )_l}$$ and $${( {u_{0i}^z} )_l}$$ represent the x- and z-components, respectively, of the difference in displacement between particle 0 and the corresponding particle in the lth layer that lies in the direction of r0i relative to particle 0. For this Nth-order lattice unit, we have   $${{{\bf L}}^{( N )}} = \sum\limits_{l = 1}^N {{h_l}\left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\cos ( {al{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 al{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {al{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 al{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right)} .$$ (35) Figure 5. View largeDownload slide The configuration of the Nth-order lattice unit. Figure 5. View largeDownload slide The configuration of the Nth-order lattice unit. Similar to the second-order lattice unit, for the Nth-order lattice unit, the dispersion also comes from the difference between L(N) and $$\mathop {\lim }\limits_{| {a{{\bf k}}} | \to 0} {{\bf L}}$$, that is, d(N). Substituting the Taylor expansion of L(N) at |ak| = 0 into this difference yields   $$\begin{array}{@{}c@{}} {d^{\left( N \right)}} = \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{1}{2}{{\cos }^2}{\theta _1} - \displaystyle\frac{{\sum\limits_{l = 1}^N {{l^2}{h_l}} }}{{2!}}{{\cos }^2}{\theta _1}}\\ {{{\cos }^2}{\theta _2} - \displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^2}{h_l}} }}{{2!}}{{\cos }^2}{\theta _2}}\\ {\displaystyle\frac{1}{2}{{\cos }^2}{\theta _3} - \displaystyle\frac{{\sum\limits_{l = 1}^N {{l^2}{h_l}} }}{{2!}}{{\cos }^2}{\theta _3}}\\ {{{\cos }^2}{\theta _4} - \displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^2}{h_l}} }}{{2!}}{{\cos }^2}{\theta _4}} \end{array}} \right) + \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\sum\limits_{l = 1}^N {{l^4}{h_l}} }}{{4!}}{{\cos }^4}{\theta _1} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^4}{h_l}} }}{{4!}}{{\cos }^4}{\theta _2} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{\sum\limits_{l = 1}^N {{l^4}{h_l}} }}{{4!}}{{\cos }^4}{\theta _3} \cdot {{\left| {a{{\bf k}}} \right|}^2}}\\ {\displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^4}{h_l}} }}{{4!}}{{\cos }^4}{\theta _4} \cdot {{\left| {a{{\bf k}}} \right|}^2}} \end{array}} \right) \cdots \cdots - \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} }}{{\left( {2N} \right)!}}{{\cos }^{2N}}{\theta _1} \cdot {{\left| {a{{\bf k}}} \right|}^{2N - 2}}}\\ {\displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} }}{{\left( {2N} \right)!}}{{\cos }^{2N}}{\theta _2} \cdot {{\left| {a{{\bf k}}} \right|}^{2N - 2}}}\\ {\displaystyle\frac{{\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} }}{{\left( {2N} \right)!}}{{\cos }^{2N}}{\theta _3} \cdot {{\left| {a{{\bf k}}} \right|}^{2N - 2}}}\\ {\displaystyle\frac{{2\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} }}{{\left( {2N} \right)!}}{{\cos }^{2N}}{\theta _4} \cdot {{\left| {a{{\bf k}}} \right|}^{2N - 2}}} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^{2N}}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^{2N}}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^{2N}}} \right)}\\ {{\rm{Order}}\left( {{{\left| {a{{\bf k}}} \right|}^{2N}}} \right)} \end{array}} \right). \end{array}$$ (36) To reduce the dispersion and maintain the elastic properties of the Nth-order lattice unit, hl is given by   $$\left\{ {\begin{array}{@{}*{1}{c}@{}} {\sum\limits_{l = 1}^N {{l^2}{h_l}} = 1}\\ {\sum\limits_{l = 1}^N {{l^4}{h_l}} = 0}\\ {\sum\limits_{l = 1}^N {{l^6}{h_l}} = 0}\\ \vdots \\ {\sum\limits_{l = 1}^N {{l^{2N}}{h_l}} = 0} \end{array}} \right..$$ (37) Under this relation, the dispersion is reduced to Order((ak)2N). This means that the high-order DLM using an Nth-order lattice unit will have 2Nth-order accuracy in terms of its lattice spacing. Based on the dispersion analysis, we can precisely choose the high-order lattice unit based on the accuracy requirement in terms of seismic wave simulation. Since the high-order DLM still focuses on the particle interactions that depend on the bonds that connect these particles, discontinuities in geological materials, such as fractures and faults, can also be interpreted by the breaking or weakening of these bonds. Therefore, the high-order DLM retains the advantage of using the discrete particle scheme in simulating complex geological materials, and it can be used to study both the static and dynamic deformations of the material. On the other hand, due to the particle feature of high-order DLMs, we form an absorbing pad by simply adding a friction force in the opposite direction of the particle velocity on each particle in the absorbing pad to gradually dampen their movements (see Appendix  B). In the aspect of 3-D seismic wave simulation, one can build a 3-D high-order lattice unit by combining three 2-D high-order lattice units in three orthogonal planes; then, the mathematical derivation of the 2-D high-order DLM can be easily implemented for the 3-D case. STABILITY CONDITION In regard to wave simulation in the time domain, due to the time iteration, the stability of the numerical simulation is also an important issue. We analyse the stability condition based on Von Neumann's method, which is described by Press et al. (1986). The numerical error is assumed as   $$\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _i^x\left( t \right)}\\ {\delta _i^z\left( t \right)} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {{U_x}\left( t \right)}\\ {{U_z}\left( t \right)} \end{array}} \right)\exp \left( {j{{\bf k}} \cdot {{\bf r}}_{0i}^0} \right),$$ (38)in which $$\delta _i^x( t )$$ and $$\delta _i^z( t )$$ are the x- and z-components of the numerical error for the displacement with respect to particle i, respectively. For this numerical error, the Verlet accumulation takes the form of   $$\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _i^x\left( {t + \Delta t} \right)}\\ {\delta _i^z\left( {t + \Delta t} \right)} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {\ddot{\delta }_i^x\left( t \right)}\\ {\ddot{\delta }_i^z\left( t \right)} \end{array}} \right)\Delta {t^2} + 2\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _i^x\left( t \right)}\\ {\delta _i^z\left( t \right)} \end{array}} \right) - \left( {\begin{array}{@{}*{1}{c}@{}} {\delta _i^x\left( {t - \Delta t} \right)}\\ {\delta _i^z\left( {t - \Delta t} \right)} \end{array}} \right).$$ (39) According to eq. (33), we have   $$\left( {\begin{array}{c} {\ddot{\delta }_0^x\left( t \right)}\\ {\ddot{\delta }_0^z\left( t \right)} \end{array}} \right) = \displaystyle\frac{{{\bf G}}}{\rho }\left( {\begin{array}{c} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right),$$ (40a)where   $${{\bf G}} = \left( {\begin{array}{c{@}{\quad }c} {g_{11}^{\left( N \right)}}&{g_{12}^{\left( N \right)}}\\ {g_{21}^{\left( N \right)}}&{g_{22}^{\left( N \right)}} \end{array}} \right),$$ (40b)in which   $$g_{11}^{\left( N \right)} = {{{\bf A}}_{11}} \times {{{\bf L}}^{\left( N \right)}},$$ (40c)  $$g_{12}^{\left( N \right)} = g_{21}^{\left( N \right)} = {{{\bf A}}_{12}} \times {{{\bf L}}^{\left( N \right)}},$$ (40d)and   $$g_{22}^{\left( N \right)} = {{{\bf A}}_{23}} \times {{{\bf L}}^{\left( N \right)}}.$$ (40e) For the Nth-order lattice unit,   $${{{\bf L}}^{( N )}} = \sum\limits_{l = 1}^N {{h_l} \left( {\begin{array}{@{}*{1}{c}@{}} {\displaystyle\frac{{\cos ( {la{{\bf k}} \cdot {{{\bf e}}_{01}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 la{{\bf k}} \cdot {{{\bf e}}_{02}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {la{{\bf k}} \cdot {{{\bf e}}_{03}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}}\\ {\displaystyle\frac{{\cos ( {\sqrt 2 la{{\bf k}} \cdot {{{\bf e}}_{04}}} ) - 1}}{{{a^2}{{{\bf k}}^2}}}} \end{array}} \right)} .$$ (40f) Substituting eq. (40) into eq. (39), we have   $$\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( {t + \Delta t} \right)}\\ {\delta _0^z\left( {t + \Delta t} \right)} \end{array}} \right) - 2\left( {{{\bf I}} - \displaystyle\frac{{\Delta {t^2}}}{{2\rho }}{{\bf G}}} \right)\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right) + \left( {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( {t - \Delta t} \right)}\\ {\delta _0^z\left( {t - \Delta t} \right)} \end{array}} \right) = 0,$$ (41) We define the growth rate of the numerical error by   $$\gamma = \displaystyle\frac{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( {t + \Delta t} \right)}\\ {\delta _0^z\left( {t + \Delta t} \right)} \end{array}} \right\|}}{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right\|}} = \displaystyle\frac{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right\|}}{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( {t - \Delta t} \right)}\\ {\delta _0^z\left( {t - \Delta t} \right)} \end{array}} \right\|}}.$$ (42) Therefore, eq. (41) can be written as   $${\gamma ^2} - 2\eta \gamma + 1 = 0,$$ (43a)where   $$\eta = \displaystyle\frac{{\left\| {\left( {{{\bf I}}{\rm{ + }}\displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}{{\bf G}}} \right)\left( {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right)} \right\|}}{{\left\| {\begin{array}{@{}*{1}{c}@{}} {\delta _0^x\left( t \right)}\\ {\delta _0^z\left( t \right)} \end{array}} \right\|}}.$$ (43b) The stability condition requires   $$\left| \gamma \right| \le 1.$$ (44)  This leads to the restriction on η given by   $$\left| \eta \right| \le 1.$$ (45)  Since   $$\left| \eta \right| \le \left| {1 + \displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}\mu } \right|,$$ (46)in which μ is the eigenvalue of matrix G, given by   $$\mu = \displaystyle\frac{1}{2}\left( {g_{11}^{\left( N \right)} + g_{22}^{\left( N \right)} \pm \sqrt {{{\left( {g_{11}^{\left( N \right)} - g_{22}^{\left( N \right)}} \right)}^2} + 4g_{12}^{\left( N \right)}g_{21}^{\left( N \right)}} } \right),$$ (47)eq. (45) can be rewritten as   $$- 2 \le \displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}\mu \le 0.$$ (48) The maximum values of μ are reached when   $${{\bf k}} = \left( {\begin{array}{@{}*{1}{c}@{}} {{k_x}}\\ {{k_z}} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {{\pi \mathord{\left/ {\vphantom {\pi a}} \right.} a}}\\ 0 \end{array}} \right),$$ (49a)  $${{\bf k}} = \left( {\begin{array}{@{}*{1}{c}@{}} {{k_x}}\\ {{k_z}} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} 0\\ {{\pi \mathord{\left/ {\vphantom {\pi a}} \right.} a}} \end{array}} \right),$$ (49b)or   $${{\bf k}} = \left( {\begin{array}{@{}*{1}{c}@{}} {{k_x}}\\ {{k_z}} \end{array}} \right) = \left( {\begin{array}{@{}*{1}{c}@{}} {{\pi \mathord{\left/ {\vphantom {\pi a}} \right.} a}}\\ {{\pi \mathord{\left/ {\vphantom {\pi a}} \right.} a}} \end{array}} \right).$$ (49c) Substituting eq. (49) into eq. (48) yields the general stability condition for the Nth-order lattice given by   $$0 \le \displaystyle\frac{{\Delta {t^2}}}{{\rho {a^2}}}\left( {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{{a^2}}} + {k_1}} \right)\sum\limits_{m = 1}^{{N \mathord{\left/ {\vphantom {N 2}} \right.} 2}} {{h_{2m - 1}}} \le 1,$$ (50a)  $$0 \le \displaystyle\frac{{\Delta {t^2}}}{{\rho {a^2}}}\left( {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{{a^2}}} + {k_3}} \right)\sum\limits_{m = 1}^{{N \mathord{\left/ {\vphantom {N 2}} \right.} 2}} {{h_{2m - 1}}} \le 1,$$ (50b)  $$0 \le \displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}\left( {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{{a^2}}} + {k_1} + {k_2} + {k_4} \pm \sqrt {k_1^2 + {{\left( {\displaystyle\frac{{{b_{23}} - {b_{12}}}}{{{a^2}}} + {k_4} - {k_2}} \right)}^2}} } \right)\sum\limits_{m = 1}^{{N \mathord{\left/ {\vphantom {N 2}} \right.} 2}} {{h_{2m - 1}}} \le 1,$$ (50c)and   $$0 \le \displaystyle\frac{{\Delta {t^2}}}{{2\rho {a^2}}}\left( {\displaystyle\frac{{{b_{12}} + {b_{23}}}}{{{a^2}}} + {k_3} + {k_2} + {k_4} \pm \sqrt {k_3^2 + {{\left( {\displaystyle\frac{{{b_{23}} - {b_{12}}}}{{{a^2}}} + {k_4} - {k_2}} \right)}^2}} } \right)\sum\limits_{m = 1}^{{N \mathord{\left/ {\vphantom {N 2}} \right.} 2}} {{h_{2m - 1}}} \le 1,$$ (50d)where $${N \mathord{/ {\vphantom {N 2}}} 2}$$ represents the division of the integer N by 2. The stability analysis shows that the high-order dynamic lattice not only reduces the numerical dispersion but also changes the stability condition in seismic wave simulation. NUMERICAL EXAMPLES Numerical tests are implemented to validate the high-order DLM. We apply the first-, second-, third- and fourth-order DLMs in elastic wave simulation on the TI media shown in Fig. 6, where we load a single force source in the vertical direction at x = 2.25 km and z = 2.25 km to generate elastic waves, and use the absorbing boundary condition (see Appendix  B) at all boundaries. The locations of the receivers are listed in Table 1. To compare the spatial dispersions of different DLMs, we set the dominant frequency of the Ricker source wavelet as 50 Hz and the lattice spacing as 5 m.The lattice network for this simulation includes 900 × 900 particles. The time interval is 0.5 ms, and the simulation goes through 1300 time steps. Table 2 shows the CPU times for the different DLMs when simulating the elastic wave propagation. The snapshots of the x-component of the displacement fields from the DLMs are plotted in Fig. 7. These snapshots reveal that the spatial dispersion decreases as the order of the lattice unit increases. We also compare the seismograms from the high-order DLMs with those from the DRP (dispersion relation preserving)/opt MacCormack FDM (Zhu et al. 2009; Zhang et al. 2012) in Fig. 8. This particular FDM solves the first-order partial differential velocity–stress equations with low dispersion or dissipation; it is non-staggered and fourth-order accurate in space. The seismogram comparison is also conducted on the model given in Fig. 6. The lattice spacing for the DLMs is 1 m, whereas the grid spacing for the FDM is 0.5 m. The dominant frequency of the Ricker wavelet is 20 Hz, and the time interval, Δt, is 0.1 ms. Since the maximum frequency of the Ricker wavelet is estimated to be 2.5 times the dominant frequency, fc, the relation between the grid spacing and point-per-wavelength (PPW) is given by   $$\Delta h \approx \displaystyle\frac{{{V_S}}}{{2.5{f_c}{\rm{PPW}}}},$$ (51)where Δh is the grid spacing for the FDM or the lattice spacing for the DLM, fc is the dominant frequency and VS is the shear wave velocity. For the TTI media, we substitute VS with quasi-S-wave velocity, VS0, which, according to Thomsen (1986), is given by   $${V_{S{\rm{0}}}} = \sqrt {\displaystyle\frac{{{c_{44}}}}{\rho }} ,$$ (52)where c44 is the component of the elasticity tensor, and ρ represents the density. Therefore, the PPW for the DLM simulations is approximately 33, whereas the PPW for the FDM simulation is approximately 66. Seismograms are recorded at receivers 1–4. The comparison indicates that when the lattice spacing is much smaller than the wavelength, the seismograms from different DLMs and those from the FDM coincide with each other. On the other hand, it can also be witnessed that the differences between the high-order DLMs and the FDM are smaller than those between the first-order DLM and the FDM. Figure 6. View largeDownload slide Transversely isotropic (TI) model used for the seismic wave simulations. Figure 6. View largeDownload slide Transversely isotropic (TI) model used for the seismic wave simulations. Figure 7. View largeDownload slide Snapshots of the x-components of the displacement fields generated by different DLMs. Parts of the displacement fields in the red frames are enlarged for detailed comparison. Figure 7. View largeDownload slide Snapshots of the x-components of the displacement fields generated by different DLMs. Parts of the displacement fields in the red frames are enlarged for detailed comparison. Figure 8. View largeDownload slide Seismograms recorded at receivers 1–4. The differences (black lines) between the seismograms from the DLMs (red solid lines) and those from the DRP/opt MacCormack FDM (blue dashed lines) are amplified to be significant. Figure 8. View largeDownload slide Seismograms recorded at receivers 1–4. The differences (black lines) between the seismograms from the DLMs (red solid lines) and those from the DRP/opt MacCormack FDM (blue dashed lines) are amplified to be significant. Table 1. Locations of the receivers. Receiver  1  2  3  4  5  6  7  x (km)  2.70  2.60  2.45  2.20  2.25  2.49  2.49  z (km)  2.25  2.60  2.65  2.65  2.49  2.49  2.25  Receiver  1  2  3  4  5  6  7  x (km)  2.70  2.60  2.45  2.20  2.25  2.49  2.49  z (km)  2.25  2.60  2.65  2.65  2.49  2.49  2.25  View Large Table 2. CPU times for different DLMs Order  1  2  3  4  t (s)  248.0  593.9  857.6  1214.2  Order  1  2  3  4  t (s)  248.0  593.9  857.6  1214.2  View Large To investigate the spatial dispersion for different methods, we carry out elastic wave simulations on the same TTI model shown in Fig. 6 with different spacings using different DLMs The dominant frequency for the source wavelet is 30 Hz, and the time interval is 0.1 ms. Fig. 9 shows the vertical displacements recorded at receiver 6 from different numerical methods with lattice spacing varying from 2 to 8 m. The seismogram generated by the DRP/opt MacCormack FDM with a grid spacing of 0.5 m is used as a reference. The seismograms in Fig. 9 show the effect of changing the spacing on different methods. The first-order DLM appears to be most affected by spatial dispersion when the spacing varies. When the spacing gets bigger, the seismograms from high-order DLMs, especially those from the third- and fourth-order DLMs, still fit the referencing seismogram well. We also compare the convergence for different DLMs as the spacing decreases in Fig. 10, where the corresponding PPW is used as the second horizontal axis and the misfit error is given by   $${\rm{error}} = \displaystyle\frac{{\sum\limits_i^N {{{\left[ {S\left( {{t_i}} \right) - {S^0}\left( {{t_i}} \right)} \right]}^2}} }}{{\sum\limits_i^N {{{\left[ {{S^0}\left( {{t_i}} \right)} \right]}^2}} }},$$ (53)in which S(ti) is the seismogram generated by a certain numerical scheme with a given spacing from a certain receiver, and S0(ti) is the referencing seismogram generated by the DRP/opt FDM with a grid spacing of 0.5 m from the same receiver. Through this comparison, we can clearly see that the seismograms from the high-order DLMs quickly converge to the referencing seismogram and have less misfit error at a large lattice spacing than the first-order DLM. Figure 9. View largeDownload slide Seismograms from DLMs when the spacing is 2, 5 and 8 m (from top to bottom). The reference seismogram is obtained by the DRP/opt MacCormack FDM with a grid spacing of 0.5 m. Figure 9. View largeDownload slide Seismograms from DLMs when the spacing is 2, 5 and 8 m (from top to bottom). The reference seismogram is obtained by the DRP/opt MacCormack FDM with a grid spacing of 0.5 m. Figure 10. View largeDownload slide Convergence curves for the DLMs in terms of the misfit error with respect to the seismograms recorded at different receivers with decreasing spacing. Figure 10. View largeDownload slide Convergence curves for the DLMs in terms of the misfit error with respect to the seismograms recorded at different receivers with decreasing spacing. We test different DLMs on a more realistic model shown in Fig. 11. A single force source in the vertical direction is located at x = 2 km and z = 0.6 km. The dominant frequency of the Ricker wavelet is 20 Hz, and the time interval is 0.5625 ms. The spacing is set to be 2.5 m. Receivers are arranged in three seismic profiles. The locations of the receivers are listed in Tables 3 –5. Snapshots of the horizontal displacements obtained by different methods are shown in Fig. 12. Detailed seismogram comparisons are shown in Figs 13 and 14. The results obtained from different methods coincide, demonstrating the ability of the high-order DLM to handle realistic media. Figure 11. View largeDownload slide The TTI overthrust model. The elastic properties of the TTI overthrust are represented by the Thomsen moduli. Each of the four thrust blocks has its own tilted angle. The location of the single force source is marked by the red star. Receivers are arranged in three profiles to record the seismograms. Figure 11. View largeDownload slide The TTI overthrust model. The elastic properties of the TTI overthrust are represented by the Thomsen moduli. Each of the four thrust blocks has its own tilted angle. The location of the single force source is marked by the red star. Receivers are arranged in three profiles to record the seismograms. Figure 12. View largeDownload slide Snapshots of the horizontal displacement wavefields simulated by the DRP/opt MacCormack FDM and different DLMs. Figure 12. View largeDownload slide Snapshots of the horizontal displacement wavefields simulated by the DRP/opt MacCormack FDM and different DLMs. Figure 13. View largeDownload slide Seismogram comparison of horizontal displacements between different DLMs (red solid lines) and the DRP/opt MacCormack FDM (blue dashed lines). Records from the receivers in profiles 1, 2 and 3 are plotted in the first, second and third columns, respectively. Amplified differences between these two kinds of methods are represented by the black dotted lines. Figure 13. View largeDownload slide Seismogram comparison of horizontal displacements between different DLMs (red solid lines) and the DRP/opt MacCormack FDM (blue dashed lines). Records from the receivers in profiles 1, 2 and 3 are plotted in the first, second and third columns, respectively. Amplified differences between these two kinds of methods are represented by the black dotted lines. Figure 14. View largeDownload slide Seismogram comparison of vertical displacements between different DLMs (red solid lines) and the DRP/opt MacCormack FDM (blue dashed lines). Records from the receivers in profiles 1, 2 and 3 are plotted in the first, second and third columns, respectively. Amplified differences between the DLMs and the FDM are represented by the black dotted lines. Figure 14. View largeDownload slide Seismogram comparison of vertical displacements between different DLMs (red solid lines) and the DRP/opt MacCormack FDM (blue dashed lines). Records from the receivers in profiles 1, 2 and 3 are plotted in the first, second and third columns, respectively. Amplified differences between the DLMs and the FDM are represented by the black dotted lines. Table 3. Locations of the receivers in seismic profile 1. Receiver  1  2  3  4  5  6  x (km)  1.80  1.80  1.80  1.80  1.80  1.80  z (km)  0.20  0.35  0.50  0.65  0.80  0.95  Receiver  1  2  3  4  5  6  x (km)  1.80  1.80  1.80  1.80  1.80  1.80  z (km)  0.20  0.35  0.50  0.65  0.80  0.95  View Large Table 4. Locations of the receivers in seismic profile 2. Receiver  1  2  3  4  5  6  x (km)  2.10  2.25  2.40  2.55  2.70  2.85  z (km)  0.80  0.80  0.80  0.80  0.80  0.80  Receiver  1  2  3  4  5  6  x (km)  2.10  2.25  2.40  2.55  2.70  2.85  z (km)  0.80  0.80  0.80  0.80  0.80  0.80  View Large Table 5. Locations of the receivers in seismic profile 3. Receiver  1  2  3  4  5  6  x (km)  3.40  3.40  3.40  3.40  3.40  3.40  z (km)  0.20  0.35  0.50  0.65  0.80  0.95  Receiver  1  2  3  4  5  6  x (km)  3.40  3.40  3.40  3.40  3.40  3.40  z (km)  0.20  0.35  0.50  0.65  0.80  0.95  View Large On the other hand, it is important to discuss the computational cost of this high-order numerical scheme. However, the exact computational cost depends on many factors, such as the efficiency of the code, the architecture of the hardware and the compiler. Therefore, discussions of the computational cost of the high-order DLM depend on the amount of required variables and multiplications because they can be regarded as the theoretical measurements of the computational cost that eliminate other influences, such as coding and hardware. For a n × n lattice network that contains n2 particles, the Nth-order DLM needs 32N × n2 variables to store the coefficients of the linear combination that are used in eq. (33) to calculate the interaction force on each particle, and another 6n2 variables are needed to store the horizontal and vertical components of the particle displacements in three consecutive time steps. In terms of the CPU time consumption, the high-order DLM spends most of the CPU time calculating the interaction force acting upon the particles through eq. (33). For the Nth-order DLM, it takes 32N × n2 multiplications to obtain the interaction forces on every particle at each time step. Therefore, the CPU time spent calculating the interaction force for the Nth-order DLM is N times that of the first-order DLM. This trend is witnessed in the CPU times for different DLMs listed in Table 2. Because the high-order DLM significantly reduces the dispersion, it allows us to use a larger lattice spacing. This will considerably reduce the particle density as well as the computational cost. CONCLUSIONS In an effort to improve the accuracy and suppress the spatial dispersion of the original DLM, we analyse the properties obtained using the different layouts of lattice units. Due to the special arrangements of the crisscross- and x-shaped lattice units, these two types of basic lattice units are limited to certain elastic media. As a result, the basic square lattice is used to develop a high-order DLM. Our research indicates that the spatial dispersion decreases as we add more particles to the lattice unit. Based on the numerical dispersion analysis, we build the theory for a high-order DLM. According to our theory, the accuracy of the high-order DLM depends on the order of the lattice unit. The higher the order of the lattice unit is, the more layers of particles are included in the lattice unit. The high-order DLM using the Nth-order lattice unit is of 2Nth-order accuracy in terms of lattice spacing. On the other hand, to better elaborate the high-order DLM, we also derive the stability condition for the method. Numerical tests show that our high-order DLM significantly reduces the spatial dispersion and improves the accuracy as well as the convergence rate for seismic wave simulation in TI media. The detailed comparison of the accuracy between the referencing FDM and different DLMs validates the theory of the high-order DLM. Based on our theory, one can precisely choose the order of the DLM based on the accuracy requirement. The high-order DLM allows us to use a larger lattice spacing, which can greatly reduce the computational cost. This improvement will greatly benefit the application of the DLM in seismic wave simulations. ACKNOWLEDGEMENTS The authors thank W. Zhang for kindly providing the FD program and gratefully acknowledge X. Qiang, Q.-L. Li, Y. Jiang, L. Yang, C. Hu and Q.-H. Li for their fruitful discussions. This study received support from National Key Research and Development Program of China (no. 2017YFB0202903) and the National Natural Science Foundation of China (41774121 and 41374006). REFERENCES Aki K., Richards P.G., 2002. Quantitative Seismology , University Science Books. Cundall P.A., Strack O.D., 1979. A discrete numerical model for granular assemblies, Géotechnique , 29, 47– 65. https://doi.org/10.1680/geot.1979.29.1.47 Google Scholar CrossRef Search ADS   Del Valle-García R., Sánchez-Sesma F.J., 2003. Rayleigh waves modeling using an elastic lattice model, Geophys. Res. Lett. , 30, 1866, doi:10.1029/2003GL017600. https://doi.org/10.1029/2003GL017600 Google Scholar CrossRef Search ADS   Hixon R., 2000. Prefactored small-stencil compact schemes, J. Comput. Phys. , 165, 522– 541. https://doi.org/10.1006/jcph.2000.6631 Google Scholar CrossRef Search ADS   Hoover W.G., Ashurst W.T., Olness R.J., 1974. Two-dimensional computer studies of crystal stability and fluid viscosity, J. Chem. Phys. , 60, 4043– 4047. https://doi.org/10.1063/1.1680855 Google Scholar CrossRef Search ADS   Hu X., Jia X., 2016. A dynamic lattice method for elastic seismic modeling in anisotropic media, Geophysics , 81( 3), T131– T143. https://doi.org/10.1190/geo2015-0511.1 Google Scholar CrossRef Search ADS   Kelly K.R., Ward R.W., Treitel S., Alford R.M., 1976. Synthetic seismograms—a finite difference approach, Geophysics , 41, 2– 27. https://doi.org/10.1190/1.1440605 Google Scholar CrossRef Search ADS   Levander A.R., 1988. Fourth-order finite-difference P-SV seismograms, Geophysics , 53( 11), 1425– 1436. https://doi.org/10.1190/1.1442422 Google Scholar CrossRef Search ADS   Lubbe R., Worthington M.H., 2006. A field investigation of fracture compliance, Geophys. Prospect. , 54( 3), 319– 331. https://doi.org/10.1111/j.1365-2478.2006.00530.x Google Scholar CrossRef Search ADS   Madariaga R., 1976. Dynamics of an expanding circular fault, Bull. seism. Soc. Am. , 66, 163– 182. Moczo P., Kristek J., Vavrycuk V., Archuleta R., Halda L., 2002. 3D heterogeneous staggered-grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities, Bull. seism. Soc. Am. , 92, 3042– 3066. https://doi.org/10.1785/0120010167 Google Scholar CrossRef Search ADS   Monette L., Anderson M.P., 1994. Elastic and fracture properties of the two-dimensional triangular and square lattices, Model. Simul. Mater. Sci. Eng. , 2, 53– 66. https://doi.org/10.1088/0965-0393/2/1/004 Google Scholar CrossRef Search ADS   Mora P., Place D., 1998. Numerical simulation of earthquake faults with gouge: toward a comprehensive explanation for the heat flow paradox, J. geophys. Res. , 103, 21 067– 21 089. https://doi.org/10.1029/98JB01490 Google Scholar CrossRef Search ADS   Möllhoff M., Bean C.J., 2008. Validation of elastic wave measurements of rock fracture compliance using numerical discrete particle simulations, Geophys. Prospect. , 7( 5), 883– 895. O'Brien G.S., Bean C.J., 2004. A 3D discrete numerical elastic lattice method for seismic wave propagation in heterogeneous media with topography, Geophys. Res. Lett. , 31, L14608, doi:10.1029/2004GL020069. https://doi.org/10.1029/2004GL020069 Google Scholar CrossRef Search ADS   O'Brien G.S., Bean C.J., 2009. Volcano topography, structure and intrinsic attenuation: their relative influences on a simulated 3D visco-elastic wavefield, J. Volc. Geotherm. Res. , 183, 122– 136. https://doi.org/10.1016/j.jvolgeores.2009.03.004 Google Scholar CrossRef Search ADS   O'Brien G.S, Bean C.J., 2011. An irregular lattice method for elastic wave propagation, Geophys. J. Int. , 187( 3), 1699– 1707. https://doi.org/10.1111/j.1365-246X.2011.05229.x Google Scholar CrossRef Search ADS   O'Brien G.S., Bean C.J., Tapamo H., 2009. Dispersion analysis and computational efficiency of elastic lattice methods for seismic wave propagation, Comput. Geosci. , 35, 1768– 1775. https://doi.org/10.1016/j.cageo.2008.12.004 Google Scholar CrossRef Search ADS   Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1986. Numerical Recipes in FORTRAN: The Art of Scientific Computing , Cambridge Univ. Press. Saenger E.H., Bohlen T., 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid, Geophysics , 69, 583– 591. https://doi.org/10.1190/1.1707078 Google Scholar CrossRef Search ADS   Saenger E.H., Gold N., Shapiro S.A., 2000. Modeling the propagation of elastic waves using a modified finite-difference grid, Wave Motion , 31, 77– 92. https://doi.org/10.1016/S0165-2125(99)00023-2 Google Scholar CrossRef Search ADS   Saltzer S.D., Pollard D.D., 1992. Distinct element modeling of structures formed in sedimentary overburden by extensional reactivation of basement normal faults, Tectonics , 11, 165– 174. https://doi.org/10.1029/91TC02462 Google Scholar CrossRef Search ADS   Thomsen L., 1986. Weak elastic anisotropy, Geophysics , 51, 1954– 1966. Google Scholar CrossRef Search ADS   Toomey A., Bean C.J., 2000. Numerical simulation of seismic waves using a discrete particle scheme, Geophys. J. Int. , 141, 595– 604. https://doi.org/10.1046/j.1365-246x.2000.00094.x Google Scholar CrossRef Search ADS   Toomey A., Bean C.J., Scotti O., 2002. Fracture properties from seismic data—a numerical investigation, Geophys. Res. Lett. , 29 4, 1050, doi:10.1029/2001GL013867. https://doi.org/10.1029/2001GL013867 Google Scholar CrossRef Search ADS   Virieux J., 1984. SH-wave propagation in heterogeneous media—velocity-stress finite-difference method, Geophysics , 49, 1933– 1942. https://doi.org/10.1190/1.1441605 Google Scholar CrossRef Search ADS   Virieux J., 1986. P-SV wave propagation in heterogeneous media—velocity-stress finite-difference method, Geophysics , 51, 889– 901. https://doi.org/10.1190/1.1442147 Google Scholar CrossRef Search ADS   Zhang W., Shen Y., Zhao L., 2012. Three-dimensional anisotropic seismic wave modelling in spherical coordinates by a collocated-grid finite-difference method, Geophys. J. Int. , 188, 1359– 1381. https://doi.org/10.1111/j.1365-246X.2011.05331.x Google Scholar CrossRef Search ADS   Zhu H.J., Zhang W., Chen X.F., 2009. Two-dimensional seismic wave simulation in anisotropic media by non-staggered finite difference method (in Chinese), Chin. J. Geophys. , 52, 1536– 1546. Google Scholar CrossRef Search ADS   APPENDIX A: DERIVATION OF THE RELATION BETWEEN STIFFNESS COEFFICIENTS AND THE ELASTICITY TENSOR The elastic energy in the distorted lattice unit is defined by eq. (1). Assuming θij ≈ θ0, that is, the distortion is relatively much smaller than the lattice spacing, yields   $$\cos {\theta _{ij}} - \cos {\theta _0} \approx - \sin {\theta _0}\left( {{\theta _i} - {\theta _j}} \right),$$ (A1)where θi and θj are the angles by which the bonds 0–i and 0–j rotate from their original positions, respectively. θi can be expressed as   $${\theta _i} = \displaystyle\frac{{{{{\bf u}}_{0i}} \cdot {{\bf s}}_{0i}^0}}{{{{\left| {{{\bf r}}_{0i}^0} \right|}^2}}},$$ (A2)in which u0i is the difference in displacement between particle i and 0, $${{\bf s}}_{0i}^0$$ is the vector perpendicular to $${{\bf r}}_{0i}^0$$ in a right-handed system and $$| {{{\bf s}}_{0i}^0} | = | {{{\bf r}}_{0i}^0} |$$. On the other hand, when u0i is much smaller than $${{\bf r}}_{0i}^0$$, we also have   $$\left| {{{{\bf r}}_{0i}}} \right| - \left| {{{\bf r}}_{0i}^0} \right| = \displaystyle\frac{{{{{\bf u}}_{0i}} \cdot {{\bf r}}_{0i}^0}}{{\left| {{{\bf r}}_{0i}^0} \right|}}.$$ (A3) Thus, for the elastic energy given by eq. (1), we have   $$E \approx \displaystyle\frac{1}{2}{\sum\limits_{i = 1}^{n/2} {{k_i}\left( {\displaystyle\frac{{{{{\bf u}}_{0i}} \cdot {{\bf r}}_{0i}^0}}{{\left| {{{\bf r}}_{0i}^0} \right|}}} \right)} ^2} + \displaystyle\frac{1}{2}{\sum\limits_{i = 1,j = {\bmod}\!({i,n}){\rm{ + }}1}^n {{b_{ij}}{{\sin }^2}{\theta _0}\left( {\displaystyle\frac{{{{{\bf u}}_{0i}} \cdot {{\bf s}}_{0i}^0}}{{{{\left| {{{\bf r}}_{0i}^0} \right|}^2}}} - \displaystyle\frac{{{{{\bf u}}_{0j}} \cdot {{\bf s}}_{0j}^0}}{{{{\left| {{{\bf r}}_{0j}^0} \right|}^2}}}} \right)} ^2}.$$ (A4) To find the connection between the elastic energy in the lattice unit and the elastic energy in the continuum theory given by eq. (2), we substitute u0i with the strain tensor. According to the continuum theory, u0i can be expressed as   $${{{\bf u}}_{0i}} = \left( {\displaystyle\frac{{\partial {u_x}}}{{\partial x}}{{\left( {r_{0i}^0} \right)}_x} + \displaystyle\frac{{\partial {u_x}}}{{\partial z}}{{\left( {r_{0i}^0} \right)}_z}} \right){{{\bf e}}_x} + \left( {\displaystyle\frac{{\partial {u_z}}}{{\partial x}}{{\left( {r_{0i}^0} \right)}_x} + \displaystyle\frac{{\partial {u_z}}}{{\partial z}}{{\left( {r_{0i}^0} \right)}_z}} \right){{{\bf e}}_z},$$ (A5)where $${( {r_{0i}^0} )_x}$$ and $${( {r_{0i}^0} )_z}$$ are the x- and z-components of $${{\bf r}}_{0i}^0$$, respectively; meanwhile, the strain tensor is defined by   $${\varepsilon _{xx}} = \displaystyle\frac{{\partial {u_x}}}{{\partial x}},$$ (A6a)  $${\varepsilon _{zz}} = \displaystyle\frac{{\partial {u_z}}}{{\partial z}},$$ (A6b)and   $${\varepsilon _{xz}} = {\varepsilon _{zx}} = \displaystyle\frac{1}{2}\left( {\displaystyle\frac{{\partial {u_x}}}{{\partial z}} + \displaystyle\frac{{\partial {u_z}}}{{\partial x}}} \right),$$ (A6c)where $$\displaystyle\frac{{\partial {u_x}}}{{\partial x}}$$, $$\displaystyle\frac{{\partial {u_z}}}{{\partial z}}$$, $$\displaystyle\frac{{\partial {u_x}}}{{\partial z}}$$ and $$\displaystyle\frac{{\partial {u_z}}}{{\partial x}}$$ are the first-order derivatives of the displacement field. For the crisscrossed lattice unit (Fig. 1a), the x- and z-components of $${{\bf r}}_{0i}^0$$ are given by   $$\left( {\begin{array}{c@{\quad}c} {{{\left( {r_{01}^0} \right)}_x}}&{{{\left( {r_{01}^0} \right)}_z}}\\ {{{\left( {r_{02}^0} \right)}_x}}&{{{\left( {r_{02}^0} \right)}_z}}\\ {{{\left( {r_{03}^0} \right)}_x}}&{{{\left( {r_{03}^0} \right)}_z}}\\ {{{\left( {r_{04}^0} \right)}_x}}&{{{\left( {r_{04}^0} \right)}_z}} \end{array}} \right) = \left( {\begin{array}{c@{\quad}c} a&0\\ 0&a\\ { - a}&0\\ 0&{ - a} \end{array}} \right),$$ (A7)where a is the lattice spacing. Substituting eqs (A7), (A6) and (A5) into eq. (A4) yields   $$E = \displaystyle\frac{1}{2}{k_1}{a^2}{\varepsilon _{xx}}^2 + \displaystyle\frac{1}{2}{k_2}{a^2}{\varepsilon _{zz}}^2 + 2\left( {{b_{12}} + {b_{23}} + {b_{34}} + {b_{41}}} \right){\varepsilon _{xz}}^2.$$ (A8) Eq. (3) is obtained by comparing eq. (A8) with eq. (2) through   $$E = {a^2}\varphi ,$$ (A9)where φ is the elastic energy density given by eq. (2). In the case of the x-shaped lattice unit (Fig. 1b), the x- and z-components of $${{\bf r}}_{0i}^0$$ are given by   $$\left( {\begin{array}{c@{\quad}c} {{{\left( {r_{01}^0} \right)}_x}}&{{{\left( {r_{01}^0} \right)}_z}}\\ {{{\left( {r_{02}^0} \right)}_x}}&{{{\left( {r_{02}^0} \right)}_z}}\\ {{{\left( {r_{03}^0} \right)}_x}}&{{{\left( {r_{03}^0} \right)}_z}}\\ {{{\left( {r_{04}^0} \right)}_x}}&{{{\left( {r_{04}^0} \right)}_z}} \end{array}} \right) = \left( {\begin{array}{c@{\quad}c} {\displaystyle\frac{{\sqrt 2 }}{2}a}&{\displaystyle\frac{{\sqrt 2 }}{2}a}\\ { - \displaystyle\frac{{\sqrt 2 }}{2}a}&{\displaystyle\frac{{\sqrt 2 }}{2}a}\\ { - \displaystyle\frac{{\sqrt 2 }}{2}a}&{ - \displaystyle\frac{{\sqrt 2 }}{2}a}\\ {\displaystyle\frac{{\sqrt 2 }}{2}a}&{ - \displaystyle\frac{{\sqrt 2 }}{2}a} \end{array}} \right).$$ (A10) Similar to the crisscross-shaped lattice unit, the elastic energy for the x-shaped lattice unit is expressed as   $$\begin{array}{@{}*{3}{l}@{}} E& = &{\left( {{b_{12}} + {b_{23}} + {b_{34}} + {b_{41}} + \displaystyle\frac{{{a^2}{k_1}}}{4} + \displaystyle\frac{{{a^2}{k_2}}}{4}} \right)\left( {\displaystyle\frac{1}{2}{\varepsilon _{xx}}^2 + \displaystyle\frac{1}{2}{\varepsilon _{zz}}^2} \right)} + {\left( {\displaystyle\frac{{{a^2}{k_1}}}{4} + \displaystyle\frac{{{a^2}{k_2}}}{4} - {b_{12}} - {b_{23}} - {b_{34}} - {b_{41}}} \right){\varepsilon _{xx}}{\varepsilon _{zz}}}\\ && +\,{\displaystyle\frac{{{a^2}}}{2}\left( {{k_1} + {k_2}} \right){\varepsilon _{xz}}^2 + \displaystyle\frac{{{a^2}}}{2}\left( {{k_1} - {k_2}} \right)\left( {{\varepsilon _{xx}} + {\varepsilon _{zz}}} \right){\varepsilon _{xz}}} \end{array}.$$ (A11) Substituting eqs (2) and (A11) into eq. (A9) yields eq. (4). APPENDIX B: ABSORBING BOUNDARY CONDITION We use an absorbing pad to realize the absorbing boundary condition. Since the DLM is particle-based, we add a friction force in the opposite direction of the particle velocity on each particle in the absorbing pad. Therefore, the total force acting upon the particle in the absorbing pad is given by   $${{{\bf f}}_{{\rm{total}}}} = {{{\bf f}}^0} + {{{\bf f}}_{{\rm{friction}}}},$$ (B1)where ftotal is the total force, f0 is the interaction force given by eq. (33) and ffriction is the friction force, which takes the form of   $${{{\bf f}}_{{\rm{friction}}}} = - \mu {{\bf \dot{u}}},$$ (B2)where μ is the friction factor and $${{\bf \dot{u}}}$$ is the particle velocity. This friction force dampens the kinetic energy of the particle. In practice, we increase the friction factor as the absorbing pad expands to better attenuate the motions of the particles. Fig. B1 shows the effect of this absorbing pad. Figure B1. View largeDownload slide Snapshots of the wavefield at different time steps. Absorbing layers are added on the left and on the top of the area. The interfaces between the absorbing layers and the full elastic area are marked by the black lines. The boundary conditions on the right and bottom represent the fixed boundary conditions. Figure B1. View largeDownload slide Snapshots of the wavefield at different time steps. Absorbing layers are added on the left and on the top of the area. The interfaces between the absorbing layers and the full elastic area are marked by the black lines. The boundary conditions on the right and bottom represent the fixed boundary conditions. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations