An improved wave equation of fractured-porous media for predicting reservoir permeability

The wave characteristics of fractured-porous media can be utilized for permeability identification; however, further research is necessary to enhance the accuracy of this identification. A novel wave equation for fractured-porous media is formulated, and theoretical analysis suggests its effectiveness in accurately identifying reservoir permeability. The proposed methodology establishes a wave equation for fractured-porous media using the volume averaging method and employs finite difference method on staggered grids to calculate wave field dispersion and attenuation, exploring the influence of fracture network structure and confining pressure on the solution of the wave equation. By analyzing the wave equation under various aspect ratios and confining pressure of fractures, it is observed that these factors significantly affect velocity and attenuation, providing valuable insights into seismic response in fractured-porous media. Furthermore, the research findings reveal promising potential in utilizing the new wave equations specific to fractured-porous media for permeability identification purposes. By constructing a three-dimensional fractured-porous network model, the wave equation for permeability identification can examine the correlation between the parameters of the equation and permeability, and establishes an association between fracture parameters and permeability, paving the way for a novel approach to permeability identification.

1 Introduction

The permeability is one of the crucial parameters for evaluating reservoir storage performance. Previous research has demonstrated that analyzing wave behavior in fractured-porous media can improve the accuracy of identifying reservoir permeability. The prediction of reservoir permeability using the wave equation in fractured-porous media can be categorized into three distinct stages. In the initial stage, Biot (1956a, b), based on the homogeneous porous media model and discovered that variations in permeability can lead to frequency dispersion and attenuation of seismic wave, resulting in changes in seismic velocity, which provided the foundation for subsequent related studies. In the subsequent stage, Numerous fractured-porous media models (Schoenberg, 1980; Schoenberg, 1983; Schoenberg and Douma, 1988; Schoenberg and Sayers, 1995; Gurevich et al., 2009; Tang, 2011; Tang et al., 2012) on Biot’s theory have been developed to investigate the relationship between permeability and seismic velocity in fractured-porous media. For instance, Chapman (Chapman et al., 2002; Chapman et al., 2016) established a microscale porous model, yet failed to provide an expression for the permeability of the porous media. Chichinina (Chichinina et al., 2007; Chichinina et al., 2009) proposed an anisotropic model for horizontally layered fractures that focused on P/S wave attenuation but did not address the permeability in fracture system. Gurevich (Gurevich et al., 2009), using a constructed fractured-porous media model, observed significant variations in velocity due to changes in the structures of pore and fracture, while changes in permeability resulted in relatively minor alterations in wavefield velocity. However, these studies did not explicitly explore the relationship between velocity and microstructure or permeability of fractured-porous media due to limitations in computer technology at that time. In the third stage, researchers initiated the construction of fractured-porous models incorporating microstructure, studying the relationship between fractures and wavefield velocity as well as reservoir permeability. For example, Du et al. (2011) studied an equivalent media model for fractured porous rocks. Fractures are modeled by constitutive relationship in terms of fracture-induced anisotropy. Guo et al. (2018) used the branch function method to construct finite-thickness fractures and examined the dispersion and attenuation of P-wave propagation perpendicular to the fracture surface. Lissa et al. (2019) explored the impact of fractures with varying widths on seismic attenuation and velocity dispersion. Song et al. (2020) investigated dispersion and attenuation of P-wave in heterogenous porous media containing oriented fractures, revealing that factors such as mechanical form of fractures and fluid flow properties significantly influence attenuation sensitivity. Wei et al. (2021) conducted research on a 3D fractured-porous network model, presenting a calculation approach for permeability while uncovering notable influences of aspect ratio on P-wave velocity and characteristic frequency. Wang et al. (2022) studied the propagation law for complex fractures in two-phase media through linear slip theory.

The prediction of reservoir permeability in fractured-porous media using the wave equation still encounters numerous challenges and difficulties. Firstly, it is imperative to enhance the microstructure and multiscale characteristics of fractured-porous media (Bai et al., 2021). Secondly, it is crucial to investigate the influences of factors such as fractures, confining pressure, permeability, and porosity on velocity as well as explore their interactions in fractured-porous media. Only through these efforts can we elucidate the relationship between wave properties and permeability in fractured-porous media.

The present work proposes a novel fractured-porous model and derives its wave equation based on volume averaging method, investigating the influence of fracture network structure on the solution of the wave equation and the relationship between fracture aspect ratio and wave velocity as well as permeability. This research possesses several distinctive features compared to previous studies: firstly, the solutions of the wave equation under different fracture aspect ratios and confining pressure conditions are studied; secondly, by analyzing various parameter combinations, such as velocity and attenuation, significant effects of fracture aspect ratio and confining pressure on wave behavior are revealed; thirdly, a proposed permeability identification method based on the curves of velocity-fracture parameters and velocity-permeability establishes correlations between fracture parameters and permeability.

2 Elliptic cylindrical fracture model and its fluid motion equation

Fractures exhibit a predominant extension direction, thus the elliptic cylindrical model can be used to describe their extensibility. The aspect ratio, defined as the ratio between the short and long radii of the elliptical cylinder cross-section, governs the structure of this model. As the aspect ratio approaches unity, it degenerates into a conventional porous model; whereas for very small values of aspect ratio, it accentuates more features of fracture structures. This flexible elliptical cylinder model enables us to depict various characteristics of fractures in reservoirs under different aspect ratios and provides insights into fracture structures with varying extensibility.

The derivation of the fluid motion equations in elliptic cylindrical fracture structure is presented below. The incompressible Newtonian fluid within the microtubule space of the elliptical cylinder satisfies the following equations (Landau and Lifshitz, 1987).

The equation of fluid mass conservation:

$\frac{\partial {\rho }_f}{\partial t}+\nabla \cdot \left({\rho }_f\dot{\mathbf{U}}\right)=0(1)$

The equation of fluid momentum conservation:

$\nabla \cdot \mathbf{s}={\rho }_f\frac{\partial \mathbf{U}}{\partial t}+\nabla p(2)$

Constitutive relationship of fluid mechanics:

$\mathbf{s}=\eta \left[\nabla \dot{\mathbf{U}}+{\left(\nabla \dot{\mathbf{U}}\right)}^{\mathrm{T}}-\frac{2}{3}\left(\nabla \cdot \dot{\mathbf{U}}\right)\mathbf{I}\right](3)$

where ${\rho }_{\mathit{f}}$ represents the fluid density, $\mathbf{s}$ denotes the fluid stress tensor, and $p$ signifies the fluid pressure, $\mathbf{U}$ is the displacement of fluid, $\dot{\mathbf{U}}$ represents the flow velocity, $\eta$ stands for the fluid viscosity coefficient, and $D\dot{\mathbf{U}}/Dt$ corresponds to the material derivative of fluid velocity field with respect to time. Divergence is taken on both sides of the fluid constitutive relation equation to obtain the fluid dynamic equation: By applying divergence operator on both sides of the constitutive equation for fluid, the equation of fluid dynamics is yielded:

${\rho }_f\frac{D\dot{\mathbf{U}}}{Dt}=-\nabla p+\eta \left({\nabla }^2\dot{\mathbf{U}}+\frac{1}{3}\nabla \left(\nabla \cdot \dot{\mathbf{U}}\right)\right)(4)$

The fluid velocity field is expressed as the combination of steady-state and unsteady-state fields:

$\dot{\mathbf{U}}\left(x,y,t\right)={\dot{\mathbf{U}}}_0\left(x,y,z\right)+{\dot{\mathbf{U}}}_1\left(x,y,z,t\right)(5)$

Under low-speed flow conditions, the compressibility of pore fluid can be neglected when subjected to elastic waves (denoted as $\nabla \cdot \dot{\mathbf{U}}=0$). Considering that the only non-zero velocity component is along the axis of the elliptical cylinder and the flow velocity depends solely on radial coordinate, therefore, $D\dot{\mathbf{U}}/Dt=\partial \dot{\mathbf{U}}/\partial t$ is obtained.

The fluid dynamics equations mentioned above can be solved by the series expansion of Mathieu function in the elliptic cylindrical coordinate system, thereby obtaining the velocity field and flow rate in elliptic cylindrical space (Xiong et al., 2021). The expression for the steady-state flow rate at both ends of the elliptical cylinder is:

$\left\{\begin{array}{c}{Q}_U=\frac{-\mathrm{i}\pi R^2}{{\rho }_{f}c}\left(\frac{P_U\cos \left(\omega L/c\right)-P_D}{\sin \left(\omega L/c\right)}\right)\left(\frac{2J_1\left(KR\right)}{KR{J}_0\left(KR\right)}-1\right)\\ Q_D=\frac{-\mathrm{i}\pi R^2}{{\rho }_{f}c}\left(\frac{P_U-P_D\cos \left(\omega L/c\right)}{\sin \left(\omega L/c\right)}\right)\left(\frac{2J_1\left(KR\right)}{KR{J}_0\left(KR\right)}-1\right)\end{array}\right.(6)$

The expression describing the flow rate of the pulsating flow field is as follows:

$\left\{\begin{array}{c}{Q}_{Ue}=2\left(\frac{a^3}{1+a^2}\right)\frac{-\mathrm{i}\pi R^2}{{\rho }_{f}c}\left(\frac{P_U\cos \left(\omega L/c\right)-P_D}{\sin \left(\omega L/c\right)}\right)\left(\frac{2J_1\left(KR\right)}{KR{J}_0\left(KR\right)}-1\right)\\ Q_{De}=2\left(\frac{a^3}{1+a^2}\right)\frac{-\mathrm{i}\pi R^2}{{\rho }_{f}c}\left(\frac{P_U-P_D\cos \left(\omega L/c\right)}{\sin \left(\omega L/c\right)}\right)\left(\frac{2J_1\left(KR\right)}{KR{J}_0\left(KR\right)}-1\right)\end{array}\right.(7)$

The flow rate of the elliptic cylindrical fractures is observed to be dependent on various factors, including the length $L$ of the elliptical cylinder, pressures $P_U$ and $P_D$ at both ends, major radius $R$ and aspect ratio $\alpha$ of the elliptical cross-section, characteristic frequency $\omega$ of non-steady flow field, and fluid density ${\rho }_f$. In Equation 7, $J_0$ and $J_1$ represent zero- and first-order Bessel functions respectively, $K=\sqrt{\mathrm{i}\omega {\rho }_f/\eta }$, and $c$ denotes the wave velocity in fluid.

The elliptical curve representation is inadequate for accurately depicting the tip of a genuine fracture. In order to precisely depict the variation in the major axis direction radius $R_1$ of actual fracture under confining pressure, an improved method for calculating the fracture permeability that not only considers the deformation of wedge-shaped fractures under confining pressure (Mavko and Nur, 1978), but also incorporates the relationship between fracture radius and confining pressure (Xiong et al., 2021), thereby reflecting the variation in flow rate of wedge-shaped fractures with changing confining pressure.

Confining pressure signifies the force exerted by overlying geological materials on a specific point within the Earth’s crust. This pressure significantly affects permeability and elasticity of porous media. Firstly, under elevated confining pressure, the closure of pores and fractures within porous media leads to a significant reduction in permeability, diminishing the capacity for fluid flow. Secondly, confining pressure alters the elasticity of the porous material, specifically impacting its bulk modulus. This shift in elasticity has ramifications for the mechanical behavior of the material and the propagation of stress-related phenomena, including seismic waves.

In this study, it is assumed that flow conservation holds for each node in the fractured-porous network, ensuring that inflow equals outflow at every node. A flow conservation equation is formulated for each node, resulting in a system of linear equations with unknown variables representing the pressures at individual nodes. The pressures at the inlet and outlet are considered as non-homogeneous terms. By solving this system of equations to determine the pressures at each node, the flow within the 3D porous network can be accurately calculated. Furthermore, by incorporating Darcy’s law, dynamic permeability prediction for the 3D microtubule network in porous media can be effectively conducted (Xiong et al., 2021).

3 Wave equation for fractured-porous media

The volume averaging theorem, established by researchers like Whitaker (Whitaker, 1966; Whitaker, 1967), Slattery (1967), and Gray and Lee (1977), plays a pivotal role in connecting microscale parameters to macroscale behaviors within porous media. It has been widely applied in fields ranging from hydrogeology to geophysics, providing a consistent framework for modeling complex porous media systems. This theorem’s mathematical rigor and practical relevance have been extensively validated through real-world applications, underlining its significance. The use of the volume averaging theorem is firmly grounded in this well-established foundation of literature and research, reinforcing its vital role in this work and aligning with the reviewer’s emphasis on its importance. By introducing the volume average theorem, a connection between the macroscopic behavior of fractured-porous media and the microscopic parameters of the media is established, which derives a macroscopic wave equation that describes the wave behavior based on this theorem.

3.1 Volume average theorem of porous media

The unit cell region $\mathrm{\Omega }$, filled with fluid in the porous media, is considered. Its position $\mathbf{x}$ represents the “center” or “centroid”. In the fluid-filled porous media of the 3D elliptical microtubule network, the characteristic length of region $\mathrm{\Omega }$ 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\ll H_h\ll 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 $\mathrm{\Omega }$ 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 ${\psi }_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 ${\psi }_f$ within the entire region $\mathrm{\Omega }$, as illustrated below.

$〈{\psi }_f〉=\frac{1}{V}{\int }_{\mathrm{\Omega }}{\psi }_f\left(\mathbf{x}\right)\mathrm{d}V(9)$

where $〈{\psi }_f〉$ is characterized by the smoothness with respect to the centroid $\mathbf{x}$ of region $\mathrm{\Omega }$. Another variable related to the volume average is the volume-averaging eigenvalue, defined as follows:

${\overline{\psi }}_f=\frac{1}{V_f}{\int }_{\mathrm{\Omega }}{\psi }_f\left(\mathbf{x}\right)\mathrm{d}V(10)$

According to the above definition, there is $\phi =〈{\psi }_{\mathit{f}}〉/{\overline{\psi }}_{\mathit{f}}$.

According to the research conducted by Slattery (Slattery, 1967) and Whitaker (Whitaker, 1966) on the volume average theorem, considering the physical variable ${\psi }_f$ as an illustrative example, the volume average term including the spatial derivative is as follows:

${\int }_{\mathrm{\Omega }}\nabla \cdot {\mathbf{\psi }}_f\mathrm{d}V=\nabla \cdot {\int }_{\mathrm{\Omega }}{\mathbf{\psi }}_f\mathrm{d}V+{\int }_{\mathrm{\Sigma }}{\mathbf{\psi }}_f\cdot \mathbf{n}\mathrm{d}S(11)$

The volume average term including time derivative is as follows:

${\int }_{\mathrm{\Omega }}{\partial }_t{\psi }_f\mathrm{d}V={\partial }_t{\int }_{\mathrm{\Omega }}{\psi }_f\mathrm{d}V-{\int }_{\mathrm{\Sigma }}{\psi }_f\mathbf{w}\cdot \mathbf{n}\mathrm{d}S(12)$

where $\mathrm{\Sigma }$ represents the interface where solid and fluid come into contact within the region Ω. It is typically specified that the normal vector $\mathbf{n}$ points outward from fluid, and $\mathbf{w}$ denotes the velocity vector of the interface between solid and fluid, which is generally taken as zero vector.

3.2 Microscopic equation for fractured-porous media

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.

${\rho }_s\frac{{\partial }^2\mathbf{u}}{\partial t^2}=\nabla \cdot \mathbf{\sigma }(13)$

where ${\rho }_{\mathit{s}}$ denotes the mass density of solid skeleton, and $\mathbf{u}$ represents the displacement of solid skeleton, $t$ refers to time and $\mathbf{\sigma }$ signifies solid stress tensor.

The microscopic constitutive equation of solid is:

$\mathbf{\sigma }=K_{s}e\cdot \mathbf{I}+G\left[\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\mathrm{T}}-\frac{2}{3}e\cdot \mathbf{I}\right](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 $\mathrm{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:

$\frac{1}{K_f}\frac{\mathrm{d}p}{\mathrm{d}t}=\frac{1}{{\rho }_f}\frac{\mathrm{d}{\rho }_f}{\mathrm{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\nabla \cdot \mathbf{u}(16)$

Assuming the presence of a porous media comprising a solid skeleton and a fluid, with only one existing interface denoted as $\mathrm{\Sigma }$, 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:

$\begin{array}{l}\dot{\mathbf{U}}\cdot \mathbf{n}=\dot{\mathbf{u}}\cdot \mathbf{n}\\ -p\cdot \mathbf{n}+\mathbf{s}\cdot \mathbf{n}=\mathbf{\sigma }\cdot \mathbf{n}\end{array}(17)$

where $\mathbf{n}$ denotes the outward normal vector of contact interface, and s is the deviatoric stress in fluid.

3.3 Macroscopic integral equation of fractured-porous media

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:

$\frac{{\partial }^2}{\partial t^2}〈{\rho }_f\mathbf{U}〉=-\nabla 〈p〉+\nabla \cdot 〈\mathbf{s}〉+\overline{p}\nabla \varphi -\frac{\eta {\varphi }^2}{\kappa }\left(\overline{\dot{\mathbf{U}}}-\overline{\dot{\mathbf{u}}}\right)(18)$

where $\kappa$ represents the permeability of porous media. Since fluid does not undergo shear deformation, the fluid stress is equivalent to the normal stress $\mathbf{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:

$\frac{{\partial }^2}{\partial t^2}〈{\rho }_s\mathbf{u}〉=-\nabla \cdot 〈\mathbf{\sigma }〉-\overline{p}\nabla \varphi +\frac{\eta {\varphi }^2}{\kappa }\left(\overline{\dot{\mathbf{U}}}-\overline{\dot{\mathbf{u}}}\right)(19)$

The stress $\mathbf{\sigma }$ of solid skeleton can be decomposed into two components: normal stress ${\mathbf{\sigma }}_s$ and shear stress ${\mathbf{\sigma }}_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:

${\mathbf{\sigma }}_s=-\left(1-\phi \right){\overline{p}}_s=a_{11}\nabla \cdot \overline{\mathbf{u}}+a_{12}\nabla \cdot \overline{\mathbf{U}}(20)$

$\mathbf{s}=-\phi \overline{p}=a_{21}\nabla \cdot \overline{\mathbf{u}}+a_{22}\nabla \cdot \overline{\mathbf{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}=-\frac{{\left(1-\phi \right)}^2{K}_f}{\phi L_1}-\frac{{\left(1-\phi \right)}^2{K}_f}{L_1\left[\left(1-\phi \right)K_s-K_b\right]}(22)$

$a_{12}=-\frac{\left(1-\phi \right)K_f}{L_1}(23)$

$a_{21}=-\frac{\phi K_s}{L_2}+\frac{\phi K_b}{\left(1-\phi \right)L_2}(24)$

$a_{22}=-\frac{{\phi }^2{K}_s}{\left(1-\phi \right)L_2}(25)$

$L_1=-\frac{\left(1-\phi \right)K_f}{\phi K_s}-\frac{\left(1-\phi \right)K_s}{\left(1-\phi \right)K_s-K_b}(26)$

$L_2=-\frac{\phi K_s}{\left(1-\phi \right)K_f}-1+\frac{K_b}{\left(1-\phi \right)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 $\nu$ can be determined (Mavko and Nur, 1978):

$\frac{1}{K^{\prime }}=\frac{1}{K_s}\left[1+\frac{2\pi }{3}\frac{\left(1-{\nu }^2\right)}{\left(1-2\nu \right)}{\displaystyle \sum _{i=1}^{N_c}}\frac{{R_1}_i^2{d}_i}{V}\right](28)$

where ${{\mathit{R}}_1}_{\mathit{i}}$ refers to the radius of the elliptical principal axis of the $i-\text{th}$ fracture, and $d_i$ stands for fracture extension length. In a cubic unit cell with side length $l$, let $M\times N\times L$ represent the spatial dimensions of the fracture network, where $M$ fractures are distributed in the $\mathrm{X}$ direction, $N$ fractures in the $\mathrm{Y}$ direction and $L$ fractures in the $\mathrm{Z}$ direction. The relationship among the main radius of fracture network model, porosity $\varphi$, aspect ratio of fractures $a$, network scale and sample side length shows as follows (Xiong et al., 2021):

$R_1=\sqrt{\frac{\left(MN+NL+ML\right)\phi l^2}{\pi a}}(29)$

Assuming that the shear stress in the porous media is solely borne by the solid skeleton, disregarding any fluid-induced shear deformation:

${\mathbf{\sigma }}_e=G\left(\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^T-\frac{2}{3}\left(\nabla \cdot \mathbf{u}\right)\mathbf{I}\right)(30)$

Where $G$ denotes the shear modulus of solid skeleton and $\mathbf{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:

$\left\{\begin{array}{c}〈\mathbf{\sigma }〉=\left(a_{11}\nabla \cdot \overline{\mathbf{u}}+a_{12}\nabla \cdot \overline{\mathbf{U}}\right)\mathbf{I}+G\left(\nabla \overline{\mathbf{u}}+{\left(\nabla \overline{\mathbf{u}}\right)}^{\mathrm{T}}-\frac{2}{3}\left(\nabla \cdot \overline{\mathbf{u}}\right)\mathbf{I}\right)\\ 〈\mathbf{s}〉=\left(a_{21}\nabla \cdot \overline{\mathbf{u}}+a_{22}\nabla \cdot \overline{\mathbf{U}}\right)\mathbf{I}\end{array}\right.(31)$

The final form of the wave equation is derived by combining the constitutive equation with the volume-averaged momentum conservation equation:

$\left\{\begin{array}{c}〈{\rho }_s〉\ddot{\mathbf{u}}=G{\nabla }^2\mathbf{u}+a_{11}^{\prime }\nabla e+a_{12}\nabla \xi +C\left(\dot{\mathbf{U}}-\dot{\mathbf{u}}\right)\\ 〈{\rho }_f〉\ddot{\mathbf{U}}=a_{21}\nabla e+a_{22}\nabla \xi -C\left(\dot{\mathbf{U}}-\dot{\mathbf{u}}\right)\end{array}\right.(32)$

where, $C=\eta {\varphi }^2/\kappa$, $a_{11}^{\prime }=a_{11}+4G/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 $\mathbf{u}$ and $\mathbf{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 ${\mathbf{u}}_P=\nabla \cdot \mathbf{u}$ and ${\mathbf{U}}_P=\nabla \cdot \mathbf{U}$, the P-wave equations is obtained:

$\left\{\begin{array}{c}〈{\rho }_s〉\frac{{\partial }^2{\mathbf{u}}_P}{\partial t^2}=a_{11}^{\prime }{\nabla }^2{\mathbf{u}}_P+a_{12}{\nabla }^2{\mathbf{U}}_P+C\left(\frac{\partial {\mathbf{U}}_P}{\partial t}-\frac{\partial {\mathbf{u}}_P}{\partial t}\right)\\ 〈{\rho }_f〉\frac{{\partial }^2{\mathbf{U}}_P}{\partial t^2}=a_{21}{\nabla }^2{\mathbf{u}}_P+a_{22}{\nabla }^2{\mathbf{U}}_P-C\left(\frac{\partial {\mathbf{U}}_P}{\partial t}-\frac{\partial {\mathbf{u}}_P}{\partial t}\right)\end{array}\right.(33)$

By applying the curl operator to both sides of the wave Eq. 32 and introducing ${\mathbf{u}}_S=\nabla \times \mathbf{u}$ and ${\mathbf{U}}_S=\nabla \times \mathbf{U}$, the s-wave equations is derived:

$\left\{\begin{array}{c}〈{\rho }_s〉\frac{{\partial }^2{\mathbf{u}}_S}{\partial t^2}=G{\nabla }^2{\mathbf{u}}_S+C\left(\frac{\partial {\mathbf{U}}_S}{\partial t}-\frac{\partial {\mathbf{u}}_S}{\partial t}\right)\\ 〈{\rho }_f〉\frac{{\partial }^2{\mathbf{U}}_S}{\partial t^2}=-C\left(\frac{\partial {\mathbf{U}}_S}{\partial t}-\frac{\partial {\mathbf{u}}_S}{\partial t}\right)\end{array}\right.(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:

$\left\{\begin{array}{c}{\mathbf{u}}_P={\mathbf{u}}_{P0}{\mathrm{e}}^{i\left(\omega t-\mathbf{k}\cdot \mathbf{x}\right)}\\ {\mathbf{U}}_P={\mathbf{U}}_{P0}{\mathrm{e}}^{i\left(\omega t-\mathbf{k}\cdot \mathbf{x}\right)}\end{array}\right.(35)$

where ${\mathbf{u}}_{P0}$ and ${\mathbf{U}}_{P0}$ refer to P-wave amplitude, and $\mathbf{k}$ denotes wave number, which is typically a complex quantity in nature; $\omega$ signifies circular frequency, and $i$ represents imaginary unit. After substituting the above equations into the P-wave equation, the subsequent expression can be derived:

$\left\{-{\omega }^2\left[\begin{array}{cc}〈{\rho }_s〉& 0\\ 0& 〈{\rho }_f〉\end{array}\right]+k^2\left[\begin{array}{cc}{a}_{11}^{\prime }& a_{12}\\ a_{21}& a_{22}\end{array}\right]+\mathrm{i}\omega \left[\begin{array}{cc}-C& C\\ C& -C\end{array}\right]\right\}\left[\begin{array}{c}{\mathbf{u}}_{P0}\\ {\mathbf{U}}_{P0}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right](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=\left(〈{\rho }_s〉+\frac{\mathrm{i}C}{\omega }\right)\left(〈{\rho }_f〉+\frac{\mathrm{i}C}{\omega }\right)+\frac{C^2}{{\omega }^2}(38)$

$Z_2=-a_{11}^{\prime }〈{\rho }_f〉-a_{22}〈{\rho }_s〉-\frac{\mathrm{i}C\left(a_{11}^{\prime }+a_{12}+a_{21}+a_{22}\right)}{\omega }(39)$

$Z_3=a_{11}^{\prime }{a}_{22}-a_{12}{a}_{21}(40)$

where $X={\omega }^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={\left(E^2+F^2\right)}^{1/4}\cdot sign\left(\cos \frac{\alpha }{2}\right)\cdot \sqrt{\frac{1}{2}+\frac{E}{2\sqrt{E^2+F^2}}}(42)$

$V_P^I={\left(E^2+F^2\right)}^{1/4}\cdot sign\left(\sin \frac{\alpha }{2}\right)\cdot \sqrt{\frac{1}{2}-\frac{E}{2\sqrt{E^2+F^2}}}(43)$

$\tan \alpha =\frac{F}{E}(44)$

where the coefficients $E$ and $F$ are

$E=\frac{2\left(a_{11}^{\prime }{\overline{\rho }}_f+a_{22}{\overline{\rho }}_s\right){\overline{\rho }}_s{\overline{\rho }}_f+2\left(a_{11}^{\prime }+a_{12}+a_{21}+a_{22}\right)\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)\frac{C^2}{{\omega }^2}\pm 2{\overline{\rho }}_s{\overline{\rho }}_{f}A\pm 2\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)\frac{C}{\omega }B}{4{\overline{\rho }}_s^2{\overline{\rho }}_f^2+4{\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)}^2\frac{C^2}{{\omega }^2}}(45)$

$F=\frac{2\left(a_{11}^{\prime }+a_{12}+a_{21}+a_{22}\right){\overline{\rho }}_s{\overline{\rho }}_f\frac{C}{\omega }-\left[2\left(a_{11}^{\prime }{\overline{\rho }}_f+a_{22}{\overline{\rho }}_s\right)\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)\frac{C}{\omega }\pm 2\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)\frac{C}{\omega }A\right]\pm 2{\overline{\rho }}_s{\overline{\rho }}_{f}B}{4{\overline{\rho }}_s^2{\overline{\rho }}_f^2+4{\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)}^2\frac{C^2}{{\omega }^2}}(46)$

$A={\left(x^2+y^2\right)}^{\frac{1}{4}}\cdot sign\left(\cos \frac{\theta }{2}\right)\cdot \sqrt{\frac{1}{2}+\frac{x}{2\sqrt{x^2+y^2}}}(47)$

$B={\left(x^2+y^2\right)}^{\frac{1}{4}}\cdot sign\left(\sin \frac{\theta }{2}\right)\cdot \sqrt{\frac{1}{2}-\frac{x}{2\sqrt{x^2+y^2}}}(48)$

$\tan \theta =\frac{y}{x}(49)$

$x={\left(a_{11}^{\prime }{\overline{\rho }}_f+a_{22}{\overline{\rho }}_s\right)}^2-4\left(a_{11}^{\prime }{a}_{22}-a_{12}{a}_{21}\right){\overline{\rho }}_f{\overline{\rho }}_s-{\left(a_{11}^{\prime }+a_{12}+a_{21}+a_{22}\right)}^2\frac{C^2}{{\omega }^2}(50)$

$y=2\left(a_{11}^{\prime }{\overline{\rho }}_f+a_{22}{\overline{\rho }}_s\right)\left(a_{11}^{\prime }+a_{12}+a_{21}+a_{22}\right)\frac{C}{\omega }-4\left(a_{11}^{\prime }{a}_{22}-a_{12}{a}_{21}\right)\left({\overline{\rho }}_f+{\overline{\rho }}_s\right)\frac{C}{\omega }(51)$

where $sign\left(x\right)$ refers to the application of the sign of $x$, and ${\overline{\rho }}_s=<{\rho }_s>$, ${\overline{\rho }}_f=<{\rho }_f>$.

In general, for a given frequency $\omega$, 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 $\mathrm{z}$-direction can be expressed as:

$\left\{\begin{array}{c}{\mathbf{u}}_S={\mathbf{u}}_{S0}{\mathrm{e}}^{i\left(\omega t-\mathbf{k}\cdot \mathbf{x}\right)}\\ {\mathbf{U}}_S={\mathbf{U}}_{S0}{\mathrm{e}}^{i\left(\omega t-\mathbf{k}\cdot \mathbf{x}\right)}\end{array}\right.(52)$

where ${\mathbf{u}}_{S0}$ and ${\mathbf{U}}_{S0}$ refer to S-wave amplitude. After substituting the above equations into the S-wave equation, the subsequent expression can be derived:

$\left\{-{\omega }^2\left[\begin{array}{cc}〈{\rho }_s〉& 0\\ 0& 〈{\rho }_f〉\end{array}\right]+k^2\left[\begin{array}{cc}G& 0\\ 0& 0\end{array}\right]+\mathrm{i}\omega \left[\begin{array}{cc}C& -C\\ -C& C\end{array}\right]\right\}\left[\begin{array}{c}{\mathbf{u}}_{S0}\\ {\mathbf{U}}_{S0}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right](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=〈{\rho }_s〉〈{\rho }_f〉-\frac{iC}{\omega }\left(〈{\rho }_s〉+〈{\rho }_f〉\right)(55)$

$Y_2=-\left(\frac{\mathrm{i}C}{\omega }-〈{\rho }_f〉\right)G(56)$

where $X={\omega }^2/k^2$. Thus the analytical expression for S-wave velocity can be solved as:

$X=E^{\prime }+F^{\prime }i(57)$

$E^{\prime }=-\frac{{\overline{\rho }}_f^2{\overline{\rho }}_{s}G+\frac{C^2}{{\omega }^2}\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)G}{{\left({\overline{\rho }}_s{\overline{\rho }}_f\right)}^2+\frac{C^2}{{\omega }^2}{\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)}^2}(58)$

$F^{\prime }=\frac{{\overline{\rho }}_f{\overline{\rho }}_{s}G-{\overline{\rho }}_f\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)G}{{\left({\overline{\rho }}_s{\overline{\rho }}_f\right)}^2+\frac{C^2}{{\omega }^2}{\left({\overline{\rho }}_s+{\overline{\rho }}_f\right)}^2}\cdot \frac{C}{\omega }(59)$

In general, for a given frequency $\omega$, 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.

4 Finite difference numerical method for wave equation of fractured-porous media

4.1 Staggered grid difference method for the first-order wave equation

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 ${\mathrm{v}}^f=\overline{\dot{\mathrm{U}}}$ and ${\mathrm{v}}^s=\overline{\dot{\mathrm{u}}}$, the macroscopic equation for fluids and solid are rewritten as

$〈{\rho }_f〉{\dot{\mathbf{v}}}^f=\nabla \cdot 〈\mathbf{s}〉-\frac{\eta {\varphi }^2}{\kappa }\left({\mathbf{v}}^f-{\mathbf{v}}^s\right)(60)$

$〈{\rho }_f〉{\dot{\mathbf{v}}}^s=\nabla \cdot 〈\mathbf{\sigma }〉+\frac{\eta {\varphi }^2}{\kappa }\left({\mathbf{v}}^f-{\mathbf{v}}^s\right)(61)$

In the equation provided above, the brackets in $〈\cdot 〉$ 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:

$\left\{\begin{array}{c}{\dot{\mathbf{v}}}^s=\frac{1}{{\rho }_s}\nabla \cdot \mathbf{\sigma }+\frac{C}{{\rho }_s}\left({\mathbf{v}}^f-{\mathbf{v}}^s\right)\\ {\dot{\mathbf{v}}}^f=\frac{1}{{\rho }_f}\nabla \cdot \mathbf{s}-\frac{C}{{\rho }_f}\left({\mathbf{v}}^f-{\mathbf{v}}^s\right)\end{array}\right.(62)$

Here $C=\frac{\eta {\varphi }^2}{\kappa }$. By applying the time derivative to both sides of the constitutive equation, the governing equation for the fractured-porous model can be obtained:

$\left\{\begin{array}{c}{\dot{\sigma }}_{ij}=\left(a_{11}\dot{e}+a_{12}\dot{\xi }\right){\delta }_{ij}+2G\left({\dot{\epsilon }}_{ij}-\frac{1}{3}\dot{e}\cdot {\delta }_{ij}\right)\\ {\dot{s}}_{ij}=\left(a_{21}\dot{e}+a_{22}\dot{\xi }\right){\delta }_{ij}\end{array}\right.(63)$

It should be noted that the stress tensor exhibits symmetry, resulting in a total of 9 unknowns, namely, ${\sigma }_{xx}$, ${\sigma }_{xy}$, ${\sigma }_{xz}$, ${\sigma }_{yy}$, ${\sigma }_{yz}$, ${\sigma }_{zz}$, $s_{xx}$, $s_{yy}$ and $s_{zz}$, which are all associated with stress and therefore necessitate updates using constitutive relations. In the case of fluids, it satisfies the condition of $s_{xx}=s_{yy}=s_{zz}$, 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_{xx}$, $s_{yy}$ and $s_{zz}$ 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\Delta x$, $y=j\Delta y$, $z=k\Delta z$ and $t=n\Delta t$, where $\Delta x$, $\Delta y$ and $\Delta z$ denote spatial step lengths and $\Delta 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 $\mathrm{x}$-direction velocity components for solid and fluid particles respectively, with analogous formulations for the $\mathrm{y}$ and $\mathrm{z}$ direction components.

FIGURE 1

FIGURE 1. The schematic diagram of finite difference staggered grid.

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=\frac{{\theta }_1{f}_{i+1,j,k}^n-{\theta }_2{f}_{i,j,k}^n+{\theta }_3{f}_{i+2,j,k}^n-{\theta }_4{f}_{i-1,j,k}^n}{\Delta h}(64)$

$D_t{f}_{i,j,k}^{n+1/2}=\frac{f_{i,j,k}^{n+1}-f_{i,j,k}^n}{\Delta 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$. ${\theta }_i$ ($i=$ 1,2,3,4) are the coefficients. $\Delta h$ and $\Delta t$ represent the spatial and time grid sizes, respectively. For uniformly spaced grid, $\Delta x=\Delta y=\Delta z=\Delta h$, ${\theta }_1={\theta }_2=9/8$, ${\theta }_3={\theta }_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}=\frac{1}{{\rho }_s}\left(D_x{{\sigma }_{xx}}_{i+1/2,j+1/2,k}^n+D_y{{\sigma }_{xy}}_{i,j+1,k}^n+D_z{{\sigma }_{xz}}_{i,j+1/2,k+1/2}^n\right)+\frac{C}{{\rho }_s}\left({v_x^f}_{i,j+1/2,k}^{n-1/2}-{v_x^s}_{i,j+1/2,k}^{n-1/2}\right)(66)$

$D_t{v_y^s}_{i+1,j,k}^{n-1/2}=\frac{1}{{\rho }_s}\left(D_x{{\sigma }_{xy}}_{i+1,j,k}^n+D_y{{\sigma }_{yy}}_{i+1/2,j+1/2,k}^n+D_z{{\sigma }_{yz}}_{i+1/2,j,k+1/2}^n\right)+\frac{C}{{\rho }_s}\left({v_y^f}_{i+1/2,j,k}^{n-1/2}-{v_y^s}_{i+1/2,j,k}^{n-1/2}\right)(67)$

$\begin{array}{l}{D}_t{v_z^s}_{i+1/2,j+1/2,k+1/2}^{n-1/2}=\frac{1}{{\rho }_s}\left(D_x{{\sigma }_{xz}}_{i+1,j+1/2,k+1/2}^n+D_y{{\sigma }_{yz}}_{i+1/2,j+1,k}^n+D_z{{\sigma }_{xz}}_{i+1/2,j+1/2,k+1}^n\right)\\ +\frac{C}{{\rho }_s}\left({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}\right)\end{array}(68)$

$D_t{v_x^f}_{i,j+1/2,k}^{n-1/2}=\frac{1}{{\rho }_f}{D}_x{s}_{i+1/2,j+1/2,k}^n-\frac{C}{{\rho }_f}\left({v_x^f}_{i,j+1/2,k}^{n-1/2}-{v_x^s}_{i,j+1/2,k}^{n-1/2}\right)(69)$

$D_t{v_y^f}_{i+1/2,j,k}^{n-1/2}=\frac{1}{{\rho }_f}{\beta }_2{D}_y{s}_{i+1/2,j+1/2,k}^n-\frac{C}{{\rho }_f}\left({v_y^f}_{i+1/2,j,k}^{n-1/2}-{v_y^s}_{i+1/2,j,k}^{n-1/2}\right)(70)$

$D_t{v_z^f}_{i+1/2,j+1/2,k+1/2}^{n-1/2}=\frac{1}{{\rho }_f}{D}_z{s}_{i+1/2,j+1/2,k+1}^n-\frac{C}{{\rho }_f}\left({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}\right)(71)$

The finite difference equations for stress update are as follows:

$\begin{array}{l}{D}_t{{\sigma }_{xx}}_{i+1/2,j+1/2,k}^n=a_{11}^{\prime }{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}\left(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}\right)\end{array}(72)$

$\begin{array}{l}{D}_t{{\sigma }_{yy}}_{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}^{\prime }{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}\left(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}\right)\end{array}(73)$

$\begin{array}{l}{D}_t{{\sigma }_{zz}}_{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}^{\prime }{D}_z{v_z^s}_{i+1/2,j+1/2,k+1/2}^{n+1/2}\\ +a_{12}\left(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}\right)\end{array}(74)$

$D_t{{\sigma }_{xy}}_{i,j,k}^n=G\left(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}\right)(75)$

$D_t{{\sigma }_{xz}}_{i,j+1/2,k+1/2}^n=G\left(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}\right)(76)$

$D_t{{\sigma }_{yz}}_{i+1/2,j,k+1/2}^n=G\left(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}\right)(77)$

$\begin{array}{l}{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}\left(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}\right)\end{array}(78)$

where $a_{11}^{*}=a_{11}-2G/3$, $a_{11}^{\prime }=a_{11}+4G/3$.The notations ${{\sigma }_{xx}}_{i,j,k}^n$, ${{\sigma }_{yy}}_{i,j,k}^n$, ${{\sigma }_{zz}}_{i,j,k}^n$, ${{\sigma }_{xy}}_{i,j,k}^n$, ${{\sigma }_{yz}}_{i,j,k}^n$, ${{\sigma }_{xz}}_{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 ${\sigma }_{xx}$, ${\sigma }_{yy}$ and ${\sigma }_{zz}$ in three directions of the solid phase, as well as the fluid stress $s$. The source function is based on Gaussian curve:

$f\left(t\right)=-2\xi \left(1-2\xi T^2\right)e^{-\xi T^2}(79)$

where $\xi =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:

$\left({\displaystyle \prod _{j=1}^2}\left(\cos {\alpha }_j\frac{\partial }{\partial t}-c\frac{\partial }{\partial x}\right)\right)u=0(80)$

where ${\alpha }_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 $\pm {\alpha }_j$. Following discretization on the difference grid, the differential operator form of the absorbing boundary condition is presented as follows:

$B\left(E_x,E_t^{-1}\right)={\displaystyle \prod _{j=1}^2}\left\{{\beta }_j\left(\frac{I-E_t^{-1}}{\Delta t}\right)\left[\left(1-a\right)I+a{E}_x\right]-c\left(\frac{E_x-I}{\Delta x}\right)\left[\left(1-b\right)I+b{E}_t^{-1}\right]\right\}(81)$

where ${\beta }_j=\cos {\alpha }_j$ represents a positive number, and $\Delta x$ signifies the spatial step size in the $\mathrm{x}$ direction on the boundary grid, and $I$ denotes the identity operator, and $E_{x}I=E_x$, $aI=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\Delta t$ and at the point $x=i\Delta x$ and $i=0$ (i.e., the left boundary in the $\mathrm{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=-\frac{1}{B_9}\left(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\right)(82)$

Given the initial conditions, the wavefield’s physical variables at each node are known at time $\mathbf{t}=0$. Since only the wave field values of two time steps are available, first-order absorbing boundary conditions can be employed to address the boundary reflection at time $t=\Delta t$. Following simplification, the solution can be obtained

$f_{i,j,k}^n=-\frac{1}{A_1^4}\left(A_1^1{f}_{i+1,j,k}^{n-1}+A_1^2{f}_{i,j,k}^{n-1}+A_1^3{f}_{i+1,j,k}^n\right)(83)$

where each coefficient refers to as follows:

$\begin{array}{c}{p}_j=c\Delta t/{\beta }_j\Delta x,A_j^1=-\left(p_{j}b+a\right),A_j^2=p_{j}b+a-1,A_j^3=a-p_j\left(1-b\right),\hfill \\ A_j^4=\left(1-a\right)+p_j\left(1-b\right),B_1=A_1^1{A}_2^1,B_2=A_1^1{A}_2^2+A_1^2{A}_2^1,\hfill \\ B_3=A_1^2{A}_2^2,B_4=A_1^1{A}_2^3+A_1^3{A}_2^1,B_5=A_1^1{A}_2^4+A_1^2{A}_2^3+A_1^3{A}_2^2+A_1^4{A}_2^1,\hfill \\ B_6=A_1^2{A}_2^4+A_1^4{A}_2^2,B_7=A_1^3{A}_2^3\text{\hspace{0.17em}and\hspace{0.17em}}{B}_8=A_1^3{A}_2^4+A_1^4{A}_2^3.\hfill \end{array}$

Therefore, the calculation of the wavefield value at the node $x=i\Delta x$ and $i=0$ (i.e., the left boundary in the $\mathrm{x}$ direction) with first-order absorbing boundary requires three layers of data: the wavefield value $F_1=f_{i=0,k}^{n=0}$ of the previous time step’s boundary layer (i.e., $t=0$ and $x=0$); the wavefield value $F_2=f_{i=1,k}^{n=0}$ of the first layer inside the previous time step’s boundary (i.e., $t=0$ and $x=\Delta x$); and the wavefield value $F_3=f_{i=1,k}^{n=1}$ of the first layer inside the current time step’s boundary (i.e., $t=\Delta t$ and $x=\Delta x$).

The proposed method can be employed to implement absorbing boundary processing for other physical variables on boundaries within the wavefield region. Due to the directional nature of wave velocity on different boundaries, the wave velocity $c$ in the absorbing boundary operator can have either a positive or negative sign. For instance, on the left and upper boundaries, reflected waves propagate in the positive direction along the coordinate axis, thus requiring a positive sign for $c$; otherwise, it necessitates a negative sign.

4.2 Comparison of wave velocity and attenuation between finite difference and plane wave analysis

To validate the effectiveness of the finite difference method for the wave equation in fractured-porous network, the core sample saturated oil shown in Table 1 is selected. Based on the finite difference in the staggered grid, the numerical solutions of the wave equation is obtained to generate wavefield snapshots, followed by the analysis of wavefield velocity and attenuation.

TABLE 1

TABLE 1. Parameter table of porous media.

The parameters of the differential grid are presented in Table 2.

TABLE 2

TABLE 2. Parameter table of finite difference simulation.

The calculation space range is $250\mathrm{m}\times 250\mathrm{m}$, and the wave propagation time is measured to be 0.055 s. The seismic source is positioned at the center of the calculation area. Figure 2 displays the calculated snapshot of solid particle velocity, with spatial coordinates representing its extent. The analysis findings reveal distinct fast P-wave within this snapshot. Moreover, due to waves propagating in various directions, a phase reversal in velocity occurs at the midpoint of the wavefield. Additionally, faint outlines of S-wave are discernible but exhibit weaker amplitude; however, slow P-wave are challenging to observe due to significant dissipation.

FIGURE 2

FIGURE 2. The wave field snapshot of oil-bearing fractured-porous media (The velocity of solid particle).

Wavefield information is computed using parameters corresponding to various aspect ratios, and the obtained results are compared with the dispersion and attenuation outcomes derived from plane wave analysis. The aspect ratio “a” ranges from 0.1 to 1.0 in increments of 0.1, resulting in a total of 10 values. In this study, the proposed method is employed to calculate the dry skeleton modulus and permeability (under steady flow conditions), which are subsequently substituted into both the numerical solution of the wave equation and the analytical solution of the plane wave analysis. The values of dry skeleton modulus and permeability for different aspect ratios are shown in Figure 3.

FIGURE 3

FIGURE 3. Dry skeleton modulus (K_b) and permeability (κ) corresponding to different aspect ratios (a).

The procedure for computing the dispersion and attenuation of fast P-wave based on seismic data is outlined as follows: Firstly, in the calculation process of the wave field snapshot, two receivers are strategically positioned to measure particle velocity and stress amplitude at distinct locations. Let $A\left(x_1\right)$ denote the recorded amplitude value at position $x_1$, while $A\left(x_2\right)$ represents the corresponding value at position $x_2$. The separation between these receivers is denoted as $\Delta s=\left|x_2-x_1\right|$. Secondly, obtain multiple sets of waveform graphs, correlate the waveforms at these two positions, determine the time difference $\Delta t$, and thereby derive the velocity value $v=\Delta \mathrm{s}/\Delta \mathrm{t}$. Additionally, the quality factor $Q=\frac{\pi f\left(x_2-x_1\right)}{v}{\left[\ln \left(\frac{A\left(x_1\right)}{A\left(x_2\right)}\right)\right]}^{-1}$ is determined by employing the amplitude attenuation method (Gong, et al., 2009), which calculates by the ratio of amplitudes at two different distances (or different times), where $f$ represents the source frequency. The two receivers in this example are positioned horizontally at the same elevation as the seismic source. The first receiver is located 25 m to the left of the source, with a separation distance of 10 grid cells or 25 m between them.

The theoretical solutions of dispersion and attenuation, along with the numerical results obtained from the finite difference method, are compared in Figures 4, 5. It can be observed from these figures that the fast P-wave velocity calculated using the proposed finite difference scheme exhibits excellent agreement with the theoretical solution based on plane wave analysis. Additionally, a consistent trend is observed for the inverse quality factor. These computational findings provide validation for the efficacy of the proposed method.

FIGURE 4

FIGURE 4. Comparison between theoretical solution and numerical solution of fast P-wave velocity under different aspect ratios.

FIGURE 5

FIGURE 5. Comparison between theoretical solution and numerical solution of fast P-wave attenuation under different aspect ratios.

4.3 Influence of confining pressure on P-wave velocity

The confining pressure exerted on the porous media has a direct impact on $K_b$, thereby influencing the values of wave velocity. For quantitatively investigating the influence of confining pressure on wave velocity, the rock physics parameters used in the calculation are the same as those in Table 1, in which the seismic wave frequency is fixed at 100 Hz, the aspect ratio is 0.15, the confining pressure changes from 10³ Pa to 10⁹ Pa. The seismic velocities can be obtained by substituting the corresponding $K_b$ under different confining pressures, and the results show in Figures 6, 7.

FIGURE 6

FIGURE 6. Variation of fast P wave velocity with confining pressure.

FIGURE 7

FIGURE 7. Variation of slow P wave velocity with confining pressure.

The calculation results demonstrate that the wave velocity increases with the applied confining pressure, owing to the concurrent rise in dry skeleton modulus $K_b$. Additionally, a similar relationship is observed between wave velocity and confining pressure, as well as between dry skeleton modulus $K_b$ and confining pressure, with their respective curves exhibiting analogous trends. Moreover, the changing trend of slow P-wave is similar to that of fast P-wave. The crucial point to emphasize here is that the frequency of seismic waves exerts minimal influence on the relationship curve between confining pressure $p_c$ and wave velocity, particularly within the relatively low-frequency range, as exemplified by this case where the frequency is set at 100 Hz. If any value within the frequency range of 10⁰–10⁴ Hz is considered, it can be observed that the curves remain the same, which can be inferred from the dispersion curve maintaining a horizontal shape in the low frequency range.

Based on this case study, it can be inferred that the velocities of fast P-waves and slow P-waves exhibit a gradual increase as the confining pressure increases. This observed trend remains consistent across different frequencies of seismic waves, particularly at relatively low frequencies (<10⁴ Hz).

5 Permeability prediction based on wave equation of fractured-porous media

The propagation of waves in porous media containing fluids gives rise to phenomena such as dispersion and attenuation, wherein both the wave velocity and amplitude exhibit frequency-dependent changes. The dispersion and attenuation curves of the wavefield are influenced by rock physical properties including porosity, permeability, fluid properties, and solid skeleton bulk modulus. Among these factors, we specifically investigated the sensitivity of dispersion and attenuation curves to variations in permeability. If alterations in permeability can be discerned on dispersion/attenuation curves, it becomes feasible to estimate these changes through observation and calculation of these curves or even determination of numerical values for permeability using a template method.

5.1 Permeability prediction of permeability based on dispersion and attenuation

According to the wave equation and plane wave analysis, expressions for dispersion/attenuation of the wave field has been yielded. The dispersion and attenuation curves of the wave field are then calculated based on rock physics parameters under different permeability. A set of solid skeleton parameters and fluid parameters is selected as presented in Table 3 (Johnson, 2001; Lo, et al., 2005).

TABLE 3

TABLE 3. Rock physical parameters.

The calculated dispersion/attenuation curves are depicted in Figures 8, 9.

FIGURE 8

FIGURE 8. Variation of fast P-wave dispersion curves with permeability.

FIGURE 9

FIGURE 9. Variation of fast P-wave attenuation curves with permeability.

In the example, various porous network structures of different scales ($4\times 4\times 4$, $5\times 5\times 5$, $6\times 6\times 6$, $7\times 7\times 7$ and $8\times 8\times 8$) were employed within a unit cube with a side length of $l$. As the density of fractures and pores increased, the connectivity between pores became denser, leading to a gradual reduction in permeability (from 3.5888 to 0.22417 Darcy). Different permeability corresponded to distinct shapes of fast P-wave dispersion curves for diverse-scale porous network structures. Figure 8 demonstrates that as permeability changes, the dispersion curves shift, with a discrepancy in P-wave velocity reaching up to 40 m/s at identical frequency in this case study. Figure 9 illustrates how fast P-waves’ attenuation curve varies with permeability, indicating an observable shift in the position of attenuation peaks.

In conclusion, the permeability can be modified by adjusting the porous network structure (pore density), based on rock physical parameters such as porosity, fluid parameters, and solid skeleton bulk modulus. Different permeability conditions result in noticeable and regular changes in the wavefield curves of dispersion and attenuation. Therefore, it is possible to calculate the variation pattern of dispersion and attenuation with respect to permeability when rock physical parameters are known. Consequently, changes in permeability can be inferred by measuring the movement of dispersion/attenuation curves.

5.2 Influence of fracture aspect ratio on permeability

By utilizing the curves of permeability and bulk modulus for different fracture aspect ratios within a fully saturated elastic model, dispersion curves for P- and S-wave velocities, along with prediction for attenuation, are successfully derived. Subsequently, by obtaining the permeability and P/S-wave velocities in porous media with varying fracture aspect ratios, the plot can be drawn, illustrating the relationship between P/S-wave velocity and permeability.

The graphical analysis of seismic attributes is presented in Figure 10, illustrating the calculation results. For sandstone, the calculated parameters are as follows: density of 2,650 kg/m³, bulk modulus of 37 GPa, shear modulus of 44 GPa, Poisson’s ratio of 0.08, and Young’s modulus of 94.5 GPa. When determining wave velocity, a fixed frequency of 30 Hz is employed.

FIGURE 10

FIGURE 10. Relationship between P-wave velocity and permeability of 3D water-bearing fracture network.

The permeability gradually decreases as the aspect ratio of fractures decreases under a given porosity condition, as observed from Figure 10. The data points on the graph move from bottom right to top left, and with approaching zero aspect ratio, there is a rapid increase in the distance between data points.

From the calculation results, it is evident that there exists a significant variation in the velocity ratio of P- and S-wave under different fracture/pore aspect ratios, indicating a substantial impact of aspect ratio on this velocity ratio. The influence of aspect ratio on the velocity ratio arises from two aspects: 1) the effect of varying aspect ratios on the bulk modulus of rock skeleton; 2) the effect of varying aspect ratios on permeability. These two influences are closely intertwined, and when assessing the feasibility of predicting permeability using seismic attributes, both aspects need to be simultaneously considered rather than isolating one-sided effects.

5.3 Identification of permeability change based on velocity ratio-impedance template

As aspect ratio have effects on both bulk modulus and permeability, their contributions to the velocity ratio are intertwined, making it impossible to separately examine the influence of permeability on wave velocity. To further analyze the relationship between wave velocity ratio and permeability, the fracture aspect ratio parameter is held constant while varying only the numerical value of permeability, then calculate data for the velocity ratio-impedance template.

To facilitate comparison with previous examples, the aspect ratio is standardized to 1 and the corresponding permeability values from those examples (1.2 $\times {10}^{-5}$ D, 1.8 $\times {10}^{-4}$ D, 0.036D, 0.2D, 0.47D) are adopted. Figure 11 illustrates the velocity ratio-impedance template.

FIGURE 11

FIGURE 11. Variation of velocity ratio with impedance in 3D water-bearing fracture network.

From Figure 11, it is evident that when the aspect ratio remains constant (in this case, 1) and only the permeability varies, there is a significant reduction in the range of variation observed in the velocity ratio. Simultaneously, there is also a substantial decrease in the range of impedance changes. Within the same range of permeability variations as considered in previous calculations, discerning changes in permeability solely from the velocity ratio becomes nearly impossible. Hence, it can be concluded that alterations in bulk modulus resulting from variations in aspect ratio are primarily responsible for considerable changes observed in velocity ratio, while modifications in permeability due to changes in aspect ratio have relatively minor impacts on velocity ratio.

In order to ascertain the extent of permeability changes that can be discerned in the velocity ratio-impedance template, the variation pattern of velocity ratio as permeability ranged from 1 millidarcy to 100 darcies has been investigated. Figure 12 illustrates that when permeability exhibits a wide range of variability (0.001–100 Darcy), there is a significant amplification in the magnitude of impedance (from 0.32 to 130,000 kg/s.m²). Hence, during substantial fluctuations in permeability, corresponding variations in impedance become more pronounced. However, it should be noted that the numerical value for velocity ratio change remains exceedingly small (0.00041), rendering it inadequate for independently distinguishing alterations in permeability.

FIGURE 12

FIGURE 12. Variation of velocity ratio with impedance in 3D water-bearing fracture network (the permeability changing from 0.001 to 100 Darcy and the aspect ratio of 1).

The following analysis investigates the patterns observed under different aspect ratios, replicating the aforementioned calculations for aspect ratios of 0.5 and 0.05. Figures 13, 14 demonstrate that as the aspect ratio decreases, there is a significant reduction in the range of variation in wave impedance, decreasing from 130,000 kg/s.m² at $a=1$ to 11 kg/s.m² at $a=0.05$. However, concurrently with the decrease in aspect ratio, there is a continuous increase in the wave velocity ratio, escalating from 0.00041 to 0.0023.

FIGURE 13

FIGURE 13. Variation of velocity ratio with impedance in 3D water-bearing fracture network (the permeability changing from 0.001 to 100 Darcy and the aspect ratio of 0.5).

FIGURE 14

FIGURE 14. Variation of velocity ratio with impedance in 3D water-bearing fracture network (the permeability changing from 0.001 to 100 Darcy and the aspect ratio of 0.05).

From the analysis of the calculation results above, it is evident that when considering rocks with identical physical parameters and varying only in permeability, there is minimal variations observed in the velocity ratio of P- to S-wave at the frequency of 30 Hz. Consequently, identifying changes in permeability becomes exceedingly challenging. In scenarios where there exists a significant alteration in permeability (e.g., from 1 millidarcy to 100 Darcy), porous networks with the large aspect ratio exhibit a relatively extensive range of impedance variation. However, even within low porosity and low permeability rocks, the dispersion degree of data points on the velocity ratio-impedance template remains insufficient for isolating the impact of permeability changes, thereby presenting substantial challenges.

5.4 Identification of permeability change based on attenuation-young’s modulus

The variation patterns of P-wave attenuation and Young’s modulus in relation to different permeability have been carried out, revealing a significant range of relative changes in attenuation. Figure 15 illustrates the variations in P-wave attenuation-Young’s modulus resulting from alterations in permeability within water-bearing porous sandstone.

FIGURE 15

FIGURE 15. Template of P-wave velocity attenuation-Young’s modulus in 3D water-bearing fracture network (The permeability changing from 0.001 to 100 Darcy and the aspect ratio of 1.

The following conclusions can be inferred from the attenuation-Young’s modulus template:

(1) The relative change in P-wave attenuation is quite noticeable. As the permeability increases, the attenuation value gradually increases, exhibiting an approximately linear growth trend.

(2) With the increase in permeability, Young’s modulus also gradually increases, and it grows more rapidly at higher permeability.

(3) If accurate data on P-wave attenuation can be obtained, the velocity attenuation-Young’s modulus template can be utilized to analyze and identify changes in permeability within reservoir when permeability changes significantly in a wide range.

(4) In cases of low porosity and low permeability where there is limited variation in permeability, the distribution range of data points on the velocity attenuation-Young’s modulus template is also narrow, posing challenges for identifying changes in permeability.

5.5 Sensitivity analysis of seismic parameters changing with permeability

In order to identify seismic attributes that exhibit high sensitivity to variations in permeability, the comprehensive sensitivity analysis on various seismic attribute parameters has been conducted. During the calculations, five permeability values were uniformly selected ranging from 0.001 Darcy to 100 Darcy, while keeping other rock physics parameters constant. The individual sensitivities of the velocity ratio of P- to S-wave, P-wave impedance, wave velocity, shear modulus, Young’s modulus, and attenuation of P- and S-wave velocities towards changes in permeability were separately evaluated. The parameters investigated in our study have distinct physical meanings in the context of fractured-porous media. The velocity ratio of P- to S-wave (V_p/V_s) serves as an indicator of fluid-filled zone. P-wave impedance (Z) and wave velocities (V_p and V_s) provide insights into rock density and elastic properties. Shear modulus (μ) and Young’s modulus (E) are mechanical properties that describe rock resistance to deformation, with lower values indicating increased deformability under stress. Lastly, attenuation parameters (Q_p and Q_s) quantify the energy loss of seismic waves, with higher permeability associated with greater attenuation, crucial for detecting fluid-saturated zones. These parameters collectively help in characterizing subsurface formations, fluid dynamics, and rock behavior in fractured-porous media, aiding in hydrocarbon exploration and reservoir assessment.

Herein, the sensitivity $s$ is defined as $s\left(x\right)=\left(x_{\max }-x_{\min }\right)/\overline{x}$ , $x_{\max }$ and $x_{\min }$ represent the maximum and minimum calculated values of the seismic attributes respectively, while $\overline{x}$ denotes their average value. Figures 16–18 present comparative charts illustrating the sensitivity of each attribute for sandstone containing water, oil or gas.

FIGURE 16

FIGURE 16. Sensitivity comparison of seismic attributes changing with permeability in 3D water-bearing fracture network (The permeability changing from 0.001 to 100 Darcy, the aspect ratio of 1 and the frequency of 30 Hz).

FIGURE 17

FIGURE 17. Sensitivity comparison of seismic attributes changing with permeability in 3D oil-bearing fracture network (The permeability changing from 0.001 to 100 Darcy, the aspect ratio of 1 and the frequency of 30 Hz).

FIGURE 18

FIGURE 18. Sensitivity comparison of seismic attributes changing with permeability in 3D gas-bearing fracture network (The permeability changing from 0.001 to 100 Darcy, the aspect ratio of 1 and the frequency of 30 Hz).

The comparative analysis reveals that the sensitivity of seismic attributes, such as velocity ratio, P-wave impedance, wave velocity, shear modulus, and Young’s modulus, to changes in permeability is relatively low. Conversely, the attenuation of P-wave and S-wave velocities exhibits a comparatively high sensitivity. The additional calculation and analysis for the aspect ratio of 0.5 are conducted, finding that the sensitivity of seismic attributes to changes in permeability remains consistent with the case where the aspect ratio is 1.

6 Conclusion

In this paper, an improved wave equation of fractured-porous media is proposed. Through this research, the significant influence of fracture aspect ratio and confining pressure on wave properties is unveiled, and a novel method for permeability identification based on the curves of velocity-fracture parameters and velocity-permeability is proposed, while numerical simulation methods are refined. The research findings demonstrate that fracture aspect ratio and confining pressure exert a substantial impact on wave properties such as velocity and attenuation; there exists a correlation between fracture parameters and permeability, which can be utilized to predict permeability by utilizing the curves of either velocity-fracture parameters or velocity-permeability. These discoveries offer fresh insights into the wave behavior of fractured-porous media, providing a theoretical foundation and innovative approaches for permeability prediction. Moreover, the methodologies and concepts presented in this study can also be applied to seismic rock physics research in various fields, thereby further advancing related areas’ development. Nevertheless, certain limitations persist in this study including the complexity of fluid flow within porous media as well as challenges associated with accurately describing the actual motion state of pore fluid. In future research, these challenges will continue to be addressed to enhance models and methods, enhancing accuracy and applicability in permeability identification.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

WZ: Data curation, Funding acquisition, Investigation, Writing–original draft. WS: Conceptualization, Formal Analysis, Methodology, Writing–original draft. JG: Project administration, Supervision, Validation, Writing–review and editing. RH: Investigation, Visualization, Writing–review and editing.

Funding

The authors declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Natural Science Foundation of China (Grant Nos 42074144 and 41874137).

Conflict of interest

Authors WZ, JG, and RH were employed by PetroChina.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The authors declare that this study received funding from CNPC Science Research and Technology Development Project under Grant 2021DJ3701. The funder had the following involvement in the study: study design, data collection and decision to publish.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Bai, B., Jiang, S., Liu, L., Li, X., and Wu, H. (2021). The transport of silica powders and lead ions under unsteady flow and variable injection concentrations. Powder Technol. 387, 22–30. doi:10.1016/j.powtec.2021.04.014

CrossRef Full Text | Google Scholar

Biot, M. A. (1956a). Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range. J. Acoust. Soc. Am. 28, 168–178. doi:10.1121/1.1908239

CrossRef Full Text | Google Scholar

Biot, M. A. (1956b). Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range. J. Acoust. Soc. Am. 28, 179–191. doi:10.1121/1.1908241

CrossRef Full Text | Google Scholar

Chapman, M., Zatsepin, S. V., and Crampin, S. (2002). Derivation of a microstructural poroelastic model. Geophys. J. Int. 151, 427–451. doi:10.1046/j.1365-246X.2002.01769.x

CrossRef Full Text | Google Scholar

Chapman, S., Tisato, N., Quintal, B., and Holliger, K. (2016). Seismic attenuation in partially saturated Berea sandstone submitted to a range of confining pressures. J. Geophys. Res. Solid Earth 121, 1664–1676. doi:10.1002/2015jb012575

CrossRef Full Text | Google Scholar

Chichinina, T., Obolentseva, I. R., Ronquillo-Jarillo, G., Sabinin, V. I, Gik, L. D, Bobrov, B. A, et al. (2007). Attenuation anisotropy of P- and S- waves: theory and laboratory experiment. J. Seismic Explor. 16, 235–264. doi:10.1093/petrology/egm066

CrossRef Full Text | Google Scholar

Chichinina, T. I., Obolentseva, I. R., and Ronquillo-Jarillo, G. (2009). Anisotropy of seismic attenuation in fractured media: theory and ultrasonic experiment. Transp. Porous Media 79, 1–14. doi:10.1007/s11242-008-9233-9

CrossRef Full Text | Google Scholar

Du, Q.-Z., Wang, X.-M., Ba, J., and Zhang, Q. (2011). An equivalent medium model for wave simulation in fractured porous rocks. Geophys. Prospect. 60 (5), 940–956. doi:10.1111/j.1365-2478.2011.01027.x

CrossRef Full Text | Google Scholar

Gong, T. J., Sun, C. Y., and Peng, H. C. (2009). Comparison of several computational methods of quality factor. Prog. Explor. Geophys. 32 (4), 252–256.

Google Scholar

Gray, W. G., and Lee, P. C. Y. (1977). On the theorems for local volume averaging of multiphase systems. Int. J. Multiph. Flow 3, 333–340. doi:10.1016/0301-9322(77)90013-1

CrossRef Full Text | Google Scholar

Guo, J., Rubino, J. G., Barbosa, N. D., Glubokovskikh, S., and Gurevich, B. (2018). Seismic dispersion and attenuation in saturated porous rocks with aligned fractures of finite thickness: theory and numerical simulations — Part 1: P-wave perpendicular to the fracture plane. Geophysis 83, WA49-WA62. doi:10.1190/geo2017-0065.1

CrossRef Full Text | Google Scholar

Gurevich, B., Brajanovski, M., Galvin, R. J., Müller, T. M., and Toms-Stewart, J. (2009). P-wave dispersion and attenuation in fractured and porous reservoirs – poroelasticity approach. Geophys. Prospect. 57, 225–237. doi:10.1111/j.1365-2478.2009.00785.x

CrossRef Full Text | Google Scholar

Higdon, R. L. (1986). Absorbing boundary conditions for difference approximations to the multi-dimensional wave equation. Math. Comput. 47, 437–459. doi:10.2307/2008166

CrossRef Full Text | Google Scholar

Higdon, R. L. (1987). Numerical absorbing boundary conditions for the wave equation. Math. Comput. 49, 65–90. doi:10.1090/s0025-5718-1987-0890254-1

CrossRef Full Text | Google Scholar

Johnson, D. L. (2001). Theory of frequency dependent acoustics in patchy-saturated porous media. J. Acoust. Soc. Am. 110, 682–694. doi:10.1121/1.1381021

CrossRef Full Text | Google Scholar

Landau, L. D., and Lifshitz, E. M. (1987). Fluid mechanics. 2nd Edition. Pergamon.

Google Scholar

Lissa, S., Barbosa, N. D., Rubino, J. G., and Quintal, B. (2019). Seismic attenuation and dispersion in poroelastic media with fractures of variable aperture distributions. Solid earth. 10, 1321–1336. doi:10.5194/se-2019-18

CrossRef Full Text | Google Scholar

Lo, W.-C., Sposito, G., and Majer, E. (2005). Wave propagation through elastic porous media containing two immiscible fluids. Water Resour. Res. 41, W02025. doi:10.1029/2004wr003162

CrossRef Full Text | Google Scholar

Mavko, G. M., and Nur, A. (1978). The effect of nonelliptical cracks on the compressibility of rocks. J. Geophys. Res. Solid Earth 83 (B9), 4459–4468. doi:10.1029/JB083iB09p04459

CrossRef Full Text | Google Scholar

Schoenberg, M. (1980). Elastic wave behavior across linear slip interfaces. J. Acoust. Soc. Am. 68, 1516–1521. doi:10.1121/1.385077

CrossRef Full Text | Google Scholar

Schoenberg, M. (1983). REFLECTION OF ELASTIC WAVES FROM PERIODICALLY STRATIFIED MEDIA WITH INTERFACIAL SLIP. Geophys. Prospect. 31, 265–292. doi:10.1111/j.1365-2478.1983.tb01054.x

CrossRef Full Text | Google Scholar

Schoenberg, M., and Douma, J. (1988). ELASTIC WAVE PROPAGATION IN MEDIA WITH PARALLEL FRACTURES AND ALIGNED CRACKS1. Geophys. Prospect. 36, 571–590. doi:10.1111/j.1365-2478.1988.tb02181.x

CrossRef Full Text | Google Scholar

Schoenberg, M., and Sayers, C. M. (1995). Seismic anisotropy of fractured rock. GEOPHYSICS 60, 204–211. doi:10.1190/1.1443748

CrossRef Full Text | Google Scholar

Slattery, J. C. (1967). Flow of viscoelastic fluids through porous media. Aiche J. 13, 1066–1071. doi:10.1002/aic.690130606

CrossRef Full Text | Google Scholar

Song, Y., Hu, H., and Han, B. (2020). Seismic attenuation and dispersion in a cracked porous medium: an effective medium model based on poroelastic linear slip conditions. Mech. Mater. 140, 103229. doi:10.1016/j.mechmat.2019.103229

CrossRef Full Text | Google Scholar

Sun, W. T. (2009). Finite difference method for elastic wave equation. Beijing, China: Tsinghua University Press.

Google Scholar

Tang, X. (2011). A unified theory for elastic wave propagation through porous media containing cracks—an extension of Biot’s poroelastic wave theory. Sci. China Earth Sci. 54, 1441–1452. doi:10.1007/s11430-011-4245-7

CrossRef Full Text | Google Scholar

Tang, X., Chen, X., and Xu, X. (2012). A cracked porous medium elastic wave theory and its application to interpreting acoustic data from tight formations. GEOPHYSICS 77, D245–D252. doi:10.1190/geo2012-0091.1

CrossRef Full Text | Google Scholar

Tuncay, K., and Corapcioglu, M. Y. (1996). Body waves in poroelastic media saturated by two immiscible fluids. J. Geophys. Res. Solid Earth 101 (B11), 25149–25159. doi:10.1029/96jb02297

CrossRef Full Text | Google Scholar

Tuncay, K., and Corapcioglu, M. Y. (1997). Wave propagation in poroelastic media saturated by two fluids. J. Appl. Mech. 64 (2), 313–320. doi:10.1115/1.2787309

CrossRef Full Text | Google Scholar

Wang, K., Peng, S., Lu, Y., and Cui, X. (2022). Wavefield simulation of fractured porous media and propagation characteristics analysis. Geophys. Prospect. 70, 886–903. doi:10.1111/1365-2478.13198

CrossRef Full Text | Google Scholar

Wei, L. L., Gan, L. D., Xiong, F. S., Sun, W. T., Ding, Q., Yang, H., et al. (2021). Attenuation and dispersion characteristics of P-waves based on the three-dimensional pore network model. Chin. J. Geophys. 64 (12), 4618–4628. doi:10.6038/cjg2021P0496

CrossRef Full Text | Google Scholar

Whitaker, S. (1966). The equations of motion in porous media. Chem. Eng. Sci. 21, 291–300. doi:10.1016/0009-2509(66)85020-0

CrossRef Full Text | Google Scholar

Whitaker, S. (1967). Diffusion and dispersion in porous media. AIChE J. 13, 420–427. doi:10.1002/aic.690130308

CrossRef Full Text | Google Scholar

Xiong, F. S., Gan, L. D., and Sun, W. T. (2021). Characterization of reservoir permeability and analysis of influencing factors in fracture-pore media. Chin. J. Geophys. 64 (1), 279–288. doi:10.6038/cjg2021N0175

CrossRef Full Text | Google Scholar