Abstract

Hydrocarbon reservoirs often contain partially gas-saturated rocks that have attracted the attention of exploration geophysicists and geologists for many years. Wave-induced fluid flow (WIFF) is an effective mechanism to quantify seismic wave dispersion and attenuation in partially gas-saturated rocks. In this study, we focus on the local fluid flow induced by variations in fluids in different regions and present a new model that describes seismic wave propagation in partially gas-saturated rocks, namely, the cylindrical patchy-saturation model. Because the seismic wave velocity and attenuation oscillate at high frequencies, it is not ideal for studying dispersion and attenuation caused by WIFF. To avoid the high-frequency oscillation in the cylindrical patchy-saturated model, we use an approximation to the Newman function instead of the full Newman function to calculate the effective bulk modulus. We then calculate the P-wave velocity and attenuation of the proposed model and interpret the lab-measured data. The proposed model is an alternative patchy-saturation model that can explain the problem of high-frequency oscillation and low-frequency attenuation.

1. Introduction

It is widely accepted that the presence of fluids plays a significant role in seismic wave attenuation and dispersion. Also, wave-induced fluid flow is an effective mechanism to quantify the attenuation and dispersion [1]. In this study, we concentrate on the attenuation and dispersion that occurs at the mesoscopic scale and use the WIFF mechanism to interpret the attenuation and dispersion of seismic waves.

WIFF that occurs on a scale much larger than the typical pore size but smaller than the wavelength is known as mesoscopic flow. Two important models that quantify the wave attenuation caused by mesoscopic flow have been proposed, namely, the dual-porosity model and the patchy-saturated model. The dual-porosity model focuses on the mesoscopic flow that occurs between soft pores and stiff pores, while the patchy-saturation model focuses on the mesoscopic flow that occurs between two different fluids, like water and gas.

The spatial variations in rock compliance in porous media may lead to local fluid flow, which can be described by the dual-porosity model. Cracked porous media have been analyzed by using the dual-porosity model (e.g., [2]). The dual-porosity model was proposed by Pride and Berryman [3, 4] and extended to Biot-Rayleigh theory by Ba et al. [5] by combining the conventional dual-porosity model and the bubble oscillation equation [6] to describe the seismic wave attenuation and dispersion characteristics. Additionally, several models have been developed by analyzing different inclusions or patches under the framework of the Biot-Rayleigh theory [711].

The spatial variations of fluid in porous media can lead to local fluid flow as described by the patchy-saturation model. The influence of partial saturation on sonic waves has been proven by laboratory observation [1214]. A theoretical model was first proposed by White [15] and White et al. [16] and has been developed since then by Dutta and Odé [17, 18] who calculated the attenuation and dispersion of patchy-saturated rocks based on Biot theory and by Johnson [19] who used branching functions to calculate bulk modulus attenuation and dispersion. Subsequently, the decoupling method proposed by Dutta and Odé has been extended to a full frequency range [2022] and an exact solution [22]. The patch distribution demonstrates important progress of the patchy-saturation model [2326]. It should be noted that the solutions derived by White [15] and Johnson [19] are approximate expressions for patchy-saturation rocks, and the exact solutions derived by Dutta and Odé [17, 18] and Vogelaar et al. [22] do not estimate P-wave propagation at high frequencies. Thus, it is necessary to find an exact solution without high-frequency oscillations.

While the presence of complex patch distributions is common, theories that consider a complex patch distribution are still challenging. Different patch distributions have a significant influence on seismic wave propagation implying that researchers should try out various patch distributions in their modeling effort. Thus, based on the discussion above and the premise that the presence of cracks and special inclusions makes the existence of a cylindrical patch possible, we present a cylindrical patchy-saturation model and derive its exact solution without high-frequency oscillations.

This paper is organized as follows: First, we propose the concept of a cylindrical patchy saturation model, and then, we derive an exact analytical solution for the cylindrical patchy-saturation model based on Biot theory [27]. Our results indicate that the P-wave velocity and attenuation curves oscillate at high frequencies for the method proposed by Vogelaar et al. [22]. To avoid these high-frequency oscillations, we propose an approximate solution for cylindrical patchy saturated rock and compute the seismic wave attenuation and dispersion curves. Finally, we conclude that the proposed model can explain the low-frequency attenuation very well when compared with laboratory data.

2. Exact Expression of Cylindrical Patchy-Saturated Rock

2.1. Concentric Cylinder Geometry

Since being first proposed by White [15] and White et al. [16], the patchy-saturated rock theory has undergone considerable development [1624, 2833]. Periodic slabs and concentric spheres are two well-known patch distributions. We present another patchy-saturation model with different patch distribution (cylindrical patchy-saturation model in Figure 1) model and analyze its properties. Unlike the traditional patchy-saturation model, it can be approximated by concentric cylinder geometry. As for other physical properties, the cylindrical patchy-saturation model and the traditional patchy-saturation model are almost identical.

Firstly, we solve the Biot [27, 34, 35] equation, whose cylindrical symmetric solutions are as follows:
(1)rP+Qur+ur+R+QUr+Ur=0,iωηϕ2κUu=rQur+ur+RUr+Ur.
P, Q, and R are the moduli and are expressed as
(2)P=Kd+4μ/3+αϕ2/β,Q=ϕαϕβ,R=ϕ2β.
Following Johnson [19], the solutions to the above equation are the linear combinations of r, r1, Bessel functions J1qr, and H1qr. Therefore, the general solution of the above equation could be expressed as
(3)ur=Ar+Br+Q+RFJ1qr+GH1qr,Ur=Ar+BrP+QFJ1qr+GH1qr,pr=2Q+RϕA+PRQ2ϕqFJ0qr+GH0qr,τr=2HμA2μr2B2μQ+RrFJ1qr+GH1qr,
with D=κ/ηϕ2PRQ2/H, q=iω/D, and H=P+2Q+R, where r is the radius of the cylinder; u, U, p, and τ are the solid displacement, fluid displacement, fluid pressure, and total stress, respectively; ϕ, η, and κ are the porosity, fluid viscosity, and permeability, respectively; J0 and H0 are the zero-order Bessel function and Neumann function, respectively; and J1 and H1 are the first-order Bessel function and Neumann function, respectively. The value of parameters P, Q, R, D, and q is different in two cylinders because the properties of the fluid and rock are various in the two regions.

2.2. Exact Analytic Solution

The exact analytic solution for the cylindrical patchy-saturation model can be found by applying boundary conditions to solve for the unknown parameters A, B, F, and G in equation (3). The boundary condition has been discussed in Dutta and Odé [17], Johnson [19], and Vogelaar et al. [22].

At the center of the cylinder, the solution for solid displacement and fluid displacement should be finite, which implies
(4)ua0=0,(5)Ua0=0.
At the interface of the cylinder, the solution of solid displacement, fluid displacement, fluid pressure, and total stress should be continuous, which implies
(6)uaa=uba,(7)Uaa=Uba,(8)paa=pba,(9)τaa=τba.
At the surface of the cylinder, the total stress is equal to the pressure generated by the seismic wave excitation, and the fluid does not flow out, which implies
(10)ubb=Ubb,(11)τbb=τ0.
At the center of the cylinder, as described by equations (4) and (5), there is no displacement in the solid or fluid, which implies Ba=0 and Ga=0. Since the parameters Ba and Ga are zero, after substituting them into equation (10), we obtain the relation between the Ga and Fa
(12)Gb=FbJ1qbbH1qbb,
For the sake of simplicity, we define auxiliary parameters fi as
(13)f0=J0qbaH0qbaJ1qbbH1qbb,f1=J1qbaH1qbaJ1qbbH1qbb.
Subtraction of the equations (6) and (7) gives
(14)FaJ1qaa=Fbf1HbHa.
Consequently, equation (6) be simplified as
(15)Aaa=Aba+Bba2+Fbf1N,
where
(16)N=Rb+QbRa+QaHbHa.
At the interface of two cylinders, the solid pressure and fluid pressure are equal, which indicates
(17)aμAaHaμ=aμAbHbμBbaFbf1N,(18)2AaRa+Qa+2AbRb+Qb=FbPbRbQb2qbf01h,
where the auxiliary parameter h is
(19)h=ηaηbDaDbqaqbJ0qaaJ1qaaf1f0.
Combining of equations (15), (17), and (18), the value of Bb then becomes
(20)Bb=b2Abg,
with
(21)g=sa1HbHa+2N2PbRbQb2f1f01qba11h.
Finally, at the boundary of the outer cylinder, the equation (11) indicates
(22)2AbKBGb2μb2Bb=τ0.
Substituting the solid displacement of the outer cylinder, ub=Abb1g, into equation (22) and with effective modulus b/2ubτ0, the quasi-static bulk modulus Kω could be derived as
(23)Kω=KBGb+4/3μg1g.

We know that the dispersion and attenuation of the quasi-static bulk modulus are controlled by parameter g, which indicates that gas saturation sa plays an important role.

2.3. Plane Wave Velocity and Q

The complex P-wave velocity and inverse quality factor 1/Q of cylindrical patchy-saturated rock can be calculated following Carcione [36] as
(24)Vp=ReKω+4μ/3ρ,(25)Q1=ImVp2ReVp2,
where ρ=1ϕρs+ϕsbρfb+saρfa is the density and μ is the shear modulus.

Figure 2 shows the P-wave velocity and attenuation for the cylindrical patchy-saturation model in the frequency range of 100 to 106 Hz on a logarithmic scale. As seen in Figure 2(a), there is an obvious dispersion of the P-wave velocity. It can also be seen from Figure 2(b) that the P-wave attenuation curve has a peak, which is caused by mesoscopic flow. We note, however, that the attenuation and dispersion curves oscillate at high frequencies, which is a problem that is necessary to address.

3. Exact Expression for Cylindrical Patchy-Saturated Rock without High-Frequency Oscillation

High-frequency oscillations result in the improper simulation of seismic wave propagation in underground media; therefore, it is necessary to propose a solution that avoids high-frequency oscillations. Through numerical calculations, we observe that the value of the auxiliary parameter fi is unstable in high frequencies. After a series of numerical experiments, we find that if we use the Newman function approximation at high frequencies, the value of the auxiliary function fi tends to be stable rather than oscillatory. This then motivates us to derive an approximation for the general solution using the Newman function approximation Hnx2/πxsinxπ/4nπ/2 to resolve the problem. This is a simple and computationally inexpensive method. The solution for low and high frequencies can be expressed analytically as follows:
(26)ur=Ar+Br+Q+RFJ1qr+GH1qr,Ur=Ar+BrP+QFJ1qr+GH1qr,pr=2Q+RϕA+PRQ2ϕqFJ0qr+GH0qr,τr=2HμA2μr2B2μQ+RrFJ1qr+GH1qr,ωω0,ur=Ar+Br+Q+RFJ1qrG2πqrsinqr+π4,Ur=Ar+BrP+QFJ1qrG2πqrsinqr+π4,pr=2Q+RϕA+PRQ2ϕqFJ0qr+G2πqrsinqrπ4,τr=2HμA2μr2B2μQ+RrFJ1qrG2πqrsinqr+π4,ωω0,
where ω0 is the reference frequency, which could be determined when the misfit between the Neumann function and its approximation is less than 102.
Similarly, applying boundary conditions to equation (26), the approximate quasi-static bulk modulus Kω is
(27)Kω=KBGb+4/3μg1g,
with
(28)g=sa1HbHa+2N2PbRbQb2f1f01k2ba11h,h=ηaηbDaDbqaqbJ0qaaJ1qaaf1f0,f0=J0qbaH0qbaJ1qbbH1qbb,f1=J1qbaH1qbaJ1qbbH1qbbω<ω0,f0=J0qbabasinqbaπ/4sinqbb3π/4J1qbb,f1=J0qbabasinqba3π/4sinqbb3π/4J1qbbω>ω0.

The main difference between the approximate bulk modulus and its exact value is that the auxiliary parameter fi is calculated in different ways. Using the approximate quasi-static bulk modulus above, we can obtain the P-wave velocity and attenuation curve without high-frequency oscillations, which will be discussed below.

4. Numerical Calculation Results and Analysis

In this section, we compute the seismic wave velocity and attenuation for the cylindrical patchy-saturation model using our proposed theory and method. The model parameters used for the calculations are shown in Table 1.

Figure 3 gives the P-wave velocity and attenuation curves for the proposed model in the frequency range of 100 to 106 Hz. The value of the P-wave velocity agrees well with the low-frequency limit (static limit) and the high-frequency limit (no-flow limit) (the low- and high-frequency limits are derived in Appendix A). Note that the P-wave velocity and attenuation curves of the proposed model do not oscillate at high frequencies; therefore, the proposed method is an effective way to avoid high-frequency oscillations for cylindrical patchy-saturated rock.

4.1. Comparison between the Proposed and Previous Models

Figure 4 gives the P-wave velocity and attenuation curves for the White, Johnson, and modified Vogelaar et al. models and the proposed model for the frequency range of 100 to 106 Hz. The descriptions of these models can be found in Appendices B, C, and D. The P-wave velocity value from all the models coincides at the low and high frequencies (we assume that the parameter gas saturation sa of all models is equal). As shown in Figure 4, we find that the error in the P-wave velocity of the Johnson [19] model is the largest, where the error is defined as the difference between the approximate and exact model. The maximum underestimation of the velocity from the Johnson model is about 4 m/s, while that from the White model is about 2 m/s. This is compatible with Vogelaar et al.’s conclusion and implies that the solution of the cylindrical patchy-saturation model derived by the same method is the most exact way to simulate seismic wave propagation. The attenuation peak of the proposed model moves toward low frequencies, as compared with the modified Vogelaar model. The displaced attenuation peak is similar to that of the attenuation peak under a different cylinder radius.

4.2. Wave Propagation Frequency-Dependent Parameter

To describe more clearly the propagation characteristics of a P wave in the cylindrical patchy-saturation model, we compute the P-wave velocity and attenuation curves for different cylinder radii (Figure 5), porosity (Figure 6), and fluid viscosity (Figure 7).

Figure 5 shows the P-wave velocity and attenuation curves for different cylinder radii in the frequency range of 100 to 106 Hz. Figure 5 indicates that the attenuation and dispersion curves of the proposed model do not oscillate at high frequencies. As shown in Figure 5, the cylinder radius affects the seismic attenuation and dispersion caused by mesoscopic flow. The attenuation peak caused by mesoscopic flow moves toward low frequencies with increasing cylinder radius, and the amplitude of 1/Q value decreases with increasing radius. These results are consistent with the model proposed by Ba et al. [5] and Pride and Berryman [3]. In summary, the cylinder radius is a frequency-dependent parameter.

Figure 6 gives the P-wave velocity and attenuation curves for different porosities in the frequency range of 100 to 106 Hz. As shown in Figure 6, the attenuation peak induced by mesoscopic flow does not move with the increasing porosity while the amplitude of the attenuation value (1/Q) decreases slightly. If the value of porosity increases by 5%, the amplitude of the 1/Q value increases by about 0.08. This property is different from other frequency-dependent parameters, which may help us to determine the relative value of porosity in the reservoir.

Figure 7 gives the P-wave velocity and attenuation curves for different values of fluid viscosity in the frequency range of 100 to 106 Hz. It is shown from Figure 7 that the attenuation peak induced by mesoscopic flow moves toward low frequencies with the increasing fluid viscosity while the amplitude of the 1/Q value does not change. In summary, fluid viscosity is a frequency-dependent parameter.

5. Interpreting Lab-Measured Data

In this section, we prove that the concept of cylindrical patchy-saturation model can interpret the P-wave attenuation induced by WIFF in low frequency. Yao and Han [37] and Yao et al. [38] published a recent low-frequency measurement. We will use a sample of Berea sandstone for the analysis of the P-wave attenuation in the range of 100 to 102: the porosity and permeability values used here are 22.3% and 4.021013m2(407 mD), respectively, and the viscosity of water is 2.98103Pas. If the rock sample can be considered as a homogeneous medium, we can apply the proposed model to the measured data to invert for the radius of the concentric cylinder under the condition of gas saturation.

Figure 8 shows P-wave attenuation from low-frequency lab-measured data and theoretical simulation. We can infer that the inner cylinder radius a is about 0.13 m, and the outer radius b is about 0.40 m under the condition of gas saturation, Sa=1/9. In this example, we demonstrate that the cylindrical patchy-saturation model effectively explains the low-frequency attenuation.

6. Conclusion

We present a cylindrical patchy-saturation model and derive the quasi-static bulk modulus to study the properties of seismic wave propagation. The P-wave velocity and attenuation, however, oscillate at high frequencies, which has a deleterious effect on the conclusions about velocity dispersion and attenuation. To avoid this phenomenon, after studying the exact equation for cylindrical patchy-saturated rocks, we find the factors that cause high-frequency oscillation and derive a theory of cylindrical patchy-saturated rocks to avoid the high-frequency oscillation. Compared with the exact equation of cylindrical patchy-saturation model derived by Vogelaar et al. [22], the proposed model is consistent, though with some misfit, with exact equation at low frequencies and avoids the high-frequency oscillations. This indicates that the proposed model is better at simulating the effect of mesoscopic flow caused by a seismic wave. In summary, the proposed method is an effective solution for avoiding high-frequency oscillations for cylindrical patchy-saturated rock and simulating seismic wave propagation in underground media.

Appendix

A. Low-Frequency and High-Frequency Limits

If ω0, the approximation of the Bessel functions is limn0J0n=1, limn0J1n=n/2, and limn0H1n=2/πn. When we use the approximate forms of the Bessel functions, the low-frequency limit of the auxiliary parameter h can be calculated as
(A.1)limω0h=ηaηbDaDbsasb.
Substituting the auxiliary parameter g into equation (A-1), the expression of the auxiliary parameter g becomes
(A.2)limω0g=sa1HbHaκϕ2sbHbN2saDbηb+sbDaηa.
After substituting equation (A-2) into modulus, the low-frequency modulus KBGW can be rearranged in an accessible form
(A.3)KBGW=Kd+α2αϕKs+ϕKw1,
where 1/Kw=sa/Kfa+sb/Kfb. The modulus KBGWis known as the Biot-Gassmann-Wood modulus.
If ω, the high-frequency limit of the auxiliary parameter g can be calculated as
(A.4)limω0g=sa1HbHa.
Therefore, the high-frequency modulus KBGH can be written as
(A.5)1KBGH+4/3μ=saKBGHa+4/3μ+sbKBGHb+4/3μ=saHa+sbHb.

The modulus KBGH is known as the Biot-Gassmann-Hill modulus.

B. White Model

As proposed by White [15], the complex bulk modulus K as a function of angular frequencyωis given by
(B.1)K=KBGH1KBGHW,
where
(B.2)W=3a2RaRbQa+QbiωZa+Zbb3,Ra=KaKd1Kd/Ks3Kb+4μKb3Ka+4μ+4μKaKbSa,Rb=KbKd1Kd/Ks3Ka+4μKb3Ka+4μ+4μKaKbSa,Za=ηaaκ1e2α1aα1a1+α1a+1e2α1a,Zb=ηbaκα2b+1+α2b1e2α2baα2b+1α2a1α2b1α1a+1e2α2ba,α1=iωηaκKEa1/2,α2=iωηbκKEb1/2,KEa=1Kfa1Ka/Ks1Kd/KsϕKa1Kfa/KsKAa,KEb=1Kfb1Kb/Ks1Kd/KsϕKb1Kfb/KsKAb,KAa=ϕKfa+1ϕKsKdKs21,KAb=ϕKfb+1ϕKsKdKs21,Qa=1Kd/KsKAaKa,Qb=1Kd/KsKAbKb.

a and b denote the two different regions; η and Kf are the fluid viscosity and bulk modulus of fluid, respectively; ϕ is the porosity; κ is the permeability; Ks is the bulk modulus of the solid grains; μ is the shear modulus.

C. Johnson Model

As proposed by Johnson [19], the complex bulk modulus K as a function of angular frequency ω is given by
(C.1)Kω=KBGHKBGHKBGW1ζ+ζ1iωτ/ζ2,
where
(C.2)τ=KBGHKBGWKBGHG,ζ=KBGHKBGW2KBGWτT,T=KBGWϕ230κb33ηbgb2+5ηaηbgagb3ηaga2a515ηbgbgbgaa3b2+5gb3ηbgb2ηb+ηagaa2b33ηbgb2b5,ga=1Kd/Ks1/Kw1/Kfa1Kd/KsϕKd/Ks+ϕKd/Kw,gb=1Kd/Ks1/Kw1/Kfb1Kd/KsϕKd/Ks+ϕKd/Kw,ΔpfPe=Rb+QbKa+4μ/3Ra+QaKb+4μ/3ϕSaKaKb+4μ/3+ϕSbKbKa+4μ/3,G=ΔpfPe23a2b3iqiω0.5,q=iω/D,D=κKBGHηaDa+ηbDb,Da=κηaϕ2PaRaQa2Ha,Db=κηbϕ2PbRbQb2Hb,Kw=SaKfa+SbKfb1,P=Kd+4μ/3+αϕ2β,Q=ϕαϕβ,R=ϕ2β.

KBGH and KBGW are discussed in Appendix A. a and b denote inner sphere radius and outer sphere radius; η and Kf are the fluid viscosity and bulk modulus of fluid, respectively; ϕ is the porosity; κ is the permeability; Ks is the bulk modulus of the solid grains.

D. Modified Vogelaar et al.’s Model

Following the idea proposed in this paper, we make an approximation to the general solution of Vogelaar et al. [22] model using the spherical Newman function approximation nnx1/πxsinxnπ/2. The approximate quasi-static bulk modulus Kω then becomes
(D.1)Kω=KBGb+4/3μg1g,
where
(D.2)g=sa1HbHa+3N2PbRbQb2f1f01k2ba11h,h=ηaηbDaDbqaqbj0qaaj1qaaf1f0,f0=j0qban0qbaj1qbbn1qbb,f1=j1qban1qbaj1qbbn1qbbω<ω0,f0=j0qbabasinqbaπ/2sinqbbπj1qbb,f1=j0qbabasinqbaπsinqbbπj1qbbω>ω0,D=κηϕ2PRQ2H,q=iωD,H=P+2Q+R.

η is the fluid viscosity; ϕ is the porosity; κ is the permeability; j0 and n0 are the zero-order spherical Bessel function and spherical Neumann function, respectively; j1 and n1 are the first-order spherical Bessel function and spherical Neumann function, respectively.

Data Availability

Previously reported P-wave attenuation data were used to support this study and are available at https://library.seg.org/doi/epub/10.1190/geo2013-0410.1. These prior studies (and datasets) are cited at relevant places within the text as references Yao et al. [37].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

We would like to acknowledge the sponsorship of the National Natural Science Foundation of China (42174139, 41974119, 42030103) and the Science Foundation from Innovation and Technology Support Program for Young Scientists in Colleges of Shandong province and Ministry of Science and Technology of China.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).