# Fluid Discrimination Based on Frequency-Dependent AVO Inversion with the Elastic Parameter Sensitivity Analysis

Authors: Pu Wang; Jingye Li; Xiaohong Chen; Kedong Wang; Benfeng Wang
Journal: Geofluids
 (2019)

---

## Abstract

Fluid discrimination is an extremely important part of seismic data interpretation. It plays an important role in the refined description of hydrocarbon-bearing reservoirs. The conventional AVO inversion based on Zoeppritz’s equation shows potential in lithology prediction and fluid discrimination; however, the dispersion and attenuation induced by pore fluid are not fully considered. The relationship between dispersion terms in different frequency-dependent AVO equations has not yet been discussed. Following the arguments of Chapman, the influence of pore fluid on elastic parameters is analyzed in detail. We find that the dispersion and attenuation of Russell fluid factor, Lamé parameter, and bulk modulus are more pronounced than those of P-wave modulus. The Russell fluid factor is most prominent among them. Based on frequency-dependent AVO inversion, the uniform expression of different dispersion terms of these parameters is derived. Then, incorporating the P-wave difference with the dispersion terms, we obtain new P-wave difference dispersion factors which can identify the gas-bearing reservoir location better compared with the dispersion terms. Field data application also shows that the dispersion term of Russell fluid factor is optimal in identifying fluid. However, the dispersion term of Russell fluid factor could be unsatisfactory, if the value of the weighting parameter associated with dry rock is improper. Then, this parameter is studied to propose a reasonable setting range. The results given by this paper are helpful for the fluid discrimination in hydrocarbon-bearing rocks.


---

## Body

## 1. Introduction

Rock physics is an effective tool for studying reservoir petrophysical properties from elastic and anelastic parameters [1]. Effective elastic media theory and poroelasticity theory are very important components of rock physics theories. The difference is that whether the mobility of pore fluid is considered or not [1]. The former considers the effects of pore geometry [2–4], instead of the effects of fluid flow on velocity. The Biot-Gassmann theory [5] is the research basis of most poroelasticity theories which are used to study the wave propagation in fluid-saturated rocks. The Gassmann theorem is proved to be valid for many types of rocks [6, 7], and it is used to predict P- and S-wave velocities in the low-frequency band [8, 9]. Due to the influence of pore fluid, dispersion and attenuation occur in the elastic parameters of rocks [10–12]. It is the same in the seismic frequency band, and some researchers have described the dispersion and attenuation by rock physics modeling. White proposed a model for seismic wave attenuation which is caused by the patchy saturation of two kinds of pore fluid [13]. The White model is extended afterward to be incorporated with Biot’s theory [14] or to characterize the complex pore structure [15]. Chapman et al. [16] presented a model which considers the squirt flow related to cracks. This model is consistent with the Gassmann’s relations at low frequency. Considering the anisotropy caused by the parallel arrangement of the cracks, Chapman [17] derived an anisotropic attenuation model. Wu et al. [18] gave three key assumptions: the model is consistent with the Gassmann-Wood prediction, partial saturation brings about higher attenuation value than full saturation, and the characteristic frequency on fluid mobility follows the theory of squirt flow. Based on these arguments, the Chapman model [16] is used to analyze the dispersion and attenuation of elastic parameters in the seismic frequency band in this paper, and the cracks are considered to be randomly arranged.Rock physics plays an important role in reservoir prediction and can be used to obtain subsurface petrophysical properties [19, 20]. The amplitude variation with offset (AVO) inversion is an effective method in reservoir interpretation and prediction [21–23]. The exact and approximate equations of Zoeppritz’s equation are adopted for the AVO inversion [24–26]. The approximate equations have good linearity and clear physical meaning [27–30] and can be more effectively applied to reservoir prediction. Exact Zoeppritz’s equation is a nonlinear forward operator, which may increase computation cost, while it is more suitable for far offset information [31]. With the development of exploration and exploitation, the lithologic trap reservoir has drawn much attention, and accurate description for hydrocarbon-bearing rocks is required. It brings huge challenges to fluid discrimination. There are mainly two methods for fluid discrimination by incorporating rock physics and seismic inversion [32]. One is to predict reservoir petrophysical properties directly by using the empirical or theoretical rock physics model [33]. The other is to identify fluid indirectly by using elastic parameters [34, 35]. The prediction precision of the former method is based on the accuracy of the rock physics model, which is the quantitative interpretation method. The latter depends on the understanding of the relationship between petrophysical properties and elastic parameters, which is the qualitative interpretation method. The AVO inversion shows a certain ability in lithology prediction and fluid discrimination [36, 37]. However, the elastic parameter inversion based on Zoeppritz’s equation does not consider the effect of dispersion and attenuation caused by pore fluid. Wilson et al. proposed frequency-dependent AVO (FAVO) following the approximation of Smith and Gidlow [38]. And the dispersion term used for fluid discrimination can be obtained by using FAVO inversion. The resolution of spectral decomposition can affect the accuracy of FAVO inversion. Wu et al. [39] and Luo et al. [40] used the smooth pseudo Wigner-Ville distribution and inverse spectral decomposition for the time-frequency analysis, respectively, to enhance the resolution and accuracy of FAVO inversion results. In order to identify the hydrocarbon-bearing reservoir accurately, more sensitive frequency-dependent fluid factor needs to be studied [41]. Zhang et al. proposed a new fluid factor following the approximation of Russell et al. [42]. Chen et al. [43] and Wang et al. [44] suggested P-wave difference dispersion factor (PDDF) and bulk modulus difference dispersion factor (MDDF) to identify pore fluid, respectively.In this paper, the attenuation and dispersion of elastic parameters are analyzed in detail to identify which is more dispersion sensitive. The results show that the pore fluid can bring more attenuation to Russell fluid factor, Lamé parameter, and bulk modulus, compared with P-wave velocity. This is critical for us to decide which dispersion term to be used in the FAVO inversion. In addition, the applicability of Russell AVO approximation is analyzed, and a more reasonable range of the pending parameter is given, which can avoid a large deviation. Then, based on FAVO inversion, the dispersion terms and PDDFs related to Russell fluid factor, Lamé parameter, and bulk modulus are derived. The PDDFs are superior to dispersion terms in identifying the gas-bearing reservoir. In order to improve the antinoise ability of the inversion,ℓ1 norm regularization term [45] is applied.


## 2. Theory and Method

The current study of fluid identification is based on the dispersion and attenuation of elastic parameters. Therefore, the Chapman model [16] is used to analyze which parameter is more sensitive to the dispersion and attenuation induced by the pore fluid. Then, the theories of sparse constrained inversion spectral decomposition and frequency-dependent AVO inversion are summarized. In this section, the formula of frequency-dependent AVO is rewritten with respect to other elastic parameters and the applicability of Russell AVO approximation is analyzed.


### 2.1. Dispersion and Attenuation Analysis of Elastic Parameters

Chapman et al. [16] presented a frequency-dependent model which considers the wave-induced exchange of fluids between pores and cracks, as well as between cracks. The cracks in this model are randomly arranged without causing anisotropy. In the low-frequency limit, this model is consistent with the predicted result by using the Gassmann equation. By correcting the isotropic elastic tensor of the material with the perturbation due to the presence of cracks and pores, it defines a valid expression:
(1)Cijkl=Cijkl0λ,μ−Cijkl1λ,μ,ω,τ,ϕ,r,ε,κf,where Cijkl0• is the isotropic elastic tensor without inclusions and Cijkl1• is the term affected by cracks and pores; λ, μ, ω, τ, ϕ, r, ε, and κf are the Lamé parameter, shear modulus, angular frequency, timescale parameter, porosity, crack aspect ratio, crack density, and fluid bulk modulus, respectively. The detailed expression of Cijkl is given by Chapman et al. [46].According to equation (1), we can get the elastic parameters influenced by frequency and fluid. Following the arguments of Chapman et al., the S-wave velocity also has attenuation and dispersion. This is because some of the cracks will necessarily have an orientation such that they will be compressed due to applied shear stress [47]. The attenuation of S-wave velocity is smaller than that of P-wave velocity in the Chapman model. To obtain attenuation characteristics of different elastic parameters, other fluid-sensitive frequency-dependent elastic parameters are studied by removing the shear modulus from the P-wave modulus:
(2)λ=M−2μ,K=M−43μ,f=M−γdry2μ,where K, M, and f are the bulk modulus, P-wave modulus, and Russell fluid factor, respectively; γdry2 is the weighting parameter associated with the dry rock and γdry2=VP/VSdry2 [35]; VP indicates P-wave velocity, and VS indicates S-wave velocity; subscript dry indicates the dry rock. In the FAVO inversion, the dispersion terms which are induced by the pore fluid are related to these parameters. Thus, these parameters are analyzed before the introduction of FAVO inversion.The dispersion characteristics of elastic parameters are shown in Figure1. It shows that the shear modulus is affected by the frequency component quite weakly in both the sandstone (SH=0.2) and the shale (SH=0.8). For a particular lithology, the change of pore fluid has a great influence on Russell fluid factor (f), Lamé parameter (λ), bulk modulus (K), and P-wave modulus (M), instead of the shear modulus (μ). The dispersion brought by gas (SG=0.6) is much larger than that of oil (SO=0.6). The value of γdry2 is discussed by Russell et al. [35]. We choose an intermediate value 2.3 of the values they listed. This parameter is also discussed later in the paper.Figure 1
Dispersion analysis of different elastic parameters. SH indicates shale content,SH=0.2 corresponds to sandstone reservoir and SH=0.8 corresponds to shale reservoir; SG indicates gas saturation in gas-water mixed fluid and SO indicates oil saturation in oil-water mixed fluid; K, μ, M, λ, and f indicate the bulk modulus, shear modulus, P-wave modulus, Lamé parameter, and Russell fluid factor, respectively; the value of γdry2 is given as 2.3. (a) Gas-bearing shale reservoir; (b) oil-bearing shale reservoir; (c) gas-bearing sandstone reservoir; (d) oil-bearing sandstone reservoir.
(a)
(b)
(c)
(d)Figure2 shows the attenuation of each elastic parameter under different reservoir petrophysical properties. The used parameter values of this model are summarized in Table 1. Since the dispersion and attenuation described by this rock physics model are caused by the effect of the pore fluid, the value of attenuation here corresponds to the dispersion magnitude. Figure 2 indicates that the attenuation values of Russell fluid factor, Lamé parameter (λ), bulk modulus (K), P-wave modulus (M), and shear modulus (μ) of the hydrocarbon-bearing reservoir are distributed from large to small. The elastic parameter attenuation value of the shale reservoir (SH=0.8) is larger than that of the sandstone reservoir (SH=0.2). The attenuation of the gas-bearing reservoir (SG=0.6) is larger than that of the oil-bearing reservoir (SO=0.6). It also shows that the attenuation of shear modulus is not affected by the pore fluid at low frequency. The attenuation of P-wave modulus is the same as that of P-wave velocity. Thus, the attenuation of Russell fluid factor, Lamé parameter, and bulk modulus is larger than that of P-wave velocity and P-wave modulus. Finally, we can obtain the dispersion terms of Russell fluid factor, Lamé parameter, or bulk modulus to improve the precision of FAVO inversion.Figure 2
Attenuation analysis of different elastic parameters. SH indicates shale content,SH=0.2 corresponds to sandstone reservoir and SH=0.8 corresponds to shale reservoir; SG indicates gas saturation in gas-water mixed fluid and SO indicates oil saturation in oil-water mixed fluid; K, μ, M, λ, and f indicate the bulk modulus, shear modulus, P-wave modulus, Lamé parameter, and Russell fluid factor, respectively; the value of γdry2 is given as 2.3. (a) Gas-bearing shale reservoir; (b) oil-bearing shale reservoir; (c) gas-bearing sandstone reservoir; (d) oil-bearing sandstone reservoir.
(a)
(b)
(c)
(d)Table 1
The used parameter values in the Chapman model.
Mineral/fluid
K (GPa)
μ (GPa)
ρ (g/cm3)
Parameters
Value
Sand
37
44
2.65
Time scale (τ)
5×10−3
Shale
21
7
2.60
Crack density (ε)
0.1
Water
2.62
0
1.01
Porosity (ϕ)
0.25
Oil
1.44
0
0.85
Crack aspect ratio (r)
0.001
Gas
0.05
0
0.21


### 2.2. Sparse Constrained Inversion Spectral Decomposition

In order to use the dispersion terms to identify the hydrocarbon-bearing rocks, the time-frequency distribution should be obtained as a prerequisite. According to the previous study [40, 44], inverse spectral decomposition has a higher time-frequency resolution. It has higher accuracy compared with other methods including short-time Fourier transform, wavelet transform, generalized S transform and Wigner-Ville distribution. Spectral decomposition based on sparse constrained inversion is thus used to obtain the high-resolution time-frequency distribution.Using a dictionary of Ricker wavelets with different central frequencies, the seismic data can be decomposed into corresponding frequency-dependent pseudoreflectivity. The convolutional model can be represented by a linear system of equations [48], and the seismic trace s can be expressed as
(3)s=W1W2⋯Wlr1r2⋮rl=Rm,where Wi refers to the Ricker wavelet with central frequency fi and ri refers to the corresponding pseudoreflectivity sequence; R indicates the dictionary of Ricker wavelets and m indicates the vector containing the pseudoreflectivity sequences [40].Theℓ1 norm regularization term can be applied to obtain a high-resolution result [45, 49] and it can promote a sparser solution than ℓ2 norm. The objective function J can be established based on ℓ1 norm constraint as shown in equation (4). It can be solved by utilizing SPGL1 algorithm, which is used for one-norm regularized least squares and is suitable for problems that are large scale [50]:
(4)J=s−Rm2+ζm1,where ζ is the trade-off parameter which balances the sparsity and the misfit.


### 2.3. Frequency-Dependent AVO Inversion

FAVO is based on the approximation of Zoeppritz’s equation [38, 39]. Aki and Richards [27] derived a simplified form of P-wave reflectivity in terms of density, P- and S-wave velocities. Gray et al. [51] provided two approximation equations about Lamé parameter, bulk modulus, shear modulus, and density. Another approximation considering fluid factor is derived by Russell et al. [35]. From the appendix, the approximate form proposed by Russell et al. has good general applicability, which can reduce to other approximations by modifying the parameter γdry2.Figure3 illustrates the AVO features of these four AVO approximations and the value of γdry2 is given as 2.3. The relevant upper and lower elastic parameters of the interface are shown in Table 2. Although they have different expressions, the AVO results are almost identical. Therefore, the FAVO study conducted below is based on the same accuracy of different AVO approximations. Since γdry2 is usually not available, it should be set to a fixed value by users. Equation (A.6) will cause a large deviation when the value of γdry2 is very close to VP/VSsat2, which implies that there is no fluid component in the rock. According to Table 2, the value of VP/VSsat2 is 2.6213. Figure 4 shows that if the value of γdry2 is very close to 2.6213, the AVO curve deviates from the curve of Aki-Richards AVO approximation (f‐μ‐ρ‐2.622 and f‐μ‐ρ‐2.612). When the values of γdry2 and VP/VSsat2 are not close, as shown in Figure 4, no matter the value of γdry2 is smaller or larger than VP/VSsat2, the AVO curve can be consistent with that of Aki-Richards AVO approximation (f‐μ‐ρ‐2.100 and f‐μ‐ρ‐3.000). Therefore, we need to ensure that the value of γdry2 cannot be chosen as any possible value of VP/VSsat2 in the reservoir during the AVO inversion to prevent the values of γdry2 and VP/VSsat2 from being very close. Since Poisson’s ratio of rocks must be larger than 0 and VP/VSsat2 must be larger than 2, therefore, equations (A.1), (A.4), and (A.5) are applicable to all the reservoirs without this deviation. An alternative way to avoid this deviation is to obtain all the potential values of VP/VSsat2 from the well logging data or the laboratory measurement, and then, take a value that is smaller or larger than them. However, by the application analysis of dispersion term later, only the smaller value is reasonable.Figure 3
Reflection amplitude analysis with incident angle of different AVO approximations.VP‐VS‐ρ indicates the Aki-Richards AVO approximation as equation (A.1) shows; λ‐μ‐ρ indicates the AVO approximation as equation (A.4) shows; K‐μ‐ρ indicates the AVO approximation as equation (A.5) shows; f‐μ‐ρ indicates the Russell et al. AVO approximation as equation (A.6) shows.Table 2
Parameters of the upper and lower layers.
P-wave velocity (m/s)
S-wave velocity (m/s)
Density (g/cm3)
Upper layer
3300
2000
2.2
Lower layer
3500
2200
2.3Figure 4
Reflection amplitude analysis with incident angle of differentγdry2. VP‐VS‐ρ indicates the Aki-Richards AVO approximation; f‐μ‐ρ‐2.622, f‐μ‐ρ‐2.612, f‐μ‐ρ‐2.100, and f‐μ‐ρ‐3.000 correspond to the Russell et al. AVO approximation with γdry2 equaling to 2.622, 2.612, 2.100, and 3.000, respectively.Wilson et al. derived a frequency-dependent AVO formula following the approximation of Smith and Gidlow [38]. This approximation eliminates the density term by using the Gardner formula. The density is not affected by the frequency, and the velocity is reversed. Thus, it is not recommended to use the Gardner formula in frequency-dependent AVO inversion. Other AVO approximations containing density term, such as equations (A.1), (A.4), (A.5), and (A.6), are utilized here. Besides, the dispersion term of Russell fluid factor has been studied by some researchers [42]. These AVO approximations can be rewritten in a uniform formula:
(5)RPPθ=AθΔXX+BθΔYY+CθΔρρ,where Aθ, Bθ, and Cθ are coefficients related to θ; ΔX/X and ΔY/Y are elastic parameter reflectivity of X and Y, respectively.Following the arguments of Wilson, we allowX and Y to vary with frequency ω:
(6)RPPθ,ω=AθΔXXω+BθΔYYω+CθΔρρ.Due to the narrow frequency band in seismic data, this equation can be expanded as a Taylor series around a representative frequencyω0:
(7)RPPθ,ω−RPPθ,ω0=ω−ω0AθIa+ω−ω0BθIb,where Ia=d/dωΔX/X and Ib=d/dωΔY/Y.By using sparse constrained inversion spectral decomposition, the time-frequency spectrumSt,θ,ω of seismic data Dt,θ can be obtained. The spectral decomposition result S can be affected by the overprint of source wavelet [52]. Therefore, spectral balancing should be performed on the spectral amplitude [53] to obtain the true spectral behaviour coming from the effect of geology and saturating fluid [18] by designing a suitable weight function wω:
(8)BRt,θ,ω=St,θ,ωwω.According to Wu et al. [18], wω is given as follows:
(9)wω=∑kS2t,θ,ωdom∑kS2t,θ,ω,where ωdom denotes the dominant frequency of the seismic data; k is the number of sampling points in a defined window.For the sample pointt0, the dispersion terms Ia and Ib can be obtained by solving the following equations:
(10)BRt0,θ,ω−RPPt0,θ,ωdom=ω−ωdomAt0,θIa+ω−ωdomBt0,θIb,RPPt0,θ,ωdom=At0,θΔXXt0,ωdom+Bt0,θΔYYt0,ωdom+Ct0,θΔρρ.A simplified matrix form can be characterized by(11)T=DIaIb.The matricesT and D are given as follows:
(12)T=BRt0,θ1,ω1−RPPt0,θ1,ωdom⋮BRt0,θ1,ωm−RPPt0,θ1,ωdom⋮BRt0,θn,ω1−RPPt0,θn,ωdom⋮BRt0,θn,ωm−RPPt0,θn,ωdom,D=ω1−ωdomAt0,θ1ω1−ωdomBt0,θ1⋮⋮ωm−ωdomAt0,θ1ωm−ωdomBt0,θ1⋮⋮ω1−ωdomAt0,θnω1−ωdomBt0,θn⋮⋮ωm−ωdomAt0,θnωm−ωdomBt0,θn,where n is the number of angle gathers and m is the number of frequencies after the spectral decomposition.


## 2.1. Dispersion and Attenuation Analysis of Elastic Parameters

Chapman et al. [16] presented a frequency-dependent model which considers the wave-induced exchange of fluids between pores and cracks, as well as between cracks. The cracks in this model are randomly arranged without causing anisotropy. In the low-frequency limit, this model is consistent with the predicted result by using the Gassmann equation. By correcting the isotropic elastic tensor of the material with the perturbation due to the presence of cracks and pores, it defines a valid expression:
(1)Cijkl=Cijkl0λ,μ−Cijkl1λ,μ,ω,τ,ϕ,r,ε,κf,where Cijkl0• is the isotropic elastic tensor without inclusions and Cijkl1• is the term affected by cracks and pores; λ, μ, ω, τ, ϕ, r, ε, and κf are the Lamé parameter, shear modulus, angular frequency, timescale parameter, porosity, crack aspect ratio, crack density, and fluid bulk modulus, respectively. The detailed expression of Cijkl is given by Chapman et al. [46].According to equation (1), we can get the elastic parameters influenced by frequency and fluid. Following the arguments of Chapman et al., the S-wave velocity also has attenuation and dispersion. This is because some of the cracks will necessarily have an orientation such that they will be compressed due to applied shear stress [47]. The attenuation of S-wave velocity is smaller than that of P-wave velocity in the Chapman model. To obtain attenuation characteristics of different elastic parameters, other fluid-sensitive frequency-dependent elastic parameters are studied by removing the shear modulus from the P-wave modulus:
(2)λ=M−2μ,K=M−43μ,f=M−γdry2μ,where K, M, and f are the bulk modulus, P-wave modulus, and Russell fluid factor, respectively; γdry2 is the weighting parameter associated with the dry rock and γdry2=VP/VSdry2 [35]; VP indicates P-wave velocity, and VS indicates S-wave velocity; subscript dry indicates the dry rock. In the FAVO inversion, the dispersion terms which are induced by the pore fluid are related to these parameters. Thus, these parameters are analyzed before the introduction of FAVO inversion.The dispersion characteristics of elastic parameters are shown in Figure1. It shows that the shear modulus is affected by the frequency component quite weakly in both the sandstone (SH=0.2) and the shale (SH=0.8). For a particular lithology, the change of pore fluid has a great influence on Russell fluid factor (f), Lamé parameter (λ), bulk modulus (K), and P-wave modulus (M), instead of the shear modulus (μ). The dispersion brought by gas (SG=0.6) is much larger than that of oil (SO=0.6). The value of γdry2 is discussed by Russell et al. [35]. We choose an intermediate value 2.3 of the values they listed. This parameter is also discussed later in the paper.Figure 1
Dispersion analysis of different elastic parameters. SH indicates shale content,SH=0.2 corresponds to sandstone reservoir and SH=0.8 corresponds to shale reservoir; SG indicates gas saturation in gas-water mixed fluid and SO indicates oil saturation in oil-water mixed fluid; K, μ, M, λ, and f indicate the bulk modulus, shear modulus, P-wave modulus, Lamé parameter, and Russell fluid factor, respectively; the value of γdry2 is given as 2.3. (a) Gas-bearing shale reservoir; (b) oil-bearing shale reservoir; (c) gas-bearing sandstone reservoir; (d) oil-bearing sandstone reservoir.
(a)
(b)
(c)
(d)Figure2 shows the attenuation of each elastic parameter under different reservoir petrophysical properties. The used parameter values of this model are summarized in Table 1. Since the dispersion and attenuation described by this rock physics model are caused by the effect of the pore fluid, the value of attenuation here corresponds to the dispersion magnitude. Figure 2 indicates that the attenuation values of Russell fluid factor, Lamé parameter (λ), bulk modulus (K), P-wave modulus (M), and shear modulus (μ) of the hydrocarbon-bearing reservoir are distributed from large to small. The elastic parameter attenuation value of the shale reservoir (SH=0.8) is larger than that of the sandstone reservoir (SH=0.2). The attenuation of the gas-bearing reservoir (SG=0.6) is larger than that of the oil-bearing reservoir (SO=0.6). It also shows that the attenuation of shear modulus is not affected by the pore fluid at low frequency. The attenuation of P-wave modulus is the same as that of P-wave velocity. Thus, the attenuation of Russell fluid factor, Lamé parameter, and bulk modulus is larger than that of P-wave velocity and P-wave modulus. Finally, we can obtain the dispersion terms of Russell fluid factor, Lamé parameter, or bulk modulus to improve the precision of FAVO inversion.Figure 2
Attenuation analysis of different elastic parameters. SH indicates shale content,SH=0.2 corresponds to sandstone reservoir and SH=0.8 corresponds to shale reservoir; SG indicates gas saturation in gas-water mixed fluid and SO indicates oil saturation in oil-water mixed fluid; K, μ, M, λ, and f indicate the bulk modulus, shear modulus, P-wave modulus, Lamé parameter, and Russell fluid factor, respectively; the value of γdry2 is given as 2.3. (a) Gas-bearing shale reservoir; (b) oil-bearing shale reservoir; (c) gas-bearing sandstone reservoir; (d) oil-bearing sandstone reservoir.
(a)
(b)
(c)
(d)Table 1
The used parameter values in the Chapman model.
Mineral/fluid
K (GPa)
μ (GPa)
ρ (g/cm3)
Parameters
Value
Sand
37
44
2.65
Time scale (τ)
5×10−3
Shale
21
7
2.60
Crack density (ε)
0.1
Water
2.62
0
1.01
Porosity (ϕ)
0.25
Oil
1.44
0
0.85
Crack aspect ratio (r)
0.001
Gas
0.05
0
0.21


## 2.2. Sparse Constrained Inversion Spectral Decomposition

In order to use the dispersion terms to identify the hydrocarbon-bearing rocks, the time-frequency distribution should be obtained as a prerequisite. According to the previous study [40, 44], inverse spectral decomposition has a higher time-frequency resolution. It has higher accuracy compared with other methods including short-time Fourier transform, wavelet transform, generalized S transform and Wigner-Ville distribution. Spectral decomposition based on sparse constrained inversion is thus used to obtain the high-resolution time-frequency distribution.Using a dictionary of Ricker wavelets with different central frequencies, the seismic data can be decomposed into corresponding frequency-dependent pseudoreflectivity. The convolutional model can be represented by a linear system of equations [48], and the seismic trace s can be expressed as
(3)s=W1W2⋯Wlr1r2⋮rl=Rm,where Wi refers to the Ricker wavelet with central frequency fi and ri refers to the corresponding pseudoreflectivity sequence; R indicates the dictionary of Ricker wavelets and m indicates the vector containing the pseudoreflectivity sequences [40].Theℓ1 norm regularization term can be applied to obtain a high-resolution result [45, 49] and it can promote a sparser solution than ℓ2 norm. The objective function J can be established based on ℓ1 norm constraint as shown in equation (4). It can be solved by utilizing SPGL1 algorithm, which is used for one-norm regularized least squares and is suitable for problems that are large scale [50]:
(4)J=s−Rm2+ζm1,where ζ is the trade-off parameter which balances the sparsity and the misfit.


## 2.3. Frequency-Dependent AVO Inversion

FAVO is based on the approximation of Zoeppritz’s equation [38, 39]. Aki and Richards [27] derived a simplified form of P-wave reflectivity in terms of density, P- and S-wave velocities. Gray et al. [51] provided two approximation equations about Lamé parameter, bulk modulus, shear modulus, and density. Another approximation considering fluid factor is derived by Russell et al. [35]. From the appendix, the approximate form proposed by Russell et al. has good general applicability, which can reduce to other approximations by modifying the parameter γdry2.Figure3 illustrates the AVO features of these four AVO approximations and the value of γdry2 is given as 2.3. The relevant upper and lower elastic parameters of the interface are shown in Table 2. Although they have different expressions, the AVO results are almost identical. Therefore, the FAVO study conducted below is based on the same accuracy of different AVO approximations. Since γdry2 is usually not available, it should be set to a fixed value by users. Equation (A.6) will cause a large deviation when the value of γdry2 is very close to VP/VSsat2, which implies that there is no fluid component in the rock. According to Table 2, the value of VP/VSsat2 is 2.6213. Figure 4 shows that if the value of γdry2 is very close to 2.6213, the AVO curve deviates from the curve of Aki-Richards AVO approximation (f‐μ‐ρ‐2.622 and f‐μ‐ρ‐2.612). When the values of γdry2 and VP/VSsat2 are not close, as shown in Figure 4, no matter the value of γdry2 is smaller or larger than VP/VSsat2, the AVO curve can be consistent with that of Aki-Richards AVO approximation (f‐μ‐ρ‐2.100 and f‐μ‐ρ‐3.000). Therefore, we need to ensure that the value of γdry2 cannot be chosen as any possible value of VP/VSsat2 in the reservoir during the AVO inversion to prevent the values of γdry2 and VP/VSsat2 from being very close. Since Poisson’s ratio of rocks must be larger than 0 and VP/VSsat2 must be larger than 2, therefore, equations (A.1), (A.4), and (A.5) are applicable to all the reservoirs without this deviation. An alternative way to avoid this deviation is to obtain all the potential values of VP/VSsat2 from the well logging data or the laboratory measurement, and then, take a value that is smaller or larger than them. However, by the application analysis of dispersion term later, only the smaller value is reasonable.Figure 3
Reflection amplitude analysis with incident angle of different AVO approximations.VP‐VS‐ρ indicates the Aki-Richards AVO approximation as equation (A.1) shows; λ‐μ‐ρ indicates the AVO approximation as equation (A.4) shows; K‐μ‐ρ indicates the AVO approximation as equation (A.5) shows; f‐μ‐ρ indicates the Russell et al. AVO approximation as equation (A.6) shows.Table 2
Parameters of the upper and lower layers.
P-wave velocity (m/s)
S-wave velocity (m/s)
Density (g/cm3)
Upper layer
3300
2000
2.2
Lower layer
3500
2200
2.3Figure 4
Reflection amplitude analysis with incident angle of differentγdry2. VP‐VS‐ρ indicates the Aki-Richards AVO approximation; f‐μ‐ρ‐2.622, f‐μ‐ρ‐2.612, f‐μ‐ρ‐2.100, and f‐μ‐ρ‐3.000 correspond to the Russell et al. AVO approximation with γdry2 equaling to 2.622, 2.612, 2.100, and 3.000, respectively.Wilson et al. derived a frequency-dependent AVO formula following the approximation of Smith and Gidlow [38]. This approximation eliminates the density term by using the Gardner formula. The density is not affected by the frequency, and the velocity is reversed. Thus, it is not recommended to use the Gardner formula in frequency-dependent AVO inversion. Other AVO approximations containing density term, such as equations (A.1), (A.4), (A.5), and (A.6), are utilized here. Besides, the dispersion term of Russell fluid factor has been studied by some researchers [42]. These AVO approximations can be rewritten in a uniform formula:
(5)RPPθ=AθΔXX+BθΔYY+CθΔρρ,where Aθ, Bθ, and Cθ are coefficients related to θ; ΔX/X and ΔY/Y are elastic parameter reflectivity of X and Y, respectively.Following the arguments of Wilson, we allowX and Y to vary with frequency ω:
(6)RPPθ,ω=AθΔXXω+BθΔYYω+CθΔρρ.Due to the narrow frequency band in seismic data, this equation can be expanded as a Taylor series around a representative frequencyω0:
(7)RPPθ,ω−RPPθ,ω0=ω−ω0AθIa+ω−ω0BθIb,where Ia=d/dωΔX/X and Ib=d/dωΔY/Y.By using sparse constrained inversion spectral decomposition, the time-frequency spectrumSt,θ,ω of seismic data Dt,θ can be obtained. The spectral decomposition result S can be affected by the overprint of source wavelet [52]. Therefore, spectral balancing should be performed on the spectral amplitude [53] to obtain the true spectral behaviour coming from the effect of geology and saturating fluid [18] by designing a suitable weight function wω:
(8)BRt,θ,ω=St,θ,ωwω.According to Wu et al. [18], wω is given as follows:
(9)wω=∑kS2t,θ,ωdom∑kS2t,θ,ω,where ωdom denotes the dominant frequency of the seismic data; k is the number of sampling points in a defined window.For the sample pointt0, the dispersion terms Ia and Ib can be obtained by solving the following equations:
(10)BRt0,θ,ω−RPPt0,θ,ωdom=ω−ωdomAt0,θIa+ω−ωdomBt0,θIb,RPPt0,θ,ωdom=At0,θΔXXt0,ωdom+Bt0,θΔYYt0,ωdom+Ct0,θΔρρ.A simplified matrix form can be characterized by(11)T=DIaIb.The matricesT and D are given as follows:
(12)T=BRt0,θ1,ω1−RPPt0,θ1,ωdom⋮BRt0,θ1,ωm−RPPt0,θ1,ωdom⋮BRt0,θn,ω1−RPPt0,θn,ωdom⋮BRt0,θn,ωm−RPPt0,θn,ωdom,D=ω1−ωdomAt0,θ1ω1−ωdomBt0,θ1⋮⋮ωm−ωdomAt0,θ1ωm−ωdomBt0,θ1⋮⋮ω1−ωdomAt0,θnω1−ωdomBt0,θn⋮⋮ωm−ωdomAt0,θnωm−ωdomBt0,θn,where n is the number of angle gathers and m is the number of frequencies after the spectral decomposition.


## 3. Application

The used field seismic data set is from HampsonRussell software, recorded over the shallow Cretaceous. Figure5 shows the stack section with gathers from CDP 260 to 390. The angles vary from 3° to 24°. The dominant frequency is about 35 Hz. The red circle in Figure 5 indicates the location of gas reservoir with the strong reflector. A real gas well is located at CDP 330. From the well logging data, the lithology of the target reservoir is gas-bearing sandstone and the gas saturation is about 50%. This is a typical gas reservoir, which can be used to test whether our method is valid or not. This set of data has been used as an application study by other researchers [54].Figure 5
Stack section of a field dataset from HampsonRussell software. The seismic traces are shown with the intervals of 3. The red circle highlights the location of the gas reservoir.In order to obtainVS/VPsat2 in the FAVO inversion, the P- and S-wave velocities are estimated based on the exact Zoeppritz’s equation. The Bayesian framework is used to improve the stability of the inversion. For simplicity, VS/VPsat2 can be approximated to a fixed value with some accuracy loss. Figure 6 shows the inversion results. It shows that the accuracy of P-wave velocity around the red circle is suppressed due to the existence of the gas reservoir compared with the overlying strata. Figure 7(d) shows relatively weaker dispersion of P-wave velocity in the gas reservoir. There is also a deviation in the position of the strong dispersion which cannot match the reservoir location very well. Figures 7(a)–7(c) (γdry2=2.3) have the same dispersion characteristics and can better identify the fluid. The dispersion of Figure 7(b) is slightly weaker than that of Figures 7(a) and 7(c). Figure 7 shows that the dispersion terms related to Russell fluid factor, Lamé parameter, and bulk modulus are superior to that of P-wave velocity, and among them, Russell fluid factor is slightly more prominent, followed by Lamé parameter.Figure 6
Inversion results based on the exact Zoeppritz’s equation. (a) is the section of P-wave velocity, (b) is the section of S-wave velocity, and (c) is the section of density. The real well logs are inserted at CDP 330. The red circle highlights the location of the gas reservoir.
(a)
(b)
(c)Figure 7
The results of FAVO inversion.Ia indicates the dispersion term. (a) is the section of dispersion term of Lamé parameter, (b) is the section of dispersion term of bulk modulus, (c) is the section of dispersion term of Russell fluid factor, and (d) is the section of dispersion term of P-wave velocity. The red circle highlights the location of the gas reservoir. The dispersion terms are normalized.
(a)
(b)
(c)
(d)The value ofγdry2 is determined from the well logging data. Figure 8 shows that all values of VP/VSsat2 are larger than 2.3. If we let γdry2=4.1 or 8.9, the dispersion term of Russell fluid factor is unsatisfactory (Figures 9(a) and 9(b)). These two values of γdry2 are within the range of VP/VSsat2, and they can cause big deviations to the dispersion term. According to the previous analysis, the AVO approximation is accurate when the value of γdry2 is much larger than VP/VSsat2. However, this does not mean that the dispersion term is correct. Figure 9(c) shows poor dispersion term even for a large value of γdry2 (γdry2=15). Therefore, γdry2 should be smaller than any possibility of VP/VSsat2. In addition, from equations (2), if the value of γdry2 is 2, 4/3, or 0, Russell fluid factor f can be changed into Lamé parameter λ, bulk modulus K, or P-wave modulus M. As depicted in Figure 2, with a larger value of γdry2, the Russell fluid factor can have more attenuation induced by the pore fluid.Figure 8
Well logging curves. (a) indicates P-wave velocity; (b) indicates S-wave velocity; (c) indicatesVP/VSsat2.Figure 9
Analysis of dispersion terms of Russell fluid factor with differentγdry2. (a) corresponds to the dispersion term with γdry2=4.1; (b) corresponds to the dispersion term with γdry2=8.9; (c) corresponds to the dispersion term with γdry2=15; the red circle highlights the location of the gas reservoir.
(a)
(b)
(c)According to Chen et al. [43], we can derive new P-wave difference dispersion factor (PDDF) combining with the P-wave difference. The new PDDFs are related to Russell fluid factor, Lamé parameter, and bulk modulus by replacing Ia in the following equation:
(13)PDDF=ΔVP⋅Ia,where ΔVP is obtained by subtracting the P-wave velocity of the upper layer from that of the lower layer. Figure 10 shows the prediction process of dispersion term and P-wave difference dispersion factor.Figure 10
The prediction process of dispersion term and P-wave difference dispersion factor.Figure11 indicates that the PDDF is superior to dispersion term Ia (Figure 7) in fluid discrimination. The section of PDDF has better accuracy and better continuity. Due to the difference in dispersion terms, the PDDF results are different obviously. The PDDFs related to Russell fluid factor, Lamé parameter, and bulk modulus are almost similar, and they are all better than that of P-wave velocity in identifying the gas-bearing rocks. By comparison, Russell fluid factor is slightly more obvious.Figure 11
The P-wave difference dispersion factors of different dispersion terms. PDDF indicates the P-wave difference dispersion factor; (a) is the section of PDDF related to Lamé parameter, (b) is the section of PDDF related to bulk modulus, (c) is the section of PDDF related to Russell fluid factor (γdry2=2.3), and (d) is the section of PDDF related to P-wave velocity. The red circle highlights the location of the gas reservoir. The PDDFs are normalized.
(a)
(b)
(c)
(d)The various parameter curves at borehole-side trace are displayed in Figure12. In Figure 12, there are two peaks on the curves of dispersion term and PDDF at the location of gas reservoir. It is consistent with the analysis of HampsonRussell software about the two sets of gas layer: top gas and base gas. The dispersion terms and PDDFs related to Russell fluid factor (f), Lamé parameter (λ), and bulk modulus (K) are more pronounced. Due to the attenuation and dispersion caused by the gas reservoir, the inverted P-wave velocity (Figure 12(b)) based on the elastic theory of Zoeppritz’s equation is affected and the prediction accuracy of the lower layer is unsatisfactory. This also explains why FAVO inversion can be used to identify oil and gas.Figure 12
Display of various parameters at the well location. (a) is the poststack seismic trace; (b) is the P-wave velocity, the red curve is the real well log after smoothing, and the blue curve is the inverted result; (c) is the dispersion term of Lamé parameter (λ); (d) is the dispersion term of bulk modulus (K); (e) is the dispersion term of Russell fluid factor (f) (γdry2=2.3); (f) is the dispersion term of P-wave velocity (VP); (g) is the PDDF of Lamé parameter (λ); (h) is the PDDF of bulk modulus (K); (i) is the PDDF of Russell fluid factor (f) (γdry2=2.3); (j) is the PDDF of P-wave velocity (VP). The yellow dotted rectangle highlights the location of the gas reservoir.


## 4. Discussion and Conclusions

Pore fluid can cause dispersion and attenuation of elastic parameters in the seismic frequency range. Therefore, the study of dispersion term of elastic parameter can be beneficial for fluid discrimination. By analysis, the shear modulus is almost unaffected by the pore fluid; however, the dispersion and attenuation of Russell fluid factor, Lamé parameter, and bulk modulus are stronger and decrease in turn. Uniform expression of dispersion terms of these parameters is derived, and P-wave difference dispersion factors are studied by incorporating with P-wave difference.The applicability of Russell AVO approximation is analyzed, and an inappropriate weighting parameter associated with the dry rock can cause a large deviation. In order to ensure that the AVO attribute and the dispersion term are correct and reasonable, the valueγdry2 should be smaller than any possibility of VP/VSsat2 from the well logging data or the laboratory measurement.Field data application shows that the dispersion terms related to Russell fluid factor with proper weighting parameter, Lamé parameter, and bulk modulus are preferable, and among them, Russell fluid factor is slightly more prominent, followed by Lamé parameter, which is consistent with the Chapman model. The P-wave difference dispersion factors related to Russell fluid factor, Lamé parameter, and bulk modulus are almost similar and better than those of P-wave velocity, and Russell fluid factor has a slight advantage in identifying pore fluid.



---

*Source: 8750127-2019-12-18.xml*