Wave propagation across fluid-solid interfaces with LBM-LSM coupling schemes

Abstract

Seismic wave propagation in fluid-solid coupled media is currently a popular topic. However, traditional wave equation-based simulation methods have to consider complex boundary conditions at the fluid-solid interface. To address this challenge, we propose a novel numerical scheme that integrates the lattice Boltzmann method (LBM) and lattice spring model (LSM). In this scheme, LBM simulates viscoacoustic wave propagation in the fluid area and LSM simulates elastic wave propagation in the solid area. We also introduce three different LBM-LSM coupling strategies, a standard bounce back scheme, a specular reflection scheme, and a hybrid scheme, to describe wave propagation across fluid-solid boundaries. To demonstrate the accuracy of these LBM-LSM coupling schemes, we simulate wave propagation in a two-layer model containing a fluid-solid interface. We place excitation sources in the fluid layer and the solid layer respectively, to observe the wave phenomena when seismic waves propagate to interface from different sides. The simulated results by LBM-LSM are compared with the reference wavefields obtained by the finite difference method (FDM) and the analytical solution (ANA). Our LBM-LSM coupling scheme was verified effective, as the relative errors between the LBM-LSM solutions and reference solutions were within an acceptable range, sometimes around 1.00%. The coupled LBM-LSM scheme is further used to model seismic wavefields across a more realistic rugged seabed, which reveals the potential applications of the coupled LBM-LSM scheme in marine seismic imaging techniques, such as reverse-time migration and full-waveform inversion. The method also has potential applications in simulating wave propagation in complex twoand multi-phase media.

Headnote

ABSTRACT

Seismic wave propagation in fluid-solid coupled media is currently a popular topic. However, traditional wave equation-based simulation methods have to consider complex boundary conditions at the fluid-solid interface. To address this challenge, we propose a novel numerical scheme that integrates the lattice Boltzmann method (LBM) and lattice spring model (LSM). In this scheme, LBM simulates viscoacoustic wave propagation in the fluid area and LSM simulates elastic wave propagation in the solid area. We also introduce three different LBM-LSM coupling strategies, a standard bounce back scheme, a specular reflection scheme, and a hybrid scheme, to describe wave propagation across fluid-solid boundaries. To demonstrate the accuracy of these LBM-LSM coupling schemes, we simulate wave propagation in a two-layer model containing a fluid-solid interface. We place excitation sources in the fluid layer and the solid layer respectively, to observe the wave phenomena when seismic waves propagate to interface from different sides. The simulated results by LBM-LSM are compared with the reference wavefields obtained by the finite difference method (FDM) and the analytical solution (ANA). Our LBM-LSM coupling scheme was verified effective, as the relative errors between the LBM-LSM solutions and reference solutions were within an acceptable range, sometimes around 1.00%. The coupled LBM-LSM scheme is further used to model seismic wavefields across a more realistic rugged seabed, which reveals the potential applications of the coupled LBM-LSM scheme in marine seismic imaging techniques, such as reverse-time migration and full-waveform inversion. The method also has potential applications in simulating wave propagation in complex twoand multi-phase media.

Keywords:

Lattice Boltzmann method (LBM)

Lattice spring model (LSM)

LBM-LSM coupling

Finite difference method (FDM)

Fluid-solid configuration

Seismic wave simulation

(ProQuest: ... denotes formulae omitted.)

1. Introduction

Wave propagation across the fluid-solid interfaces is a common and vital topic in many engineering application fields (e.g., marine seismic survey, borehole sonic log, non-destructive testing, rock physics; Zhang, 2004; Zhan et al., 2020), and can be used in full-waveform inversion (FWI) (e.g., Liu et al., 2013; Shi et al., 2014; Liu et al., 2019; Wang et al., 2019; Qu et al., 2020; Aghamiry et al., 2020; Yao et al., 2020; Zhang et al., 2020; Wang et al., 2022b; Fang et al., 2023; Li et al., 2023) and reverse-time migration (RTM) (e.g., Yu et al., 2016; Li et al., 2019, 2021; Mao et al., 2021; Li and Qu, 2022; Zhou et al., 2022; Wang et al., 2022d; Qu et al., 2023; Mu et al., 2023; Li et al., 2024). To address these problems, a lot of simulation approaches have been proposed, such as finite-difference method (FDM) (Virieux, 1986; Graves, 1996; van Vossen et al., 2002; Zhou et al., 2005; Qu et al., 2016; Chen et al., 2017, 2019; Wu et al., 2018; Zou et al., 2020; Wang et al., 2022a, 2022c; Ren et al., 2022), and finite-element method (Lysmer and Drake, 1972; Serón et al., 1990, 1996). Two common computational methods are frequently employed to calculate the wavefields in fluid-solid coupled media. The partitioned approach uses one numerical calculation method to solve the acoustic wave in the fluid region and another simulation method to calculate the elastic wave in the solid region. The coupling of the fluid-solid boundary is achieved by introducing appropriate explicit boundary conditions that enable seismic waves to propagate throughout the computational region. However, this approach usually requires detailed information such as the specific location and morphology of the fluid-solid interface to be known (Yu et al., 2016). van Vossen et al. (2002) proposed a mirroring method to deal with the fluid-solid boundary for the wavefield error problem caused by the deviation of the interface position in the FDM simulation. It is shown that this method can suppress the numerical dispersion caused by the tilted interface to a certain extent, but it is not suitable for irregular fluid-solid interfaces. Carcione and Helle (2004) simulated the seismic wave crossing through a fluid and viscoelastic solid coupling media based on domain decomposition techniques and also analyzed the waves associated with the fluid-solid interface using amplitude variation with offset (AVO). Wilcox et al. (2010) introduced a high-order discontinuous Galerkin scheme for the numerical solution of three-dimensional wave propagation problems in coupled elastic-acoustic media. Soares (2011) summarized the common partitioned coupling methods for computing fluid-solid coupling problems and found that such methods usually require the use of domain decomposition techniques where either common interface nodes or partially or fully overlapping computational grids are used at the interface, and then the direct transfer process of the wavefield in two different regions is realized with the help of some explicit direct coupling method or implicit iterative coupling method. To better fit with the topographic interface on the sea bottom, Sun et al. (2021) proposed a curvilinear-grid FDM to implement the fluid-solid boundary conditions. Lecoulant et al. (2022) employed the spectral finite-element scheme to simulate the low-frequency seismo-acoustic waves generated by submarine earthquakes in the ocean and found they are comparable to the theoretical ones.

While conventional numerical schemes can enhance the computational accuracy of wave propagation near interfaces by incorporating boundary conditions (such as normal displacement and traction continuity for horizontal interfaces) and averaging elastic moduli and densities, they fall short in providing an intuitive depiction of wave scattering at irregular interfaces. Moreover, these methods rely on continuity assumptions about the wave equation and may impede the flexibility in modeling interfacial wave phenomena.

The lattice Boltzmann method (LBM) is a microscopic dynamic simulation method based on statistical physics. It is suitable for modeling various physical phenomena in fluids without solving the traditional wave equations or the Navier-Stokes equations. LBM achieves the purpose of simulating macroscopic complex physical phenomena by tracking the simple collision and streaming process between microscopic particles. LBM originated from the Navier-Stokes equation and was first proposed by McNamara and Zanetti (1988) to realize the numerical simulation of flow field. Due to its flexible boundary processing, simple calculation, easy parallelism, and complete discrete characteristics, it has a certain development in the field of wave propagation. Chopard et al. (1998) studied the wave and fracture evolution processes in solid systems based on the LBM model with the Bhatnagar-Gross-Krook collision operator. Xu and Sagaut (2011) proposed an ideal strategy for reducing dispersion and dissipation errors of the LBM model for simulating aeroacoustics problems. Li and Shan (2011) studied the propagation and attenuation laws of adiabatic acoustic waves based on the LBM theory. Frantziskonis (2011) and O'Brien et al. (2012) successively tried to use LBM to simulate the propagation of elastic waves and discussed the dispersion and computational stability issues in numerical simulations, respectively. Xia et al. (2017a, 2022) modified LBM to simulate wave propagation in viscous media, revealing a quantitative relationship between the relaxation time in LBM and the quality factor in the Kelvin-Voigt-based wave equation. Escande et al. (2020) proposed von Neumann stability analysis and Chapman-Enskog analysis for LBM for wave propagation in isotropic linear elastic solids, and successfully calculated the surface wave in elastic solids by LBM for the first time. Jiang et al. (2022) proposed a set of stability models for the multiple-relaxation-time lattice Boltzmann model and verified by the seismic wavefield simulation experiments. Although some experts have applied LBM to the numerical simulation of elastic waves, but the relevant research is not yet mature, so LBM is more suitable for the simulation of waves in fluids.

The lattice spring model (LSM) is a non-wave equation modeling scheme used for numerical simulations of dynamic wave propagation or quasi-static deformation in solid media. Toomey and Bean (2000) proposed a discrete particle model where a continuous solid media is represented as closely spaced microscopic spheres arranged in a hexagonal structure. The numerical simulation of elastic waves in Poisson solid was achieved through Hooke's law and the Verlet algorithm, which is known as the discrete particle scheme (DPS) or elastic lattice method (ELM). However, the Poisson's ratio of the media described by the DPS model is fixed at 0.25, which limits its application in certain fields. O'Brien and Bean (2004) added a non-central force term to the ELM model proposed by Toomey and Bean (2000) to get rid of the restriction on the Poisson's ratio in the ELM model. Pazdniakou and Adler (2012) introduced two groups of angular springs (i.e., 45° and 60°) into the original LSM model to improve its flexibility and adaptability, and the idea of elastic elements was also incorporated in the new model to realize wavefield simulation in complex media like porous media. Hu and Jia (2016, 2018) modified the ELM model and termed it as dynamic lattice method (DLM) by making use of the orientations of the bonds in the lattice and by including more adjacent nodes or bonds to increase the calculation accuracy, and the proposed DLM scheme was employed to simulate wave propagation in heterogeneous TI media and further to be tested in seismic imaging. Xia et al. (2015, 2017b, 2018) combined LSM with the Verlet algorithm and proposed a new method to simulate elastic wave propagation in subsurface media. The proposed method makes the mesh generation more flexible, which means a wide potential application. As a promising simulation technique in complex solid/porous media, the non-wave equation-based LSM characterizes the elastoplasticity from a micromechanics perspective thus enjoys good elastic system dynamics and can effectively simulate waveform propagation in complex media. More recently, Liu et al. (2020) proposed a modified LSM that can be used to analyze and characterize the size-dependent effect on wave propagation in an elastic media with a higher Poisson's ratio by comparison with the modified couple stress theory. Wei (2021) studied fatigue cracking in homogeneous and composite materials by using a similar lattice particle method and achieved reasonable computational accuracy and efficiency. Tang et al. (2022) proposed a novel perfectly matched layer (PML) boundary condition for LSM to model wave propagation in complex media. It can be seen that LSM or its derivative model is very suitable for simulating the wave problem in solid media.

The collaboration between LBM and LSM holds promise in providing a comprehensive and versatile approach to address a wide range of wave simulation challenges across diverse media and scenarios. To simulate seismic wave propagation in geological models with significant variations in media properties, the development of accurate and stable methods for solving coupled fluid-solid boundary models is a crucial research area in exploration seismology and other related fields. Ladd and Verberg (2001) proposed a numerical calculation method that can simulate the hydrodynamic interaction between suspended solid particles by modifying the classical lattice Boltzmann equation and using the idea of momentum exchange to integrate the forces between solids and fluids into the updated formulation of the LBM. O'Brien and Bean (2004) implemented the coupling between the LBM and LSM considering only the action of the wire spring with the help of momentum exchange algorithm. This technique was used to solve the mechanical coupling problem of the deformation of the fluid-rock system, and the effectiveness of the coupling method was verified by using the reciprocity theorem. Pazdniakou and Adler (2013) developed a novel LBM-LSM coupling model based on O'Brien and Bean (2004), in which the LSM not only considers the conventional linear spring but also incorporates angular springs, allowing the Poisson's ratio of the simulated solid media to be variable within a certain range. Adler and Pazdniakou (2018) further employed this LBM-LSM coupling method to study the corresponding seismic wavefields of complex geological models containing pores and compared the differences between the corresponding wavefields with and without pores under different model assumptions. Jiang and Zhao (2019) and Li et al. (2022) developed a similar fluid-solid coupling scheme based on the LBM and a distinct LSM, in which the coupled model shows its ability to mimic hydraulic fracturing in formations with complex fracture networks. Nguyen et al. (2021) numerically investigated the acoustic properties of dry or saturated Fontainebleau sandstone porous media using a combination of LSM and LBM. Wang et al. (2023) reviewed the recent progress of LBM-solid coupling schemes, which provide new directions for wave simulation with advanced LBM-LSM coupling schemes. As far as we know, the commonly used LBM-LSM-based fluid-solid boundary coupling schemes for wave simulations are mainly carried out through combination of the standard bounce back strategy and the momentum exchange algorithm. Under this assumption, the fluid particles will be bounced back along the original path when they hit the solid boundary, while the momentum carried by the fluid particles will be transferred to the solid nodes at the interface, thus triggering the deformation of the solid interface, so this type of fluid-solid boundary coupling method is called momentum exchange algorithm.

Several coupling techniques between LBM and LSM have been proposed, but their applications have mainly been in fluid dynamics and not directly in numerical simulations of wavefields. It is essential to note that the two methods are not interchangeable; instead, they complement each other and can work in tandem to tackle different aspects of wave propagation problems effectively. In this paper, we focus on simulating wavefields in fluid and solid media using LBM and LSM, respectively, and explore a feasible and effective mutual coupling technique for studying seismic wave propagation in two-phase media with fluid-solid interfaces.

The rest parts of this paper are as follows. In Section 2, we give the basic theory of the LBM method, LSM method, and the LBM-LSM coupling method in detail. In Section 3, we verify the effectiveness of the proposed fluid-solid coupling schemes in two-layer fluid-solid coupled media and apply it to a complex seabed model. At last, we make discussions and reach some conclusions.

2. LBM-LSM coupling schemes

The lattice Boltzmann method (LBM) is a computational fluid dynamics technique used to simulate and model fluid flow and other complex physical phenomena. It is a mesoscopic method, meaning it operates at an intermediate level between the macroscopic continuum and the microscopic particle-based approaches. The lattice spring model (LSM) is a discrete model used to study the mechanical behavior of materials, particularly in the context of solid mechanics and material science. It is a mesoscopic approach that bridges the gap between the macroscopic continuum and the microscopic atomistic models. The frameworks of LBM and LSM are shown in Appendix A and Appendix B, respectively.

Based on the LBM-LSM coupling method to simulate wavefields in complex media with fluid-solid dual phases, the most critical aspect is to handle the propagation of seismic waves at the fluid-solid boundary. The most commonly employed strategy for coupling LBM and LSM at the fluid-solid boundary is to discretize the entire computational domain using the same set of spatial grids and to differentiate which nodes represent solids and which represent fluids by defining Boolean variables or phase functions. The schematic diagram of LBM-LSM fluid-solid coupling method is shown in Fig. 1, where the red solid dots represent LSM solid nodes, and the blue hollow circles represent LBM fluid nodes. The fluid-solid boundary is generally located at the midpoint of the line connecting the solid node and the adjacent fluid node in the vicinity of the boundaries, which is labeled with the green dashed line in Fig. 1. Suppose x represents the fluid boundary node coordinate, and ... represents the adjacent solid boundary node coordinate, which are located on both sides of the boundary. When the fluid-solid boundary is fixed, the standard bounce back scheme is expressed as

... (1)

where -i is the opposite direction of i. On the other hand, when the fluid-solid boundary is moving and both the boundary velocity and solid node velocity are v , the standard bounce back scheme is expressed as (O'Brien and Bean, 2004; Buxton et al., 2005)

... (2)

where, ... is a constant. In this case, due to the presence of the fluid node, the total force acting on the LSM solid nodes is

... (3)

where

... (4)

N in the summation sign of Eq. (3) is the summation over the adjacent linear springs, and M in the summation sign is the summation over the adjacent fluid nodes. In other words, the part before the plus sign in Eq. (3) is the force exerted by the neighboring solid nodes, and the part after the plus sign is the force exerted by the neighboring fluid. K is the spring elastic constant and gn is the coupling constant. The final elastic wave equation can be ensured to be isotropic by choosing the value of gn reasonably. uj and un are the displacements of nodes xj and xn , respectively, and ... is the unit vector pointing from node xj to node xn.

The standard bounce back scheme is commonly used in literature when fluid particles come into contact with the fluid-solid interface. However, this scheme may not be appropriate for seismic waves that cross the wave impedance interface, as the reflected wave is usually reflected back to the fluid region as a specular reflection. Furthermore, actual subsurface strata are not smooth and regular partition interfaces, so seismic waves are likely to be diffusely reflected when they cross such uneven wave impedance interfaces. Therefore, two different LBM-LSM fluid-solid coupling schemes, namely the specular reflection scheme and the hybrid scheme, have been developed and utilized in this work.

The sketch maps of the two typical LBM-LSM boundary conditions, standard bounce back scheme and specular reflection scheme, are shown in Fig. 2. For simplicity, we assume that xf represents the coordinates of the fluid boundary node above the interface and xs represents the coordinates of the solid boundary node below the interface, which are located on both sides of the boundary. Suppose at a certain moment, a fluid particle streaming along direction i at node xf will encounter a solid node xs . If the particle has a standard bounce back, the particle will stream along direction -i at node xf at the next moment, and if the particle has a specular reflection, the particle will stream along direction +i (see Fig. 2). The following describes the basic principles for LBM-LSM coupling schemes in three cases:

1) Standard bounce back scheme

Assuming that the motion velocity of the solid node xs at the fluid-solid boundary is vs , when the fluid particle hits the fluid-solid boundary with velocity cf and a standard bounce back occurs, consider that the streaming velocity of the LBM fluid particle does not change, and according to the momentum exchange algorithm, the particle distribution function along the standard bounce back direction -i at node xf after the collision can be expressed as

... (5)

where ... is the particle distribution function along the direction i after the streaming step of grid node xf, ci is the particle streaming velocity along the i direction, wi is the LBM discrete weighting factor, and cs is the LBM lattice speed of sound. Since the forces are reciprocal, the force exerted by a fluid particle on a neighboring solid node xs can be expressed as

... (6)

in which, * and * are the time and space sampling steps, respectively.

2) Specular reflection scheme

Similarly, according to the momentum exchange algorithm, the particle distribution function along the specular reflection direction +i at the node xf after the collision can be expressed as

... (7)

Accordingly, the force exerted by the fluid particles on the adjacent solid node xs can be expressed as

... (8)

3) Hybrid scheme

According to the fluid-solid boundary momentum exchange algorithm, assuming that part of the fluid particles on the node xf after encountering the fluid-solid boundary (accounting for a ) moves along the standard bounce back direction -i, and the other part (accounting for β ) moves along the specular reflection direction +i . In this case, the particle distribution function on the node xf after the collision step can be expressed as

... (9)

where the scaling factors a and β are required to satisfy the following conditions

... (10)

Consequentially, the force of the fluid particles at the boundary on the neighboring solid particles can be expressed as

... (11)

With such LSM and LBM fluid-solid coupling boundary conditions, the seismic wavefield can propagate in both directions between the fluid and solid regions.

In summary, the simulation of fluid-solid coupling wavefields with the LBM-LSM scheme can be carried out in three steps as shown in Fig. 3 and described below.

1) LSM iteration. The first step involves the LSM iteration, which includes locating all solid grid points and calculating linear and angular spring forces using Eq. (B-7), and updating relevant displacements, velocities, and accelerations using the velocity Verlet algorithm in Eq. (B-8).

2) LBM iteration. The second step involves the LBM iteration, which includes locating all fluid grid points, performing the collision step using the collision operator in Eq. (A-3), and performing the streaming step using the improved streaming step in Eq. (A-8). This step is carried out differently depending on whether the next node is fluid or solid. If it is fluid, then streaming is performed directly, while if it is solid, then the bounce-back occurs.

3) LBM and LSM coupling calculation. The third step involves the coupling calculation of LBM and LSM. In this step, the vibration velocity of a solid node is assigned to a neighboring fluid node using interpolation or a weighting factor multiplication if the neighboring node is fluid. The fluid forces on solids are transferred to neighboring solid nodes through momentum exchange or interpolation algorithms if the neighboring node is solid after the LBM calculation.

3. Numerical examples

3.1. 2D layered models

To examine whether aforementioned fluid-solid coupling method is feasible, this paper will first adopt a relatively simple two-layer fluid-solid coupling model to conduct wave simulation tests, with the upper layer of the model being the fluid region and the lower layer the solid region. The P-wave velocity in the fluid region is 1155 m/s and the density is 1000 kg/m3; the P-wave and S-wave velocities in the solid region are 1443 m/s and 962 m/s, respectively, and the density is 2000 kg/m3. The dominant frequency of the Ricker wavelet source used in the simulation is 15 Hz, and the delay time is 0.08 s. The position of the source is (x = 500 m, z = 450 m), when the source is added on a node in the fluid region, the depth of the fluid-solid interface is set to 550 m; when the source is added to the solid region, the depth of the fluid-solid interface is set to 300 m (see Fig. 4). The spatial sampling interval is 1.0 m, the time sampling interval for the FDM results of the Kelvin-Voigt model (denoted as FDM-KV) (Madja et al., 1985) is 0.125 ms, and the time sampling interval for the LBM-LSM scheme is 0.50 ms. The number of discretization grids used in the entire computing area is 1001 × 1001. In the numerical simulation tests in this section, the LSM uses a square grid considering both linear and angular spring forces, while the LBM uses a relatively simple but typical single-relaxation-time LBM model with the relaxation time to be 0.6 or 1.0. Three distinct LBM-LSM coupling schemes are tested as follows.

3.1.1. Standard bounce back scheme

First, consider the case where a fluid particle encounters the fluid-solid boundary with a standard bounce back, i.e., the fluid particle bounces back to the fluid region along the original direction of incidence (see Fig. 4). In order to examine the application of this fluid-solid coupling method more comprehensively, the cases of seismic waves propagating from the fluid region to the fluid-solid interface and then incident into the solid region are simulated (the simulation results are shown in Fig. 5, Fig. 6, Fig. 7), and the cases of seismic waves propagating from the solid region to the fluid-solid interface and then incident into the fluid region (the simulation results are shown in Fig. 8, Fig. 9). The wave profiles calculated by the LBM-LSM coupling method along the depth direction at x = 400 m in Fig. 6, Fig. 9 are extracted in Fig. 5, Fig. 8, respectively, for comparison with the FDM-KV solution and an analytical solution for the wavefield reflected at a fluid-solid boundary (denoted as ANA) (de Hoop and van der Hijden, 1983). The relative errors are calculated with the least-squares misfit function defined in Xia et al. (2018).

According to Fig. 5, Fig. 6, Fig. 7, Fig. 8, Fig. 9, it can be seen that the seismic wavefields solved by the LBM-LSM coupling method match well with those calculated by FDM-KV with the same simulation parameter settings when the standard bounce back scheme is used for the fluid-solid coupling interface of LBM-LSM, except for some minor differences in the details of the wavefields at the interface as shown in the residual wavefields between LBM-LSM and FDM-KV. The main morphology and polarity characteristics of the wavefields simulated by both methods are in agreement. Both methods agree well regardless of whether the source is loaded on a fluid node or a solid node.

3.1.2. Specular reflection scheme

Next, consider the case where a fluid particle hits the fluid-solid boundary and undergoes a specular reflection, i.e., the fluid particle bounces back to the fluid region along the mirror-symmetric direction. The simulation parameters are consistent with subsection 3.1.1, except for the interface coupling scheme. Without loss of generality, Fig. 10, Fig. 13 contrast the snapshots calculated by LBM-LSM and FDM-KV for the seismic source added on the fluid node and the solid node, respectively. Fig. 11, Fig. 14 abstract the corresponding wave profiles in the depth direction. Moreover, the seismic traces for LBM-LSM, FDM-KV and ANA recorded at (900 m, 700 m) for the two-layer fluid-solid coupling model with the source added on fluid LBM node are compared in Fig. 12, which shows the difference between LBM-LSM/FDM-KV and the analytical solution is a little large because the analytical solution is fully elastic, while the adopted LBM-LSM coupling and FDM-KV schemes take into account a certain viscosity (i.e., τ = 1.00). It is easy to conclude from Fig. 10, Fig. 11, Fig. 12, Fig. 13, Fig. 14 that the wavefields of the z-component of the vibration velocity calculated by the LBM-LSM agree well with the FDM-KV results, while the simulated results of the x-component differ more significantly from the FDM-KV reference solution. The reason why there is a good match of z-component of the vibration velocity and a bad match of x-component of the vibration velocity is that the velocity of the z-component of the bounced particle is almost constant and the velocity of the x-component of the bounced particle is reversed when changing from the standard bounce back scheme to the specular reflection scheme. By comparing the results in this subsection with those in the previous subsection, it is easy to find that the difference in the motion of fluid particles at the interface when implementing LBM-LSM fluid-solid coupling scheme can directly lead to the difference between the simulated seismic wavefields.

3.1.3. The hybrid scheme

In this subsection, a hybrid scheme is employed to simulate the motions of fluid particles at the fluid-solid interface, and for simplicity, the scaling factors a and β in Eq. (9) are both chosen to be 0.50. Similarly, the snapshots and depth profiles of the LBM-LSM and FDM-KV are also compared here. As can be seen from the snapshots in Fig. 15, the snapshots of the LBM-LSM using the hybrid method have similar transmitted/reflected Pand S-wave patterns with those of FDM-KV, but still there are deviations between the LBM-LSM calculated wave profiles and the FDM-KV results. Similar conclusions can be drawn from the results of the seismic profiles in Fig. 16.

This section highlights the differences between wavefields obtained from simulations using the three different LBM-LSM fluid-solid interface treatment strategies. Table 1 presents the relative errors between the profiles or snapshots of LBM-LSM and FDM-KV. The specular reflection scheme shows the best agreement with FDM-KV, while the standard bounce back scheme shows the worst agreement. The conclusion cannot be drawn that the standard bounce back scheme is necessarily incorrect, as it may be necessary to use both standard bounce back and specular reflection particles at an unsmooth interface, where a full reflection scenario is more probable. Besides, one can find from Table 1 that the relative errors of Vx are generally larger than those of Vz, that's because the variables for z component are less affected by the different rebound methods, while the variables for x component are more affected. Therefore, LBM-LSM offers flexibility in coupling the fluid-solid interface. In the following application examples, unless otherwise specified, the standard bounce back strategy is used for numerical calculations of LBM-LSM for simplicity.

3.1.4. Scholte wave simulation

The numerical simulation of seismic wavefields in fluid-solid coupled media is full of challenges. In addition to accurately reproducing the complex wavefield characteristics near the fluid-solid interface, such as the propagation law of behavior of the classical P-wave and S-wave, the generation of the Scholte wave is also a crucial yet difficult problem that cannot be overlooked. The conditions under which Scholte waves are generated have been summarized through numerical studies: Scholte wave can emerge at the fluid-solid interface when the P-wave velocity of the fluid medium lies between the P-wave velocity and the S-wave velocity of the solid medium (Glorieux et al., 2001; Carcione et al., 2018). To enhance the clarity of the Scholte wave demonstration, we employed a small relaxation time (i.e., τ = 0.65) and designed a fluid-solid interface located close to the source at a depth of 500 m, as depicted in Fig. 4. The source location, the source wavelet function, and parameters of the upperand lower-layer media remain consistent with the previous examples. For comparison, Fig. 17 illustrates the FDM-KV solution as well as snapshots of the three different LBM-LSM coupling schemes. Observations from the simulation data reveal that while the Scholte waves generated using LBM-LSM schemes exhibit a close resemblance to those derived from the conventional FDM scheme, nuanced variations are evident among the outcomes of the different LBM-LSM coupling methods. Consequently, it can be inferred that LBM-LSM schemes are efficacious for simulating Scholte waves, with the stipulation that varying fluid-solid boundary conditions exert a direct influence on the wave characteristics at the fluid-solid interface.

3.2. The seabed model

Since the media parameters are homogeneous in either the fluid or solid phase in the previous simulation examples, a more general and realistic complex geological model of the undulating seabed interface encountered in marine seismic exploration will be tested, and the specific velocity model is shown in Fig. 18. In which the depth of the water is around 450 m.

By applying a single-shot source (its horizontal distance to be 500 m) and multi-shot seismic sources (its horizontal distance to be 350, 450, 550, and 650 m) to the vertical component the vibration velocity at the depth of 637 m, the corresponding snapshots of the two components of the vibration velocity are displayed in Fig. 19. Since the sea floor is a relatively strong wave impedance interface, the reflected, transmitted, and converted waves at the sea floor can be clearly seen in Fig. 19. One can find from Fig. 19 that the simulation results of LBM-LSM closely resemble those of FDM-KV, the wave patterns near the seafloor appear more natural, and the snapshots produced by LBM-LSM exhibit almost no artifacts. The minor discrepancy in the simulation results between the two methods primarily stems from their inherent differences in solving the wave equations and the way they handle the fluid-solid boundary. Therefore, the combined LBM and LSM approach can be employed to solve complex seismic wave propagation problems in the fluid-solid two-phase coupled media.

4. Discussions

LBM is ideal for modeling waves in fluids, while LSM is suited for modeling waves in solids. Thus, handling the bidirectional propagation mechanism of waves at the fluid-solid boundary is crucial to solving wave problems in complex media. Based on our literature research, we found that the standard bounce back scheme of LBM-LSM is often used in fluid-solid coupling mechanics applications for treating the fluid-solid boundary. However, this scheme may not be suitable for simulating seismic wavefields in fluid-solid coupling media due to the propagation characteristics of seismic waves. When seismic waves cross the wave impedance interface, the reflected waves are typically reflected to the fluid region as specular reflections. Furthermore, due to the irregularity of subsurface strata, seismic waves are often diffusely reflected when crossing unsmooth wave impedance interfaces. The LBM-LSM coupling method allows for a hybrid approach of standard bounce back and specular reflection to simulate this phenomenon, with a weighting factor used to control the proportion of the two bounce back boundaries. The hybrid approach, which accounts for the actual motion of seismic waves at interfaces, is more valuable due to its enhanced accuracy and comprehensive insights into seismic activities. In addition, the proposed LBM-LSM scheme has also been employed to simulate the Scholte wave, and reasonable simulation results have been obtained. Furthermore, the fluid-solid boundary coupling method studied in this paper, through further research and development, is expected to take into account the boundary deformation, fluid transport, energy dissipation and other factors caused by the wavefield passing through the fluid-solid interface, and thus can better solve the wave propagation problem in real fluid saturating porous media. The energy decay or dissipation is directly related to the magnitude of the LBM relaxation time and the coupling strategy of the fluid-solid interface. This will provide a novel numerical simulation tool for the seismic wave propagation phenomenon in fluid-solid coupled media models (Adler and Pazdniakou, 2018; Ding et al., 2018, 2021).

However, the current fluid-solid coupling approach in this paper employs only the simplest momentum exchange algorithm. To improve simulation accuracy for practical wavefield simulations, future work may incorporate interpolation techniques into the momentum exchange algorithm. Additionally, the weighting coefficients of the proposed hybrid scheme for standard bounce back and specular reflection particles can be adjusted freely beyond the current 50% limitation. While this coupled method can solve complex seismic wavefield simulation problems in various two-phase or multi-phase media, it is still necessary to validate against physical experimental results. When dealing with the problem of seismic wave propagation at the fluid-solid interface in real geological models, it is important to choose the appropriate fluid-solid coupling treatment based on the physical characteristics of the interface. Compared to traditional forward simulation methods based on wave equations, LBM-LSM has several advantages, including independence from wave equations, flexibility in handling irregular fluid-solid boundaries, and higher simulation accuracy. This method has promising applications in seismic rock physics, reservoir fluid prediction, numerical simulation of acoustic logging, and velocity dispersion studies. It is worth noting that this paper only considers two-dimensional simulations, but the proposed method and strategy can be easily extended to three-dimensional cases.

5. Conclusions

We propose a novel method for simulating wavefields in fluid-solid coupling media using the LBM-LSM coupling approach with the standard bounce back scheme and the innovative hybrid method of standard bounce back and specular reflection fluid-solid coupling boundary conditions. This method takes into account the complex and realistic wavefield propagation effects at the fluid-solid boundary, and can be applied to seismic exploration-scale geological models with high accuracy. We validate the proposed fluid-solid coupling scheme by comparing the results of the three different boundary conditions used in seismic forward simulations for a fluid-solid two-layer model with FDM and ANA solutions. The observed differences between the proposed scheme and the two reference solutions, FDM and ANA, primarily arise from the handling of the solid-fluid coupling boundary conditions, the characteristics of the algorithms employed for wavefield tracking, and the limited elasticity capacity of the LBM algorithm. The LBM-LSM fluid-solid coupling method is then applied to simulate wave propagation in marine seismic exploration models, producing reliable results.

The LBM-LSM coupling method used in this study is flexible, and can be extended to three-dimensional cases. The numerical dispersion characteristics and stability conditions of the coupling scheme are important topics worthy of further study. The LBM-LSM coupling method can be employed in various applications including wavefield simulation of porous media, rock physics modeling, seismic data interpretation, reservoir fluid prediction, hydraulic fracturing simulation, seismic monitoring, and geotechnical engineering. With further development, it also has great potential for applications in seismic rock physics, reservoir fluid prediction, numerical simulation of acoustic logging, and velocity dispersion studies of fluid saturated porous media, etc.

CRediT authorship contribution statement

Mu-Ming Xia: Writing – original draft, Software, Methodology, Investigation. Hui Zhou: Writing – review & editing, Supervision, Methodology. Chun-Tao Jiang: Writing – original draft, Resources, Conceptualization. Han-Ming Chen: Writing – review & editing, Supervision, Formal analysis. Jin-Ming Cui: Software, Investigation. Can-Yun Wang: Supervision, Project administration. Chang-Chun Yang: Supervision, Project administration.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported in part by R & D Department of China National Petroleum Corporation (2022DQ0604-01), National Natural Science Foundation of China (42204132), the China Postdoctoral Science Foundations (2020M680667, 2021T140661), Harvard-CUP Joint Laboratory on Petroleum Science and "111" project (B13010). Mu-Ming Xia greatly acknowledges the financial support from the CAS Special Research Assistant Project.

Sidebar

ARTICLE INFO

Article history:

Received 6 May 2023

Received in revised form 4 March 2024

Accepted 16 May 2024

Available online 18 May 2024

Edited by Jie Hao and Meng-Jiao Zhou

Footnote

* Corresponding author.

E-mail address: [email protected] (H.-M. Chen).

References

References

Adler, P.M., Pazdniakou, A., 2018. Wave propagation in heterogeneous media with cavities by an LSM/LBM-coupled model. Rock Mech. Rock Eng. 51, 3907-3923. https://doi.org/10.1007/s00603-018-1602-2.

Aghamiry, H.S., Gholami, A., Operto, S., 2020, Compound regularization of fullwaveform inversion for imaging piecewise media. IEEE Trans. Geosci. Rem. Sens. 58 (2), 1192-1204. https://doi.org/10.1109/TGRS.2019.2944464.

Buxton, G.A., Verberg, R., Jasnow, D., et al., 2005. Newtonian fluid meets an elastic solid: coupling lattice Boltzmann and lattice-spring models. Phys. Rev. E 71 (5). https://doi.org/10.1103/PhysRevE. 71.056707, 056707-1-056707-16.

Carcione, J.M., Helle, H.B., 2004. The physics and simulation of wave propagation at the ocean bottom. Geophysics 69 (3), 825-839. https://doi.org/10.1190/ 1.1759469.

Carcione, J.M., Bagaini, C., Ba, J., et al., 2018. Waves at fluid-solid interfaces: explicit versus implicit formulation of the boundary condition. Geophys. J. Int. 215, 37-48. https://doi.org/10.1093/gji/ggy262.

Chen, H.M., Zhou, H., Zhang, Q.C., et al., 2017. Modeling elastic wave propagation using K-space operator-based temporal high-order staggered-grid finitedifference method. IEEE Trans. Geosci. Electron. 55 (2), 801-815. https:// doi.org/10.1109/TGRS.2016.2615330.

Chen, H.M., Zhou, H., Jiang, S.Q., et al., 2019. Fractional Laplacian viscoacoustic wave equation low-rank temporal extrapolation. IEEE Access 7, 93187-93197. https:// doi.org/10.1109/ACCESS.2019.2927760.

Chopard, B., Luthi, P., Marconi, S., 1998. A lattice Boltzmann model for wave and fracture phenomena, arXiv. Prepr. cond-mat/9812220 1-4. https://doi.org/ 10.48550/arXiv.cond-mat/9812220,

de Hoop, A.T., van der Hijden, J.H.M.T., 1983. Generation of acoustic waves by an impulsive line source in a fluid/solid configuration with a plane boundary. J. Acoust. Soc. Am. 74, 333-342, https://doi.org/10.1121/1.390970.

Ding, P.B., Di, B.R., Wang, D., et al., 2018. Pand S-wave velocity and anisotropy in saturated rocks with aligned cracks. Wave Motion 81, 1-14. https://doi.org/ 10.1016/j.wavemoti. 2018.05.001,

Ding, P.B., Gong, F., Zhang, F., et al., 2021. A physical model study of shale seismic responses and anisotropic inversion. Petrol. Sci. 18 (4), 1059-1068, https:// doi.org/10.1016/j.petsci.2021.01.001.

Escande, M., Kolluru, P.K., Cléon, L.M., et al., 2020. Lattice Boltzmann method for wave propagation in elastic solids with a regular lattice: theoretical analysis and validation. https://doi.org/10.48550/arXiv.2009.06404 arXiv preprint arXiv: 2009.06404, 1-27.

Fang. J.W., Zhou, H., Li, Y.Y.E., et al., 2023. Time-domain elastic full waveform inversion with frequency normalization, IEEE Trans. Geosci. Rem. Sens. 61 (1-17), 5903917. https://doi.org/10.1109/TGRS.2023.3246772.

Frantziskonis, G.N., 2011. Lattice Boltzmann method for multimode wave propagation in viscoelastic media and in elastic solids. Phys. Rev. E 83 (6). https:// doi.org/10.1103/PhysRevE.83.066703, 066703(1-8).

Glorieux, C., Van de Rostyne, K., Nelson, K., et al., 2001. On the character of acoustic waves at the interface between hard and soft solids and liquids. Geophys. J. Int. 110 (3), 1299-1306. https://doi.org/10.1121/1.1396333.

Graves, R.W., 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. Bull. Seismol. Soc. Am. 86 (4), 1091-1106. https://doi.org/10.1785/BSSA0860041091.

Hu, X.L., Jia, X.F., 2016. A dynamic lattice method for elastic seismic modeling in anisotropic media. Geophysics 81 (3), T131-T143. https://doi.org/10.1190/ geo2015-0511.1.

Hu, X.L., Jia, X.F., 2018. High-order dynamic lattice method for seismic simulation in anisotropic media, Geophys. J. Int. 212 (3), 1868-1889. https://doi.org/10.1093/ gji/ggx508.

Jiang, C., Zhao, G.F., 2019. A coupling model of distinct lattice spring model and lattice Boltzmann method for hydraulic fracturing. Rock Mech. Rock Eng. 52. 3675-3690. https://doi.org/10.1007/s00603-019-01819-3.

Jiang, C.T., Zhou, H., Xia, M.M., et al., 2022. Stability conditions of multiplerelaxation-time lattice Boltzmann model for seismic wavefield modeling. J. Appl. Geophys. 204 j.jappgeo.2022.104742. (1-18), 104742. https://doi.org/10.1016/j.jappgeo.2022.104742

Jiang, C.T., Zhou, H., Xia, M.M., et al., 2023. A joint absorbing boundary for the multiple-relaxation-time lattice Boltzmann method in seismic acoustic wavefield modeling. Petrol. Sci. 20 (2023), 2113-2126. https://doi.org/10.1016/ J.petsci.2023.02.019.

Ladd, A.J.C., Verberg, R., 2001. Lattice-Boltzmann simulations of particle-fluid suspensions. J. Stat. Phys. 104, 1191-1251. https://doi.org/10.1023/A: 1010414013942.

Landau, J.D., Lifshitz, E.M., 1959. Theory Elasticity, vol. 1. Pergamon, London, pp. 1-42.

Lecoulant, J., Guennou, C., Guillon, L, et al., 2022. Numerical modeling and observations of seismo-acoustic waves propagating as modes in a fluid-solid waveguide. J. Acoust. Soc. Am. 151 (5), 3437-3447. https://doi.org/10.1121/ 10.0010529.

Li, Z.C., Qu, Y.M., 2022. Research progress on seismic imaging technology. Petrol. Sci. 19 (1), 128-146. https://doi.org/10.1016/j.petsci.2022.01.015.

Li, Y.B., Shan, X.W., 2011. Lattice Boltzmann method for adiabatic acoustics. Philos. Trans. R. Soc. 369 (1944), 2371-2380. https://doi.org/10.1098/rsta.2011.0109.

Li, Q.Q., Fu, LY., Wei, W., et al., 2019. Stable and high-efficiency attenuation compensation in reverse-time migration using wavefield decomposition algorithm. IEEE Geosci. Rem. Sens. Lett. 16 (10), 1615-1619. https://doi.org/10.1109/ LGRS.2019.2903955.

Li, Q.Q., Fu, L.Y., Wu, R.S., et al., 2021. A novel wavefield-reconstruction algorithm for RTM in attenuating media. IEEE Geosci. Rem. Sens. Lett. 18 (4), 731-735. https:// doi.org/10.1109/LGRS.2020.2982038.

Li, Z., Zhao, G.F., Deng, X.F., et al., 2022. Further development of distinct lattice spring model for stability and collapse analysis of deep foundation pit excсаvation. Comput. Geotech. 144 (1-19), 104619. https://doi.org/10.1016/ j.compgeo.2021.104619.

Li, Q.Y., Wu, G.C., Jia, Z.F., et al., 2023. Full-waveform inversion in acoustic-elastic coupled media with irregular seafloor based on the generalized finitedifference method. Geophysics 88 (2), T83-T100. https://doi.org/10.1190/ GEO2022-0408.1.

Li, J.L., Qu, Y.M., Li, Y.F., et al., 2024. Ocean bottom dual-sensor Q-compensated elastic least-squares reverse time migration based on viscoacoustic and separated-viscoelastic coupled equations. Geophysics 89 (3), 1-65. https:// doi.org/10.1190/GEO2023-0433.1.

Liu, L., Liu, H., Zhang, H., et al., 2013. Full waveform inversion based on modified quasi-Newton equation. Chin. J. Geophys. 56 (7), 2447-2451. https://doi.org/ 10.6038/cig20130730 (in Chinese).

Liu, Y.Z., Huang, X.Q., Wan, X.W., et al., 2019. Elastic multi-parameter full-waveform inversion for anisotropic media. Chin. J. Geophys. 62 (5), 1809-1823. https:// doi.org/10.6038/cig2019M0176 (in Chinese).

Liu, N., Fu, LY., Tang, G., et al., 2020. Modified LSM for size-dependent wave propagation: comparison with modified couple stress theory. Acta Mech. 231, 1285-1304, https://doi.org/10.1007/s00707-019-02580-y.

Lysmer, J., Drake, LA., 1972. A finite element method for seismology. Methods Comput. Phys.: Advances in Research and Applications 11, 181-216. https:// doi.org/10.1016/8978-0-12-460811-5,50009-Χ.

Madja, G., Chin, R.C., Followill, F.E., 1985. A perturbation theory for love waves in anelastic media. Geophys. J. Int. 80 (1), 1-34. https://doi.org/10.1111/j.1365246X.1985.tb05076.x.

Mao, X.R., Zhou, Y.G., Lei, F., et al., 2021. Reverse-time migration by combining Laplacian filtering with wavefield decomposition for high-resolution subsurface imaging. IEEE J. Multiscale Multiphys. Comput. Tech. 6, 73-80. https://doi.org/ 10.1109/JMMCT.2021.3086833,

McNamara, G.R., Zanetti, G., 1988. Use of the Boltzmann equation to simulate lattice-gas automata. Phys. Rev. Lett. 61 (20), 2332-2335. https://doi.org/ 10.1103/PhysRevLett.61.2332.

Mu, X.R., Mao, Q., Huang, J.P., 2023. Stable attenuation-compensated reverse time migration and its application to land seismic data. Petrol. Sci. 20 (5), 2784-2795. https://doi.org/10.1016/j.petsci.2023.03.014.

Nguyen, T.A., Pazdniakou, A., Adler, P.M., 2021. Acoustic properties of Fontainebleau samples by lattice models. Transport Porous Media 137, 683-706. https:// doi.org/10.1007/s11242-021-01583-5.

O'Brien, G.S., Bean, C.J., 2004. A discrete numerical method for modeling volcanic earthquake source mechanisms. J. Geophys. Res. Solid Earth 109 (B9). https:// doi.org/10.1029/2004JB003023. B09301(1-7).

O'Brien, G.S., Nissen-Meyer, T., Bean, C.J., 2012. A lattice Boltzmann method for elastic wave propagation in a Poisson solid. Bull. Seismol. Soc. Am. 102 (3), 1224-1234. https://doi.org/10.1785/0120110191.

Pazdniakou, A., Adler, P.M., 2012. Lattice spring models. Transport Porous Media 93 (2), 243-262. https://doi.org/10.1007/s11242-012-9955-6.

Pazdniakou, A., Adler, P.M., 2013. Acoustic wave propagatsion in saturated porous media. In: ASME Int. Mech. Eng. Congr. Expo, vol. 2013, 56437, https://doi.org/ 10.1115/IMECE2013-62884. V014T15A001(1-11).

Qian, Y.H., Deng, Y.F., 1997. A lattice BGK model for viscoelastic media. Phys. Rev. Lett. 79 (14), 2742-2745. https://doi.org/10.1103/PhysRevLett.79.2742.

Qian, Y.H., d'Humières, D., Lallemand, P., 1992. Lattice BGK models for Navier-Stokes equation. Europhys. Lett. 17 (6), 479-484. https://doi.org/10.1209/0295-5075/ 17/6/001.

Qu, Y.M., Huang, J.P., Li, Z.C., et al., 2016. A hybrid grid method in an auxiliary coordinate system for irregular fluid-solid interface modeling. Geophys. J. Int. 208, 1540-1556, https://doi.org/10.1093/gji/ggw429.

Qu, Y.M., Guan, Z., Li, J.L., et al., 2020. Fluid-solid coupled full-waveform inversion in the curvilinear coordinates for ocean-bottom cable data. Geophysics 85 (3), R113-R133. https://doi.org/10.1190/GEO2018-0743.1.

Qu, Y.M., Li, J.L., Li, Y.F., et al., 2023. Joint acoustic and decoupled-elastic leastsquares reverse time migration for simultaneously using water-land dual-detector data. IEEE Trans. Geosci. Rem. Sens. 61 (1-11), 5909511. https://doi.org/ 10.1109/TGRS.2023.3270930.

Ren, Z.M., Dai, X., Bao, Q.Z., 2022. Source wavefield reconstruction based on an implicit staggered-grid finite-difference operator for seismic imaging. Petrol. Sci. 19 (5), 2095-2106. https://doi.org/10.1016/j.petsci.2022.05.008.

Serón, F.J., Sanz, F.J., Kindelan, M., et al., 1990. Finite-element method for elastic wave propagation. Commun. Appl. Numer. Methods 6 (5), 359-368. https:// doi.org/10.1002/cnm.1630060505.

Serón, F.J., Badal, J., Sabadell, F.J., 1996. A numerical laboratory for simulation and visualization of seismic wavefields. Geophys. Prospect. 44 (4), 603-642. https:// doi.org/10.1111/j.1365-2478.1996.tb00168.x.

Shi, Y.M., Zhang, Y., Yao, F.C., et al., 2014. Methodology of seismic imaging for hydrocarbon reservoirs based on acoustic full waveform inversion. Chin. J. Geophys. 57 (2), 607-617. https://doi.org/10.6038/cig20140224 (in Chinese).

Soares, D., 2011. Coupled numerical methods to analyze interacting acousticdynamic models by multidomain decomposition techniques. Math. Probl Eng. 2011 (1-28), 245170. https://doi.org/10.1155/2011/245170.

Sun, Y.C., Zhang. W., Ren, H.X., et al., 2021. 3D Seismic-wave modeling with a topographic fluid-solid interface at the sea bottom by the curvilinear-grid finite-difference method. Bull. Seismol. Soc. Am. 111 (5), 2753-2779, https:// doi.org/10.1785/0120200363.

Tang, J.X., Zhou, H., Jiang, C.T., et al., 2022. A perfectly matched layer technique applied to lattice spring model in seismic wavefield forward modeling for Poisson's solids. Bull. Seismol. Soc. Am. 112 (2), 608-621. https://doi.org/ 10,1785/0120210166.

Toomey, A., Bean, C.J., 2000. Numerical simulation of seismic waves using a discrete particle scheme. Geophys. J. Int. 141 (3), 595-604. https://doi.org/10.1046/ j.1365-246x.2000.00094.x.

van Vossen, R., Robertsson, J.O.A., Chapman, C.H., 2002. Finite-difference modeling of wave propagation in a fluid-solid configuration. Geophysics 67 (2), 618-624. https://doi.org/10.1190/1.1468623.

Verlet, L, 1967. Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev. 159 (1), 98-103. https:// doi.org/10.1103/PhysRev.159.98.

Virieux, J., 1986. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics 51 (4), 889-901. https://doi.org/10.1190/ 1.1442147.

Wang, Z.Y., Huang, J.P., Liu, D.J., et al., 2019. 3D variablegrid fullwaveform inver sion on GPU. Petrol. Sci. 16 (2019), 10011014. https://doi.org/10.1007/s12182019-00368-2..

Wang. N., Xing, G.C., Zhu, T.Y., et al., 2022a. Propagating seismic waves in VII attenuating media using fractional viscoelastic wave equation. J. Geophys. Res. Solid Earth 127 (4), 1-31. https://doi.org/10.1029/2021JB023280.

Wang. T.F., Cheng. J.B., Geng. J.H., 2022b. Reflection-based traveltime and waveform inversion with second-order optimization. Petrol. Sci. 19 (4), 1582-1591. https://doi.org/10.1016/j.petsci.2022.02.003.

Wang, Y.F., Bai, M., Yang, LQ., et al., 2022c. An unsplit CFS-PML scheme for the second-order wave equation with its application in fractional viscoacoustic simulation. IEEE Trans. Geosci. Rem. Sens. 60 (1-11), 5905211. https://doi.org/ 10.1109/TGRS.2021.3092714.

Wang, Y.F., Hu, J.X.Y., Harris, M., et al., 2022d. Crosswell seismic imaging using Qcompensated viscoelastic reverse time migration with explicit stabilization. IEEE Trans. Geosci. Rem. Sens. 60 (1-11), 5915611. https://doi.org/10.1109/ TGRS.2022 3176749.

Wang, L., Liu, Z.L., Rajamuni, M., 2023. Recent progress of lattice Boltzmann method and its applications in fluid-structure interaction, Proc. Inst. Mech. Eng., Part C: J. Mech. Eng. Sci. 237 (11), 2461-2484. https://doi.org/10.1177/ 09544062221077583.

Wei, H.Y., 2021. Material Failure Simulation with Random Microstructure Using Lattice Particle Method and Neural Network, vol. 6. Arizona State University. pp. 112-137.

Wilcox, LC, Stadler, G., Burstedde, C., et al., 2010. A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media. J. Comput. Phys. 229 (24), 9373-9396. https://doi.org/10.1016/j.jcp.2010.09.008.

Wu, G.C., Li, Q.Y., Wu, J.L., et al., 2018. High-order finite-difference seismic forward modeling method for fluid-solid boundary coupling media. Acta Seismol, Sin. (Chin. Ed.) 40 (1), 32-44. https://doi.org/10.11939/iass.20170010 (in Chinese).

Xia, M.M., Zhou, H., Zhang, Q.C., et al., 2015. Modeling Elastic Waves with Lattice Spring Model. SEG Technical Program Expanded Abstracts 2015. Society of Exploration Geophysicists, pp. 3734-3738. https://doi.org/10.1190/segam20155839756.1.

Xia, M.M., Wang, S.C., Zhou, H., et al., 2017a, Modelling viscoacoustic wave propagation with the lattice Boltzmann method. Sci. Rep. 7 (1-9), 10169. https:// doi.org/10.1038/s41598-017-10833-w.

Xia, M.M., Zhou, H., Li, Q.Q., et al., 2017b. A general 3D lattice spring model for modeling elastic waves. Bull. Seismol. Soc. Am. 107 (5), 2194-2212. https:// doi.org/10.1785/0120170024.

Xia, M.M., Zhou, H., Chen, H.M., et al., 2018. A rectangular-grid lattice spring model for modeling elastic waves in Poisson's solids. Geophysics 83 (2), T69-T86. https://doi.org/10.1190/geo2016-0414.1.

Xia, M.M., Zhou, H., Jiang, C.T., et al., 2022. Viscoacoustic wave simulation with the lattice Boltzmann method. Geophysics 87 (6), T403-T416. https://doi.org/ 10.1190/geo2021-0663.1.

Xu, H., Sagaut, P., 2011. Optimal low-dispersion low-dissipation LBM schemes for computational aeroacoustics. J. Comput. Phys. 230 (13), 5353-5382. https:// doi.org/10.1016/j.jcp.2011.03.040.

Yao, G., Wu, D., Wang, S.X., 2020. A review on reflection-waveform inversion. Petrol. Sci. 17 (2020), 334-351. https://doi.org/10.1007/s11430-016-9091-9.

Yu, P.F., Geng, J.H., Li, X.B., et al., 2016, Acoustic-elastic coupled equation for ocean bottom seismic data elastic reverse time migration. Geophysics 81 (5). 333-345. https://doi.org/10.1190/geo2015-0535.1.

Zhan, Q.W., Zhuang, M.W., Mao, Y.Q., et al., 2020. Unified Riemann solution for multi-physics coupling: anisotropic poroelastic/elastic/fluid interfaces. J. Comput. Phys. 402, 108961. https://doi.org/10.1016/j.jcp.2019.108961.

Zhang, J.F., 2004. Wave propagation across fluid-solid interfaces: a grid method approach. Geophys. J. Int. 159, 240252. https://doi.org/10.1111/j.1365246X.2004.02372.x.

Zhang, Q.C., Mao, W.J., Fang, J.W., 2020. Elastic full waveform inversion with sourceindependent crosstalk-free source-encoding algorithm. IEEE Trans. Geosci. Rem. Sens. 58 (4), 2915-2927. https://doi.org/10.1109/TGRS.2019.2957829.

Zhou, H., Sato, M., Liu, H., 2005. Migration velocity analysis and prestack migration of common-transmitter GPR data. IEEE Trans. Geosci. Rem. Sens. 43 (1), 86-91. https://doi.org/10.1109/TGRS.2004.839920.

Zhou, X.Y., Chang, X., Wang, Y.B., et al., 2022. Non-artifact vector Pand S-wave separation for elastic reverse time migration. Petrol. Sci. 19 (6), 2695-2710. https://doi.org/10.1016/j.petsci.2022.08.006.

Zou, Q.. Huang. J.P., Yong, P., et al., 2020. 3D elastic waveform modeling with an optimized equivalent staggered-grid finite-difference method. Petrol. Sci. 17 (2020), 967-989. https://doi.org/10.1007/s12182-020-00477-3.

Lattice Boltzmann method The lattice Boltzmann method (LBM) is a numerical modeling technique that simulates complex physical phenomena at the macroscopic level by depicting the moving state and interaction mechanism of microscopic particles at the mesoscopic physical level. Its governing equation is the lattice Boltzmann equation shown as (Qian et al., 1992)

... (A-1)

where, ... is the particle distribution function in the i direction, ci is the corresponding discrete velocity vector, and ... is the collision term, which represents the increment of the particle distribution function along the i direction due to collisions between particles. In the LBM system, the fluid density r and fluid velocity v in the macroscopic sense can be easily obtained by a simple weighted summation using the particle distribution function ... and the discrete velocity ci.

... (A-2)

In Eq. (A-1), the most commonly used collision operator in LBM is the simple linearized Bhatnagar-Gross-Krook collision operator (short as LBM-BGK) (Qian et al., 1992):

... (A-3)

such a collision operator represents a relaxation process in which the particle distribution function approaches its equilibrium state, where r is the relaxation time and f eq i is the equilibrium distribution function, which is derived from the MaxwellBoltzmann velocity distribution function in statistical mechanics (Qian and Deng, 1997; Xia et al., 2022):

... (A-4)

The choice of the discrete velocity model is of vital importance when using LBM for numerical calculations. An LBM model with too few discrete velocities may result in some physical quantities that should be conserved not satisfying the conservation law, while a model with too many discrete velocities may result in computational waste. The term "DdQq" is commonly used in LBM to label different discrete velocity models, where "d" represents the number of spatial dimensions and "q" represents the number of discrete velocities. The discrete velocity sets for the D2Q9 and D3Q19 models are shown in Figs. A-1(a) and (b), respectively. In this paper, the classical D2Q9 model is employed to conduct 2D LBM simulations. The discrete velocity sets for D2Q9 model and its corresponding weighting factors are as follows:

... (A-5)

... (A-6)

and the lattice speed of sound for the same model is ... .

To better describe the motion of the wave propagation to the interface in the inner region, the reflection coefficient (R) and transmission coefficient (T) are defined as (Xia et al., 2022)

... (A-7)

where p1, p2 and n1; n2 are the densities and streaming velocities on both sides of the interface. The improved LBM streaming step after taking the reflection and transmission effects into account is

... (A-8)

where the subscript –i represents the opposite lattice direction of i. The collision step remains unchanged. For the wavefield simulation of a viscous absorbing medium, it is particularly important to match the relaxation time t of the LBM with the quality factor Q of the media.

Through a mathematical derivation, Xia et al. (2022) obtain

... (A-9)

where Ct is the time conversion factor, and fm is the dominant frequency. With this quantitative relationship, we can adjust the different relaxation time to suit the attenuation media with different quality factors according to the actual needs. For the inner boundary condition treatment, there are common schemes such as standard bounce back, mirror bounce back, and extrapolation, etc. When it comes to wave problems, artificial reflections from truncated boundaries seriously interfere with the effective signal. By combining the viscous absorbing boundary and the PML absorbing boundary based on LBM-BGK, Jiang et al. (2023) proposed a joint absorbing boundary that better solves the truncated boundary reflection problem. In the two-dimensional case, its form is given as

... (A-10)

where t* is the modified relaxation time, h is an attenuation function, h is the perturbation component of feq, and the auxiliary quantity

APPENDIX B

Lattice spring model

Lattice spring model (LSM) is a mesoscopic-level physical model used for studying the elasticity of solid media. It consists of a grid of mass points connected by linear and angular springs, treating a macroscopic continuous medium as a series of microscopic discrete elastomeric elements. Such a discrete model can simulate physical processes such as deformation, crack propagation, and wave propagation.

In 3D space, an elastic solid can be assumed to consist of a series of interconnected networks of springs distributed in a cubic grid. The discrete models of LSM are denoted as "DdQq", in which the number "d" after "D" represents the spatial dimension, and the number "q" after "Q" represents the number of neighboring nodes connected to the central node by springs. Figs. B-1(a) and (b) show the diagram for the D2Q8 model and the D3Q18 model, respectively. For an isotropic medium, the elastic energy contained in the elastomeric element shown in Fig. B-1(b) with node i as the central node can be expressed by the following equation

... (B-1)

where K and c are the elastic constants of the linear spring and the angular spring, respectively, ... is the angle with node i as the vertex and the bonds ji and ki together form an angle, ui is the displacement of node i, ... is the unit vector pointing from node i to neighboring node j, and xij is the vector pointing from node i to neighboring node j. The first term of Eq. (B-1) is the total contribution to the elastic energy from the central force action or linear spring action of the 18 neighboring points j on the central node i, and the second term is the total contribution to the elastic energy from the bond bending action or angular spring action on all angles with the central node i as the vertex. In the case of a square mesh, the energy density * of the elastomeric element can be expressed as

... (B-2)

where a is the grid side length. By linearizing Eq. (B-2) and substituting it into Eq. (B-1), we obtain (O'Brien and Bean, 2004)

... (B-3)

where, the signs " * " and " * " represent dot product and tensor product, respectively

Using the Taylor expansion for Eq. (B-3) and expressing the displacement in terms of strain, we obtain the following equation by omitting the higher-order terms:

... (B-4)

Comparing Eq. (B-4) with the classical expression for the elastic energy of a 3D isotropic continuum (Landau and Lifshitz, 1959) yields

... (B-5)

Therefore, LSM can be employed to simulate elastic waves in a continuous medium, and the corresponding Pand S-wave velocities of computational equation have the following form (Xia et al., 2017b):

... (B-6)

where p is the bulk density of the solid media. By taking the partial derivative of the expression for the elastic energy density * with respect to ui, the total spring force at the center node i of the lattice spring model is given by the following expression for the total spring force on the neighboring nodes:

... (B-7)

LSM needs to be used in conjunction with the Verlet algorithm to perform quasi-static deformation or dynamic wavefield numerical simulations. The Verlet algorithm appears in molecular dynamics (Verlet, 1967) and it is a time integration algorithm derived from Newton's law. First of all, we adopt the Verlet algorithm to calculate the displacement u0, velocity v0, and acceleration a0 at a given time t and a given point 0. Since the position 0 can be arbitrary chosen, we replace node index 0 by (x, z), then we have

... (B-8a)

... (B-8b)

... (B-8c)

... (B-8d)

where * denotes a time step and m is the mass of the media.

In summary, the implementation process of the Verlet algorithm can be roughly divided into four steps

1) Let each node of LSM allocate the mass m.

2) Add a constant x similar to viscosity, introduce it into the vibration evolution calculation process to improve the stability of the calculation.

3) The displacement u, velocity v, and acceleration a of each node corresponding to the initial time t need to be set in advance to apply the Verlet algorithm for numerical calculation.

4) The Verlet algorithm is used to calculate the corresponding values of each physical quantity u, v, and a at time t + *.

The general LSM template can be a regular square, so as to adapt to the discretization of different geologic models in a more flexible way. Here the D2Q8 model is employed as an example to give the elastic constants in different directions. By step-by-step mathematical derivation, the isotropic spring constants are (Xia et al., 2018)

... (B-9)

where *x, *z are the spatial sampling intervals. In general, the linear springs of the same length have the same elastic constants, so we have

... (B-10)

By the way, the truncated boundary reflection problem encountered when using LSM for wavefield simulation needs to be handled by loading the perfectly matched layer absorbing boundary was developed by Tang et al. (2022).