The modeling of elastic wave equations considered for porous media saturated with two immiscible fluids is derived using volume averaging theory by Tuncay and Corpcioglu (1997). They combined the macroscopic momentum balance equation with the constitutive relation and incorporated the capillary effects into the bulk modulus. Based on this framework and the corresponding mathematical stability analysis, Xiong et al. (2021) simplified the system by transforming the two saturated fluids into an equivalent fluid, considering the case in which one fluid phase is entirely wrapped by the other. Consequently, the wrapped fluid is assumed to have no contact with the boundary of the representative elementary volume (REV) selected by the volume averaging theory. In this way, only one fluid phase, which affects the whole unit, needs to be considered.

By averaging the microscopic relations over the REV, the macroscopic momentum balance equations for a porous medium with an effective fluid can be represented as follows: ρ s u ¯ ¨ s = ∇ · σ + η ϕ 2 κ u ¯ ˙ f − u ¯ ˙ s ρ f u ¯ ¨ f = ∇ · s + η ϕ 2 κ u ¯ ˙ f − u ¯ ˙ s (1) where the angle brackets and overbar denote operators of taking the volume average and intrinsic average of the corresponding microscopic quantity, respectively; the subscripts s and f respectively represent the solid phase and equivalent fluid; u and U are the average displacements of solid and fluid, respectively. The dot above the quantity indicates the time derivative of the quantity. ϕ is the total porosity which can be divided into ϕ ₁ = S ₁ ϕ and ϕ ₂ = S ₂ ϕ, where S ₁ and S ₂ are the saturations of the two fluids, with the relation of S ₁ + S ₂ = 1; η and κ are the viscosity and intrinsic permeability of the effective fluid, respectively. σ and s respectively represent the stress components of solid and fluid, which need to be determined by the constitutive relation of the solid and the effective fluid phases further.

Based on the relations of the disturbances δ(1-ϕ) and δϕ, the pressures of solid (p _s ) and fluid (p _f ) in terms of ∇ · u ¯ s and ∇ · u ¯ f , the macroscopic constitutive equations of the effective fluid model can be expressed as follows: σ = a ∼ 11 ∇ · u ¯ s + a ∼ 12 ∇ · u ¯ f I + G ∇ u ¯ + ∇ u ¯ T − 2 3 ∇ · u ¯ I s = a ∼ 21 ∇ · u ¯ s + a ∼ 22 ∇ · u ¯ f I (2) where the elastic coefficients a _ij (i, j = 1, 2) are determined with the following Equation 3: a ∼ 11 E 3 = K s E 1 1 − ϕ K 1 E 2 S 1 + K 1 K 2 + K 2 A 2 S 2 + K s K b ϕ K 1 S 2 + K 2 S 1 + E 2 , a ∼ 12 E 3 = K s K 1 E 1 S 1 ϕ E 2 + K 2 , a ∼ 21 = a ∼ 12 , a ∼ 22 E 3 = K 1 S 1 ϕ ϕ K s 2 K 2 S 1 + E 2 + K 2 E 1 E 2 S 2 . (3) where E 1 = 1 − ϕ K s − K b and E 2 = α 100 1 − S 1 − n / n − 1 − 1 1 − n / n 1 − S 1 − 2 n − 1 / n − 1 , α and n are the fitting parameters needed to calculate the capillary force (Van Genuchten, 1980). E 3 = E 1 K 1 E 2 S 1 + K 1 K 2 + K 2 E 2 S 2 + K s 2 ϕ K 1 S 2 + E 2 + K 2 S 1 .

By substituting the macroscopic constitutive Equation 2 into the macroscopic momentum balance Equation 1, we obtain the system of poro-elastic wave equations for an isotropic porous medium saturated with two fluids. This equation system can be expressed in terms of displacements of solid and fluid phases as follows: ρ s ∂ 2 u ∂ t 2 = ∇ a ′ ∼ 11 ∇ · u + a ∼ 12 ∇ · U + ∇ · G ∇ u + b ∂ U ∂ t − ∂ u ∂ t , ρ f ∂ 2 U ∂ t 2 = ∇ a ∼ 21 ∇ · u + a ∼ 22 ∇ · U − b ∂ U ∂ t − ∂ u ∂ t , (4)

The governing equations for P and S waves can be derived by taking the divergence and curl to both sides of Equation 4, respectively.

Here we introduce another wave propagation model for porous media containing two fluids. To ensure the simplicity of the equations form, we adopt the framework of Biot’s theory and then introduce the effects of two fluids on wave propagation. The wave equations established by Biot (1956) can be expressed as follows: N ∇ 2 u i + A + N ∇ e + Q ∇ ξ = ρ 11 u ¨ i + ρ 12 U ¨ i + b u ˙ i − U ˙ i , Q ∇ e + R ∇ ξ = ρ 12 u ¨ i + ρ 22 U ¨ i − b u ˙ i − U ˙ i , (5) where u ˙ i = ∂ u i / ∂ t , U ˙ i = ∂ U i / ∂ t , and i = 1, 2, 3 represents the three directions in space, respectively. Based on Gedanken experiments conducted under specific conditions, the elastic responses of the fluid phase are characterized by elastic coefficients expressed as follows (Biot and Willis, 1957; Johnson, 1986; Williams, 1992): A = 1 − ϕ 1 − ϕ − K b / K s K s + ϕ K s / K f K b 1 − ϕ − K b / K s + ϕ K s / K f − 2 μ b 3 , N = G = μ b , Q = 1 − ϕ − K b / K s ϕ K s 1 − ϕ − K b / K s + ϕ K s / K f , R = ϕ 2 K s 1 − ϕ − K b / K s + ϕ K s / K f . (6)

Here K _f denotes the bulk modulus of the fluid phase. With the help of Gassmann’s fluid substitution theory (Gassmann, 1951), Biot’s theory is capable of describing complex fluid flow patterns resulting from the existence of two fluids by considering the equivalent bulk modulus of their effective fluid. This can be expressed through the arithmetic mean and harmonic mean as follows (Wollner and Dvorkin, 2018): K ¯ f = X · K f , AR + 1 − X · K f , HR (7) with Equation 8 K f , AR = S 1 K f 1 + S 2 K f 2 , K f , HR = S 1 K f 1 + S 2 K f 2 − 1 , (8) where X can be determined by Zong et al. (2022) Equation 9: X = ϕ M dry ϕ M dry + 1 − 1 + ϕ K dry / K s + G 1 − ϕ − K dry / K s · 4 / 3 · S 1 K f 1 + S 2 K f 2 , (9)

Here, M dry = K dry + G · 4 / 3 denotes the wave modulus of the dry rock skeleton. This implies that by substituting Equation 7 into Equations 5 and 6, Biot’s equations can effectively describe the wave propagation through porous rock saturated with two fluids.

We have introduced two sets of wave propagation models above, which did not consider the influence of relatively complex fluid flow mechanisms, i.e., fluid flow occurred below the macroscopic scale, on wave propagation. To accurately depict the intricate dispersions and attenuations of seismic wave over a wide frequency range, we are motivated to incorporate the squirt flow on a microscopic scale, which simultaneously exists and interacts with the global flow described by the above-mentioned SEF and Biot’s models. Following the approach outlined by Shi et al. (2024), based on the model proposed by Gurevich et al. (2010) for evaluating the stiffening effect of the combined squirt flow, the frequency-dependent elastic modulus corresponding to multiscale fluid flow patterns can be expressed as follows: 1 K m f P , ω = 1 K h + 1 1 1 K m P − 1 K h + 1 1 K f * P , ω − 1 K s ϕ c P , 1 N m f P , ω = 1 N m P − 4 15 1 K m P − 1 K m f P , ω , (10) where K f * P , ω = 1 − 2 J 1 k a k a J 0 k a K e f , k a = 1 α − 3 i ω η K e f 1 / 2 , and ω is the angular frequency, K _h denotes the original dry bulk modulus of the rock skeleton, J ₀ and J ₁ are respectively the zero-order and first-order Bessel functions, and K _ef means the bulk modulus of the effective fluid, which can be calculated using Equation 7. More details on the derivation of Equation 10 and its required parameters can be found in Shi et al. (2024).

Therefore, by substituting K _mf and N _mf in Equation 10 into Equation 4 of SEF model and Equation 5 of Biot’s equations with the equivalent fluid K _f , we can incorporate the squirt flow mechanism into the above macroscopic wave propagation theories.

In addition to the plane wave solutions, the Biot’s model and SEF model can be universally reformulated as a first-order velocity-stress formulation as Equation 11: ∂ v ∂ t − ∂ q D ∂ t = 1 ρ ρ a − ρ f 2 ρ a ρ f ρ f ρ ∇ j − p δ i j + σ i j η ϕ 2 κ q D + ∇ i p f (11) where ⟨ρ⟩ represents the total density of porous medium, and ⟨ρ _f ⟩ = ϕρ _f corresponds to the fluid density ρ _f scaled by the porosity ϕ, ⟨ρ _a ⟩ denotes the coupling coefficient between the solid and fluid phases, which is set to zero for the SEF model. We solved this system by using a conservative staggered space-time grid discretization, which is mathematically equivalent to a finite-volume (FV) method. This approach was developed by Alkhimenkov et al. (2021) for Biot’s equations on multi-graphical processing unit (GPU) architectures, along with a semi-implicit discretization for the Darcy’s flux q ^D (Alkhimenkov et al., 2021).

The discretized rectangular time-space grids with sides parallel to the coordinate axes is assumed to be superimposed on the domain {x, y, z ≥ 0} × {t ≥ 0}, with the spatial grid sizes {Δx, Δy, Δz}> 0 and time discretization Δt > 0. The time and spatial coordinates are discretized as t m = m Δ t , x i = i Δ x , y j = j Δ y , z k = k Δ z .

In this framework, the velocity field variables of both the solid and fluid phases are located at half-integer spatial nodes and inter time nodes, and can be expressed as v x s i + 1 / 2 , j , k m , v y s i , j + 1 / 2 , k m , v z s i , j , k + 1 / 2 m . Moreover, the pressure scalar field variables are allocated to the integer spatial nodes and half-integer time nodes, they are represented as p i , j , k m + 1 / 2 and p f i , j , k m + 1 / 2 , respectively.

A similar staggered discretization is employed to the stress deviator tensor fields and Darcy’s flux, constructing a conservative staggered-grid formulation with second-order accuracy in both space and time. The one-dimensional discrete representation of this system can be written as Equation 12: v i + 1 / 2 m + 1 − v i + 1 / 2 m Δ t = ρ a ρ ρ a − ρ f 2 − p + σ x x i + 1 m + 1 / 2 − − p + σ x x i l + 1 / 2 Δ x + ρ f ρ ρ a − ρ f 2 · η ϕ 2 κ χ q D i + 1 / 2 m + 1 + 1 − χ q D i + 1 / 2 m + p f i + 1 m + 1 / 2 − p f i m + 1 / 2 Δ x , − q D i + 1 / 2 m + 1 + q D i + 1 / 2 m Δ t = ρ f ρ ρ a − ρ f 2 − p + σ x x i + 1 m + 1 / 2 − − p + σ x x i m + 1 / 2 Δ x + ρ ρ ρ a − ρ f 2 · η ϕ 2 κ χ q D i + 1 / 2 m + 1 + 1 − χ q D i + 1 / 2 m + p f i + 1 m + 1 / 2 − p f i m + 1 / 2 Δ x . (12)

Herein, the weight parameter χ is introduced to control the stability of the equation system and the convergence rate of the numerical scheme, χ = 0 and 1 correspond to a fully explicit scheme and an implicit scheme, respectively.

For numerical experiments, we consider a two-dimensional domain of 9.35 × 9.35 m, discretized into 256 × 256 grid cells. Based on the von-Neumann stability analysis by Alkhimenkov et al. (2021) for linear schemes with constant coefficients, the time step Δt is determined with the given spatial grid sizes {Δx, Δy, Δz}, following the Courant–Friedrichs–Lewy (CFL) condition. In addition, it has been pointed out that this CFL condition can be extended to a fourth-order accurate in space, second-order accurate in time numerical scheme in the case that the fourth-order approximation to the first derivative is applied (Masson et al., 2006). This extension can be further explored in the construction of coefficients.

Furthermore, considering the configurations of SEF model and its discretization, we imposed reflecting boundary conditions and implemented a Ricker wavelet as the source function, which is centered within the spatial domain and can be expressed as Equation 13: s t = 1 − 2 π 2 f c 2 t − t 0 2 e − π 2 f c 2 t − t 0 2 , (13) where f _c is the source peak frequency and t ₀ is the wavelet delay.

We adapt and extend the GPU-based finite-volume staggered-grid scheme of Alkhimenkov et al. (2021) to solve the SEF-SQ and Biot-SQ models efficiently. The solver adopts second-order spatial and temporal discretization on a staggered grid, with improved stability and significant computational speed-up over CPU-based implementations.

In this section, we present an investigation of the seismic wave propagation characteristics in the porous solid media saturated with two immiscible fluids based on the two models established above, which includes the variation of P-wave velocities, the inverse quality factors with frequency, and the time-domain seismic wavefield. The rock physics parameters of the reservoir considered are enumerated in Table 1 (Tuncay and Corapcioglu, 1996; Santos et al., 2004).

The dispersion and attenuation analysis of seismic waves play crucial roles in the field of seismic rock physics. Here we examine the relationship between the wave velocity (V _P ) and the inverse attenuation quality factor (Q ⁻¹) varies with the frequency (ω) and the water saturation (S ₁).

To investigate the impact of incorporating the squirt flow mechanism, we compare the original macroscopic models, which include the SEF model and Biot’s equations with effective fluid, with their corresponding multiscale extensions, referred to as the SEF-SQ and Biot-SQ models in the following. Figure 1 illustrates the predicted V _P and Q ⁻¹ as functions of frequency for the porous media partially saturated with water.

From the left column of Figure 1, it can be observed that compared with the original (single-scale) models, both the SEF-SQ model and the Biot-SQ model with K _f exhibit an additional increase stage in attenuation. The first peaks, present in the single- and multiscale models, occur in the same frequency region. This observation is consistent with the fact that the attenuation characteristic occurring in the lower frequency range is primarily governed by the global wave-induced fluid flow. Meanwhile, the second peak, which appears at a higher frequency (approximately 10⁷ Hz), is exclusively predicted by the multiscale models. This phenomenon is attributed to the incorporated squirt flow effects, which are introduced into the two macroscopic models through K _mf and N _mf in Equation 10. As illustrated in the right column of Figure 1, a similar pattern can be observed in the dispersion curves. The curves by the single-scale models and their corresponding multiscale models overlapped well within the lower frequency range. As frequency increases, all models predict a general increase, but only the SEF-SQ and Biot-SQ models can capture the wave dispersion characteristics related to the squirt flow at higher frequencies, as evidenced by the presence of the second peaks. From a broad perspective, the multiscale wave propagation theories provide a more comprehensive characterization of P-wave attenuation and dispersion across a wide frequency range. These models not only accurately depict the attenuation and dissipation characteristics captured by the original single-scale models, but also capture additional attenuation peaks and enhance the representation of wave dispersion in the high-frequency band, which benefits the interpretation of seismic wave observations as frequency increases.

In addition to frequency-dependent effects, the multiscale Biot-SQ model includes an additional factor K _f , allowing for considering partially saturated cases. However, it can be revealed from the comparisons between the first and second rows in Figure 2 that, while the Biot’s equations with K _f is applicable for predicting the attenuation and dispersion in the partially saturated porous medium, its predictions are consistently lower than those of the SEF and SEF-SQ models.

To further validate the combination of the Biot’s equations with K _f and to investigate the influence of saturation on wave dispersion and attenuation, we analyze the variations of V _p and Q ⁻¹ as functions of water saturation S ₁ predicted by the single-scale and multiscale models in Figure 2. As saturation increases, all models exhibit similar trends in the inverse quality factor curves: their change is smooth but then decreases suddenly as water saturation approaches 1. However, even though both the Biot-SQ and SEF-SQ models are multiscale models, the predictions of the Biot-SQ model consistently appear to be lower than those of the SEF-SQ model at low water saturation. A similar discrepancy can be observed in comparison to P-wave velocity curves. Although V _P curves from these two multiscale models show sudden increases as S _w approaches 1, distinct variations occur when water saturation is approximately less than 0.4. This can be explained by the fact that, in the Biot-SQ model, as water saturation (S1) increases from 0 to approximately 0.5, the effective fluid bulk modulus increases significantly due to the replacement of compressible gas by less compressible water in the pore space, resulting in an increase in the composite P-wave velocity. Beyond this range, capillary and squirt flow effects become more significant, which tend to stabilize or slightly reduce the velocity as saturation approaches unity. This discrepancy can be attributed to the fact that the SEF model accounts for capillary force effects on the bulk modulus of both fluids (as detailed in the above Methodology section). As a result, the SEF-SQ model demonstrates superior predictive capabilities over the Biot-SQ model, providing a more reliable theoretical framework for interpreting seismic responses, particularly in cases of relatively low-saturation conditions.

The developed integrated SEF-SQ and Biot-SQ models complement each other. While both account for multiscale WIFF effects, SEF-SQ explicitly incorporates capillary effects on bulk modulus, providing superior predictions at low saturations. Biot-SQ, on the other hand, extends the Biot framework by including squirt flow, facilitating its application within a familiar theoretical structure.

In this section, we investigate the time-domain seismic wavefields to evaluate the two established models. A conservative staggered space-time grid discretization, developed by Alkhimenkov et al. (2021) specifically designed for multi-graphical processing units (GPU) numerical applications, is employed. This numerical scheme, equivalent to a finite volume approach, focuses on effectively solving the first-order velocity-stress formulation of the framework of Biot’s equations using GPU. Since the form of the equations in the SEF model can be considered as a special case of Biot’s equations, the scheme remains applicable to obtaining the wavefields based on the required relationship between the elastic parameters of the two models as Equation 14: Q = ρ 11 a 12 ρ s + ρ 12 a 22 ρ f , R = ρ 12 a 12 ρ s + ρ 22 a 22 ρ f , N = ρ 11 ρ 12 − ρ 12 2 Q ρ 12 ρ 22 + ρ 11 ρ 12 − ρ 12 2 a 21 ρ 12 ρ f − ρ 11 ρ 12 − ρ 12 2 a 11 ρ 22 ρ s , A = 2 ρ 12 2 − ρ 11 ρ 12 Q ρ 12 ρ 22 + 2 ρ 11 ρ 12 − ρ 12 2 a 11 ρ 22 ρ s + ρ 11 ρ 12 − ρ 12 2 a 21 ρ 12 ρ f . (14)

We utilize a rectangular time-space grid with time discretization (Δt) and spatial grid (Δx, Δy, Δz) for the following analysis. The velocity field variables of solid ( v ) and fluids are situated at half-integer spatial nodes and integer time nodes, while the pressure scalar field variables of solid ( σ ) and fluids (s) are defined at integer spatial nodes and half-integer time nodes. This configuration establishes a conservative staggered grid formulation with second-order accuracy in spatial and temporal fields. The numerical experiments are conducted on a two-dimensional domain of 9.35 m × 9.35 m with 256 × 256 grid cells. The time step Δt is determined according to the Courant-Friedrichs-Lew (CFL) condition proposed by Alkhimenkov et al. (2021), and 300 time-steps are executed. The reflecting boundary conditions are employed, and the seismic source function is set as a Morlet wavelet centered at the midpoint of the spatial domain.

Figures 3, 4 depict the snapshots of the time-domain seismic wavefield in the x-z spatial domain for fluid-saturated porous medium with different water saturations, as outlined in Table 1. The components {v _x , v _y , σ _xx , σ _xy , s} are plotted for comparison. Differences across the entire numerically computed fields are calculated using the following Equation 15: E r r o r = ζ S E F − S Q − ζ B i o t − S Q max ζ B i o t − S Q . (15)

Here, the symbol ζ denotes the component under consideration for estimation, including the solid velocity fields (v _x , v _y ), along with the stress fields of the solid phase (σ _xx , σ _xy ) and that of the fluid. The operator max(·) represents taking the maximum value. The subscripts ‘SEF-SQ’ and ‘Biot-SQ’ respectively represent the simulation results from the two models.

Figure 3 illustrates a comparative analysis involving the Biot-SQ model (Figure 3a), the SEF-SQ model (Figure 3b), and their difference (Figure 3c) with water saturation (S ₁) of 0.9. It is observed that both models are capable of accurately simulating the wavefields of components {v _x , v _y , σ _xx , σ _xy , s} across the entire spatial domain for the partially saturated porous medium. As depicted in Figure 3c, the relative errors of the solid velocity fields {v _x , v _y } and the solid stress fields {σ _xx , σ _xy } are consistently maintained below 2%, while the error of the fluid stress field {s} is relatively significant. It can be observed that the relative discrepancy related to the fluid stress calculated by the two models can reach up to 40%, which may be attributed to the different effective fluid modeling approaches for the two immiscible fluids applied in the models.

Figure 4 depicts a similar comparative analysis involving the Biot-SQ model (Figure 4a), the SEF-SQ model (Figure 4b), and their difference (Figure 4c) with water saturation (S ₁) of 0.1. We can observe the similar wavefields to those in Figure 3, albeit with smaller magnitudes of the calculated fluid stress fields due to the lower water saturation. It is noted that the errors of the fluid stress field {s} between the two models remain almost the same for different water saturations, further confirming that this discrepancy essentially originates from the models rather than the parameters.