ABSTRACT

Linearized approximations to the P-wave reflectivity as a function of the incidence angle (called amplitude variation with offset) involve the extraction of band-limited reflectivity terms that are a function of changes in the elastic constants of the earth across each lithologic interface. The most common of these extracted reflectivities are the intercept and gradient, usually labeled A and B, respectively. The extended elastic impedance (EEI) method uses a rotation angle χ to map A and B into a new reflectivity corresponding to a particular elastic parameter. The success of EEI depends on finding an optimum value for the angle χ. This value is usually calculated by correlating the EEI result over a range of χ angles with various elastic parameters and then finding the best correlation coefficient. We have developed a new approach for the interpretation of the EEI method, which incorporates the Biot-Gassmann poroelastic theory and attaches a physical meaning to the χ angle. We call this method extended poroelastic impedance (EPI). The main advantage of the EPI method is that the χ angle is now interpreted as a parameter that is dependent on the dry-rock properties of the reservoir, rather than a parameter whose value is estimated empirically. The method is evaluated by numerical and synthetic seismic examples and by application to field data from a gas sand reservoir.

INTRODUCTION

Over the past 30 years, numerous techniques have been developed to estimate fluid and lithology effects in the earth’s subsurface using prestack seismic data. These include the amplitude variation with offset (AVO) method (Gidlow et al., 1992), prestack simultaneous inversion for P-impedance, S-impedance, and density (Hampson et al., 2005), the lambda-mu-rho (LMR) method (Goodway et al., 1997), and the elastic impedance (EI) (Connolly, 1999) and extended EI (EEI) methods (Whitcombe et al., 2002). In addition to the deterministic methods just described, much work has been done on stochastic and probabilistic approaches to inversion (Buland and Omre, 2003; Spikes et al., 2008; Alemie and Sacchi, 2011; Grana, 2014; Connolly and Hughes, 2016), which assign reasonable error bounds to the estimates of rock properties. Recent work also extends the geostatistical inversion approach to geologic facies and petrophysics models in conjunction with impedances (Lang and Grana, 2017; Fjeldstad and Grana, 2018). All of this work is related to the development of empirical rock-physics models, which is well-described in many publications (Avseth et al., 2005; Mavko et al., 2009; Dvorkin et al., 2014; Vernik, 2016). We will use several of these rock-physics models to develop a synthetic seismic data set to test our method.

Poroelasticity theory (Biot, 1941; Gassmann, 1951) was incorporated into the theory of prestack inversion (Hedlin, 2000; Russell et al., 2003), which led to a generalization of the LMR method. Gray et al. (1999) derive three-term AVO equations for LMR and bulk modulus, mu and rho, and these equations were generalized by Russell et al. (2011) into a three-term poroelastic AVO equation. This work has been extended to other parameters such as the Young’s modulus and Poisson’s ratio and by the use of a Bayesian probabilistic approach by Zong et al. (2012, 2013) and Yin and Zhang (2014). In this paper, the use of poroelasticity is extended to a generalization of the theory of EEI, called extended poroelastic impedance (EPI). Before discussing the theory of this new approach, the basic theory of prestack linearized amplitude analysis will be reviewed.

For a plane wave of incidence angle θi striking the boundary between two elastic layers of P-wave velocities VPi and VPi+1, S-wave velocities VSi and VSi+1, and densities ρi and ρi+1, the reflectivity as a function of the angle can be written as the linearized sum of three terms (Richards and Frasier, 1976; Aki and Richards, 1980) given by  
RPP(θ)=aΔVP2VP+bΔVS2VS+cΔρ2ρ,
(1)
where a=1/cos2θ, b=8ksatsin2θ, c=14ksatsin2θ, ΔVP=VPi+1VPi, VP=(VPi+1+VPi)/2, ΔVS=VSi+1VSi, VS=(VSi+1+VSi)/2, Δρ=ρi+1ρi, ρ=(ρi+1+ρi)/2, ksat=(VS/VP)sat2, and θ=(θi+1+θi)/2, where θi+1=sin1(sinθi(VPi+1/VPi)). In the term ksat, the subscript sat stands for the saturated, or in situ, velocities. In some formulations, the factor two in the denominator of the three terms ΔVP/2VP, ΔVS/2VS, and Δρ/2ρ is incorporated into the terms a, b, and c. However, writing these terms with a denominator of two indicates that these are elastic reflectivity terms found by dividing the difference between the elastic parameters below and above the interface by their sum. For example, the three reflectivity terms in equation 1 can be written as RVP=(ΔVP/2VP)=(VPi+1VPi)/(VPi+1+VPi), RVS=(ΔVS/2VS)=(VSi+1VSi)/(VSi+1+VSi), and Rρ=(Δρ/2ρ)=(ρi+1ρi)/(ρi+1+ρi), where the subscripts indicate which elastic parameter is being used in the calculation. We will use this subscript notation in subsequent forms of the reflectivity.
Equation 1 is often referred to as the Aki-Richards equation and can be rearranged as  
R(θ)=A+Bsin2θ+Csin2θtan2θ,
(2)
where A=(ΔVP/2VP)+(Δρ/2ρ), B=(ΔVP/2VP)8ksat(ΔVS/2VS)4ksat(Δρ/2ρ), and C=(ΔVP/2VP). In this formulation, considered to be the standard AVO, equation, A is called the intercept, B the gradient, and C the curvature. Wiggins et al. (1983) are the first to derive this form of the AVO equation, although their paper included only the A and B terms and had a sign error on the ΔVP/2VP term in B.

Note that A, B, and C can also be thought of as reflectivity terms. This is obvious for the C term, which is equal to RVP. However, by using the logarithmic differential approximation from calculus, we find that (ΔP/2P)(ΔlnP/2) for any parameter P. Using this formulation, A can be thought of as the acoustic impedance reflectivity, or A=RAI=(ΔAI/2AI)=(ΔVP/2VP)+(Δρ/2ρ), where AI=ρVP is the acoustic impedance. In addition, if we make the assumption that (VS/VP)sat0.5, or ksat=0.25, we find that BRAI2RSI=(ΔAI/2AI)2(ΔSI/2SI), where RSI=(ΔSI/2SI)=(ΔVS/2VS)+(Δρ/2ρ) is the shear impedance reflectivity.

Using this same logarithmic approximation, Connolly (1999) shows that the Aki-Richards equation (equation 1) could be derived from what he called EI, using the following approximation:  
REI(θ)=ΔEI(θ)2EI(θ)ΔlnEI(θ)2,
(3)
where  
EI(θ)=VPaVSbρc,
(4)
where a, b, and c are the scaling terms in the original Aki-Richards equation (equation 1) given above.
EEI (Whitcombe et al., 2002) is based on the rearrangement of the Aki-Richards equation given earlier in equation 2. To extend EI to EEI, Whitcombe et al. (2002) drop the third term in equation 2 and then transform sin2θ to tan χ to get  
R(θ)A+Btanχ,90o<χ<90o.
(5)
Finally, they divide by cosχ to get the final form of the EEI reflectivity, which we will call the scaled EEI reflectivity, or  
REEI(θ)=Acosχ+Bsinχ.
(6)
The impedance form of EEI is given by Whitcombe et al. (2002) by scaling equation 4 by reference values (which are computed by averaging over the zone of interest) and then substituting the χ angle formulation into the a, b, and c terms, or  
EEI(χ)=VP0ρ0[(VPVP0)p(VSVS0)q(ρρ0)r],
(7)
where p=cosχ+sinχ, q=8ksatsinχ, and r=cosχ4ksatsinχ. Whitcombe et al. (2002) caution that χ is a rotation angle that maps A and B into a new reflectivity corresponding to a particular elastic parameter rather than a physical angle such as the angle of incidence θ. The value of χ is found by rotating equation 7 through a range of χ angles, correlating each rotation with a chosen elastic parameter, and choosing the χ angle that gives the best correlation.

EPI THEORY

Using theory developed by Dong (1996), explained in detail in Appendix A, Whitcombe et al. (2002) develop explicit reflectivity terms for K (bulk modulus), λ (the first Lamé parameter), and μ (shear modulus), which they write as  
RK=ΔK2K=(A+B3+2f)(3+2f34ksat),
(8)
 
Rλ=Δλ2λ=(A+B2+f)(2+f24ksat),
(9)
 
Rμ=Δμ2μ=(ABf)(f4ksat),
(10)
where f=C/A. In each reflectivity, the term under B can be interpreted as the inverse of tanχ, and the second term in brackets can be interpreted as a scaling term. To find a value for f, Whitcombe et al. (2002) use Gardner’s equation (Gardner et al., 1974), given by  
ρ=aVP0.25.
(11)
Gardner’s equation is an empirical density-velocity relationship that was derived in the Gulf Coast and is intended for the description of clastic rocks. However, it has been found to give a reasonable fit to a wide range of lithologies. To explicitly derive a value for f, note first that differentiation of equation 11 leads to the differential form (Δρ/ΔVP)=(a/4)VP0.75. This then gives us the relationship (Δρ/ρ)=1/4(ΔVP/VP), which upon substitution into the intercept term A from equation 2 gives A=(5ΔVP/8VP). Using the expression for C in equation 2, this finally gives a value of f=0.8 and leads to χ angles of  
χK=tan1(13+1.6)=tan1(0.217)=12.3o
(12)
for the bulk modulus reflectivity,  
χλ=tan1(12+0.8)=tan1(0.357)=19.7o
(13)
for the first Lamé constant reflectivity, and  
χμ=tan1(10.8)=tan1(1.25)=51.3o
(14)
for the shear modulus reflectivity. As pointed out by Whitcombe et al. (2002), the first two angles give real values of θ, but the third value for the shear modulus does not. In the above three equations, the arguments inside the arctangent function have been written in two different ways, as a ratio of two numbers and as the result of doing the division. The choice between a ratio and a single value will have an effect on the final calculation. In the cases shown above, the single-valued argument was used to compute the angles that correspond exactly to the angles given by Whitcombe et al. (2002).
The methods described so far depend on elastic constants that were estimated from the in situ reservoir rocks saturated with the reservoir fluid. In the saturated case, the P- and S-wave velocities can be written as  
VPsat=Ksat+(4/3)μsatρsat=λsat+2μsatρsat
(15)
and  
VSsat=μsatρsat,
(16)
where Ksat is the saturated bulk modulus, λsat is the saturated first Lamé parameter, and μsat is the saturated shear modulus (or the first Lamé parameter).
Biot (1941) and Gassmann (1951) extend the concept of elasticity to poroelasticity, in which the effects of the porosity and fluid saturation of the reservoir rocks are separated from the dry rock frame. Two main results came out of this work, the first being that the shear modulus is independent of the saturation, or  
μsat=μdry,
(17)
and the second being that the dry bulk modulus and second Lamé parameter are related to their saturated terms by the addition of a term (which we will call p) that is related to the porosity and fluid saturation of the reservoir. We can write these two expressions as  
Ksat=Kdry+p
(18)
and  
λsat=λdry+p,
(19)
where Kdry is the dry bulk modulus, λdry is the dry second Lamé parameter, and p=α2M is a fluid/porosity term, where α is the Biot coefficient and M is the Biot fluid modulus. The Biot coefficient and fluid modulus that were introduced by Biot (1941) can be related to Gassmann (1951) formulation by the equations  
α=1KdryKm
(20)
and  
M=(ϕKf+αϕKm)1,
(21)
where Km is the mineral bulk modulus, Kf is the fluid bulk modulus, and ϕ is the porosity. This leads to the following equation for the P-wave velocity (equation 16 stays the same due to the relationship in equation 17):  
VPsat=Kdry+p+(4/3)μρsat=λdry+p+2μρsat.
(22)
Using equations 16 and 22, Russell et al. (2003) derive a generalized form of the λρ equation, written as  
pρsat=IPsat2dISsat2,
(23)
where IPsat = saturated P-wave impedance, or ρsatVPsat, ISsat = saturated S-wave impedance, or ρsatVSsat, and  
d=(VPVS)dry2=Kdryμ+43=λdryμ+2
(24)
is the dry-rock velocity ratio squared. Various values of the dry-rock velocity ratios are summarized in Table 1.
The Dong (1996) formulation can be used to combine the concept of poroelasticity with the EEI technique, which will be called EPI. The mathematical derivation is given in Appendix A, and it shows that we can write the generalized EPI equation using the dry rock P- to S-wave velocity ratio squared d as  
Rp=(A+d4+(4d)fB)(4+(4d)f44dksat).
(25)
To be more consistent with EEI, equation 25 can be reformulated using kdry, the inverse of d, giving  
Rp=(A+B4kdry+(4kdry1)f)(4kdry+(4kdry1)f4kdry4ksat),
(26)
where kdry=(1/d)=(VS/VP)dry2. We will now show that equations 25 and 26 are a generalization of equations 810.
The second term in brackets in equations 25 and 26 can be thought of as a scaling term, and the first term can be thought of as equivalent to equation 5, meaning that tanχ can be interpreted using the squared dry rock VP-to-VS ratio in equation 25 with the equation  
tanχ=sin2θ=d4+(4d)f,
(27)
or, using the inverse of d in equation 26 as  
tanχ=sin2θ=14kdry+(4kdry1)f.
(28)
Equations 27 and 28 give us a way to compute the angle χ using the squared dry rock velocity ratio and the ratio f=C/A rather than empirically finding this value by correlating well logs transformed using equation 7 with actual well-log values.

NUMERICAL RESULTS

Table 1 shows a range of values for the dry rock VP-to-VS ratio squared d, the inverse of d, or the VS-to-VP ratio squared kdry, the dry rock Poisson’s ratio νdry, the dry bulk modulus over shear modulus ratio Kdry/μ, and the dry first Lamé parameter over shear modulus ratio λ/μ. Using this table gives us a way to physically interpret equations 2328.

In Table 1, two values are of particular importance: d=4/3 and d=2. When d is equal to 4/3 (or kdry=0.75), this implies that the Kdry/μ ratio equals zero and the term p is equivalent to the bulk modulus Ksat. Substitution into equation 25 or 26 gives  
Rp(d=4/3)=(A+B3+2f)(3+2f34ksat),
(29)
which is identical to equation 8 or RK. When d is equal to 2 (or kdry=0.5), this implies that the λdry/μ ratio equals zero and p is equivalent to λsat. Substitution into equation 25 or 26 gives  
Rp(d=2)=(A+B2+f)(2+f24ksat),
(30)
which is identical to equation 9 or Rλ. Finally, consider the limiting value as d approaches (or, equivalently, as kdry approaches zero). Substitution of kdry=0 into equation 26 gives  
Rp(kdry=0)=(ABf)(f4ksat),
(31)
which is identical to equation 10 or Rμ. Thus, equations 25 and 26 reduce to the three reflectivities from Whitcombe et al. (2002).

In Table 1, the value of d=2.333 (or kdry=0.429), gives a Kdry/μ ratio of 1.0 and is often used to model clean sands, and the value of d=3 (or kdry=0.333) has been observed in deeper sediments offshore Brazil (Dillon et al., 2003). Finally, as discussed by Hedlin (2000), a value of d=2.233 (or kdry=0.448) corresponds to a Kdry/μ ratio of 0.9, and fits the experimental work of Murphy et al. (1993). Hedlin (2000) also suggests that this was a better value to use in equation 21 than 2 (the accepted LMR value) because a value of 2 implies a dry rock Poisson’s ratio of zero, as shown in Table 1.

To explore the full relationship between the dry rock velocity ratio and χ, note that Figure 1 shows a plot of χ angle (in °) versus kdry for three different values of the C/A ratio f: 0.5, 0.8, and 1.1. A value of f=0.8 is obtained from Gardner et al. (1974) relationship, and it is the value used by Whitcombe et al. (2002). By varying this value over this large range, it can be seen that the value of f has very little impact on the result.

In Figure 1, notice that kdry values of 0, 0.5, and 0.75 correspond to χ rotation angles of 51.3°, 12.4°, and 19.8°, respectively, which are the same values computed by Whitcombe et al. (2002) for the χμ, χK, and χλ reflectivities, respectively. These values can be computed from equations 25 or 26, as can the values of incident angle θ. Specifically, kdry values of 0.5 and 0.75 correspond to θ values of 28° and 37°, respectively, whereas a kdry value of zero does not give us a real value of θ. This partially explains the +90° to 90° apparent discontinuity that occurs between 12.4° and 51.3° in Figure 1. This apparent discontinuity appears at different angles for the three values of f, and it is the value at which the denominator in equation 26 equals zero. For f=0.8, it occurs when kdry=1/9. Also of note in Figure 1 and equation 26 is that the χ rotation angle decays asymptotically to zero as the value of kdry increases and that the curve is independent of f when kdry=0.25.

However, the discontinuity seen in Figure 1 is not a physical phenomenon but rather a mathematical artifact. The plot in Figure 1 has been constrained to lie within the angle range of +90° and 90°, so the values “wrap” around whenever they exceed either of those values. As was pointed out in the previous section, there are two ways of computing the arctangent, one that uses a single value z, where z=y/x, and the other that uses y and x separately. Figure 1 was created using the tan1(z) function (written atan(z) in computer languages), and it gives the values that were computed in the previous section for the three reflectivities in Whitcombe et al. (2002). A more correct calculation of the arctangent that is found in most computer languages uses the itan2(x,y) function, which uses all four quadrants of the x and y plane in the calculation. Figure 2 shows the result of using the atan2 function, in which the curve is now “unwrapped.” Note that χμ=131.3°, which is equal to 51.3°+180°, and that the curves now monotonically decrease between kdry values of zero and one.

As a summary of these results, Table 2 shows a subset of the dry rock d=(VP/VS)2 and kdry=(VS/VP)2 values shown in Table 1, but now with the χ angle included in the right column. On this table, it should be noted that a χ angle of 22.4°, corresponding to a d value of 2.233 (or a kdry value of 0.448) ties together work done in a rock-physics laboratory (Murphy et al., 1993) to the theory of EEI (Whitcombe et al., 2002) and to the Biot-Gassmann theory.

SYNTHETIC EXAMPLE

To test the EPI method on a synthetic seismic example, we created a four-layer well log of alternating shale and sandstone layers. For each layer, we built a petroelastic model with realistic petrophysical parameters. Because our model for the shale was unique, we will describe it in detail here. The input variables for the shale model were volume of shale (Vsh) and total porosity (ϕt). The effective bulk and shear moduli of the shale were computed using linear fits to ϕt, which were developed from a North Sea field, and given by  
Keff=Km37ϕt
(32)
and  
μeff=μm19ϕt,
(33)
where the mineral moduli were computed from a Reuss average of the shale and quartz constituents given by  
Km=(VshKsh+1VshKqz)1
(34)
and  
μm=(Vshμsh+1Vshμqz)1,
(35)
with Kqz=37  GPa, μqz=45, Ksh=21  GPa, and μsh=5.3. The bulk density was computed from a Voigt average of the brine and mineral components given by  
ρb=ρwϕt+ρm(1ϕt),
(36)
where  
ρm=ρshVsh+ρqz(1Vsh),
(37)
with ρqz=2.65  g/cm3 and ρsh=2.650  g/cm3. The brine densities were computed from the Batzle-Wang equations (Batzle and Wang, 1992) using a pore pressure of 20 MPa, temperature of 35°C, and salinity of 0.01. The compressional and shear velocities VP and VS of the shale were then computed from the effective bulk and shear moduli and densities by  
VP=Keff+43μeffρb
(38)
and  
VS=μeffρb.
(39)
The elastic parameters for the dry rock frame of the sandstone layers were computed using the high-porosity sandstone model of Dvorkin and Nur (1996), which combines the Hertz-Mindlin theory with the modified lower Hashin-Shtrikman bound. Because the theory is fully described by Dvorkin and Nur (1996), we will not rewrite it here. However, there are several parameters in this model that were not described above, the effective pressure (Peff), coordination number (C), and critical porosity (ϕ0). For this model, we used Peff=25  MPa, C=9, and ϕ0=0.4. The mineral modulii and bulk-density values were computed as in equations 3437, with the values of the shale and quartz components as given earlier. The Batzle-Wang equations were used to calculate the properties of the gas and brine fluid components. The effective bulk modulus of the porous, fluid-filled sandstone was computed from the dry-rock bulk modulus using the Gassmann equation, which can be found by substituting equations 20 and 21 into equation 18. The compressional and shear velocities VP and VS of the shale are then computed using equations 15 and 16. The input parameters for our four layers are shown in Table 3, in which the shales both have identical parameters, the top sand consists of 100% gas and has no shale, and the bottom sand consists of 100% brine and has 20% shale.

The resulting synthetic well-log curves are shown in Figure 3. In both shales, the P-wave velocity is 3335  m/s, the S-wave velocity is 1460  m/s (giving a VP/VS ratio of 2.28 and a Poisson’s ratio of 0.38), and the density is 2.28  g/cm3. In the gas sand, the P-wave velocity is 2982  m/s, the S-wave velocity is 2022  m/s (giving a VP/VS ratio of 1.48 and a Poisson’s ratio of 0.08), and the density is 2.03  g/cm3. In the basal wet sand, the P-wave velocity is 2679  m/s, the S-wave velocity is 1431  m/s (giving a VP/VS ratio of 1.87 and a Poisson’s ratio of 0.3), and the density is 2.15  g/cm3. The low P-impedance in the wet sandstone layer suggests a gas sand, but this low impedance is caused by the presence of 20% shale (recall that in our petroelastic model, the bulk modulus of shale was much lower than that of sandstone), and it is clear from the VP/VS ratio of this sand that it is not hydrocarbon bearing.

Figure 4 shows the synthetic angle gather computed from the elastic logs shown in Figure 3, computed using a far angle of 45° (in Figure 4 and the subsequent seismic displays in our synthetic study, the amplitude colorbar has been normalized to standard deviations away from the mean, with red for positive amplitude and blue for negative amplitude, where the brightest colors show one standard deviation away from the mean). This synthetic was generated using the full Zoeppritz equations at each layer interface and a 30 Hz Ricker wavelet. Notice the strong increasing AVO response at the top and base of the gas sand, indicative of a class 3 AVO anomaly (Rutherford and Williams, 1989), and the lack of an amplitude response at the top of the basal wet sand. This is quantified in Figure 5, which shows the picked amplitudes from the angle gather shown in Figure 4 for the three interfaces.

Figure 6 shows a 50-trace synthetic anticlinal model that was created by inserting the well shown in Figure 3 at common depth point (CDP) 25, replicating the four layers at each CDP location and stretching, squeezing, and truncating the log values to conform to the picked horizons for the anticline. The top event defines the top of the anticline, and the other two events were kept flat and truncated at the anticlinal surface. Synthetic gathers were then created at each of the 50 trace locations, again using the complete Zoeppritz equations and a 30 Hz Ricker wavelet, but this time with a far angle of 30°. The resulting synthetic gathers are shown in Figure 7. Note the amplitude “tuning” events at the edges of the anticlinal structure. Figure 8 shows the CDP stack of the synthetic gathers shown in Figure 7, with the P-wave sonic log inserted at CDP 25. From the stack, it is impossible to infer whether the bright amplitudes are created from hydrocarbon-bearing or nonhydrocarbon-bearing sands.

Figure 9 shows the AVO intercept (A) and AVO gradient (B) of the synthetic gathers shown in Figure 7. The top of the gas sand shows negative amplitudes on the intercept and gradient, and the base of the gas sand shows positive amplitudes on the intercept and gradient, which is indicative of a class 3 AVO sand. However, the top of the wet sand shows negative amplitudes on the intercept and almost no amplitude effect on the gradient. This is to be expected from the AVO curves shown in Figure 5.

Finally, Figure 10 shows the EEI/EPI results at 22° and 51° using the intercept and gradient results shown in Figure 7. These results were computed using equation 6. The 22° result shows a clear class 3 AVO anomaly from the top and base of the gas sand, and this result is reversed on the 51° result, indicating that the matrix of the sand is stiffer than the matrix of the shale when the fluid effects are removed. That is, the 22° result is responding to the fluid and the 51° result is responding to the matrix. At the top of the wet sand, we see the same negative response on the 22° and 51° results, indicating that there is no hydrocarbon anomaly present.

CASE STUDY

In our case study, the EPI method is applied to a shallow, highly porous gas sand play from southern Alberta. Figure 11 shows a display of a set of angle gathers over this gas sand reservoir. As indicated on the display, the angle range is from 0° to 30°. An integrated P-wave sonic log is inserted at common midpoint (CMP) 330, the gas well location. The gas sand is centered at 630 ms and extends from roughly CMP 322 to 340. This feature is called a class 3 AVO anomaly (Rutherford and Williams, 1989), and it manifests as a strong negative amplitude event overlying a strong positive amplitude, both with increasing absolute amplitudes from near to far angles. In Figure 11 and all subsequent seismic displays, the amplitude color bar has been normalized to standard deviations away from the mean, with red for positive amplitude and blue for negative amplitude, where the brightest colors show two standard deviations away from the mean. Although an angle range of 0°–30° is not large enough to adequately capture the AVO response in some areas of the world, it has done a good job for this shallow, highly porous gas sand.

Figure 12 then shows the P-wave velocity, density, VP/VS ratio, resistivity, and gamma ray logs for the gas well, along with the correlated synthetic trace (in red) and the extracted CMP gather from trace 330. All of the log curves were measured in situ, except for the VS log that was used to create the VP/VS ratio log. The VS log was created using the Arco “mudrock” line (Castagna et al., 1985) outside the gas zone and Biot-Gassmann modeling inside the gas zone. This modeling was done by assuming a 50% gas saturation over the zone between the markers’ top gas and base gas shown in Figure 12 and by computing the porosity using the density log. As expected, the VP/VS ratio shows a large drop over the gas sand zone, which has created the class 3 AVO anomaly seen on the seismic data.

However, there is another set of geologic events seen on the well-log suite, which contributed to the drilling of “false” anomalies throughout this area before the use of the AVO method. This consists of the hard streaks caused by carbonate stringers, which are seen below the gas sand on the P-wave and density logs. These hard streaks have been indicated in Figure 12. The hard streaks can create the appearance of large amplitude “bright spots” on the seismic stack, which was the main technique used to identify gas sands before the use of AVO and is sometimes still used. One of the purposes of the use of EEI/EPI in this area is to differentiate the hard streaks from the gas sand.

Figure 13 shows the CMP stack of the gathers shown in Figure 3, again with the integrated P-wave sonic log overlain at CMP 330. For comparison purposes, the color bar is again normalized to standard deviations away from the mean. The gas sand on the sonic log is centered at a time of 630 ms. Notice on the CMP stack in Figure 13 that the base of the gas sand (the large positive amplitudes in red at just greater than 630 ms) is quite visible but that the top of the gas (the large negative amplitudes in blue at just less than 630 ms) does not stand out from the continuous reflection event. In addition, there is a very strong negative event below the base of the gas sand, which is due to the hard streaks shown on the well log in Figure 12. As pointed out earlier, this event does not contain hydrocarbons, and it is often misinterpreted if only the stack is being used as an exploration tool.

Figure 14 shows the intercept and gradient attributes (A and B from equation 2) computed from the gathers in Figure 11. When comparing these attributes with the CMP stack shown in Figure 13, notice that the two strongest events on the intercept are the base of gas sand and the base of the hard streak below the gas sand, whereas the top of the gas sand is slightly weaker than on the stack. On the other hand, the strongest events on the gradient are the top of the gas sand and the base of the hard streak. Using the EEI/EPI technique, the optimum combination of intercept and gradient attributes shown in Figure 14 that will best separate the lithologic features, e.g., carbonate stringers, from the gas sand, can be found.

Figure 15 shows the application of the EPI/EEI technique to our gas sand example from Alberta. In Figure 15a, the EEI stack on the right was created by applying the EEI equation (equation 6) to the intercept and gradient attributes (A and B from equation 1) using a rotation angle of χ=22°. This is called the fluid stack. The EPI interpretation of this result is that kdry=0.45, which, as noted in Table 1, is a measured value for clean, high-porosity sands.

In Figure 15b, the EEI stack was created by applying an angle of χ=51° and is called the lithologic stack. The lithologic stack assumes that kdry=0, which agrees with EPI interpretation of pure shear modulus reflectivity. The fluid stack in Figure 15a clearly shows a gas sand centered at the gas well and a time of 630 ms, whereas the lithologic stack in Figure 15b largely ignores the gas sand and shows the set of carbonate stringers centered on 650 ms to the right of the gas well.

In summary, the application of the EEI/EPI technique in this area has allowed us to clearly differentiate the gas sand from the nonhydrocarbon-bearing hard streaks.

CONCLUSION

A new interpretation of the EEI method has been presented, which incorporates Biot-Gassmann poroelastic theory. This method is called EPI. The main advantage of the EPI method is that the χ term is recast as a continuous parameter that is dependent on the dry-rock properties of the reservoir, rather than a parameter whose value is estimated by correlation with other elastic parameters.

A discussion of the theory of the EPI method was first given, and then a table of representative dry-rock velocity ratios was used to give a physical basis for the method. These velocity ratios were then related to the χ angle, and it was found that a χ angle of 22.4°, corresponding to a dry rock S- to P-wave velocity ratio squared of 0.448, was optimum for the identification of a hydrocarbon anomaly. This observation tied together work that had been done in a rock-physics laboratory to the theory of EEI and to the Biot-Gassmann theory.

Next, a synthetic seismic example was used to illustrate the method. In this example, we created a four-layer earth consisting of overlying layers of shale and sandstone in which realistic petroelastic models were built for both lithologies. The two shales had identical parameters, but the top sandstone was 100% gas filled and had no shale content, whereas the basal sandstone was 100% brine filled and had 20% shale content. An anticlinal model was then created using the synthetic well log as the center of the anticline. From this model, synthetic seismic gathers were created using the Zoeppritz equations and analyzed using a simple CMP stack and the intercept/gradient approach followed by the EEI/EPI method. Because of the shale content of the basal sand, its P-wave velocity was actually lower than that of the gas-filled sand; thus, these reflection events looked equivalent on the stacked synthetic seismic section. But by using χ angles of 22 and −51, we were able to differentiate the wet sand from the gas sand.

Finally, the EPI method was applied to a gas sand field example from Alberta. The use of this data set indicated that, as in our numerical and synthetic studies, a dry rock S- to P-wave velocity ratio squared of 0.45, corresponding to χ=22°, did the best job of revealing the gas sand. Conversely, a dry rock S- to P-wave velocity ratio of zero, corresponding to χ=51° did the best job of revealing the lithologic anomalies, which, in this case, were thin carbonate stringers.

Although the term impedance has been used in extended elastic and poroelastic impedance, a more exact term would have been band-limited extended elastic reflectivity and poroelastic reflectivity. The reason that the term impedance is used is that the process is based on an underlying impedance concept, as discussed in the “EPI theory” section. Also, the band-limited reflectivity can be inverted to impedance in a straightforward way by the application of a modified poststack inversion scheme.

One caution in using this method is that it is not always clear what is meant by the term “dry-rock” when applied to in situ reservoir rocks. For clean sandstones and vugular carbonates, the meaning is unambiguous, but for shales, shaley sands, and fractured carbonates, there is no clear definition. In reservoirs with these types of rocks, the term “dry-rock velocity ratio” may be more of a tuning parameter for the process, rather than an actual measurable quantity.

ACKNOWLEDGMENTS

We would like to thank the anonymous reviewer who suggested that we needed a synthetic seismic example and Grace Ng at CGG, whose programming skills made the work of creating this example much easier than anticipated.

DATA AND MATERIALS AVAILABILITY

Data associated with this research are available and can be obtained by contacting the corresponding author.

DERIVATION OF EPI THEORY

Dong (1996) differentiates the P-wave velocity VP (as given in the first form on the right side of equation 15) with respect to bulk modulus K to obtain  
dVPdK=12ρVP.
(A-1)
Equating delta terms to derivatives gives  
ΔVPVP=ΔK2ρVP2.
(A-2)
Using the same approach, Dong (1996) also obtains the relationships among VP, μ, and ρ, as  
ΔVPVP=2Δμ3ρVP2
(A-3)
and  
ΔVPVP=Δρ2ρ.
(A-4)
Similarly, Dong (1996) differentiates the S-wave velocity VS with respect to μ and ρ to obtain  
ΔVSVS=Δμ2ρVS2
(A-5)
and  
ΔVSVS=Δρ2ρ.
(A-6)
Finally, we have the trivial result:  
Δρρ=1.
(A-7)
Putting this in a matrix form, Dong (1996) finds that  
[ΔVPVPΔVSVSΔρρ]=1ρ[12VP223VP212012VS212001][ΔKΔμΔρ].
(A-8)
If we recall the definitions of A, B, and C given in equation 2, this allowed Dong (1996) to rewrite equation A-8 as  
(ABC)=1ρ(14VP223VP21414VP253VP21414VP213VP214)(ΔKΔμΔρ).
(A-9)
Dong (1996) then rearranges equation A-9 to derive an optimum combination of intercept and gradient as follows:  
(5A+BCBAC)=1ρ(32VP20102VP200012)(ΔKΔμΔρ).
(A-10)
Whitcombe et al. (2002) use the inverse of equation A-10 to compute the elastic coefficients as a combination of A, B, and C. Without performing the full inversion, let us specifically look at the bulk modulus and note that the first and third equations can be written as  
5A+B=3ΔK2ρVP2+Δρρ
(A-11)
and  
2A2C=Δρρ,
(A-12)
respectively. Subtracting these two equations and rearranging gives  
ΔK=(3A+B+2C)ρVP21.5.
(A-13)
Noting that K=ρVP2(4/3)ρVS2, we can compute the bulk modulus reflectivity by dividing equation A-13 by this value, to get  
ΔK2K=3A+B+2C34k,
(A-14)
where k=(VS/VP)2. Whitcombe et al. (2002) then algebraically rearrange equation A-14 using the term f=C/A, giving  
RK=ΔK2K=(A+B3+2f)(3+2f34k).
(A-15)
By a similar analysis, Whitcombe et al. (2002) find that  
Rλ=Δλ2λ=(A+B2+f)(2+f24ksat)
(A-16)
and  
Rμ=Δμ2μ=(ABf)(f4ksat).
(A-17)
To derive poroelastic equations 19 and 20 given earlier, we start by noting that the P-wave velocity in equation 15 can be written as follows:  
VP=μd+pρ,
(A-18)
where d=(Kdry/μ)+4/3 and p=α2M. This allows us to rewrite equation A-8 as  
(ΔVPVPΔVSVSΔρρ)=1ρ(12VP2d3VP212012VS212001)(ΔpΔμΔρ).
(A-19)
Notice in equation A-19 that if we let d=4/3, the second term in the first row of the matrix on the right side becomes 2/3VP2, and we get exactly the same formulation as in equation A-8. That is, Δp becomes ΔK. Also, if we let d=2, the second term in the first row of the matrix on the right side becomes 1/VP2, and we get the λμρ formulation and Δp becomes Δλ. Then, we rederive the second and third steps of Dong (1996), as shown in equations A-9 and A-10. That is, equation A-9 becomes  
(ABC)=1ρ(14VP2d4VP21414VP28+d4VP21414VP2d4VP214)(ΔpΔμΔρ),
(A-20)
and equation A-10 becomes  
((8dd)A+BCBAC)=1ρ(2dVP204d2d02VP200012)(ΔpΔμΔρ).
(A-21)
Using the same approach as before, the porosity/fluid reflectivity term is computed, starting with the first and third equations in matrix equation A-21, or  
(8dd)A+B=2ΔpρVP2+(4d2d)Δρρ
(A-22)
and  
2A2C=Δρρ.
(A-23)
Eliminating the density term between equations A-22 and A-23 and rearranging gives  
Δp=(4A+dB+(4d)C)ρVP22.
(A-24)
Noting that we can rewrite p as  
p=ρVP2dρVS2.
(A-25)
Allows us to divide Δp by 2p to find the p reflectivity term given by  
Δp2p=(4A+dB+(4d)C)ρVP24(ρVP2dρVS2)=4A+dB+(4d)C44dksat,
(A-26)
where ksat=(VS/VP)sat2. Again, eliminating C using f=C/A, we then get  
RP=Δp2p=(A+d4+(4d)fB)(4+(4d)f44dksat).
(A-27)
Note that this is equation 18 from the main text. Using the fact that d=1/kdry allows us to derive the final form of equation (equation 19 in the main text):  
Rp=(A+B4kdry+(4kdry1)f)(4kdry+(4kdry1)f4kdry4ksat).
(A-28)
Equation A-28 is valid for the full range of kdry values, but note that if kdry=3/4, then equation A-28 reduces to Whitcombe et al.’s (2002),RK, if kdry=1/2 equation A-28 reduces to Whitcombe et al.’s (2002),Rλ, and if kdry=0 equation A-28 reduces to Whitcombe et al.’s (2002)Rμ.
Freely available online through the SEG open-access option.