An Exact Expression for the Effective Bulk Modulus for Acoustic Wave Propagation in Cylindrical Patchy-Saturation Rocks


 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.


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 [7][8][9][10][11].
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 [12][13][14]. 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 [20][21][22] and an exact solution [22]. The patch distribution demonstrates important progress of the patchy-saturation model [23][24][25][26]. It should be noted that the solutions derived by White [15] and Johnson [19] are approximate expressions for patchysaturation 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 highfrequency 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.

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 patchysaturated rock theory has undergone considerable development [16][17][18][19][20][21][22][23][24][28][29][30][31][32][33]. 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 patchysaturation model, it can be approximated by concentric cylinder geometry. As for other physical properties, the cylindrical patchy-saturation model and the traditional patchysaturation model are almost identical. Firstly, we solve the Biot [27,34,35] equation, whose cylindrical symmetric solutions are as follows: ð1Þ P, Q, and R are the moduli and are expressed as Following Johnson [19], the solutions to the above equation are the linear combinations of r, r −1 , Bessel functions J 1 ðqrÞ, and H 1 ðqrÞ. Therefore, the general solution of the above equation could be expressed as with D = ðκ/ηϕ 2 ÞððPR − Q 2 Þ/HÞ, q = ffiffiffiffiffiffiffiffiffiffiffiffi ffi −iω/D p , 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; J 0 and H 0 are the zero-order Bessel function and Neumann function, respectively; and J 1 and H 1 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.

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].  Lithosphere At the center of the cylinder, the solution for solid displacement and fluid displacement should be finite, which implies At the interface of the cylinder, the solution of solid displacement, fluid displacement, fluid pressure, and total stress should be continuous, which implies 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 At the center of the cylinder, as described by equations (4) and (5), there is no displacement in the solid or fluid, which implies B a = 0 and G a = 0. Since the parameters B a and G a are zero, after substituting them into equation (10), we obtain the relation between the G a and F a For the sake of simplicity, we define auxiliary parameters f i as Subtraction of the equations (6) and (7) gives Consequently, equation (6) be simplified as where At the interface of two cylinders, the solid pressure and fluid pressure are equal, which indicates where the auxiliary parameter h is Combining of equations (15), (17), and (18), the value of B b then becomes with Finally, at the boundary of the outer cylinder, the equation (11) indicates Substituting the solid displacement of the outer cylinder, uðbÞ = A b bð1 − gÞ, into equation (22) and with effective modulus −ðb/2uðbÞÞτ 0 , the quasi-static bulk modulus KðωÞ could be derived as We know that the dispersion and attenuation of the quasi-static bulk modulus are controlled by parameter g, which indicates that gas saturation s a plays an important role.

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 where ρ = ð1 − ϕÞρ s + ϕ½s b ρ f b + s a ρ f a 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 10 0 to 10 6 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.

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 f i 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 f i tends to be stable rather than oscillatory. This then motivates us to derive an approximation for the general solution using the Newman function approximation H n ðxÞ ⟶ ffiffiffiffiffiffiffiffiffiffiffiffiffi ð2/πxÞ p sin ðx − ðπ/4Þ − ðnπ/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: where ω 0 is the reference frequency, which could be determined when the misfit between the Neumann function and its approximation is less than 10 -2 . Similarly, applying boundary conditions to equation (26), the approximate quasi-static bulk modulus KðωÞ is The main difference between the approximate bulk modulus and its exact value is that the auxiliary parameter f i 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.

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 10 0 to 10 6 Hz. The value of the P-wave velocity agrees well with the low-frequency limit (static limit) and the highfrequency limit (no-flow limit) (the low-and highfrequency 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 pro-posed method is an effective way to avoid high-frequency oscillations for cylindrical patchy-saturated rock. 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 10 0 to 10 6 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 s a 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     5 Lithosphere 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.

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 10 0 to 10 6 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  6 Lithosphere 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 frequencydependent parameter. Figure 6 gives the P-wave velocity and attenuation curves for different porosities in the frequency range of 10 0 to 10 6 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 7 Lithosphere of 10 0 to 10 6 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.

Interpreting Lab-Measured Data
In this section, we prove that the concept of cylindrical patchy-saturation model can interpret the P-wave attenua-tion 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 10 0 to 10 2 : the porosity and permeability values used here are 22.3% and 4:02 ⋅ 10 −13 m 2 (407 mD), respectively, and the viscosity of water is 2:98 ⋅ 10 -3 Pa ⋅ s. 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, S a = 1/9. In this example, we demonstrate that the cylindrical patchy-saturation model effectively explains the low-frequency attenuation.

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, Theoretical data Lab-measured data Figure 8: Low-frequency lab-measured attenuation (circles) for a sample of Berea sandstone and the theoretical P-wave attenuation (line). 8 Lithosphere 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 patchysaturated rock and simulating seismic wave propagation in underground media.
η is the fluid viscosity; ϕ is the porosity; κ is the permeability; j 0 and n 0 are the zero-order spherical Bessel function and spherical Neumann function, respectively; j 1 and n 1 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.