## ABSTRACT

Accurate inversion of time-lapse amplitude variation with offset (AVO) to reservoir pressure change requires knowledge of strain/stress changes in the reservoir and overburden. These geomechanical effects have not been taken into account in previous studies. We have determined that there are three contributions of importance that must be included in the AVO equations: a density term in the reservoir, the seismic anisotropy generated by the reservoir deformations, and the impact of the changing overburden properties. To calculate the error in neglecting these contributions, we derived analytic expressions for the 4D intercept and gradient as functions of the strains developed in the reservoir and overburden during production. The 4D AVO behavior is found to be controlled by three $R$-factors that relate to time-lapse variations in P- and S-wave velocities. These factors can be estimated from rock-physics measurements in the laboratory or calibrated to field data measurements of time shifts via the well-known Hatchell-Bourne-Røste $R$-factor. Numerical computation indicates that the error in not accounting for the anisotropic contribution to the P-P reflectivity is 50%–300% for the 4D zero-offset intercept and 62%–155% for the 4D gradient. Thus, in most applications, the time-lapse anisotropy in the reservoir and the overburden cannot be neglected when estimating pressure changes using 4D AVO.

## INTRODUCTION

The use of time-lapse amplitude variation with offset (AVO) has gained popularity as a possible approach for obtaining quantitative estimates of reservoir pressure and saturation changes from 4D seismic data. A theoretical framework for this method is published by Landrø (2001), who develops expressions for the interpretation of intercept and gradient at a top reservoir interface. Like all AVO approaches, the method relies upon a calibration of the governing rock-physics coefficients for the reservoir followed by a nonlinear or linear inversion step. Over the years, there have been several incremental improvements to this 4D AVO theory and its implementation. Examples are the use of a quadratic pressure term (Meadows et al., 2001), multicomponent data (Stovas and Landrø, 2004), and an engineering-consistent multiattribute approach (MacBeth et al., 2006). Despite the theoretical expectations, there are relatively few practical case studies that convincingly validate the theory per se, although anecdotally it is understood that qualitative AVO predictions roughly hold (see, e.g., Corte et al. [2019] who, as with several others, observe that the saturation changes are more pronounced on the far-offset stacked data in accordance with the theoretical gradients). Several reasons may be envisaged for the lack of good case studies; however, we consider an overriding factor to be a difficulty in obtaining clean and trustworthy 3D AVO signatures. This current paper contributes to this discussion by demonstrating that geomechanical effects in the reservoir and the overburden must also play a role in the quantitative interpretation. Because the reservoir undergoes triaxial deformation in general (Fjaer et al., 2008), we must necessarily introduce production-induced elastic property changes that are anisotropic (e.g., Nur and Simmons, 1969). The overburden properties are also dependent on the stress and strain changes transferred to these nonreservoir rocks. These overburden changes due to production-induced reservoir deformations are therefore also anisotropic. The reservoir and overburden anisotropy must be taken into account to ensure the physical correctness of an AVO inversion (Herwanger and Horne, 2005). To obtain an expression for the error in not taking these effects into account, we redevelop the time-lapse AVO equations in a framework better suited for direct geomechanical interpretation using 4D seismic data. Following developments in time-lapse time shifts, we choose to relate the wave velocity and density changes to anisotropic strain changes in the reservoir and overburden. Therefore, we follow the paradigm of Hatchell and Bourne (2005a) and Røste et al. (2005) and the subsequent development of Rodriguez-Herrera et al. (2015). The resultant strain-dependent equations are found to be simple in form, easy to interpret, and also compatible with the rock-physics driven approach by Landrø (2001). The equations provide a means of reexpressing the 4D AVO in a geomechanics friendly form that captures the strain-induced anisotropy.

## THEORY

*’*s before and after production and then we subtract to obtain the time-lapse change. This calculation is performed for the sand using equation 4 and the shale using equation 5, with the relevant initial moduli and third-order coefficients inserted. To further simplify the calculation, we introduce “$R$” factors following the original proposition by Hatchell and Bourne (2005a) and Røste et al. (2005) and the development of Herwanger (2008) and Rodriguez-Herrera et al. (2015). Specifically, the latter two authors define three $R$ factors $R1$, $R2$, and $R3$ as

### Geomechanical changes only in reservoir rocks

The above description is satisfactory if no corrections are necessary for the geomechanics, but we have shown that this is required. Therefore, for accurate application of the 4D AVO, the two error terms in equation 17 must be applied: one for density and one for the anisotropy. To determine how big a problem this might be, we consider that $R1res$ is typically approximately 2 (MacBeth et al., 2018) — in which case, the density term can represent a 50% error. Assuming $R1>R2$ for the anisotropic error term $(R1res\u2212R2res)/{\u2212R1res+2(R1res\u2212R2res)}$ we arrive at an estimate of the error of 100%. Both errors are therefore not trivial: Adjustments must be made, and equations 12 and 14 must be used in the analysis. In practical terms, the 4D AVO gradient should be larger than originally estimated by the isotropic equation and the contribution of pressure change to the far offsets will diminish faster with offset than anticipated.

### Contribution from geomechanical changes in the shale

A rough estimate of $\epsilon 33ob/\epsilon 33res$ can be determined by considering the analytic solution of Geertsma (1973), which provides the simplest and most popular route to assess how the overburden deforms in the presence of a compacting reservoir. The solution assumes a nucleus of strain embedded in a uniform elastic half-space, a 2D flat geometry for the reservoir, and no mechanical contrast between the reservoir and the overburden (Fjaer et al., 2008). By summing this solution over the reservoir volume (Hodgson, 2009), it is possible to predict the overburden strain for a variety of reservoir geometries. Considering only the overburden strain along a vertical axis through the center of a disc-shaped reservoir, it is found that Geertsma predicts $\epsilon 33ob/\epsilon 33res=\u22120.5\xd7(h/D)$, where $h$ is the reservoir thickness and $D$ is the reservoir depth (see Appendix B and Hatchell and Bourne, 2005b). Taking, as an example, a reservoir of thickness 40 m lying at a depth of 2 km, this gives a strain ratio of 1/50 and a negligible error. For typical overburden and reservoir HBR $R$-factors of 5 and 2, respectively, this results in an intercept error of only 3.3%. The strain ratio does, however, vary with the reservoir depth and lateral dimension.

In reality, there are many situations for which the vertical strain ratio $\epsilon 33ob/\epsilon 33res$ can be significant:

- •
One situation may be if there is pressure communication between the reservoir and the immediate surroundings. In this case, the reservoir and the overburden suffer from depletion with similar strains; therefore, the ratio is large and positive.

- •
An important controlling factor is the mechanical contrast between the reservoir and the overburden layers. For example, if the reservoir is composed of soft sands, then stress is concentrated in the stiffer surrounding layers and helps shield the depleting zone from the weight of the overburden. Large strains are developed in the reservoir compared to the overburden; thus, we might expect the strain ratio to be very small.

- •
Conversely, if the reservoir is composed of stiff cemented sandstones, then the strain development in the overburden may be enhanced. A stiff layer lying above the reservoir shields the displacement from propagating upward to the free surface. In this case, there would be a higher than normal strain development between the stiff layer and the reservoir; thus, the vertical strain in the overburden would contribute significantly to the intercept.

- •
For a compact region of reservoir depletion, stress arching occurs, which shields the reservoir from the overburden weight and limits the strain development.

- •
A thin soft layer, such as a shale in a stacked reservoir interval, absorbs the strain and prevents propagation upward into the overburden (Rangel et al., 2016). This is known to occur in the shales overlying a stiff chalk (e.g., De Gennaro et al., 2008).

To capture the full diversity of the above effects, the vertical strain ratio may be considered to vary between $+1$ and $\u22121$. For the extremes of this phenomenon, the error in the intercept is now observed to be sizable and up to 300%.

### Discussion of the approximations

^{o}, 45

^{o}, and 90

^{o}to the symmetry axis of a sample. The inverted TOE coefficients can be converted to the $R$-factors using equation 6. As an example, Table 1 provides a compilation of such values from the work of Winkler and Liu (1996), Sarkar et al. (2003), Prioul and Lebrat (2004), and Prioul et al. (2004). It is immediately noticeable that these laboratory values are large in magnitude and somewhat scattered with no distinct trend. However, we recognize that $|R1|\u226b|R2|>|R3|$, and this may suggest that further approximations/assumptions can be made. Because our study is focused on solutions relating to seismic data acquired in the field, we also consider the field observations. In particular, we consider the understanding based on Røste et al. (2005) and Hatchell and Bourne (2005a), who originally define an $R$-factor (which we term as the original $R$-factor, $R$) to determine the link between the fractional vertical velocity change and the vertical strain

^{1}compares the AVO behavior for the four levels of approximation as defined by the equations above, and Figure

^{2}shows the errors as a function of the ratio of vertical strains for the overburden and reservoir. For a dynamically active overburden, we distinguish between overburden strains equally as big as those in the reservoir but with opposite polarities, i.e., $\epsilon 33ob=\xb1|\epsilon 33res|$, but we also calculate the intercept and gradient for more typical values of $10\u22124$. Although the trend of decreasing amplitude with offset is not changed between the different approximations, there is significant difference in the overall error made by the original equations. The original equation is in error relative to the fully corrected equation by at least 50% for the intercept, and this can reach 300% for some values of the strain ratio. The gradient errors are at least 62.5%, rising to a possible 155%. These are large errors and therefore cannot be ignored in a quantitative analysis for pressure estimation using 4D seismic data.

## DISCUSSION

The AVO equations resulting from our development are based on several assumptions and approximations, which could limit the widespread application of this method. These fall into two groups.

### Model for strain dependence

The TOE model used in this work is that of Prioul et al. (2004), who assume an isotropic symmetry for the third-order coefficients. A more complete approach would be to consider an anisotropic symmetry and thus work with a greater number of coefficients (see Fuck et al., 2009). It is understood (Bakk et al., 2020) that, for shales, this more generalized form may be necessary. Our justification for choosing the smaller number of coefficients is that the laboratory data of Prioul et al. (2004) validate the isotropic assumption to within a few percent using a combination of hydrostatic and biaxial data. Therefore, we are currently satisfied that the reduction in the number of parameters from the isotropic assumption outweighs the need for completeness, even for the anisotropic overburden. There is also some opinion on the strain-polarity asymmetry of the $R$-factor model (Hatchell and Bourne, 2005a), which the TOE model does not address because it is linear in strain. However, the arguments for this asymmetry are difficult to separate from the lithology dependence (Røste et al., 2015; MacBeth et al., 2018). The laboratory measurements of Winkler and Liu (1996), Prioul et al. (2004), and Sarkar et al. (2003) do consider nonhydrostatic paths, but not explicitly unloading. Unloading studies of shales (Holt et al., 2018) and sandstones (Holt and Stenebråten, 2013) do indicate hysteresis effects and possible asymmetries in the loading/unloading response. These results suggest a modification for the TOE model may be required in practice. Finally, of concern also is that the TOE model does not contain an explicit dependence on prior stress/strain history, particularly for the shale overburden. There has been substantial discussion in the literature regarding stress state and stress path being important for the stress/strain dependence for shales (Holt et al., 2018) and to a lesser extent for sandstones (Holt and Stenebråten, 2013). The adaptation of the TOE model to such situations is still under investigation (Bakk et al., 2020). We propose that the general conclusions of our study are not dependent on the strain model that we choose, and will hold in principle for most models of production-induced anisotropy. However, specific details and parameterizations may vary if, for example, a crack model is used.

### Geomechanical assumptions

The understanding developed in the current work is based on two commonly held assumptions: zero volumetric strain in the overburden and uniaxial strain in the reservoir. These lead to simplified orthorhombic overburden changes and VTI reservoir changes due to pore pressure decrease or increase when third-order elasticity theory is applied. Based on this theory, we developed an analytic expression for 4D AVO to gain insight into the seismic response due to geomechanical deformations. Reality is, as expected, more complicated than the assumptions that we have made as can be attested by running a geomechanical simulator for any field-scale model. Although we find the uniaxial strain condition holds well in our experience, the assumption of zero volumetric strain (or mean stress) in the overburden is much more easily violated. Recall that the assumption is based on Geertsma (1973), who assumes no contrast between the reservoir and the surroundings. We know that in practice that stress/strain paths in the overburden may significantly deviate from this ideal depending on the property contrast between the reservoir and overburden, geometry, and lateral variations. Therefore, our theory is strictly limited to reservoirs that lie in uniform, flat-lying geology for which the geomechanical property contrast between the reservoir and the overburden is small. A more complete generalization of our result, which captures the field architecture and mechanical property variations, requires the use of a full numerical simulator (see, e.g., the work of Jaramillo et al., 2019). Geomechanical simulation would provide access to a more comprehensive set of anisotropic strain components, including the distribution and orientation of the principal strains. The approximations used here to build our conceptual understanding would not then be applicable, but numerical calculation via the TOE model would allow the reflection coefficients to be calculated using full anisotropic reflectivity equations. Thus one might proceed as in Fuck et al. (2009), to calculate the reflectivity from an anisotropic-anisotropic interface, with anisotropy of arbitrary symmetry in the overburden. Although this remains a compelling step, we find that it is of considerable practical interest and insight to reduce the number of parameters describing the AVO behavior to a smaller number. Indeed, our formulation has shown that major effects can be approximated by the two parameters $R1$ and $R2$ linked to well-known empirical parameters from time-shift analysis. Thus, a future direction for this work is to validate this simpler approach on synthetic seismic in the presence of heterogeneity and real data. Another avenue to be assessed is the effect of the improved formulation on the joint inversion of pressure and saturation. Finally, it is recognized that horizontal strain is an important component to understand in practice. The relations developed here can also be used to invert for vertical strain and the ratio of horizontal strains in the overburden from 4D AVO.

## CONCLUSION

All hydrocarbon reservoirs undergoing production (and hence recovery) are geomechanically active to some degree or other because the evolution of reservoir and overburden strains with field production life is inevitable. Thus, the problems highlighted above are present for time-lapse AVO because production-induced anisotropy is always generated and cannot be dismissed. We show that for compacting (or expanding) reservoirs, this production-induced anisotropy arises mainly through changes in Thomsen’s $\delta $ influencing the AVO gradient (assuming here that two-term AVO is a valid approximation for our data), and it is sufficient to lead to a sizable error in pressure estimation. Additional factors not previously taken into account include a density effect in the reservoir that cannot be neglected, and also an overburden contribution due to deformation of the layers overlying the reservoir. Although this overburden strain is possibly a smaller factor than the other two sources of error, it may become important under specific field conditions. The geomechanical response of the reservoir and overburden is complex and is highly dependent on the reservoir thickness, the geometry, the mechanical contrast between the reservoir and overburden, and faulting, in addition to general heterogeneity, and it requires numerical computation and evaluation on a case-by-case basis. We have outlined extreme end-member cases to this behavior in our work.

## ACKNOWLEDGMENT

The authors thank the sponsors of the Edinburgh Time-Lapse Project, Phase VII: AkerBP, BP, CGG, Chevron/Ithaca Energy, CNOOC, Equinor, ConocoPhillips, ENI, Petrobras, Norsar, Woodside, Taqa, Halliburton, ExxonMobil, OMV and Shell for their financial support and technical discussion. CM also thanks Joerg Herwanger for initially posing this interesting geomechanical challenge, and for fruitful and interactive discussions during development.

## DATA AND MATERIALS AVAILABILITY

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

### UNIAXIAL COMPRESSIBILITY VERSUS HYDROSTATIC COMPRESSIBILITY

*ν*)/3(1 − ν). The compressibility $Cpp$ is known to be a function of the pore pressure, porosity, effective pressure, and lithology. It increases with the increasing pore pressure, and it decreases with the increasing effective confining pressure. Typical pore volume compressibilities at low effective pressures (less than 6.90 MPa or 1000 psi) for clastic reservoirs are in the range of 4.4 to 11.6 × 10

^{−4}MPa

^{−1}(3 to 8 × 10

^{−6}psi

^{−1}). However, for strongly compacting chalk reservoirs (such as Ekofisk and Valhall), compressibility values can reach 72.5 to 145 × 10

^{−4}MPa

^{−1}(50 to 100 × 10

^{−6}psi

^{−1}).

### CONVERTING HYDROSTATIC VELOCITY PERTURBATIONS TO UNIAXIAL VELOCITY PERTURBATIONS

### CURVATURE TERM IN THE 4D AVO

Biographies and photographs of the authors are not available.