A new lattice spring model for seismic wave simulation in heterogeneous anisotropic media

SUMMARY

Lattice spring model (LSM) provides an alternative numerical approach for simulating seismic wave propagation in heterogeneous media. This method has gained great popularity in fractured media due to its intuitive physical representation. Originating from the discrete element method, the LSM allows particles to achieve micromechanical interactions through springs rather than directly solving the differential equation. The most important issue in the LSM is calibrating the spring coefficients, which can be derived through experiments or physical principles. By simply removing the springs that exceed their strength, the LSM can easily simulate the entire failure process of materials, a task that is challenging for continuum-based methods such as the finite difference method (FDM) and finite element method. In this paper, we propose a new LSM for seismic wave simulation in heterogeneous anisotropic media, which yields more accurate results compared to the regular particle-based methods. Unlike the conventional LSM, which calibrates spring coefficients using the wave equation with an implicit homogeneous approximation, our new LSM calibrates the coefficients using a modified wave equation in heterogeneous media. Compared with the conventional LSM, whose spring coefficients only contain the elasticity tensor itself, our new LSM additionally takes the first derivative terms of the elasticity tensor into account, and thus can accurately handle the scattering waves in seismic wave simulation. We investigate the spring coefficients of the two LSMs and derive the numerical dispersion and stability condition. To validate the accuracy of the new LSM, we test several scattering, layered and complex heterogeneous anisotropic models, respectively, comparing their results with those obtained using the high-accuracy FDM. Numerical experiments demonstrate the high quality of the new LSM in complex media compared with the conventional LSM. Finally, two fracture models are simulated to illustrate the new LSM’s capability in modelling the complex failure process.

Seismic anisotropy, Wave propagation, Wave scattering and diffraction

1 INTRODUCTION

Seismic wave forward modelling is crucial in various geophysical applications, including fault rupture tracking, strong ground motion simulation, subsurface structure inversion and resource exploration. Various numerical tools are employed to simulate the wavefield propagation depending on the material’s rheology. The most commonly used methods include the finite difference method (FDM; Alterman & Karal Jr. 1968; Kelly et al. 1976; Virieux 1984, 1986; Graves 1996; Zhang & Chen 2006; Zhang & Zhang 2022), the finite element method (FEM; Lysmer & Drake 1972; Marfurt 1984; Zhu & Dorman 2000; Qu et al. 2022; Xu et al. 2022) and the spectral element method (SEM; Komatitsch & Vilotte 1998; Komatitsch & Tromp 1999; Tromp et al. 2008; Zhang et al. 2020). Based on continuum mechanics, these numerically discretized methods can offer accurate solutions of the wave equation under the assumption of smoothly varying media. Therefore, in the naturally occurring discontinuous features in actual media, such as faults and fractures, the model must be smoothed to satisfy the continuum differential equation. The discontinuous Galerkin mehod (DGM) allows for discontinuous wavefields and is adept at handling such cases (Virieux et al. 2011). Recently, Masson (2023) and Masson & Virieux (2023) proposed the distributional finite difference method (DFDM) that integrates some essential features of FDM, FEM, SEM and DGM. The DFDM enables discontinuous wavefields between neighbouring elements, making it a promising tool for capturing the discontinuous features (Lyu et al. 2024; Masson et al. 2024).

In the meantime, to avoid the continuum limitations of the differential equation, another category of discrete particle models has been developed in simulating dynamic fractures and fault propagation. Cundall & Strack (1979) adopted number of discs in the distinct element method to study the deformation and fracturing behaviours of granular assemblies. This idea was then applied to study material fractures (Monette & Anderson 1994), earthquake dynamics (Mora & Place 1998) and tectonic processes (Saltzer & Pollard 1992). The discrete element method (DEM) focuses on contact mechanics and particle interactions, excelling in scenarios involving granular and particulate materials. Traditional DEM models for continuum behaviour often require a calibration process to determine the correct interparticle parameters by experiments and involve time-consuming contact detections and interactions (Cheng et al. 2009). Nevertheless, their accuracy is typically lower than that of continuum-based methods. Consequently, combined methods such as DEM–FDM (Yin et al. 2020) have been studied, employing the DEM in severely damaged areas and the FDM in the continuum. However, the boundary coupling of these two methods introduces complexity. Therefore, the emergence of some special DEM models such as the lattice spring method (LSM), offering good accuracy comparable to traditional continuum-based numerical methods, has drawn significant attention. The LSM can be considered as a specialized form of DEM with a compact arrangement of particles connected by springs, which yields highly accurate results for wave propagation in continuum media (Liu & Liu 2006). For wave problems, the spring coefficients in the LSM can be calibrated using wave equations. In this regard, the LSM is similar to the FDM and the FEM, as all three approaches rely on the local approximation of constitutive laws. The LSM does not discretize the wave equation directly; instead, it discretizes the geological structure as discrete particles and attempt to investigate the microscopic interactions through springs. Therefore, the discrete particle methods present great flexibility for solving different wave equations. Moreover, inheriting some advantages from DEM, the LSM excels in capturing fracture mechanics. LSM has been widely applied in simulations for fracture initiation and propagation by controlling the springs break and removal (Ostoja-Starzewski 2002; Zhao et al. 2011; Chen et al. 2014, 2016; Zhao 2015; Buxton 2022).

In the field of seismic wave modelling, Toomey & Bean (2000) presented a particle-based model with hexagonal lattices to simulate elastic wave propagation in 2-D heterogeneous media. The Poisson's ratio was restricted to 0.25, as only the central force between a pair of particles was considered. Del Valle‐García & Sánchez‐Sesma (2003) added the bond-bending force to remove the restrictions of fixed Poisson's ratio and simulated Rayleigh waves using an elastic lattice model (ELM). Based on their work, O'Brien & Bean (2004) extended the ELM (known as LSM) to 3-D cases and then studied the effects of the volcanic topography. O'Brien et al. (2009) analysed the dispersion and computational efficiency of the LSM and added a new force term to reduce the numerical dispersion. O′Brien (2014) provided a new implementation for the free surface and improved the accuracy of modelling Rayleigh waves. Hu & Jia (2016) extended the LSM from isotropic media to transversely isotropic media and accelerated the calculation using a linear transformation scheme. Xia et al. (2017) compared three types of 3-D LSMs with different patterns of springs and studied their dispersion features. Xia et al. (2018) developed a rectangular-grid LSM to simulate elastic wave propagation, which allowed for easy and flexible mesh generation. Hu & Jia (2018) added more layers of particles to build high-order LSM schemes, significantly reducing the spatial dispersion and improving the accuracy. O'Brien (2021) extended the LSM to an isotropic nonlinear viscoelastic medium and studied the dynamic and static deformation. Tang et al. (2022) implemented the perfectly matched layer (PML) in the LSM for Poisson solids and reduced boundary reflections. Tang et al. (2023) proposed an adaptive location strategy in the LSM to improve the simulation accuracy for irregular interfaces. Liu & Fu (2021) combined the modified LSM with discrete fracture networks to study seismic responses in complex fractured media. The conventional LSMs are essentially based on an implicit homogeneous approximation, as the spring coefficients are calibrated using the wave equation in homogeneous media instead of heterogeneous media. As a result, it leads to some errors for the wavefields in heterogeneous media, especially in strongly varied media. Indeed, some authors have compared the influence of these two types of wave equations on the FDM (Di Bartolo et al. 2012; Li et al. 2019). In this paper, we calibrate spring coefficients using the modified wave equation in heterogeneous media and develop a new LSM for seismic wave simulation in heterogeneous anisotropic media. First, we derive the body force using Hooke's law for a network of springs and then calibrate the spring coefficients of the conventional LSM and the new LSM using the second-order displacement equation in the homogeneous and heterogeneous anisotropic media, respectively. Meanwhile, we theoretically analyse the differences of the spring coefficients between the conventional LSM and our new LSM. The dispersion property and stability condition are derived to determine proper parameters in simulations. By comparing three numerical results with those obtained by the high-accuracy FDM, the accuracy of the conventional LSM and our new LSM is examined. Numerical results show that, for the seismic wave simulation in the heterogeneous media, the conventional LSM produces errors for the scattering and the reflected waves; on the contrary, our new LSM has good accuracy when dealing with scatterings in complex media. Finally, two fractured models are employed to verify the new LSM's capability in simulating complex failure processes.

2 LATTICE SPRING MODEL

The LSM is a numerical method that provides a mesoscopic perspective for describing subsurface media. It employs a particle-based model to represent actual media, treating the medium as a collection of lattices with interconnected springs between discrete particles. By directly converting the medium into numerous particles with specific physical properties and interactions, the lattice model serves as a replacement for the medium. Seismic wave propagation can be simulated by analysing the force fields generated by the interconnected springs on particles.

2.1 Basic lattice units in the LSM

Similar to the molecular structure in chemistry, the LSM is comprised of numerous stacked lattice units. We employ the square D2Q8 lattice model, which implies a 2-D square lattice model with eight neighbour particles (Xia et al. 2018). Fig. 1(a) illustrates a D2Q8 basic unit of the LSM. Each particle can interact with its neighbours through springs. Fig. 1(b) depicts the basic types of springs used in the LSM. The linear springs between two adjacent particles provide central interactions and the angular springs among three adjacent particles provide shearing angular interactions. Correspondingly, k and b are introduced to represent the spring coefficients, which reflect the level of deformation difficulty when the springs are stretched or bent. Large coefficients imply that it is hard to stretch the linear spring or bend the angular spring. According to Monette & Anderson (1994), the elastic energy density in a distorted lattice model (Fig. 1c) can be expressed as

$$\begin{eqnarray} {\emptyset }_0 = \frac{1}{{2{S}_0{\rm{\ }}}}\mathop \sum \limits_{i = 1}^N {k}_{0i}{\left( {|{{\boldsymbol{r}}}_{0i}\left| - \right|{\boldsymbol{r}}_{0i}^0|} \right)}^2 + \frac{1}{{2{S}_0}}\mathop \sum \limits_{i = 1,{\rm{\ }}j = \mathrm{ mod}\left( {i,N} \right) + 1}^N {b}_{i0j}{\left( {\mathrm{ cos}{\theta }_{i0j} - \mathrm{ cos}\theta _{i0j}^0} \right)}^2, \end{eqnarray}$$

(1)

in which |${\emptyset }_0$| is the elastic energy density for the central particle, |${S}_0$| is the corresponding effective area that equals |${\rm{\Delta }}x{\rm{\Delta }}z$|⁠, and |${\rm{\Delta }}x$| and |${\rm{\Delta }}z$| are the intervals in the x- and z-directions, respectively. N is the total number of the neighbouring particles which equals 8 in Fig. 1. |$| * |$| denotes the absolute value of *. |${{\boldsymbol{r}}}_{0i}$| and |${\boldsymbol{r}}_{0i}^0$| are the vectors for the linear spring (referred to as spring 0-i for short) that connects particles 0 and i in the distorted and undistorted lattice, respectively. Similarly, |${\theta }_{i0j}$| and |$\theta _{i0j}^0$| are the angles for the angular spring (referred to as spring i-0-j for short) between the two linear springs 0-i and 0-j in the distorted and undistorted lattice, respectively. |${k}_{0i}$| and |${b}_{i0j}$| are the coefficients for the linear spring 0-i and angular spring i-0-j, respectively. The function mod(i, N) denotes the remainder function after |$\ i$| is divided by N. The first term on the right-hand side of eq. (1) represents the elastic energy density due to the two-point central interactions related to the stretching of the linear spring, and the second term represents the elastic energy density due to the three-point angular interactions associated with the bending of the angular springs. Therefore, the spring force acting on particle 0 in the distorted lattice (Fig. 1c) can be calculated by

$$\begin{eqnarray} {{\boldsymbol{F}}}_0 = - \frac{{\partial {\emptyset }_0{\rm{\ }}}}{{\partial {{\boldsymbol{u}}}_0}} = - \frac{1}{{{S}_0{\rm{\ }}}}\mathop \sum \limits_{i = 1}^8 {k}_{0i}\frac{{\partial (|{{\boldsymbol{r}}}_{0i}\left| - \right|{\boldsymbol{r}}_{0i}^0|)}}{{\partial {{\boldsymbol{u}}}_0}} - \frac{1}{{{S}_0}}\mathop \sum \limits_{i = 1,{\rm{\ }}j = \mathrm{ mod}\left( {i,8} \right) + 1}^8 {b}_{i0j}\left( {\mathrm{ cos}{\theta }_{i0j} - \mathrm{ cos}\theta _{i0j}^0} \right)\frac{{\partial \left( {\mathrm{ cos}{\theta }_{i0j} - \mathrm{ cos}\theta _{i0j}^0} \right)}}{{\partial {{\boldsymbol{u}}}_0}}, \end{eqnarray}$$

(2)

where |${{\boldsymbol{u}}}_{0i} = {{\boldsymbol{u}}}_i - {{\boldsymbol{u}}}_0$|⁠, representing the difference between the displacements of particle |$i\ $|and particle 0.

$Diagram for the D2Q8 lattice spring model. (a) An undistorted square lattice. $k_{0i}$ represents the linear spring and $b_{i0j}$ represents the angular spring; ${\boldsymbol{k}}$ and ${\boldsymbol{b}}$ represent the coefficients of the linear springs and the angular springs, respectively. (b) Diagram for stretching of the linear spring and bending of the angular spring. (c) A distorted lattice from its balanced position. The dotted model represents the original position of the undistorted lattice and the dashed model represents the current position of the distorted lattice.$

Figure 1.

Diagram for the D2Q8 lattice spring model. (a) An undistorted square lattice. |$k_{0i}$| represents the linear spring and |$b_{i0j}$| represents the angular spring; |${\boldsymbol{k}}$| and |${\boldsymbol{b}}$| represent the coefficients of the linear springs and the angular springs, respectively. (b) Diagram for stretching of the linear spring and bending of the angular spring. (c) A distorted lattice from its balanced position. The dotted model represents the original position of the undistorted lattice and the dashed model represents the current position of the distorted lattice.

Open in new tab Download slide

2.2 Calibrated continuum elastic theory

On the one hand, the elastic properties of the LSM are characterized by spring coefficients. On the other hand, the properties of underground heterogeneous media are typically defined by an elasticity tensor. Therefore, it is essential to establish the relationship between the micromechanical spring coefficients and the macroscopic elasticity tensor. To calibrate the spring coefficients, we compare the body force of a distorted lattice, given by eq. (2), with the 2-D heterogeneous anisotropic elastic equations in the form of

$$\begin{eqnarray} {F}_x = \rho \frac{{{\partial }^2{u}_x}}{{\partial {t}^2}} = \frac{\partial }{{\partial x}}\left[ {{c}_{11}\frac{{\partial {u}_x}}{{\partial x}} + {c}_{13}\frac{{\partial {u}_z}}{{\partial z}} + {c}_{15}\left( {\frac{{\partial {u}_x}}{{\partial z}} + \frac{{\partial {u}_z}}{{\partial x}}} \right)} \right] + \frac{\partial }{{\partial z}}\left[ {{c}_{15}\frac{{\partial {u}_x}}{{\partial x}} + {c}_{35}\frac{{\partial {u}_z}}{{\partial z}} + {c}_{55}\left( {\frac{{\partial {u}_x}}{{\partial z}} + \frac{{\partial {u}_z}}{{\partial x}}} \right)} \right] \end{eqnarray}$$

(3a)

and

$$\begin{eqnarray} {F}_z = \rho \frac{{{\partial }^2{u}_z}}{{\partial {t}^2}} = \frac{\partial }{{\partial x}}\left[ {{c}_{15}\frac{{\partial {u}_x}}{{\partial x}} + {c}_{35}\frac{{\partial {u}_z}}{{\partial z}} + {c}_{55}\left( {\frac{{\partial {u}_x}}{{\partial z}} + \frac{{\partial {u}_z}}{{\partial x}}} \right)} \right] + \frac{\partial }{{\partial z}}\left[ {{c}_{13}\frac{{\partial {u}_x}}{{\partial x}} + {c}_{33}\frac{{\partial {u}_z}}{{\partial z}} + {c}_{35}\left( {\frac{{\partial {u}_x}}{{\partial z}} + \frac{{\partial {u}_z}}{{\partial x}}} \right)} \right], \end{eqnarray}$$

(3b)

where |$\rho $| is the density, |${c}_{ij}$| are the components of the elasticity tensor, and |${u}_x$| and |${u}_z$| are the x- and z-components of the displacement, respectively.

When the medium is homogeneous or varies smoothly, the derivative terms of the medium properties can be neglected; thus, eq. (3) is simplified as

$$\begin{eqnarray} {F}_x = {c}_{11}\frac{{{\partial }^2{u}_x}}{{\partial {x}^2}} + 2{c}_{15}\frac{{{\partial }^2{u}_x}}{{\partial x\partial z}} + {c}_{55}\frac{{{\partial }^2{u}_x}}{{\partial {z}^2}} + {c}_{15}\frac{{{\partial }^2{u}_z}}{{\partial {x}^2}} + \left( {{c}_{13} + {c}_{55}} \right)\frac{{{\partial }^2{u}_z}}{{\partial x\partial z}} + {c}_{35}\frac{{{\partial }^2{u}_z}}{{\partial {z}^2}} \end{eqnarray}$$

(4a)

and

$$\begin{eqnarray} {F}_z = {c}_{33}\frac{{{\partial }^2{u}_z}}{{\partial {z}^2}} + 2{c}_{35}\frac{{{\partial }^2{u}_z}}{{\partial x\partial z}} + {c}_{55}\frac{{{\partial }^2{u}_z}}{{\partial {x}^2}} + {c}_{35}\frac{{{\partial }^2{u}_x}}{{\partial {z}^2}} + \left( {{c}_{13} + {c}_{55}} \right)\frac{{{\partial }^2{u}_x}}{{\partial x\partial z}} + {c}_{15}\frac{{{\partial }^2{u}_x}}{{\partial {x}^2}}. \end{eqnarray}$$

(4b)

In the conventional LSM, spring coefficients are essentially calibrated based on the homogeneous wave equation of eq. (4) (Toomey & Bean 2000; Del Valle‐García & Sánchez‐Sesma 2003; O'Brien & Bean 2004; O'Brien et al. 2009; O′Brien 2014; Xia et al. 2017, 2018; Buxton 2022; Tang et al. 2022, 2023). In this paper, we propose a new LSM, in which spring coefficients are calibrated based on the heterogeneous wave equation of eq. (3). We refer to the new LSM as the medium's derivative LSM (MDLSM) considering that its spring coefficients incorporate the medium's derivative terms while the conventional LSM not. From eqs (3) and (4), the forces used to calibrate the spring coefficients of the conventional LSM and MDLSM are different. To see the difference clearly, we subtract eq. (4) from eq. (3) as

$$\begin{eqnarray} F_x^{\textrm{MDLSM}} - F_x^{\mathrm{ LSM}} = \left( {\frac{{\partial {c}_{11}}}{{\partial x}} + \frac{{\partial {c}_{15}}}{{\partial z}}} \right){\rm{\ }}\frac{{\partial {u}_x}}{{\partial x}} + \left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{35}}}{{\partial z}}} \right){\rm{\ }}\frac{{\partial {u}_z}}{{\partial z}} + \left( {\frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right){\rm{\ }}\left( {\frac{{\partial {u}_x}}{{\partial z}} + \frac{{\partial {u}_z}}{{\partial x}}} \right) \end{eqnarray}$$

(5a)

and

$$\begin{eqnarray} F_z^{\textrm{MDLSM}} - F_z^{\mathrm{ LSM}} = \left( {\frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}}} \right){\rm{\ }}\frac{{\partial {u}_x}}{{\partial x}} + \left( {\frac{{\partial {c}_{35}}}{{\partial x}} + \frac{{\partial {c}_{33}}}{{\partial z}}} \right){\rm{\ }}\frac{{\partial {u}_z}}{{\partial z}} + \left( {\frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{35}}}{{\partial z}}} \right){\rm{\ }}\left( {\frac{{\partial {u}_x}}{{\partial z}} + \frac{{\partial {u}_z}}{{\partial x}}} \right), \end{eqnarray}$$

(5b)

where |${F}^{{\rm{MDLSM}}}$| and |${F}^{{\rm{LSM}}}$| are the forces used to calibrate the spring coefficients of the MDLSM and the conventional LSM, respectively. Eq. (5) indicates that the MDLSM calculates the forces employed by the conventional LSM, as well as the forces with the products of the medium's derivative terms and the displacement's derivative terms. These extra derivative terms play an important role in handling the problem of scattering waves in seismic wave simulation.

2.3 Parameter transformation

Eq. (2) defines the unit force in a distorted lattice. Under the assumption that the displacements of particles are much smaller than their intervals, we have

$$\begin{eqnarray} {\rm{cos}}{\theta }_{i0j} - {\rm{cos}}\theta _{i0j}^0 = - {\rm{sin}}\theta _{i0j}^0\left( {{\theta }_i - {\theta }_j} \right), \end{eqnarray}$$

(6)

$$\begin{eqnarray} {\theta }_i = \frac{{{{\boldsymbol{u}}}_{0i} \cdot {\boldsymbol{s}}_{0i}^0}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}},{\rm{\ }}{\theta }_j = \frac{{{{\boldsymbol{u}}}_{0j} \cdot {\boldsymbol{s}}_{0j}^0}}{{{{\left| {{\boldsymbol{r}}_{0j}^0} \right|}}^2}} \end{eqnarray}$$

(7)

and

$$\begin{eqnarray} |{{\boldsymbol{r}}}_{0i}\left| { - |{\boldsymbol{r}}_{0i}^0} \right| = \frac{{{{\boldsymbol{u}}}_{0i} \cdot {\boldsymbol{r}}_{0i}^0}}{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}, \end{eqnarray}$$

(8)

where |${\theta }_i$| and |${\theta }_j$| are the angles that linear spring 0-i and spring 0-j rotate from their equilibrium positions shown in Fig. 2, respectively. The vector |${\boldsymbol{s}}_{0i}^0$| is perpendicular to the vector |${\boldsymbol{r}}_{0i}^0$| in the counter-clockwise direction and |$|{\boldsymbol{s}}_{0i}^0| = $|||${\boldsymbol{r}}_{0i}^0|$|⁠.

$Geometric relations between the distorted and undistorted lattices. Particles $0,{\boldsymbol{i}},{\boldsymbol{j}}$ are at the equilibrium position in the undistorted lattice; particles 0, i′, j′ are the current position of particles $0,{\boldsymbol{\ i}},{\boldsymbol{\ j}}$ in the distorted lattice;$

Figure 2.

Geometric relations between the distorted and undistorted lattices. Particles |$0,{\boldsymbol{i}},{\boldsymbol{j}}$| are at the equilibrium position in the undistorted lattice; particles 0, i′, j′ are the current position of particles |$0,{\boldsymbol{\ i}},{\boldsymbol{\ j}}$| in the distorted lattice;

Open in new tab Download slide

Substituting eqs (6)–(8) into eq. (2), we have the two components of |${{\boldsymbol{F}}}_0$| as

$$\begin{eqnarray} {F}_x = \frac{1}{{\Delta x\Delta z}}\mathop \sum \limits_{i = 1}^8 {k}_{0i}\frac{{u_{0i}^xr_{0i}^{0x} + u_{0i}^zr_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}r_{0i}^{0x} + \frac{1}{{\Delta x\Delta z}}{\sin }^2\theta _{i0j}^0\mathop \sum \limits_{i = 1,{\rm{\ }}j = \,\mathrm{ mod}\left( {i,8} \right) + 1}^8 {b}_{i0j}\left( {\frac{{u_{0i}^xs_{0i}^{0x} + u_{0i}^zs_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}} - \frac{{u_{0j}^xs_{0j}^{0x} + u_{0j}^zs_{0j}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0j}^0} \right|}}^2}}} \right)\left( {\frac{{s_{0i}^{0x}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}} - \frac{{s_{0j}^{0x}}}{{{{\left| {{\boldsymbol{r}}_{0j}^0} \right|}}^2}}} \right) \end{eqnarray}$$

(9a)

and

$$\begin{eqnarray} {F}_z = \frac{1}{{\Delta x\Delta z}}\mathop \sum \limits_{i = 1}^8 {k}_{0i}\frac{{u_{0i}^xr_{0i}^{0x} + u_{0i}^zr_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}r_{0i}^{0z} + \frac{1}{{\Delta x\Delta z}}{\sin }^2\theta _{i0j}^0\mathop \sum \limits_{i = 1,\ j = \mathrm{ mod}\left( {i,8} \right) + 1}^8 {b}_{i0j}\left( {\frac{{u_{0i}^xs_{0i}^{0x} + u_{0i}^zs_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}} - \frac{{u_{0j}^xs_{0j}^{0x} + u_{0j}^zs_{0j}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0j}^0} \right|}}^2}}} \right)\left( {\frac{{s_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}} - \frac{{s_{0j}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0j}^0} \right|}}^2}}} \right). \end{eqnarray}$$

(9b)

|${\boldsymbol{r}}_{0i}^0$| is defined as

$$\begin{eqnarray} {\boldsymbol{r}}_{0i}^0 = r_{0i}^{0x}{{\boldsymbol{e}}}_x + r_{0i}^{0z}{{\boldsymbol{e}}}_z \end{eqnarray}$$

(10a)

and

$$\begin{eqnarray} {\boldsymbol{s}}_{0i}^0 = s_{0i}^{0x}{{\boldsymbol{e}}}_x + s_{0i}^{0z}{{\boldsymbol{e}}}_z, \end{eqnarray}$$

(10b)

where |${{\boldsymbol{e}}}_x$| and |${{\boldsymbol{e}}}_z$| are the unit vectors in the x- and z–directions, respectively, and |$r_{0i}^{0x}$| and |$r_{0i}^{0z}$| can be given by a matrix form of

$$\begin{eqnarray} \left[ {\begin{array}{@{}*{2}{c}@{}} {r_{01}^{0x}}&{r_{01}^{0z}}\\ {\begin{array}{@{}*{1}{c}@{}} {r_{02}^{0x}}\\ {r_{03}^{0x}}\\ {\begin{array}{@{}*{1}{c}@{}} {r_{04}^{0x}}\\ {r_{05}^{0x}}\\ {\begin{array}{@{}*{1}{c}@{}} {r_{06}^{0x}}\\ {\begin{array}{@{}*{1}{c}@{}} {r_{07}^{0x}}\\ {r_{08}^{0x}} \end{array}} \end{array}} \end{array}} \end{array}}&{\begin{array}{@{}*{1}{c}@{}} {r_{02}^{0z}}\\ {r_{03}^{0z}}\\ {\begin{array}{@{}*{1}{c}@{}} {r_{04}^{0z}}\\ {r_{05}^{0z}}\\ {\begin{array}{@{}*{1}{c}@{}} {r_{06}^{0z}}\\ {\begin{array}{@{}*{1}{c}@{}} {r_{07}^{0z}}\\ {r_{08}^{0z}} \end{array}} \end{array}} \end{array}} \end{array}} \end{array}} \right] = \left[ {\begin{array}{@{}*{2}{c}@{}} h&0\\ {\begin{array}{@{}*{1}{c}@{}} h\\ 0\\ {\begin{array}{@{}*{1}{c}@{}} { - h}\\ { - h}\\ {\begin{array}{@{}*{1}{c}@{}} { - h}\\ {\begin{array}{@{}*{1}{c}@{}} 0\\ h \end{array}} \end{array}} \end{array}} \end{array}}&{\begin{array}{@{}*{1}{c}@{}} h\\ {\begin{array}{@{}*{1}{c}@{}} h\\ h\\ {\begin{array}{@{}*{1}{c}@{}} 0\\ { - h} \end{array}} \end{array}}\\ {\begin{array}{@{}*{1}{c}@{}} { - h}\\ { - h} \end{array}} \end{array}} \end{array}} \right], \end{eqnarray}$$

(11)

in which we set |${\rm{\Delta }}x = {\rm{\Delta }}z = h$| for the uniform lattice pattern. |$s_{0i}^{0x}$| and |$s_{0i}^{0z}$| are presented by

$$\begin{eqnarray} \left[ {\begin{array}{@{}*{1}{c}@{}} {s_{0i}^{0x}}\\ {s_{0i}^{0z}} \end{array}} \right] = \left( {\begin{array}{@{}*{2}{c}@{}} 0&{ - 1}\\ 1&0 \end{array}} \right)\left[ {\begin{array}{@{}*{1}{c}@{}} {r_{0i}^{0x}}\\ {r_{0i}^{0z}} \end{array}} \right]. \end{eqnarray}$$

(12)

Then, eq. (9) can be simplified as

$$\begin{eqnarray} {F}_x = \frac{1}{{{h}^2}}\left\{ \mathop \sum \limits_{\begin{array}{@{}*{1}{c}@{}} {i = 1\ j = \mathrm{ mod}\left( {i,8} \right) + 1,}\\ {k = \mathrm{ mod}\left( {i,8} \right) - 1} \end{array}}^8 \left[\frac{{r_{0i}^{0x}r_{0i}^{0x}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {\sin }^2\theta _{i0j}^0(g_{ij}^{xx}{b}_{i0j} + g_{ik}^{xx}{b}_{i0k})\right] {u_{0i}^x + } \left[\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {\sin }^2\theta _{i0j}^0\left( {g_{ij}^{xz}{b}_{i0j} + g_{ik}^{xz}{b}_{i0k}} \right)\right]u_{0i}^z\right\} \end{eqnarray}$$

(13a)

and

$$\begin{eqnarray} {F}_z = \frac{1}{{{h}^2}}\Bigg\{ \mathop \sum \limits_{\begin{array}{@{}*{1}{c}@{}} {i = 1\ j = \mathrm{ mod}\left( {i,8} \right) + 1,}\\ {k = \mathrm{ mod}\left( {i,8} \right) - 1} \end{array}}^8 \left[\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {\sin }^2\theta _{i0j}^0(g_{ij}^{zx}{b}_{i0j} + g_{ik}^{zx}{b}_{i0k})\right] {u_{0i}^x + } \left[\frac{{r_{0i}^{0z}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {\sin }^2\theta _{i0j}^0\left( {g_{ij}^{zz}{b}_{i0j} + g_{ik}^{zz}{b}_{i0k}} \right)\right]u_{0i}^z, \end{eqnarray}$$

(13b)

in which

$$\begin{eqnarray} g_{ij}^{mn} = \frac{{s_{0i}^{0m}p_{ij}^n}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}\,\,\,\,\,\,g_{ik}^{mn} = \frac{{s_{0i}^{0m}p_{ik}^n}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}},{\rm{\ }}m = x,z,{\rm{\ \ \ }}n = x,z, \end{eqnarray}$$

(14)

$$\begin{eqnarray} p_{ij}^m = \frac{{s_{0i}^{0m}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}} - \frac{{s_{0j}^{0m}}}{{{{\left| {{\boldsymbol{r}}_{0j}^0} \right|}}^2}},{\rm{\ }}m = x,z. \end{eqnarray}$$

(15)

The summations in eq. (13) are implemented over all the neighbours with |$i = 1,\ 2,\ 3,...,\ 8$|⁠. Specially, when |$i = 1$|⁠, we set |$k = 8$| instead of |$k = {\rm{mod}}( {i,8} ) - 1 = 0$|⁠.

To calibrate the spring coefficients of the MDLSM, we expand |$u_{0i}^p$| (⁠|$p = x,z$|⁠) at particle 0 to a second-order accuracy as

$$\begin{eqnarray} u_{0i}^p \approx \frac{{\partial {u}_p}}{{\partial x}}r_{0i}^{0x} + \frac{{\partial {u}_p}}{{\partial z}}r_{0i}^{0z} + \frac{1}{2}\frac{{{\partial }^2{u}_p}}{{\partial {x}^2}}{(r_{0i}^{0x})}^2 + \frac{1}{2}\frac{{{\partial }^2{u}_p}}{{\partial {z}^2}}{(r_{0i}^{0z})}^2 + \frac{{{\partial }^2{u}_p}}{{\partial x\partial z}}r_{0i}^{0x}r_{0i}^{0z},{\rm{\ }}p = x,z, \end{eqnarray}$$

(16)

and substitute eq. (16) into eq. (13). Eqs (3) and (13) are the force equations expressed by the medium parameters and the spring coefficients, respectively. By combining the two equations, the relationship between the spring coefficients of the MDLSM and the elasticity tensor can be obtained in Appendix A. Similarly, the relationship between the spring coefficients of the conventional LSM and the elasticity tensor is obtained by combining eqs (4) and (13). Note that due to disparities in the calibrated equations of eqs (3) and (4), the spring coefficients of the conventional LSM and the MDLSM are different. We have the spring coefficients of the conventional LSM as

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} \displaystyle {k_{01}^{\mathrm{ LSM}} = k_{05}^{\mathrm{ LSM}} = {c}_{11} - \frac{{{c}_{13}}}{2} - \frac{{{c}_{55}}}{2}}\\ \displaystyle {k_{02}^{\mathrm{ LSM}} = k_{06}^{\mathrm{ LSM}} = \frac{{{c}_{13}}}{2} + {c}_{15} + \frac{{{c}_{55}}}{2}}\\ {\begin{array}{@{}*{1}{c}@{}} \displaystyle {k_{03}^{\mathrm{ LSM}} = k_{07}^{\mathrm{ LSM}} = {c}_{33} - \frac{{{c}_{13}}}{2} - \frac{{{c}_{55}}}{2}}\\ \displaystyle {k_{04}^{\mathrm{ LSM}} = k_{08}^{\mathrm{ LSM}} = \frac{{{c}_{13}}}{2} - {c}_{15} + \frac{{{c}_{55}}}{2}} \end{array}} \end{array}} \right. \end{eqnarray}$$

(17a)

and

$$\begin{eqnarray} b_{102}^{\mathrm{ LSM}} = b_{203}^{\mathrm{ LSM}} = b_{304}^{\mathrm{ LSM}} = b_{405}^{\mathrm{ LSM}} = b_{506}^{\mathrm{ LSM}} = b_{607}^{\mathrm{ LSM}} = b_{708}^{\mathrm{ LSM}} = b_{801}^{\mathrm{ LSM}} = - {h}^2\left( {{c}_{13} - {c}_{55}} \right) \end{eqnarray}$$

(17b)

The spring coefficients of the MDLSM are obtained as

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{l}@{}} \displaystyle {k_{01}^{\textrm{MDLSM}} = k_{01}^{\mathrm{ LSM}} + \frac{h}{2}\left( {\frac{{\partial {c}_{11}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial x}} - 2\frac{{\partial {c}_{55}}}{{\partial x}}} \right)}\\ \displaystyle {k_{02}^{\textrm{MDLSM}} = k_{02}^{\mathrm{ LSM}} + \frac{h}{2}\left( {\frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{15}}}{{\partial z}} + \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {k_{03}^{\textrm{MDLSM}} = k_{03}^{\mathrm{ LSM}} + \frac{h}{2}\left( {\frac{{\partial {c}_{13}}}{{\partial z}} + \frac{{\partial {c}_{33}}}{{\partial z}} - 2\frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {k_{04}^{\textrm{MDLSM}} = k_{04}^{\mathrm{ LSM}} + \frac{h}{2}\left( {\frac{{\partial {c}_{15}}}{{\partial x}} - \frac{{\partial {c}_{15}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {k_{05}^{\textrm{MDLSM}} = k_{05}^{\mathrm{ LSM}} - \frac{h}{2}\left( {\frac{{\partial {c}_{11}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial x}} - 2\frac{{\partial {c}_{55}}}{{\partial x}}} \right)}\\ \displaystyle {k_{06}^{\textrm{MDLSM}} = k_{06}^{\mathrm{ LSM}} - \frac{h}{2}\left( {\frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{15}}}{{\partial z}} + \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {k_{07}^{\textrm{MDLSM}} = k_{07}^{\mathrm{ LSM}} - \frac{h}{2}\left( {\frac{{\partial {c}_{13}}}{{\partial z}} + \frac{{\partial {c}_{33}}}{{\partial z}} - 2\frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {k_{08}^{\textrm{MDLSM}} = k_{08}^{\mathrm{ LSM}} - \frac{h}{2}\left( {\frac{{\partial {c}_{15}}}{{\partial x}} - \frac{{\partial {c}_{15}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \end{array}} \right. \end{eqnarray}$$

(18a)

and

$$\begin{eqnarray} \left\{{\begin{array}{@{}*{1}{l}@{}} \displaystyle {b_{102}^{\textrm{MDLSM}} = b_{102}^{\mathrm{ LSM}} - {h}^3\left( {\frac{{\partial {c}_{13}}}{{\partial x}} - \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle{b_{203}^{\textrm{MDLSM}} = b_{203}^{\mathrm{ LSM}} + {h}^3\left( {\frac{{\partial {c}_{13}}}{{\partial x}} - \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {b_{304}^{\textrm{MDLSM}} = b_{304}^{\mathrm{ LSM}} - {h}^3\left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} - \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ {b_{405}^{\textrm{MDLSM}} = b_{405}^{\mathrm{ LSM}} + {h}^3 \displaystyle \left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} - \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {b_{506}^{\textrm{MDLSM}} = b_{506}^{\mathrm{ LSM}} + {h}^3\left( {\frac{{\partial {c}_{13}}}{{\partial x}} - \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {b_{607}^{\textrm{MDLSM}} = b_{607}^{\mathrm{ LSM}} - {h}^3\left( {\frac{{\partial {c}_{13}}}{{\partial x}} - \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {b_{708}^{\textrm{MDLSM}} = b_{708}^{\mathrm{ LSM}} + {h}^3\left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} - \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {b_{801}^{\textrm{MDLSM}} = b_{801}^{\mathrm{ LSM}} - {h}^3\left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} - \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \end{array}.} \right. \end{eqnarray}$$

(18b)

The complete deduction of the spring coefficients can be referred to in Appendix A. From eqs (17) and (18), when the medium is homogeneous or varies smoothly, the MDLSM will convert into the conventional LSM; otherwise, in heterogeneous media, the conventional LSM will result in errors caused by the lack of the derivative terms of the medium in the expressions of spring coefficients.

2.4 Verlet algorithm

The Verlet algorithm is devised for molecular dynamics (Verlet 1967), which has a second-order accuracy in time. To apply the Verlet algorithm, the displacement |${{\boldsymbol{u}}}_i$|⁠, the velocity |${{\boldsymbol{v}}}_i$| and the acceleration |${{\boldsymbol{a}}}_i$| of the particle i are required. The iterative equations are given by

$$\begin{eqnarray} {{\boldsymbol{u}}}_i\left( {t + \Delta t} \right) = {{\boldsymbol{u}}}_i\left( t \right) + {{\boldsymbol{v}}}_i\left( t \right)\Delta t + \frac{1}{2}{{\boldsymbol{a}}}_i\left( t \right)\Delta {t}^2, \end{eqnarray}$$

(19a)

$$\begin{eqnarray} {{\boldsymbol{v}}}_i\left( {t + \frac{1}{2}\Delta t} \right) = {{\boldsymbol{v}}}_i\left( t \right) + \frac{1}{2}{{\boldsymbol{a}}}_i\left( t \right)\Delta t, \end{eqnarray}$$

(19b)

$$\begin{eqnarray} {{\boldsymbol{a}}}_i\left( {t + \Delta t} \right) = \frac{{{{\boldsymbol{F}}}_i\left( {t + \Delta t} \right)}}{{{\rho }_i}} \end{eqnarray}$$

(19c)

and

$$\begin{eqnarray} {{\boldsymbol{v}}}_i\left( {t + \Delta t} \right) = {{\boldsymbol{v}}}_i\left( {t + \frac{1}{2}\Delta t} \right) + \frac{1}{2}{{\boldsymbol{a}}}_i\left( {t + \Delta t} \right)\Delta t, \end{eqnarray}$$

(19d)

where |${{\boldsymbol{F}}}_i$| is the total body force on particle |$\ i$| calculated using eq. (13), |${\rm{\Delta }}t$| is the time interval.

3 DISPERSION ANALYSIS

Generally, the dispersion analysis is derived in homogeneous media, and thus the dispersion relations of the MDLSM and the conventional LSM are the same. We assume the plane wave solutions at |$xoz$| plane as

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} {u = {u}_x{{\boldsymbol{e}}}_x + {u}_z{{\boldsymbol{e}}}_z}\\ {{u}_x = u_x^0\mathrm{ exp}\left( { - \beth \omega t + \beth {k}_xx + \beth {k}_zz} \right)}\\ {{u}_z = u_z^0\mathrm{ exp}\left( { - \beth \omega t + \beth {k}_xx + \beth {k}_zz} \right)} \end{array}} \right., \end{eqnarray}$$

(20)

where |$u_x^0$| and |$u_z^0$| are the amplitude of the displacement in x- and z-directions, respectively; |${k}_x = k\,{\rm{cos}}\,\theta $| and |${k}_z = \ k\,{\rm{sin}}\,\theta $| are the wavenumbers in x- and z-directions, respectively. |${\boldsymbol{k}}$| is the wavenumber and makes a wave-propagation angle |$\theta $| with x-axis; |$\omega $| is the angular frequency; |$\beth = \sqrt { - 1} $| is the imaginary unit. For particle 0, according to eq. (20), we have

$$\begin{eqnarray} {{\boldsymbol{u}}}_{0i} = {{\boldsymbol{u}}}_i - {{\boldsymbol{u}}}_0 = {\boldsymbol{u}}\left[ {\mathrm{ exp}\left( {\beth {k}_xr_{0i}^{0x} + \beth {k}_zr_{0i}^{0z}} \right) - 1} \right]. \end{eqnarray}$$

(21)

Substituting eq. (21) into eq. (13), we obtain

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} {{F}_x = {m}_{11}{u}_x + {m}_{12}{u}_z}\\ {{F}_z = {m}_{21}{u}_x + {m}_{22}{u}_z} \end{array}} \right., \end{eqnarray}$$

(22)

where

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{l}@{}} \displaystyle {{m}_{11} = \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0x}}}{{{{\left| {r_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xx}{b}_{i0j} + g_{ik}^{xx}{b}_{i0k}} \right)} \right]\left[ {\mathrm{ exp}\left( {\beth {k}_xr_{0i}^{0x} + \beth {k}_zr_{0i}^{0z}} \right) - 1} \right]}\\ \displaystyle {{m}_{12} = \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xz}{b}_{i0j} + g_{ik}^{xz}{b}_{i0k}} \right)} \right]\left[ {\mathrm{ exp}\left( {\beth {k}_xr_{0i}^{0x} + \beth {k}_zr_{0i}^{0z}} \right) - 1} \right]}\\ \displaystyle {{m}_{21} = \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zx}{b}_{i0j} + g_{ik}^{zx}{b}_{i0k}} \right)} \right]\left[ {\mathrm{ exp}\left( {\beth {k}_xr_{0i}^{0x} + \beth {k}_zr_{0i}^{0z}} \right) - 1} \right]}\\ \displaystyle {{m}_{22} = \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0z}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zz}{b}_{i0j} + g_{ik}^{zz}{b}_{i0k}} \right)} \right]\left[ {\mathrm{ exp}\left( {\beth {k}_xr_{0i}^{0x} + \beth {k}_zr_{0i}^{0z}} \right) - 1} \right]} \end{array}.} \right. \end{eqnarray}$$

(23)

On the other hand, based on the Verlet algorithm, we have

$$\begin{eqnarray} a\left( t \right) = \frac{{{\partial }^2u\left( t \right)}}{{\partial {t}^2}} = \frac{{u\left( {t + \Delta t} \right) - 2u\left( t \right) + u\left( {t - \Delta t} \right)}}{{\Delta {t}^2}}, \end{eqnarray}$$

(24)

and thus,

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} \displaystyle {{F}_x = \rho \frac{{{\partial }^2{u}_x}}{{\partial {t}^2}} = - \frac{{4\rho }}{{\Delta {t}^2}}{{\sin }}^2\left( {\frac{{\omega \Delta t}}{2}} \right){u}_x}\\ \displaystyle {{F}_z = \rho \frac{{{\partial }^2{u}_z}}{{\partial {t}^2}} = - \frac{{4\rho }}{{\Delta {t}^2}}{{\sin }}^2\left( {\frac{{\omega \Delta t}}{2}} \right){u}_z} \end{array}} \right.. \end{eqnarray}$$

(25)

Combining eqs (22) and (25), we gain

$$\begin{eqnarray} \left[ {\begin{array}{@{}*{2}{c}@{}} {{m}_{11}}&{{m}_{12}}\\ {{m}_{21}}&{{m}_{22}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {{u}_x}\\ {{u}_z} \end{array}} \right] + \frac{{4\rho }}{{\Delta {t}^2}}{\sin }^2\left( {\frac{{\omega \Delta t}}{2}} \right)\left[ {\begin{array}{@{}*{1}{c}@{}} {{u}_x}\\ {{u}_z} \end{array}} \right] = 0. \end{eqnarray}$$

(26)

Hence, the phase velocity |$v( {{k}_x,{k}_z} )$| can be expressed as

$$\begin{eqnarray} v\left( {{k}_x,{k}_z} \right) = \frac{\omega }{k} = \frac{2}{{k\Delta t}}\textrm{arcsin}\Bigg[\frac{{\Delta t}}{{2\sqrt {2\rho } }}\sqrt { - \left( {{m}_{11} + {m}_{22}} \right) \pm \sqrt {{{\left( {{m}_{11} - {m}_{22}} \right)}}^2 + 4{m}_{12}{m}_{21}} } {\rm{\ }}\Bigg], \end{eqnarray}$$

(27)

where ‘|$\pm $|’ represents two modes of waves. For isotropic media, ‘|$+ $|’ and ‘|$- $|’ represent P wave and S wave, respectively.

The number of particles per wavelength is defined by

$$\begin{eqnarray} \mathrm{ PPW} = \frac{\lambda }{h}, \end{eqnarray}$$

(28)

where |$\lambda $| represents the wavelength, typically calculated using the shortest wavelength. To analyse the numerical dispersion errors of the MDLSM, we choose an isotropic model with a Poisson's ratio of 0.25 (i.e. |${v}_p/{v}_s = \sqrt 3 $| ). In this section, PPW is calculated using the shortest S-wave wavelength. Fig. 3 illustrates the phase-velocity numerical errors with PPW = 10. It shows that the relative error of the phase velocity for the P wave is slightly larger than that for the S wave overall. To keep the phase velocity within |$\pm $|2 per cent for the MDLSM, PPW |$\ge 10$| is recommended in numerical experiments. Fig. 4 shows normalized phase-velocity curves versus the PPW for different propagation angles. The propagation angle is set from |$0^\circ $| to |$45^\circ $| due to symmetry. From Fig. 4, as the PPW increases, the P-wave and S-wave phase velocities gradually approach to their respective true velocities, which indicates more accurate results can be obtained by increasing the PPW.

$The spatial variation of the relative errors for the phase velocities with PPW = 10. $\frac{{{{\boldsymbol{v}}}_{\boldsymbol{p}}{\boldsymbol{\mathrm{ d}t}}}}{{\boldsymbol{h}}} = 0.5$. The distance from the origin to the curve denotes the value of relative error.$

Figure 3.

The spatial variation of the relative errors for the phase velocities with PPW = 10. |$\frac{{{{\boldsymbol{v}}}_{\boldsymbol{p}}{\boldsymbol{\mathrm{ d}t}}}}{{\boldsymbol{h}}} = 0.5$|⁠. The distance from the origin to the curve denotes the value of relative error.

Open in new tab Download slide

Figure 4.

The normalized phase-velocity curves versus the PPW for different propagation angles. The top and bottom panels show the curves of P and S waves, respectively. Curves of different colours represent results for different propagation angles.

Open in new tab Download slide

4 STABILITY CONDITION

We derive the stability condition of the MDLSM by Von Neumann's method (Press et al. 1986). Assuming that the coefficients of eq. (3) vary slowly in space and time, we have the numerical error in the form of

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} {{\delta }_x\left( t \right) = \delta _x^0\left( t \right)\mathrm{ exp}\left( {\beth {k}_{x}x + \beth {k}_{z}z} \right)}\\ {{\delta }_z\left( t \right) = \delta _z^0\left( t \right)\mathrm{ exp}\left( {\beth {k}_{x}x + \beth {k}_{z}z} \right)} \end{array}} \right., \end{eqnarray}$$

(29)

where |${\delta }_x( t )$| and |${\delta }_z( t )$| are the numerical errors for the horizontal and vertical displacements, respectively. A heuristic procedure to test the stability is assuming that the error functions’ evolution satisfies the same law as the original wave equation. Substituting eq. (29) into eq. (13) and combining with the Verlet algorithm, we have

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} \displaystyle {\frac{{\delta _x^0{\rm{\ }}\left( {t + \Delta t} \right) - 2\delta _x^0\left( t \right) + \delta _x^0\left( {t - \Delta t} \right)}}{{\Delta {t}^2}} = \frac{1}{\rho }\left[ {{m}_{11}\delta _x^0\left( t \right) + {m}_{12}\delta _z^0\left( t \right)} \right]}\\ \displaystyle {\frac{{\delta _z^0\left( {t + \Delta t} \right) - 2\delta _z^0\left( t \right) + \delta _z^0\left( {t - \Delta t} \right)}}{{\Delta {t}^2}} = \frac{1}{\rho }\left[ {{m}_{21}\delta _x^0\left( t \right) + {m}_{22}\delta _z^0\left( t \right)} \right]} \end{array}} \right.. \end{eqnarray}$$

(30)

Here, the growth rates of the errors can be defined as constants:

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} {{\gamma }_{xx} = \displaystyle \frac{{\delta _x^0\left( {t + \Delta t} \right)}}{{\delta _x^0\left( t \right)}} = \frac{{\delta _x^0\left( t \right)}}{{\delta _x^0\left( {t - \Delta t} \right)}}}\\ \displaystyle {{\gamma }_{zz} = \frac{{\delta _z^0\left( {t + \Delta t} \right)}}{{\delta _z^0\left( t \right)}} = \frac{{\delta _z^0\left( t \right)}}{{\delta _z^0\left( {t - \Delta t} \right)}}} \end{array}} \right. \end{eqnarray}$$

(31)

and

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} \displaystyle {{\gamma }_{xz} = \frac{{\delta _x^0\left( t \right)}}{{\delta _z^0\left( t \right)}}}\\ \displaystyle {{\gamma }_{zx} = \frac{{\delta _z^0\left( t \right)}}{{\delta _x^0\left( t \right)}}} \end{array}} \right.. \end{eqnarray}$$

(32)

Jointly using eqs (30), (31) and (32), we have

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} {\gamma _{xx}^2 - 2{\alpha }_x{\gamma }_{xx} + 1 = 0}\\ {\gamma _{zz}^2 - 2{\alpha }_z{\gamma }_{zz} + 1 = 0} \end{array}} \right., \end{eqnarray}$$

(33)

where

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} \displaystyle {{\alpha }_x = 1 + \frac{{\Delta {t}^2}}{{2\rho }}\left( {{m}_{11} + {\gamma }_{zx}{m}_{12}} \right)}\\ \displaystyle {{\alpha }_z = 1 + \frac{{\Delta {t}^2}}{{2\rho }}\left( {{\gamma }_{xz}{m}_{21} + {m}_{22}} \right)} \end{array}} \right.. \end{eqnarray}$$

(34)

The stability condition requires |$| {{\gamma }_{xx}} |,\ | {{\gamma }_{zz}} | \le 1$|⁠, which leads to |$| {{\alpha }_x} |,\ | {{\alpha }_z} | \le 1$|⁠. In practice, without loss of generality, we assume |${\gamma }_{xx} = {\gamma }_{zz}$|⁠. Therefore, based on eq. (34), the stability condition is obtained as

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{c}@{}} { - 8\rho \le \Delta {t}^2\bigg[({m}_{11} + {m}_{12}) + \sqrt {{{\left( {{m}_{11} - {m}_{22}} \right)}}^2 + 4{m}_{12}{m}_{21}} \bigg] \le 0}\\ { - 8\rho \le \Delta {t}^2\bigg[({m}_{21} + {m}_{22}) - \sqrt {{{\left( {{m}_{11} - {m}_{22}} \right)}}^2 + 4{m}_{12}{m}_{21}} \bigg] \le 0} \end{array}} \right.. \end{eqnarray}$$

(35)

In the case of isotropic media, the stability condition becomes simpler (O'Brien & Bean 2004; O'Brien et al. 2009) in the form of

$$\begin{eqnarray} \Delta t < h/{v}_p, \end{eqnarray}$$

(36)

where |${v}_p$| stands for the phase velocity of the P wave.

5 NUMERICAL EXAMPLES

In this section, we demonstrate the accuracy of the MDLSM by comparing it with the conventional LSM and a high-accuracy collocated-grid FDM ( Zhang & Chen 2006; Sun et al. 2016) in anisotropic media. This FDM solves the hyperbolic velocity–stress equations and employs a dispersion relation preserving/optimization MacCormack scheme to minimize the dispersion and dissipation. The high-accuracy FDM has been extensively validated using the generalized reflection–transmission coefficient method and SEM. According to Zhang & Chen (2006), 10 grid points per shortest wavelength are sufficient to achieve global accuracy for complex models with topography. Scattering, layered and complex models are used to validate the MDLSM's capability in wave propagation. Subsequently, the MDLSM's ability to predict fracture propagation paths is demonstrated using two models with pre-existing fractures. Table 1 shows a summary of simulation parameters for numerical experiments. From Table 1, we specifically used large PPW in our experiments to minimize numerical dispersion errors. The source is applied on the |${v}_z$| component as a single force source; 10 layers of standard PMLs (Zhang & Shen 2010) are set to eliminate the boundary reflections in seismic wave propagation.

Table 1.

Open in new tab

The summary of simulation parameters. The maximum frequency of the Ricker wavelet is calculated by 2.5|${f}_\mathrm{ c}$|⁠, in which |${f}_\mathrm{ c}$| is the dominant frequency. Numerical parameters in the scattering, layered and SEAM models for the FDM, the MDLSM and the conventional LSM are the same.

Numerical model	Scattering	Layered	SEAM	Fracture
Model size (x-direction)	1995 m	1995 m	4345 m	0.0106 m
Model size (z-direction)	1995 m	1995 m	3720 m	0.0306 m
Spatial interval	5 m	5 m	5 m	2 \|$\ \times $\| 10^–4 m
Maximum PPW for P wave	28.8	28.8	44.1	–
Minimum PPW for P wave	21.6	20.0	11.9	–
Maximum PPW for S wave	16	16	24.9	–
Minimum PPW for S wave	14.4	12.8	8.0	–
Time step	1 ms	1 ms	0.5 ms	6.67 \|$\ \times $\| 10^–7 ms
Dominant frequency	10 Hz	10 Hz	10 Hz	–
Time-shift	0.15 s	0.15 s	0.15 s	–
Number of time steps	1001	501	1001	–

Numerical model	Scattering	Layered	SEAM	Fracture
Model size (x-direction)	1995 m	1995 m	4345 m	0.0106 m
Model size (z-direction)	1995 m	1995 m	3720 m	0.0306 m
Spatial interval	5 m	5 m	5 m	2 \|$\ \times $\| 10^–4 m
Maximum PPW for P wave	28.8	28.8	44.1	–
Minimum PPW for P wave	21.6	20.0	11.9	–
Maximum PPW for S wave	16	16	24.9	–
Minimum PPW for S wave	14.4	12.8	8.0	–
Time step	1 ms	1 ms	0.5 ms	6.67 \|$\ \times $\| 10^–7 ms
Dominant frequency	10 Hz	10 Hz	10 Hz	–
Time-shift	0.15 s	0.15 s	0.15 s	–
Number of time steps	1001	501	1001	–

Table 1.

Open in new tab

Numerical model	Scattering	Layered	SEAM	Fracture
Model size (x-direction)	1995 m	1995 m	4345 m	0.0106 m
Model size (z-direction)	1995 m	1995 m	3720 m	0.0306 m
Spatial interval	5 m	5 m	5 m	2 \|$\ \times $\| 10^–4 m
Maximum PPW for P wave	28.8	28.8	44.1	–
Minimum PPW for P wave	21.6	20.0	11.9	–
Maximum PPW for S wave	16	16	24.9	–
Minimum PPW for S wave	14.4	12.8	8.0	–
Time step	1 ms	1 ms	0.5 ms	6.67 \|$\ \times $\| 10^–7 ms
Dominant frequency	10 Hz	10 Hz	10 Hz	–
Time-shift	0.15 s	0.15 s	0.15 s	–
Number of time steps	1001	501	1001	–

Numerical model	Scattering	Layered	SEAM	Fracture
Model size (x-direction)	1995 m	1995 m	4345 m	0.0106 m
Model size (z-direction)	1995 m	1995 m	3720 m	0.0306 m
Spatial interval	5 m	5 m	5 m	2 \|$\ \times $\| 10^–4 m
Maximum PPW for P wave	28.8	28.8	44.1	–
Minimum PPW for P wave	21.6	20.0	11.9	–
Maximum PPW for S wave	16	16	24.9	–
Minimum PPW for S wave	14.4	12.8	8.0	–
Time step	1 ms	1 ms	0.5 ms	6.67 \|$\ \times $\| 10^–7 ms
Dominant frequency	10 Hz	10 Hz	10 Hz	–
Time-shift	0.15 s	0.15 s	0.15 s	–
Number of time steps	1001	501	1001	–

5.1 Scattering model

The background medium is characterized by |$\rho = 1500\ {\rm{kg\ }}{{\rm{m}}}^{ - 3},$| and Thomsen parameters |${v}_{p0} = 3000{\rm{\ m\ }}{{\rm{s}}}^{ - 1}$|⁠, |${v}_{s0} = 2000\ {\rm{m\ }}{{\rm{s}}}^{ - 1}$| and |$\delta = \varepsilon = 0.2$|⁠. The relationship between Thomsen parameters and the elasticity tensor (Thomsen 1986) is given by

$$\begin{eqnarray} {c}_{11} = \rho \left( {1 + 2\varepsilon } \right)v_{p0}^2, \end{eqnarray}$$

(37a)

$$\begin{eqnarray} {c}_{33} = \rho v_{p0}^2, \end{eqnarray}$$

(37b)

$$\begin{eqnarray} {c}_{13} = \rho v_{p0}^2\sqrt {\left( {1 - \frac{{v_{s0}^2}}{{v_{p0}^2}}} \right)\left( {1 - \frac{{v_{s0}^2}}{{v_{p0}^2}} + 2\delta } \right)} - \rho v_{s0}^2 \end{eqnarray}$$

(37c)

and

$$\begin{eqnarray} {c}_{55} = {\rm{\ }}\rho v_{s0}^2, \end{eqnarray}$$

(37d)

where |${v}_{p0}$|⁠, |${v}_{s0}$|⁠, |$\delta $| and |$\varepsilon $| are the Thomsen parameters.

The scattering body is placed at (1145 m, 850 m) with a size of 5 m |$\ \times $| 5 m, which is characterized by |$\rho = 1400{\rm{\ kg\ }}{{\rm{m}}}^{ - 3}$|⁠, |${v}_{p0} = 2700\ {\rm{m\ }}{{\rm{s}}}^{ - 1}$|⁠, |${v}_{s0} = 1800\ {\rm{m\ }}{{\rm{s}}}^{ - 1}$| and |$\delta = \varepsilon = 0.1$|⁠. Fig. 5 shows the acquisition geometry for the scattering model. The source is located at (745 m, 1250 m) and six receivers are laid for the observation of scattering waves from different directions.

$Diagram of the scattering model. The geometry is shown on the ${{\boldsymbol{v}}}_{{\boldsymbol{p}}0}$ model. The star represents the location of the source and the inverted triangles represent receivers. The dashed box enlarges the region containing the scatter body.$

Figure 5.

Diagram of the scattering model. The geometry is shown on the |${{\boldsymbol{v}}}_{{\boldsymbol{p}}0}$| model. The star represents the location of the source and the inverted triangles represent receivers. The dashed box enlarges the region containing the scatter body.

Open in new tab Download slide

Fig. 6 shows snapshots of |${v}_z$| for the FDM, the conventional LSM, the MDLSM and their differences for the scattering model, respectively. Given the relatively small size of the scattering body, the scattering waves are mingled with the directed waves. Notably, the snapshots obtained through the three methods exhibit high similarity. From Figs 6(d) and (e), the differences between the two LSMs and the FDM are mainly caused by numerical dispersion of the directed S waves. These differences generally remain below 5 per cent and conceal the scattering waves. Fig. 7 shows the seismograms of |${v}_z$| for the six receivers shown in Fig. 5. The seismograms from the conventional LSM, the MDLSM and the FDM exhibit an overall agreement, with discrepancies remaining below 5 per cent. The numerical dispersion of directed waves makes it challenging to discern scattering waves in the observations.

$Snapshots of ${v}_z$ for the scattering model. (a) FDM. (b) Conventional LSM. (c) MDLSM. (d) Differences between the conventional LSM and the FDM. (e) Differences between the MDLSM and the FDM.$

Figure 6.

Snapshots of |${v}_z$| for the scattering model. (a) FDM. (b) Conventional LSM. (c) MDLSM. (d) Differences between the conventional LSM and the FDM. (e) Differences between the MDLSM and the FDM.

Open in new tab Download slide

$Seismograms of ${v}_z$ at the six receivers. In Fig. 7 (a) and (b), the solid lines represent the seismograms calculated by the FDM. The triangles and circles represent the seismograms calculated by the conventional LSM and the MDLSM, respectively. In Fig. 7 (c) and (d), the two solid lines represent the differences of seismograms for the conventional LSM and the MDLSM compared with the FDM, respectively. The differences are magnified by 20 times. Panels (a) and (c) are ${v}_z$ at receivers 1, 2 and 3. Panels (b) and (d) are ${v}_z$ at receivers 4, 5 and 6.$

Figure 7.

Seismograms of |${v}_z$| at the six receivers. In Fig. 7 (a) and (b), the solid lines represent the seismograms calculated by the FDM. The triangles and circles represent the seismograms calculated by the conventional LSM and the MDLSM, respectively. In Fig. 7 (c) and (d), the two solid lines represent the differences of seismograms for the conventional LSM and the MDLSM compared with the FDM, respectively. The differences are magnified by 20 times. Panels (a) and (c) are |${v}_z$| at receivers 1, 2 and 3. Panels (b) and (d) are |${v}_z$| at receivers 4, 5 and 6.

Open in new tab Download slide

To further investigate scattering effects, we subtract the reference wavefields calculated in the background homogeneous model from the total wavefields. Fig. 8 shows snapshots of |${v}_z$| for the scattering waves calculated by the FDM, the conventional LSM, and the MDLSM, respectively. From Fig. 8, the amplitudes of the scattering waves are relatively small, constituting only about 0.1 per cent of the directed waves. It is clearly seen that the scattering waves calculated by the MDLSM and the FDM are similar, whereas those calculated by the conventional LSM and the FDM differ significantly. To better observe the differences among them, Fig. 9 shows the seismograms of |${v}_z$| for the scattering waves at six receivers. Receivers 1, 2 and 3 record the backward scattering waves and receivers 4, 5 and 6 record the forward scattering waves. In Fig. 9(a), the backward scattering waves calculated by the MDLSM agree well with those of the FDM and only have slight differences in amplitudes, likely attributable to the numerical dispersion. In Fig. 9(b), the forward scattering waves between the MDLSM and the FDM show higher consistency in amplitudes compared to the backward scattering waves. On the contrary, the backward scattering waves calculated by the conventional LSM markedly differ from those obtained by the FDM, displaying significant differences in amplitudes and even a reversal in polarities. Although the forward scattering waves with the conventional LSM provide relatively accurate phases, there are notable errors in amplitudes. This illustrates that the extra medium derivative terms in the spring coefficients of the MDLSM play a crucial role in accurately simulating the scattering waves.

$Snapshots of ${v}_z$ for the scattering waves. (a) FDM. (b) Conventional LSM. (c) MDLSM.$

Figure 8.

Snapshots of |${v}_z$| for the scattering waves. (a) FDM. (b) Conventional LSM. (c) MDLSM.

Open in new tab Download slide

$Seismograms of ${v}_z$ for the scattering waves. (a) Backward scattering waves. (b) Forward scattering waves.$

Figure 9.

Seismograms of |${v}_z$| for the scattering waves. (a) Backward scattering waves. (b) Forward scattering waves.

Open in new tab Download slide

The results above indicate that the MDLSM offers high accuracy and can produce nearly identical wavefields to the FDM. Meanwhile, the conventional LSM accurately captures the directed waves but yields unreliable results for scattering waves, especially for the backward scattering waves. Therefore, for cases involving scattering effects, the MDLSM is more preferred than the conventional method.

5.2 Layered model

Table 2 shows the anisotropic properties of the medium. Fig. 10 depicts the acquisition geometry and the velocity model. The source is located at (995 m, 1000 m) and five receivers are laid for observation of the wavefields.

$Diagram of the two-layered model. The geometry is shown on the ${{\boldsymbol{v}}}_{{\boldsymbol{p}}0}$ model.$

Figure 10.

Diagram of the two-layered model. The geometry is shown on the |${{\boldsymbol{v}}}_{{\boldsymbol{p}}0}$| model.

Open in new tab Download slide

Table 2.

Open in new tab

The anisotropic properties of the medium for the two-layered model.

Depth (km)	\|${v}_{p0}$\| (⁠\|${\rm{m\ }}{{\rm{s}}}^{ - 1}$\|⁠)	\|${v}_{s0}$\| (⁠\|${\rm{m\ }}{{\rm{s}}}^{ - 1}$\|⁠)	\|$\delta $\|	\|$\varepsilon $\|	\|$\rho \ ( {{\rm{kg\ }}{{\rm{m}}}^{ - 3}} )$\|
0–1.25	2500	1600	0.1	0.1	1000
1.25–1.995	3000	2000	0.2	0.2	1500

Depth (km)	\|${v}_{p0}$\| (⁠\|${\rm{m\ }}{{\rm{s}}}^{ - 1}$\|⁠)	\|${v}_{s0}$\| (⁠\|${\rm{m\ }}{{\rm{s}}}^{ - 1}$\|⁠)	\|$\delta $\|	\|$\varepsilon $\|	\|$\rho \ ( {{\rm{kg\ }}{{\rm{m}}}^{ - 3}} )$\|
0–1.25	2500	1600	0.1	0.1	1000
1.25–1.995	3000	2000	0.2	0.2	1500

Table 2.

Open in new tab

The anisotropic properties of the medium for the two-layered model.

Depth (km)	\|${v}_{p0}$\| (⁠\|${\rm{m\ }}{{\rm{s}}}^{ - 1}$\|⁠)	\|${v}_{s0}$\| (⁠\|${\rm{m\ }}{{\rm{s}}}^{ - 1}$\|⁠)	\|$\delta $\|	\|$\varepsilon $\|	\|$\rho \ ( {{\rm{kg\ }}{{\rm{m}}}^{ - 3}} )$\|
0–1.25	2500	1600	0.1	0.1	1000
1.25–1.995	3000	2000	0.2	0.2	1500

Depth (km)	\|${v}_{p0}$\| (⁠\|${\rm{m\ }}{{\rm{s}}}^{ - 1}$\|⁠)	\|${v}_{s0}$\| (⁠\|${\rm{m\ }}{{\rm{s}}}^{ - 1}$\|⁠)	\|$\delta $\|	\|$\varepsilon $\|	\|$\rho \ ( {{\rm{kg\ }}{{\rm{m}}}^{ - 3}} )$\|
0–1.25	2500	1600	0.1	0.1	1000
1.25–1.995	3000	2000	0.2	0.2	1500

Fig. 11 shows snapshots of |${v}_z$| for the FDM, the conventional LSM, the MDLSM and their differences for the two-layered model, respectively. In Figs 11(a)–(c), the directed waves, reflected waves and transmitted waves are distinctly visible. Figs 11(a) and (b) reveal a notable difference in the amplitudes of the reflected waves between the conventional LSM and the FDM, while the transmitted waves exhibit identical patterns. Comparatively, Figs 11(a) and (c) demonstrate that wavefields obtained by the MDLSM and the FDM have identical phase positions and similar amplitudes for the whole snapshot. Figs 11(d) and (e) provide a clearer depiction, highlighting significant errors in the conventional LSM. In contrast, the MDLSM exhibits much smaller errors overall. Fig. 12 shows the seismograms of |${v}_z$| for the five receivers shown in Fig. 10. Notably, the differences between the MDLSM and the FDM are consistently smaller than those between the conventional LSM and the FDM. This reveals that the MDLSM enhances the accuracy of the conventional LSM for seismic wave simulation in heterogeneous media.

$Snapshots of ${v}_z$ for the two-layered model. (a) FDM. (b) Conventional LSM. (c) MDLSM. (d) Differences between the conventional LSM and the FDM. (e) Differences between the MDLSM and the FDM.$

Figure 11.

Snapshots of |${v}_z$| for the two-layered model. (a) FDM. (b) Conventional LSM. (c) MDLSM. (d) Differences between the conventional LSM and the FDM. (e) Differences between the MDLSM and the FDM.

Open in new tab Download slide

$Seismograms of ${v}_z$ at the five receivers. The differences are magnified by five times. Panels (a) and (c) are ${v}_z$ at receivers 1, 2 and 3. Panels (b) and (d) are ${v}_z$ at receivers 4 and 5.$

Figure 12.

Seismograms of |${v}_z$| at the five receivers. The differences are magnified by five times. Panels (a) and (c) are |${v}_z$| at receivers 1, 2 and 3. Panels (b) and (d) are |${v}_z$| at receivers 4 and 5.

Open in new tab Download slide

Fig. 13 provides snapshots of the reflected waves and transmitted waves of |${v}_z$| for the FDM, the conventional LSM and the MDLSM, respectively. The reflected waves are obtained by subtracting the reference wavefields calculated in the homogeneous background model from the total wavefields calculated in the two-layered model. From Figs 13(a)–(c), the reflected waves calculated by the MDLSM coincide well with those calculated by the FDM. In comparison, the reflected waves calculated by the conventional LSM have low reliability, as their snapshots significantly deviate from those calculated by the other two methods. From Figs 13(d)–(f), the three methods show good consistency in the phases of the transmitted waves. Fig. 14 shows the seismograms of |${v}_z$| for the reflected waves and transmitted waves. Receivers 1, 2 and 3, positioned in the top layer, record the reflected waves, while receivers 4 and 5 in the bottom layer record the transmitted waves. Fig. 14(a) demonstrates that the reflected waves calculated by the MDLSM and the FDM agree well with each other; in comparison, the reflected waves calculated by the conventional LSM exhibit substantial differences in amplitudes and phases. For the transmitted waves, Fig. 14(b) demonstrates good consistency between the MDLSM and the FDM, displaying little differences in amplitudes and phases. In contrast, despite relatively accurate phases, the conventional LSM has notable errors in amplitudes compared with the FDM. The inaccuracies in wavefields calculated by the conventional LSM may have significant consequences, especially when processing reflected waves.

$Snapshots of ${v}_z$ for the reflected waves and the transmitted waves. (a) Reflected waves by the FDM. (b) Reflected waves by the conventional LSM. (c) Reflected waves by the MDLSM. (d) Transmitted waves by the FDM. (e) Transmitted waves by the conventional LSM. (f) Transmitted waves by the MDLSM.$

Figure 13.

Snapshots of |${v}_z$| for the reflected waves and the transmitted waves. (a) Reflected waves by the FDM. (b) Reflected waves by the conventional LSM. (c) Reflected waves by the MDLSM. (d) Transmitted waves by the FDM. (e) Transmitted waves by the conventional LSM. (f) Transmitted waves by the MDLSM.

Open in new tab Download slide

$Seismograms of ${v}_z$ for the reflected waves and the transmitted waves. (a) Reflected waves. (b) Transmitted waves.$

Figure 14.

Seismograms of |${v}_z$| for the reflected waves and the transmitted waves. (a) Reflected waves. (b) Transmitted waves.

Open in new tab Download slide

5.3 SEAM model

To further validate the accuracy of the MDLSM in modelling seismic wave propagation in complex media, the SEAM model is tested in this experiment. Fig. 15 depicts the acquisition geometry and the velocity structure for the SEAM model. The source is located at (2190 m, 1865 m) and two receivers are laid for the observation of the wavefields.

$Diagram of a vertical 2D slice of the SEAM model. The geometry is shown on the ${v}_{p0}$ model. The model contains a salt body at shallow depth. The sedimentary layers are folded around the salt and segmented by some faults. (a) ${v}_{p0}$. (b) ${v}_{s0}$. (c) $\delta $. (d) $\varepsilon $. (e) $\ \rho $.$

Figure 15.

Diagram of a vertical 2D slice of the SEAM model. The geometry is shown on the |${v}_{p0}$| model. The model contains a salt body at shallow depth. The sedimentary layers are folded around the salt and segmented by some faults. (a) |${v}_{p0}$|⁠. (b) |${v}_{s0}$|⁠. (c) |$\delta $|⁠. (d) |$\varepsilon $|⁠. (e) |$\ \rho $|⁠.

Open in new tab Download slide

To illustrate the differences of the spring coefficients between the conventional LSM and the MDLSM, the horizontal |${k}_1$| spring is shown in Fig. 16. Note that the differences are mainly at the interface, indicating impacts on scattering waves. Fig. 17 shows snapshots of |${v}_z$| for the FDM, the conventional LSM, the MDLSM and their differences for the SEAM model, respectively. From Figs 17(a)–(c), differences in amplitudes and phases between the conventional LSM and the FDM are notable, for example, at locations (2.3 km, 0.8 km) and (2.1 km, 1.6 km). Fig. 17(d) highlights significant differences in wavefields between the conventional LSM and the FDM, indicating substantial errors in the wavefields calculated by the conventional LSM within this complex model. In contrast, Fig. 17(e) shows that the differences in wavefields between the MDLSM and the FDM are small, demonstrating the high accuracy of the MDLSM in complex media. Fig. 18 shows seismograms of |${v}_z$| for the two receivers shown in Fig. 15(a). Notably, the difference between the MDLSM and the FDM is significantly smaller overall than that between the conventional LSM and the FDM.

$The spring coefficient ${k}_1$ for the SEAM model. (a) MDLSM results. (b) Conventional LSM results. (c) Differences between the MDLSM and the conventional LSM results.$

Figure 16.

The spring coefficient |${k}_1$| for the SEAM model. (a) MDLSM results. (b) Conventional LSM results. (c) Differences between the MDLSM and the conventional LSM results.

Open in new tab Download slide

$Snapshots of ${v}_z$ for the SEAM model. (a) FDM. (b) Conventional LSM. (c) MDLSM. (d) Differences between the conventional LSM and the FDM. (e) Differences between the MDLSM and the FDM.$

Figure 17.

Snapshots of |${v}_z$| for the SEAM model. (a) FDM. (b) Conventional LSM. (c) MDLSM. (d) Differences between the conventional LSM and the FDM. (e) Differences between the MDLSM and the FDM.

Open in new tab Download slide

$Seismograms of ${v}_z$ at the two receivers. The differences are magnified by five times. Panels (a) and (b) are ${v}_z$ at receivers 1 and 2.$

Figure 18.

Seismograms of |${v}_z$| at the two receivers. The differences are magnified by five times. Panels (a) and (b) are |${v}_z$| at receivers 1 and 2.

Open in new tab Download slide

5.4 Fracture model

Naturally fractured structures are ubiquitous, and the LSM is suitable for modelling such cases because it can simulate fractures by simply removing springs that exceed the strength. In this section, to clarify the motivation for extending the LSM to fracture propagation, we conduct two simple fracture failure simulations using the proposed MDLSM. The first fracture model was previously studied by Chen et al. (2016). The geometry and the loading setup for the model are shown in Fig. 19. The horizontal fracture is positioned approximately at the centre of the model, which measures 4 mm in length and 0.2 mm in width. The pre-existing fracture is defined by simply removing the corresponding springs. Elastic constants for the model are

$$\begin{eqnarray} C = \begin{bmatrix}281 & 133 & 0 \\133 & 200 & 0 \\ 0 & 0 & 114 \end{bmatrix} {\rm{GPa}}. \end{eqnarray}$$

(38)

$Diagram of the horizontal fracture model. The white area at the centre of the model represents the pre-existing fracture.$

Figure 19.

Diagram of the horizontal fracture model. The white area at the centre of the model represents the pre-existing fracture.

Open in new tab Download slide

The top boundary is fixed in the vertical direction. The left and the right boundaries are implemented as free surfaces by removing any springs outside the model. The bottom boundary is applied with a very small force to achieve the uniaxial tensile loading.

We adopt a critical elongation criterion similar to that used in Chen et al. (2016). The critical elongation |$s_j^{{\rm{critical}}\ }$| in each direction, that is, j = 1, 2, 3, depends on the corresponding length of the spring and material strength. It can be written as

$$\begin{eqnarray} s_j^{\textrm{critical}{\rm{\ }}} = {\alpha }_j \cdot {R}_j, \end{eqnarray}$$

(39)

where |${\alpha }_j$| is calibrated from experimental tests, and |${R}_j$| is the spring length. In the lattice model shown in Fig. 1, tensile failure occurs when

$$\begin{eqnarray} \delta {l}_j > {\rm{\ }}s_j^{\textrm{critical}{\rm{\ }}}, \end{eqnarray}$$

(40)

where |$\delta {l}_j$| is the elongation of the|${\rm{\ }}$|linear spring |$0$|-j. With the assumption of critical elongations, allowing for the material's symmetry, we set

$$\begin{eqnarray} {\alpha }_1 = {\alpha }_5 = 0.024, \end{eqnarray}$$

(41a)

$$\begin{eqnarray} {\alpha }_3 = {\alpha }_7 = 0.016, \end{eqnarray}$$

(41b)

$$\begin{eqnarray} {\alpha }_2 = {\alpha }_4 = {\rm{\ }}{\alpha }_6 = {\alpha }_8 = 0.008. \end{eqnarray}$$

(41c)

In the failure process, when |$\delta {l}_j$| reaches the critical elongation, the linear spring between particle 0 and particle j is removed, and so are the angular springs containing particle j. For instance, when the linear spring 0-2 is removed, the angular springs 1-0-2 and 2-0-3 are also eliminated.

Fig. 20 shows the failure paths for the horizontal fracture model. The linear springs are marked by lines when they are removed. As time goes, new fractures initiate from the ends of the pre-existing fracture and gradually extend horizontally until they reach the boundary. The final fracture path aligns with findings of Chen et al. (2016). This test demonstrates the MDLSM's capability for fracture failure simulation at the rock scale. To demonstrate the MDLSM's capability of simulating complex fracture propagation, we analysed a model with a pre-existing inclined fracture. The geometry and loading setup for this model are shown in Fig. 21, with the same modelling parameters and failure criteria applied as in the horizontal fracture model (Fig. 19). Fig. 22 shows the failure paths for the inclined fracture model. As the simulation progresses, new fractures initiate at the ends of the original fracture and gradually extend horizontally until reaching the boundary. Compared to the horizontal fracture model, new fracture formation in the inclined model occurs later. In addition, new fractures in the horizontal fracture model propagate symmetrically, with nearly identical failure velocities on both the left and right sides, whereas in the inclined fracture model, new fractures propagate faster on the left side than on the right.

$Dynamic failure path of the fracture model. Panels (a)–(f) represent the results at the 100th, 5000th, 5500th, 6000th, 7000th and 9000th time steps, respectively.$

Figure 20.

Dynamic failure path of the fracture model. Panels (a)–(f) represent the results at the 100th, 5000th, 5500th, 6000th, 7000th and 9000th time steps, respectively.

Open in new tab Download slide

$Diagram of the inclined fracture model. The white area at the centre of the model represents the pre-existing fracture.$

Figure 21.

Diagram of the inclined fracture model. The white area at the centre of the model represents the pre-existing fracture.

Open in new tab Download slide

$Dynamic failure path of the inclined fracture model. Panels (a)–(f) represent the results at the 100th, 9500th, 10 000th, 11 000th, 12 000th and 15 000th time steps, respectively.$

Figure 22.

Dynamic failure path of the inclined fracture model. Panels (a)–(f) represent the results at the 100th, 9500th, 10 000th, 11 000th, 12 000th and 15 000th time steps, respectively.

Open in new tab Download slide

It is noted that, in the context of quasi-static or static loading conditions used in rock failure simulations, the free surface caused by fractures can be effectively represented by simply removing the springs. The LSM has already been widely used in fracture simulation for brittle or quasi-brittle materials, such as concrete, rock and wood (Zhao et al. 2011; Chen et al. 2014, 2016; Zhao 2015; Buxton 2022). However, for dynamic seismic wave modelling, the free surface caused by fractures may not be accurately represented by simply removing the springs, especially when the Poisson’s ratio of the medium deviates from 0.25 (O'Brien 2014). O’Brien (2014) provided an empirical correction term for isotropic media to modify the spring coefficients beneath the free surface, which helps reduce the mismatch with the free surface boundary condition, thereby improving the accuracy of the LSM for arbitrary Poisson's ratios. Our future research will focus on how to modify the spring coefficients in complex anisotropic media to more accurately implement the free surface. Additionally, different challenges also arise in seismic wave modelling, such as modelling the infinite underground space within a finite computational domain, which requires the implementation of absorbing boundaries, or dealing with fractures smaller than the particle intervals.

6 CONCLUSIONS

We propose a new MDLSM for seismic wave simulation in heterogeneous anisotropic media. Unlike conventional LSMs that rely on the wave equation with an implicit homogeneous approximation for calibrating spring coefficients, our MDLSM calibrates these coefficients using the modified wave equation for heterogeneous media. We investigate the differences of spring coefficients between the MDLSM and the conventional LSM. The dispersion property and stability conditions are derived for proper parameters selections in simulations. To validate the accuracy of the MDLSM, three anisotropic models are employed, that is, a scattering model, a two-layered model and a complex SEAM model. Results from these experiments reveal that the conventional LSM exhibits relatively large errors in heterogeneous media, producing incorrect amplitudes and even inverted polarities for the backward scattering waves and reflected waves. In contrast, the MDLSM agrees well with the FDM and exhibits effectiveness in handling complex media. Consequently, the MDLSM proves to be a more robust tool for seismic wave simulation in complex heterogeneous media. Additionally, we conduct simulations in two fracture models using the MDLSM, demonstrating its capability to capture complex failure process.

To expand the application of the MDLSM in seismic exploration, future efforts should focus on several aspects. The free surface caused by fractures may not be accurately captured by simply removing the springs, particularly when the Poisson’s ratio differs from 0.25. Future work will explore how to adjust the spring coefficients beneath the free surface in anisotropic media to enhance the accuracy of its representation. Implementing absorbing boundaries and undulating free surfaces will be critical for enhancing simulation accuracy. Irregular MDLSM models are required for accurately capturing scattering and fracture shapes; quantitative validation of fracture models through experimental testing is necessary. Moreover, the complete earthquake process can be simulated by employing the DEM in large deformation areas and the MDLSM in small deformation areas.

ACKNOWLEDGMENTS

This study received support from the National Natural Science Foundation of China (42074125).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

REFERENCES

Alterman

Karal

F.C.

Jr.,

1968

Propagation of elastic waves in layered media by finite difference methods

Bull. seism. Soc. Am.

(

367

–

398

10.1785/BSSA0580010367

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Buxton

G.A.

2022

An irregular lattice spring model: uniform elasticity, grid refinement and isotropic crack propagation

Model. Simul. Mater. Sci. Eng.

(

055 002

doi:10.1088/1361-651X/ac6c43

10.1088/1361-651X/ac6c43

Chen

Jiao

Liu

2016

A nonlocal lattice particle model for fracture simulation of anisotropic materials

Compos. B: Eng.

141

–

151

10.1016/j.compositesb.2015.12.028

Chen

Lin

Jiao

Liu

2014

A generalized 2D non-local lattice spring model for fracture simulation

Comput. Mech.

1541

–

1558

10.1007/s00466-014-1075-4

Cheng

Liu

2009

New discrete element models for elastoplastic problems

Acta Mech. Sin.

629

–

637

10.1007/s10409-009-0271-5

Cundall

P.A.

Strack

O.D.

1979

A discrete numerical model for granular assemblies

Géotechnique

(

–

10.1680/geot.1979.29.1.47

Del Valle-García

Sánchez-Sesma

F.J.

2003

Rayleigh waves modeling using an elastic lattice model

Geophys. Res. Lett.

1866

Google Scholar

OpenURL Placeholder Text

WorldCat

Di Bartolo

Dors

Mansur

W.J.

2012

A new family of finite-difference schemes to solve the heterogeneous acoustic wave equation

Geophysics

(

T187

–

T199

10.1190/geo2011-0345.1

Graves

R.W.

1996

Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences

Bull. seism. Soc. Am.

(

1091

–

1106

10.1785/BSSA0860041091

Jia

2016

A dynamic lattice method for elastic seismic modeling in anisotropic media

Geophysics

(

T131

–

T143

10.1190/geo2015-0511.1

Jia

2018

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

Geophys. J. Int.

212

(

1868

–

1889

Kelly

K.R.

Ward

R.W.

Treitel

Alford

R.M.

1976

Synthetic seismograms: a finite-difference approach

Geophysics

(

–

Komatitsch

Tromp

1999

Introduction to the spectral element method for three-dimensional seismic wave propagation

Geophys. J. Int.

139

(

806

–

822

10.1046/j.1365-246x.1999.00967.x

Komatitsch

Vilotte

J.-P.

1998

The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures

Bull. seism. Soc. Am.

(

368

–

392

10.1785/BSSA0880020368

Duan

2019

Elastic wavefield forward modeling in heterogeneous media based on the quasi-regular grid high-order finite difference

Oil Geophys. Prospect.

(

539

–

550

Google Scholar

OpenURL Placeholder Text

WorldCat

Liu

2006

Application of discrete element method for continuum dynamic problems

Arch. Appl. Mech.

(

3-4

229

–

243

10.1007/s00419-006-0018-8

Liu

2021

Modeling seismic responses in complex fractured media using the modified lattice spring model coupled with discrete fracture networks

J. Nat. Gas Sci. Eng.

104 206

doi:10.1016/j.jngse.2021.104206

10.1016/j.jngse.2021.104206

Lysmer

Drake

L.A.

1972

A Finite Element Method for Seismology. Methods of Computational Physics

Vol. 11

Academic Press

New York

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Lyu

Masson

Romanowicz

Zhao

2024

Introduction to the distributional finite difference method for 3D seismic wave propagation and comparison with the spectral element method

J. Geophys. Res.: Solid Earth

129

(

doi:10.1029/2023JB027576

10.1029/2023JB027576

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Marfurt

K.J.

1984

Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations

Geophysics

(

533

–

549

Masson

2023

Distributional finite-difference modelling of seismic waves

Geophys. J. Int.

233

(

264

–

296

Masson

Lyu

Moczo

Capdeville

Romanowicz

Virieux

2024

2-D seismic wave propagation using the distributional finite-difference method: further developments and potential for global seismology

Geophys. J. Int.

237

(

339

–

363

Masson

Virieux

2023

P-SV-wave propagation in heterogeneous media: velocity-stress distributional finite-difference method

Geophysics

(

T165

–

T183

10.1190/geo2022-0118.1

Monette

Anderson

M.P.

1994

Elastic and fracture properties of the two-dimensional triangular and square lattices

Model. Simul. Mater. Sci. Eng.

(

doi:10.1088/0965-0393/2/1/004

10.1088/0965-0393/2/1/004

Mora

Place

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

O'Brien

G.S.

2014

Elastic lattice modelling of seismic waves including a free surface

Comput. Geosci.

(

117

–

124

10.1016/j.cageo.2014.03.011

O'Brien

G.S.

2021

A lattice method for seismic wave propagation in nonlinear viscoelastic media

Geophys. J. Int.

224

(

1572

–

1587

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.

(

L14608

doi:10.1029/2004GL020069

O'Brien

G.S.

Bean

C.J.

Tapamo

2009

Dispersion analysis and computational efficiency of elastic lattice methods for seismic wave propagation

Comput. Geosci.

(

1768

–

1775

10.1016/j.cageo.2008.12.004

Ostoja-Starzewski

2002

Lattice models in micromechanics

Appl. Mech. Rev.

(

–

Press

W.H.

Teukolsky

S.A.

Flannery

B.P.

Vetterling

W.T.

1986

Numerical Recipes in Fortran 77: The Art of Scientific Computing

, 2nd edn,

Cambridge Univ. Press

Cambridge

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Zhang

Eisenträger

Song

2022

A time-domain approach for the simulation of three-dimensional seismic wave propagation using the scaled boundary finite element method

Soil Dyn. Earthq. Eng.

152

107 011

doi:10.1016/j.soildyn.2021.107011

10.1016/j.soildyn.2021.107011

Saltzer

S.D.

Pollard

D.D.

1992

Distinct element modeling of structures formed in sedimentary overburden by extensional reactivation of basement normal faults

Tectonics

(

165

–

174

Sun

Zhang

Chen

2016

Seismic-wave modeling in the presence of surface topography in 2D general anisotropic media by a curvilinear grid finite-difference method

Bull. seism. Soc. Am.

106

(

1036

–

1054

Tang

Xia

Zhou

Jiang

Zheng

2023

Lattice spring model for irregular interface based on an adaptive location strategy

IEEE Trans. Geosci. Remote Sens.

–

Google Scholar

OpenURL Placeholder Text

WorldCat

Tang

Zhou

Jiang

Xia

Chen

Zheng

2022

A perfectly matched layer technique applied to lattice spring model in seismic wavefield forward modeling for Poisson's solids

Bull. seism. Soc. Am.

112

(

608

–

621

Thomsen

1986

Weak elastic anisotropy

Geophysics

(

1954

–

1966

Toomey

Bean

C.J.

2000

Numerical simulation of seismic waves using a discrete particle scheme

Geophys. J. Int.

141

(

595

–

604

10.1046/j.1365-246x.2000.00094.x

Tromp

Komatitsch

Liu

2008

Spectral-element and adjoint methods in seismology

Commun. Comput. Phys.

(

–

Google Scholar

OpenURL Placeholder Text

WorldCat

Verlet

1967

Computer ‘experiments’ on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules

Phys. Rev.

159

(

doi:10.1103/PhysRev.159.98

10.1103/PhysRev.159.98

Virieux

1984

SH-wave propagation in heterogeneous media: velocity-stress finite-difference method

Geophysics

(

1933

–

1942

Virieux

1986

P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method

Geophysics

(

889

–

901

Virieux

Calandra

Plessix

2011

A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging

Geophys. Prospect.

794

–

813

10.1111/j.1365-2478.2011.00967.x

Xia

Zhou

Chen

Zhang

2018

A rectangular-grid lattice spring model for modeling elastic waves in Poisson's solids

Geophysics

(

T69

–

T86

10.1190/geo2016-0414.1

Xia

Zhou

Chen

Wang

2017

A general 3D lattice spring model for modeling elastic waves

Bull. seism. Soc. Am.

107

(

2194

–

2212

Chen

Zhang

Pan

2022

An adaptive modal discontinuous Galerkin finite element parallel method using unsplit multi-axial perfectly matched layer for seismic wave modeling

Commun. Comput. Phys.

(

1083

–

1113

10.4208/cicp.OA-2021-0118

Yin

Wang

Zhang

2020

Effect of particle shape on the progressive failure of shield tunnel face in granular soils by coupled FDM–DEM method

Tunn. Undergr. Space Technol.

100

103 394

10.1016/j.tust.2020.103394

Zhang

2022

Efficient 2D acoustic wave finite-difference numerical simulation in strongly heterogeneous media using the adaptive mesh refinement technique

Geophysics

(

T29

–

T42

10.1190/geo2020-0801.1

Zhang

Wang

Zhang

2020

A procedure for 3D seismic simulation from rupture to structures by coupling SEM and FEM

Bull. seism. Soc. Am.

110

(

1134

–

1148

Zhang

Chen

2006

Traction image method for irregular free surface boundaries in finite difference seismic wave simulation

Geophys. J. Int.

167

(

337

–

353

10.1111/j.1365-246X.2006.03113.x

Zhang

Shen

2010

Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling

Geophysics

(

T141

–

T154

Zhao

2015

Modelling 3D jointed rock masses using a lattice spring model

Int. J. Rock Mech. Min. Sci.

–

10.1016/j.ijrmms.2015.05.011

Zhao

Fang

Zhao

2011

A 3D distinct lattice spring model for elasticity and dynamic failure

Int. J. Numer. Anal. Meth. Geomech.

(

859

–

885

Zhu

Dorman

2000

Two-dimensional, three-component wave propagation in a transversely isotropic medium with arbitrary-orientation—finite-element modeling

Geophysics

(

934

–

942

APPENDIX A: DERIVATIONS OF CALIBRATING THE SPRING COEFFICIENTS THROUGH THE ELASTICITY TENSOR

Comparing eq. (3) with eq. (13), we have

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0x}}}{{{{\left| {r_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xx}{b}_{i0j} + g_{ik}^{xx}{b}_{i0k}} \right)} \right]r_{0i}^{0x} = \frac{{\partial {c}_{11}}}{{\partial x}} + \frac{{\partial {c}_{15}}}{{\partial z}}, \end{eqnarray}$$

(A1)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0x}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xx}{b}_{i0j} + g_{ik}^{xx}{b}_{i0k}} \right)} \right]r_{0i}^{0z} = \frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}, \end{eqnarray}$$

(A2)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0x}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xx}{b}_{i0j} + g_{ik}^{xx}{b}_{i0k}} \right)} \right]\frac{1}{2}{(r_{0i}^{0x})}^2 = {c}_{11}, \end{eqnarray}$$

(A3)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0x}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xx}{b}_{i0j} + g_{ik}^{xx}{b}_{i0k}} \right)} \right]\frac{1}{2}{(r_{0i}^{0z})}^2 = {c}_{55}, \end{eqnarray}$$

(A4)

(A5)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xz}{b}_{i0j} + g_{ik}^{xz}{b}_{i0k}} \right)} \right]r_{0i}^{0x} = \frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}, \end{eqnarray}$$

(A6)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xz}{b}_{i0j} + g_{ik}^{xz}{b}_{i0k}} \right)} \right]r_{0i}^{0z} = \frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{35}}}{{\partial z}}, \end{eqnarray}$$

(A7)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xz}{b}_{i0j} + g_{ik}^{xz}{b}_{i0k}} \right)} \right]\frac{1}{2}{(r_{0i}^{0x})}^2 = {c}_{15}, \end{eqnarray}$$

(A8)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xz}{b}_{i0j} + g_{ik}^{xz}{b}_{i0k}} \right)} \right]\frac{1}{2}{(r_{0i}^{0z})}^2 = {c}_{35}, \end{eqnarray}$$

(A9)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{xz}{b}_{i0j} + g_{ik}^{xz}{b}_{i0k}} \right)} \right]r_{0i}^{0x}r_{0i}^{0z} = {c}_{13} + {c}_{55}, \end{eqnarray}$$

(A10)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zx}{b}_{i0j} + g_{ik}^{zx}{b}_{i0k}} \right)} \right]r_{0i}^{0x} = \frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}}, \end{eqnarray}$$

(A11)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zx}{b}_{i0j} + g_{ik}^{zx}{b}_{i0k}} \right)} \right]r_{0i}^{0z} = \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{35}}}{{\partial z}}, \end{eqnarray}$$

(A12)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zx}{b}_{i0j} + g_{ik}^{zx}{b}_{i0k}} \right)} \right]\frac{1}{2}{(r_{0i}^{0x})}^2 = {c}_{15}, \end{eqnarray}$$

(A13)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zx}{b}_{i0j} + g_{ik}^{zx}{b}_{i0k}} \right)} \right]\frac{1}{2}{(r_{0i}^{0z})}^2 = {c}_{35}, \end{eqnarray}$$

(A14)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0x}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zx}{b}_{i0j} + g_{ik}^{zx}{b}_{i0k}} \right)} \right]r_{0i}^{0x}r_{0i}^{0z} = {c}_{13} + {c}_{55}, \end{eqnarray}$$

(A15)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0z}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zz}{b}_{i0j} + g_{ik}^{zz}{b}_{i0k}} \right)} \right]r_{0i}^{0x} = \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{35}}}{{\partial z}}, \end{eqnarray}$$

(A16)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0z}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zz}{b}_{i0j} + g_{ik}^{zz}{b}_{i0k}} \right)} \right]r_{0i}^{0z} = \frac{{\partial {c}_{35}}}{{\partial x}} + \frac{{\partial {c}_{33}}}{{\partial z}}, \end{eqnarray}$$

(A17)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0z}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zz}{b}_{i0j} + g_{ik}^{zz}{b}_{i0k}} \right)} \right]\frac{1}{2}{(r_{0i}^{0x})}^2 = {c}_{55}, \end{eqnarray}$$

(A18)

$$\begin{eqnarray} \frac{1}{{{h}^2}}\sum {\rm{\ }}\left[ {\frac{{r_{0i}^{0z}r_{0i}^{0z}}}{{{{\left| {{\boldsymbol{r}}_{0i}^0} \right|}}^2}}{k}_{0i} + {{\sin }}^2\theta _{i0j}^0\left( {g_{ij}^{zz}{b}_{i0j} + g_{ik}^{zz}{b}_{i0k}} \right)} \right]\frac{1}{2}{(r_{0i}^{0z})}^2 = {c}_{33} \end{eqnarray}$$

(A19)

and

(A20)

The above equations must be solved at all particles to obtain spring coefficients from the elasticity tensor. When solving eqs (A1)–(A20), we find that |${c}_{15}$| must be equal to |${c}_{35}$|⁠, indicating that this model is not suitable for general 2D anisotropic media. Liu & Liu (2006) also observed a similar phenomenon that the nine-disc model with hexagon lattices can solve orthotropic media, while the seven-disc model with square lattices can solve general anisotropic media. In the meantime, the number of the unknown spring coefficients is 16 (i.e. 8 linear springs and 8 angular springs), while the number of the linearly independent equations above is 14, which means the parameter options are not unique. To keep their physical meaning, Xia et al. (2017) chose the spring coefficients that are larger than zero. Drawing on the previous research works (Hu & Jia 2016; Xia et al. 2017), according to the central symmetry, we add two equations to make the solution unique, that is,

$$\begin{eqnarray} {b}_{102} + {b}_{203} = {b}_{506} + {b}_{607} \end{eqnarray}$$

(A21)

and

$$\begin{eqnarray} {b}_{304} + {b}_{405} = {b}_{708} + {b}_{801}. \end{eqnarray}$$

(A22)

Finally, the spring coefficients are obtained as

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{l}@{}} \displaystyle {{k}_{01} = {c}_{11} - \frac{{{c}_{13}}}{2} - \frac{{{c}_{55}}}{2} + \frac{h}{2}\left( {\frac{{\partial {c}_{11}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial x}} - 2\frac{{\partial {c}_{55}}}{{\partial x}}} \right)}\\ \displaystyle {{k}_{02} = \frac{{{c}_{13}}}{2} + {c}_{15} + \frac{{{c}_{55}}}{2} + \frac{h}{2}\left( {\frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{15}}}{{\partial z}} + \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {{k}_{03} = {c}_{33} - \frac{{{c}_{13}}}{2} - \frac{{{c}_{55}}}{2} + \frac{h}{2}\left( {\frac{{\partial {c}_{13}}}{{\partial z}} + \frac{{\partial {c}_{33}}}{{\partial z}} - 2\frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {{k}_{04} = \frac{{{c}_{13}}}{2} - {c}_{15} + \frac{{{c}_{55}}}{2} + \frac{h}{2}\left( {\frac{{\partial {c}_{15}}}{{\partial x}} - \frac{{\partial {c}_{15}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {{k}_{05} = {c}_{11} - \frac{{{c}_{13}}}{2} - \frac{{{c}_{55}}}{2} - \frac{h}{2}\left( {\frac{{\partial {c}_{11}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial x}} - 2\frac{{\partial {c}_{55}}}{{\partial x}}} \right)}\\ \displaystyle {{k}_{06} = \frac{{{c}_{13}}}{2} + {c}_{15} + \frac{{{c}_{55}}}{2} - \frac{h}{2}\left( {\frac{{\partial {c}_{15}}}{{\partial x}} + \frac{{\partial {c}_{15}}}{{\partial z}} + \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {{k}_{07} = {c}_{33} - \frac{{{c}_{13}}}{2} - \frac{{{c}_{55}}}{2} - \frac{h}{2}\left( {\frac{{\partial {c}_{13}}}{{\partial z}} + \frac{{\partial {c}_{33}}}{{\partial z}} - 2\frac{{\partial {c}_{55}}}{{\partial z}}} \right)}\\ \displaystyle {{k}_{08} = \frac{{{c}_{13}}}{2} - {c}_{15} + \frac{{{c}_{55}}}{2} - \frac{h}{2}\left( {\frac{{\partial {c}_{15}}}{{\partial x}} - \frac{{\partial {c}_{15}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \end{array}} \right. \end{eqnarray}$$

(A23)

and

$$\begin{eqnarray} \left\{ {\begin{array}{@{}*{1}{l}@{}} \displaystyle {{b}_{102} = - {h}^2\left[ {{c}_{13} - {c}_{55} + h\left( {\frac{{\partial {c}_{13}}}{{\partial x}} - \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \right]}\\ \displaystyle {{b}_{203} = - {h}^2\left[ {{c}_{13} - {c}_{55} - h\left( {\frac{{\partial {c}_{13}}}{{\partial x}} - \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \right]}\\ \displaystyle {{b}_{304} = - {h}^2\left[ {{c}_{13} - {c}_{55} + h\left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} - \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \right]}\\ \displaystyle {{b}_{405} = - {h}^2\left[ {{c}_{13} - {c}_{55} - h\left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} - \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \right]}\\ \displaystyle {{b}_{506} = - {h}^2\left[ {{c}_{13} - {c}_{55} - h\left( {\frac{{\partial {c}_{13}}}{{\partial x}} - \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \right]}\\ \displaystyle {{b}_{607} = - {h}^2\left[ {{c}_{13} - {c}_{55} + h\left( {\frac{{\partial {c}_{13}}}{{\partial x}} - \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} + \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \right]}\\ \displaystyle {{b}_{708} = - {h}^2\left[ {{c}_{13} - {c}_{55} - h\left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} - \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \right]}\\ \displaystyle {{b}_{801} = - {h}^2\left[ {{c}_{13} - {c}_{55} + h\left( {\frac{{\partial {c}_{13}}}{{\partial x}} + \frac{{\partial {c}_{13}}}{{\partial z}} - \frac{{\partial {c}_{55}}}{{\partial x}} - \frac{{\partial {c}_{55}}}{{\partial z}}} \right)} \right]} \end{array}} \right., \end{eqnarray}$$

(A24)

where a second-order central FD operator is used to calculate the first spatial derivative of the medium.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Download all slides

Month:	Total Views:
December 2024	39
January 2025	162
February 2025	67
March 2025	52
April 2025	34

Article Contents

A new lattice spring model for seismic wave simulation in heterogeneous anisotropic media

SUMMARY

1 INTRODUCTION

2 LATTICE SPRING MODEL

2.1 Basic lattice units in the LSM

2.2 Calibrated continuum elastic theory

2.3 Parameter transformation

2.4 Verlet algorithm

3 DISPERSION ANALYSIS

4 STABILITY CONDITION

5 NUMERICAL EXAMPLES

5.1 Scattering model

5.2 Layered model

5.3 SEAM model

5.4 Fracture model

6 CONCLUSIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

REFERENCES

APPENDIX A: DERIVATIONS OF CALIBRATING THE SPRING COEFFICIENTS THROUGH THE ELASTICITY TENSOR

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

A new lattice spring model for seismic wave simulation in heterogeneous anisotropic media

SUMMARY

1 INTRODUCTION

2 LATTICE SPRING MODEL

2.1 Basic lattice units in the LSM

2.2 Calibrated continuum elastic theory

2.3 Parameter transformation

2.4 Verlet algorithm

3 DISPERSION ANALYSIS

4 STABILITY CONDITION

5 NUMERICAL EXAMPLES

5.1 Scattering model

5.2 Layered model

5.3 SEAM model

5.4 Fracture model

6 CONCLUSIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

REFERENCES

APPENDIX A: DERIVATIONS OF CALIBRATING THE SPRING COEFFICIENTS THROUGH THE ELASTICITY TENSOR

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only