A Pressure- and Frequency-Dependent Multiscale Model for Describing the Wave Propagation Characteristics of Fluid-Saturated Porous Media

Abstract

The mechanism of wave propagation in fluid-saturated porous media is influenced by pressure and frequency. Pressure dependence is mainly dominated by the opening and closing of compliant and stiff pores in rocks, as well as nonlinear deformation respect to high-order elastic constants. Frequency dependence is mainly reflected in the dispersion and attenuation caused by wave-induced fluid flow (WIFF). Therefore, the propagation characteristics of seismic waves in subsurface rocks when pressure and frequency are coupled have broad practical significance, such as geofluid discrimination and in situ stress detection. A new equivalent elastic modulus applicable to fluid-saturated porous media has been established, which simultaneously considers the effects of pressure and WIFF. First, the dual-porosity model is incorporated to account for the changes in rock porosity under pressure and corresponding linear and nonlinear deformations. Then, based on the heterogeneity of rock at the mesocale and microscale, a unified pressure- and frequency-dependent elastic modulus over a wide frequency band is established using the Zener model. The wave equation of fluid-saturated porous media is constructed using the new model, and the pressure- and frequency-dependent phase velocities are derived. Rock physics and digital simulation experiments are applied to analyze the variation of elastic parameters and velocity with pressure and frequency. Comparison with experimental measurement data shows that the new model has higher accuracy than traditional models, especially in the low effective pressure region and the frequency band respect to seismic exploration.

1. Introduction

Unconventional oil and gas have become a hotspot in the field of oil and gas exploration and development, attracting attention to the knowledge of subsurface pressure fields. The internal pore structure and shape of subsurface rocks are closely related to pore pressure. Under effective pressure, linear or nonlinear deformation occurs (depending on whether they are stiff or compliant pores) due to compression, which greatly affects the elastic modulus and elastic wave velocity of the rock. The traditional linear elastic theory, widely used in elastic wave dynamics, lacks an accurate description of nonlinear deformation of compliant pores caused by pressure, so there may be deviations in predicting the elastic modulus or velocity. Furthermore, when reservoir fluids are considered, the dispersion and attenuation characteristics of seismic waves propagation in subsurface media cannot be ignored. Due to the heterogeneity of the medium, the relative motion between pore fluid and rock skeleton during wave propagation is called the wave-induced flow mechanism (WIFF). The energy dissipation caused by WIFF leads to the phenomenon of seismic wave dispersion and attenuation. In general, WIFF is coupled with the effective pressure, and the equivalent elastic modulus and velocity of rock are the results of their joint influence. Therefore, clarifying the seismic response mechanism of fluid-saturated porous media under the influence of pressure and frequency is of great significance for the precise exploration and efficient development of reservoirs.

Pressure-induced rock elastic modulus and velocity variations have been studied for a long time. Based on finite deformation, the acoustoelastic theory that describes the influence of stress on rock mechanical properties is more in line with the deformation characteristics of deep rocks and therefore has received widespread attention. The acoustoelastic theory suggests that the rock’s deformation is the result of static strain caused by initial stress combined with small perturbations of seismic waves [1, 2]. On this basis, the stiffness matrix of the rock not only depends on its second-order elastic constant but also on its third-order elastic constant and initial stress–strain state, which is called the equivalent stiffness matrix [3-7]. Kim and Hong [8] studied the relation between stress and seismic reflection coefficient by analyzing the anisotropy vary with uniaxial stress based on Ruger approximation. Pistre et al. [9] proposed the concept of elastic constant difference constant and established the relation between principal stress difference and elastic constant difference. Ba et al. [10] derived the relation equations between stress, velocity, and frequency of rocks under hydrostatic stress and uniaxial stress conditions. Liu et al. [11] conducted acoustoelastic analysis on saturated porous rocks, and based on this, obtained rock elastic modulus and velocities that are more consistent with experimental measurement results. However, when large magnitude of pressure is applied, there is a significant deviation between the predicted results of traditional acoustoelastic theory and the measured values, mainly caused by the pressure-induced deformation of rock pores and cracks. Based on the dual-porosity model, an improved porous acoustoelastic model is proposed, which can account for the linear or nonlinear deformation of the internal pores of rocks [12-14].

Some researches show that the equivalent medium model of acoustoelastic theory in dry rock is in good agreement with the real measurement results, but there are still obvious deviations in saturated rock [15, 16]. In order to study the elastic properties of saturated rocks, the effect of WIFF is considered, which is more suitable for the real situation of subsurface media. Biot [17] pioneered the establishment of an elastic wave theoretical model for fluid-saturated porous media and started the research on the fluid dispersion and attenuation effect caused by WIFF. Mavko and Jizba [18] proposed the squirt-flow theory, which improved the underestimation of high-frequency attenuation by Biot’s macroscopic fluid flow. And the wave-induced pressure gradients from a microscopic perspective is considered [19]. It is worth noting that the obvious frequency dispersion and attenuation phenomenon in these studies mostly occurs in the ultrasonic frequency band, lacking analysis of WIFF effects in the seismic frequency band. The recent laboratory measurements have showed that the squirt flow can also apparent in seismic frequency band, when the viscosity of fluid in rocks is high [20, 21]. But, another dispersion and attenuation mechanism caused by the heterogeneity of rocks, which often occurs in frequency band, is not involved in these theory yet [22]. Pride et al. [23] proposed a mesoscopic DPDP model, which predicted the WIFF in the seismic frequency band based on the heterogeneity characteristics of rocks or pore fluids and establishes the relationship between seismic data and dispersion and attenuation characteristics of rock. Then, Zhang et al. [24] used Zener model to couple this mesoscopic phenomenon with the microscopic one and proposed a unified bulk modulus that can describe the dispersion and attenuation due to heterogeneity of rocks. Zhao et al. [25] used dynamic volumetric strain to analyzed WIFF of heterogeneous rocks and successfully capture the dispersion and attenuation characteristics for porous rock. Li and Yan [26] also considered this heterogeneity through crack reservoirs and derived a multiscale elastic modulus to characterize the WIFF in porous media. However, in the classic WIFF model mentioned above, the influence of pressure on rock mechanical properties or pore’s deformation was not considered. Gurevich et al. [27] established a frequency- and pressure-dependent rock elastic model based on the pressure relaxation in the intergranular contact area. The model is consistent with Gassmann model and Mavko-Jizba equation in both high- and low-frequency bands. Then, some generalized and simplified equations based on this model are proposed to adapt to more practical situations [28-30]. Chen et al. [31] obtained an updated Shapiro’s model considering the nonlinear change of stiff pores under pressure and combined with Gurevich’s theory to establish a more accurate pressure- and frequency-dependent elastic modulus equation. Moreover, the effect of crack closure, squirt flow, and acoustoelastic are considered together to establish a pressure-dependent model [32]. Then the differential effective medium and dual-porosity theory are used to derive a broadband model [33, 34].

At present, the pressure- and frequency-coupled elastic modulus equation is mainly based on the dual-porosity model, which is a semiempirical relationship and lacks clear physical significance. Moreover, the frequency variation characteristics in the seismic frequency band are usually ignored, which limits its application in practical production. In this study, by analyzing the dispersion and attenuation mechanism of fluid-saturated porous media at mesoscale and microscale simultaneously, a new model which can accurately describe the frequency variation characteristics of rocks in the full frequency band is proposed. On this basis, combined with the pressure-dependent dual-porosity model, a unified pressure- and frequency-coupled elastic modulus equation is established. According to the poro-acoustoelastic theory, the energy change caused by nonlinear deformation of rock (including the nonlinear deformation of compliant pores under compression) is analyzed. The new unified elastic modulus is substituted into it, and the pressure- and frequency-depended wave equation is constructed. As a result, the mechanism of wave propagation in fluid-saturated porous media can be fully understood by solving the complex velocities through wave equation. The feasibility of these new models and velocities is verified by comparing with the measured experimental data.

2. Methodology

The elastic modulus is an important parameter that represents the physical properties of rocks. It serves as a bridge between the microstructure of rocks such as mineral composition and distribution and macroscopic elastic properties such as rock velocity. Thus, correct characterization of elastic modulus plays an important role in the construction of rock constitutive equations and the solution of wave equations. For porous media, except physical properties of the constituent minerals, there are two aspects that cannot be ignored: pressure effect and frequency dependence.

2.1 Pressure Effect

Shapiro [12] divided the total porosity $φ$ of porous rocks into two parts, the stiff porosity part $φ_{s}$ and the compliant porosity part $φ_{c}$ ⁠:

φ = φ_{c} + φ_{s} .

(1)

Then, the aspect ratio $γ$ was proposed as a parameter for evaluating the relation of minimum and maximum dimensions of a pore; it is a dimensionless parameter closely related to pressure. It is generally believed that excessive pressure causes the closure of internal compliant pores, and the aspect ratio of rock is mainly supported by stiff pores (supported by more-or-less isometric pores), which is greater than 0.1. Otherwise, the aspect ratio of rock is dominated by compliant pores (supported by thin cracks and grain contact vicinities), which is less than 0.01. The pressure-dependent aspect ratio $γ (P)$ be represented as [13]:

\begin{matrix} γ (P) = γ_{i} - γ_{m i n} \\ γ_{m i n} = \frac{4 (1 - o_{d r s}^{2}) P}{π E_{d r s}} \end{matrix},

(2)

where P represents effective pressure, $γ_{i}, γ_{m i n}$ represent initial aspect ratio and minimum aspect ratio of the unclosed compliant pores at effective pressure P, $o_{d r s}, E_{d r s}$ represent Poisson’s ratio and Young’s modulus of rock skeleton when $φ_{c} \to 0$ ⁠, respectively. The aspect ratio decreases with increasing effective pressure.

Under dual-porosity theory, the compressibility of rock $C_{d r}$ can be provided by Taylor expansion:

C_{d r} (φ_{s}, φ_{c}) = C_{d r s} [1 + θ_{s} φ_{s} + θ_{c} φ_{c}],

(3)

where $θ_{s}, θ_{c}$ are elastic bulk strain of $φ_{s}$ and $φ_{c}$ ⁠, respectively. $C_{d r s}$ is compressibility of a hypothetical medium without compliant pores. The relation between rock compressibility and bulk modulus is:

K_{d r} = 1 / C_{d r}, K_{d r s} = 1 / C_{d r s} .

(4)

In equation (3), $φ_{s}, φ_{c}$ are both pressure-dependent and can be shown as

{\begin{cases} \begin{matrix} φ_{s} = P (C_{g} - C_{d r s}) \\ φ_{c} = φ_{c 0} e^{(- θ_{c} P C_{d r s})} \end{matrix} \end{cases}

(5)

substituting equation (5) into equations (3) and (4) yields

C_{d r} (P) = C_{d r s} [1 - θ_{s} (C_{d r s} - C_{g}) P + θ_{c} φ_{c 0} e^{(- θ_{c} P C_{d r s})}],

(6)

and

K_{d r} (P) = K_{d r s} [1 + θ_{s} (C_{d r s} - C_{g}) P - θ_{c} φ_{c 0} e^{(- θ_{c} P C_{d r s})}] .

(7)

The quantities θ can be approximately estimated using the upper Hashin–Strikma bound:

{\begin{cases} θ_{s} \approx 1 + \frac{3 K_{g}}{4 μ_{g}} \\ θ_{c} \approx \frac{K_{g} (3 K_{g} + 4 μ_{g})}{π γ_{m i n} μ_{g} (3 K_{g} + μ_{g})} \end{cases},

(8)

where $C_{g}$ represent compressibility of rock matrix, $K_{g}, μ_{g}$ represent bulk and shear modulus of rock matrix, $φ_{c 0}$ represent initial compliant porosity in the case of P $\to$ 0.

From equation (7), we get the equivalent bulk modulus of porous rocks under pressure (⁠ $K_{d r s}$ is bulk modulus when $φ_{c} \to 0$ ⁠). Similar to equation (7), the equivalent shear modulus of porous rocks can be expressed as

μ_{d r} (P) = μ_{d r s} [1 + θ_{s μ} (C_{d r s} - C_{g}) P - θ_{c μ} φ_{c 0} e^{(- θ_{c} P C_{d r s})}],

(9)

where

{\begin{cases} θ_{s μ} \approx 1 + \frac{6 K_{g} + 2 μ_{g}}{9 K_{g} + 8 μ_{g}} \\ θ_{c μ} \approx \frac{1}{5} [1 + \frac{4 (3 K_{g} + 4 μ_{g}) (9 K_{g} + 4 μ_{g})}{3 π γ_{m i n} (3 K_{g} + μ_{g}) (3 K_{g} + 2 μ_{g})}] \end{cases}

(10)

$θ_{s μ}, θ_{c μ}$ represent elastic shear strain of $φ_{s}$ and $φ_{c}$ ⁠.

The pressure-dependent bulk and shear modulus of porous media obtained from equations (7) and (9) provide a basis for studying the pressure effects on the elastic modulus and velocities of rocks, and details of the derivation process are provided in Shapiro’s work [12].

2.2 Frequency Dependence

When elastic waves propagate in fluid-saturated porous media, a pressure gradient is formed within the pores, leading to relative motion between the pore fluid and the rock skeleton. This relative motion is reflected as dispersion and attenuation of bulk modulus, which is happened in macroscale (Biot model), mesoscale (DPDP model), and microscale (squirt-flow theory). Establishing a unified bulk modulus model $K_{u n i f i e d}$ across full frequency bands is beneficial for more accurate characterization of the frequency dependence of rocks.

The dispersion characteristics of saturated medium can be divided into two regimes and a transition. In low-frequency regime, fluid and skeleton move freely, and the equivalent bulk modulus at this time can be regarded as bulk modulus of dry rock, known as the relaxed modulus. In high-frequency regime, the fluid pressure between the stiff and compliant pores does not have enough time to reach equilibrium, and the fluids in stiff pores are isolated and can be considered as part of the rock matrix. At this point, the equivalent bulk modulus is equal to the unrelaxed modulus when $φ_{c}$ = 0. There is a transition between these two regimes, in which the rock bulk modulus changes rapidly, and the frequency corresponding to this change is called the characteristic frequency. Overall, the change of bulk modulus in the full frequency range $K_{u n i f i e d}$ can be described by Zener model:

K_{u n i f i e d} (ω) = K_{l} \frac{(1 + i \frac{ω}{ω_{c}} \sqrt{\frac{K_{h}}{K_{l}}})}{(1 + i \frac{ω}{ω_{c}} \sqrt{\frac{K_{l}}{K_{h}}})},

(11)

where $K_{l}, K_{h}, ω, a n d ω_{c}$ represent relaxed modulus, unrelaxed modulus, frequency, and characteristic frequency. In this equation, when $ω ≪ ω_{c}$ ⁠, $K_{u n i f i e d}$ → $K_{l}$ ⁠,when $ω ≫ ω_{c}$ ⁠, $K_{u n i f i e d}$ → $K_{h}$ ⁠.

When the mesoscopic heterogeneity is considered (patchy fluids, cracks, and heterogeneity elastic modulus), it is generally believed that there will be a significant dispersion and attenuation within the seismic frequency band, so another modified Zener model is proposed [24]. When $ω_{c}^{m i}$ is microscopic characteristic frequency, $ω_{c}^{m e}$ is mesoscopic characteristic frequency, $K_{h}^{m i}$ and $K_{h}^{m e}$ are corresponding microscopic and mesoscopic unrelaxed modulus of rocks, the modified Zener model is

K_{u n i f i e d} (ω) = K_{l} \frac{(1 + i \frac{ω}{ω_{c}^{m e}} \sqrt{\frac{K_{h}^{m e}}{K_{l}}})}{(1 + i \frac{ω}{ω_{c}^{m e}} \sqrt{\frac{K_{l}}{K_{h}^{m e}}})} \cdot \frac{(1 + i \frac{ω}{ω_{c}^{m i}} \sqrt{\frac{K_{h}^{m i}}{K_{h}^{m e}}})}{(1 + i \frac{ω}{ω_{c}^{m i}} \sqrt{\frac{K_{h}^{m e}}{K_{h}^{m i}}})}

(12)

These characteristic frequencies can be expressed as [35, 36]:

{\begin{matrix} ω_{c}^{m e} = 2 π \frac{M κ}{η L^{2}} \\ ω_{c}^{m i} = 2 π \frac{γ^{3} K_{g}}{η L^{2}} \end{matrix},

(13)

where $κ, η, a n d L$ are permeability, fluid viscosity, and the scale of the mesoscopic heterogeneity, respectively. And the fluid storage coefficient M is equal to

M = {[\frac{φ}{K_{f}} + \frac{α - φ}{K_{l}}]}^{- 1,}

(14)

where $K_{f}$ is bulk modulus of fluid, $α$ is Biot coefficient.

Besides, the unrelaxed moduli of microscopic and mesoscopic are as follows:

{\begin{cases} K_{h}^{m e} = K_{l} \frac{a_{33}}{a_{33} - a_{13}^{2} K_{l}} \\ K_{h}^{m i} = K_{d r s} \end{cases},

(15)

where $a_{13} a n d a_{33}$ are both relative coefficients in DPDP model, their specific form can be found in online supplementary Appendix 1.

From equations (12), (13), and (15), the change of the bulk modulus of fluid-saturated porous media in the full frequency range can be divided into three frequency-independent stages, corresponding to $K_{l}, K_{h}^{m e}, K_{h}^{m i}$ ⁠, and two frequency-dependent transform regions when frequency is close to the characteristic frequency $ω_{c}^{m e}, ω_{c}^{m i}$ ⁠, as shown in Figure 1. When $ω ≪ ω_{c}^{m e}$ ⁠, equivalent bulk modulus is equal to the bulk modulus of dry rock. When $ω ≫ ω_{c}^{m i}$ ⁠, it is equal to the bulk modulus when the compliant pores are closed. When $ω_{c}^{m e} ≪ ω ≪ ω_{c}^{m i}$ ⁠, equivalent bulk modulus is a constant, corresponding to the intermediate state in which the fluid pressure in mesoscale cannot be balanced, but the microfluid can still reach equilibrium.

According to the pressure effect in Section 2.1, the bulk modulus of dry rock is a function of pressure:

K_{l} = K_{d r} (P) .

(16)

Then we have obtained a unified bulk modulus of rocks, which is related to pressure and frequency:

K_{u n i f i e d} (P, ω) = K_{d r} (P) \cdot \frac{(1 + i \frac{ω}{ω_{c}^{m e}} \sqrt{\frac{K_{h}^{m e}}{K_{d r} (P)}})}{(1 + i \frac{ω}{ω_{c}^{m e}} \sqrt{\frac{K_{d r} (P)}{K_{h}^{m e}}})} \cdot \frac{(1 + i \frac{ω}{ω_{c}^{m i}} \sqrt{\frac{K_{h}^{m i}}{K_{h}^{m e}}})}{(1 + i \frac{ω}{ω_{c}^{m i}} \sqrt{\frac{K_{h}^{m e}}{K_{h}^{m i}}})}

(17)

Hence, according to equation (20) in Gurevich’s research [27] and equation (17), the equivalent shear modulus of rocks can be expressed as

\begin{aligned} μ_{u n i f i e d} {(P, ω)}^{- 1} = μ_{d r} {(P)}^{- 1} - 4 \\ [K_{d r} {(P)}^{- 1} - K_{u n i f i e d} {(P, ω)}^{- 1}] / 15 \end{aligned}

(18)

From equations (17) and (18), we derived a new unified modulus which is a function of pressure and frequency. Different from some existing works, this new model simultaneously considers both the effect of mesoscopic flow (caused by the heterogeneity of fluid) and squirt flow. Thus, it can be applied to a wider frequency band to characterize the pressure effect and dispersion within the seismic frequency band. Then this modulus is applied to wave equation to calculate phase velocity of fluid-saturated porous media later.

2.3 Velocities

Fluid-saturated porous rock exhibits significant nonlinear characteristics when subjected to a certain degree of pressure. At this time, the solution of dynamic equation under linear elasticity theory is not suitable. Combining the acoustoelastic theory, a more accurate equation of fluid-saturated porous media is established. The acoustoelastic theory divided the deformed rock into three configurations, the unstressed configuration (or natural configuration), initial configuration, and final configuration. When stress applied to rocks, they transform from unstressed configuration to initial configuration, and consider the elastic waves, rocks transform from initial configuration to final configuration. Considering that the state of rock in unstressed configuration is often unknown, it is more reasonable to use the initial configuration to describe the motion of rock.

Based on Biot’s theory, the energy equation for saturated porous rocks is [37]:

\begin{matrix} \frac{d}{d t} (\frac{\partial T}{\partial u_{i, t}}) + \frac{\partial D}{\partial u_{i, t}} = \frac{d}{d X_{j}} \frac{\partial W}{\partial u_{i, j}}, (i = 1, 2, 3) \\ \frac{d}{d t} (\frac{\partial T}{\partial w_{i, t}}) + \frac{\partial D}{\partial w_{i, t}} = \frac{d}{d X_{j}} \frac{\partial W}{\partial w_{i, j}}, (i = 1, 2, 3), \end{matrix}

(19)

where W represents the potential of rock, T represents kinetic energy, D represents viscous dissipation, $u_{i}$ represents the displacements of solid skeleton, $w_{i}$ represents the fluid flow relative to the solid phase, $w_{i} = φ (U_{i} - u_{i})$ ⁠, $U_{i}$ represents the displacements of fluid and $X_{j}$ represents Lagrangian initial coordinates.

T and D can be shown as [17]:

\begin{matrix} 2 T = ρ u_{i, t}^{2} + 2 ρ_{f} u_{i, t} w_{i, t} + ρ_{c} w_{i, t}^{2} \\ 2 D = \frac{η}{κ} w_{i, t}^{2,} \end{matrix}

(20)

where $ρ = (1 - φ) ρ_{s}$ + $φ ρ_{f}$ , $ρ_{c} = τ ρ_{f} / φ$ ⁠. $ρ_{s}, ρ_{f}$ represent density of solid and fluid respectively and $τ$ is tortuosity of the porous space.

The Lagrangian strains of solids $E$ and fluids $ζ$ relative to solids are expressed as [10]:

\begin{matrix} E_{i j} = \frac{1}{2} (\frac{\partial u_{i}}{\partial X_{j}} + \frac{\partial u_{j}}{\partial X_{i}} + \frac{\partial u_{k}}{\partial X_{i}} \frac{\partial u_{k}}{\partial X_{j}}) \\ ζ = - (\frac{\partial w_{i}}{\partial X_{i}} + \frac{1}{2} \frac{\partial w_{j}}{\partial X_{j}} \frac{\partial w_{k}}{\partial X_{k}}) \end{matrix}

(21)

Substituting equations (20) and (21) into equation (19), yields:

\begin{matrix} \frac{\partial}{\partial t^{2}} (ρ u_{i} + ρ_{f} w_{i}) = \frac{\partial}{\partial X_{j}} [(δ_{k i} + \frac{\partial u_{i}}{\partial X_{k}}) \frac{\partial W}{\partial E_{k j}}] \\ \frac{\partial}{\partial t^{2}} [ρ_{f} u_{i} + ρ_{+} w_{i}] = - \frac{\partial}{\partial X_{i}} [(1 + \frac{\partial w_{k}}{\partial X_{k}}) \frac{\partial W}{\partial ζ}], \end{matrix}

(22)

Where

ρ_{+} = ρ_{c} + η / κ

According to Grinfeld and Norris’s work, the energy of isotropic saturated rocks can be expressed as [38]:

\begin{matrix} W = \frac{1}{2} (\bar{K} - \frac{2}{3} \bar{u} + α^{2} M) E^{2} + \bar{u} E_{i j}^{2} - α M E ζ + \frac{1}{2} M ζ^{2} \\ + \frac{1}{6} ψ_{1} E^{3} + ψ_{2} E_{i j}^{2} E + \frac{4}{3} ψ_{3} E_{i k} E_{i l} E_{k l} + ψ_{4} ζ^{3} \\ + ψ_{5} E^{2} ζ + Ψ_{6} E ζ^{2} + \frac{1}{2} ψ_{7} (E - E_{i j}^{2}) ζ, \end{matrix}

(23)

where $Ψ_{1,2, 3}$ are third-ordered elastic constants without fluid, which are equal to the elastic constants $ν_{1, 2, 3}$ proposed [3]. $Ψ_{4,5, 6,7}$ represent elastic constants for characterizing fluid effects in media, corresponding to $γ_{1, 2, 3}$ and $γ$ in Liu’s study [11]. $\bar{K} a n d \bar{μ}$ are functions of effective pressure and frequency, corresponding to $K_{u n i f i e d}, μ_{u n i f i e d}$ mentioned above, which is function of effective pressure and frequency.

When triaxial stress and pore pressure are determined, the constitutive equation of saturated rocks can be expressed as

\begin{matrix} σ_{i j} = (\bar{K} - \frac{2}{3} \bar{u} + α^{2} M) e^{s} δ_{i j} + 2 \bar{u} e_{i j}^{s} - α M ζ^{s} δ_{i j} \\ P_{p} = M (ζ^{s} - α e^{s}) \end{matrix}

(24)

The effective pressure can be expressed as

P = P_{c} - P_{p,}

(25)

where $σ_{i j}, P_{c}, P_{p}$ are triaxial stress, confining pressure, and pore pressure of rocks, $P_{c} = σ / 3$ 。

Assuming that the propagation of waves in rocks is a result of static strain superimposed with small dynamic fields of wave, $u^{s}$ and $w^{s}$ represent the static displacement from the unstressed state to the initial state, $u^{d}$ and $w^{d}$ represent the displacement from the initial state to the final state. The total displacement is $u = u^{s} + u^{d}$ ⁠, the same as $w$ ⁠. According to acoustoelastic theory, static strain is much greater than dynamic strain and ${\partial^{2} u}^{s} / \partial X_{j}^{2} = 0, {\partial^{2} w}^{s} / \partial X_{j}^{2} = 0$ ⁠. From equations (21)–(25), we can derive wave equation of saturated rocks [11]:

\begin{matrix} ρ \frac{\partial^{2} u_{i}^{d}}{\partial t^{2}} + ρ_{f} \frac{\partial^{2} w_{i}^{d}}{\partial t^{2}} = A_{i j k l} u_{k, l j}^{d} + B_{i j} w_{k, k j}^{d} \\ ρ_{f} \frac{\partial^{2} u_{i}^{d}}{\partial t^{2}} + ρ_{+} \frac{\partial^{2} w_{i}^{d}}{\partial t^{2}} = B_{k j} u_{k, l j}^{d} + B w_{k, i k}^{d,} \end{matrix}

(26)

where A_ijkl,B_ij,B_kj, and B are functions of static stress and strain, bulk and shear modulus of rock and third-ordered elastic constants; the specific form is shown in online supplementary Appendix 2.

The two wave functions of solid phase and liquid phase are represented as

\begin{matrix} u^{d} = u_{0} e^{j ω (\frac{n_{i} \cdot X_{i}}{V} - t)} \\ w^{d} = w_{0} e^{j ω (\frac{n_{i} \cdot X_{i}}{V} - t),} \end{matrix}

(27)

$u_{0} a n d w_{0}$ are wave amplitudes, $n_{i}$ is a wave propagation vector, $X_{i}$ is a position vector, t is time, and V is wave phase velocity.

Substituting equation (27) into (26):

\begin{matrix} (Z_{i k} - ρ V^{2} δ_{i k}) u_{k 0} + (N_{i k} - ρ_{f} V^{2} δ_{i k}) w_{k 0} = 0 \\ (Z_{k i} - ρ_{f} V^{2} δ_{i k}) u_{k 0} + (M_{k i} - ρ_{+} V^{2} δ_{i k}) w_{k 0} = 0, \end{matrix}

(28)

Where

Z_{i k} = A_{i j k l} n_{l} n_{j}, N_{i k} = B_{i k} n_{k} n_{j}, N_{k i} = B_{j k} n_{i} n_{j}, M_{k i} = B n_{i} n_{k}

To make the equation have nontrivial solutions, there are:

| \begin{matrix} Z_{i k} - ρ V^{2} δ_{i k} & N_{i k} - ρ_{f} V^{2} δ_{i k} \\ N_{k i} - ρ_{f} V^{2} δ_{i k} & N_{k i} - ρ_{+} V^{2} δ_{i k} \end{matrix} | = 0 (i, k = 1, 2, 3)

(29)

Assuming n = (0,0,1), then two transverse shear velocities (V_s1,s2) can ultimately be obtained:

\begin{matrix} V_{S 1}^{2} = \frac{A_{1313}}{ρ `} \\ V_{S 2}^{2} = \frac{A_{2323}}{ρ `,} \end{matrix}

(30)

Where,

ρ^{'} = ρ - ρ_{f}^{2} / ρ_{+}

And, two complex longitude velocities (V_p1,p2):

V_{P 1, P 2}^{2} = \frac{Λ}{2 (ρ ρ_{+} - ρ_{f}^{2})} [1 \pm \sqrt{1 - \frac{Ξ}{Λ^{2}}}],

(31)

where

\begin{array}{ll} Λ = ρ B + A_{3333} ρ_{+} - 2 ρ_{f} B_{33} \\ Ξ = 4 {(ρ ρ_{+} - ρ_{f}^{2})}^{A_{3333}} \end{array} .

(32)

In most porous media, there is $Λ ≫ Ξ$ [39]. And in general, A_ijkl is much larger than B_ij,B_kj,B, these longitude velocities could be simplified as

\begin{matrix} V_{p 1} \approx \frac{A_{3333}}{ρ^{'}} \\ V_{p 2} \approx \frac{B}{ρ_{+}} \end{matrix} .

(33)

Then corresponding phase velocities $ν_{i}$ and attenuation factor ${1 / Q}_{i}$ are:

v_{i} = {[Re (V_{i}^{- 1})]}^{- 1} and 1/ Q_{i} = Im (V_{i}) / Re (V_{i}), (i = P 1, P 2, S 1, S 2),

(34)

where Re and Im represent the real and imaginary part of a complex number.

3. Model Experiments

3.1 Analysis of Pressure Effect

Some models are applied to test the accuracy of the proposed unified elastic modulus and analyze the dispersion characteristics of the new model in the full frequency range under pressure. First, the rock physics model of Liu et al. [11] was applied for analysis, as shown in Table 1. Seven third-order elastic constants $Ψ_{1, 2, 3, 4, 5, 6 a n d 7}$ are 3378, −1310, −345, 1.93, 220, −110, and −1150 GPa, respectively. Since the experimental data do not contain the pressure-dependent aspect ratio $γ$ ⁠, we choose the value of $γ_{i}$ is 0.005. Then the pressure-dependent $γ$ (P) can be calculated through equation (2). The variations of aspect ratio $γ$ with pore pressure can be seen in Figure 2. It should be noted that some related parameters in this model were obtained under ultrasonic background (approximately 1 MHz). In this part of analysis, the pressure dependence of elastic modulus is the dominate factor, and the mesoscopic dispersion property is ignored (due to it mainly appears in the low-frequency range).

According to equations (5), (8), and (10), we tested the variation of the relevant parameters that form the elastic modulus with pore pressure. Figure 3 shows the variations of strains and pore porosity of compliant pores with pressure when the confining pressure is fixed at 65 MPa. Both the volumetric strain and shear strain increase with the decrease of effective pressure, and the volumetric strain of the bulk modulus is slightly greater than that of the shear modulus. The same variation pattern is also reflected in the compliant porosity.

Usually, because of $θ_{c}$ is much greater than $θ_{s}$ ⁠, $φ_{s} θ_{s}$ is on the order of 0.01, while $φ_{c} θ_{c}$ is on the order of 0.1, $φ_{s} θ_{s}$ in equations (7) and (9) can be omitted. The approximate equations for bulk modulus and shear modulus are rewritten as follows:

\begin{matrix} K_{d r} (P) \approx K_{d r s} [1 - θ_{c} φ_{c 0} e^{(- θ_{c} P C_{d r s})}] \\ μ_{d r} (P) \approx μ_{d r s} [1 - θ_{c μ} φ_{c 0} e^{(- θ_{c} P C_{d r s})}] \end{matrix} .

(35)

According to the model and the calculated parameter information in Figures 2 and 3, the variation of the unified bulk modulus and shear modulus with pressure is shown in Figure 4. Overall, the modulus decreases with the decrease of effective pressure. When the pore pressure is less than 35 MPa, the compliant pores of the rock have closed and the modulus change is gentle. When it is greater than 35 MPa, the compliant pores dominate the modulus of rocks, and the modulus decreases sharply with the increase of pore pressure.

Substitute equation (35) into (17) and determine the frequency as a constant of 1 MHz. According to equations (23)–(25) and (30), (33), the variation of V_p and V_s with pressure is shown in Figure 5. In this figure, the black curve represents the velocity under the influence of the new unified pressure-dependent elastic model in this article, while the blue curve represents the velocity obtained by traditional porous acoustoelastic theory (ignore pressure influence on compliant pores), calculated by Liu et al. [11]. The red scatter represents the measured value. It can be seen that when the pore pressure is less than 25 MPa, corresponding to the case of high effective pressure, both the two methods are consistent with the actual measured values. This situation is due to the closure of compliant pores when pressure is high enough, resulting in the change of $K_{u n i f i e d}$ and $μ_{u n i f i e d}$ is not significant. While, when pressure is below a certain level, the compliant pores in the rock open, and $K_{u n i f i e d}$ and $μ_{u n i f i e d}$ changes significantly with pressure. As a result, there is a significant deviation between the theoretical value of porous acoustoelastic theory and measured values. However, the method proposed in this article is closer to the measured values, as shown in Figure 5. Variation of characteristic frequency $ω_{c}^{m i}$ under different pressure conditions is shown in Figure 6. When the pressure is around 25 MPa, and the rock characteristic frequency is around the measurement frequency of 1 MHz. When the pore pressure is greater than 25 MPa, the characteristic frequency is higher than the measurement frequency. According to equation (17), $K_{u n i f i e d}$ is approximately equal to $K_{d r s}$ when frequency is lower than $ω_{c}^{m i}$ ⁠. It is known that $K_{d r s}$ is measured in a high effective pressure (small pore pressure) situation. So when pressure less than 25 MPa, the velocity calculated by the method in this article is basically consistent with the measured velocity. On the contrary, when the pore pressure is high, the value of $K_{d r s}$ will deviate from the true situation. Meanwhile, the squirt-flow influence when the frequency not reach $ω_{c}^{m i}$ may also lead to errors in the predicted results. These may explain the error between this article’s method and the measured value in the situation when pore pressure is higher than 25 MPa in Figure 5.

3.2.Analysis of Frequency Effect

To analyze the dispersion characteristics of the new unified model in the full frequency range under pressure, the data in Table 1 are used again, but the frequency is no longer a constant and the mesoscopic dispersion is no longer ignored. The mesoscopic heterogeneity is incorporated by a framework of DPDP model, which divides this part of medium into two categories based on different physical properties: phase 1 and phase 2. The relevant mesoscopic parameters are: the scale of the mesoscopic heterogeneity L = 0.4 mm, bulk modulus of heterogeneous phase 1 $K_{d r 1}$ = 1.05 $K_{d r}$ ⁠, bulk modulus of heterogeneous phase 2 $K_{d r 2}$ = 0.5 $K_{d r}$ ⁠, their volume fraction and Skempton’s coefficient are v₁ = 94.55%, v₂ = 5.45%, B₁ = 0.626, B₂ =0.587, the confining pressure is 65 MPa. After knowing these parameters, $K_{u n i f i e d}$ and $μ_{u n i f i e d}$ can be calculated through equations (13)–(15). Then, the velocity and attenuation factor changing with frequency are shown in Figure 7. From the figure, it can be seen that there are two distinct transition regions throughout the full frequency band, corresponding to mesoscopic heterogeneity and microscopic heterogeneity, respectively.

The different colored curves in Figure 7 represent the dispersion characteristics of the model under different pressures states. It can be clearly seen that the characteristic frequency and attenuation at mesoscopic are less affected by pressure. At microscopic, as the pressure increases, the characteristic frequency and attenuation decrease. The influence of pressure on the frequency characteristics of the model is mainly reflected in the high-frequency range.

The mesoscopic dispersion characteristics of heterogeneous rocks are mainly determined by mesoscopic scale parameters. By changing the relative parameters of the new unified model, their effect on the dispersion and attenuation characteristics of the model is analyzed. The confining pressure is 65 MPa and pore pressure is 50 MPa. Figure 8 shows the influence of mesoscale L on the model, and it can be seen that as L continues to increase, the characteristic frequency $ω_{c}^{m e}$ gradually moves toward low frequencies, and the attenuation is not affected by L. Figures 9 and 10 show the influence of bulk modulus of heterogeneous phases 1 and phase 2 on mesoscopic dispersion and attenuation, respectively. As can be seen, the influence of bulk modulus respect to mesoscopic heterogeneity in both figures is similar: at the mesoscopic, the attenuation increases with the increase of modulus, while the characteristic frequency $ω_{c}^{m e}$ remains unchanged.

In addition, to verify the accuracy of the new unified model, a rock physics model is used for test and its data were obtained through experimental measurements by Borgomano et al. [40], covering a wide frequency range from low to high frequency. The parameters are shown in Table 2. The applied pressure is: $P_{c}$ = 12 MPa, $P_{p}$ = 2 MPa. Since the experimental data do not provide relevant parameters at mesoscopic, we refer to the work of Zhang et al. [24]. In their study, they assume the mesoscopic heterogeneity is caused by the flow squirts from the clay phase to the stiff phase. And in their experiment, they defined volume fraction of clay (5.45%) and skempton’s coefficient are 0.625 of stiff phase and 0.587 of clay phase. The rock matrix bulk moduli of stiff phase to the clay phase are 13.6 GPa and 6 GPa. The representative volume element of are the harmonic means of those of the two phases and equal to 12.8 GPa. Therefore, we have chosen the values of these parameters within a reasonable range: L = 0.7 mm, $K_{d r 1}$ = 1.05 $K_{d r}$ ⁠, $K_{d r 2}$ = 0.5 $K_{d r}$ ⁠, v₁ = 95%, v₂ = 5%, B₁ = 0.63, B₂ = 0.58.

Figure 11 shows the comparison between the bulk modulus and attenuation calculated using the unified pressure- and frequency-dependent elastic model proposed in this article and the experimental measurement results. Figure 11(a) shows the variation of bulk modulus with frequency, and Figure 11(b) shows the attenuation factor of bulk modulus $Q_{K}^{- 1}$ = Im(K)/Re(K). The red curve in these figure represents the calculation results of the method proposed in this article, while the blue curve represents the results obtained using the modified Shapiro model, which is also a pressure- and frequency-dependent elastic model but does not take the mesoscopic dispersion of the medium into account [31]; the black scatter is the experimental measurement data. From the distribution of measurement data points in Figure 11(a), it can be seen that there are two characteristic frequencies of the rock, corresponding to the two attenuated peaks in Figure 11(b). The Shapiro’s model only considers the heterogeneity of rock at the microscopic, so it only matches the measurement data well in the high-frequency part and lacks accurate characterization of the low-frequency part. The method presented in this study can better reflect the phenomenon of dispersion and attenuation caused by mesoscopic heterogeneity of rocks, thus having great advantages in describing the nonlinear heterogeneous properties in fluid-saturated porous media.

4. Discussion

Propagation characteristics in fluid-saturated porous rocks under hydrostatic environment were discussed above. When external pressures or stresses are spatially nonequilibrium, which are generally exist in subsurface media, the rock exhibits anisotropy [41]. In this case, the rock’s velocities along different propagation directions exhibit different dispersion and attenuation characteristics. Assuming that parameters of the rock are isotropic, only the stress applied to rock in the three axes are different and $P_{c}$ = (⁠ $σ_{11}$ + $σ_{22}$ + $σ_{33}$ ⁠)/3, $σ_{i i}$ represents stress in different directions, and the variation of rock velocity and attenuation factor in each direction with frequency is shown in Figure 12. The blue curve in the figure represents direction n = (1, 0, 0) with $σ_{11}$ = 65 MPa, the red curve represents direction n = (0, 1, 0) with $σ_{22}$ = 55 MPa, and the black curve represents direction n = (0, 0, 1) with $σ_{33}$ = 60 MPa. By comparing the velocity and attenuation factor of three curves across the full frequency range, it can be seen that the greater the stress, the higher the rock velocity along this direction, however, the less significant the changes in velocity and attenuation factor. Figure 13 shows a polarization diagram of the variation of rock phase velocities with orientation at a specific frequency (with a certain incident angle of 90^o), and the triaxial stress applied to rock is the same as the situation in Figure 12. The P-wave velocity is maximum along the direction of maximum horizontal stress, while the two s-wave velocities are maximum along the direction of maximum horizontal stress and minimum horizontal stress, respectively. From Figures 12 and 13, the stress and frequency-dependent velocity distribution characteristics of rocks are preliminary analyzed. And these characteristics related to principle stressed may provide a basis for in situ stress field detection of subsurface reservoirs.

5. Conclusions

A new unified elastic model is proposed to describe the dispersion and attenuation phenomenon of fluid-saturated porous rocks under pressure. New model incorporates a dual-porosity theory into a modified Zener model to accurately reflect the influence of rock properties under pressure and frequency coupling conditions. The nonlinear deformation of compliant pores in the dual-porosity model reflects the pressure dependence. The modified Zener model unifies the frequency variation characteristics of rocks at both mesoscale and microscales, reflecting the frequency dependence of rocks throughout full frequency band. The Lagrangian function considering rock kinetic energy and strain energy is applied to establish the dynamic equation of fluid-saturated porous rocks using porous acoustoelastic theory. Substitute the new unified elastic model into the dynamic equation to solve for the various complex phase velocities. Under this new framework, the variation of rock velocities (or elastic modulus) and attenuation factor with pressure and frequency can be more accurately obtained, and the wave propagation mechanism in fluid-saturated porous media under pressure influence can be more comprehensively understood. The velocities or modulus calculated from the new unified model was compared with traditional methods and measured experimental data. The comparison results show that in terms of pressure effect, compared with the porous acoustoelastic model, the new method has a higher degree of agreement with the measured data, especially in the low effective pressure region; in terms of frequency variation characteristics, compared with the improved Shapiro’s model, the new method can accurately reflect the dispersion and attenuation phenomenon within the seismic frequency band while ensuring the accuracy of high-frequency areas. These experimental analyses demonstrate the superiority of the new method in describing the propagation mechanism of fluid-saturated porous media under pressure–frequency coupling effects.

Data Availability

The data used in this study are included or cited within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

This research was supported by the Laoshan National Laboratory of Science and Technology Foundation (No. LSKJ202203400) and the National Natural Science Foundation Project of China (No. 42030103).

Acknowledgments

Thanks to the two anonymous reviewers for their constructive input and the associate editor Gang Fang for the efficient handling of this paper.

Supplementary Materials

Appendix A: DPDP model and relative parameters.

Appendix B: Elastic tensor in acoustoelastic theory.

Properties	Values
Porosity $φ$	0.197
Tortuosity $τ$	2.848
Fluid viscosity (Pa.s) $η$	0.00106
Rock density (kg/m³) $ρ_{s}$	2270
Water density (kg/m³) $ρ_{f}$	1000
Permeability (md) $κ$	43.50
Bulk moduli of rock matrix (GPa) $K_{g}$	37
Shear moduli of rock matrix (GPa) $μ_{g}$	44
Bulk moduli of fluid (GPa) $K_{f}$	2.25
Bulk moduli of rock skeleton (GPa) $K_{d r s}$	13.2
Shear moduli of rock skeleton (GPa) $μ_{d r s}$	5.4

Properties	Values
Porosity $φ$	0.197
Tortuosity $τ$	2.848
Fluid viscosity (Pa.s) $η$	0.00106
Rock density (kg/m³) $ρ_{s}$	2270
Water density (kg/m³) $ρ_{f}$	1000
Permeability (md) $κ$	43.50
Bulk moduli of rock matrix (GPa) $K_{g}$	37
Shear moduli of rock matrix (GPa) $μ_{g}$	44
Bulk moduli of fluid (GPa) $K_{f}$	2.25
Bulk moduli of rock skeleton (GPa) $K_{d r s}$	13.2
Shear moduli of rock skeleton (GPa) $μ_{d r s}$	5.4

Porperties	Values
Porosity $φ$	0.075
Crack porosity $φ_{c}$	0.0015
Fluid viscosity (Pa.s) $η$	1
Rock density (kg/m³) $ρ_{s}$	2540
Crack density (kg/m³) $ρ_{f}$	1300
Permeability (m²) $κ$	5 × 10^-17
Aspect ratio $γ$	2.75 × 10^-4
Bulk moduli of rock matrix (GPa) $K_{g}$	77
Shear moduli of rock matrix (GPa) $μ_{g}$	84
Bulk moduli of fluid (GPa) $K_{f}$	4.36
Bulk moduli of rock skeleton (GPa) $K_{d r s}$	58

Porperties	Values
Porosity $φ$	0.075
Crack porosity $φ_{c}$	0.0015
Fluid viscosity (Pa.s) $η$	1
Rock density (kg/m³) $ρ_{s}$	2540
Crack density (kg/m³) $ρ_{f}$	1300
Permeability (m²) $κ$	5 × 10^-17
Aspect ratio $γ$	2.75 × 10^-4
Bulk moduli of rock matrix (GPa) $K_{g}$	77
Shear moduli of rock matrix (GPa) $μ_{g}$	84
Bulk moduli of fluid (GPa) $K_{f}$	4.36
Bulk moduli of rock skeleton (GPa) $K_{d r s}$	58

Abstract

1. Introduction