A New Seismic Inversion Scheme Using Fluid Dispersion Attribute for Direct Gas Identification in Tight Sandstone Reservoirs

Abstract

Sufficient gas accumulation is an essential factor that controls the effective development of tight sandstone gas reservoirs that are generally characterized by low porosity and permeability. Seismic methods are important for predicting potential gas areas in tight sandstones. However, the complex relationships between rock physical properties and gas saturation make gas enrichment estimation with seismic methods challenging. Nonetheless, seismic velocity dispersion using a wave-induced fluid flow mechanism can enable gas identification by utilizing the associated dispersion attributes. This paper proposes a method for improved gas identification using a new fluid dispersion attribute obtained by incorporating the decoupled fluid-solid seismic amplitude variation with offset representation into the frequency-dependent inversion scheme. Numerical analyses and synthetic data tests confirmed the enhanced sensitivity of the fluid dispersion attribute to gas saturation compared to the conventionally used compressional wave velocity dispersion attribute. Field data applications further validated the ability of the proposed fluid dispersion attribute to improve gas prediction in tight sandstone reservoirs. The results of the measurements enable rational interpretation of the geological significance of assessments of reservoir properties from gas-producing wellbores. The proposed fluid dispersion attribute is a reliable indicator for gas prediction and represents a useful tool for characterizing tight sandstone reservoirs.

1. Introduction

Effective development of tight sandstone gas reservoirs requires identifying areas with sufficient gas accumulation in tight sandstones, which are generally characterized by low porosity (<10%) and permeability (<1 mD). Many previous studies have attempted to characterize tight sandstone gas reservoirs using seismic methods [1,2,3,4,5]. However, because the elastic properties, and their corresponding seismic attributes, may be insensitive to pore fluids in tight sandstone reservoirs, it is challenging to predict gas saturation using conventional methods based on the elastic assumption.

The inelastic properties of tight sandstones produced by poroelastic behavior associated with fluid flow mechanisms provide an opportunity to use seismic attributes associated with velocity dispersion and attenuation. Xue et al. [6] used seismic attenuation to predict gas-bearing tight sandstones. Experimental measurements [7] and theoretical modeling [8] indicate that the elastic parameters are frequency-dependent owing to fluid mobility in tight sandstones. In addition, Jin et al. [9] characterized tight gas sandstones using compressional wave velocity dispersion. Various poroelastic theories have been proposed to characterize the dissipation mechanisms in heterogeneous reservoir rocks resulting from wave-induced fluid flow [10,11,12,13,14]. Developing new seismic methods for improved gas identification is feasible by considering the inelastic mechanisms of seismic waves in tight sandstones.

The frequency-dependent amplitude variation with offset (FD-AVO) approach is a practical method for fluid identification which utilizes velocity dispersion in hydrocarbon resources. Wilson et al. [15] extended the linearized approximation of the PP-wave reflection coefficient by Smith and Gidlow [16] to the frequency-dependent form, assuming that variation in elastic properties with frequency causes the seismic reflection coefficient to change with the frequency. The authors proposed an FD-AVO inversion scheme to estimate P-wave velocity dispersion. Subsequently, the dispersion attributes associated with P-wave velocity were applied to fluid detection in various reservoirs [3,17,18]. Some researchers have also focused on improving the FD-AVO inversion scheme by incorporating advanced spectral decomposition techniques with high resolution [19] but at the cost of reducing computational efficiency.

In principle, AVO expressions can be extended to a frequency-dependent form for dispersion attribute inversion. As Wang et al. [4] reported, most existing FD-AVO methods estimate dispersion attributes associated with the properties of the overall rock. However, with complex correlations between physical properties and gas saturation due to the low porosity and permeability of tight sandstones, it is challenging to ensure reliable gas identification using the properties of the entire rock. To address this challenge, in elastic seismic inversion, Yin et al. [20] and Zhang et al. [21] decoupled the fluid and solid terms in the AVO expression based on poroelastic theory, proposing the effective fluid bulk modulus as an indicator for hydrocarbon identification. This decoupled fluid-solid AVO expression aims to achieve more accurate fluid identification by directly estimating the effective fluid bulk modulus.

Recognizing the inelastic properties due to fluid flow, we sought to extend the decoupled fluid-solid AVO representation and to develop a new FD-AVO method to obtain the fluid dispersion attribute for improved gas prediction in tight sandstones. First, we briefly reviewed the decoupled fluid-solid AVO representation and constructed a fluid dispersion inversion scheme by extending the decoupled fluid-solid representation to the frequency-dependent form. Next, the proposed method was tested using synthetic data of realistic models from well-log data. The sensitivity of the proposed fluid dispersion attribute to gas saturation was analyzed and compared with that of the conventionally used P-wave velocity dispersion. Finally, the proposed method was applied to seismic field data and the results were validated using the reservoir properties of gas-producing wells.

2. Methodology

2.1. Decoupled Fluid-Solid AVO Linearized Approximation

Russell et al. [22] incorporated poroelastic theory into the AVO representation by Aki and Richards [23] to derive an expression in terms of the Gassmann fluid term f, the shear modulus μ, and the density ρ as follows:

$$R_{PP}\left(\theta \right)=\left[\left(1-\frac{{\gamma }_{dry}^2}{{\gamma }_{sat}^2}\right)\frac{{\sec }^2\theta }{4}\right]\frac{\Delta f}{f}+\left(\frac{{\gamma }_{dry}^2}{4{\gamma }_{sat}^2}{\sec }^2\theta -\frac{2}{{\gamma }_{sat}^2}{\sin }^2\theta \right)\frac{\Delta \mu }{\mu }+\left(\frac{1}{2}-\frac{{\sec }^2\theta }{4}\right)\frac{\Delta \rho }{\rho },$$

(1)

where θ is the incidence angle, γ_dry and γ_sat are the P-to-S-wave velocity ratios for dry and saturated rocks, respectively, and the symbol Δ denotes the difference of the corresponding properties across the interface. For simplicity, the denominators of the reflectivity terms in Equation (1) represent the average values of related properties. This simplified expression is utilized in the relevant equations in this paper.

The Gassmann fluid term f in Equation (1) has the form:

$$f={V_P}^2\rho -{\gamma }_{dry}^2{V_S}^2\rho .$$

(2)

where V_P and V_S denote the P- and S-wave velocities of fluid-saturated rock, respectively.

Han and Batzle [24] decoupled the fluid item f as the product of the gain function G(φ) and the effective fluid bulk modulus K_f:

$$f\approx G\left(\phi \right)K_f,$$

(3)

where $G\left(\phi \right)={\left(1-K_{dry}/K_m\right)}^2/\phi$, K_dry and K_m represent the bulk moduli of the dry rock and solid matrix, respectively, and φ denotes the porosity.

By introducing the critical porosity model [25], the gain function G(φ) in Equation (3) can be related to the porosity φ and the critical porosity φ_c:

$$G\left(\phi \right)=\phi /{{\phi }_c}^2.$$

(4)

Then, after substituting Equations (3) and (4) into Equation (1), the decoupled fluid-solid AVO linearized approximation can be obtained according to Yin et al. [20]:

$$\begin{array}{c}{R}_{PP}\left(\theta \right)=\left[\left(1-\frac{{\gamma }_{dry}^2}{{\gamma }_{sat}^2}\right)\frac{{\sec }^2\theta }{4}\right]\frac{\Delta K_f}{K_f}+\left(\frac{{\gamma }_{dry}^2}{4{\gamma }_{sat}^2}{\sec }^2\theta -\frac{2}{{\gamma }_{sat}^2}{\sin }^2\theta \right)\frac{\Delta f_m}{f_m}\\ \hspace{1em}\hspace{1em} +\left(\frac{1}{2}-\frac{{\sec }^2\theta }{4}\right)\frac{\Delta \rho }{\rho }+\left(\frac{{\sec }^2\theta }{4}-\frac{{\gamma }_{dry}^2}{2{\gamma }_{sat}^2}{\sec }^2\theta +\frac{2}{{\gamma }_{sat}^2}{\sin }^2\theta \right)\frac{\Delta \phi }{\phi },\end{array}$$

(5)

where f_m = μφ represents the dry rock matrix term. The effective fluid bulk moduli (K_f) in Equation (5) can be used to evaluate the fluid properties directly.

2.2. A New Fluid Dispersion Attribute Inversion Scheme Based on Decoupled Fluid-Solid AVO

It is rational to assume that the density (ρ) and the porosity (φ) in Equation (5) do not change with frequency and that K_f and f_m are frequency-dependent. Therefore, the frequency dependence of K_f and f_m causes the reflection coefficients to vary with frequency ω. Consequently, a new dispersion inversion scheme can be established by extending the decoupled fluid-solid AVO approximation in Equation (5) to the frequency-dependent form:

$$\begin{array}{c}{R}_{PP}\left(\theta ,\omega \right)=\left[\left(1-\frac{{\gamma }_{dry}^2}{{\gamma }_{sat}^2}\right)\frac{{\sec }^2\theta }{4}\right]\frac{\Delta K_f}{K_f}\left(\omega \right)+\left(\frac{{\gamma }_{dry}^2}{4{\gamma }_{sat}^2}{\sec }^2\theta -\frac{2}{{\gamma }_{sat}^2}{\sin }^2\theta \right)\frac{\Delta f_m}{f_m}\left(\omega \right)\\ \hspace{1em} +\left(\frac{1}{2}-\frac{{\sec }^2\theta }{4}\right)\frac{\Delta \rho }{\rho }+\left(\frac{{\sec }^2\theta }{4}-\frac{{\gamma }_{dry}^2}{2{\gamma }_{sat}^2}{\sec }^2\theta +\frac{2}{{\gamma }_{sat}^2}{\sin }^2\theta \right)\frac{\Delta \phi }{\phi }.\end{array}$$

(6)

By assuming

$$A\left(\theta \right)=\left(1/4\right)\left(1-{\gamma }_{dry}^2/{\gamma }_{sat}^2\right){\sec }^2\theta ,$$

$$B\left(\theta \right)=\left({\gamma }_{dry}^2/4{\gamma }_{sat}^2\right){\sec }^2\theta -\left(2/{\gamma }_{sat}^2\right){\sin }^2\theta ,$$

$$C\left(\theta \right)=1/2-\left(1/4\right){\sec }^2\theta ,$$

$$D\left(\theta \right)=\left(1/4\right){\sec }^2\theta -\left({\gamma }_{dry}^2/2{\gamma }_{sat}^2\right){\sec }^2\theta -\left(2/{\gamma }_{sat}^2\right){\sin }^2\theta ,$$

Equation (6) was expanded using the Taylor series at a reference frequency ω₀ and a first-order term was retained by ignoring the high-order terms:

$$R_{PP}\left(\theta ,\omega \right)\approx R_{PP}\left(\theta ,{\omega }_0\right)+\left(\omega -{\omega }_0\right)A\left(\theta \right)\frac{\partial }{\partial \omega }\left(\frac{\Delta K_f}{K_f}\right)+\left(\omega -{\omega }_0\right)B\left(\theta \right)\frac{\partial }{\partial \omega }\left(\frac{\Delta f_m}{f_m}\right),$$

(7)

where

$$R_{PP}\left(\theta ,{\omega }_0\right)=A\left(\theta \right)\frac{\Delta K_f}{K_f}\left({\omega }_0\right)+B\left(\theta \right)\frac{\Delta f_m}{f_m}+C\left(\theta \right)\frac{\Delta \rho }{\rho }\left({\omega }_0\right)+D\left(\theta \right)\frac{\Delta \phi }{\phi }.$$

(8)

Here, D_Kf and D_fm were defined as the dispersion attributes associated with K_f and f_m:

$$D_{K_f}=\frac{\partial }{\partial \omega }\left(\frac{\Delta K_f}{K_f}\right),$$

(9)

$$D_{f_m}=\frac{\partial }{\partial \omega }\left(\frac{\Delta f_m}{f_m}\right).$$

(10)

Furthermore, Equation (7) was expressed in a simplified form:

$$\begin{array}{l}\Delta R\left(\theta ,\omega \right)=R_{PP}\left(\theta ,\omega \right)-R_{PP}\left(\theta ,{\omega }_0\right)\\ \hspace{1em}\hspace{1em}\hspace{1em} =\left(\omega -{\omega }_0\right)\left[A\left(\theta \right)D_{K_f}+B\left(\theta \right)D_{f_m}\right],\end{array}$$

(11)

In applications to field data, a wavelet with a spectrum of W(ω) = [W(ω₁), W(ω₂), …, W(ω_n)] was introduced, and Equation (11) was transformed into the time-frequency domain (S), using spectra decomposition methods to generate the following equation:

$$\begin{array}{l}\Delta S\left(\theta ,\omega \right)=S_{PP}\left(\theta ,\omega \right)-S_{PP}\left(\theta ,{\omega }_0\right)\\ \hspace{1em}\hspace{1em}\hspace{1em} =\left(\omega -{\omega }_0\right)W\left(\omega \right)\left[A\left(\theta \right)D_{K_f}+B\left(\theta \right)D_{f_m}\right],\end{array}$$

(12)

Based on Equation (12), the inversion for the dispersion attributes D_Kf and D_fm can be represented in terms of the following matrix for field data application:

$$d=G\left[\begin{array}{l}{D}_{K_f}\\ D_{f_m}\end{array}\right],$$

(13)

where matrices d and G have the following form:

$$d=\left[\begin{array}{c}\begin{array}{c}\Delta S\left({\theta }_1,{\omega }_1\right)\\ ⋮\\ \Delta S\left({\theta }_1,{\omega }_n\right)\end{array}\\ \begin{array}{c}\begin{array}{l}⋮\\ ⋮\end{array}\\ \Delta S\left({\theta }_m,{\omega }_1\right)\\ ⋮\\ \Delta S\left({\theta }_m,{\omega }_n\right)\end{array}\end{array}\right],$$

(14)

$$G=\left[\begin{array}{cc}W\left({\omega }_1\right)\left({\omega }_1-{\omega }_0\right)A\left({\theta }_1\right)& W\left({\omega }_1\right)\left({\omega }_1-{\omega }_0\right)B\left({\theta }_1\right)\\ ⋮& ⋮\\ W\left({\omega }_n\right)\left({\omega }_n-{\omega }_0\right)A\left({\theta }_1\right)& W\left({\omega }_n\right)\left({\omega }_n-{\omega }_0\right)B\left({\theta }_1\right)\\ ⋮& ⋮\\ ⋮& ⋮\\ W\left({\omega }_1\right)\left({\omega }_1-{\omega }_0\right)A\left({\theta }_m\right)& W\left({\omega }_1\right)\left({\omega }_1-{\omega }_0\right)B\left({\theta }_m\right)\\ ⋮& ⋮\\ W\left({\omega }_n\right)\left({\omega }_n-{\omega }_0\right)A\left({\theta }_m\right)& W\left({\omega }_n\right)\left({\omega }_n-{\omega }_0\right)B\left({\theta }_m\right)\end{array}\right].$$

(15)

For the spectral decomposition of the seismic AVO data with n incidence angles and m frequencies at each sampling time, the Levenberg–Marquardt algorithm was used to calculate the dispersion attributes D_Kf and D_fm:

$$\left[\begin{array}{l}{D}_{K_f}\\ D_{f_m}\end{array}\right]={\left(G^{\mathrm{T}}G+kI\right)}^{-1}{G}^{\mathrm{T}}d,$$

(16)

where I denotes the identity matrix and k denotes the damping factor.

2.3. Inverse Spectral Decomposition Technique

Inverse spectral decomposition (ISD) is an effective technique for obtaining time-frequency spectra with high temporal and frequency resolutions, which can facilitate the dispersion attribute inversion presented in Equations (13) to (16). As suggested by Han et al. [26], if the frequency of a signal s(t) varies over time, it can be expressed by extending the traditional convolution of the wavelet w(t) and reflection coefficient r(t) to the corresponding frequency-dependent form:

$$s\left(t\right)={\displaystyle \sum _{n=1}^N\left[w_n\left(t\right)\ast r_n\left(t\right)\right]},$$

(17)

where * indicates the convolution operation, N denotes the number of frequencies, n is the dominant frequency of the n-th Ricker wavelet w_n, and r_n represents the reflection coefficient for the n-th Ricker wavelet.

Equation (17) can be expressed in a matrix form as follows:

$$s={\displaystyle \sum _{n=1}^N\left(W_n{R}_n\right)},$$

(18)

where W_n = [W₁, W₂, …, W_N] denotes the wavelet matrix consisting of the frequency-dependent Ricker wavelets, and R_n = [R₁, R₂, …, R_N]^T is the corresponding reflection coefficient varying with frequency.

Equation (18) can be represented in a simplified form as:

$$s=WR,$$

(19)

by neglecting the subscript “n” for simplicity.

The method for calculating the time-frequency spectrum R in Equation (19) is called ISD. Han et al. [26] showed that to reduce the ambiguity of the inverse problem, R can be sparsely constrained by the L1 norm with regularization represented as follows:

$$J=\frac{1}{2}{\Vert WR-s\Vert }_2^2+\lambda {\Vert R\Vert }_1, \lambda >0,$$

(20)

where J denotes the cost function, ||R||₁ indicates the constraint term, and λ is the parameter that controls sparseness. Equation (20) can be solved using a fast iterative soft-thresholding algorithm.

In Figure 1, spectral decompositions obtained using ISD were compared with those using the continuous wavelet transform (CWT) based on synthetic data. A synthetic signal containing three Ricker wavelets, with different dominant frequencies and center times, is displayed in Figure 1a. A Ricker wavelet with 20 Hz dominant frequency was centered at 100 ms, and the other two Ricker wavelets with 40 Hz dominant frequency were superimposed and located around 300 ms. The time-frequency spectra of the signal were computed using the CWT and ISD, as shown in Figure 1b,c. The results indicated that the ISD method provided superior temporal and frequency resolutions compared to the CWT.

2.4. Rock Physics Model for Gas-Bearing Tight Sandstones

Previous studies suggest that tight sandstones exhibit poroelastic behavior with dissipation mechanisms related to wave-induced fluid flow [7,8]. Reasonable interpretations of inverted dispersion attributes require appropriate rock physics modeling for tight sandstones. We developed a rock physics modeling workflow to describe dissipation mechanisms in tight sandstones, as shown in Figure 2.

First, the Hashin–Shtrikman bounds (HSB) were used to obtain the elastic modulus of the rock matrix, which mainly consisted of quartz and clay. Then, the dual-porosity model, dividing the pore space into soft and stiff pores, was applied to characterize the micropore structure of the tight sandstone [27]. Based on the self-consistent approximation (SCA), the two types of pores were added to the solid matrix to obtain the elastic modulus of the dry rock frame. Finally, fluids were filled into the pore spaces using patchy saturation theory that describes the dissipation mechanism caused by mesoscale wave-induced flow within the seismic frequency band. A detailed description of patchy saturation theory is provided in Appendix A.

3. Theoretical Model Analyses

3.1. Rock Physics Modeling for Gas-Bearing Tight Sandstone

The properties of the minerals and fluids used in rock physics’ modeling for tight sandstones are listed in

Table 1. Properties for rock physics modeling of tight sandstone [28,29].

Properties	Clay	Quartz	Water	Gas
K (GPa)	21.00	36.60	2.25	0.012
μ (GPa)	7.00	45.00	0.00	0.000
ρ (g/cm3)	2.60	2.65	1.04	0.078
η (cP)	-	-	3.00	0.015

. Based on the work of Guo et al. [27] for the study area, the volumetric fractions of quartz and clay were set to 0.87 and 0.13, respectively. The porosity was assumed to have a constant value of 0.1. The soft pores or microcracks were assumed to have an aspect ratio of 0.01, whereas the stiff pores were assumed to have an aspect ratio of 1. The volumetric proportions of the stiff and soft pores were set to 9:1.

Figure 3a shows the frequency-dependent P-wave velocity (V_P) computed with the proposed model (Figure 2). As the gas saturation (S_g) varies from 0.1 to 0.7, the magnitude of V_P decreases, and the corresponding P-wave velocity dispersion varies within the seismic band below 100 Hz. In addition, as indicated by Equation (2), the velocity dispersion of the compressional wave can lead to the frequency-dependence of the Gassmann fluid term f. As indicated by Equation (3), f can be represented as the product of G(φ) and K_f. Because G(φ) is only related to the porosity and the properties of the solid matrix and the dry frame, it is rational to assume that G(φ) does not change with frequency. The frequency dependence of f is attributed to K_f varying with frequency. Thus, we extended Equation (3) to its frequency-dependent from:

$$f\left(\omega \right)\approx G\left(\phi \right)K_f\left(\omega \right).$$

(21)

Because K_f is only associated with the fluid properties and is not affected by the rock matrix, it is predicted to exhibit greater sensitivity to S_g than to V_P of rock. The frequency-dependent K_f value was obtained using Equation (21); the frequency-dependent f value was computed with Equation (2) and G(φ) was determined with our proposed rock physics model (Figure 2). Then, we obtained the frequency-dependent V_P and K_f values for various S_g values (Figure 3a,b). The results indicated that the variation rate in V_P at 40 Hz was ~4% for S_g values ranging from 0.1 to 0.7, while K_f varied by ~70% for the same S_g range.

The reflectivity terms of ΔV_P/V_P and ΔK_f/K_f versus frequency for different S_g values were computed for an interface model separating mudstone and tight sandstone (Figure 4a,b). In the modeling, the V_P value of the overlaid mudstone was assumed to be 3900 m/s according to logging data; the V_P values of the tight sandstone for various S_g models are given in Figure 3a. The K_f values of tight sandstone for different S_g values are given in Figure 3b. The overlaid mudstone was assumed to be 100% water-saturated with the corresponding K_f value equal to that of water. The results presented in Figure 4 indicate that ΔK_f/K_f was an order of magnitude higher than ΔV_P/V_P in terms of the frequency-dependent variation for different S_g models. Therefore, the enhanced sensitivity can be predicted from the fluid dispersion attribute D_Kf, defined by the relationship of ΔK_f/K_f to frequency (Equation (9)).

In this study, we assumed that the fluid dispersion attribute was associated with the wave-induced fluid flow described by patchy saturation theory [30]. However, depending on the characteristics of tight sandstones in specific areas, complex poroelastic behaviors in tight sandstones may be described by different theories, including the multi-scale fracture model [10,11] and the double-porosity model [12,13,14].

3.2. Synthetic Data Test Using a Realistic Model

To accurately evaluate our proposed D_Kf for gas prediction, we utilized synthetic data for a realistic model based on logging data from the study area. The main challenge faced in applying the realistic model is rock physical and seismic modeling that simultaneously considers the inelasticity and heterogeneity of the target reservoir. In this study, we addressed this challenge by integrating the model for the inelastic pathy-saturation mechanism (Figure 2) and the dual-porosity model applied to well-log data [27]. Then, synthetic data for the inelastic layered model were computed with a seismic modeling tool based on the propagator matrix method (PMM) [31,32].

Logging data from a borehole in the study area in the Ordos basin suggested that the gas-bearing tight sandstone was located at a depth of ~1335–1664 m (Figure 5). The porosity (φ), permeability (κ), gas saturation (S_g), and volumetric mineral fraction curves, providing realistic model properties utilized for modeling, are shown.

Figure 6a shows the S_g curves of the target tight sandstone within the depth range indicated in Figure 5. For the tight sandstone within this depth range, we considered three models where S_g had change rates of −50%, 0%, and +50%, corresponding to the original S_g curve in Figure 5. Figure 6b–d illustrates the frequency-dependent V_P values of the three S_g models in Figure 6a, which were calculated based on the abovementioned method.

Using the PMM as the seismic modeling tool and the seismic wavelets (Figure 7) extracted from real seismic data, we calculated the synthetic AVO data from the three realistic tight sandstone models with different S_g values (Figure 6). The stratigraphic structure surrounding the target sandstones for seismic modeling can be extracted from Figure 5. Figure 8a–c show the computed AVO responses of the three S_g models in Figure 6b–d. Variations in seismic reflections are apparent for the three S_g models. In Figure 8, actual reflections from the bottom of the tight sandstone correspond to peak amplitudes at ~55 ms. The time window of the selected seismic reflections has a slightly wider range than the actual reflections from the tight sandstone models to maintain consistency with the interpreted horizon in the field data illustrated in the next section.

Using the synthetic data, we compared the conventionally used D_P and our proposed D_Kf (Figure 8d,e). The results indicated that D_Kf exhibited more distinct anomalies in the target tight sandstone with less ambiguity, suggesting the superiority of D_Kf over D_P.

In addition, we computed D_Kf and D_P for the models with continuously varied S_g, with a change rate from −50% to 50% with respect to the original S_g curve measured in the borehole (Figure 5). Then, the root-mean-square values of D_Kf and D_P for the target layer in Figure 8 were estimated and normalized to a range of 0–1 for comparison (Figure 9). As shown in Figure 9, the magnitude of D_Kf was much larger than that of D_P, where D_P was comparable to D_Kf after approximately 100× magnification. The results indicate the enhanced sensitivity of D_Kf to variation in S_g.

Although D_Kf and 100 × D_P exhibited similar responses to S_g with change rates between −50% and −20%, D_Kf exhibited more apparent sensitivity to S_g than D_P for a relatively higher S_g with a change rate from −20% to +50%. This result confirms the improved ability of the proposed D_Kf compared to D_P for discriminating potential areas with high gas enrichment in tight sandstones.

4. Field Data Applications

4.1. Dataset

Our proposed method was applied to real seismic data for the tight sandstone gas reservoirs from the Ordos Basin, China. Figure 10 shows the two-way travel time of reflections from the target reservoir which exhibit a relatively smooth tectonic relief. The locations of four wells are indicated; wells A and B are gas-producing and wells C and D are dry. The seismic profile across the four wells is shown in Figure 11. The seismic reflections exhibit continuous events around the gas-producing wells A and B in the target layer. However, reliably locating gas zones using seismic amplitudes can be challenging.

Figure 12 and Figure 13 show the decomposed spectra computed with the CWT and ISD methods at frequencies of 10, 20, 30, 40, 50, and 60 Hz, respectively, for the post-stack seismic profile (Figure 11). The spectral decompositions obtained with the ISD method exhibited higher temporal resolutions than those using the CWT method.

4.2. Calibration of D_Kf Using Actual Seismic Data near the Gas-Producing Well

The proposed D_Kf was calibrated for improved fluid identification in tight gas sandstones using actual AVO data collected adjacent to the gas-producing well A. As shown in Figure 14a, the AVO data had incidence angles of 5°, 15°, and 25°. The decomposed spectrum of the corresponding stacked trace was calculated with the ISD method (Figure 14b).

As shown in Figure 14c, strong D_P anomalies, unrelated to the target gas layer, can cause ambiguity in gas identification in tight sandstone. Notably, the anomaly response of D_Kf (Figure 14d) was more apparent than that of D_P (Figure 14c) in the gas-bearing tight sandstone. Moreover, the anomaly responses unrelated to the target gas layer were suppressed by D_Kf compared to D_P. The results support the proposed D_Kf for improved gas identification in tight sandstones. The results of the synthetic data tests for the realistic models (Figure 8d,e) coincided with the field data results (Figure 14c,d). The good correlations obtained justify the rock physical and seismic modeling method used for the realistic model test, and confirm the reliability of the results presented in Figure 9.

4.3. Applications to Three-Dimensional (3D) Seismic Data

The proposed inversion method for the fluid dispersion attribute was applied to 3D seismic data acquired in the study area. The sections of the inverted dispersion attributes D_P and D_Kf are shown in Figure 15a,b. High D_P and D_Kf value anomalies were observed in the tight sandstones in the gas-producing wells A and B, indicating potential gas zones in the target reservoirs. In contrast, no D_P or D_Kf responses were observed in the target layer in wells C and D, which exhibited no gas yield. Furthermore, the results for D_Kf indicated that the gas zone around well A exhibited a more distinct anomaly and covered a broader range than the zone adjacent to well B (Figure 15b). In contrast, D_P predicted a more distinct gas zone in well B than in well A (Figure 15a).

Figure 16 shows horizontal slices of D_P and D_Kf for the target tight sandstone, computed using 3D seismic data from the study area. The attribute D_P (Figure 16a) indicated that the gas zone around well B covered a relatively larger area than the zone around well A. In contrast, the dispersion attribute D_Kf (Figure 16b) indicated that the gas zone around well A covered a much larger area than the zone adjacent to well B.

The results in Figure 15 and Figure 16 can be calibrated with the properties of the gas-bearing tight sandstone reservoirs measured in wells A and B (

Table 2. Tight gas sandstone reservoir properties measured in two gas-producing wells.

Well	h (m)	φ	κ (mD)	Sg
Well A	11.90	0.136	2.70	0.72
Well B	2.80	0.131	1.77	0.51

). As shown in Table 2, the porosities φ of the two wells were comparable. The permeability κ was relatively higher in well A than in well B but of the same order of magnitude in the two wells. However, the gas-bearing tight sandstone exhibited a higher S_g in well A than in well B. Therefore, it is reasonable to expect that D_Kf would exhibit more distinct anomaly responses compared to D_P for the gas zone around well A with a higher S_g, as shown in Figure 15b and Figure 16b.

As shown in Table 2, the gas-bearing tight sandstone was thicker in well A than in well B. It can be inferred that the thicker gas-bearing sandstone around well A covered a broader area than that around well B, as shown in Figure 16b, confirming the reliability of the results presented in Figure 15b and Figure 16b. Therefore, the proposed attribute D_Kf represents a preferred indicator for reliable gas identification. The results in Figure 16b provide reasonable estimates of the distribution of potential gas zones for future prospecting.

5. Discussion

5.1. Advantages of the Proposed Fluid Dispersion Attribute

In principle, various AVO representations can be extended to obtain corresponding seismic dispersion attributes based on the frequency-dependent inversion scheme [4]. However, most existing approaches generally estimate dispersion attributes associated with the overall properties of the rocks, such as the commonly used P-wave velocity dispersion attribute (D_P) [3,9,15,17,19]. In comparison, this study proposed the fluid dispersion attribute (D_Kf) for direct gas identification, which was obtained by incorporating the decoupled fluid-solid AVO equation into the frequency-dependent inversion scheme. Accordingly, the proposed D_Kf can be anticipated to exhibit enhanced fluid sensitivity for gas prediction in tight sandstones.

Numerical analyses (Figure 3 and Figure 4) and synthetic data tests (Figure 8 and Figure 9) validated the enhanced sensitivity of our proposed D_Kf to gas saturation (S_g) in tight sandstones. Specifically, the model test results suggested that D_Kf was of a much higher order of magnitude than D_P. Moreover, the anomaly response of D_Kf showed an enhanced sensitivity over that of D_P for higher S_g (Figure 9), indicating that D_Kf is a preferable indicator for gas identification in tight gas sandstones.

Furthermore, field data applications indicated that the anomaly response of D_Kf (Figure 14d) was more apparent than that of D_P (Figure 14c) in gas prediction, with the response being unrelated to the target gas layer well suppressed by D_Kf. The results suggest that D_Kf can discriminate gas zones more accurately than D_P. Moreover, the distribution of potential gas areas estimated by D_Kf had more rational geological significance with respect to reservoir properties measured in the gas-producing wellbores (Figure 15b and Figure 16b). The results support the use of the proposed D_Kf for improved gas identification in tight sandstones.

5.2. Rock Physics Modeling of Gas-Bearing Tight Sandstones

Extending Equation (5) to the corresponding frequency-dependent form, as given in Equation (6), does not entail any particular inelastic mechanisms of the seismic wave propagating in tight gas sandstones. Therefore, appropriate interpretations of inverted dispersion attributes require rational rock physics modeling for tight sandstones. To quantify the velocity dispersion associated with varied S_g for interpreting the obtained D_Kf, we attributed the inelastic properties of tight gas sandstones to the wave-induced fluid flow described by patchy saturation theory [8,30]; a flowchart of the modeling method is illustrated in Figure 2. However, complex poroelastic behavior in tight sandstones may be associated with various inelastic mechanisms that are described by different theories, including the multi-scale fracture model [9,11] and the double-porosity model [12,13,14]. In addition, the shear-modulus-related dispersion attribute obtained in Equation (10) was not discussed in the present study due to a lack of knowledge of inelastic shear wave behaviors using current rock physics methods. Future laboratory measurements and rock physics modeling will enable comprehensive interpretations to be provided for the obtained dispersion attributes.

To accurately simulate the seismic signatures of gas-bearing tight sandstones, the proposed model (Figure 2) was extended, in realistic modeling tests, to the case of layered structures by incorporating our previous modeling method [27]. The modeling results for frequency-dependent velocities with layered structures (Figure 6) provided more accurate models of the realistic tight gas sandstones for testing the different frequency-dependent inversion methods.

5.3. Seismic Modeling for Inelastic Tight Gas Sandstone

An advanced modeling method is essential to accurately model the seismic signatures of inelastic layered models. The PMM method [31,32,33] provides a powerful tool for efficiently simulating seismic signatures for the established models (Figure 6). As illustrated in Figure 8 and Figure 14, the modeled seismic responses for the realistic models were consistent with actual seismic data adjacent to the corresponding borehole, leading to the characteristics of D_Kf and D_P for realistic models similar to those obtained from real seismic data. Therefore, the results obtained confirm the applicability of the proposed rock physics modeling and seismic simulation methods for generating reliable synthetic data. Furthermore, the methods can be generalized and applied to other hydrocarbon reservoirs with specific properties and geological structures.

5.4. Spectral Decomposition Techniques Used in the Dispersion Attribute Inversion

Previous studies have mainly focused on improving the FD-AVO inversion by incorporating advanced spectral decomposition techniques with high temporal and frequency resolutions [18,19]. This paper proposed that the FD-AVO applied for gas prediction can be significantly improved by introducing an effective dispersion attribute (i.e., D_Kf). Specifically, the results of D_Kf resulted in more reliable estimates of potential gas areas than D_P when the two attributes were obtained by conducting frequency-dependent inversion using the same spectral decomposition method (e.g., ISD in this paper). In future research, other advanced time-frequency analysis techniques may be employed to obtain spectra with higher temporal and frequency resolutions. Wavelet analysis methods can be used to capture the time-frequency features of seismic signatures [34,35]. However, the trade-off between computational efficiency, robustness, and temporal and frequency resolutions should be properly evaluated.

At specific frequencies, the decomposed spectra (Figure 12 and Figure 13) may exhibit similar features to the dispersion attributes (Figure 15a,b). However, the physical implications of the time-frequency analysis results have not been explicitly delineated, making the straightforward application of time-frequency analysis techniques for gas prediction challenging. In comparison, the dispersion attributes obtained with FD-AVO are physics-based and can significantly improve gas identification by employing a more efficient dispersion attribute (i.e., D_Kf).

6. Conclusions

Effective seismic methods for reliable fluid identification are essential for predicting favorable gas-bearing areas in tight sandstone reservoirs. Seismic velocity dispersion due to poroelastic behaviors of tight sandstone enables gas prediction by employing the associated dispersion attributes. Accordingly, this paper developed a new FD-AVO method to estimate D_Kf for improved gas prediction. The proposed method was tested using synthetic data from realistic models and applied to field data with predicted results validated by reservoir properties measured in gas-producing wells. The main conclusions are as follows:

• The fluid dispersion attribute (D_Kf) obtained by incorporating the decoupled fluid-solid AVO equation into the frequency-dependent inversion scheme provides a preferable gas indicator that improves gas prediction over conventional D_P in tight sandstones.
• Numerical analyses and synthetic data tests validated the enhanced sensitivity of D_Kf to S_g. The response of D_Kf was of a much higher order of magnitude than D_P and showed more distinct anomalies than D_P for higher S_g, confirming that D_Kf represents a superior indicator for gas prediction identification in tight sandstones.
• Field data applications indicated that the results of D_Kf led to more reliable estimates of potential gas-bearing areas. The distribution of gas zones estimated by D_Kf were more meaningful geologically as revealed by the reservoir properties from gas-producing boreholes, supporting the use of D_Kf for improved gas identification in tight sandstones.
• The rock physics modeling and seismic simulation methods proposed offer powerful tools for testing models of the fluid dispersion attribute inversion. The synthetic data computed using the established realistic models provide reliable representations of the seismic signatures of actual tight sandstone reservoirs.

Future research should include exploration of other potential inelastic mechanisms of seismic waves propagating in tight sandstones based on rock physics modeling and laboratory measurements. The shear wave-related dispersion attribute also deserves further investigation to enable comprehensive interpretation of the obtained dispersion attributes. The proposed rock physics modeling and seismic simulation methods can be generalized to other reservoirs with particular properties and geological structures.