Extended poroelastic impedance

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 in-version 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. (2012Zong et al. ( , 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 V Pi and V Piþ1 , S-wave velocities V Si and V Siþ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 R PP ðθÞ ¼ a (1) , where θ iþ1 ¼ sin −1 ðsin θ i ðV Piþ1 ∕V Pi ÞÞ.In the term k sat , the subscript sat stands for the saturated, or in situ, velocities.In some formulations, the factor two in the denominator of the three terms ΔV P ∕2V P , ΔV S ∕2V S , 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 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 (2) , and C ¼ ðΔV P ∕2V P Þ.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 ΔV P ∕2V P 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 R VP .However, by using the logarithmic differential approximation from calculus, we find that ðΔP∕2PÞ ≈ ðΔ ln P∕2Þ for any parameter P. Using this formulation, A can be thought of as the acoustic impedance reflectivity, or A ¼ R AI ¼ ðΔAI∕2AIÞ ¼ ðΔV P ∕2V P Þ þ ðΔρ∕2ρÞ, where AI ¼ ρV P is the acoustic impedance.In addition, if we make the assumption that ðV S ∕V P Þ sat ≈ 0.5, or k sat ¼ 0.25, we find that B≈R AI −2R SI ¼ ðΔAI∕2AIÞ−2ðΔSI∕2SIÞ, where R SI ¼ ðΔSI∕2SIÞ ¼ ðΔV S ∕2V S ÞþðΔρ∕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: R EI ðθÞ ¼ ΔEIðθÞ 2EIðθÞ ≈ Δ ln EIðθÞ 2 ; (3) where 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 sin 2 θ to tan χ to get (5) Finally, they divide by cos χ to get the final form of the EEI reflectivity, which we will call the scaled EEI reflectivity, or 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 where p ¼ cos χ þ sin χ, q ¼ −8k sat sin χ, and r ¼ cos χ− 4k sat sin χ.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 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 ρ ¼ aV 0.25 P : 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 ðΔρ∕ΔV P Þ ¼ ða∕4ÞV −0.75 P .This then gives us the relationship ðΔρ∕ρÞ ¼ 1∕4ðΔV P ∕V P Þ, which upon substitution into the intercept term A from equation 2 gives A ¼ ð5ΔV P ∕8V P Þ.Using the expression for C in equation 2, this finally gives a value of f ¼ 0.8 and leads to χ angles of for the bulk modulus reflectivity, for the first Lamé constant reflectivity, and 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 singlevalued 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 and where K sat 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 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 and where K dry is the dry bulk modulus, λ dry is the dry second Lamé parameter, and p ¼ α 2 M 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 where K m is the mineral bulk modulus, K f 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): Using equations 16 and 22, Russell et al. (2003) derive a generalized form of the λρ equation, written as where I Psat = saturated P-wave impedance, or ρ sat V Psat , I Ssat = saturated S-wave impedance, or ρ sat V Ssat , and 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 To be more consistent with EEI, equation 25 can be reformulated using k dry , the inverse of d, giving where k dry ¼ ð1∕dÞ ¼ ðV S ∕V P Þ 2 dry .We will now show that equations 25 and 26 are a generalization of equations 8-10.
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 V P -to-V S ratio in equation 25 with the equation or, using the inverse of d in equation 26 as 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 V P -to-V S ratio squared d, the inverse of d, or the V S -to-V P ratio squared k dry , the dry rock Poisson's ratio ν dry , the dry bulk modulus over shear modulus ratio K dry ∕μ, and the dry first Lamé parameter over shear modulus ratio λ∕μ.Using this table gives us a way to physically interpret equations 23-28.
In Table 1, two values are of particular importance: d ¼ 4∕3 and d ¼ 2. When d is equal to 4/3 (or k dry ¼ 0.75), this implies that the K dry ∕μ ratio equals zero and the term p is equivalent to the bulk modulus K sat .Substitution into equation 25 or 26 gives which is identical to equation 8 or R K .When d is equal to 2 (or k dry ¼ 0.5), this implies that the λ dry ∕μ ratio equals zero and p is equivalent to λ sat : Substitution into equation 25 or 26 gives which is identical to equation 9 or R λ .Finally, consider the limiting value as d approaches ∞ (or, equivalently, as k dry approaches zero).Substitution of k dry ¼ 0 into equation 26 gives 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 k dry ¼ 0.429), gives a K dry ∕μ ratio of 1.0 and is often used to model clean sands, and the value of d ¼ 3 (or k dry ¼ 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 k dry ¼ 0.448) corresponds to a K dry ∕μ 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 k dry 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 k dry 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, k dry values of 0.5 and 0.75 correspond to θ values of 28°and 37°, respectively, whereas a k dry 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 k dry ¼ 1∕9.Also of note in Figure 1 and equation 26 is that the χ rotation angle decays asymptotically to zero as the value of k dry increases and that the curve is independent of f when k dry ¼ 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 Figure 1.A plot of dry rock ðV S ∕V P Þ 2 versus χ for three different values of the C∕A ratio f, made using the atanðzÞ function.Note the discontinuity where the denominator of equation 20 equals zero.created using the tan −1 ð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 i tan 2ð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 k dry values of zero and one.
As a summary of these results, Table 2 shows a subset of the dry rock d ¼ ðV P ∕V S Þ 2 and k dry ¼ ðV S ∕V P Þ 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 k dry value of 0.448) ties together work done in a rockphysics 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 (V sh ) 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 and where the mineral moduli were computed from a Reuss average of the shale and quartz constituents given by with K qz ¼ 37 GPa, μ qz ¼ 45, K sh ¼ 21 GPa, and μ sh ¼ 5.3.The bulk density was computed from a Voigt average of the brine and mineral components given by where with ρ qz ¼ 2.65 g∕cm 3 and ρ sh ¼ 2.650 g∕cm 3 .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 V P and V S of the shale were then computed from the effective bulk and shear moduli and densities by 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 (P eff ), coordination number (C), and critical porosity (ϕ 0 ).For this model, we used P eff ¼ 25 MPa, C ¼ 9, and ϕ 0 ¼ 0.4.The mineral modulii and bulk-density values were computed as in equations 34-37, 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 V P and V S 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 V P ∕V S ratio of 2.28 and a Poisson's ratio of 0.38), and the density is 2.28 g∕cm 3 .In the gas sand, the P-wave velocity is 2982 m∕s, the S-wave velocity is 2022 m∕s (giving a V P ∕V S ratio of 1.48 and a Poisson's ratio of 0.08), and the density is 2.03 g∕cm 3 .In the basal wet sand, the P-wave velocity is 2679 m∕s, the S-wave velocity is 1431 m∕s (giving a V P ∕V S ratio of 1.87 and a Poisson's ratio of 0.3), and the density is 2.15 g∕cm 3 .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 V P ∕V S 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 Figure 3.The synthetic sonic logs generated from the petroelastic models shown in Table 3.  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 nonhydrocarbonbearing 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°u sing 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 dis-   play 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 ampli-   tude, 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, V P ∕V S 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 V S log that was used to create the V P ∕V S ratio log.The V S 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 V P ∕V S 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 welllog 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 stan-Figure 12.The log suite for the gas well, along with the correlated synthetic and extracted gather at CMP 330 from the seismic gathers shown in Figure 11.dard 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 k dry ¼ 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 k dry ¼ 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.

andFigure 2
Figure 2. A plot of dry rock ðV S ∕V P Þ 2 versus χ for three different values of the C∕A ratio f, made using the atan2ðy; xÞ function, in which the plot monotonically decreases between k dry values of zero and one.

Figure 4 .
Figure 4.The synthetic angle gather computed from the elastic logs shown in Figure 3.

Figure 5 .
Figure 5.The picked amplitudes from the angle gather shown in Figure 4 for the three interfaces.

Figure 6 .
Figure 6.A 50-trace synthetic anticlinal model created by inserting the well shown in Figure 3 at CDP 25.

Figure 7 .
Figure 7.The 50 CMP prestack synthetic seismic anticlinal model created from the model shown in Figure 6.

Figure 8 .
Figure 8.The CMP stack of the synthetic gathers shown in Figure 7, with the P-wave sonic log inserted at CMP 25.

Figure 9 .
Figure 9.The (a) AVO intercept (A) and (b) AVO gradient (B) of the synthetic gathers shown in Figure 7.

Figure 10 .
Figure 10.The (a) EEI/EPI result at 22°and (b) EEI/EPI result at −51°found by applying equation 6 in the text to the intercept and gradient results shown in Figure 9.

Figure 11
Figure 11.A set of angle gathers over a shallow Alberta gas sand.The color bar is normalized to standard deviations away from a mean of zero.

Figure 13 .
Figure 13.The CMP stack of the angle gathers from Figure 11.For comparison purposes, the color bar is again normalized to standard deviations away from the mean.

Figure 14 .
Figure 14.(a) The intercept attribute computed from the gathers in Figure 11 and (b) the gradient attribute computed from the gathers in Figure 11.

Table 1 .
Ranges of values of d and k dry and their relationship with other elastic constant ratios.

Table 2 .
A subset of the dry rock d V P ∕V S 2 and k dry V S ∕V P 2 values shown in Table1, with the χ angle included in the right column.The input parameters for the petroelastic models used in the synthetic study discussed in the text.