Skip to main content

ORIGINAL RESEARCH article

Front. Earth Sci., 29 December 2023
Sec. Solid Earth Geophysics
This article is part of the Research Topic Advances of New Technologies in Seismic Exploration View all 22 articles

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

  • 1Research Institute of Petroleum Exploration and Development-Northwest, PetroChina, Lanzhou, China
  • 2School of Aerospace Engineering, Tsinghua University, Beijing, China

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:

ρft+ρfU˙=0(1)

The equation of fluid momentum conservation:

s=ρfUt+p(2)

Constitutive relationship of fluid mechanics:

s=ηU˙+U˙T23U˙I(3)

where ρf represents the fluid density, s denotes the fluid stress tensor, and p signifies the fluid pressure, U is the displacement of fluid, U˙ represents the flow velocity, η stands for the fluid viscosity coefficient, and DU˙/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:

ρfDU˙Dt=p+η2U˙+13U˙(4)

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

U˙x,y,t=U˙0x,y,z+U˙1x,y,z,t(5)

Under low-speed flow conditions, the compressibility of pore fluid can be neglected when subjected to elastic waves (denoted as 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, DU˙/Dt=U˙/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:

QU=iπR2ρfcPUcosωL/cPDsinωL/c2J1KRKRJ0KR1QD=iπR2ρfcPUPDcosωL/csinωL/c2J1KRKRJ0KR1(6)

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

QUe=2a31+a2iπR2ρfcPUcosωL/cPDsinωL/c2J1KRKRJ0KR1QDe=2a31+a2iπR2ρfcPUPDcosωL/csinωL/c2J1KRKRJ0KR1(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 PU and PD at both ends, major radius R and aspect ratio α of the elliptical cross-section, characteristic frequency ω of non-steady flow field, and fluid density ρf. In Equation 7, J0 and J1 represent zero- and first-order Bessel functions respectively, K=iωρf/η, 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 R1 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 Ω, filled with fluid in the porous media, is considered. Its position x represents the “center” or “centroid”. In the fluid-filled porous media of the 3D elliptical microtubule network, the characteristic length of region Ω is denoted as Hh, the micro and macro characteristic length are represented by h and H, respectively. These three lengths are governed by the following equation.

hHhH(8)

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

The volume of the unit cell Ω is denoted as V. The pore space comprises two components: the volume Vs of the solid skeleton and the volume Vf of fluid in porous media (i.e., the volume of the pore space). Thus, V=Vs+Vf. Let ψf represent a physical quantity associated with fluid in porous media, with its value outside of the fluid defined as 0. The volume averaging method is applied to ψf within the entire region Ω, as illustrated below.

ψf=1VΩψfxdV(9)

where ψf is characterized by the smoothness with respect to the centroid x of region Ω. Another variable related to the volume average is the volume-averaging eigenvalue, defined as follows:

ψ¯f=1VfΩψfxdV(10)

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

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

ΩψfdV=ΩψfdV+ΣψfndS(11)

The volume average term including time derivative is as follows:

ΩtψfdV=tΩψfdVΣψfwndS(12)

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

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.

ρs2ut2=σ(13)

where ρs denotes the mass density of solid skeleton, and u represents the displacement of solid skeleton, t refers to time and σ signifies solid stress tensor.

The microscopic constitutive equation of solid is:

σ=KseI+Gu+uT23eI(14)

where Ks stands for the bulk modulus of solid particle, G refers to the shear modulus of solid skeleton, e denotes the first strain invariant, and superscript T indicates transpose.

The previous sections present the mass conservation equation, the momentum conservation equation and the constitutive equation of fluid. The equation of fluid state is expressed in the following form:

1Kfdpdt=1ρfdρfdt(15)

where p represents fluid pressure and Kf signifies the fluid bulk modulus. By combining the continuity equation and the state equation of fluid, the expression of pressure term is obtained:

p=Kfu(16)

Assuming the presence of a porous media comprising a solid skeleton and a fluid, with only one existing interface denoted as Σ, which represents the contact between the fluid and the solid skeleton. The boundary conditions employed to describe this interface include the no-slip condition and fluid-solid equilibrium condition of normal force, expressed as follows:

U˙n=u˙npn+sn=σn(17)

where 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:

2t2ρfU=p+s+p¯ϕηϕ2κU˙¯u˙¯(18)

where κ represents the permeability of porous media. Since fluid does not undergo shear deformation, the fluid stress is equivalent to the normal stress s. By volume-averaging the solid momentum conservation equation and integrating it with the fluid-solid equilibrium condition at boundaries, it is obtained as follows:

2t2ρsu=σp¯ϕ+ηϕ2κU˙¯u˙¯(19)

The stress σ of solid skeleton can be decomposed into two components: normal stress σs and shear stress σe. In terms of the normal stress component of solid skeleton, Tuncay and Corapcioglu (1996), Tuncay and Corapcioglu (1997) investigated the constitutive relationship associated with volume change (solely induced by compressive deformation), and derived the equations of solid normal stress and fluid stress:

σs=1φp¯s=a11u¯+a12U¯(20)
s=φp¯=a21u¯+a22U¯(21)

where ps represents the applied confining pressure on porous media, and the respective equations for each constant in this expression are provided below:

a11=1φ2KfφL11φ2KfL11φKsKb(22)
a12=1φKfL1(23)
a21=φKsL2+φKb1φL2(24)
a22=φ2Ks1φL2(25)
L1=1φKfφKs1φKs1φKsKb(26)
L2=φKs1φKf1+Kb1φKs(27)

where Kb denotes the bulk modulus of dry skeleton.

The equivalent bulk modulus of the porous unit with Nc fractures, volume V, and Poisson’s ratio ν can be determined (Mavko and Nur, 1978):

1K=1Ks1+2π31ν212νi=1NcR1i2diV(28)

where R1i refers to the radius of the elliptical principal axis of the ith fracture, and di stands for fracture extension length. In a cubic unit cell with side length l, let M×N×L represent the spatial dimensions of the fracture network, where M fractures are distributed in the X direction, N fractures in the Y direction and L fractures in the Z direction. The relationship among the main radius of fracture network model, porosity ϕ, aspect ratio of fractures a, network scale and sample side length shows as follows (Xiong et al., 2021):

R1=MN+NL+MLφl2πa(29)

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

σe=Gu+uT23uI(30)

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

In conclusion, the constitutive relationship between solid and fluid can be summarized as follows:

σ=a11u¯+a12U¯I+Gu¯+u¯T23u¯Is=a21u¯+a22U¯I(31)

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

ρsu¨=G2u+a11e+a12ξ+CU˙u˙ρfU¨=a21e+a22ξCU˙u˙(32)

where, C=ηϕ2/κ, a11=a11+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 u and U. These equations exhibit hyperbolic characteristics with dissipative terms due to momentum transfer.

By applying the divergence operator to both sides of the wave Eq. 32 and introducing uP=u and UP=U, the P-wave equations is obtained:

ρs2uPt2=a112uP+a122UP+CUPtuPtρf2UPt2=a212uP+a222UPCUPtuPt(33)

By applying the curl operator to both sides of the wave Eq. 32 and introducing uS=×u and US=×U, the s-wave equations is derived:

ρs2uSt2=G2uS+CUStuStρf2USt2=CUStuSt(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:

uP=uP0eiωtkxUP=UP0eiωtkx(35)

where uP0 and UP0 refer to P-wave amplitude, and k denotes wave number, which is typically a complex quantity in nature; ω signifies circular frequency, and i represents imaginary unit. After substituting the above equations into the P-wave equation, the subsequent expression can be derived:

ω2ρs00ρf+k2a11a12a21a22+iωCCCCuP0UP0=00(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:

Z1X2+Z2X+Z3=0(37)
Z1=ρs+iCωρf+iCω+C2ω2(38)
Z2=a11ρfa22ρsiCa11+a12+a21+a22ω(39)
Z3=a11a22a12a21(40)

where X=ω2/k2. The analytical expression for P-wave velocity can be obtained by solving the above equations:

Vp=VpR+iVpI(41)
VpR=E2+F21/4signcosα212+E2E2+F2(42)
VPI=E2+F21/4signsinα212E2E2+F2(43)
tanα=FE(44)

where the coefficients E and F are

E=2a11ρ¯f+a22ρ¯sρ¯sρ¯f+2a11+a12+a21+a22ρ¯s+ρ¯fC2ω2±2ρ¯sρ¯fA±2ρ¯s+ρ¯fCωB4ρ¯s2ρ¯f2+4ρ¯s+ρ¯f2C2ω2(45)
F=2a11+a12+a21+a22ρ¯sρ¯fCω2a11ρ¯f+a22ρ¯sρ¯s+ρ¯fCω±2ρ¯s+ρ¯fCωA±2ρ¯sρ¯fB4ρ¯s2ρ¯f2+4ρ¯s+ρ¯f2C2ω2(46)
A=x2+y214signcosθ212+x2x2+y2(47)
B=x2+y214signsinθ212x2x2+y2(48)
tanθ=yx(49)
x=a11ρ¯f+a22ρ¯s24a11a22a12a21ρ¯fρ¯sa11+a12+a21+a222C2ω2(50)
y=2a11ρ¯f+a22ρ¯sa11+a12+a21+a22Cω4a11a22a12a21ρ¯f+ρ¯sCω(51)

where signx refers to the application of the sign of x, and ρ¯s=<ρs>, ρ¯f=<ρf>.

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

According to the method of plane wave analysis, the plane harmonic S-wave propagating along the z-direction can be expressed as:

uS=uS0eiωtkxUS=US0eiωtkx(52)

where uS0 and US0 refer to S-wave amplitude. After substituting the above equations into the S-wave equation, the subsequent expression can be derived:

ω2ρs00ρf+k2G000+iωCCCCuS0US0=00(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:

Y1X+Y2=0(54)
Y1=ρsρfiCωρs+ρf(55)
Y2=iCωρfG(56)

where X=ω2/k2. Thus the analytical expression for S-wave velocity can be solved as:

X=E+Fi(57)
E=ρ¯f2ρ¯sG+C2ω2ρ¯s+ρ¯fGρ¯sρ¯f2+C2ω2ρ¯s+ρ¯f2(58)
F=ρ¯fρ¯sGρ¯fρ¯s+ρ¯fGρ¯sρ¯f2+C2ω2ρ¯s+ρ¯f2Cω(59)

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

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 vf=U˙¯ and vs=u˙¯, the macroscopic equation for fluids and solid are rewritten as

ρfv˙f=sηϕ2κvfvs(60)
ρfv˙s=σ+ηϕ2κvfvs(61)

In the equation provided above, the brackets in denoting the volume average is omitted in the macro-scale equations for the sake of simplicity. Thus, the wave equation in velocity-stress format is:

v˙s=1ρsσ+Cρsvfvsv˙f=1ρfsCρfvfvs(62)

Here C=ηϕ2κ. By applying the time derivative to both sides of the constitutive equation, the governing equation for the fractured-porous model can be obtained:

σ˙ij=a11e˙+a12ξ˙δij+2Gε˙ij13e˙δijs˙ij=a21e˙+a22ξ˙δij(63)

It should be noted that the stress tensor exhibits symmetry, resulting in a total of 9 unknowns, namely, σxx, σxy, σxz, σyy, σyz, σzz, sxx, syy and szz, which are all associated with stress and therefore necessitate updates using constitutive relations. In the case of fluids, it satisfies the condition of sxx=syy=szz, 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 sxx, syy and szz instead.

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

FIGURE 1
www.frontiersin.org

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

The staggered grid differencing operators, denoted as Dh in space and Dt in time, are introduced as follows: (Sun, 2009):

Dhfi+1/2,j,kn=θ1fi+1,j,knθ2fi,j,kn+θ3fi+2,j,knθ4fi1,j,knΔh(64)
Dtfi,j,kn+1/2=fi,j,kn+1fi,j,knΔt(65)

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

Utilizing the staggered-grid finite difference, the discretization of unknown variables in the first-order velocity-stress equation occurs at distinct grid points, leading to the subsequent finite difference equations for velocity update:

Dtvxsi,j+1/2,kn1/2=1ρsDxσxxi+1/2,j+1/2,kn+Dyσxyi,j+1,kn+Dzσxzi,j+1/2,k+1/2n+Cρsvxfi,j+1/2,kn1/2vxsi,j+1/2,kn1/2(66)
Dtvysi+1,j,kn1/2=1ρsDxσxyi+1,j,kn+Dyσyyi+1/2,j+1/2,kn+Dzσyzi+1/2,j,k+1/2n+Cρsvyfi+1/2,j,kn1/2vysi+1/2,j,kn1/2(67)
Dtvzsi+1/2,j+1/2,k+1/2n1/2=1ρsDxσxzi+1,j+1/2,k+1/2n+Dyσyzi+1/2,j+1,kn+Dzσxzi+1/2,j+1/2,k+1n+Cρsvzfi+1/2,j+1/2,k+1/2n1/2vzsi+1/2,j+1/2,k+1/2n1/2(68)
Dtvxfi,j+1/2,kn1/2=1ρfDxsi+1/2,j+1/2,knCρfvxfi,j+1/2,kn1/2vxsi,j+1/2,kn1/2(69)
Dtvyfi+1/2,j,kn1/2=1ρfβ2Dysi+1/2,j+1/2,knCρfvyfi+1/2,j,kn1/2vysi+1/2,j,kn1/2(70)
Dtvzfi+1/2,j+1/2,k+1/2n1/2=1ρfDzsi+1/2,j+1/2,k+1nCρfvzfi+1/2,j+1/2,k+1/2n1/2vzsi+1/2,j+1/2,k+1/2n1/2(71)

The finite difference equations for stress update are as follows:

Dtσxxi+1/2,j+1/2,kn=a11Dxvxsi+1,j+1/2,kn+1/2+a11*Dyvysi+1/2,j+1,kn+1/2+a11*Dzvzsi+1/2,j+1/2,k+1/2n+1/2+a12Dxvxfi+1,j+1/2,kn+1/2+Dyvyfi+1/2,j+1,kn+1/2+Dzvzfi+1/2,j+1/2,k+1/2n+1/2(72)
Dtσyyi+1/2,j+1/2,kn=a11*Dxvxsi+1,j+1/2,kn+1/2+a11Dyvysi+1/2,j+1,kn+1/2+a11*Dzvzsi+1/2,j+1/2,k+1/2n+1/2+a12Dxvxfi+1,j+1/2,kn+1/2+Dyvyfi+1/2,j+1,kn+1/2+Dzvzfi+1/2,j+1/2,k+1/2n+1/2(73)
Dtσzzi+1/2,j+1/2,kn=a11*Dxvxsi+1,j+1/2,kn+1/2+a11*Dyvysi+1/2,j+1,kn+1/2+a11Dzvzsi+1/2,j+1/2,k+1/2n+1/2+a12Dxvxfi+1,j+1/2,kn+1/2+Dyvyfi+1/2,j+1,kn+1/2+Dzvzfi+1/2,j+1/2,k+1/2n+1/2(74)
Dtσxyi,j,kn=GDyvxsi,j+1/2,kn+1/2+Dxvysi+1/2,j+1,kn+1/2(75)
Dtσxzi,j+1/2,k+1/2n=GDzvxsi,j+1/2,k+1n+1/2+Dxvzsi+1/2,j+1/2,k+1/2n+1/2(76)
Dtσyzi+1/2,j,k+1/2n=GDzvysi+1/2,j,k+1n+1/2+Dyvzsi+1/2,j+1/2,k+1/2n+1/2(77)
Dtsi+1/2,j+1/2,kn=a21Dxvxsi+1,j+1/2,k+1/2n+1/2+a21Dyvysi+1/2,j+1,kn+1/2+a21Dzvzsi+1/2,j+1/2,k+1/2n+1/2+a22Dxvxfi+1,j+1/2,kn+1/2+Dyvyfi+1/2,j+1,kn+1/2+Dzvzfi+1/2,j+1/2,k+1/2n+1/2(78)

where a11*=a112G/3, a11=a11+4G/3.The notations σxxi,j,kn, σyyi,j,kn, σzzi,j,kn, σxyi,j,kn, σyzi,j,kn, σxzi,j,kn represent the solid stress at spatial grid points i,j,k and time grid n. si,j,kn is the discretized fluid stress. vxsi,j,kn, vysi,j,kn, vzsi,j,kn are discretized solid velocity. vxfi,j,kn, vyfi,j,kn, vzfi,j,kn 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 σxx, σyy and σzz in three directions of the solid phase, as well as the fluid stress s. The source function is based on Gaussian curve:

ft=2ξ12ξT2eξT2(79)

where ξ=F02/0.1512 represents pulse width, T=tts serves as time-shift parameter and ts=1.5/F0. The Higdon absorbing boundary conditions (Higdon, 1986; Higdon, 1987) are employed as the boundary conditions. Taking the second-order Higdon absorbing boundary condition as an illustration, the wave field at the boundary satisfies the following equation:

j=12cosαjtcxu=0(80)

where αj denotes the incident angle of absorption, and c refers to the propagation velocity of incident wave. The above formula is capable of effectively attenuating incident waves with angle of ±αj. Following discretization on the difference grid, the differential operator form of the absorbing boundary condition is presented as follows:

BEx,Et1=j=12βjIEt1Δt1aI+aExcExIΔx1bI+bEt1(81)

where βj=cosαj represents a positive number, and Δx signifies the spatial step size in the x direction on the boundary grid, and I denotes the identity operator, and ExI=Ex, 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, fi,j,kn represents the physical variable of the wave field at the time t=nΔt and at the point x=iΔx and i=0 (i.e., the left boundary in the x direction). Upon implementation of the boundary processing, it must satisfy the condition of zero reflection for the reflected wave, thereby enabling calculation of fi,j,kn through utilization of the absorption boundary difference formula:

fi,j,kn=1B9B1fi+2,j,kn2+B2fi+1,j,kn2+B3fi,j,kn2+B4fi+2,j,kn1+B5fi+1,j,kn1+B6fi,j,kn1+B7fi+2,j,kn+B8fi+1,j,kn(82)

Given the initial conditions, the wavefield’s physical variables at each node are known at time 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=Δt. Following simplification, the solution can be obtained

fi,j,kn=1A14A11fi+1,j,kn1+A12fi,j,kn1+A13fi+1,j,kn(83)

where each coefficient refers to as follows:

pj=cΔt/βjΔx,Aj1=pjb+a,Aj2=pjb+a1,Aj3=apj1b,Aj4=1a+pj1b,B1=A11A21,B2=A11A22+A12A21,B3=A12A22,B4=A11A23+A13A21,B5=A11A24+A12A23+A13A22+A14A21,B6=A12A24+A14A22,B7=A13A23 and B8=A13A24+A14A23.

Therefore, the calculation of the wavefield value at the node x=iΔx and i=0 (i.e., the left boundary in the x direction) with first-order absorbing boundary requires three layers of data: the wavefield value F1=fi=0,kn=0 of the previous time step’s boundary layer (i.e., t=0 and x=0); the wavefield value F2=fi=1,kn=0 of the first layer inside the previous time step’s boundary (i.e., t=0 and x=Δx); and the wavefield value F3=fi=1,kn=1 of the first layer inside the current time step’s boundary (i.e., t=Δt and x=Δ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
www.frontiersin.org

TABLE 1. Parameter table of porous media.

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

TABLE 2
www.frontiersin.org

TABLE 2. Parameter table of finite difference simulation.

The calculation space range is 250m×250m, 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
www.frontiersin.org

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
www.frontiersin.org

FIGURE 3. Dry skeleton modulus (Kb) 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 Ax1 denote the recorded amplitude value at position x1, while Ax2 represents the corresponding value at position x2. The separation between these receivers is denoted as Δs=x2x1. Secondly, obtain multiple sets of waveform graphs, correlate the waveforms at these two positions, determine the time difference Δt, and thereby derive the velocity value v=Δs/Δt. Additionally, the quality factor Q=πfx2x1vlnAx1Ax21 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
www.frontiersin.org

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

FIGURE 5
www.frontiersin.org

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 Kb, 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 103 Pa to 109 Pa. The seismic velocities can be obtained by substituting the corresponding Kb under different confining pressures, and the results show in Figures 6, 7.

FIGURE 6
www.frontiersin.org

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

FIGURE 7
www.frontiersin.org

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 Kb. Additionally, a similar relationship is observed between wave velocity and confining pressure, as well as between dry skeleton modulus Kb 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 pc 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 100–104 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 (<104 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
www.frontiersin.org

TABLE 3. Rock physical parameters.

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

FIGURE 8
www.frontiersin.org

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

FIGURE 9
www.frontiersin.org

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

In the example, various porous network structures of different scales (4×4×4, 5×5×5, 6×6×6, 7×7×7 and 8×8×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/m3, 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
www.frontiersin.org

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 ×105 D, 1.8 ×104 D, 0.036D, 0.2D, 0.47D) are adopted. Figure 11 illustrates the velocity ratio-impedance template.

FIGURE 11
www.frontiersin.org

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.m2). 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
www.frontiersin.org

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.m2 at a=1 to 11 kg/s.m2 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
www.frontiersin.org

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
www.frontiersin.org

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
www.frontiersin.org

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 (Vp/Vs) serves as an indicator of fluid-filled zone. P-wave impedance (Z) and wave velocities (Vp and Vs) 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 (Qp and Qs) 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 sx=xmaxxmin/x¯ , xmax and xmin represent the maximum and minimum calculated values of the seismic attributes respectively, while x¯ denotes their average value. Figures 1618 present comparative charts illustrating the sensitivity of each attribute for sandstone containing water, oil or gas.

FIGURE 16
www.frontiersin.org

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
www.frontiersin.org

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
www.frontiersin.org

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

Keywords: fractured-porous media, wave equation, permeability identification, seismic attributes, numerical simulation

Citation: Zhao W, Sun W, Gao J and He R (2023) An improved wave equation of fractured-porous media for predicting reservoir permeability. Front. Earth Sci. 11:1269616. doi: 10.3389/feart.2023.1269616

Received: 30 July 2023; Accepted: 13 November 2023;
Published: 29 December 2023.

Edited by:

Xintong Dong, Jilin University, China

Reviewed by:

Bing Bai, Beijing Jiaotong University, China
Danping Cao, China University of Petroleum (East China), China

Copyright © 2023 Zhao, Sun, Gao and He. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Wanjin Zhao, zhao_wj@petrochina.com.cn

Disclaimer: 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.