SUMMARY

Accurate modelling of the frequency-dependence of seismic wave velocity related to fracture system and fluid content is crucial to the quantitative interpretation of seismic data in fractured reservoirs. Both mesoscale fractures and patchy saturation effects can cause significant velocity dispersion and attenuation in the seismic frequency band due to wave-induced fluid flow (WIFF) mechanism. Considering the coupled impact of ‘mesoscale fractures’ and ‘patchy saturation’, we derive expressions for the frequency-dependent anisotropy in partially saturated porous rock containing two fracture sets with different orientations, sizes and connectivities. Especially, we simplify the rock-physics model as an orthorhombic (ORT) media by assuming the mesoscale fractures to be orthogonal and give the explicit expressions for frequency-dependent elastic constants. Finally, we give the expressions for the frequency-dependent phase velocity in patchy saturated and fractured ORT media and investigate the effect of patchy saturation on P-wave velocity at different polar and azimuth angles. In this paper, we investigate the effects of fluid saturation and fluid pressure on frequency-dependent velocities and Thomsen anisotropy parameters. Also, the effect of the relative permeability is very noticeable. The relaxation frequency can be lower in partially saturated fractured rocks compared with the fully saturated case, which makes the rock have a larger stiffness. The non-monotonic relationships between frequency-dependent anisotropy and fluid saturation add complexity to seismic forward modelling and inversion in reservoirs with complex fracture patterns.

1 INTRODUCTION

Accurate description of fracture system and fluid content in porous rock is a key problem to a variety of applications including oil/gas reservoir characterization and CO2 migration monitoring. Modelling and analysis of field observations suggest that frequency-dependent anisotropy can indicate fluid saturation and fracture size information in fractured reservoirs (Maultzsch et al. 2003; Maultzsch 2007). A proper theoretical rock-physics model is essential to extract fluid and fracture information from seismic data.

Wave-induced fluid flow (WIFF) which builds the relationship between fluid-saturated pore space and elastic properties, is widely known as a driving mechanism behind frequency-dependence of seismic velocity dispersion and attenuation in heterogeneous porous media (Müller et al. 2010). Therefore, the mesoscopic scale (larger than the pore size but smaller than the wavelength) fluid flow is increasingly believed to be the main cause of velocity dispersion and attenuation in the seismic exploration frequency band (Pride et al. 2004; Müller et al. 2010; Quintal 2012). Two most obvious scenarios of mesoscopic fluid flow that have gained substantial attention are patchy saturated rocks and fractured reservoirs.

For patchy saturation case, White (1975) and Dutta & Seriff (1979) first developed a discrete inclusion-based model in 3-D space to describe the influence of biphasic fluids on seismic waves. White et al. (1975) introduced an interlayer flow model to compute the frequency-dependent P-wave modulus for a partially saturated medium. Lots of theoretical refinements are proposed based on these two models (Norris 1993; Müller & Gurevich 2004; Müller et al. 2008). For anisotropic fracture case, Hudson et al. (1996) modelled penny-shaped cracks and quantified seismic dispersion and attenuation induced by a single crack neglecting the interactions. Gurevich (2003) gave the expressions for static elastic moduli of a fluid-saturated fractured rock with equant porosity from the properties of unfractured porous rock and fracture compliances of the dry fractured rock. Based on this model, Brajanovski et al. (2005) and Galvin & Gurevich (2006, 2007) introduced planar fractures and finite-sized penny-shaped cracks in a porous medium, respectively. Based on the unified analytical framework of Gurevich et al. (2009), Many researchers studied the effects of fracture connectivity on seismic dispersion in porous and fractured rocks (Rubino et al. 2014; Shao & Pyrak-Nolte 2016; Guo et al. 2017, 2018). Considering the coupled fluid motion on two scales: the grain scale (pores and microcracks) and fracture scale (mesoscale fractures), Chapman (2003) developed a model to describe the frequency-dependence of elastic stiffness tensor in porous media based on the approach of Chapman et al. (2002). The results of this model have reached a good agreement with measured data in laboratory (Maultzsch et al. 2003; Tillotson et al. 2011; Ding et al. 2017). Chapman (2009) extended this model to account for multiple sets of fractures with different sizes and orientations. Pang & Stovas (2020) compared the dispersion of anisotropy parameters due to layering and fracturing effect at low frequencies based on Chapman's (2003) model. According to Chapman (2009)’s model, Shuai et al. (2020) specified the expressions for frequency-dependent anisotropy due to two sets of orthogonal fractures and combined the fracture-induced with layer-induced anisotropy under the assumption of weak-contrast approximation.

To describe the coupled ‘squirt flow’ and ‘patchy saturation’ effects on seismic propagation in porous rock saturated by two immiscible fluids, Papageorgiou & Chapman (2017) combined the approaches of Papageorgiou & Chapman (2015) and Papageorgiou et al. (2016) into a consistent model valid over the full range of frequencies and capillary pressures. Jin et al. (2018) extended this isotropic model to account for fracture-induced anisotropy and pronounced the effect of relative permeability and used this model to explain laboratory measurements of seismic anisotropy in partially saturated and fractured rocks. The purpose of this paper is to develop a model to describe the frequency-dependent anisotropy in partially saturated rock with multiple sets of mesoscale fractures. We do this by extending the model of Jin et al. (2018) to contain multiple sets of fractures with different orientations, sizes and connectivities.

The paper is organized as follows. First, we give the expressions for partially saturated rock with two sets of mesoscale fractures. Then, we derive the equations for elastic constants of porous rock with two sets of orthogonal fractures and saturated by two immiscible fluids. Finally, we analyse the effects of uneven fluid pressure and fluid saturation on frequency-dependent anisotropy in our anelastic orthorhombic (ORT) model and discuss the coupled effects of ‘mesoscale fracture’ and ‘patchy saturation’ on seismic wave propagation.

2 FREQUENCY-DEPENDENT ANISOTROPY DUE TO MULTIPLE SETS OF MESOSCALE FRACTURES IN POROUS ROCK SATURATED BY TWO IMMISCIBLE FLUIDS

Following Papageorgiou & Chapman (2017), we assume that the wave-induced fluid pressure is different within each inclusion. The induced pressure of one fluid phase is proportionally to that of the other fluid phase and this relationship is described as
(1)

Here, |${P_g}$| and |${P_w}$| represent the pressure of the two fluids with the fluid modulus |${K_g} < {K_w}$|⁠. The parameter q quantifies variation of capillary pressure and lies within |$( {{q_0},1} ]$| where |${q_0} = \frac{{{K_g}}}{{{K_w}}}\ $| and the fluid modulus |${K_g}$| is smaller than |${K_w}$|⁠. |$q\ = \ 1$| corresponds to the iso-stress condition that the two kinds of fluid are mixed well. This condition is referred as ‘uniform saturation’ (Mavko & Mukerji 1998). Other values between |${q_0}$| and 1 correspond to unequilibrated pressures between the fluids which can help capture the effect of ‘patchy saturation’.

In this study, we consider the pores and two sets of mesoscale fractures with different sizes and orientations in porous rocks. The normal direction of the mth fracture set is defined by polar angle |${\theta _m}$| and azimuthal angle |${\varphi _m}$|⁠: |${n^m} = ( {\cos {\varphi _m}\sin {\theta _m},\sin {\varphi _m}\sin {\theta _m},\cos {\theta _m}} )\ ,\ m\ = \ 1,2$|⁠. With the assumption that the saturations |${S_w}\ $|are equal in each inclusion type, the fluid pressures in grain-scale pores and two sets of mesoscale fractures can be written as
(2)
where |${p^*}$|⁠, |${f_1}$| and |${f_2}$| are pore fluid pressures in pore, fracture set 1 and fracture set 2, respectively. The subscript g and w refer to the two immiscible fluids. The parameter |$\tilde{q} = {S_w}\ + q( {1 - {S_w}} )$|⁠.
For each fluid type, the fluid density |${\rho _f}$| can be estimated from the undisturbed fluid density |$\rho _f^0$|⁠, the fluid pressure |${P_f}$| and the fluid bulk modulus |${K_f}$|⁠:
(3)
The fluid mass in pore and fracture volume can be expressed by
(4)

Here, the superscript p, 1, 2 represent pore, fracture set 1 and fracture set 2, respectively. |${\phi _p}$|⁠, |${\phi _1}$| and |${\phi _2}$| corresponds to the porosity of pore, fracture set 1 and fracture set 2.

For porous rock saturated by two immiscible fluids, the fluid flow obeys the multifluid Darcy's law. The volume flux of each fluid |${Q_w}$|⁠, |${Q_g}$| depends on the fluid mobility and the pressure gradient along the x direction |${\partial _x}{P_w}$|⁠, |${\partial _x}{P_g}$|⁠:
(5)
where A is the cross-sectional area of the fluid flow, |${\eta _w}$| and |${\eta _g}$| represent the dynamic viscosity of the two fluids, |${\kappa _w}$| and |${\kappa _g}$| are the relative permeability of the two fluids and |$\kappa $| is the absolute permeability of the rock. For each fluid type, the relative permeability is usually assumed as an increasing function of its saturation. According to Jin et al. (2018), the model of relative permeability can be approximated by |${\kappa _w} = S_w^3\ $| and |${\kappa _g} = {( {1 - {S_w}} )^2}\ $|⁠, which largely fits the experimental curves for multiphase fluid flow in CO2 storage reservoirs in Benson et al. (2013).
According to the approximation in Chapman et al. (2002), the multifluid flow in pores and fractures in terms of applied stress and fluid pressure can be described by
(6)

Where |${\sigma _{ll}}$|⁠, |${\sigma _i}$| and |${\sigma _{ii}}$| are the normal stress acting on pore, fracture set 1 and fracture set 2, respectively. v is the Poisson's ratio of rock matrix and |${K_w}$|⁠, |${K_g}$| are bulk modulus of two fluid types. |$\mu $| is the shear modulus of the isotropic background medium. Also, we define the notations |$\ \ \sigma _c^1 = \frac{{\pi \mu {r_1}}}{{2( {1 - v} )}}\ $| and |$\sigma _c^2 = \frac{{\pi \mu {r_2}}}{{2( {1 - v} )}}\ $|⁠, where |${r_m}$| is the aspect ratios of the mth fracture set. In this work, we assume that the fractures are uniformly sized and neglect the dependence of velocity dispersion on fracture aspect ratio.

By extending the single-fluid model of Chapman (2009) to two-fluid case, we have the system of equations:
(7)
where |${A_1}$| and |${A_2}$| corresponds to the cross-sectional area of fracture set 1 and fracture set 2. The Darcy constants |${g_w}$| and |${g_g}$| are defined as
(8)
Here, l is a characteristic length scale and indicates the distance between the adjacent intergranular voids. |${M_w}$| and |${M_g}$| are fluid mobilities of water and gas, respectively, where
(9)
By substituting the eq. (6) into eq. (7), we can solve for the fluid pressure in each inclusion and derive the frequency-dependent elastic stiffness elements |${C_{ij}}( \omega )$| as shown in Appendix  A. The effective fluid bulk modulus |${K_f}$| and the effective fluid mobility |${M_f}$| are as follow
(10)

The effective fluid bulk modulus |${K_f}$| is a Ruess average of the two fluid moduli weighted by the patchy parameter q. The effective mobility of the fluid mixture |${M_f}$| is influenced by both patchy parameter q and the relative permeabilities |${\kappa _w}$|⁠, |${\kappa _g}$|⁠. Therefore, the characteristic frequencies of the velocity dispersion will be not only controlled by the fracture size, but also related to the patchy saturation effects and the relative permeability.

In partially saturated porous rock with two sets of mesoscale fractures, the velocity dispersion and attenuation are caused by the WIFF mechanism, which take place between fractures and their embedding background. Here, the characteristic frequencies related to fracture set 1 |${\omega _{f1}}$| and fracture set 2 |${\omega _{f2}}$| are defined by
(11)
where |$\omega _0^{\prime}$| and |$\omega _0^{^{\prime\prime}}$| are the values of |${\omega _{f1}}$| and |${\omega _{f2}}$| for full water saturation case, respectively. From eq. (11), we can find that the main difference between partially saturation and full water saturation cases comes from the effective mobility of the fluid mixture. A lower fluid mobility results in a lower characteristic frequency. This can be explained by that the fluid with lower mobility needs more time for the pressure gradient to relax and reach the pressure equilibrium, which in turn leads to a lower characteristic frequency. Therefore, the dispersion and attenuation in partially saturated porous rock with mesoscale fractures are more complex than single fluid case.

Considering the orientations of the two fractures as |${\theta _1} = \pi /3\ $|⁠, |${\varphi _1} = \pi /3\ $|⁠, |${\theta _2} = \pi /6\ $| and |${\varphi _2} = \pi /6\ $|⁠, we study the effects of patchy saturation and dispersion on elastic stiffness tensor. Table 1 shows the parameter settings in the theoretical model describing porous rock saturated by water and supercritical CO2. Fig. 1 shows the variations of |${C_{33}}$| with water saturation for a range of the patchy parameter q values. In the case of |$\omega _0^{\prime}$| frequency, the elastic stiffness coefficient |${C_{33}}$| depends on water saturation and the patchy parameter q. Both the bulk modulus and the viscosity of the water are higher than those of gas. As a result, the |${C_{33}}$| increases as the water saturation increases. What's more, as patchy parameter q increases, the |${C_{33}}$| decreases. In the case of zero-frequency limit, the |${C_{33}}$| decreases with increase in patchy parameter q. What's more, the result when |$q\ = \ 1$| is consistent with anisotropic Gassmann-wood theory (Gassmann 1951) in which the bulk and shear moduli are irrelevant to frequency.

The variations of ${C_{33}}$with water saturation for different q values at zero frequency (solid lines) and $\omega _0^{\prime}$ frequency (dashed lines).
Figure 1.

The variations of |${C_{33}}$|with water saturation for different q values at zero frequency (solid lines) and |$\omega _0^{\prime}$| frequency (dashed lines).

Table 1.

Parameters of porous rock saturated by water and gas.

Lamé parameter of rock matrix (Pa)|$\lambda $||$1.58\ \times {10^{10}}$|Aspect ratio of fracture set 2|${r_2}$||$1\ \times {10^{ - 4}}$|
Lamé parameter of rock matrix (Pa)|$\mu $||$9.31\ \times {10^9}$|Fracture density 1|${\varepsilon _{f1}}$|0.03
Solid density (kg m–3)|$\rho $|2650Fracture density 2|${\varepsilon _{f2}}$|0.02
Porosity|${\phi _p}$|10 per centCharacteristic frequency of fracture set 1 at full water-saturation (Hz)|$\omega _0^{\prime}$||$50$|
Aspect ratio of fracture set 1|${r_1}$||$1\ \times {10^{ - 4}}$|Characteristic frequency of fracture set 2 at full water-saturation (Hz)|$\omega _0^{^{\prime\prime}}$|10
Water viscosity (Pa·s)|${\eta _w}$||$6\ \times {10^{ - 4}}$|CO2 viscosity (Pa·s)|${\eta _g}$||$2.1\ \times {10^{ - 5}}$|
Water density (kg m–3)|${\rho _w}$|1000CO2 density (kg/m3)|${\rho _g}$|240
Water modulus (Pa)|${K_w}$||$2.4\ \times {10^9}$|CO2 modulus (Pa)|${K_g}$||$1.1\ \times {10^7}$|
Lamé parameter of rock matrix (Pa)|$\lambda $||$1.58\ \times {10^{10}}$|Aspect ratio of fracture set 2|${r_2}$||$1\ \times {10^{ - 4}}$|
Lamé parameter of rock matrix (Pa)|$\mu $||$9.31\ \times {10^9}$|Fracture density 1|${\varepsilon _{f1}}$|0.03
Solid density (kg m–3)|$\rho $|2650Fracture density 2|${\varepsilon _{f2}}$|0.02
Porosity|${\phi _p}$|10 per centCharacteristic frequency of fracture set 1 at full water-saturation (Hz)|$\omega _0^{\prime}$||$50$|
Aspect ratio of fracture set 1|${r_1}$||$1\ \times {10^{ - 4}}$|Characteristic frequency of fracture set 2 at full water-saturation (Hz)|$\omega _0^{^{\prime\prime}}$|10
Water viscosity (Pa·s)|${\eta _w}$||$6\ \times {10^{ - 4}}$|CO2 viscosity (Pa·s)|${\eta _g}$||$2.1\ \times {10^{ - 5}}$|
Water density (kg m–3)|${\rho _w}$|1000CO2 density (kg/m3)|${\rho _g}$|240
Water modulus (Pa)|${K_w}$||$2.4\ \times {10^9}$|CO2 modulus (Pa)|${K_g}$||$1.1\ \times {10^7}$|
Table 1.

Parameters of porous rock saturated by water and gas.

Lamé parameter of rock matrix (Pa)|$\lambda $||$1.58\ \times {10^{10}}$|Aspect ratio of fracture set 2|${r_2}$||$1\ \times {10^{ - 4}}$|
Lamé parameter of rock matrix (Pa)|$\mu $||$9.31\ \times {10^9}$|Fracture density 1|${\varepsilon _{f1}}$|0.03
Solid density (kg m–3)|$\rho $|2650Fracture density 2|${\varepsilon _{f2}}$|0.02
Porosity|${\phi _p}$|10 per centCharacteristic frequency of fracture set 1 at full water-saturation (Hz)|$\omega _0^{\prime}$||$50$|
Aspect ratio of fracture set 1|${r_1}$||$1\ \times {10^{ - 4}}$|Characteristic frequency of fracture set 2 at full water-saturation (Hz)|$\omega _0^{^{\prime\prime}}$|10
Water viscosity (Pa·s)|${\eta _w}$||$6\ \times {10^{ - 4}}$|CO2 viscosity (Pa·s)|${\eta _g}$||$2.1\ \times {10^{ - 5}}$|
Water density (kg m–3)|${\rho _w}$|1000CO2 density (kg/m3)|${\rho _g}$|240
Water modulus (Pa)|${K_w}$||$2.4\ \times {10^9}$|CO2 modulus (Pa)|${K_g}$||$1.1\ \times {10^7}$|
Lamé parameter of rock matrix (Pa)|$\lambda $||$1.58\ \times {10^{10}}$|Aspect ratio of fracture set 2|${r_2}$||$1\ \times {10^{ - 4}}$|
Lamé parameter of rock matrix (Pa)|$\mu $||$9.31\ \times {10^9}$|Fracture density 1|${\varepsilon _{f1}}$|0.03
Solid density (kg m–3)|$\rho $|2650Fracture density 2|${\varepsilon _{f2}}$|0.02
Porosity|${\phi _p}$|10 per centCharacteristic frequency of fracture set 1 at full water-saturation (Hz)|$\omega _0^{\prime}$||$50$|
Aspect ratio of fracture set 1|${r_1}$||$1\ \times {10^{ - 4}}$|Characteristic frequency of fracture set 2 at full water-saturation (Hz)|$\omega _0^{^{\prime\prime}}$|10
Water viscosity (Pa·s)|${\eta _w}$||$6\ \times {10^{ - 4}}$|CO2 viscosity (Pa·s)|${\eta _g}$||$2.1\ \times {10^{ - 5}}$|
Water density (kg m–3)|${\rho _w}$|1000CO2 density (kg/m3)|${\rho _g}$|240
Water modulus (Pa)|${K_w}$||$2.4\ \times {10^9}$|CO2 modulus (Pa)|${K_g}$||$1.1\ \times {10^7}$|

3 FREQUENCY-DEPENDENT ANISOTROPY DUE TO TWO ORTHOGONAL MESOSCALE FRACTURES IN POROUS ROCK SATURATED BY TWO IMMISCIBLE FLUIDS

Next, we consider a specific case of the theory in chapter 2 and set the two sets of fractures are orthogonal. We note the indices 1 and 2 correspond to horizontal and vertical aligned fractures and the orientation angles are set as |${\theta _1} = \ 0$|⁠, |${\varphi _1} = \ 0$|⁠, |${\theta _2} = \pi /2\ $| and |${\varphi _2} = \ 0$|⁠. The partially saturated porous rock with two sets of orthogonal fractures can be described as anelastic ORT media, as shown in Fig. 2. The effective frequency-dependent elastic stiffness tensor of a partially saturated rock with two sets of orthogonal mesoscale fractures |${C_{ij}}( \omega )$| can be calculated by following Appendix  A. What's more, we can give the specific expression for the fluid-pressure in pores in this anelastic ORT model
(12)
Schematic of porous rock with orthogonal fractures.
Figure 2.

Schematic of porous rock with orthogonal fractures.

The elastic correction of the fracture set 1 in standard form is
(13)
and the elastic correction of the fracture set 2 is
(14)
The elements of |$C_{ij}^1( \omega )$| and |$C_{ij}^2( \omega )$| are define in eq. (A17). The fluid pressure in fracture set 1 can be expressed by
(15)
(16)
(17)
(18)
Also, the fluid pressures in fracture set 2 are
(19)
(20)
(21)
(22)

4 PHASE VELOCITY IN PARTIALLY SATURATED ROCK CONTAINING TWO SETS OF ORTHOGONAL FRACTURES

The phase velocity V and the displacement vector U of plane waves in anelastic ORT media satisfy the Christoffel equation as follows
(23)
where |$\rho $| is the density of the porous rock, |${\delta _{ik}}$| is the Kroneker's symbolic unit matrix and |${G_{ik}}$| is the symmetric Christoffel matrix. For partially saturated rock with two sets of orthogonal fractures, the elements of |${G_{ik}}$| are frequency-dependent:
(24)
(25)
(26)
(27)
(28)
(29)
Here, the |${\rm{\Theta }}$| and |${\rm{\Phi }}$| are polar and azimuthal phase angles of the phase velocity vector, respectively. The notations |$\lambda $|⁠, |$\mu $|⁠, d, e, |$a_{ij}^1$| and |$a_{ij}^2$| are defined in Appendix  A. Then, by setting the determinant of |$[ {{G_{ik}} - \rho {V^2}{\delta _{ik}}} ]$| to be zero, we can obtain a cubic equation for frequency-dependent phase velocities in anelastic ORT media:
(30)
where the frequency-dependent coefficients |$a( \omega )$|⁠, |$b( \omega )$| and |$c( \omega )$| are expressed by
(31)
(32)
(33)
According to the solution of the cubic equation given by Schoenberg & Helbig (1997) and Tsvankin (1997), we only take the real part of the |${G_{ik}}$| and the frequency-dependent phase velocity of anelastic ORT media can be expressed as
(34)
Where the largest root |$( {k\ = \ 0} )\ $|corresponds to the phase velocity of the quasi P-wave, while the other two roots |$( {k\ = \ 1,2} )$| correspond to the quasi split shear waves. The directional velocities for the quasi waves will change according to the symmetry of the media. The equations for |$\upsilon ( \omega )$|⁠, |$d( \omega )$| and |$q( \omega )$| are defined by
(35)
(36)
(37)
The condition for the roots of eq. (30) to be real, the frequency-dependent coefficient |$d( \omega )$| should be negative (Schoenberg & Helbig 1997) and the following combination of |$d( \omega )$| and |$q( \omega )$| should be non-positive
(38)
If |$Q\ ( \omega ) = \ 0$| two of the roots are identical, which leads to shear wave singularities in certain directions. By introducing the weak anisotropy approximation for ORT medium (Tsvankin 1997), the P-wave velocity versus frequency, polar and azimuthal angle is given by
(39)
where |${V_{P0}}$| is the vertical P-wave velocity and coefficients |$\varepsilon $| and |$\delta $| are given by
(40)
(41)

The velocities and anisotropy parameters in vertical direction are described by equations in Appendix  B.

5 NUMERICAL EXAMPLE

To illustrate the velocity dispersion and attenuation due to WIFF which takes place between fractures and their embedding background in partially saturated rock with two orthogonal mesoscale fractures, we use parameter settings in Table 1 and investigate the vertical velocities and Thomsen-style anisotropy parameters vary with water saturation and fluid pressure. Fig. 3 shows the variations of vertical P-wave velocity and attenuation factor with water saturations and patchy parameter q at zero frequency and |$\omega _0^{\prime}$| frequency. At zero-frequency limit, the vertical P-wave velocity increases with increase in water saturation wherever decreases with increase q value. It is obvious that the relationship between P-wave velocity and water saturation is monotonic. The attenuation factor is zero which means there is no energy attenuation at zero frequency. At |$\omega _0^{\prime}$| frequency, there is a non-monotonic variation of P-wave velocity with water saturation in the case of |$q\ = {q^0}\ $|⁠. This is related to the relative permeability effect that leads to a lower fluid mobility and a higher bulk modulus of Pwave for intermediate water saturations. As water saturation increases, the attenuation factor first increases then decreases and reaches the maximum value. The characteristic water saturation corresponding to the maximum attenuation factor increases with increase in patchy parameter q.

The variations of vertical P-wave velocity: (a) and attenuation factor and (b) with water saturation for different q values at zero frequency (solid lines) and $\omega _0^{\prime}$ frequency (dashed lines).
Figure 3.

The variations of vertical P-wave velocity: (a) and attenuation factor and (b) with water saturation for different q values at zero frequency (solid lines) and |$\omega _0^{\prime}$| frequency (dashed lines).

Fig. 4 shows the variation of vertical P-wave velocity with frequencies and water saturations at different q values. The P-wave velocity increases with increase in frequency. The relationship between P-wave velocity and water is very complex and is non-monotonic when |$q\ = {q_0}\ $|⁠. As shown in Fig. 4(a), the P-wave velocity at intermediate water saturation is higher than that at fully water saturation ( |${S_w} = \ 1$|⁠) at 0–10 Hz frequency band. This may be associated with the relative permeability effect which leads to a lower the fluid mobility and a higher P-wave velocity for intermediate water saturation at some frequencies. As shown in Figs 4(a) and (b), the characteristic frequency of the P-wave velocity first decreases the increases as the water saturation. In the cases of |$q\ = \ 50{q_0}$| and |$q\ = \ 1$|⁠, the P-wave velocity increases with increase in frequency but decreases with increase water saturation. Also, the P-wave velocities for partially saturation cases lie between velocities for |${S_w} = \ 0$| and |${S_w} = \ 1$|⁠. The characteristic frequency of the P-wave velocity decreases as water saturation increases. The variations of attenuation factor with frequencies and water saturations at different q values are show in Fig. 5. The characteristic frequency depends both on water saturation and patchy parameter q. the maximum of the attenuation factor first increases then decreases as the water saturation increases and show non-monotonic trend.

The variations of P-wave velocity with frequency for different water saturation. (a) $q\ = {q_0}\ $, (b) $q\ = \ 10{q_0}$, (c) $q\ = 50{q_0}\ $ and (d) $q\ = \ 1$.
Figure 4.

The variations of P-wave velocity with frequency for different water saturation. (a) |$q\ = {q_0}\ $|⁠, (b) |$q\ = \ 10{q_0}$|⁠, (c) |$q\ = 50{q_0}\ $| and (d) |$q\ = \ 1$|⁠.

The variations of attenuation factor with frequency for different water saturation. (a) $q\ = {q_0}\ $, (b) $q\ = \ 10{q_0}$, (c) $q\ = 50{q_0}\ $ and (d) $q\ = \ 1$.
Figure 5.

The variations of attenuation factor with frequency for different water saturation. (a) |$q\ = {q_0}\ $|⁠, (b) |$q\ = \ 10{q_0}$|⁠, (c) |$q\ = 50{q_0}\ $| and (d) |$q\ = \ 1$|⁠.

Fig. 6 shows the variations of Thomsen-style anisotropy parameters with water saturations and q values at zero frequency and |$\ \omega _0^{\prime}$| frequency. As q increases, Thomsen-style anisotropy parameters increase. In |$[ {{x_1},{x_3}\ } ]$| plane (Figs 6a and b), the anisotropy parameter difference between |${S_w} = \ 0$| and |${S_w} = \ 1$| at |$\omega _0^{\prime}$| frequency is obviously larger than that at zero frequency. Both the anisotropy parameter |${\varepsilon _1}$| and |${\delta _1}$| first increases then decreases and finally increases as water saturation increases and show a non-monotonic relationship at |$\omega _0^{\prime}$| frequency. The anisotropy parameter |${\delta _1}$| can be negative at |$\omega _0^{\prime}$| frequency. For anisotropy parameters |${\varepsilon _2}$|⁠, |${\delta _2}$| and |${\delta _3}$|⁠, as shown in Figs 6(c), (d) and (e) the anisotropy parameter |${\varepsilon _2}$|⁠, |${\delta _2}$| and |${\delta _3}$| decrease with increase in frequency and water saturation. Also, the increase in q value can enlarge the anisotropy parameters.

The variations of Thomsen-style anisotropy parameters with water saturation for different q values at zero frequency (solid lines) and $\omega _0^{\prime}$ frequency (dashed lines). (a) ${\varepsilon _1}$, (b) ${\delta _1}\ $, (c) ${\varepsilon _2}\ $, (d) ${\delta _2}$ and (e) ${\delta _3}$.
Figure 6.

The variations of Thomsen-style anisotropy parameters with water saturation for different q values at zero frequency (solid lines) and |$\omega _0^{\prime}$| frequency (dashed lines). (a) |${\varepsilon _1}$|⁠, (b) |${\delta _1}\ $|⁠, (c) |${\varepsilon _2}\ $|⁠, (d) |${\delta _2}$| and (e) |${\delta _3}$|⁠.

Based on the anelastic ORT model we proposed, we investigate the influences of uneven fluid pressure and water saturation on angel-and azimuth dependent phase velocities. Fig. 7 shows vertical P-wave velocity varies with polar and azimuthal angles for different q values with water saturation is 0.5. From Fig. 7 (a), we can observe that in the case of |$q\ = {q_0}\ $|⁠, the P-wave velocity increases with increase polar angle |${\rm{\Theta }}$|⁠, and theP-wave velocity first increase then decrease as azimuth angle |${\rm{\Phi }}$| increases and the minimum P-wave velocity at |${\rm{\Phi \ }} = {\rm{\ }}45^\circ $| for each polar angle profile. The P-wave velocity firstly decreases then increases with increase in polar angle |${\rm{\Theta }}$| at |${\rm{\Phi \ }} = {\rm{\ }}90^\circ $| profile and reaches the minimum velocity at |${\rm{\Theta \ }} = {\rm{\ }}45^\circ $|⁠. For|$\ q\ = \ 10{q_0}$| and |$q\ = 50{q_0}\ $|⁠, the P-wave velocity firstly increases then decreases with increase in polar angle |${\rm{\Theta }}$| at |${\rm{\Phi \ }} = {\rm{\ }}90^\circ $| profile and reaches the trough value at |${\rm{\Theta \ }} = {\rm{\ }}45^\circ $|⁠, which corresponds to the averaged polar angle of the two fractures |${{\rm{\Theta }}_1}$| and |${{\rm{\Theta }}_2}$|⁠. For |$\ q\ = \ 1$|⁠, the P-wave velocity increases as polar angle |${\rm{\Theta }}$| increases at |${\rm{\Phi \ }} = {\rm{\ }}90^\circ $| profile.

Vertical P-wave velocity versus polar and azimuthal angles for different q values. (a) $q\ = {q_0}\ $, (b) $q\ = 10{q_0}\ $, (c) $q\ = 50{q_0}\ $ and (d) $q\ = \ 1$.
Figure 7.

Vertical P-wave velocity versus polar and azimuthal angles for different q values. (a) |$q\ = {q_0}\ $|⁠, (b) |$q\ = 10{q_0}\ $|⁠, (c) |$q\ = 50{q_0}\ $| and (d) |$q\ = \ 1$|⁠.

Fig. 8 shows vertical P-wave velocity varies with polar and azimuthal angles for different water saturations with |$\ q\ = {q_0}\ $|⁠. As water saturation increase, the P-wave velocity increases obviously at different polar and azimuth angles. For |${S_w} = \ 0$| (Fig. 8a), the P-wave velocity increase with |${\rm{\Theta \ and\ \Phi }}.$| For |${S_w} = \ 0.2$|⁠, |${S_w} = \ 0.4$| and |${S_w} = \ 0.6$|⁠, the P-wave velocity firstly decreases then increases and reaches the minimum velocity at |${\rm{\Phi \ }} = \ 45^\circ $| in |${\rm{\Theta \ }} = \ 90^\circ $| profile. For each polar angle |${\rm{\Theta }}$| profile, the P-wave velocity firstly decreases then increases and reaches the minimum velocity at |${\rm{\Theta \ }} = \ 45^\circ $|⁠. For |${S_w} = \ 0.8$|⁠, the variation of P-wave velocity with azimuth angle |${\rm{\Phi }}$| is not obvious in |${\rm{\Theta \ }} = \ 90^\circ $| profile. The P-wave velocity firstly decreases then increases but the minimum velocity is not at |${\rm{\Theta \ }} = \ 45^\circ $| in polar angle |${\rm{\Theta \ }} = {\rm{\ }}0^\circ $| profile. The P-wave velocity increases with increase in azimuth angle |${\rm{\Phi }}$| in |${\rm{\Theta \ }} = {\rm{\ }}90^\circ $|⁠. For |${S_w} = \ 1$|⁠, the P-wave velocity firstly increases then decreases and reaches the maximum velocity at |${\rm{\Phi \ }} = \ 45^\circ $|⁠. For each azimuth angle |${\rm{\Phi }}$| profile, the P-wave velocity increases with increase in polar angle |${\rm{\Theta }}$|⁠.

Vertical P-wave velocity versus polar and azimuthal angles for different water saturation. (a) ${S_w} = \ 0$, (b) ${S_w} = \ 0.2$, (c) ${S_w} = \ 0.4$, (d) ${S_w} = \ 0.6$, (e) ${S_w} = \ 0.8$ and (f) ${S_w} = \ 1$.
Figure 8.

Vertical P-wave velocity versus polar and azimuthal angles for different water saturation. (a) |${S_w} = \ 0$|⁠, (b) |${S_w} = \ 0.2$|⁠, (c) |${S_w} = \ 0.4$|⁠, (d) |${S_w} = \ 0.6$|⁠, (e) |${S_w} = \ 0.8$| and (f) |${S_w} = \ 1$|⁠.

6 DISCUSSION

Seismic monitoring of injected CO2 plumes and hydrocarbon production in fractured reservoirs depends on accurate rock physics model describing the fracture systems and fluid content. In this paper, we account for patchy saturation effects by introducing a non-dimensional parameter |$q\ $| which lies between the ratio of fluid bulk moduli and 1.

Papageorgiou & Chapman (2017) developed a model that describe the multifluid saturated rock through a unified model in which the squirt flow mechanism, the unbalanced fluid pressure and the relative permeability are considered. Jin et al. (2018) extended this theory to contain the mesoscale fractures and analyse the non-monotonic relationship between SWS with uneven fluid pressure and relative permeability. The Brine-CO2 partial saturation experiment in cracked sandstone with a well-defined fracture network indicate that P-wave velocity provides more reliable information about CO2 content than shear wave anisotropy (Falcon-Suarez et al. 2020). In this study, we give the theoretical models to describe the coupled effects of ‘mesoscale fractures’ and ‘patchy saturation’ by introducing the non-dimensional parameter q in Chapman's model (2009). This model emphasizes the importance of accurate modelling when we deal with partially saturated reservoir with complex fracture patterns.

To investigate the effects of relative permeability on vertical P-wave velocity in partially saturated rock with two orthogonal mesoscale fractures, we calculate the results from eq. (10) and the equation from Amalokwu et al. (2014):
(42)

Fig. 9 compares the P-wave velocities from the theoretical model with relative permeability effect (solid line) and without relative permeability effect (dashed line). For two frequency cases, the relationship between the P-wave velocity and water saturation changes from monotonic to non-monotonic due to relative permeability effect. However, the variation of P-wave velocity with water saturation with relative permeability effect for two frequencies are different.

The variations of P-wave velocity with water saturation for a range of q values at two frequencies. The solid lines represent results with relative permeability effects and the dashed lines represent results given by eq. (42) without relative permeability effects. (a) $\omega _0^{\prime}$ frequency case and (b) $\omega _0^{^{\prime\prime}}$ frequency case.
Figure 9.

The variations of P-wave velocity with water saturation for a range of q values at two frequencies. The solid lines represent results with relative permeability effects and the dashed lines represent results given by eq. (42) without relative permeability effects. (a) |$\omega _0^{\prime}$| frequency case and (b) |$\omega _0^{^{\prime\prime}}$| frequency case.

In the case of partially saturated rock with two sets of orthogonal fractures, setting orientations of fractures can be equivalent to changing in observe direction of seismic wave propagation. In this model, we assume that the relative permeability models are the same in all inclusions. Future work will address the geometry of each inclusions and combined multiple relative permeability model in this approach.

7 CONCLUSION

We have presented an extension on the model of Jin et al. (2018) and give the expressions for frequency-dependent anisotropy due to patchy saturation and multiple sets of mesoscale fractures. The elastic stiffness varies significantly with q and water saturation. Specially, we derive the equations in the case of the two sets of fractures are orthogonal and describe the seismic propagation from different directions. The characteristic frequency can be lowered due to the relative permeability, compared with fully water saturation case. The maximum attenuation factor depends on water saturation and patchy parameter |$\ q$|⁠. The relationship between the P-wave velocity and water saturation can be non-monotonic due to the relative permeability. The relationship of anisotropic parameters with water saturation changes with frequency. The patchy saturation effects can influence not only the magnitude of P-wave velocity, but also the variations of P-wave velocity with polar and azimuth angles. Therefore, accurate descriptions of both the fracture system and patchy effects are of great importance to seismic anisotropy in forward and inversion process in fractured reservoirs.

ACKNOWLEDGEMENTS

This research work is funded by the National Natural Science Foundation of China (No. 52074251 and No. 92058211), the Fundamental Research Funds for Central Universities of China (No. 202012003) and 111 project (No. B20048). Alexey Stovas would like to acknowledge the Petromaks 2 project for financial support.

DATA AVAILABILITY

Data and code related this research are available from the first author upon reasonable request.

REFERENCES

Amalokwu
K.
Chapman
M.
Best
A.I.
Sothcott
J.
Minshull
T.A.
Li
X.Y.
2014
Experimental observation of water saturation effects on S-wave splitting in synthetic rock with fractures aligned at oblique angles
Geophys. J. Int.
200
17
24
..

Benson
S.
Pini
R.
Reynolds
C.
Krevor
S.
2013
Relative permeability analyses to describe multi-phase flow in CO2 storge reservoirs
GCCSI Report No. 2 Relative Permeability, Global CCS Institute
.

Bond
W.L.
,
1943
.
The mathematics of the physical properties of crystals
,
Bell Syst. Tech. J.
,
22
,
1
72
.

Brajanvski
M.
Gurevich
B.
Schoenberg
M.
2005
A model for p-wave attenuation and dispersion in porous medium permeated by aligned fractures
Geophys. J. Int.
163
372
384
..

Chapman
M.
,
2003
.
Frequency-dependent anisotropy due to meso-scale fractures in the presence of equant porosity
,
Geophys. Prospect.
,
51
,
369
379
.

Chapman
M.
,
2009
.
Modeling the effect of multiple sets of mesoscale fractures in porous rock on frequency-dependent anisotropy
,
Geophysics
,
74
,
D97
D103
.

Chapman
M.
,
Zatsepin
S.V.
,
Crampin
S.
,
2002
.
Derivation of a microstructural poroelastic model
,
Geophys. J. Int.
,
151
,
427
451
.

Ding
P.
,
Di
B.
,
Wang
D.
,
Wei
J.
,
Li
X.Y.
,
2017
.
Measurements of seismic anisotropy in synthetic rocks with controlled crack geometry and different crack densities
.
Pure appl. Geophys.
,
174
.
1907
1922
.

Dutta
N.C.
,
Seriff
A.J.
,
1979
.
On White's model of attenuation in rocks with partial gas saturation
,
Geophysics
,
44
,
1806
1812
.

Falcon-Suarez
I.H.
Papageorgiou
G.
Jin
Z.
Muñoz-Ibáñez
A.
Chapman
M.
Angus
I.B.
2020
CO2-Brine substitution effects on ultrasonic wave propagation through sandstone with oblique fractures
Geophys. Res. Lett.
47
e2020GL088439
,doi:.

Galvin
R.J.
,
Gurevich
B.
,
2006
.
Interaction of an elastic wave with a circular crack in a fluid-saturated porous medium
,
Appl. Phys. Lett.
,
88
,
061918
.

Galvin
R.J.
,
Gurevich
B.
,
2007
.
Scattering of a longitudinal wave by a circular crack in a fluid-saturated porous medium
,
Int. J. Solids Struct.
,
44
,
7389
7398
.

Gassmann
F.
,
1951
.
Elastic waves through a packing of spheres
,
Geophysics
,
16
,
673
685
.

Guo
J.
Rubino
J.G.
Glubokovskikh
S.
Gurevich
B.
2017
Effects of fracture intersections on seismic dispersion: theoretical predictions versus numerical simulations
Geophys. Prospect.
65
1264
1276
..

Guo
J.
Rubino
J.G.
Glubokovskikh
S.
Gurevich
B.
2018
Dynamic seismic signatures of saturated porous rocks containing two orthogonal sets of fractures: theory versus numerical simulations
Geophys. J. Int.
213
1244
1262
.

Gurevich
B.
2003
Elastic properties of saturated porous rocks with aligned fractures
J. appl. Geophys.
54
203
218
.

Gurevich
B.
,
Brajanovski
M.
,
Galvin
R.J.
,
Müller
T.M.
,
Toms-Stewart
J.
,
2009
.
P-wave dispersion and attenuation in fractured and porous reservoirs-poroelasticity approach
,
Geophys. Prospect.
,
57
,
225
237
.

Hudson
I.A.
,
Liu
E.
,
Crampin
S.
,
1996
.
The mechanical properties of materials with interconnected cracks and pores
,
Geophys. J. Int.
,
124
,
105
112
.

Jin
Z.
,
Chapman
M.
,
Papageorgiou
G.
,
2018
.
Frequency-dependent anisotropy in a partially saturated fractured rock
.
Geophys. J. Int.
,
215
,
1985
1998
.

Maultzsch
S.
,
2007
.
Modeling and analysis of attenuation anisotropy in multiazimuth VSP data from the Clair field
,
Geophys. Prospect.
,
55
,
627
642
.

Maultzsch
S.
,
Chapman
M.
,
Liu
E.
,
Li
X.Y.
,
2003
.
Modelling frequency-dependent seismic anisotropy in fluid saturated rock with aligned fractures: implication of fracture size estimation from anisotropic measurements
,
Geophys. Prospect.
,
51
,
381
392
.

Mavko
G.
,
Mukerji
T.
,
1998
.
Bounds on low-frequency seismic velocities in partially saturated rocks
.
Geophysics
,
63
,
918
924
.

Müller
T.M.
Gurevich
B.
Lebedev
M.
2010
Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks—a review
Geophysics
75
A147
A164
.

Müller
T.M.
Toms-Stewart
J.
Wenzlau
F.
2008
Velocity-saturation relation for partially saturated rocks with fractal pore fluid distribution
Geophys. Res. Lett.
35
L09306
,doi:.

Müller
T.M.
,
Gurevich
B.
,
2004
.
One-dimensional random patchy saturation model for velocity and attenuation in porous rocks
,
Geophysics
,
69
,
1166
1172
.

Norris
A.N.
,
1993
.
Low-frequency dispersion and attenuation in partially saturated rocks
,
J. acoust. Soc. Am.
,
94
,
359
370
.

Pang
S.
,
Stovas
A.
,
2020
.
Low-frequency anisotropy in fractured and layered media
,
Geophys. Prospect.
,
68
,
353
370
.

Papageorgiou
G.
,
Amalokwu
K.
,
Chapman
M.
,
2016
.
Theoretical derivation of a Brie-like fluid mixing law
,
Geophys. Prospect.
,
64
,
1048
1053
.

Papageorgiou
G.
,
Chapman
M.
,
2015
.
Multifluid squirt flow and hysteresis effects on the bulk modulus-water saturation relationship
,
Geophys. J. Int.
,
203
,
814
817
.

Papageorgiou
G.
,
Chapman
M.
,
2017
.
Wave-propagation in rocks saturated by two immiscible fluids
,
Geophys. J. Int.
,
209
,
1761
1767
.

Pride
S.R.
,
Berryman
J.G.
,
Harris
J.M.
,
2004
.
Seismic attenuation due to wave induced flow
,
J. geophys. Res.
,
109
,
B01201
.

Quintal
B.
2012
Frequency-dependent attenuation as a potential indicator of oil saturation
J. appl. Geophys.
82
119
128
.

Rubino
J.G.
Müller
T.M.
Guarracino
L.
Milani
M.
Holliger
K.
2014
Seismoacoustic signatures of fracture connectivity
J. geophys. Res.
119
2252
2271
.

Schoenberg
M.
,
Helbig
K.
,
1997
.
Orthorhombic media: modeling elastic wave behavior in a vertically fractured earth
,
Geophysics
,
62
,
1954
1974
.

Shao
S.
,
Pyrak-Nolte
L.J.
,
2016
.
Wave propagation in isotropic media with two orthogonal fracture sets
,
Rock. Mech. Rock. Eng.
,
49
,
4033
4048
.

Shuai
D.
,
Stovas
A.
,
Wei
J.
,
Di
B.
,
2020
.
Frequency-dependent anisotropy due to two orthogonal sets of mesoscale fractures in porous media
,
Geophys. J. Int.
,
221
,
1450
1467
.

Tillotson
P.
,
Chapman
M.
,
Best
A.I.
,
Sothcott
J.
,
McCann
C.
,
Wang
S.
,
Li
X. Y.
,
2011
.
Observations of fluid-dependent shear-wave splitting in synthetic porous rocks with aligned penny-shaped fractures
,
Geophys. Prospect.
,
59
,
111
119
.

Tsvankin
I.
,
1997
.
Anisotropic parameters and P-wave velocity for orthorhombic media
,
Geophysics
,
62
,
1292
1309
.

White
J.E.
,
1975
.
Computed seismic speeds and attenuation in rocks with partial gas saturation
,
Geophysics
,
40
,
224
232
.

White
J.E.
,
Mikhaylova
N.G.
,
Lyakhovitskiv
F.M.
,
1975
.
Low-frequency seismic waves in fluid-saturated layered rocks
,
J. acoust. Soc. Am.
,
11
,
654
659
.

APPENDIX A.

Frequency-dependent anisotropic elastic constants for partially saturated porous rock with two sets of mesoscale fractures

Substituting eq. (6) into eq. (7), we can solve the system for the fluid pressure in each inclusion as a function of applied stress. The normal direction of the mth fracture with the orientation angles |${\theta _m}$| and |${\varphi _m}$| are defined by
(A1)
First, we define the notations:
(A2)
(A3)
(A4)
and
(A5)
where |${r_m}$| is the aspect ratio of the mth set of fractures.
The fluid pressure in the pores can be calculated by
(A6)
where
(A7)
Where the notation |${K_p} = \frac{{4\mu }}{{3{K_f}}}\ $|⁠. Likewise, the fluid pressures in mth sets of fractures are given by
(A8)
in which we define
(A9)
Following Chapman (2009), The effective frequency-dependent elastic stiffness tensor of a partially saturated rock with two sets of mesoscale fractures will be
(A10)
where |$C_{ij}^0$| is computed by Lamé parameter |$\lambda $| and |$\mu $|⁠. In terms of the normal direction of fractures in eq. (A1), the rotation matrix of mth set of fracture can be defined by
(A11)
The applied stress tensors are defined as
(A12)
(A13)
(A14)
and
(A15)
The forms of these tensors can isolate the various elements of fracture correction tensors easily. A set of fluid pressures with respect to the above stress tensors can be calculated:
(A16)
The elastic corrections of fractures can be computed from the fluid pressures following eq (A16):
(A17)
|$a_{ij}^m$| can be written in matrix form as follows
(A18)
and the elastic correction of the mth set of fracture with orientation angles |${\theta _m}$| and |${\varphi _m}$| can be given by
(A19)
where M is the Bond transformation matrix (Bond 1943).
The elastic correction for pores to the effective elastic parameters of the rock |$C_{ij}^p( \omega )$| is expressed by
(A20)
where
(A21)
and,
(A22)
The fluid pressure in pores is
(A23)
where
(A24)

APPENDIX B.

Vertical velocities and anisotropy parameters in anelastic ORT medium

According to Shuai et al. (2020), the vertical P- and S-wave velocities of partially saturated porous rock with orthogonal fractures are expressed by
(B1)
(B2)
The anisotropic parameters in |$[ {{x_1},\ {x_3}} ]$| symmetry plane are defined by
(B3)
(B4)
(B5)
The anisotropic parameters in |$[ {{x_2},\ {x_3}} ]$| symmetry plane are given by
(B6)
(B7)
(B8)
The anisotropic parameter in |$[ {{x_1},\ {x_2}} ]\ $| symmetry plane
(B9)
This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)