Extended poroelastic impedance

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.

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)sat≈0.5$⁠, or $ksat=0.25$⁠, we find that $B≈RAI−2RSI=(ΔAI/2AI)−2(ΔSI/2SI)$⁠, where $RSI=(ΔSI/2SI)=(ΔVS/2VS)+(Δρ/2ρ)$ is the shear impedance reflectivity.

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 23–28.

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 ¹ 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 ¹, 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 ¹. 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 ¹ 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 ¹ is not a physical phenomenon but rather a mathematical artifact. The plot in Figure ¹ 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 ¹ was 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 $itan2(x,y)$ function, which uses all four quadrants of the $x$ and $y$ plane in the calculation. Figure ² 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.

The resulting synthetic well-log curves are shown in Figure ³. In both shales, the P-wave velocity is $3335m/s$⁠, the S-wave velocity is $1460m/s$ (giving a $VP/VS$ ratio of 2.28 and a Poisson’s ratio of 0.38), and the density is $2.28g/cm3$⁠. In the gas sand, the P-wave velocity is $2982m/s$⁠, the S-wave velocity is $2022m/s$ (giving a $VP/VS$ ratio of 1.48 and a Poisson’s ratio of 0.08), and the density is $2.03g/cm3$⁠. In the basal wet sand, the P-wave velocity is $2679m/s$⁠, the S-wave velocity is $1431m/s$ (giving a $VP/VS$ ratio of 1.87 and a Poisson’s ratio of 0.3), and the density is $2.15g/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 ⁴ shows the synthetic angle gather computed from the elastic logs shown in Figure ³, computed using a far angle of 45° (in Figure ⁴ 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 ⁵, which shows the picked amplitudes from the angle gather shown in Figure ⁴ for the three interfaces.

Figure ⁶ shows a 50-trace synthetic anticlinal model that was created by inserting the well shown in Figure ³ 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 ⁷. Note the amplitude “tuning” events at the edges of the anticlinal structure. Figure ⁸ shows the CDP stack of the synthetic gathers shown in Figure ⁷, 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 ⁹ shows the AVO intercept (⁠$A$⁠) and AVO gradient (⁠$B$⁠) of the synthetic gathers shown in Figure ⁷. 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 ⁵.

Finally, Figure ¹⁰ shows the EEI/EPI results at 22° and $−51°$ using the intercept and gradient results shown in Figure ⁷. 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.

In our case study, the EPI method is applied to a shallow, highly porous gas sand play from southern Alberta. Figure ¹¹ 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 ¹¹ 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 ¹² 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 ¹² 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 ¹². 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 ¹³ shows the CMP stack of the gathers shown in Figure ³, 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 ¹³ 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 ¹². 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 ¹⁴ shows the intercept and gradient attributes (⁠$A$ and $B$ from equation 2) computed from the gathers in Figure ¹¹. When comparing these attributes with the CMP stack shown in Figure ¹³, 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 ¹⁴ that will best separate the lithologic features, e.g., carbonate stringers, from the gas sand, can be found.

Figure ¹⁵ shows the application of the EPI/EEI technique to our gas sand example from Alberta. In Figure ¹⁵, 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 ¹⁵, 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 ¹⁵ clearly shows a gas sand centered at the gas well and a time of 630 ms, whereas the lithologic stack in Figure ¹⁵ 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.

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.

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 associated with this research are available and can be obtained by contacting the corresponding author.