Rock Physical Modeling and Seismic Dispersion Attribute Inversion for the Characterization of a Tight Gas Sandstone Reservoir

Gas identification using seismic data is challenging for tight gas reservoirs with low porosity and permeability due to the complicated poroelastic behaviors of tight sandstone. In this study, the Chapman theory was used to simulate the dispersion and attenuation caused by the squirt flow of fluids in the complex pore spaces, which are assumed to consist of high aspect-ratio pores (stiff pores) and low aspect-ratio microcracks (soft pores). The rock physics modeling revealed that as the gas saturation varies, P-wave velocity dispersion and attenuation occurs at seismic frequencies, and it tends to move to high frequencies as the gas saturation increases. The velocity dispersion of the tight gas sandstone causes a frequency-dependent contrast in the P-wave impedance between the tight sandstone and the overlying mudstone, which consequently leads to frequency-dependent incidence reflection coefficients across the interface. In the synthetic seismic AVO modeling conducted by integrating the rock physics model and the propagator matrix method, the variations in the amplitudes and phases of the PP reflections can be observed for various gas saturations. The tests of the frequency-dependent AVO inversion of these synthetic data revealed that the magnitude of the inverted P-wave dispersion attribute can be used to indicate gas saturation in tight sandstone reservoirs. The applications of the frequency-dependent AVO inversion to the field pre-stacked seismic data revealed that the obtained P-wave dispersion attribute is positively correlated with the gas production from the pay zone at the well locations. Thus, the methods of the rock physics modeling and the frequency-dependent AVO inversion conducted in this study have good potential for the evaluation of the gas saturation in tight gas sandstone reservoirs.


INTRODUCTION
Tight gas sandstone is generally characterized by a strong heterogeneity, low effective porosity, and extremely low permeability. It is challenging to predict the complicated poroelastic behaviors of tight gas sandstone due to the presence of a complex pore and crack system. Correspondingly, the seismic responses from tight gas sandstone are complex because of the dispersion and attenuation caused by the squirt flow of fluids in the complex pore spaces when the seismic waves pass through the tight sandstone. Thus, in order to evaluate the fluid properties of tight sandstone, it is essential to construct valid rock physics models for the investigation of seismic dispersion and attenuation and to develop effective techniques for estimating gas saturation using such poroelastic properties.
For rock physics modeling for a poroelastic medium, Biot (1956a,b) established a macroscopic pore elasticity model in order to consider global fluid-flow in a fluid saturated medium. White (1975) proposed a model that considers the dispersion and attenuation that occurs at the seismic frequencies associated with the patchy saturation of two fluid phases. Dvorkin and Nur (1993) developed the Biot/Squirt (BISQ) model by combining the global Biot flow and the local squirt-flow mechanism. Chapman et al. (2002) and Chapman (2003Chapman ( , 2009 proposed the theory of multiscale fractures for an effective medium containing microscopic pores and cracks and with mesoscopic fractures. Ba et al. (2012) investigated the Biot-Rayleigh theory for modeling tight gas sandstone using a double-pore model. Ba et al. (2016) researched compressional wave dispersion due to rock matrix stiffening by clay squirt flow. Ba et al. (2017) investigated the wave propagating in an anelastic rock with patchy-saturation and fabric heterogeneity.
For the investigation of the seismic responses associated with fluid flow, the frequency-dependent amplitude variation with offset (AVO) method was proposed. Chapman et al. (2005Chapman et al. ( , 2006 analyzed the frequency dependency of the AVO signatures associated with the dispersion and attenuation caused by the presence of fluids. Wilson et al. (2009) and Wilson (2010) extended the conventional AVO approximation to frequency dependency and derived the P-wave dispersion attribute as an indicator for reservoir fluid characterization, in which the spectral decomposition technique was used in the implement of the dispersion attribute inversion. Wu (2010) improved the accuracy of the time-frequency analysis in the inversion in the application of field data. Zhang et al. (2014) extended the AVO Shuey approximation and introduced the concepts of a frequency-dependent intercept and gradient into the AVO inversion. Li (2013) investigated numerical and physical modeling as well as field data for the frequency-dependent AVO in order to enable the quantitative estimation of gas saturation. Wu et al. (2014) also illustrated another field case study of the frequency-dependent AVO by improving the procedure of inversion. Chen et al. (2015) quantitatively calculated the gas saturation based on the dispersion inversion using stacked seismic data. Wang et al. (2019) investigated the problem of reflection and transmission of plane elastic waves at an interface between two double-porosity media.
In this study, we developed a rock physics model for a tight gas sandstone reservoir based on the Chapman theory that considers the squirt flow of fluids in complex pore spaces. Then we used our model to investigate the influence of gas saturation on P-wave velocity dispersion and attenuation. The frequency-dependent AVO modeling and inversion were tested and analyzed for various gas saturations using synthetic data. Finally, the frequency-dependent AVO inversion was applied to field data for a tight gas sandstone reservoir, and the P-wave dispersion attributes were inverted and calibrated in order to determine the gas saturation in the tight sandstone reservoir.

Rock Physics Model of Tight Sandstone
For tight gas sandstone reservoirs, their complex pore structures result in complicated elastic properties and seismic responses associated with fluids. Smith et al. (2009) and Ruiz and Cheng (2010) investigated a rock model including soft and stiff pores, which divided the complex pores into stiff pores and soft pores. They studied the influence of the distributions of stiff and soft pores on the elastic properties of tight sandstone. Their results indicated that the dual-porosity model can well predict the elastic behaviors of tight sandstone. However, this model could not describe inelasticity associated with fluids. Chapman (2003) proposed the theory of multiscale pores and fractures and studied the dispersion and attenuation related to the squirt flow of fluids in complex pore spaces in gas sandstone reservoirs. Thus, we used the concept of stiff and soft pores in Chapman's theory to describe pore structures of tight gas sandstone reservoirs. A schematic diagram of the complex pore spaces consisting of stiff and soft pores is shown in Figure 1.
In Chapman's theory (Chapman, 2003), the expression for the anisotropic stiffness matrix of an equivalent medium is where C 0 ijkl is the stiffness of the isotropic background; C 1 ijkl , C 2 ijkl , and C 3 ijkl are the disturbances caused by the pores, microcracks, and fractures, respectively; ϕ p is the porosity; ε c is the crack density; and ε f is the fracture density.
According to geological understanding of the tight sandstone reservoir in this study area, mesoscale fractures are not developed, and the pore space is dominated by micro pores and cracks. Thus, in Eq. 1, we maintain the terms of the micro pores and cracks, while we drop the term of the mesoscale fractures. So, Eq. 1 becomes According to Chapman et al. (2002), the relaxation time τ plays an important role in determining the frequency bandwidth of the seismic dispersion and attenuation: where µ is the shear modulus; and v is the Poisson ratio of the solid matrix. a is the crack radius; r is the aspect ratio; η is the fluid viscosity; ς is the grain size; and k is the permeability. The frequency-dependent P-wave velocity V p (ω) and the attenuation factor 1/Q(ω) can be calculated from C 3333 (ω) in Eq. 2 and the density ρ s of the fluid-saturated sandstone. The equations are as follows: where Re and Im indicate the real part and imaginary part of a complex modulus, respectively.

Physical Parameters of a Mixed Fluid
The viscosity of the mixed fluid in Eq. 1 has an impact on the frequency band of the dispersion and attenuation. For a mixture of water and gas, the viscosity has the following form (Davide and Carcione, 2003): where η g is the viscosity of the gas; η w is the viscosity of water; and S g is the gas saturation. Under different gas saturations, the density of the fluid mixture and that of the fluid saturated sandstone are as follows: where ρ f is the density of the fluid mixture; ρ g is the density of the gas; ρ w is the density of water; ρ s is the density of the saturated rock; ρ m is the density of the solid matrix; and ϕ is the porosity.
The bulk module K f and the velocity V f of the fluid mixture are calculated using Wood's equations (Wood, 1995): In these equations, K g and K w are the bulk modulus of gas and water, respectively.

The Method for the Modeling of the Frequency-Dependent AVO
Based on the rock physics model, the propagator matrix theory was used to calculate the seismic responses of the tight sandstone gas reservoir. The reflection and transmission coefficient vectors in the case of the P-wave incidence are solved using the following equation (Carcione, 2001): where the matrices A 1 and A 2 are the propagation matrices related to the elastic moduli of the upper and lower media, respectively. B α =T(0)T −1 (hα)(α=1,...,N) is the propagation matrix of the middle layer, which has an N-layer structure. h α is the thickness of each individual layer. i p is the P-wave incident vector related to the elastic properties of the incidence medium. The reflection coefficient R pp at each frequency can be calculated using the propagator matrix theory. We use R f to denote R pp in the following. By multiplying the frequency reflection coefficient R f by the spectrum of the seismic wavelet W f , we can calculate the amplitude spectrum U f of the reflected wave.
After the inverse Fourier transform of U f , the reflection waveform U t in the time domain can be calculated as where f is the angular frequency; i is the unit of the imaginary number; and t is time.

The Frequency-Dependent AVO Theory
Based on the AVO approximate formula given by Shuey (1985); Chapman et al. (2006) deduced the frequency-dependent AVO formula shown as where D P represents the derivatives of the seismic wave velocities with frequency; and D G is the frequency-dependent AVO gradient.

Dispersion and Attenuation in Tight Sandstone
Based on the rock physics model discussed above, Figure 2 shows the calculated velocity dispersion and attenuation of a tight gas sandstone reservoir for different gas saturations. The properties used in the rock physics model are presented in Table 1. As can be seen from Figure 2A, the P-wave velocity Vp in the seismic frequency bandwidth decreases as the gas saturation Sg increases from 0 to 0.6 to 1. Moreover, the range of the frequency at which the velocity dispersion occurs moves to a higher frequency. Accordingly, as illustrated in Figure 2B, the attenuation peak also moves to a higher frequency as Sg increases. According to Eqs. 3, 6, Sg affects the viscosity coefficient of the fluid mixture, and thus, it changes the relaxation time parameter τ, which controls the frequency range of the velocity dispersion and attenuation. Therefore, variations in the gas saturation have a significant impact on the dispersion and attenuation of the tight gas sandstone reservoirs.

Seismic Modeling and Dispersion Attribute Inversion in the Theoretical Models
In the theoretical model shown in Figure 1, the velocity dispersion and attenuation of the tight gas sandstone are calculated using results shown in Figure 2. Given a density of 2,460 kg/m 3 for tight sandstone, the corresponding P-wave impedance can be calculated for various gas saturations as shown in Figure 3A. Also, we assume that the surrounding mudstone is elastic and has a P-wave velocity of 4,000 m/s, a shear wave velocity of 2,350 m/s, and a density of 2,650 kg/m 3 . Thus, the FIGURE 7 | (A) Stacked seismic responses for three gas saturations of 1, 0.6, and 0; (B) spectrum of the PP reflection for a gas saturation of 1; (C) spectrum of the PP reflection for a gas saturation of 0.6; (D) spectrum of the PP reflection for a gas saturation of 0; (E) inverted Dp for the three gas saturation cases. frequency-independent P-wave impedance of the mudstone is shown by the gray line in Figure 3A.
The two intersecting points in Figure 3A show the frequencies where the P-wave impedance Ip of the shale equals that of the tight sandstone with a gas saturation of 0.6 and 1, respectively. Accordingly, Figure 3B shows that normal-incidence reflection coefficients are 0 at the two frequencies, and will change from negative to positive for increasing frequencies across the two frequencies.
In the following, the frequency-dependent seismic responses of the tight sandstone gas reservoir are calculated using the propagator matrix method. We assume that the thickness of the tight sandstone layer is 5 m. Figure 4 shows the absolute values of the complex frequency-dependent PP wave reflection coefficients for gas saturations of 1, 0.6, and 0, with incidence angles ranging from 0 • to 40 • . For an incident Ricker wavelet of 45 Hz, the amplitude spectra of the reflected PP waves can be obtained by multiplying the spectrum of the incident Ricker wavelet by the frequency-dependent reflection coefficients. The results are shown in Figure 5. Furthermore, the reflected waveforms can be calculated by conducting an inverse Fourier transform of the reflected spectra. The results are shown in Figure 6. For the various gas saturations, there were significant differences in the PP wave frequency-dependent reflection coefficients, the spectrum, and the waveforms.
As can be seen from the AVO waveforms in Figure 6, the phase of the waveform for complete gas saturation is opposite to that for complete water saturation, and the variations in the reflection amplitude with incidence angle are different. When the gas saturation is 0.6, the polarity reversal can be observed as the incidence angle increases. Figure 7A shows the stacked seismic responses for various gas saturations obtained from AVO gathers in Figure 6. Figures 7B-D show the frequency spectrum of the AVO reflections corresponding to gas saturations of 1, 0.6, and 0, respectively. Figure 7E shows the P-wave dispersion attribute D P calculated using the frequency-dependent AVO inversion method. As can be seen, the case of high gas saturation corresponds to the value of the high D P attribute. Due to the change in reflection phase FIGURE 11 | Profile of the dispersion attribute Dp calculated using the frequency-dependent AVO inversion method. Frontiers in Earth Science | www.frontiersin.org for the case of the gas saturation of 0.6, there is a time shift for the inverted D P attribute compared to the cases of complete gas saturation and complete water saturation. Thus, the calculation results using synthetic data of the theoretical model verify the feasibility of identifying gas-bearing tight sandstone reservoirs using the P-wave velocity dispersion attribute.

REAL DATA APPLICAIONS
Logging and Seismic Data Figure 8 shows the logging data for the tight sandstone gas reservoir in the study area. The target layer is tight sandstone with a thickness of about 5.5 m below a depth of about 1,325 m. The target layer has a lower gamma ray value and relative higher Pand S-wave velocities than the surrounding mudstone. It has a density about 2.5 g/cm 3 and a porosity of about 0.12. The quartz content of the tight sandstone is greater than 85%. Figure 9 shows the post-stacked seismic profile across the well. The yellow curve shows the gamma ray values from the well, and the target layer can be recognized by its low gamma value of around 920 ms. Figure 10 shows the results of the time-frequency analysis of the post-stacked seismic data in Figure 9. The time-frequency analysis profiles correspond to frequencies of 10, 30, 50, and 70 Hz. As shown in Figure 10, the variations in the strength reflected from the target layer at different frequencies form the basis for the inversion of the P-wave dispersion attributes.

Inversion of P-Wave Dispersion Attribute Dp
Based on the frequency-dependent AVO inversion theory, the dispersion attribute profile of Dp is calculated and shown in Figure 11. The red curve in Figure 11 is the logging gamma ray value. Strong P-wave velocity dispersion can be observed at the position of the target layer around 920 ms, as pointed out by an arrow.
However, in Figure 11, a high value of the Dp attribute can be observed just above the target zone. A reasonable interpretation on this may be the interbedded structure consisting of thin layers of mudstone and sandstone over the position of the target zone. Such interbedded layers that can be observed on logging data in Figure 8 mean that seismic events represent more phases (peaks and troughs) as shown in Figure 9. Strong reflection energy reveals the presence of interference between reflected waves.
The frequency-dependent attributes of seismic reflections are affected by both the interference caused by stacked thin layers, and by the dispersion and attenuation related to fluid flow. Thus, when we intend to identify reservoir fluids using frequencydependent attributes, it is necessary to notice this uncertainty.
When comparing Figures 9, 11, it is meaningful that the inverted Dp attribute reveals the location of the pay zone which has high gas saturation. Such a pay zone cannot be predicted on the post-stacked seismic profile, where the seismic events show fewer obvious variations laterally. Figure 12 shows a horizontal slice of the RMS amplitude of the seismic reflection from the target zone. Figure 13 is the  corresponding slice of P-wave velocity dispersion attribute Dp. The values of the dispersion attribute Dp for wells A, B, and C in the study area range from high to low. The development of the reservoir shows that the productivities of the three wells are consistent with the inverted Dp attributes, which verifies the feasibility of identifying gas-bearing tight sandstone reservoirs based on frequency-dependent AVO inversion.

DISCUSSION AND CONCLUSION
We conducted rock physics modeling of tight gas sandstone by employing the Chapman theory to simulate the dispersion and attenuation caused by the squirt flow of fluids in complex pore spaces. The frequency dependence of the AVO signatures and the inversion of the P-wave dispersion attribute were investigated in different theoretical models with various gas saturations in tight gas sandstone. The field data applications of the frequencydependent AVO inversion indicate that the estimated P-wave dispersion attribute can be used as a reliable indicator of gas identification in tight gas sandstone reservoirs.
In the constructed rock physics model, the complex pore space is equivalent to a combination of high aspect-ratio pores (stiff pores) and low aspect-ratio microcracks (soft pores). The velocity dispersion and attenuation were simulated by the squirt flow of fluids in the complex pore spaces, which is affected by the gas saturation and the viscosity of the fluid mixtures of gas and brine. The rock physics modeling revealed that as the gas saturation increases, the P-wave velocity generally decreases. P-wave velocity dispersion and attenuation occurs at seismic frequencies, and it tends to move toward higher frequencies as the gas saturation increases.
According to the rock physics modeling results, P-wave velocity dispersion in tight sandstone results in the frequencydependence of the contrast in the elastic impedance at the interface between the tight sandstone and the overlying mudstone, and therefore, it leads to frequency-dependent reflection coefficients. The synthetic model based on the propagator matrix method revealed that the reflection coefficient varies significantly with frequency, and the variations in the amplitudes and phases of the PP reflections can be observed for various gas saturations.
The tests of the frequency-dependent AVO inversion using synthetic data revealed that the magnitude of the P-wave dispersion attribute is an effective indicator of gas saturation in the tight sandstone. Finally, the frequency-dependent AVO inversion was applied to pre-stacked field seismic data. The results indicate that the inverted P-wave dispersion attribute is well correlated with the production of the pay zone for the evaluation wells. Thus, the methods investigated in this study have good potential for the evaluation of gas saturation in tight gas sandstone reservoirs.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
ZG contributed as the corresponding author of the manuscript. HJ did a part of writing and coding works. CL, YZ, CN, DW, and YL provided some interesting ideas. All authors contributed to the article and approved the submitted version.

FUNDING APPENDIX A
The Frequency-Dependent AVO Theory Shuey (1985) gave an approximate formula for the PP reflection coefficient for isotropic media: where θ is the incident angle. V p and ρ are the mean values of the P-wave velocity and density across the interface, respectively. V p and ρ are the differences in the P-wave velocity and density across the interface, respectively. G is the AVO gradient. Based on Eq. A1, Chapman et al. (2006) concluded that the dispersion of the P-wave velocity caused the frequency dependency of the reflection coefficient. Thus, Eq. A1 can be rewritten as Using the Taylor series expansion method, we obtained the frequency-dependent R PP for a reference frequency: where D P represents the derivatives of the seismic wave velocities with frequency; and D G is the frequency-dependent AVO gradient.
In the application of the frequency-dependent AVO method to real data, the time-frequency spectra S(t,θ,f) of the pre-stacked gathers contains the information for the incidence wavelet spectrum, so it is necessary to eliminate the effect of the incidence wavelet.
M t, θ, f = S t, θ, f w f , θ (A6) where max A f =ref (θ) is the maximum value of the amplitude spectrum in the selected time window used for the calculation. According to Eq. A4, where M(t,θ,f )=M(t,θ,f )−M(t,θ,f0 ) . Eq. A8 can be rewritten in the form of a matrix: where the D P and D G values of inversion parameters can be obtained by solving the over-determined equation.