# Inversion for Geofluid Discrimination Based on Poroelasticity and AVO Inversion

Authors: Lingqian Wang; Hui Zhou; Bo Yu; Yanxin Zhou; Wenling Liu; Yukun Tian
Journal: Geofluids
 (2019)

---

## Abstract

Geofluid discrimination plays an important role in reservoir characterization and prospect identification. Compared with other fluid indicators, the effective pore-fluid bulk modulus is more sensitive to the property of fluid contained in reservoirs. We combine the empirical relations with deterministic models to form a new kind of linearized relationship between the mixed fluid/rock term and the fluid modulus. On the one hand, the linearized relationship can decouple the fluid bulk modulus from the mixed fluid/rock term; on the other hand, the decoupled terms are more stable especially in low-porosity situations compared with previous approaches. In terms of the new linearized equation of the fluid modulus, we derive a novel linearized amplitude variation with offset (AVO) approximation to avoid the complicated nonlinear relationship between the fluid modulus and the reflectivity series. Convoluting this linearized AVO approximation with seismic wavelets, the forward modeling is constructed to combine the prestack seismic records with the fluid modulus. Meanwhile, we introduce the Bayesian inference with multivariable Cauchy prior to the fluid modulus inversion for a stable and high-resolution solution. Model examples demonstrate the accuracy of the proposed linearized AVO approximation compared with the exact Zoeppritz equation and Aki-Richards approximate equation. The synthetic and field data tests illustrate the accuracy and feasibility of the proposed fluid modulus inversion approach for geofluid discrimination.


---

## Body

## 1. Introduction

Fluid discrimination plays an important role in seismic exploration and reservoir description. Qualitative interpretation [1–6] and direct quantitative estimation [7–11] for rock properties are two main methods to discriminate different pore-fluid types. However, the qualitative interpretation, such as the bright-spot, dim-spot, and flat-spot techniques, has difficulties in discriminating the geofluid in complicated lithologic reservoirs. Estimation of the rock-physical property directly from seismic records is a better way to differentiate the geofluid filled in reservoirs. The major challenge is the uncertainties in geofluid discrimination associated with two factors. First, the fluid indicators are very likely to provide ambiguous results for fluid identification due to the mixed effects of the fluid and rock porosity. Second, the instability and inaccuracy of the traditional elastic parameters, such as the P-wave velocity and S-wave velocity estimated by prestack seismic inversion, may magnify the uncertainty of the fluid indicators. A lot of efforts are taken to study various fluid indicators from prestack seismic records for reservoir prediction and fluid discrimination. Smith and Gidlow [12] initially incorporated the mudrock line of Castagna et al. [13] and the pseudo-Poisson’s ratio reflectivity to define the fluid factor in the form of the weighted difference between the P- and S-wave velocity reflectivities. Rutherford and Williams [2] determined gas sands based on the normal-incidence reflection coefficient and the contrast in Poisson’s ratio. Goodway et al. [14] converted the velocity measurements to Lame’s parameters to enhance the sensitivity to fluids. Quakenbush et al. [15] utilized the crossplot of acoustic impedance and shear impedance to investigate discrimination between different lithology and fluid types.In terms of the poroelasticity theory presented by Biot [16, 17] and Gassmann [18], Russell et al. [19] generalized the formulation of Goodway et al. [14] by defining the fluid indicator as the weighted difference between the P-wave impedance and the S-wave impedance. Zong et al. [20] defined the P- and S-wave moduli as a new fluid indicator and introduced them into the prestack seismic inversion. Russell et al. [21] and Zong et al. [22] derived a new AVO approximate equation with different approaches to estimate the fluid/porosity term directly. It may be that these geofluid indicators can estimate the fluid type effectively to a certain extent, but they are ambiguous in complicated situations due to the integrated response of the geofluid and the petrophysical properties. In theory, fluid modulus is the most sensitive parameter to discriminate the geofluid contained in reservoirs. Yin and Zhang [23] extracted the fluid modulus as the fluid factor based on the empirical critical porosity and estimated the fluid modulus directly from the prestack seismic data. Further, a rock-physical empirical formula [24, 25] and AVO theory were combined to generate a new linearized AVO approximate equation to estimate the fluid modulus directly [26]. However, one of the inversion terms containing the porosity in the denominator means that the inversion results will be unstable in low-porosity situations. Meanwhile, the hypothesis and approximation will magnify the instability and inaccuracy. In this paper, we combine the deterministic model proposed by Han and Batzle [27] with the velocity-porosity empirical relations to derive the linearized relationship between the fluid-saturation effect term in Gassmann’s equation and the fluid modulus. The stability of the inversion results can be improved and some approximations can be omitted for more accurate results.It is important to estimate the rock-physical properties for reservoir prediction and fluid discrimination. Empirical relations and deterministic models can provide the relationship between the petrophysical parameters and the elastic parameters, such as P-wave velocity (VP), S-wave velocity (VS), and density (ρ). Compared with the poststack seismic inversion, prestack seismic inversion can provide more elastic property information taking advantage of the AVO phenomenon [28]. However, the prestack inversion also suffers difficulties stemming from noise contamination, band limitation, and nonuniqueness like other inversion problems [29]. Therefore, it is important to improve the stability and accuracy of prestack seismic inversion results with different methods. AVO inversion techniques can be classified into two categories: the stochastic inversion approaches and the deterministic inversion approaches; each has its own merits and demerits. The former category is carried out with forward modeling time by time in terms of the prior knowledge derived from well logs and multivariate geostatistical modeling [30–35]. The major limitations of its widespread application are the expensive computation and the tough choice and evaluation of inversion results. In contrast, the deterministic inversion approach can provide a determinate solution with high computational efficiency [36–38]. The Bayesian approach, as a popular method for inversion, takes advantage of available prior information of inversion parameters and the assumed distribution of noise in observed data as constraints to stabilize the inversion results [39–41]. The main concern in this approach is choosing an appropriate prior probability distribution for the inversion parameters [42]. The solution of conventional Gaussian probability distribution lacks sparsity and hence cannot improve resolution. In this paper, we utilize the Bayesian linearized methods proposed by Alemie and Sacchi [43] to obtain stable and accurate inversion results. The objective function is formulated with the pore-fluid bulk modulus and other elastic parameters based on the Bayesian framework with the assumption that the likelihood model has a Gaussian probability distribution and the prior model has a Cauchy distribution for high-resolution inversion results.In this paper, the geofluid-discrimination approach is conducted by incorporating poroelasticity and AVO inversion. Firstly, we decouple the fluid modulus from the mixed fluid/rock term in Gassmann’s equation. The deterministic models provided by Han and Batzle [27] are combined with the linearized empirical velocity-porosity relations to derive the fluid modulus. Compared with the decoupled approaches in Zong et al. [26], the decoupled terms in this paper are more stable, especially in low-porosity situations. Furthermore, we derive a new AVO approximate equation to relate the reflectivity series with the fluid modulus. The new AVO approximate equation allows us to estimate the fluid modulus of the reservoir in a more direct and stable manner than previous formulations. We also test the accuracy of the AVO approximate equation compared with the exact Zoeppritz equation and the Aki-Richards approximate equation. Eventually, we invert the fluid modulus with the linearized Bayesian methods for a more stable and accurate solution. We conduct synthetic and field data case studies to illustrate the feasibility and accuracy of this method in geofluid discrimination.


## 2. AVO Approximation Equation with Fluid Modulus

The prestack seismic response is directly influenced by the P-wave velocity, S-wave velocity, density, and incident angle based on the exact Zoeppritz equation. The P-wave velocity and S-wave velocity can be expressed as(1)VP=Ks+3/4μsρ,VS=μsρ,where Ks, μs, and ρ are the bulk modulus, shear modulus, and density of the fluid-saturated rock, respectively. From equation (1), we see that the fluid effects are mainly contained in the density and the bulk modulus of the fluid-saturated rock. Han and Batzle [27] have demonstrated that the bulk modulus of a fluid-saturated rock is more sensitive to the fluid-saturation effect than the P-wave velocity. The bulk modulus of a fluid-saturated rock containing the fluid effects can be expressed as
(2)Ks=Kd+f,(3)f=1−Kd/K02ϕ/Kf+1−ϕ/K0−Kd/K02,(4)μs=μd,following from Gassmann’s equation, where Kd, Ks, K0, and Kf are the bulk modulus of the dry rock, the saturated-rock frame, mineral grain, and geofluid, respectively, ϕ is the porosity of the rock, f is the increment resulting from the fluid saturation, and μs and μd are saturated and dry rock shear moduli. These equations indicate that the fluid in pores will affect bulk modulus but not shear modulus. Russell et al. [21] and Zong et al. [22] linearized the poroelasticity to estimate the fluid properties directly from prestack seismic records. The linearized AVO approximation with fluid factor can be expressed as
(5)RPPθ=14−γdry24γsat2sec2θΔff+γdry24γsat2sec2θ−2γsat2sin2θΔμsμs+12−sec2θ4Δρρ,where θ is the average of incident and transmission angles and γdry and γsat are the ratio of P-wave velocity to S-wave velocity of the saturated rock and dry rock, respectively. From equation (3), we see that the fluid factor f has the mixed effects of the fluid component and the rock matrix component. Therefore, the fluid factor f cannot reflect the fluid property directly.In this paper, we develop a novel P-wave reflectivity equation in terms of fluid modulus. The P-wave velocity, S-wave velocity, and density are expressed as a linearized function of the porosity for dry, clean sands [44]:
(6)VP=a+bϕ,VS=c+dϕ,ρ=k1−ϕ.The bulk modulus of the dry rock frame is expressed in terms of the P-wave velocity, S-wave velocity, and density as(7)Kd=ρVP2−43VS2.Substituting equation (6) into equation (7), we obtain a function of ϕ to express the bulk modulus. Since the modulus is the product of the density and the square of velocity, we obtain an equation that is cubic in terms of porosity. In terms of Knϕ=Kd/K0 [27], the bulk modulus is written as
(8)Kd=1+Aϕ+Bϕ2+Cϕ3K0.In low-porosity situations, equation (8) is further simplified as
(9)Kd=1−Dϕ2K0,where D is an empirical parameter. Based on the derivation of equation (9), D is determined by the statistical parameters in equation (6). Further, Han and Batzle [27] illustrated that D ranges from 1.45 to slightly more than 2.0 for consolidated rocks at high differential-pressure conditions with experience, and D is associated with the lithology. In this paper, the lithology of the target stratum is sandstone; we set D as 2 in terms of the experience results provided by Han and Batzle [27]. Substituting equation (9) into the linearized equation of the fluid modulus and the mixed pore/fluid term proposed by Han and Batzle [27],
(10)f=1−Knϕ2ϕKf,the new linearized relationship of the fluid modulus and the mixed pore/fluid term is written as
(11)f=D2ϕ2−Dϕ2Kf.In practice, the rock physics model parameterD should be tested against the well log data and calibrated in the work area.Compared with equation (10), the expression in equation (11) is more stable even in low-porosity situations. For simplification, we rewrite D2ϕ2−Dϕ2 as Y relating to the lithology and porosity, then we express equation (11) as
(12)f=YKf.We incorporate equations (5) and (12) to derive a new AVO approximation to estimate the fluid modulus directly. From equation (12), the fluid/pore term reflectivity can be divided into the fluid modulus reflectivity and the porosity-related reflectivity:
(13)Δff=ΔKfKf+ΔYY.Substituting equation (13) into equation (5) and merging the shear modulus with ΔY/Y, we define the combined new parameter as a shear modulus related term without physical meaning. The combination of these two parameters leads to the simplification of the new P-wave reflectivity and the decoupling of the fluid modulus from the mixed fluid/pore term. The new shear modulus related term of the rock matrix can be written as
(14)μm=Yμs.The form of the new shear modulus related term is the same as the relationship between the fluid modulus and the mixed fluid/pore term shown in equation (12). It should be noted that the new variable μm has no physical meaning which is proposed here only for the simplification of the expression. Substituting equation (13) into equation (5), we write the new AVO approximation as
(15)RPPθ=sec2θ4−γdry24γsat2sec2θΔKfKf+sec2θ4−γdry22γsat2sec2θ+2γsat2sin2θΔYY+γdry24γsat2sec2θ−2γsat2sin2θΔμmμm+12−sec2θ4Δρρ.Ifγdry is equal to γsat, which means that there is no fluid in the reservoir, the coefficient of the fluid modulus goes to zero, and the coefficient of the second term tends to be the opposite of that of the third term. In equation (15), the ratio of the P-wave velocity to the S-wave velocity in dry rock can be estimated from the laboratory measurements and real calculation [19]. The range of γdry2 and the range of their equivalent elastic constant ratios are shown in Table 1 provided by Russell et al. [21]. We see from Table 1 that γdry2 ranges from 4 to 1.333. However, there is no physical meaning when γdry2 equals 2 or 1.333; it means that the Poisson’s ratio is 0 or -1, respectively. Assigning 2.333 to γdry2 is appropriate for clean unconsolidated sandstones in most cases. However, it is also necessary to obtain the value of γdry2 from laboratory measurements due to its dependency on the reservoir type.Table 1
The ranges ofVP/VSdry2, VP/VSdry, Poisson’s ratio νdry, Kdry/μ, and λdry/μ of dry rocks [21].
VP/VSdry2
VP/VSdry
νdry
Kdry/μ
λdry/μ
4.000
2.000
0.333
2.667
2.000
3.333
1.826
0.286
2.000
1.333
3.000
1.732
0.250
1.667
1.000
2.500
1.581
0.167
1.167
0.500
2.333
1.528
0.125
1.000
0.333
2.250
1.500
0.100
0.917
0.250
2.233
1.494
0.095
0.900
0.233
2.000
1.414
0.000
0.667
0.000
1.333
1.155
−1.000
0.000
−0.667The weighting coefficients of the four terms ofΔKf/Kf, ΔY/Y, Δμm/μm, and Δρ/ρ in equation (15) varying with incident angles are shown in Figure 1. The dry rock velocity ratios γdry2 are 1.333 (red solid line), 2.000 (blue dotted line), and 2.333 (green dashed line), and the saturated rock velocity ratio is assigned as 4. The incident angle ranges from 0° to 60°. Figure 1(a) shows the weighting coefficient of the term of ΔKf/Kf, from which we see that it becomes larger with the increase of incident angle and decreases as the dry rock velocity ratio increases. Figure 1(b) shows the weighting coefficient of the term of ΔY/Y, whose variation trend is the same as that of ΔKf/Kf. Figure 1(c) shows the coefficient of the term of Δμm/μm, which jitters around zero when the dry rock velocity ratio is relatively small and moves away from zero with the increase of γdry2. Figure 1(d) shows the weighting coefficient of the term of Δρ/ρ; as expected, the three lines of different γdry2 are overlapped since they are independent of γdry2.Figure 1
The relationships of the coefficients of the four terms in equation (15) with the incident angle ranging from 0° to 60° for three dry rock velocity ratios of γdry2=1.333 (red solid line), γdry2=2.0 (blue dotted line), and γdry2=2.333 (green dashed line). The coefficients of the term of (a) ΔKf/Kf, (b) ΔY/Y, (c) Δμm/μm, and (d) Δρ/ρ.
(a)
(b)
(c)
(d)To test the accuracy of the new AVO approximation, we compare it with the exact Zoeppritz equation and the Aki-Richards approximate equation utilizing four models consisting of two sandstone layers saturated with gas or water. Their parameters are displayed in Table2. There are three AVO curves plotted in Figure 2 at the incident angle ranging between 0° and 45° in terms of the Zoeppritz equation, the Aki-Richards approximate equation, and the fluid modulus AVO approximation. The P- and S-wave velocities are derived in terms of the Biot-Gassmann poroelastic theory. In Model 1, the upper layer is the fully gas-saturated sandstone and the underlying layer is the fully water-saturated sandstone. From Figure 2(a), we see that the reflectivity amplitudes varying with the incident angles corresponding to Model 1 are almost the same as the Zoeppritz equation (green solid line), the Aki-Richards equation (red cross), and the fluid modulus equation (blue circle). In Model 2, there are two layers of fully water-saturated sandstone with the same fluid modulus, but the bulk moduli of the dry rock, shear modulus, and density are different. The AVO curves corresponding to the Zoeppritz equation, the Aki-Richards equation, and the fluid modulus AVO approximation for Model 2 are shown in Figure 2(b), from which it is seen that the fluid modulus equation is closer to the Zoeppritz equation compared with the Aki-Richards equation. It means that the fluid modulus AVO approximation is more accurate. In Model 3, the fully gas-saturated sandstone layer is overlapped by the fully water-saturated sandstone layer. The trend and amplitude of reflectivity of all three different equations are almost the same as those shown in Figure 2(c). In Model 4, the two layers are all sandstone fully saturated with gas; the differences between them are the bulk modulus of the dry rocks, shear modulus, and density. From Figure 2(d), we see that the Aki-Richards approximation equation is closer to the Zoeppritz equation; the error between the fluid modulus equation and the Zoeppritz equation becomes larger with the increase of the incident angle. Considering the assumption of the convolutional model, the accuracy of the fluid modulus equation at small incident angles is satisfied [45]. Therefore, the fluid modulus AVO inversion can be applied for fluid discrimination effectively. Generally speaking, the fluid modulus AVO approximation equation can provide a relatively accurate relationship between the fluid modulus and the reflectivity at different incident angles.Table 2
Four sandstone models used for the forward modeling.
Sand type
Kf (GPa)
Kd (GPa)
Ks (GPa)
μ (GPa)
ρ (kg/cm3)
ϕ
Model 1
Gas
0.10
5.00
40.0
4.00
1.99
25%
Water
2.38
5.00
40.0
4.00
2.26
25%
Model 2
Water
2.38
5.00
40.0
4.00
2.26
25%
Water
2.38
10.0
40.0
10.0
2.34
20%
Model 3
Water
2.38
10.0
40.0
10.0
2.26
20%
Gas
0.10
10.0
40.0
10.0
2.12
25%
Model 4
Gas
0.10
10.0
40.0
10.0
2.12
20%
Gas
0.10
5.0
40.0
4.00
1.99
25%Figure 2
AVO curves corresponding to the Zoeppritz equation (green solid line), the Aki-Richards approximate equation (red cross), and our fluid modulus approximation (blue circle) at the incident angle ranging from 0° to 45°: (a) Model 1, (b) Model 2, (c) Model 3, and (d) Model 4 listed in Table2.
(a)
(b)
(c)
(d)


## 3. Bayesian Inversion for Fluid Modulus

We can estimate the fluid modulus directly from prestack seismic records and utilize it to discriminate the geofluid filled in reservoirs. The forward modeling can be constructed by convoluting a wavelet with the fluid modulus AVO approximation shown in equation (15), then the inversion can be conducted based on the Bayesian inference. Equation (15) can be simplified as
(16)RPPθ=AR1+BR2+CR3+DR4,where
(17)R1=ΔKfKf,R2=ΔYY,R3=Δμmμm,R4=Δρρ,A=sec2θ4−γdry24γsat2sec2θ,B=sec2θ4−γdry22γsat2sec2θ+2γsat2sin2θ,C=γdry24γsat2sec2θ−2γsat2sin2θ,D=12−sec2θ4.We expand equation (16) to a matrix form for multiple incident angles:
(18)RPPθ1RPPθ2⋮RPPθna=A1B1C1D1A2B2C2D2⋮⋮⋮⋮AnaBnaCnaDnaR1R2R3R4,where θii=1,2,⋯,na represents the ith incident angle, Ai,Bi,Ci,Dii=1,2,⋯,na stand for the weighting coefficients of the ith incident angle, and na represents the number of incident angles. The single sample point is expanded to a single trace; the matrix equation is expressed as
(19)RPPθ1RPPθ2⋮RPPθna=A1B1C1D1A2B2C2D2⋮⋮⋮⋮AnaBnaCnaDnaR1R2R3R4,where the vector RPPθi represents the reflectivity of the ith incident angle at all sample points; i.e., RPPθi=RPP1θi,RPP2θi,⋯,RPPntθiT, nt is the number of time sampling points of a single trace, Rj=Rj1,Rj2,⋯,RjntTj=1,2,3,4 represents the jth reflectivity vector corresponding to Rj in equation (16), and Ai, Bi, Ci, and Di represent the nt×nt diagonal matrixes consisting of the weighting coefficients of A, B, C, and D in equation (18). We simplify equation (19) as
(20)R=Gm,where R stands for the left-hand side of equation (19) and G and m for the matrix and vector on the right-hand side of equation (19), respectively.Considering the wavelet convolution matrix and noise, the forward modeling formula is expressed as(21)d=Lm+n,where the coefficient matrix L can be expressed as L=WG, W is the wavelet matrix, n stands for random noise, and d represents the seismic records of different incident angles. Here, we assume the wavelet is known and fixed; i.e., there is no amplitude attenuation and phase distortion of the wavelet. The wavelet is extracted from the well log data and the seismic record via the well-seismic calibration, and the seismic record for inversion should be processed with attenuation compensation and migration before the inversion is performed.We conduct the fluid modulus prestack seismic inversion based on the Bayesian inference framework proposed by Ulrych et al. [46]. In terms of the Bayesian theory, the posterior probability distribution is expressed as
(22)Pmd=PdmPmPd∝PdmPm,where Pmd is the posterior distribution, Pd is the marginal probability distribution, and Pm is the prior distribution of four AVO parameter reflectivities. In this paper, the prior distribution of the parameters to be inverted is assumed to be a multivariable Cauchy probability distribution. On the assumption that noise in seismic records obeys Gaussian distribution, the likelihood function in equation (22) is written as
(23)Pdm=P0exp−2σn2−1d−GmTd−Gm,where P0=1/2πnant/2σnnant and σn is the standard deviation of the errors. The multivariable Cauchy probability distribution Pm is given by
(24)Pm=P1exp−2∑i=1ntln1+mTΦim,where P1=1/π2ntΨnt/2, Φi=EiTΨ−1Ei, Ψ is a 4×4 scale matrix representing the correlation of the four AVO parameter reflectivities, and Ei is a 4×4nt matrix to determine the location of the AVO parameter reflectivities for the single-trace calculation:
(25)Ejki=1,ifj=1andk=i,1,ifj=2andk=i+nt,1,ifj=3andk=i+2nt,1,ifj=4andk=i+3nt,0,otherwise,where Ei can extract the four AVO parameter reflectivities from the ith sample point. Substitution of equations (23) and (24) into equation (22) yields the posterior distribution:
(26)Pmd=P0P1exp−2σn2−1Gm−dTGm−d+2∑i=1ntln1+mTΦim.After some algebraic operations, the posterior probabilityPmd shown in equation (26) for a maximal value leads to
(27)GTG+αQm=GTd,where α∝σn2 is a hyperparameter that is determined by trial and error and Q is a 4nt×4nt nondiagonal matrix with elements of
(28)Qkn=∑i=1nt2Φkni1+mTΦim,k,n=1,2,⋯,4nt.The problem to be solved shown in equation (27) is nonlinear. In this paper, we solve it with the iterative reweighted least-squares algorithm [47, 48]. Meanwhile, considering the correlation between the four inverted parameters, we utilize the preconditioned conjugate gradient method proposed by Zong and Yin [49] in the inversion method to decouple the inverted parameters completely. Subsequent to the reconstruction of the reflectivity of different AVO parameters, the rock properties of Kf, Y, μm, and ρ can be obtained by solving Fx=m, where x represents the logarithm of four inverted parameters and F denotes the first-order difference matrix. We solve it with the conjugate gradient method with the initial model of four parameters. In the synthetic and field data tests, the initial model and true model of density and porosity are obtained from the well log, and the other three parameters are obtained with well log data and rock physics theory. The well log can provide the porosity, P-wave velocity, S-wave velocity, and density. Based on the interpretation results of the well log, we can determine the empirical parameters and utilize Gassmann’s equation expressed in equations (1), (2), and (11) to calculate the Kf, Y, and μm. In addition, the well log data are converted to the time domain via the time-depth conversion and well-seismic tie with the depth, P-wave velocity, and seismic record in both the synthetic and field data tests.


## 4. Synthetic Data Test

We test the ability of our formula and method with a synthetic angle gather with incident angles ranging from 0° to 30° with a gap of 5°. The model parameters are obtained from a section of the well log data with low-pass filtering. Figure3 shows the model parameters of Kf, Y, μm, and density which are taken from the well log data and Gassmann’s equation shown in equations (1), (2), and (11) [25, 50]. Considering that the lithology of the target stratum is sandstone, we set the empirical parameters D and Kd to 2 and 37 GPa. The synthetic angle gather is generated by convoluting a Ricker wavelet of 45 Hz dominant frequency with the reflectivity series derived from equation (15). Figure 4(a) shows the noise-free synthetic angle gather and Figure 4(b) shows the synthetic angle gather contaminated with Gaussian random noise; the signal-to-noise (S/N) ratio is 5. In this paper, S/N represents the ratio of the root-mean-square amplitude of the signal to that of the noise. The four AVO parameters estimated from the noise-free angle gather are shown in Figure 5, where the blue solid line indicates the model parameters from the well log data and the red dashed line indicates the AVO parameters estimated from the synthetic angle gather. From Figure 5, we see that the inverted fluid bulk modulus, porosity-related term Y, and new shear modulus related term are close to the true values with small deviations between 645 ms and 660 ms. The deviations of the inverted density are a little bit obvious compared with the first three AVO parameters. Downton [40] pointed out that the inverted density is sensitive to the maximal offset, and the estimated density reflectivity showed more evident bias than the P-impedance and S-impedance reflectivities when the data are contaminated. Therefore, the instability of the inverted density will be magnified when the angle gathers are contaminated with noise. Figure 6 shows the inverted AVO parameters from the noisy angle gather with S/N=5. The blue solid lines represent the true AVO parameters of well log data and the red dashed lines represent the inverted AVO parameters. There is overall accuracy degradation of the estimated results; the deviations of the fluid bulk modulus, porosity-related term Y, and new shear modulus related term are acceptable. The influence of the noise on the density shows more serious misfit.Figure 3
The model parameters ofKf, Y, μm, and density of the real well log data used to test the inversion accuracy of the proposed formula and method.Figure 4
The synthetic angle gathers with incident angles ranging from 0° to 30° obtained using the true model parameters: (a) noise-free data; (b) noisy data withS/N=5.
(a)
(b)Figure 5
The estimated AVO parameters ofKf, Y, μm, and density (dashed line) from the synthetic noise-free angle gather and their true values (solid line).Figure 6
The estimated AVO parameters ofKf, Y, μm, and density (dashed line) from the noisy synthetic angle gather with S/N=5 and their true values (solid line).The synthetic data test demonstrates the feasibility and the stability of the proposed method to estimate effectively the fluid bulk modulus, porosity-related termY, and new shear modulus related term. Similarly to conventional prestack seismic inversion, the limitation of our method is the accuracy of the inverted density, which has been analyzed by Downton [40]. Further, we compare the inversion results with the results from the method proposed by Zong et al. [26] to illustrate the performance of our method in the low-porosity situation. Figure 7 shows the fluid modulus inversion results from the noise-free angle gather. The left subplot displays the variation of porosity with time, the middle subplot shows the inverted fluid modulus with our method and the true fluid modulus, and the right subplot shows the inverted fluid modulus with the method of Zong et al. [26] and the true fluid modulus. In both Figures 7(b) and 7(c), the blue solid lines denote the fluid modulus of the true model and the red dashed lines stand for the inverted results. Figure 8 shows the fluid modulus inversion results from the noisy angle gather with S/N=5. It can be seen that the inversion results from the method of Zong et al. [26] are more unstable in the relatively low-porosity locations such as the layers shallower than 650 ms in Figures 7 and 8. However, the inversion results with our method show better stability and accuracy.Figure 7
The comparison between the inverted fluid modulus with our method and the method of Zong et al. [26] from the noise-free angle gather. The left subplot (a) displays the variation of porosity with time, the middle subplot (b) shows the inverted fluid modulus with our method (red dashed line) and the true fluid modulus (blue solid line), and the right subplot (c) shows the inverted fluid modulus with the method of Zong et al. [26] (red dashed line) and the true fluid modulus (blue solid line).
(a)
(b)
(c)Figure 8
The comparison between the inverted fluid modulus with our method and Zong et al. [26] from the noisy angle gather with S/N=5. The left subplot (a) displays the variation of porosity with time, the middle subplot (b) shows the inverted fluid modulus with our method (red dashed line) and the true fluid modulus (blue solid line), and the right subplot (c) shows the inverted fluid modulus with the method of Zong et al. [26] (red dashed line) and the true fluid modulus (blue solid line).
(a)
(b)
(c)


## 5. Field Data Test

We apply the proposed fluid bulk modulus AVO approximation to prestack seismic inversion for fluid discrimination. Strata’s demo data set is applied for fluid discrimination. Figure9(a)–9(d) shows the angle gathers ranging from 5° to 25° at CDP 176, CDP 528, CDP 880, and CDP 1232, respectively. The gathers are in the time-angle domain; the maximum incident angle is around 29°. Figure 10(a)–10(c) illustrates the partially stacked angle gathers with small, middle, and large incident angles ranging from 1° to 10°, 10° to 19°, and 19° to 28°, respectively. Well log data are used to establish the low-frequency (0-0-10-15 Hz) initial P-wave velocity, S-wave velocity, and density models as shown in Figure 11. To estimate the four AVO parameters proposed in this paper, we combine the rock physics and the well log data. The well log can provide the P-wave velocity, S-wave velocity, density, and porosity. Based on the interpretation results of the well log, the empirical parameter of the statistical function can be determined and used for Gassmann’s equation to calculate the four inverted parameters. The initial models of fluid bulk modulus (Kf), the porosity-related term (Y), and the new shear modulus related term (μm) are shown in Figure 12. For this field data test, considering that the lithology of the target stratum is sandstone, we determine the parameter D as 2 and the bulk modulus of the dry rock as 37 GPa. Meanwhile, the main purpose of our method is to estimate the fluid modulus for discrimination. From equation (15), we see that the derivation of the empirical parameter mainly influences the accuracy of the second and third inverted parameters. Based on the AVO approximation in equation (15), the new four AVO parameters estimated from the prestack data are shown in Figure 13. The black rectangle on the Kf inversion profile indicates the target sandstone saturated with gas. The gas-saturated sandstone can be identified clearly with a low fluid bulk modulus of about 0.1 GPa. However, the traditional inversion results of P-wave velocity and S-wave velocity shown in Figure 14 cannot illustrate the gas-saturated sandstone obviously compared with the inverted fluid modulus due to the mixed effects of the rock matrix and fluid properties on them.Figure 9
The angle gathers with incident angles ranging from 5° to 25° at (a) CDP 176, (b) CDP 528, (c) CDP 880, and (d) CDP 1232.
(a)
(b)
(c)
(d)Figure 10
Partially stacked angle gathers of (a) small-angle range, (b) middle-angle range, and (c) large-angle range.
(a)
(b)
(c)Figure 11
The initial model of (a) P-wave velocity, (b) S-wave velocity, and (c) density from the well log data.
(a)
(b)
(c)Figure 12
The initial model of (a)Kf, (b) Y, (c) μm, and (d) density from the initial model of P-wave velocity, S-wave velocity, and density shown in Figure 11 in terms of the rock-physical model.
(a)
(b)
(c)
(d)Figure 13
The inversion results of (a)Kf, (b) Y, (c) μm, and (d) density. The black rectangle in (a) represents the target stratum.
(a)
(b)
(c)
(d)Figure 14
The traditional inversion results of (a)VP, (b) VS, and (c) density.
(a)
(b)
(c)For further clarification about the ability of our method to discriminate the geofluid, we also plot the well log at CDP 495 in Figure15, where the location of the gas-saturated sandstone is allocated by the black rectangle. We utilize the resistivity from the well log exhibiting high values in the gas-saturated sandstone to illustrate the accuracy of our method. From Figure 15(a), we see that there are high values of resistivity from 638 ms to 654 ms. The corresponding fluid bulk moduli (Figure 15(b)) has low values indicating the layers saturated with gas. The boundary of the gas-saturated sandstone layer may result in the value of the fluid bulk moduli being relatively larger. As a mixed term of fluid bulk modulus and porosity, the P-wave velocity (Figure 15(c)) cannot obviously indicate the fluid types compared with the fluid modulus (Figure 15(b)). Meanwhile, there is no obvious correlation between the S-wave velocity curve (Figure 15(d)) and the resistivity curve (Figure 15(a)). In addition, the comparison between the inverted model and well log data is shown in Figure 16 to illustrate the accuracy of the inversion results. For the log data, the density can be obtained from the well log directly; the other three inverted parameters Kf, Y, and μm are derived based on the parameters from the well log and the rock physics theory. Due to the high-frequency component of the well log data, we apply a low-pass filter with 70-80 Hz high-cut frequency to the well log to compare with the inverted model. The blue lines represent the inverted model, and the red lines denote the low-pass filtered well log data. From Figure 16, we see that the inverted model matches well with the well log data. The correlation coefficient between the inverted model and the well log data is 0.8924, 0.8567, 0.8824, and 0.8648. Though there are also some obvious mismatches of the inverted Y, μm and density, the main purpose of our method is to estimate the fluid modulus for fluid discrimination; the accuracy of the inverted Kf is acceptable.Figure 15
(a) The resistivity from the well log, and the inversion results of (b)Kf, (c) P-wave velocity, and (d) S-wave velocity. The black rectangle represents the location of the target stratum.Figure 16
The comparison between the low-pass filtered parameters and the inverted parameters of (a)Kf, (b) Y, (c) μm, and (d) density. The blue lines represent the inverted model and the red lines denote the low-pass filtered well log data with a high-cut frequency of 70-80 Hz. The correlation coefficient between the inverted model and the well log data is 0.8924, 0.8567, 0.8824, and 0.8648.
(a)
(b)
(c)
(d)


## 6. Discussion

The challenge of the fluid discrimination is the uncertainty of the fluid indicators associated with the definition of the fluid indicators and the instability of the traditional elastic parameters such as the P-wave velocity and the S-wave velocity. Recently, the fluid bulk modulus has been demonstrated to be an effective indicator to discriminate the fluid types. Meanwhile, the fluid bulk modulus has been decoupled from the traditional Aki-Richards equation [23, 26]. However, there are so many approximations to the empirical formula in Yin and Zhang [23] resulting in a larger mismatch with the Zoeppritz equation. Zong et al. [26] proposed a new AVO approximation equation containing the porosity in the denominator of the inversion term, which leads to the instability of the inversion results in low-porosity situations. The instability of the porosity-related term and density will aggravate the accuracy of the inversion results. In this paper, we propose a new AVO approximation in a more stable form based on the empirical model and Gassmann’s equations. The fluid bulk modulus can be obtained accurately. The main disadvantage of our method is the need for abundant rock-physical parameters. It may not be necessary to provide sufficient rock-physical parameters, which can be replaced by the empirical formula. The deviation of the empirical factor will cause the magnitude of the inversion results to not reflect the accurate situations. Therefore, sufficient rock-physical parameters for an accurate empirical factor are necessary to obtain a more reliable inversion result. In addition, the linearized relationship in equation (9) may limit its application in different work areas. The generalization of other empirical relations introduced into our AVO approximation needs further studies.


## 7. Conclusions

Fluid discrimination plays an important role in reservoir prediction. The fluid bulk modulus is sensitive to the layers saturated with geofluid. In this paper, we derive a new AVO approximation equation relating the fluid bulk modulus with the reflectivity series. Initially, we illustrate the accuracy of the fluid modulus AVO approximation with the model test. The forward modeling can be constructed by convoluting the fluid modulus AVO approximation with a seismic wavelet. Then, we illustrate the accuracy and stability of the fluid modulus inversion with synthetic noise-free and noisy seismic data. The model is constructed according to the well log data. Eventually, we test the fluid bulk modulus prestack inversion on Strata’s demo seismic record where the target reservoirs are sandstone layers saturated with gas. The availability to discriminate geofluid can be demonstrated in terms of the well log data. The magnitude of the inversion results may not reflect accurate situations due to the deviation from the empirical formula and parameters. However, the inversion results are acceptable for geofluid discrimination by utilizing the variation of the inverted fluid bulk modulus and the prior information.



---

*Source: 2656747-2019-11-26.xml*