Numerical Simulation of Elastic Wave Field in Viscoelastic Two-Phasic Porous Materials Based on Constant Q Fractional-Order BISQ Model

The fractional-order differential operator describes history dependence and global correlation. In this paper, we use this trait to describe the viscoelastic characteristics of the solid skeleton of a viscoelastic two-phasic porous material. Combining Kjartansson constant Q fractional order theory with the BISQ theory, a new BISQ model is proposed to simulate elastic wave propagation in a viscoelastic two-phasic porous material. The corresponding time-domain wave propagation equations are derived, and then the elastic waves are numerically simulated in different cases. The integer-order derivatives are discretised using higher-order staggered-grid finite differences, and the fractional-order time derivatives are discretised using short-time memory central differences. Numerical simulations and analysis of the wave field characterisation in different phase boundaries, different quality factor groups, and multilayered materials containing buried bodies are carried out. The simulation results show that it is feasible to combine the constant Q fractional-order derivative theory with the BISQ theory to simulate elastic waves in viscoelastic two-phasic porous materials. The combination can better describe the viscoelastic characteristics of the viscoelastic two-phasic porous materials, which is of great significance for further understanding the propagation mechanism of elastic waves in viscoelastic two-phasic porous materials and viscoelastic two-phasic porous materials containing buried bodies. This paper provides a theoretical forward simulation for fine inversion and reconstruction of layer information and buried body structure in viscoelastic two-phasic porous materials.


Introduction
Smart material is a complex material system that integrates material and structure, intelligent processing, actuation system, control system, and sensing system. Its design and synthesis involve almost all areas of high technology disciplines. The basic material components that constitute smart materials include piezoelectric materials [1], magnetostrictive materials [2], electrostrictive materials [3], and others. The emergence of smart materials will bring human civilisation to a new level. Based on dissipative structure theory, two-phasic porous materials can intelligently dissipate elastic waves. For example, the elastic waves can be attenuated by adaptively changing physical parameters such as pore fluid viscosity, porosity, and solid skeleton viscoelasticity. Therefore, the two-phasic porous material is promising to become a kind of smart material such as magnetostrictive material and electrostrictive material.
Two-phasic porous materials are the most widely available engineering materials. The propagation of bulk and interfacial waves and the response to dynamic loading in been addressed in the literature, in physics [24,25], in mechanical fluids [26], in finance models [27][28][29], in diffusion equations [30], and others. It is important to mention that there exist many fractional-order derivatives, which are not cited above, derived from the Caputo derivatives and the Riemann-Liouville derivatives. The traditional Maxwell, Voigt, Kelvin, and Zener models cannot describe such complex mechanical properties accurately and concisely, but the fractional-order calculus can compensate for the shortcomings of these models [31]. Therefore, we choose to use the constant Q fractional-order model to describe the viscoelasticity of the skeleton of two-phasic materials.
In the past decades, the wave theory related to two-phasic porous materials has been supported by knowledge in multiple fields such as petrophysics, hydrodynamics, elastodynamics, anisotropy theory, and viscoelastic theory. Scholars have developed many theoretical models of porous materials. In general, the development of theoretical models of porous materials has gone from the initial fluid replacement materials to the solid-liquid coupled two-phase (multiphase) materials and then to the fracture porous materials; from the initial uniform isotropic two-phase materials to the lateral isotropic materials and then to the anisotropic two-phase materials; from the initial completely elastic materials to the viscoelastic materials. When elastic waves propagate in the porous materials, relative motion between the solid and liquid will occur due to the pressure gradients. As the porous fluid satisfies Darcy's law, it will regain the pressure equilibrium state, during which the relative translational and rotational friction between the solid and liquid phases will convert the elastic wave energy into heat. Based on the introduction of macroscopic elastic dynamics, Biot [32,33] used Lagrange's Equation and derived the wave equation for fluid-saturated porous applicable to the entire frequency range. In this model, the two-phase porous materials contain single uniform pores, both of which are filled with liquid. Biot introduces the dissipation function containing the relative displacement of solid and fluid phases, derives the Lagrangian equation with dissipation terms, and then derives the equation of motion. According to the principle of elastodynamics and the principle of fluid mass conservation, the pressure expression is given, and then the porous elastic Biot equation is derived. Plona [34], Kelder and Smeulders [35] observed slow P waves in artificial porous materials composed of water-saturated sintered glass and natural sandstone, respectively, confirming the correctness of the Biot theory. Biot theory can well describe the attenuation and loss of elastic waves in fluid saturated porous materials and can predict the presence of slow P waves.
When the waves propagate in two-phasic porous materials, the macroscopic Biot flow mechanism and the microscopic Squirt flow mechanism act as a coupling process to influence the propagation and attenuation of elastic waves. Dvorkin and Nur [36] started from the porous isotropic one-dimensional problem and proposed a Biot-Squirt model that treats both decay mechanisms in a unified way, taking into account the free flow of the fluid in the axial direction (reflecting Biot global flow) and the flow in the radial direction (reflecting Squirt local flow), which well explains the high dispersion and strong decay observed in the experiment. Parra [37] extended the one-dimensional isotropic BISQ model to the transversely porous isotropic case starting from the frequency-domain porous elastic wave equation. Yang and Zhang [38] extend it to the high-dimensional transversely isotropic and generally anisotropic case to the high-dimensional transverse isotropic and general anisotropic case. Many scholars have further refined and developed the basic theory and application of the BISQ model.
The BISQ model can well describe the relative motion between flow-solid phases in two-phasic porous materials, including macroscopic Biot flow and microscopic Squirt flow. However, this is only one factor leading to the attenuation of elastic waves. The viscoelasticity of the solid skeleton of the two-phasic porous material also cannot be negligible. We use fractional-order derivatives to describe the viscoelasticity of the solid skeleton and innovatively introduce them into the BISQ model to establish a new fractionalorder viscoelastic BISQ model. The two main factors leading to the attenuation of elastic waves, namely the viscoelasticity of the solid skeleton and the relative flow between solid and liquid, are considered. This will be closer to the actual state than the results of previous studies. Through the numerical simulation of elastic wave propagation in viscoelastic two-phasic porous materials, we can explore the mechanism of elastic wave propagation in viscoelastic two-phasic porous materials, which helps us analyse and study the propagation and response law and characteristics of elastic waves in viscoelastic two-phasic porous materials. This new fractional BISQ model can be applied to many fields in the future. Subsurface reservoirs can be regarded as fluid-saturated two-phasic porous materials. Seismic wave simulation is the most fundamental problem in seismic exploration. Establishing a seismic wave propagation model that more closely resembles the actual subsurface medium and obtaining wavefield solutions by numerical methods are important steps in illuminating the subsurface using seismic imaging and waveform inversion. Seafloor sediments can also be regarded as a kind of two-phasic porous medium. Based on the propagation mechanism of acoustic waves in the two-phasic porous medium of the seafloor, the accepted acoustic data can be used to invert the seafloor sediment parameters, which is a more important aspect in the field of marine geophysics. Underwater anechoic materials can also be two-phasic porous materials. Therefore, the study of elastic wave propagation mechanisms in two-phasic porous materials is also crucial for designing underwater anechoic materials or even intelligent anechoic materials.

Constant Q Fractional-Order Constitutive Relationship
In terms of the stress-strain constitutive relationship, the linear viscoelasticity of the material skeleton can be expressed as the convolution of the relaxation function with the strain on the time derivative [39], where * denotes the convolution operation, X denotes the spatial location, t represents the time variable, C ijkl (t) is a fourth-order relaxation tensor function, and σ ij (X, t) and ε kl (X, t) are the stress tensor and strain tensor, respectively. For an isotropic solid skeletal, the fourth-order relaxation tensor of 81 components degenerates into a viscoelastic material stiffness matrix, which requires only two independent relaxation tensor representations.
According to the theory of fractional-order derivatives [40], the convolution operation in Equation (4) is transformed into a fractional-order time derivative, and then Equation (4) becomes ∂t 2γ is the time fractional-order derivative of 2γ order, and we have where the order of the fractional-order derivative is defined as follows where Q p and Q s are the quality factors of P-and S-wave, respectively, and Q 12 , Q 13 can be expressed respectively as [41].
In addition, the reference modulus M p , M p⊥ , M 13 , M 12 , and M µ can be expressed by the P-wave velocities v p , v p⊥ and S-wave velocities v s at the reference frequency ω γ , respectively, as follows where ρ s is the density of the solid skeleton, C 0 11 , C 0 12 , C 0 13 , C 0 33 , C 0 44 , C 0 66 are the elements of C 0 which is the elastic stiffness matrix tensor of the dry solid skeleton. For the VTI materials, we have

Constant Q Wave Propagation Equations
In this section, Kjartansson constant Q fractional-order viscoelastic constitutive relationship is used to describe the viscoelasticity of the solid skeleton. The viscoelastic skeleton constitutive relationship, geometric equations, equations of motion, and the expression of fluid pressure from the BISQ theory are given for two-phase VTI materials. Based on these equations, we derive the two-dimensional wave propagation equations.
According to [39], the viscoelastic skeleton constitutive relationship is where τ is the total stress tensor and ε(t) is the solid skeleton strain tensor, and α(t) is the poroelastic coefficient tensor of the effective stress.
For isotropic, α(t) is The expression of fluid pressure from the BISQ theory is where is the Biot flow tensor, K f is the bulk modulus of the porous fluid, S is the Squirt tensor, U denotes the displacements of the fluid, φ is the material's porosity, I is the unit matrix.
The geometric equation is The differential equations of motion are where τ ij is the component of the total stress tensor, , ρ ai is the additional mass density of the solid-fluid coupling in the x i direction, u i and U i denote the displacements of the solid and the fluid in the x i direction and x i , i = 1, 2, 3 represent the x, y, z direction respectively, η is the fluid viscosity coefficient. Fluid impedance coefficient r ij = k ij −1 , where k ij is the element of permeability tensor k, for VTI k = diag(k 11 , k 11 , k 33 ). Combining Equations (6)-(9), the wave propagation equations in VTI viscoelastic two-phase porous materials based on the fractional-order BISQ model can be derived. In the x-z plane, the wave propagation equations consist of the following five equations According to Equation (5), the above equations can be further derived and can be written as the velocity-stress equations

Finite-Difference Numerical Solution
The discretisation form of the finite-difference algorithm is intuitively simple, the mesh division form is more flexible, the computational consumption is relatively small, and it is more suitable for the simulation of wave equations of complex models. The staggered-grid finite-difference algorithm is a more advanced meshing method. Compared with the conventional finite-difference meshing method, the computational accuracy is doubled without increasing the computational consumption. In this paper, the wave propagation equations established for viscoelastic two-phasic contain both integer-order derivatives and fractional-order derivatives. The integer-order derivatives are discretised using the staggered-grid finite-difference algorithm, while the fractional-order derivatives are discretised using the central difference method.
Fractional-order differential operators can describe complex mechanical and physical processes with history dependence and spatial full-domain correlation simply and accurately. One of the main mathematical features of fractional-order derivatives is their nonlocal nature, i.e., the current state is related to all past states, and thus the fractional-order wave equations require an enormous amount of computation and storage in numerical simulations, especially for long time histories or sizeable computational domain. In order to reduce the cost of computation and improve computational efficiency, Podlubny [42] proposed the "short-time memory algorithm", which truncates the length of the operator to achieve the goal of reducing the amount of operations and storage. Only considering the finite number of operations before the current time t and the finite time interval [t − L, t] that has a large impact on the current moment, L is called the memory length.

Discretisation of Fractional Order Time Derivatives
For fractional-order derivatives, there are various definitions, and in this paper, we use the Grünwald-Letnikov definition Define when m = 0, ϕ(2γ, 0) = 1. Using the short-time memory centre difference approximation to Equation (10). When τ → ∆t , the discretisation at t = m∆t gives

Stability and Absorption Boundary Conditions
For finite-difference numerical simulations, stability should be considered in the calculation process. The stability condition of the fractional-order partial differential equations established in this paper is obtained by Fourier stability analysis. Its expression is where ∆x and ∆z are the horizontal and vertical spatial grid steps, ∆t is the time interval, δ = 0.872 is the stability constant, and c n is the corresponding staggered-grid difference operator coefficient. In order to eliminate the false reflection of the artificial boundary, the practical and straightforward Cerjan decay boundary [43] is used, i.e., expanding N grid points outward along the artificial boundary where i is the node number within the absorber layer, and a is the attenuation coefficient. The optimal value of the attenuation coefficient can be determined by several tests.

Results and Discussion
The subsurface medium can be regarded as a fluid-saturated two-phasic porous material, and seismic waves belong to elastic waves. The discussion of the factors affecting the attenuation of elastic waves in a subsurface medium and the phenomena such as scattering and transmission when encountering buried bodies are representative. Therefore, to examine the effectiveness of the Constant Q fractional-order BISQ model proposed in this paper, we select different subsurface media models for the numerical simulation study of wave fields. The grid size of the simulation area is 300 × 300, where the spatial interval is 10 m. The source of the vibration is the Ricker wavelet, which is located at (150, 150), loading in the x-direction. The time sampling interval is set to 0.5 ms.

Single Layer Model with Different Phase Boundaries
The composition of the viscoelastic two-phasic porous materials is different, and the viscosity of the fluid contained in it may also be different. In order to analyse the influence of fluid viscosity coefficient on elastic wave propagation characteristics in detail, the single-layer models with ideal phase boundary and different fluid viscosity coefficients are designed. Some basic parameters are shown in Table 1. Numerical simulations are performed for the three models in Table 1.          Figures 1 and 2 show the presence of four types of waves in the two-phasic po material, namely, fast quasi-P1 wave, fast quasi-S1 wave, slow quasi-S2 wave, and quasi-P wave from outside to inside, which are the same as the conventional Bio BISQ models. This also proves the correctness of the proposed method in this pape analysis of the figures shows that the slow quasi-P wave in the solid phase is weaker in the liquid phase. Even when the fluid viscosity coefficient is minimal, the slow qu wave cannot be observed in the solid phase. This is because the slow quasi-P wav type of highly dissipative wave. The attenuation is faster in the solid phase because s are more rigid than liquids. With the increase of the viscosity coefficient, the fast qu wave and the fast quasi-S1 wave in the solid and liquid phase all weaken, but the quasi-P wave and the slow quasi-S2 all enhance. As the fluid viscosity coefficien creases, the coupling between the solid and liquid phases is enhanced. Therefore, bot fast quasi-P1 wave and the fast quasi-S1 wave all weaken in the source loading dire However, in the direction of vertical source loading, the slow quasi-S2 wave and the quasi-P wave all enhance due to certain total energy. Due to the enhanced couplin  Figures 1 and 2 show the presence of four types of waves in the two-phasic porous material, namely, fast quasi-P1 wave, fast quasi-S1 wave, slow quasi-S2 wave, and slow quasi-P wave from outside to inside, which are the same as the conventional Biot and BISQ models. This also proves the correctness of the proposed method in this paper. An analysis of the figures shows that the slow quasi-P wave in the solid phase is weaker than in the liquid phase. Even when the fluid viscosity coefficient is minimal, the slow quasi-P wave cannot be observed in the solid phase. This is because the slow quasi-P wave is a type of highly dissipative wave. The attenuation is faster in the solid phase because solids are more rigid than liquids. With the increase of the viscosity coefficient, the fast quasi-P wave and the fast quasi-S1 wave in the solid and liquid phase all weaken, but the slow quasi-P wave and the slow quasi-S2 all enhance. As the fluid viscosity coefficient increases, the coupling between the solid and liquid phases is enhanced. Therefore, both the fast quasi-P1 wave and the fast quasi-S1 wave all weaken in the source loading direction. However, in the direction of vertical source loading, the slow quasi-S2 wave and the slow quasi-P wave all enhance due to certain total energy. Due to the enhanced coupling between the solid and liquid phases, the slow quasi-P wave in the solid phase enhances in the presence of the enhanced slow quasi-P wave in the liquid phase. In addition, the transverse waves can be seen in the liquid phase does not mean that the transverse waves can propagate in the liquid phase, but rather because the vibrations in the solid and liquid phases are coupled together.
In Figures 3 and 4, we can see that as the viscosity coefficient increases, the amplitudes of both the fast quasi-P wave and the fast quasi-S wave in the solid phase decrease and the phase is delayed. While the fast quasi-P wave amplitude in the liquid phase diminishes, the quasi-S wave amplitude increases, and the phase is delayed. For the solid phase, as the viscosity coefficient increases, the relative motion between the solid and liquid becomes weaker, and the coupling effect between the liquid phase and the solid phase becomes stronger. Therefore, the solid phase energy becomes weaker, and the phase is delayed. As the viscosity coefficient increases, the solid-phase fast quasi-P wave and the liquid phase fast quasi-P wave couple each other, resulting in the weakening of both. At certain total energy, the residual energy in the quasi-S2 wave in the liquid phase enhances in the direction of vertical source loading. Therefore, the quasi-S wave amplitude increases. The above three models are discussed under a certain defined solid skeleton quality factor. Of course, we can also discuss them under different quality factors based on our proposed new fractional-order BISQ model.

Single Layer Models with Different Quality Factor Groups
Elastic wave generates pressure in the solid skeleton, and the pressure induces fluid flow relative to the solid skeleton in a two-phase porous material. The relative flow is known as wave-induced flow. Small particles in the fluid will slowly block the pore space. This process has a relaxation property that can lead to elastic wave attenuation. In order to analyse the influence of the attenuation on elastic wave propagation in viscoelastic two-phase VTI porous materials in detail, the single-layer models with different quality factor groups are designed. Some basic parameters are shown in Table 2. Numerical simulations are performed for the three models in Table 2. Numerical simulations are performed for the three models in Table 2. Figures 5 and  6 illustrate the effects of quality factors on the solid-phase and liquid-phase wave fields, respectively. Figures 7 and 8 illustrate the effects of the fluid viscosity coefficient on the vibration of a mass in the solid-phase and liquid-phase wave fields, respectively.        In Figures 5 and 6, the amplitudes of the waves all enhance as the quality fa increase. As the quality factors increase, the viscoelasticity of the solid skeleton beco weaker, and the attenuation of elastic waves becomes weaker. As a result, the amplit of the waves enhance in both solid and liquid phases. One of the most important feat In Figures 5 and 6, the amplitudes of the waves all enhance as the quality factors increase. As the quality factors increase, the viscoelasticity of the solid skeleton becomes weaker, and the attenuation of elastic waves becomes weaker. As a result, the amplitudes of the waves enhance in both solid and liquid phases. One of the most important features of the new model proposed in this paper is the freedom to set the viscoelasticity of the solid skeleton to perform numerical simulations of the wave field in a two-phase porous material.
In Figures 7 and 8, for quasi-p wave and quasi-s wave, the mass vibration phase is delayed as the quality factors increase. This is because as the viscoelasticity of the solid skeleton becomes weaker, the rigidity becomes stronger. Masses can reach higher peaks, and reaching higher peaks takes longer. The result is the all delay phases.
In Figure 9, the location of the buried body and the snapshots of the wave field in the double-layer containing the buried body are given. Fast quasi-P, fast quasi-S1, slow quasi-S1, and slow quasi-P waves ar transmitted, and converted at the stratigraphic interface. Multiple types of w observed simultaneously in the wavefield snapshot, and the wavefield informa In addition to this, due to the difference in the quality factors of the upper and l of the solid phase, the difference in the upper and lower amplitudes of each w clearly observed. There are two reflected waves in the upper layer of the st interface and two transmitted waves in the lower layer, which come from the st interface and the upper interface of the buried body, respectively. Complex transmission, and conversion phenomena also occur at the corners of the bu These waves superimpose with the transmitted and converted waves at the st interface, which produces a more complex wave field. The study of transmissio tering mechanisms of elastic waves in two-phasic porous materials containing b ies is important for analysing the stress distribution of the buried bodies in porous materials in the elastic wave field. For example, it can be used to analy of underground buildings in seismic wave fields.

Conclusions
In this paper, the constant Q fractional-order theory is combined with th ory to establish a more generalised constant Q fractional-order wave model Fast quasi-P, fast quasi-S1, slow quasi-S1, and slow quasi-P waves are reflected, transmitted, and converted at the stratigraphic interface. Multiple types of waves can be observed simultaneously in the wavefield snapshot, and the wavefield information is rich. In addition to this, due to the difference in the quality factors of the upper and lower layers of the solid phase, the difference in the upper and lower amplitudes of each wave can be clearly observed. There are two reflected waves in the upper layer of the stratigraphic interface and two transmitted waves in the lower layer, which come from the stratigraphic interface and the upper interface of the buried body, respectively. Complex reflection, transmission, and conversion phenomena also occur at the corners of the buried body. These waves superimpose with the transmitted and converted waves at the stratigraphic interface, which produces a more complex wave field. The study of transmission and scattering mechanisms of elastic waves in two-phasic porous materials containing buried bodies is important for analysing the stress distribution of the buried bodies in two-phasic porous materials in the elastic wave field. For example, it can be used to analyse the force of underground buildings in seismic wave fields.

Conclusions
In this paper, the constant Q fractional-order theory is combined with the BISQ theory to establish a more generalised constant Q fractional-order wave model to describe the propagation of elastic waves in the viscoelastic two-phasic porous materials. The twodimensional velocity-stress wave propagation equation is derived. Fractional-order time derivatives are truncated by the short-time memory method. Numerical simulations of wave fields in the material for different cases are performed. By analysing the simulation results, we can draw several conclusions.
As the fluid viscosity coefficient increases, the pore fluid becomes more viscous. The fluid flow between the pores and fractures becomes slower, hindering the P-wave propagation. The fast quasi-P and quasi-S1 waves in the solid phase and the fast quasi-P wave energy in the liquid phase weaken. Under the condition of constant total energy and weak dissipation of relative motion, the residual energy of the slow quasi-S2 wave in the liquid phase is enhanced. As the quality factors increase, the viscoelasticity of the solid skeleton of the viscoelastic porous material decreases and the rigidity becomes stronger. The amplitudes of the waves all enhance. The time required to reach higher peaks becomes longer, resulting in the phase delay.
In a double-layer material containing buried body, complex scattering, transmission, and waveform conversion phenomena occur at the stratigraphic interface, buried body interface, and corners. In-depth study of wave scattering, transmission, and conversion phenomena in viscoelastic two-phasic porous materials or complex buried structures and analysis of various wave propagation characteristics can deepen the understanding of wave propagation mechanism in actual viscoelastic two-phasic porous materials. This may help conduct more in-depth research on viscoelastic two-phasic porous materials using multiwave elastic data, with important theoretical and practical applications.
This paper also has several shortcomings. Physical experiments should be carried out to further verify the accuracy of the model developed. The BISQ model assumes that the pore fluid flow satisfies Darcy's law. In reality, the fluid may not satisfy Darcy's law due to the semiconnectivity of the pores or the presence of particles in the fluid that can block the pores. Therefore, a non-Darcy's law BISQ model should be considered.