*Article* **A New Anelasticity Model for Wave Propagation in Partially Saturated Rocks**

**Chunfang Wu 1, Jing Ba 1,\*, Xiaoqin Zhong 2, José M. Carcione 1,3, Lin Zhang <sup>1</sup> and Chuantong Ruan 1,4**

<sup>1</sup> School of Earth Sciences and Engineering, Hohai University, Nanjing 211100, China;

	- China; zhongxiaoqin\_cq@petrochina.com.cn

<sup>4</sup> School of Mathematics and Statistics, Zhoukou Normal University, Zhoukou 466001, China

**\*** Correspondence: jba@hhu.edu.cn

**Abstract:** Elastic wave propagation in partially saturated reservoir rocks induces fluid flow in multi-scale pore spaces, leading to wave anelasticity (velocity dispersion and attenuation). The propagation characteristics cannot be described by a single-scale flow-induced dissipation mechanism. To overcome this problem, we combine the White patchy-saturation theory and the squirt flow model to obtain a new anelasticity theory for wave propagation. We consider a tight sandstone Qingyang area, Ordos Basin, and perform ultrasonic measurements at partial saturation and different confining pressures, where the rock properties are obtained at full-gas saturation. The comparison between the experimental data and the theoretical results yields a fairly good agreement, indicating the efficacy of the new theory.

**Keywords:** partial saturation; patchy saturation; squirt flow; P-wave velocity dispersion and attenuation; anelasticity; ultrasonic measurements

#### **1. Introduction**

Seismic waves induce fluid flow and anelasticity (the wave-velocity dispersion and dissipation factor) in rocks saturated with immiscible fluids [1–8]. The level of anelasticity depends on the in situ pressure, fluid content and type, and pore structure. This subject is highly relevant to petroleum exploration and production.

WIFF (wave-induced fluid flow) occurs at various spatial scales that can be categorized as macroscopic, mesoscopic, and microscopic [9]. The first is the wavelength-scale equilibration process occurring between the peaks and troughs of a P-wave, while the mesoscopic length is much larger than the typical pore size but smaller than the wavelength. The microscopic scale is of the same order of magnitude as the pore and grain sizes.

The macroscopic mechanism has been discussed by Biot [10–12] and is often referred to as the Biot relaxation peak (usually at kHz dominant frequencies). The basic assumptions are that the rock frame is homogeneous and isotropic, and the relative motion between the grains and the pore fluid is governed by Darcy's law. Local fluid flow on meso- and micro-scales are neglected, and consequently, the Biot peak cannot explain the observed wave anelasticity at all frequencies [13].

Partial saturation leads to fluid heterogeneity at the mesoscopic scale and the pressure difference between fluid phases causes wave dissipation at low frequencies [9,14–19]. White [20] proposed the first patchy-saturation model (the White model, spherical pockets). Dutta and Odé [21] reformulated this model by using the Biot theory, while Johnson [22] generalized it to patches of arbitrary geometry by using a branch function. Liu et al. [23] analyzed the effect of the fluid properties.

Moreover, dissimilar pores, with different shapes (micro-fractures and intergranular pores) and/or orientations, also cause mesoscopic pressure gradients and squirt flow,

**Citation:** Wu, C.; Ba, J.; Zhong, X.; Carcione, J.M.; Zhang, L.; Ruan, C. A New Anelasticity Model for Wave Propagation in Partially Saturated Rocks. *Energies* **2021**, *14*, 7619. https://doi.org/10.3390/ en14227619

Academic Editor: Eugen Rusu

Received: 5 October 2021 Accepted: 6 November 2021 Published: 15 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

resulting in dissipation. At the pore level, dissipation can be described with squirt flow models [24,25]. Dvorkin and Nur [26] unified the Biot and squirt flows and proposed the BISQ model (Biot/squirt), which describes anelasticity at some frequency ranges. However, the low-frequency P-wave velocity prediction from the BISQ model is smaller than the Gassmann velocity [27], while it is consistent with the Biot one at high frequencies. Dvorkin et al. [28] extended the BISQ model to partially saturated rocks by incorporating the Wood equation [29] and proposed that the squirt flow length can be related to water saturation. Dvorkin et al. [30] reformulated the BISQ model to achieve consistency with the Gassmann velocity at the low-frequency limit. However, the P-wave velocity obtained with this model is higher than the theoretical high limit at high frequencies (when all the cracks are closed, and the P-wave velocity value is determined by the Biot model) [31]. Wu et al. [31] proposed a reformulated modified frame squirt flow model (MFS) to solve the problem.

Mavko and Jizba [32] introduced a modified frame to estimate the high-frequency unrelaxed dry rock shear and bulk moduli (M-J model), where cracks are saturated and the stiff pores are drained. To obtain the wet rock properties from the M-J model, Gurevich et al. [33] used the pressure relaxation method of Murphy et al. [34]. The model can be applied in a broad frequency range. Wu et al. [31] presented a reformulated modified frame squirt flow model, but although the prediction is acceptable at ultrasonic frequencies, Pride et al. [35] showed that the attenuation is significantly underestimated at seismic frequencies because the mesoscopic mechanism is not taken into account.

Wave anelasticity is mainly due to the effect of multi-scale fluid flow [36–40]. Rubino and Holliger [41] studied the problem at the micro and meso scales, analyzing the effects of the pore aspect ratio, while Li et al. [42] studied wave velocity in fractured poroelastic media saturated with immiscible fluids. Recently, Sun [43] proposed a model which considers the three loss mechanisms, i.e., the Biot, squirt flow, and mesoscopic relaxation peaks, in the framework of a double-porosity theory. This low-velocity limit does not honor Gassmann velocity.

We briefly review the propagation models at different scales, and propose a new one, based on the White theory and a reformulated modified frame squirt flow model (see Figure 1). Based on the numerical examples, the wave propagation characteristics of the new model and the effect of permeability and the outer diameter of the patch are analyzed. The P-wave velocity and attenuation with varying saturations are measured at different confining pressures. The crack properties and squirt flow length are obtained from the experimental data. The new model is applied to ultrasonic measurements performed on tight sandstone from the Qingyang area of the Ordos basin. The comparison between the experimental data and the theoretical results are made, so as to verify the capability of the new model in the description of those wave properties.

**Figure 1.** A new model based on a reformulated modified frame squirt flow (MFS) model with the White theory. The effects of squirt flow occurring between soft and stiff pores in the water-saturated host medium are incorporated by using an equivalent host medium of a modified solid (the MFS model). On the other hand, the White theory describes the anelasticity due to the patchy saturation of the immiscible fluid mixture.

#### **2. Model**

#### *2.1. Patchy-Saturation (White) Model*

White [20] proposed a patchy saturation model, by considering flow in a concentric spherical model where the inner sphere is saturated with one fluid type (gas), and the outer shell is saturated with a liquid (water), where the frame is assumed to be homogeneous. Let *a* and *b* be the inner and outer diameters, such as (*b* > *a*), and the gas saturation is *Sg* = *a*3/*b*3. Dutta and Odé [21] modified the White model based on the Biot model, and obtained the following wet rock bulk and shear moduli:

$$K^\*(\omega) = \frac{K\_{\infty}}{1 - \mathcal{W}K\_{\infty}},\tag{1}$$

$$\mathbf{G}^\*(\omega) = \mathbf{G}\_{dry\_{\mathcal{Y}}} \tag{2}$$

respectively, where *K*<sup>∞</sup> is the bulk modulus at the high-frequency limit, *Gdry* is the dry rock shear modulus, and *W* is a complex function of porosity, permeability, fluid viscosity, etc. (see Appendix A in Carcione et al. [44], and the Section 2.3).

#### *2.2. Squirt Flow Model*

Th flow between microcrack and grain contacts back and forth to stiff (equant) pores induces dissipation even for a single saturating fluid. The microcracks are incorporated into an effective rock skeleton, containing only stiff pores.

The reformulated modified frame squirt flow model considers both the squirt and Biot flows. According to Dvorkin et al. [30] and a boundary condition given by Gurevich et al. [33], the modified bulk modulus is (Wu et al. [31]):

$$K\_{ms} = K\_{msd} + \frac{\alpha\_c \,^2 F\_c}{\phi\_c} \left[ 1 - \frac{2J\_1(\lambda R)}{\lambda R J\_0(\lambda R)} \right],\tag{3}$$

where *Kmsd* = 1/*K*<sup>0</sup> <sup>−</sup> 1/*Khp* <sup>+</sup> 1/*Kdry*−<sup>1</sup> , *K*<sup>0</sup> is the bulk modulus of the mineral mixture, *Kdry* is the dry rock bulk modulus, *Fc* = 1/*Kf l* + 1/(*φcQc*) −<sup>1</sup> , *φ<sup>c</sup>* is the microcrack porosity, *α<sup>c</sup>* = 1 − *Kmsd*/*K*0, *Qc* = *K*0/(*αc*−*φc*), *R* is characteristic squirt flow length, *ω* is the angular frequency, *λ*<sup>2</sup> = *iωηφc*/*κ* 1/*Kf l* + 1/(*φcQc*) , *η* is the fluid viscosity, *κ* is the permeability, *Kf l* is the bulk modulus of fluid, and *J*<sup>0</sup> and *J*<sup>1</sup> are the zero- and first-order Bessel functions, respectively.

The modified dry-frame bulk and shear moduli are (Wu et al. [31]):

$$\frac{1}{K\_{md}} = \frac{1}{K\_{ms}} + \frac{1}{K\_{hp}} - \frac{1}{K\_0} \, \text{} \tag{4}$$

$$\frac{1}{G\_{md}} = \frac{1}{G\_{dry}} - \frac{4}{15} \left( \frac{1}{K\_{dry}} - \frac{1}{K\_{md}} \right) \tag{5}$$

respectively, where *Khp* is the high-pressure modulus [33]. The P-wave phase velocity and attenuation can then be obtained according to Toksöz and Johnston [45] as

$$V\_{\text{ph}P1,2} = \frac{1}{\text{Re}(X\_{1,2})},\ a\_{1,2} = \omega \text{Im}(X\_{1,2}),\tag{6}$$

where

$$X\_{1,2} = \sqrt{Y\_{1,2}}\ \ Y\_{1,2} = -\frac{B}{2A} \pm \sqrt{\left(\frac{B}{2A}\right)^2 - \frac{C}{A}}\ \ A = \frac{\phi FM\_{dry}}{\rho\_{22}^2}\ \ .$$

$$B = \frac{F\left(2\alpha\_{md} - \phi - \phi\frac{\rho\_{11}}{\rho\_{22}}\right) - \left(M\_{dry} + F\frac{\alpha\_{md}^2}{\phi}\right)\left(1 + \frac{\rho\_a}{\rho\_{22}} + i\frac{\omega\_c}{\omega}\right)}{\rho\_{22}},$$

$$C = \frac{\rho\_{11}}{\rho\_{22}} + \left(1 + \frac{\rho\_{11}}{\rho\_{22}}\right)\left(\frac{\rho\_a}{\rho\_{22}} + i\frac{\omega\_c}{\omega}\right), \rho\_{11} = (1 - \phi)\rho\_{s}, \rho\_{22} = \phi\rho\_{f1},$$

$$F = \left(1/K\_{fl} + (\alpha\_{md} - \phi)/(\phi K\_0)\right)^{-1},$$

where *ρ<sup>a</sup>* is the additional coupling density, *ω<sup>c</sup>* = *ηφ*/(*κρf l*) is the characteristic frequency, *αmd* = 1 − *Kmd*/*K*0, *φ* is porosity, *Mdry* is the uniaxial modulus of the rock skeleton under drained conditions, and *ρ<sup>s</sup>* and *ρf l* are the mineral density and fluid density, respectively.

#### *2.3. Patch-Saturation and Squirt Flow Models Combined*

The White model assumes a uniform rock skeleton and that the area outside the inclusion is fully saturated with water. Therefore, the modified dry rock moduli (4) and (5) are used in the White model, thus combining the micro and meso descriptions of anelasticity. If subindices 1 and 2 refer to the gas-inclusion region and host medium (water), respectively, we have the wet rock moduli

$$K(\omega) = \frac{K\_{\infty}}{1 - WK\_{\infty}}\tag{7}$$

$$G(\omega) = G\_{\text{md}\prime} \tag{8}$$

where

$$K\_{\infty} = \frac{K\_{\rm G2}(3K\_{\rm G1} + 4G\_{md}) + 4G\_{md}(K\_{\rm G1} - K\_{\rm G2})S\_{\rm g}}{(3K\_{\rm G1} + 4G\_{md}) - 3(K\_{\rm G1} - K\_{\rm G2})S\_{\rm g}} \tag{9}$$

$$\mathcal{W} = \frac{3i\alpha\kappa (R\_1 - R\_2)(F\_1 - F\_2)}{b^3 \omega (\eta\_1 Z\_1 - \eta\_2 Z\_2)}. \tag{10}$$

Moreover,

$$K\_{\rm G1} = \frac{K\_0 - K\_{\rm md} + \phi K\_{\rm md} \left( K\_0 / K\_{f11} - 1 \right)}{1 - \phi - K\_{\rm md} / K\_0 + \phi K\_0 / K\_{f11}} \tag{11}$$

$$K\_{\rm G2} = \frac{K\_0 - K\_{md} + \phi K\_{md} \left(K\_0 / K\_{fl2} - 1\right)}{1 - \phi - K\_{md} / K\_0 + \phi K\_0 / K\_{fl2}} \tag{12}$$

are Gassmann moduli, where *Kf l*<sup>1</sup> and *Kf l*<sup>2</sup> are fluid moduli,

$$R\_1 = \frac{(K\_{G1} - K\_{md})(3K\_{G2} + 4G\_{md})}{(1 - K\_{md}/K\_0)\left[K\_{G2}(3K\_{G1} + 4G\_{md}) + 4G\_{md}(K\_{G1} - K\_{G2})S\_{\wp}\right]} \tag{13}$$

$$R\_2 = \frac{(K\_{G2} - K\_{md})(3K\_{G1} + 4G\_{md})}{(1 - K\_{md}/K\_0)\left[K\_{G2}(3K\_{G1} + 4G\_{md}) + 4G\_{md}(K\_{G1} - K\_{G2})S\_{\mathcal{S}}\right]} \tag{14}$$

$$F\_1 = \frac{(1 - K\_{md}/K\_0)K\_{A1}}{K\_{G1}} \tag{15}$$

$$F\_2 = \frac{(1 - K\_{md}/K\_0)K\_{A2}}{K\_{G2}} \tag{16}$$

$$Z\_1 = \frac{1 - \exp(-2\gamma\_1 a)}{(\gamma\_1 a - 1) + (\gamma\_1 a + 1)\exp(-2\gamma\_1 a)}\tag{17}$$

$$Z\_2 = \frac{(\gamma\_2 b + 1) + (\gamma\_2 b - 1) \exp[-2\gamma\_2(b - a)]}{(\gamma\_2 b + 1)(\gamma\_2 a - 1) - (\gamma\_2 b - 1)(\gamma\_2 a + 1) - \exp[-2\gamma\_2(b - a)]} \tag{18}$$

$$
\gamma\_1 = \sqrt{i\omega \eta\_1 / \kappa K\_{E1}} \tag{19}
$$

$$
\gamma\_2 = \sqrt{i\omega\eta\_2/\kappa K\_{E2}}\tag{20}
$$

where *η*<sup>1</sup> and *η*<sup>2</sup> are fluid viscosities, and

$$K\_{E1} = \left[1 - \frac{K\_{f11}(1 - K\_{G1}/K\_0)(1 - K\_{\rm md}/K\_0)}{\phi K\_{G1}\left(1 - K\_{f11}/K\_0\right)}\right] K\_{A1} \tag{21}$$

$$K\_{E2} = \left[1 - \frac{K\_{f12}(1 - K\_{G2}/K\_0)(1 - K\_{\rm md}/K\_0)}{\phi K\_{G2}\left(1 - K\_{f12}/K\_0\right)}\right] K\_{A2} \tag{22}$$

$$\frac{1}{K\_{A1}} = \left(\frac{\phi}{K\_{fl1}} + \frac{1-\phi}{K\_0} - \frac{K\_{md}}{K\_0^2}\right) \tag{2.3}$$

$$\frac{1}{K\_{A2}} = \left(\frac{\phi}{K\_{f12}} + \frac{1-\phi}{K\_0} - \frac{K\_{md}}{K\_0^2}\right). \tag{24}$$

According to Wood [29], the effective bulk modulus of the gas-water mixture can be calculated from

$$\frac{1}{\mathcal{K}\_{fl}} = \frac{\mathcal{S}\_{\mathcal{S}}}{\mathcal{K}\_{fl1}} + \frac{\mathcal{S}\_{w}}{\mathcal{K}\_{fl2}}\tag{25}$$

where *Sw* is the water saturation.

Finally, the P-wave phase velocity and attenuation are

$$V\_p = \sqrt{\frac{\text{Re}(\mathcal{K}(\omega) + 4G(\omega)/3)}{\rho}},\tag{26}$$

$$Q\_p^{-1} = \frac{\text{Im}(K(\omega) + 4G(\omega)/3)}{\text{Re}(K(\omega) + 4G(\omega)/3)},\tag{27}$$

respectively, where *ρ* = (1 − *φ*)*ρ<sup>s</sup>* + *φ* - *Sgρ*<sup>1</sup> + *Swρ*<sup>2</sup> is bulk density, and *ρ*<sup>1</sup> and *ρ*<sup>2</sup> are the fluid densities.

#### *2.4. Results*

The MFS model is directly applied in partially saturated reservoir rocks, where the gas–water mixture is obtained with the Wood equation (there are no gas pockets), and the properties are listed in Table 1. The numerical examples of the characteristics of wave prorogation by the proposed model are shown in Figure 2, and the effects of permeability and the outer diameter of the patch on the wave velocity and attenuation are shown in Figures 3 and 4, respectively.



**Figure 2.** P-wave velocity (**a**) and attenuation (**b**) of the present and MFS models. The number between parentheses indicates water saturation.

**Figure 3.** P-wave velocity (**a**) and attenuation (**b**) of the present model as a function of water saturation and various permeabilities.

**Figure 4.** P-wave velocity (**a**) and attenuation (**b**) of the present model as a function of water saturation for different values of the outer diameter.

Figure 2 compares the P-wave velocity (a) and attenuation (b) of the present model with those of the MFS model, where the number between parentheses indicates water saturation. The velocities coincide at low frequencies and increase with saturation, with those of the present model higher at high frequencies. Two inflection points are clearly observed, corresponding to the mesoscopic and squirt flow attenuation peaks when the saturation is 80%, the first being the stronger point. The attenuation of the present model is higher than that of the MFS one.

Figure 3 shows the effect of permeability, where we can see that attenuation has a maximum at a given saturation which increases with permeability. Figure 4 displays the same quantities as a function of saturation for different outer diameters (*b*). The velocity increases with *b*, and the attenuation decreases.

#### **3. Ultrasonic Data**

#### *3.1. Rock Specimen and Experiment*

A tight sandstone sample (S2-9) from the Qingyang area, Ordos Basin, was tested. The sample was processed into a cylinder with a diameter of 25.2 mm and a length of 50 mm. An aluminum standard with the same shape and size was processed corresponding to the specimen. The sample was composed of quartz, feldspar, and interstitial materials (mainly carbonate minerals and clay), and its porosity was 8.85%. A thin section is shown in Figure 5. The experimental set-up consisted of a pulse generator, a temperature control unit, a confining pressure control unit, a pore pressure control unit, and an ultrasonic wave test unit [46,47].

**Figure 5.** Thin section image of a tight sandstone.

The piezoelectric ultrasonic wave transducers were glued to the top and bottom of the sample, sealed with a rubber sleeve. An electrical pulse was applied to the source transducer to generate the ultrasonic P-waves. A digital oscilloscope was used to display and record the waveforms from the receiver. The temperature, pore, and confining pressures were controlled by the appropriate units [48]. The pore pressure was 15 MPa, the effective pressures were 5, 15, 25, 35 and 45 MPa, the temperature was 20 ◦C, and the waveforms were recorded after we maintained the experimental conditions for half an hour. For the partial gas–water saturation tests, the samples were first saturated with water by using the vacuum pressure saturation method and then placed in an oven to vary the saturation. The approach of Ba et al. [19] was adopted to quantify the fluid content. The sample was tested around six different water saturation conditions, 0%, 20%, 40%, 60%, 80%, and 100%. The wave velocities were obtained from the travel times and the spectral-ratio method was used to obtain the dissipation factor.

#### *3.2. Experimental Results*

Figure 6 shows the velocity as a function of water saturation and effective pressure. As expected, the P-wave velocity increases with water saturation and pressure, approaching a linear trend at high pressures [49,50], since microcracks close. In the partially saturated rock, the rock pore spaces contain air (with a lower bulk modulus and a lower P-wave velocity) and water (with a higher bulk modulus and a higher P-wave velocity). With the increase in water saturation, the volume ratio of water increases and that of air decreases while the rock skeleton stays unchanged. Generally, the P-wave velocity increases with the water saturation. The influence of effective pressure on the stiff pores is small and can be neglected [50–52]. The S-wave velocity also increases with effective pressure, but decreases as saturation increases, due to the density effect.

**Figure 6.** P-(**a**) and S-(**b**) wave velocities as a function of effective pressure at different water saturations.

The spectral ratio method is applied to calculate the dissipation factor [53,54]. We have

$$\ln\left(\frac{A\_1(f)}{A\_2(f)}\right) = -\frac{\pi x}{QV}f + \ln\frac{G\_1(x)}{G\_2(x)}\tag{28}$$

where *f* is the frequency, *A*1(*f*) and *A*2(*f*) are the amplitude spectra of the rock sample and standard, respectively, *Q* is the quality factor, *x* is the propagation distance, *V* is the wave velocity, and *G*1(*x*) and *G*2(*x*) are the sample and standard geometrical factors, respectively. As shown in Figure 7, attenuation decreases with effective pressure. Its behavior versus saturation is similar to that of Figure 4. The attenuation variations with respect to effective pressure and saturation are similar to those of the sandstone samples analyzed by Pang et al. [55] and Amalokwu et al. [56].

**Figure 7.** P-wave attenuation as a function of effective pressure and saturation.

#### *3.3. Crack Parameters and Squirt Flow Length*

The parameters of the present model can be obtained from the experimental data. They involve the skeleton bulk and shear moduli at different pressures, the dry rock bulk modulus with microcracks closed, the microcrack porosity, the squirt flow length, etc. The following steps are considered.

(1) The dry rock bulk and shear moduli are calculated from the velocities as

$$K\_{dry} = \left(V\_{pd}^2 - \frac{4}{3}V\_{sd}^2\right)\rho\tag{29}$$

$$G\_{dry} = V\_{sd}^2 \rho \tag{30}$$

where *ρ* = *ρs*(1 − *φ*), *Vpd* and *Vsd* are the P- and S-wave velocities of full gas saturation, respectively.


The DZ model is applied to calculate the microcrack density (Figure 8a) and porosity (Figure 8b) based on the experimental data at different effective pressures. Their variations are more significant in the low-pressure range. As pressure increases, both quantities decrease. The microcrack density and porosity decreases, which can be attributed to the closure of microcracks [57,58].

**Figure 8.** Microcrack density (**a**) and porosity (**b**) as a function of effective pressure.

The dry rock density of sample S2-9 is 2410 kg/m3, and the bulk modulus of the mineral mixture is 39 GPa. The fluid properties are determined from the empirical equations of Batzle and Wang [59]. Figure 9 displays the P-wave velocity as a function of the effective pressure, where the squirt flow lengths are obtained by matching the theoretical results to the experimental data. It shows that the sample can be characterized by a constant squirt flow length at different pressures [31]. The characteristic length of sample S2-9 is 0.45 mm. This quantity is not so relevant to the pressure and it can be considered as an intrinsic rock property [26].

**Figure 9.** P-wave velocity as a function of the effective pressure compared to the experimental data. Results at different squirt flow lengths are shown.

#### **4. Comparison between Theory and Experiment**

#### *4.1. Effect of Saturation*

The present model is used to calculate the P-wave velocity and attenuation of sample S2-9 at 5 MPa effective pressure. The dry rock bulk modulus is 20.5 GPa, the Poisson ratio is 0.15, the permeability is 0.177 mD, the outer diameter is 0.12 mm, the high-pressure modulus is 23 GPa, and the fluid properties are listed in Table 1. Figure 10 displays the results for different models. The Gassmann–Hill curve does not consider the wave-induced fluid flow. The MFS curve coincides with the Gassmann–Wood curve at low saturations (Figure 10a), while the present model has a velocity of the order of the Gassmann–Hill curve. When saturation increases, the rock is stiffened by the microscopic fluid flow, resulting in a velocity increase.

**Figure 10.** P-wave velocity (**a**) and attenuation (**b**) as a function of water saturation at 5 MPa. The open circles correspond to the experimental data.

The influence of fluid flow is determined by fluid pressure gradients at the interface between different fluid phases, where the fast P-wave converts to the slow (diffusive) Biot wave (mesoscopic loss) [15,60]. Attenuation has a maximum at a given saturation due to the mesoscopic loss mechanism, absent in the MFS model, whereas at full gas or water saturation, the results are similar. The attenuation curves are similar to those of the sandstone samples analyzed by Amalokwu et al. [56].

Compared with the simplified models, the present model provides a good match between the theory and the ultrasonic data for a tight sandstone, mainly the P-wave velocity

as a function of saturation at an effective pressure of 5 MPa, showing the effectiveness of the squirt flow model combined with the White theory. However, attenuation is underestimated by the model due to the fact that the spatial variations in mineral grain and porosity are not considered.

#### *4.2. Effect of Effective Pressure*

In this example, the outer diameters at effective pressures of 15, 25, 35, and 45 MPa are 0.14, 0.16, 0.18 and 0.2 mm, respectively. Figures 11–14 display the P-wave velocity and attenuation as a function of water saturation at these pressures. The overall trend is similar to that at 5 MPa. As pressure increases, the MFS velocities approach the Gassmann–Wood velocities, and the P-wave velocity predictions from the new model increases; however, the attenuation decreases where most microcracks close, and the squirt flow effects are inhibited. Therefore, the characteristics of wave propagation from the new model are similar to those of the experimental data in Figures 11–14.

**Figure 11.** P-wave velocity (**a**) and attenuation (**b**) as a function of water saturation at 15 MPa. The open circles correspond to the experimental data.

**Figure 12.** P-wave velocity (**a**) and attenuation (**b**) as a function of water saturation at 25 MPa. The open circles correspond to the experimental data.

**Figure 13.** P-wave velocity (**a**) and attenuation (**b**) as a function of water saturation at 35 MPa. The open circles correspond to the experimental data.

**Figure 14.** P-wave velocity (**a**) and attenuation (**b**) as a function of water saturation at 45 MPa. The open circles correspond to the experimental data.

#### *4.3. Crossplots*

Figure 15 shows crossplots of the measured and theoretical velocities, showing a good agreement.

Attenuation crossplots are displayed in Figure 16. The attenuation prediction from the present model is less than the experimental one, particularly at low effective pressures. This can be due to the fact that the model considers only the mesoscopic loss caused by partial saturation. Additional attenuation may be due to the presence of many minerals, and the microcrack content, shape, and distribution [15,51,61].

**Figure 15.** Crossplots of the measured and theoretical velocities at different water saturations (**a**) and effective pressures (**b**).

**Figure 16.** Crossplots of the measured and theoretical dissipation factors at different water saturations (**a**) and effective pressures (**b**).

#### **5. Conclusions**

We combine the reformulated modified frame squirt flow model and White's mesoscopic loss theory (spherical gas pockets) to develop a new model for describing wave anelasticity in partially saturated rocks. The microcrack properties and characteristic squirt flow lengths are obtained from experimental data at different effective pressures and water saturations. We compare the results with those of simplified models, showing that the present model provides a good match between the theory and the ultrasonic data for a tight sandstone, mainly the P-wave velocity as a function of saturation and pressure. Attenuation is underestimated by the model due to the fact that mesoscopic loss (fast P- to slow P-wave conversion) due to spatial variations in mineral grain and porosity are not considered. The new model can be used to predict the characteristics of wave propagation in partially saturated tight sandstones, mainly the P-wave velocity. Moreover, a better description at high frequencies (from tens of kHz) should consider the Biot attenuation peak. The generalization of the model will be the task of a future paper.

Due to the complex microstructures and fabric heterogeneity of tight sandstone, the proposed model cannot fully describe the experimental measurement data at low effective pressures. The theories cannot perfectly match real rocks, and there might be errors/defects in the experiment measurements. In the related engineering applications of hydrocarbon reservoir exploration, the methods of big data analytics and machine learning may be applied in combination with the theoretical models, so as to improve the applicability of the model and the accuracy of reservoir property prediction or interpretation.

**Author Contributions:** Conceptualization, J.B. and C.W.; Data curation, J.B. and X.Z.; Formal analysis, C.W., J.B., X.Z., J.M.C., L.Z. and C.R.; Funding acquisition, J.B.; Investigation, C.W., J.B., J.M.C. and L.Z.; Methodology, C.W., J.B. and J.M.C.; Project administration, J.B. and X.Z.; Resources, J.B. and X.Z.; Supervision, J.B.; Validation, C.W. and J.M.C.; Writing—original draft, C.W. and J.B.; Writing—review & editing, J.B., J.M.C. and L.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Natural Science Foundation of China, grant number 41974123 and 42174161, the Jiangsu Province Outstanding Youth Fund Project, grant number BK20200021, the National Science and Technology Major Project of China, grant number 2017ZX05069- 002, and Sinopec Key Laboratory of Geophysics.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data relevant with this study can be accessed by contacting the corresponding author.

**Acknowledgments:** This work is supported by the National Natural Science Foundation of China (Grant no. 41974123, 42174161), the Jiangsu Province Outstanding Youth Fund Project (Grant no. BK20200021), the National Science and Technology Major Project of China (Grant no. 2017ZX05069- 002) and Sinopec Key Laboratory of Geophysics.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Symbols**

**Table A1.** List of symbols.


**Table A1.** *Cont.*


Microcrack Porosity Estimation at Different Effective Pressures


$$\frac{1}{K\_{stiff}^{MT}} = \frac{1}{K\_0} \left( 1 + \frac{\phi\_s}{1 - \phi\_s} P \right) \tag{A1}$$

$$\frac{1}{G\_{stiff}^{MT}} = \frac{1}{G\_0} \left( 1 + \frac{\phi\_s}{1 - \phi\_s} \mathcal{Q} \right),\tag{A2}$$

respectively; *G*<sup>0</sup> is the shear modulus and *φ<sup>s</sup>* is the stiff porosity.

$$P = \frac{(1 - \upsilon)}{6(1 - 2\upsilon)} \times \frac{4(1 + \upsilon) + 2\gamma^2(7 - 2\upsilon) - \left[3(1 + 4\upsilon) + 12\gamma^2(2 - \gamma)\right]g}{2\gamma^2 + (1 - 4\gamma^2)g + (\gamma^2 - 1)(1 + \gamma)g^2} \tag{A3}$$

$$\begin{array}{c} Q = & \frac{4(\gamma^2 - 1)(1 - \gamma)}{15\left\{8(\gamma - 1) + 2\gamma^2(3 - 4v) + \left[(7 - 8v) - 4\gamma^2(1 - 2v)\right]g\right\}}\\ & \times \left\{\frac{8(1 - v) + 2\gamma^2(3 + 4v) + \left[(8v - 1) - 4\gamma^2(5 + 2v)\right]g + 6\left(\gamma^2 - 1\right)(1 + v)g^2}{2\gamma^2 + (1 - 4\gamma^2)g + (\gamma^2 - 1)(1 + \gamma)g^2} \\ & - 3\left[\frac{8(v - 1) + 2\gamma^2(5 - 4v) + \left[3(1 - 2v) + 6\gamma^2(v - 1)\right]g}{-2\gamma^2 + ((2 - \gamma) + \gamma^2(1 + v))g}\right] \end{array} \tag{A4}$$

where *γ* is the spheroidal aspect ratio, and *υ* is the Poisson ratio of the grains, i.e.,

$$g = \begin{cases} \upsilon = (3\mathcal{K}\_0 - 2\mathcal{G}\_0) / (6\mathcal{K}\_0 + 2\mathcal{G}\_0), \\ \frac{\gamma}{(1-\gamma^2)^{3/2}} \left( \arccos \gamma - \gamma \sqrt{1-\gamma^2} \right) (\gamma < 1) \\\ \frac{\gamma}{(1-\gamma^2)^{3/2}} \left( \gamma \sqrt{1-\gamma^2} - \arccos \eta \right) (\gamma > 1) \end{cases} \tag{A5}$$

Microcracks are included into the host material by neglecting the interactions between cracks and pores. The effective moduli (host with cracks) are

$$\frac{1}{K\_{eff}^{MT}} = \frac{1}{K\_{stiff}^{MT}} \left( 1 + \frac{16\left(1 - \left(\upsilon\_{stiff}^{MT}\right)^2\right) \Gamma}{9\left(1 - 2\upsilon\_{stiff}^{MT}\right)} \right) \tag{A6}$$

$$\frac{1}{G\_{eff}^{MT}} = \frac{1}{G\_{stiff}^{MT}} \left( 1 + \frac{32 \left( 1 - \upsilon\_{stiff}^{MT} \right) \left( 5 - \upsilon\_{stiff}^{MT} \right) \Gamma}{45 \left( 2 - \upsilon\_{stiff}^{MT} \right)} \right) \tag{A7}$$

where *υMT sti f f* = 3*KMT sti f f* − <sup>2</sup>*μMT sti f f* / 6*KMT sti f f* + <sup>2</sup>*μMT sti f f* and Γ is the microcrack density. When all the microcracks close at high pressures, a least square method is used to obtain the optimal aspect ratio of the stiff pores by using Equations (A1) and (A2).



$$
\Gamma = \Gamma^i e^{-p/\not p} \tag{A8}
$$

where Γ*<sup>i</sup>* is the initial value when the effective pressure is zero, *p* is the effective pressure, and *p*ˆ is a constant.


$$\gamma\_p^i = \frac{3}{4\pi} \int\_{\Gamma^i}^{\Gamma} \frac{\left(1/K(\Gamma) - 1/K\_{eff}^{hp}\right)}{\Gamma} \frac{dp}{d\Gamma} d\Gamma \tag{A9}$$

where *K*(Γ) is the effective bulk modulus which can be obtained from Equation (A1). Substituting Equation (A8) into (A9), we obtain

$$\gamma\_p^i = \frac{3}{4\pi} \int\_{\Gamma}^{\Gamma^i} \frac{\left(1/K(\Gamma) - 1/K\_{eff}^{hp}\right)\mathcal{P}}{\Gamma^2} d\Gamma \tag{A10}$$

and by integrating Equation (A10) from Γ to Γ*<sup>i</sup>* ,

$$\gamma\_p^i = \frac{4\dot{\rho} \left[1 - \left(\upsilon\_{eff}^{hp}\right)^2 \ln\left(\frac{\Gamma^i}{\Gamma}\right)\right]}{3\pi K\_{eff}^{hp} \left[1 - 2\upsilon\_{eff}^{hp}\right]} \tag{A11}$$

where *υ hp eff* is the effective Poisson ratio at high pressures, i.e., *υ hp eff* = 3*Khp eff* <sup>−</sup> <sup>2</sup>*Ghp eff* / 6*Khp eff* <sup>+</sup> <sup>2</sup>*Ghp eff* .

Combining Equations (A8) and (A11), the relation between the minimum initial aspect ratio and the effective pressure can be obtained as

$$\gamma\_p^i = \frac{4\left[1 - \left(\upsilon\_{eff}^{hp}\right)^2\right]p}{\pi E\_{eff}^{hp}}\tag{A12}$$

where *Ehp eff* <sup>=</sup> <sup>3</sup>*Khp eff* & 1 − 2*υ hp eff* ' is the effective Young modulus at high pressures. The cumulative microcrack density decreases with pressure. If pressure changes from zero to *dp*, the corresponding reduction of the cumulative microcrack density is *d*Γ. When the pressure increment is small enough, it can be considered that the decrease of microcrack density is mainly due to the closure of microcracks with an aspect ratio less than the minimum initial aspect ratio. David and Zimmerman [52] related the microcrack porosity and density as

$$\phi\_{\mathbb{C}} = \frac{4\pi\gamma}{3}\Gamma \tag{A13}$$

Therefore, the microcrack properties can be obtained from the acoustic wave velocities as a function of the effective pressure.

#### **References**

