Tight sandstones have low porosity and permeability and strong heterogeneities with microcracks, resulting in small wave impedance contrasts with the surrounding rock and weak fluid-induced seismic effects, which make the seismic characterization for fluid detection and identification difficult. For this purpose, we propose a reformulated modified frame squirt-flow (MFS) model to describe wave attenuation and velocity dispersion. The squirt-flow length (R) is an important parameter of the model, and, at present, no direct method has been reported to determine it. We obtain the crack properties and R based on the DZ (David-Zimmerman) model and MFS model, and how these properties affect the wave propagation, considering ultrasonic experimental data of the Sichuan Basin. The new model can be useful in practical applications related to exploration areas.

Tight (oil and gas) reservoirs are becoming important in seismic exploration [13], with tight sandstones playing a dominant role in China, since they are widely distributed. Unlike conventional reservoirs, tight sandstone ones have a complex geological origin, low porosity and permeability, and strong heterogeneities with a complex pore structure. In recent years, a series of studies have been performed related to the characterization of these rocks [48].

The pore structure of the reservoirs affects the oil and gas distribution [9]. Cracks not only affect the elasticity of the rocks, but also the fluid flow [1013]. Zimmerman [14] presented an empirical relation between the compressibility and effective pressure of sandstones and David and Zimmerman [15] extended the method to predict the crack distribution from dry-rock wave velocities (DZ model) on the basis of the MT [16] and DEM [17] models. Deng et al. [18] investigated the pore structure of conventional and tight sandstones based on the DZ model, and Wei et al. [19] studied the influence of effective pressure and porosity on the crack properties.

Wave propagation in earth crust is often accompanied by significant attenuation and dispersion [2023], and many theoretical models have been proposed [11, 2427]. With new developments and analyses on field measured data, it has been realized that local fluid flow induced by heterogeneities is the main cause of dissipation [28]. Based on the Biot theory [2931], the BISQ model [32] introduced a new property called “squirt-flow length” to describe the squirt flow and combined the effects at macroscopic (Biot’s) and microscopic scales. The fast P-wave velocity dispersion predicted by the model is significantly higher than that of the Biot model. However, the fast P-wave velocity is lower than the Gassmann velocity [33] at low frequencies and is consistent with the high-frequency one predicted by the Biot theory. Mavko and Jizba [34] proposed an equation (M-J model) based on a modified rock skeleton to compute the elastic modulus of the unrelaxed wet-rock skeleton at high frequencies. The P-wave velocity predicted by the M-J model is greater than that of the Biot model at high frequencies and also greater than the Gassmann velocity at low frequencies. This equation is not applicable to the case of gas-saturated or dry rocks. In view of the limitations of the M-J model, Gurevich et al. [35] extended it based on the theory of Murphy et al. [36], which can be applied at all frequencies. Li et al. [37] pointed out that the characteristic frequency of squirt flow is related to the crack aspect ratio, fluid viscosity, and bulk moduli of the minerals. Carcione and Gurevich [38] unified the squirt-flow and Biot theories and performed numerical simulations based on the Zener mechanical model [39].

Dvorkin et al. [40] obtained the complex modulus of fluid-saturated rocks based on the BISQ mode by considering a one-dimensional radial flow. This theory is consistent with the Gassmann velocity at the low-frequency limit. However, the model does not consider the influence of the Biot flow. In addition, the P-wave velocity obtained with this model is higher than the theoretical maximum value at high frequencies (when all the cracks are closed, and the P-wave velocity value is determined by the Biot model) [41]. To overcome these problems, Wu et al. [41] presented the MFS model.

The parameter R is crucial in the squirt-flow model, but to date, it has not been obtained by direct experiments, and this restricts the use of the theory [42]. Dvorkin et al. [43] showed that this parameter represents the radius of a cylinder with its axis parallel to the direction of wave propagation. Based on ultrasonic experimental data of Best [44], Marketos and Best [45] found that the relation between R and viscosity follows a power law. Since the BISQ model was proposed, other studies have been performed, which mainly focused on extending the model [40, 41, 4648], considering anisotropy [28, 49] and wave simulation [5052]. However, the approach about how to determine R remains unclear.

First, we use ultrasonic experimental data of tight sandstones and apply the DZ model to analyze the crack distribution. Then, R is obtained with the MFS model. We analyze the relations between R and crack density, crack aspect ratio, permeability, and other properties. A semiempirical formula of R for tight sandstones is proposed.

Tutuncu et al. [53] showed that there is a close contact between adjacent grains in tight sandstone, and a large number of cracks are formed at the edges of grains. We select twelve tight sandstone samples, collected from the Upper Triassic Xujiahe Formation in the Guang’an gasfield of the central Sichuan Basin. This sandstone is one of the greatest potential strata in this basin, whose reservoirs are characterized by low porosity and permeability, well-developed cracks, and high-water saturation. The samples are mainly lithic quartz sandstone, with a small amount of siliceous quartz sandstone, which is composed of quartz, feldspar, lithic debris, and cements. The porosity ranges from 3% to 14%, and the bulk modulus of the mineral is 39 GPa. The rock samples were processed into cylinders with a diameter of 25 mm and a height of 25 to 50 mm, and both ends of the sample were polished. The experimental set-up consists of a pulse generator, a temperature control unit, a confining pressure control unit, a pore pressure control unit, and an acoustic wave test unit [54, 55].

For the tests with full gas saturation (nitrogen), the samples were sealed with a rubber sleeve. The pore pressure was 10 MPa, and effective pressures of 5, 10, 15, 20, 25, 30, and 35 MPa were applied. The transmitted waveforms were recorded at 80°C, by maintaining the experimental conditions for half an hour. Then, the samples were saturated with brine (146,377 ppm NaCl) by vacuum pumping pressurization. The velocities at different effective pressures were obtained from the first arrivals of the extracted waveforms. The errors of the velocity measurement are estimated according to Yurikov et al. [56], and they are in range of 0.242-0.465%. The rock properties are given in Table 1.

A thin-section analysis of the sample is shown in Figure 1, where blue-dyed resin indicates porosity and microcracks/contacts are observed.

Figures 2 and 3 show the S- and P-wave velocities, respectively, as a function of the effective pressure at full gas and water saturations. The P- and S-wave velocities are higher in the second case, increasing with pressure. Some water-saturated tight samples show a greater S velocity than the gas-saturated state, which is related to the heterogeneous microcracks and microstructures. For tight rocks, the stiffening effect of local fluid flow on rock skeleton, which is induced by S-waves at high frequencies, may also lead to a higher S-wave velocity at the water-saturated state [55]. In the low-pressure range, the rate of velocity variation with pressure is significant and then decreases with increasing pressure (e.g., [18, 57], and [58]). The main reason is that the cracks close when effective pressure increases. The influence of pressure on the stiff pores is small and can be neglected [15, 59], so that the velocities as a function of the pressure can be used to obtain the crack properties.

The different pressures cause the porosity changes, and the change of crack porosity is dominant. The porosities of the sandstones as a function of the effective pressure are shown in Figure 4, where we can see that porosity decreases with increasing pressure. The porosity variation is more significant in the low-pressure range, which is caused by the closure of cracks when pressure increases. At high pressures, the porosity variation diminishes.

3.1. Methodology

Based on the experimental data, parameters such as crack density and aspect ratio are obtained with the DZ model [15]. These parameters can be used to characterize the pore distribution. The specific steps are given in Appendix A, and the flow chart is shown in Figure 5.

3.2. Crack Properties

The crack properties can be obtained by following the methodology presented in Section 3.1.

Figure 6 shows the crack density as a function of the effective pressure for twelve tight sandstones. The crack density gradually decreases with increasing pressure. When the pressure is low, the variation in crack density is significant.

The crack aspect ratio is not a constant, but a continuous distribution within a certain range [18]. However, the main crack aspect ratio is constant at a certain pressure, which is the aspect ratio corresponding to the peak of the curve [60]. The main crack aspect ratio of sample GAR7 for different pressures is shown in Figure 7, where we observe that the crack porosity decreases with increasing pressure, more pronounced at low pressures [18, 42, 61, 62]. The main reason is the closure of cracks.

4.1. MFS Model

The pore space can be classified into compliant pores (cracks and grain contacts) and stiff pores (intergranular voids), where the latter correspond to the main porosity [34, 59, 63]. A modified frame squirt-flow model is proposed according to the characteristics of the microscopic pore structure. The cracks are incorporated into an effective rock skeleton as shown in Figure 8, containing only stiff pores.

The MFS model is based on the theory of Dvorkin et al. [40], by applying the boundary condition of Gurevich et al. [35] (the boundary pressure at the contact between cracks and stiff pores is constant). The P-wave velocity predicted by this model is consistent with that from the Gassmann equation at low frequencies and approaches the high limit value (the rock without compliant pores) at high frequencies [41].

The modified-frame bulk modulus is
(1)Kms=Kmsd+αc2Fcϕc12J1λRλRJ0λR,
where Kmsd=1/K01/Khp+1/Kdry1, Kdry is the dry-rock bulk modulus, Fc=1/Kfl+1/ϕcQc1, ϕc is crack porosity, αc=1Kmsd/K0, Qc=K0/αmdϕc, αmd=1Kmd/K0 is the poroelasticity coefficient, J0 and J1 are the zero- and first-order Bessel functions, respectively, λ2=iωηϕc/κ1/Kfl+1/ϕcQc, ω is angular frequency, η is the fluid viscosity, κ is permeability, and Kfl is the fluid bulk modulus.
Then, the modified dry-rock bulk and shear moduli are
(2)1Kmd=1Kms+1Khp1K0,1μmd=1μdry4151Kdry1Kmd,
respectively, where Khp is the bulk modulus of the dry rock when all the cracks are closed, which can be obtained by fitting the dry-rock velocity with effective pressure [35], and μdry is the dry-rock shear modulus. The P-wave velocity and attenuation of the fluid-saturated rock can be obtained according to Toksöz and Johnston [64]:
(3)VphP1,2=1ReX1,2,a1,2=ωImX1,2,
where
(4)X1,2=Y1,2,Y1,2=B2A±B2A2CA,A=ϕFMdryρ22,B=F2αmdϕϕρ1/ρ2Mdry+Fαmd2/ϕ1+ρa/ρ2+iωc/ωρ2,C=ρ1ρ2+1+ρ1ρ2ρaρ2+iωcω,ρ1=1ϕρs,ρ2=ϕρfl,F=1Kfl+αmdϕϕK01,
where ρa is the additional coupling density, ωc=ηϕ/κρfl is the characteristic frequency, αmd=1Kmd/K0, ϕ the is porosity, Mdry is the uniaxial modulus of the rock skeleton at drained conditions, ρs is the mineral density, and ρfl is fluid density.

Figure 9 shows the P-wave velocity of sample GAR7 as a function of frequency and different pressures, where R is computed with the least-square method (see Table 2 in Appendix B), by minimizing square of difference between the experimental data and P-wave velocity predicted by the MFS model. From 5 to 35 MPa, R takes the following values: 0.055, 0.051, 0.049, 0.044, 0.041, 0.038, and 0.03 (in mm). The fluid properties of the fluid are those at the measurement conditions according to Batzle and Wang [65]. The bulk and shear moduli of the dry rock are obtained from the velocity and density. The P-wave velocity dispersion decreases with increasing pressure.

4.2. Squirt-Flow Length

Figure 10 shows the P-wave velocity as a function of the effective pressure, where R is obtained as described above (see Table 3 in Appendix B).

It is seen that the P-wave velocity predicted by the model increases with the squirt-flow length. For most of the samples (GAR11, GAR6, GA8, GAR7, GAR12, GA1, GAR8, and GA2), R decreases as pressure increases. The main reason is that the cracks with smaller aspect ratios tend to close. In samples GA3, GAR1, GA10, and GA6, R increases with increasing pressure. The pressure is not the only factor affecting the squirt-flow length. These tight samples show different characteristics of fabric heterogeneity. The mineral composition, pore/microcrack shape/scale, grain scale, and their distributions all effect on the acoustic wave velocities. We adopt a single R in the model predictions for each sample, which in fact assumes a unique microcrack scale for the proposed idealized rock model. However, for those samples with a higher degree of fabric heterogeneity (i.e., containing more types of different microcracks, specifically with different scales), the model prediction may deviate more from the experimental measurement than those with a lower heterogeneity degree. Figure 10 shows the velocity as a function of frequency for the different samples at different pressures and the optimal squirt-flow length. It shows that each sample is approximately characterized by a constant R at different pressures. Thus, R can be considered as an intrinsic property of the rock [32]. The appropriate squirt-flow lengths of the twelve samples are 0.015, 0.035, 0.02, 0.023, 0.15, 0.12, 0.22, 0.045, 0.25, 0.35, 0.25, and 0.28 (in mm) and mainly occur between 10 and 25 MPa. In order to improve the estimation, R can be calculated by using the P- and S-wave velocities in this pressure range.

The comparison between the MFS model results and experimental data is shown in Figure 11, based on R reported in Figure 10. In samples GAR11, GAR6, GA8, GAR7, GAR12, GA1, GAR8, and GA2, at low frequencies, the predicted P-wave velocity is lower than the experimental data. With increasing pressure, the prediction gradually approaches the experimental data. For samples GA3, GAR1, GA10, and GA6, the result is higher than the experimental data in the low-pressure range. As pressure increases, the result approaches the experimental data.

4.3. Squirt-Flow Length and Crack Properties

Based on the experimental data, the factors influencing the squirt flow are analyzed with the MFS model. Figures 1214 show R as a function of crack density, main crack aspect ratio, and permeability, respectively, for the twelve rock samples. The results show that R increases with increasing crack density and permeability and decreases with increasing crack aspect ratio. Cheng et al. [66] showed that the crack fraction and its radius increase with total porosity. Consequently, a higher crack fraction and flow length enhance the effect of the local fluid flow induced by the waves.

Crack density and aspect ratio are the main parameters to represent the effect of cracks, since they affect R [67]. Basically, the analysis shows that there is an optimal value of R within a certain pressure range for a given sample. Permeability also affects R (see Figure 15).

The characteristic frequency of the squirt flow and the critical diffusion length can be used to determine the squirt-flow length. The first is given by [35, 68, 69], and [37].
(5)fc=γ3K0η.
Mavko and Mukerji [70] report a critical diffusion length.
(6)lc=κKflηfc.
According to Equations (5) and (6), we have
(7)lc=κKflK0γ3.

When the diffusion length is less than the critical one, the pore-fluid pressure has sufficient time to reach an equilibrium (relaxed state); otherwise, the pore fluid pressure cannot reach equilibrium, and the regime is unrelaxed.

According to Pride [71] and Dvorkin and Nur [32] and the previous analysis, we find that R is proportional to permeability and crack density and inversely proportional to the main crack aspect ratio.

In this work, we assume a semiempirical equation.
(8)R=aκKflK0γ3Γ+b.
The parameters a and b can be obtained by fitting experimental data and can be different for different lithologies and facies. Based on the twelve tight sandstones, we obtain
(9)R=0.03856κKflK0γ3Γ+63.42.

Equation (9) is represented in Figure 16, where we can see that R increases with the increase of crack density and permeability, when the crack density and permeability increase, and decreases with the increase of main crack aspect ratio, and the characteristic frequency of squirt flow moves to low frequencies (relaxed state) [32, 35, 40, 67, 71].

R is considered an intrinsic property, and Equation (9) provides a quantitative way for its determination. The property cannot be measured directly. However, it can be obtained by experiment using a suitable rock-physics model.

The sandstones in Guang’an area have low porosity and permeability and strong heterogeneities with microcracks and can be considered as a typical type of tight sandstone. We have proposed a rock-physics model to estimate the squirt-flow length (R) of tight sandstones, based on experiment data. The crack density and aspect ratio and frame permeability are obtained with effective-medium theories. We establish a semiempirical equation relating these quantities to the squirt-flow length, based on ultrasonic experimental data. The results show that R increases with increasing crack density and permeability and decreases with increasing crack aspect ratio. The methodology provides a new method for the calculation of R in tight sandstones, which can also be applied to other lithologies, e.g., carbonates. In a seismic exploration context, even an approximate estimation of R from seismic velocities measured at a single frequency is useful, as this allows one to infer derived parameters such as crack density, aspect ratio, and permeability, by using rock physics models.

Appendix

Figure 9 shows the P-wave velocity of sample GAR7 as a function of frequency and different pressures, where R is computed with the least-square method (see Table 2 in Appendix B, which gives the minimum square of difference between P-wave velocity prediction of the MFS model and experimental data of sample (GAR7) under different effective pressures), by minimizing square of difference between the experimental data and P-wave velocity predicted by the MFS model.

Figure 10 shows the P-wave velocity as a function of the effective pressure, where R is obtained as described above (see Table 3 in Appendix B, which gives the minimum square of difference between experimental data and P-wave velocity predicted by the MFS model and the relative difference between experimental data and model prediction for each rock in Figure 10).

A. Estimation of the Rock Properties

Step 1. Calculate the aspect ratio of the stiff pores. Based on the MT model, the quantitative relationship between elastic moduli and stiff porosity is established. The effective moduli are
(A.1)1KstiffMT=1K01+ϕs1ϕsP,(A.2)1μstiffMT=1μ01+ϕs1ϕsQ,
where KstiffMT and μstiffMT are the bulk and shear moduli of the rock only containing stiff pores, K0 and μ0 are the bulk and shear moduli of the mineral mixture, respectively, ϕs is the stiff porosity, and
(A.3)P=1ν612ν×41+ν+2γ272ν31+4ν+12γ22γg2γ2+14γ2g+γ211+γg2,Q=4γ211γ158γ1+2γ234ν+78ν4γ212νg×81ν+2γ23+4ν+8ν14γ25+2νg+6γ211+νg22γ2+14γ2g+γ211+γg238ν1+2γ254ν+312ν+6γ2ν1g2γ2+2γ+γ21+νg,g=γ1γ23/2arccosγγ1γ2γ<1,γ1γ23/2γ1γ2arccoshγγ>1,
where γ is the spheroidal aspect ratio and ν is the Poisson ratio of the grains, i.e., ν=3K02μ0/6K0+2μ0.
Cracks are introduced into the host material, by neglecting the interactions between cracks and pores. The effective moduli are
(A.4)1KeffMT=1KstiffMT1+161νstiffMT2Γ912νstiffMT,(A.5)1μeffMT=1μstiffMT1+321νstiffMT5νstiffMTΓ452νstiffMT,
where νstiffMT is the Poisson ratio of the rock with only stiff pores, i.e., νstiffMT=3KstiffMT2μstiffMT/6KstiffMT+2μstiffMT, KeffMT and μeffMT are the effective bulk and shear moduli, respectively, when the rock contains cracks and stiff pores, and Γ is the crack density. All the cracks close at high effective pressure, so that there are only stiff pores. The least-square method is used to obtain the optimal aspect ratio of the stiff pores by using Equations (A.1) and (A.2) (see Table 4 in Appendix B).

Step 2. Calculate the cumulative crack density at different pressures. It can be obtained with Equations (A.4) and (A.5) by a least-square method (see Table 5 in Appendix B). When this is known, the moduli can be obtained.

Step 3. Establish the relation between effective pressure and crack density, which can be determined by using the relation between the crack density and effective pressure [58]:
(A.6)Γ=Γiep/p,
where Γi is the initial value when the effective pressure is zero and p^ is a constant.
Step 4. Calculate the crack aspect ratio distribution. When the effective pressure increases, cracks gradually tend to close. The minimum initial aspect ratio of the unclosed cracks under each effective pressure can be obtained as [14]
(A.7)γi=34πΓiΓ1/KΓ1/KeffhpΓdpdΓdΓ,
where KΓ is the effective bulk modulus at different pressures, which can be computed from Equation (A.1).
Substituting Equation (A.8) into (A.9), we obtain
(A.8)γi=34πΓγΓi1/KΓ1/Keffhpp^Γ2dΓ,
and integrating Equation (A.10) from Γ to Γi,
(A.9)γi=4p^1νeffhp2lnΓi/Γ3πKeffhp12νeffhp,
whereνeffhp is the effective Poisson ratio at high pressures, i.e., νeffhp=3Keffhp2μeffhp/6Keffhp+2μeffhp.
Combined with Equations (A.6) and (A.9), the relation between the minimum initial aspect ratio and the effective pressure can be obtained:
(A.10)γi=41νeffhp2pπEeffhp,
where Eeffhp is the effective Young modulus at high pressures, i.e., Eeffhp=3Keffhp12νeffhp. The cumulative crack density decreases when the effective pressure increases. If the effective pressure changes from zero to dp, the corresponding reduction of cumulative crack density is dΓ. When the effective-pressure increment is small enough, the reduction of crack density can be considered to be caused by the closure of the cracks with aspect ratio less than the minimum initial aspect ratio. David and Zimmerman [14] relate crack porosity and crack density as
(A.11)ϕc=4πγ3Γ.

Therefore, the crack properties in rock can be obtained from the wave velocities as a function of the pressure from experimental measurements (e.g., [25]).

B. The Results of Least Square Method and Misfit for Each Rock over the Whole Pressure Range

The data of the modeling results can be accessed by contacting the corresponding author.

The authors declare that they have no conflicts of interest.

This work is supported by the National Natural Science Foundation of China (grant nos. 41974123 and 42174161), the Jiangsu Province Science Fund for Distinguished Young Scholars (grant no. BK20200021), and the Sinopec Key Laboratory of Geophysics.

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