The unit cell region Ω , filled with fluid in the porous media, is considered. Its position x represents the “center” or “centroid”. In the fluid-filled porous media of the 3D elliptical microtubule network, the characteristic length of region Ω is denoted as H h , the micro and macro characteristic length are represented by h and H , respectively. These three lengths are governed by the following equation. h ≪ H h ≪ H (8)

Here, the micro characteristic length h signifies the scale of mineral grains or solid constituents within porous media, typically in the micrometer range or smaller. Conversely, the macro characteristic length H represents larger-scale features, often on the order of hundreds of meters, and it characterizes the wavelength scale. This hierarchical relationship, grounded in multiscale analysis principles, allows us to effectively link micro-scale properties to macro-scale behavior, forming a foundational element in our modeling of fluid flow in porous media.

The volume of the unit cell Ω is denoted as V . The pore space comprises two components: the volume V s of the solid skeleton and the volume V f of fluid in porous media (i.e., the volume of the pore space). Thus, V = V s + V f . Let ψ f represent a physical quantity associated with fluid in porous media, with its value outside of the fluid defined as 0. The volume averaging method is applied to ψ f within the entire region Ω , as illustrated below. ψ f = 1 V ∫ Ω ψ f x d V (9) where ψ f is characterized by the smoothness with respect to the centroid x of region Ω . Another variable related to the volume average is the volume-averaging eigenvalue, defined as follows: ψ ¯ f = 1 V f ∫ Ω ψ f x d V (10)

According to the above definition, there is φ = ψ f / ψ ¯ f .

According to the research conducted by Slattery (Slattery, 1967) and Whitaker (Whitaker, 1966) on the volume average theorem, considering the physical variable ψ f as an illustrative example, the volume average term including the spatial derivative is as follows: ∫ Ω ∇ ⋅ ψ f d V = ∇ ⋅ ∫ Ω ψ f d V + ∫ Σ ψ f ⋅ n d S (11)

The volume average term including time derivative is as follows: ∫ Ω ∂ t ψ f d V = ∂ t ∫ Ω ψ f d V − ∫ Σ ψ f w ⋅ n d S (12) where Σ represents the interface where solid and fluid come into contact within the region Ω. It is typically specified that the normal vector n points outward from fluid, and w denotes the velocity vector of the interface between solid and fluid, which is generally taken as zero vector.

Assuming the skeleton is initially static and satisfies the assumptions of linear elasticity and isotropy, its equation of microscale momentum conservation can be expressed as follows. ρ s ∂ 2 u ∂ t 2 = ∇ ⋅ σ (13) where ρ s denotes the mass density of solid skeleton, and u represents the displacement of solid skeleton, t refers to time and σ signifies solid stress tensor.

The microscopic constitutive equation of solid is: σ = K s e ⋅ I + G ∇ u + ∇ u T − 2 3 e ⋅ I (14) where K s stands for the bulk modulus of solid particle, G refers to the shear modulus of solid skeleton, e denotes the first strain invariant, and superscript T indicates transpose.

The previous sections present the mass conservation equation, the momentum conservation equation and the constitutive equation of fluid. The equation of fluid state is expressed in the following form: 1 K f d p d t = 1 ρ f d ρ f d t (15) where p represents fluid pressure and K f signifies the fluid bulk modulus. By combining the continuity equation and the state equation of fluid, the expression of pressure term is obtained: − p = K f ∇ ⋅ u (16)

Assuming the presence of a porous media comprising a solid skeleton and a fluid, with only one existing interface denoted as Σ , which represents the contact between the fluid and the solid skeleton. The boundary conditions employed to describe this interface include the no-slip condition and fluid-solid equilibrium condition of normal force, expressed as follows: U ˙ ⋅ n = u ˙ ⋅ n − p ⋅ n + s ⋅ n = σ ⋅ n (17) where n denotes the outward normal vector of contact interface, and s is the deviatoric stress in fluid.

Based on the volume integral theorem, the macroscopic equation can be obtained by averaging the microscopic momentum conservation equation for fluids over the entire volume: ∂ 2 ∂ t 2 ρ f U = − ∇ p + ∇ ⋅ s + p ¯ ∇ ϕ − η ϕ 2 κ U ˙ ¯ − u ˙ ¯ (18) where κ represents the permeability of porous media. Since fluid does not undergo shear deformation, the fluid stress is equivalent to the normal stress s . By volume-averaging the solid momentum conservation equation and integrating it with the fluid-solid equilibrium condition at boundaries, it is obtained as follows: ∂ 2 ∂ t 2 ρ s u = − ∇ ⋅ σ − p ¯ ∇ ϕ + η ϕ 2 κ U ˙ ¯ − u ˙ ¯ (19)

The stress σ of solid skeleton can be decomposed into two components: normal stress σ s and shear stress σ e . In terms of the normal stress component of solid skeleton, Tuncay and Corapcioglu (1996), Tuncay and Corapcioglu (1997) investigated the constitutive relationship associated with volume change (solely induced by compressive deformation), and derived the equations of solid normal stress and fluid stress: σ s = − 1 − φ p ¯ s = a 11 ∇ ⋅ u ¯ + a 12 ∇ ⋅ U ¯ (20) s = − φ p ¯ = a 21 ∇ ⋅ u ¯ + a 22 ∇ ⋅ U ¯ (21) where p s represents the applied confining pressure on porous media, and the respective equations for each constant in this expression are provided below: a 11 = − 1 − φ 2 K f φ L 1 − 1 − φ 2 K f L 1 1 − φ K s − K b (22) a 12 = − 1 − φ K f L 1 (23) a 21 = − φ K s L 2 + φ K b 1 − φ L 2 (24) a 22 = − φ 2 K s 1 − φ L 2 (25) L 1 = − 1 − φ K f φ K s − 1 − φ K s 1 − φ K s − K b (26) L 2 = − φ K s 1 − φ K f − 1 + K b 1 − φ K s (27) where K b denotes the bulk modulus of dry skeleton.

The equivalent bulk modulus of the porous unit with N c fractures, volume V , and Poisson’s ratio ν can be determined (Mavko and Nur, 1978): 1 K ′ = 1 K s 1 + 2 π 3 1 − ν 2 1 − 2 ν ∑ i = 1 N c R 1 i 2 d i V (28) where R 1 i refers to the radius of the elliptical principal axis of the i − th fracture, and d i stands for fracture extension length. In a cubic unit cell with side length l , let M × N × L represent the spatial dimensions of the fracture network, where M fractures are distributed in the X direction, N fractures in the Y direction and L fractures in the Z direction. The relationship among the main radius of fracture network model, porosity ϕ , aspect ratio of fractures a , network scale and sample side length shows as follows (Xiong et al., 2021): R 1 = M N + N L + M L φ l 2 π a (29)

Assuming that the shear stress in the porous media is solely borne by the solid skeleton, disregarding any fluid-induced shear deformation: σ e = G ∇ u + ∇ u T − 2 3 ∇ ⋅ u I (30)

Where G denotes the shear modulus of solid skeleton and I represents the unit tensor. In the context of momentum transfer between phases, it is essential to consider the fluid’s viscosity due to its contribution to energy dissipation within the system. While fluids exhibit viscosity, the mechanical shear response of porous media primarily arises from the solid matrix.

In conclusion, the constitutive relationship between solid and fluid can be summarized as follows: σ = a 11 ∇ ⋅ u ¯ + a 12 ∇ ⋅ U ¯ I + G ∇ u ¯ + ∇ u ¯ T − 2 3 ∇ ⋅ u ¯ I s = a 21 ∇ ⋅ u ¯ + a 22 ∇ ⋅ U ¯ I (31)

The final form of the wave equation is derived by combining the constitutive equation with the volume-averaged momentum conservation equation: ρ s u ¨ = G ∇ 2 u + a 11 ′ ∇ e + a 12 ∇ ξ + C U ˙ − u ˙ ρ f U ¨ = a 21 ∇ e + a 22 ∇ ξ − C U ˙ − u ˙ (32) where, C = η ϕ 2 / κ , a 11 ′ = a 11 + 4 G / 3 , G is the skeleton shear modulus.

The system of equations above represents the wave propagation control equation for porous media saturated with single fluid at low frequencies, wherein the unknowns are the displacements of solid and fluid particles denoted as u and U . These equations exhibit hyperbolic characteristics with dissipative terms due to momentum transfer.

By applying the divergence operator to both sides of the wave Eq. 32 and introducing u P = ∇ ⋅ u and U P = ∇ ⋅ U , the P-wave equations is obtained: ρ s ∂ 2 u P ∂ t 2 = a 11 ′ ∇ 2 u P + a 12 ∇ 2 U P + C ∂ U P ∂ t − ∂ u P ∂ t ρ f ∂ 2 U P ∂ t 2 = a 21 ∇ 2 u P + a 22 ∇ 2 U P − C ∂ U P ∂ t − ∂ u P ∂ t (33)

By applying the curl operator to both sides of the wave Eq. 32 and introducing u S = ∇ × u and U S = ∇ × U , the s-wave equations is derived: ρ s ∂ 2 u S ∂ t 2 = G ∇ 2 u S + C ∂ U S ∂ t − ∂ u S ∂ t ρ f ∂ 2 U S ∂ t 2 = − C ∂ U S ∂ t − ∂ u S ∂ t (34)

The propagation of elastic waves in isotropic porous media is characterized by uniformity in all directions. Let the solution of the wave equation be expressed in the following prescribed form: u P = u P 0 e i ω t − k ⋅ x U P = U P 0 e i ω t − k ⋅ x (35) where u P 0 and U P 0 refer to P-wave amplitude, and k denotes wave number, which is typically a complex quantity in nature; ω signifies circular frequency, and i represents imaginary unit. After substituting the above equations into the P-wave equation, the subsequent expression can be derived: − ω 2 ρ s 0 0 ρ f + k 2 a 11 ′ a 12 a 21 a 22 + i ω − C C C − C u P 0 U P 0 = 0 0 (36)

The system of equations will have non-zero solutions only if the determinant of the coefficient matrix is equal to zero, resulting in the dispersion equation: Z 1 X 2 + Z 2 X + Z 3 = 0 (37) Z 1 = ρ s + i C ω ρ f + i C ω + C 2 ω 2 (38) Z 2 = − a 11 ′ ρ f − a 22 ρ s − i C a 11 ′ + a 12 + a 21 + a 22 ω (39) Z 3 = a 11 ′ a 22 − a 12 a 21 (40) where X = ω 2 / k 2 . The analytical expression for P-wave velocity can be obtained by solving the above equations: V p = V p R + i V p I (41) V p R = E 2 + F 2 1 / 4 ⋅ s i g n cos α 2 ⋅ 1 2 + E 2 E 2 + F 2 (42) V P I = E 2 + F 2 1 / 4 ⋅ s i g n sin α 2 ⋅ 1 2 − E 2 E 2 + F 2 (43) tan ⁡ α = F E (44) where the coefficients E and F are E = 2 a 11 ′ ρ ¯ f + a 22 ρ ¯ s ρ ¯ s ρ ¯ f + 2 a 11 ′ + a 12 + a 21 + a 22 ρ ¯ s + ρ ¯ f C 2 ω 2 ± 2 ρ ¯ s ρ ¯ f A ± 2 ρ ¯ s + ρ ¯ f C ω B 4 ρ ¯ s 2 ρ ¯ f 2 + 4 ρ ¯ s + ρ ¯ f 2 C 2 ω 2 (45) F = 2 a 11 ′ + a 12 + a 21 + a 22 ρ ¯ s ρ ¯ f C ω − 2 a 11 ′ ρ ¯ f + a 22 ρ ¯ s ρ ¯ s + ρ ¯ f C ω ± 2 ρ ¯ s + ρ ¯ f C ω A ± 2 ρ ¯ s ρ ¯ f B 4 ρ ¯ s 2 ρ ¯ f 2 + 4 ρ ¯ s + ρ ¯ f 2 C 2 ω 2 (46) A = x 2 + y 2 1 4 ⋅ s i g n cos θ 2 ⋅ 1 2 + x 2 x 2 + y 2 (47) B = x 2 + y 2 1 4 ⋅ s i g n sin θ 2 ⋅ 1 2 − x 2 x 2 + y 2 (48) tan ⁡ θ = y x (49) x = a 11 ′ ρ ¯ f + a 22 ρ ¯ s 2 − 4 a 11 ′ a 22 − a 12 a 21 ρ ¯ f ρ ¯ s − a 11 ′ + a 12 + a 21 + a 22 2 C 2 ω 2 (50) y = 2 a 11 ′ ρ ¯ f + a 22 ρ ¯ s a 11 ′ + a 12 + a 21 + a 22 C ω − 4 a 11 ′ a 22 − a 12 a 21 ρ ¯ f + ρ ¯ s C ω (51) where s i g n x refers to the application of the sign of x , and ρ ¯ s = < ρ s > , ρ ¯ f = < ρ f > .

In general, for a given frequency ω , Eq. 37 yields two complex roots (for wave number k there are four roots, but only two possess physical significance due to the requirement of continuous decrease in wave amplitude during propagation; thus, the imaginary part of k must be greater than 0), which implies that within a porous elastic media containing two immiscible fluids, two types of P-waves will propagate.

According to the method of plane wave analysis, the plane harmonic S-wave propagating along the z -direction can be expressed as: u S = u S 0 e i ω t − k ⋅ x U S = U S 0 e i ω t − k ⋅ x (52) where u S 0 and U S 0 refer to S-wave amplitude. After substituting the above equations into the S-wave equation, the subsequent expression can be derived: − ω 2 ρ s 0 0 ρ f + k 2 G 0 0 0 + i ω C − C − C C u S 0 U S 0 = 0 0 (53)

The system of equations will have non-zero solutions only if the determinant of the coefficient matrix is equal to zero, resulting in the dispersion equation: Y 1 X + Y 2 = 0 (54) Y 1 = ρ s ρ f − i C ω ρ s + ρ f (55) Y 2 = − i C ω − ρ f G (56) where X = ω 2 / k 2 . Thus the analytical expression for S-wave velocity can be solved as: X = E ′ + F ′ i (57) E ′ = − ρ ¯ f 2 ρ ¯ s G + C 2 ω 2 ρ ¯ s + ρ ¯ f G ρ ¯ s ρ ¯ f 2 + C 2 ω 2 ρ ¯ s + ρ ¯ f 2 (58) F ′ = ρ ¯ f ρ ¯ s G − ρ ¯ f ρ ¯ s + ρ ¯ f G ρ ¯ s ρ ¯ f 2 + C 2 ω 2 ρ ¯ s + ρ ¯ f 2 ⋅ C ω (59)

In general, for a given frequency ω , Eq. 54 yields only one complex root (for wave number k there are two roots, but only one possess physical significance due to the requirement of continuous decrease in wave amplitude during propagation; thus, the imaginary part of k must be greater than 0), which implies that the single S-wave form will propagate within porous elastic media containing a single fluid.

Firstly, the wave equation in terms of displacement is transformed into the velocity-stress formulation. Subsequently, finite difference scheme is implemented on the staggered grid. Finally, by considering seismic source and boundary conditions, snapshots of wave field in porous media with non-uniform saturation distribution are computed.

By using v f = U ˙ ¯ and v s = u ˙ ¯ , the macroscopic equation for fluids and solid are rewritten as ρ f v ˙ f = ∇ ⋅ s − η ϕ 2 κ v f − v s (60) ρ f v ˙ s = ∇ ⋅ σ + η ϕ 2 κ v f − v s (61)

In the equation provided above, the brackets in ⋅ denoting the volume average is omitted in the macro-scale equations for the sake of simplicity. Thus, the wave equation in velocity-stress format is: v ˙ s = 1 ρ s ∇ ⋅ σ + C ρ s v f − v s v ˙ f = 1 ρ f ∇ ⋅ s − C ρ f v f − v s (62)

Here C = η ϕ 2 κ . By applying the time derivative to both sides of the constitutive equation, the governing equation for the fractured-porous model can be obtained: σ ˙ i j = a 11 e ˙ + a 12 ξ ˙ δ i j + 2 G ε ˙ i j − 1 3 e ˙ ⋅ δ i j s ˙ i j = a 21 e ˙ + a 22 ξ ˙ δ i j (63)

It should be noted that the stress tensor exhibits symmetry, resulting in a total of 9 unknowns, namely, σ x x , σ x y , σ x z , σ y y , σ y z , σ z z , s x x , s y y and s z z , which are all associated with stress and therefore necessitate updates using constitutive relations. In the case of fluids, it satisfies the condition of s x x = s y y = s z z , thereby reducing the number of unknowns to 7. Henceforth, the components of fluid stress will no longer be differentiated, and s is employed to represent s x x , s y y and s z z instead.

The wave equation is represented in its component form by a set of 6 equations, while the constitutive equation is expressed through 7 equations, resulting in a total of 13 equations that correspond to the number of unknowns. It is assumed that both stress and velocity are initialized as zero at t=0. In the Cartesian coordinates, set x = i Δ x , y = j Δ y , z = k Δ z and t = n Δ t , where Δ x , Δ y and Δ z denote spatial step lengths and Δ t represents the time step length. As depicted in Figure 1, the staggered grid configuration is utilized wherein v x s and v x f represent the x -direction velocity components for solid and fluid particles respectively, with analogous formulations for the y and z direction components.

The staggered grid differencing operators, denoted as D h in space and D t in time, are introduced as follows: (Sun, 2009): D h f i + 1 / 2 , j , k n = θ 1 f i + 1 , j , k n − θ 2 f i , j , k n + θ 3 f i + 2 , j , k n − θ 4 f i − 1 , j , k n Δ h (64) D t f i , j , k n + 1 / 2 = f i , j , k n + 1 − f i , j , k n Δ t (65)

Here, the notation f i , j , k n represents the function to be discretized at spatial grid points i , j , k and time grid n . θ i ( i = 1,2,3,4) are the coefficients. Δ h and Δ t represent the spatial and time grid sizes, respectively. For uniformly spaced grid, Δ x = Δ y = Δ z = Δ h , θ 1 = θ 2 = 9 / 8 , θ 3 = θ 4 = − 1 / 24 , and thus the difference format with fourth-order accuracy in space and second-order accuracy in time is obtained.

Utilizing the staggered-grid finite difference, the discretization of unknown variables in the first-order velocity-stress equation occurs at distinct grid points, leading to the subsequent finite difference equations for velocity update: D t v x s i , j + 1 / 2 , k n − 1 / 2 = 1 ρ s D x σ x x i + 1 / 2 , j + 1 / 2 , k n + D y σ x y i , j + 1 , k n + D z σ x z i , j + 1 / 2 , k + 1 / 2 n + C ρ s v x f i , j + 1 / 2 , k n − 1 / 2 − v x s i , j + 1 / 2 , k n − 1 / 2 (66) D t v y s i + 1 , j , k n − 1 / 2 = 1 ρ s D x σ x y i + 1 , j , k n + D y σ y y i + 1 / 2 , j + 1 / 2 , k n + D z σ y z i + 1 / 2 , j , k + 1 / 2 n + C ρ s v y f i + 1 / 2 , j , k n − 1 / 2 − v y s i + 1 / 2 , j , k n − 1 / 2 (67) D t v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n − 1 / 2 = 1 ρ s D x σ x z i + 1 , j + 1 / 2 , k + 1 / 2 n + D y σ y z i + 1 / 2 , j + 1 , k n + D z σ x z i + 1 / 2 , j + 1 / 2 , k + 1 n + C ρ s v z f i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n − 1 / 2 − v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n − 1 / 2 (68) D t v x f i , j + 1 / 2 , k n − 1 / 2 = 1 ρ f D x s i + 1 / 2 , j + 1 / 2 , k n − C ρ f v x f i , j + 1 / 2 , k n − 1 / 2 − v x s i , j + 1 / 2 , k n − 1 / 2 (69) D t v y f i + 1 / 2 , j , k n − 1 / 2 = 1 ρ f β 2 D y s i + 1 / 2 , j + 1 / 2 , k n − C ρ f v y f i + 1 / 2 , j , k n − 1 / 2 − v y s i + 1 / 2 , j , k n − 1 / 2 (70) D t v z f i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n − 1 / 2 = 1 ρ f D z s i + 1 / 2 , j + 1 / 2 , k + 1 n − C ρ f v z f i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n − 1 / 2 − v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n − 1 / 2 (71)

The finite difference equations for stress update are as follows: D t σ x x i + 1 / 2 , j + 1 / 2 , k n = a 11 ′ D x v x s i + 1 , j + 1 / 2 , k n + 1 / 2 + a 11 * D y v y s i + 1 / 2 , j + 1 , k n + 1 / 2 + a 11 * D z v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 + a 12 D x v x f i + 1 , j + 1 / 2 , k n + 1 / 2 + D y v y f i + 1 / 2 , j + 1 , k n + 1 / 2 + D z v z f i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 (72) D t σ y y i + 1 / 2 , j + 1 / 2 , k n = a 11 * D x v x s i + 1 , j + 1 / 2 , k n + 1 / 2 + a 11 ′ D y v y s i + 1 / 2 , j + 1 , k n + 1 / 2 + a 11 * D z v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 + a 12 D x v x f i + 1 , j + 1 / 2 , k n + 1 / 2 + D y v y f i + 1 / 2 , j + 1 , k n + 1 / 2 + D z v z f i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 (73) D t σ z z i + 1 / 2 , j + 1 / 2 , k n = a 11 * D x v x s i + 1 , j + 1 / 2 , k n + 1 / 2 + a 11 * D y v y s i + 1 / 2 , j + 1 , k n + 1 / 2 + a 11 ′ D z v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 + a 12 D x v x f i + 1 , j + 1 / 2 , k n + 1 / 2 + D y v y f i + 1 / 2 , j + 1 , k n + 1 / 2 + D z v z f i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 (74) D t σ x y i , j , k n = G D y v x s i , j + 1 / 2 , k n + 1 / 2 + D x v y s i + 1 / 2 , j + 1 , k n + 1 / 2 (75) D t σ x z i , j + 1 / 2 , k + 1 / 2 n = G D z v x s i , j + 1 / 2 , k + 1 n + 1 / 2 + D x v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 (76) D t σ y z i + 1 / 2 , j , k + 1 / 2 n = G D z v y s i + 1 / 2 , j , k + 1 n + 1 / 2 + D y v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 (77) D t s i + 1 / 2 , j + 1 / 2 , k n = a 21 D x v x s i + 1 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 + a 21 D y v y s i + 1 / 2 , j + 1 , k n + 1 / 2 + a 21 D z v z s i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 + a 22 D x v x f i + 1 , j + 1 / 2 , k n + 1 / 2 + D y v y f i + 1 / 2 , j + 1 , k n + 1 / 2 + D z v z f i + 1 / 2 , j + 1 / 2 , k + 1 / 2 n + 1 / 2 (78) where a 11 * = a 11 − 2 G / 3 , a 11 ′ = a 11 + 4 G / 3 .The notations σ x x i , j , k n , σ y y i , j , k n , σ z z i , j , k n , σ x y i , j , k n , σ y z i , j , k n , σ x z i , j , k n represent the solid stress at spatial grid points i , j , k and time grid n . s i , j , k n is the discretized fluid stress. v x s i , j , k n , v y s i , j , k n , v z s i , j , k n are discretized solid velocity. v x f i , j , k n , v y f i , j , k n , v z f i , j , k n are discretized fluid velocity.

The selection of spatial step size should be based on the dispersion curve in order to determine the propagation velocity of seismic waves at a specific frequency. By calculating different types of wavelengths given the frequency, it is necessary to have at least 2–4 spatial grids within each wavelength, thus determining the appropriate spatial step size. The time step size should be determined by starting from the stability condition of the difference scheme and selecting a relatively small value to ensure that an unstable solution does not occur.

The seismic source is applied to the principal stresses σ x x , σ y y and σ z z in three directions of the solid phase, as well as the fluid stress s . The source function is based on Gaussian curve: f t = − 2 ξ 1 − 2 ξ T 2 e − ξ T 2 (79) where ξ = F 0 2 / 0.1512 represents pulse width, T = t − t s serves as time-shift parameter and t s = 1.5 / F 0 . The Higdon absorbing boundary conditions (Higdon, 1986; Higdon, 1987) are employed as the boundary conditions. Taking the second-order Higdon absorbing boundary condition as an illustration, the wave field at the boundary satisfies the following equation: ∏ j = 1 2 cos ⁡ α j ∂ ∂ t − c ∂ ∂ x u = 0 (80) where α j denotes the incident angle of absorption, and c refers to the propagation velocity of incident wave. The above formula is capable of effectively attenuating incident waves with angle of ± α j . Following discretization on the difference grid, the differential operator form of the absorbing boundary condition is presented as follows: B E x , E t − 1 = ∏ j = 1 2 β j I − E t − 1 Δ t 1 − a I + a E x − c E x − I Δ x 1 − b I + b E t − 1 (81) where β j = cos ⁡ α j represents a positive number, and Δ x signifies the spatial step size in the x direction on the boundary grid, and I denotes the identity operator, and E x I = E x , a I = a . The parameters a and b represent the weighted average coefficients for spatial and temporal dimensions, respectively. Different values of a and b correspond to distinct differential weighting methods.

Taking the left boundary as an example, f i , j , k n represents the physical variable of the wave field at the time t = n Δ t and at the point x = i Δ x and i = 0 (i.e., the left boundary in the x direction). Upon implementation of the boundary processing, it must satisfy the condition of zero reflection for the reflected wave, thereby enabling calculation of f i , j , k n through utilization of the absorption boundary difference formula: f i , j , k n = − 1 B 9 B 1 f i + 2 , j , k n − 2 + B 2 f i + 1 , j , k n − 2 + B 3 f i , j , k n − 2 + B 4 f i + 2 , j , k n − 1 + B 5 f i + 1 , j , k n − 1 + B 6 f i , j , k n − 1 + B 7 f i + 2 , j , k n + B 8 f i + 1 , j , k n (82)