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.

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.

Consider the case of a shale overburden with initial (preproduction) vertical transverse isotropy (VTI) wave properties overlying a VTI sandstone reservoir (anisotropic due to shale and clay deposits). In this problem, the elastic properties of the shale and sandstone are defined by the density ρ, the vertical P- and S-wave velocities α and β, and the Thomsen parameters ε and δ (because we are considering P-P waves only). Following Rüger (1997), the P-P reflection coefficient at this discrete sand-shale interface is
RpVTI(θ)=12ΔcZZ¯+12{Δcαα¯(2β¯α¯)2Δcμμ¯+Δcδ}sin2θ+12{Δcαα¯+Δcε}sin2θtan2θ,
(1a)
where μ is the shear modulus, Z is the P-wave impedance, α is the vertical P-wave velocity, β is the vertical S-wave velocity, θ is the approximate incidence angle, and we define the operator Δc to represent a small parameter contrast across the interface such that, for example, Δcα=α(sand)α(shale). We break with conventional notation here because we wish to retain the use of Δ for time-lapse changes in the development below. In addition, the overscore refers to an average of the properties either side of the interface. If the anisotropic parameters ε and δ are zero, it can be shown, by rearrangement to highlight the three perturbation terms Δα/α, Δβ/β, and Δρ/ρ, that equation 1a agrees with the isotropic expression of Aki and Richards (1980). Note that Rüger (1997) equation is organized as a normal incidence term, a “gradient” term of sin2θ, and a higher order term of sin2θtan2θ. This is different from Landrø (2001) who follows the first Smith and Gidlow (1987) rearrangement of Aki and Richards (1980) (specifically their equation 2) by formulating the AVO equation in terms of normal incidence, sin2θ and tan2θ. For comparison purposes with the Landrø publication, we also show the Smith and Gidlow (1987) version of Rüger (1997):
RpVTI(θ)=12ΔcZZ¯12{4(β¯α¯)2Δcμμ¯Δcδ+Δcε}sin2θ+12{Δcαα¯+Δcε}tan2θ.
(1b)
From these equations, we immediately observe the impact of the anisotropy in the reservoir and overburden as an additional parameter contrast in the AVO coefficients. These terms will be important when we consider time-lapse changes. To determine whether the additional factors from anisotropy are of sufficient magnitude to be important to the pressure inversion problem, we must first determine the time-lapse version of this AVO equation. For this development, we maintain the use of Rüger (1997) expression because it expresses the progression to higher order terms in the incidence angle. Our assumption is that the VTI equations remain valid for wave propagation before and after production (we will later examine and validate this assumption). To derive the time-lapse reflectivity ΔR, from the AVO equation in equation 1a, we assign the corresponding anisotropic elastic properties to the sand and shale before and after production and then we subtract to obtain the time-lapse change. Our convention is to difference the monitor (postproduction) reflectivity from the baseline (preproduction) reflectivity. Taking this time-lapse difference and noting that the operators Δ and Δc are linear and interchangeable allows us to readily transform equation 1a into a time-lapse equation, arriving at the general result (in Rüger [1997] format)
ΔR(θ)=ΔI+(ΔG+Δδ)sin2θ+(ΔH+Δε)sin2θtan2θ,
(2)
where ΔI is the time-lapse change of the AVO intercept, ΔG is the change of the nonanisotropic gradient term in equation 1a, and ΔH is the higher order curvature term. These three terms compare directly with Landrø (2001) (taking into account the rearrangement) if we replace time-lapse perturbations in the isotropic velocities by time-lapse perturbations in the anisotropic velocities. The anisotropic terms Δδ=(δresafterδresbefore)(δobafterδobbefore) and Δε=(εresafterεresbefore)(εobafterεobbefore) (where the subscript ob refers to the overburden and res refers to the reservoir) are additional contributions that relate directly to the anisotropy induced when the sands and shales change under stress (and strain) due to production and cannot be neglected. Even if the overburden changes could be ignored, the production-related change in Thomsen’s parameter δ for the reservoir must still be taken into account. To develop this observation further, we must establish a framework for evaluating the stress (or strain) induced anisotropy and the time-lapse perturbations: Δα/α, Δρ/ρ, and Δμ/μ.
Consider a time-lapse change in the reservoir sand induced by the extraction of hydrocarbon volumes. We neglect changes in fluid saturation to focus only on the geomechanical effects, on the understanding that the saturation effects are better determined and can be readily included later if desired. Note that in some fields such as those in the Southern Gas Basin or Northern North Sea, United Kingdom Continental Shelf, the assumption of geomechanics as the sole contributor is valid (Brain et al., 2017) because it may be a gas reservoir with no upward movement of the aquifer. As anticipated from geomechanical principles, the reservoir compacts in response to fluid production, whereas the overburden will sympathetically extend. If the reservoir is narrow in the horizontal dimension or the pressure depletion is locally confined, stress arching will also be very visible (see, e.g., Herwanger, 2008). Thus, unlike the case of Landrø (2001) for which pressure changes are considered to occur only in the reservoir itself and time-lapse changes are limited to sandstone property changes, we also consider the possibility of stress/strain changes in the surrounding shale. To calculate the desired time-lapse change in the AVO signature, we describe the shale and sand wave properties by a density ρ and stiffness tensor represented as a 6×6 matrix of elastic constants cij by using the Voigt notation (where i=1,6 and j=1,6). We seek to express changes in these elastic constants as a function of principal strains created by the applied stress field. We assume these to be oriented horizontally (ε11 and ε22) and vertically (ε33). For this purpose, we use the third-order elasticity model that relates postproduction moduli cij to preproduction moduli cij0 and the production-induced strains ε11,ε22, and ε33 via third-order coefficients cijk that characterize the nonlinear properties of the material (Prioul and Lebrat, 2004; Prioul et al., 2004). This model uses approximations to nonlinear elasticity theory for anisotropic media, and specifically for initially VTI media. It has been shown by these authors to be satisfactory for orthorhombic media. The model is simpler to understand and implement as a function of strain rather than stress. Importantly, in their experimental studies, Prioul et al. (2004) and Prioul and Lebrat (2004) validate the underlying assumption of an isotropic sixth-order tensor, and hence their desired approach. Explicitly, this model can be written as
c11=c110+c111ε11+c112(ε22+ε33),c22=c110+c111ε22+c112(ε11+ε33),c33=c330+c111ε33+c112(ε11+ε22),c12=c120+c123ε33+c112(ε11+ε22),c13=c130+c123ε22+c112(ε11+ε33),c23=c130+c123ε11+c112(ε22+ε33),c44=c440+c144ε11+c155(ε22+ε33),c55=c440+c144ε22+c155(ε11+ε33),c66=c660+c144ε33+c155(ε11+ε22).
(3)
In these equations, there are only three independent coefficients c111, c112, and c123 as c144=(c112c123)/2 and c155=(c111c112)/4. According to the relations in equation 3, the initial state defined by the preproduction moduli cij0 may be isotropic or anisotropic. It should be noted that because the original approximations of Prioul et al. (2004) are linear in strain, we choose to develop our subsequent equations in strain because conversion to stress via Hooke’s law would result in a more complex set of final equations (Prioul et al. [2004] show how this conversion may be accomplished). Another reason for our choice is to provide compatibility with 4D time-shift theory, which is also written in terms of strains. To simplify equation 3 for the purposes of deriving a time-lapse AVO expression, we make two assumptions. First, in the reservoir, we assume a stress state that creates uniaxial strain (i.e., with no lateral displacement, such that ε11=ε22=0,ε330). This condition is valid when the thickness of the producing unit is much smaller than its lateral extent. Because the reservoir is typically long, thin, and at depth, we consider the assumption of uniaxial deformation to be a good approximation. Under this assumption, displacements inside the reservoir are purely vertical in response to a drop, or rise, in fluid pore pressure such that there is no lateral deformation, no stress arching, and the full weight of the overburden is experienced by the reservoir. For these circumstances, equation 3 reduces to a more compact form
c11=c110+c112ε33,c22=c110+c112ε33,c33=c330+c111ε33,c12=c120+c123ε33,c13=c130+c112ε33,c23=c130+c112ε33,c44=c440+c155ε33,c55=c440+c155ε33,c66=c660+c144ε33.
(4)
From inspection of equation 4, it may be readily observed that an initially VTI sand remains VTI under uniaxial strain. The corresponding density change for uniaxial deformation is given by Δρ/ρ=ε33. In contrast, the overburden in the immediate vicinity of the top reservoir will experience a stress state in which the total volumetric strain is zero (i.e., ε11+ε22+ε33=0) (see, e.g., Fjaer et al. [2008] — this assumes a Geertsma [1973] limit with no elastic contrast between the reservoir and the surroundings). Clearly, there is no density change in this case because there is no change in volume. The elastic stiffnesses in equation 3 are modified specifically for the overburden as
c11=c110+(c111c112)ε11,c22=c110+(c111c112)ε22,c33=c330+(c111c112)ε33,c12=c120+(c123c112)ε33,c13=c130+(c123c112)ε22,c23=c130+(c123c112)ε11,c44=c440+(c144c155)ε11,c55=c440+(c144c155)ε22,c66=c660+(c144c155)ε33.
(5)
Here, it may be observed by inspection that the initially VTI shale in its preproduction state now has orthorhombic symmetry (nine independent elastic constants). This system possesses two vertical planes of symmetry aligned along the principal strain directions: the vertical [x1, x3] and [x2, x3] planes. In each of these planes, the wave behavior may be regarded as VTI (Tsvankin, 1997; Rüger, 1998). In our work, we follow Rüger (1998) in the construction of VTI AVO equations for the two vertical symmetry planes [x1, x3] and [x2, x3].
To derive the time-lapse reflectivity ΔR from the AVO equation in equation 1a, we calculate the corresponding ρ, α, β, ε, and δ from the cijs 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
R1=12c111C330,
(6a)
R2=12c112C330,
(6b)
and
R3=12c123C330.
(6c)
This, of course, implies that by definition the R-factors are a function of the elastic constants at the initial (reference baseline) state of the rock. Note that, for completeness, by application of symmetry rules
R1R24=12c155C330
(7a)
and
R2R32=12c144C330.
(7b)
The third-order coefficients can be readily calculated in the laboratory if time-lapse variations in the vertical and horizontal P- and S-wave velocities are measured along different directions through a sample undergoing hydrostatic deformation (e.g., Winkler and Liu, 1996; Prioul et al., 2004), or nonhydrostatic deformation (e.g., Winkler and Liu, 1996; Sarkar et al., 2003). Furthermore, the third-order coefficients estimated from hydrostatic deformation have been validated by a comparison of predicted and measured velocity for nonhydrostatic deformation (Prioul et al., 2004). Table 1 shows a selection of laboratory results in which the R-factors have been calculated from several published experiments. These indicate that R1 and R2 are predominantly positive, whereas R3 may be positive or negative. MacBeth et al. (2018) observe that the absolute values of R-factors from laboratory measurements are much higher than the field values. Bakk et al. (2020) and MacBeth et al. (2018) suggest reasons for why this mismatch exists. For example, the data considered by Prioul et al. (2004) were obtained under constant (drained) pore pressure and were measured for predominantly hydrostatic stress paths. This would clearly not be the case in an actual field depletion scenario; therefore, any third-order elasticity (TOE) values derived from such experiments must be considered biased. Here, we decide to use field values rather than laboratory values because the field values are more consistent with the time shifts that we observe from seismic data. Note also that the original R-factor of Hatchell and Bourne (2005a) and Røste et al. (2005) is defined in the uniaxially deforming reservoir as Runi=R1res and in the overburden (with a net volumetric strain of zero) as Rzvs=R1obR2ob, where the superscript res indicates the reservoir properties and ob indicates an overburden property. These factors are originally defined only for the vertical P-wave velocity and vertical strain and the framework introduced in our study may be regarded as a generalization of this concept. Substituting the appropriate cij’s from equations 3 and 4 into the original definitions for the Thomsen parameters (Thomsen, 1986) and retaining only first-order terms in strain, for the reservoir sandstone, we find the elliptical anisotropy (Helbig, 1983) correction
Δδ=Δε=(R1resR2res)ε33res.
(8)
For wave propagation in the [x1, x3] plane of the overburden, we find
Δε=(R1obR2ob)(ε33obε11ob)Δδ=(R1obR2ob)(2ε33ob+ε22ob),
(9)
whereas for wave propagation in the [x2, x3] plane, we find
Δε=(R1obR2ob)(ε33obε22ob)Δδ=(R1obR2ob)(2ε33ob+ε11ob).
(10)
In the above, the difference R1R2 is fundamentally the Hatchell–Bourne–Roste (HBR) R-factor (for zero volumetric strain); thus, it may be connected to another set of fundamental constants, the Thomsen parameters. The overburden equations can be further reduced to functions of vertical strain by applying the zero volumetric strain constraint described above. In the case of equal horizontal strains, both vertical symmetry planes give Δδ=Δε=3/2(R1obR2ob)ε33ob. To return to the objective of our study, we consider 4D AVO for two main situations: with and without time-lapse changes in the shale. To retain consistency with previous studies, we consider only two terms for the AVO, the zero-offset intercept and the gradient.

Geomechanical changes only in reservoir rocks

To draw a comparison between the rock-physics-focused approach of Landrø (2001) and our approach, we first consider the amplitude variation with incidence angle for the case of no time-lapse changes in the overburden shales, and for the intercept and gradient AVO terms only. This particular situation is closest to the case studied by Landrø (2001). We will find that even here there is an error because, according to geomechanical principles, the production-induced perturbations in the elastic moduli have VTI symmetry. Applying the same reasoning as in the previous section to equation 2, the time-lapse changes for the zero-offset (or intercept) term, ΔI, are written
ΔI12(Δρρ+Δαα)res.
(11)
To compute the relevant time-lapse perturbations, we recognize that the vertical P-wave velocity for the initial preproduction case is given by αi=c330/ρ. Similarly, the vertical P-wave velocity postproduction is αf=(c330(12R1resε33res))/ραi(1R1resε33res). Thus, the relative velocity change to first order in time-lapse perturbation is the familiar (Δα/α)res=R1resε33res. The density terms have already been discussed previously (Δρ/ρ)res=ε33res. Substituting these expressions into equation 11 yields
ΔI=12(1+R1res)ε33res.
(12)
By a similar development, it can be shown that the AVO gradient term ΔG in the [x1, x3] and [x2, x3] planes is identical (because this is a VTI medium); thus,
ΔG12{Δαα4(βα)2Δμμ+Δδ}res.
(13)
Substituting the proposed strain model, we find
ΔG=12{R1res+2(R1resR2res)+(R1resR2res)}ε33res
(14a)
or
ΔG=12(2R1res3R2res)ε33res.
(14b)
The first two terms in equation 14aR1res and +2(R1resR2res) may be equated to the isotropic terms of the original 4D AVO equation. However, the third term (R1resR2res) is the desired time-lapse anisotropic correction. The error in estimating the gradient and failing to take anisotropy into account is therefore (R1resR2res)/{R1res+2(R1resR2res)}.
In a reservoir for which the uniaxial assumption is valid (lateral extent of the reservoir much larger than its thickness, homogeneous initially isotropic rock, and linear poroelasticity), the vertical strain ε33res resulting from a pore pressure change ΔP is (Fjaer et al., 2008)
ε33res=Δhh=CuniΔP,
(15)
where Cuni is the uniaxial pore volume compressibility and Δh is the reduction in reservoir thickness h. Thus, the pressure decrease leads to a negative (or compressive) strain. The uniaxial compressibility coefficient is defined as the fractional change in volume with the pore pressure change when the lateral strain is constrained to be zero (for further details and values of this quantity, see Appendix  A). Thus, equations 12 and 14 may be easily rewritten as a function of reservoir pressure change replacing ε33res by CuniΔP. The equations predict that the intercept and gradient terms depend linearly on the reservoir pore pressure change — as is expected from a first-order perturbation theory. As a cross-reference with Landrø (2001), we quote his equation for the first-order pressure dependence and the approximation tan2θsin2θ
ΔR=12lαΔP+12{lα8(βα)2lβ}ΔPsin2θ,
(16)
where the fractional change in velocity per unit pressure for pressure-only changes is given as lα=(1/ΔP)Δα/α and lβ=(1/ΔP)Δβ/β. Equation 16 is directly comparable with our versions in equations 12 and 14, if the anisotropic corrections are excluded. We note that these equations are identical in form, but with differently labeled controlling factors. We can use this similarity to link the empirically derived lα and lβ values to the geomechanical parameters R1 and R2 and the uniaxial compressibility Cuni. If we ignore the density term in the intercept and the anisotropic term in the gradient,
lα=R1resCuni,
(17a)
lβ=14(αβ)2(R1resR2res)Cuni.
(17b)
These equations may also be inferred from the TOE equation 11. Interestingly, for an α/β of 2, equation 17b becomes lβ=(R1resR2res)Cuni. If we apply the approximation R1>R2, which appears to be the case from the laboratory (Prioul et al., 2004), then lαlβ. Equation 17 provides a useful quality check between R1res and R2res parameters (or the two HBR R-factors) and lα and lβ. To calculate one set of parameters from the other, we must know the uniaxial compressibility. The term Cuni can be measured in the laboratory (Fjaer et al., 2008); however, this practice is less common than laboratory measurement of pore volume compressibility. Thus, if required, Cuni can be determined directly from the pore volume compressibility Cpp (MacBeth et al., 2018). The Cpp values lie in the range of 0.22.9×104  MPa1 (1.520×106  psi1) for consolidated sandstones and 0.812.3×104  MPa1 (5.585×106  psi1) for unconsolidated sandstones. As a guide to computation for a given porosity, Hall (1953) empirical regression to laboratory data has been used for many years within many commercial softwares without question (Jalalh, 2006). Although some other regressions have been published, the results are largely dependent on the quality and quantity of the laboratory data available. Because the improvements are subject to ongoing discussion, we opt to use Hall (1953) to generate the compressibilities: Cpp=0.259/ϕ0.4358×104  MPa1 (or Cpp=1.782/ϕ0.4358×105  psi1). Table 2 provides the results of converting lα and lβ to R1res and R2res for the Gullfaks field and North Sea fields that correspond to three different sandstones from the pressure sensitivity database of MacBeth (2004). To obtain the estimates, we have also assumed α/β=2. The estimated R1 values are an order of magnitude higher than those recorded in the field (R12 or 3), but they are lower than the values measured in laboratory studies, which are several orders of magnitude higher than the field values (see Table 1).

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 (R1resR2res)/{R1res+2(R1resR2res)} 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

This case has not been studied by previous authors and reveals an additional term due to the overburden geomechanics. The overburden shale will strain in response to induced stresses caused by pore pressure changes in the reservoir. According to equation 5, an originally VTI shale will assume an orthorhombic symmetry after these changes. Following Rüger (1998), to analyze this efficiently, we consider the two orthogonal planes [x1, x3] and [x2, x3], where the wave behavior is effectively VTI. Applying equation 1a to either vertical plane of symmetry, the time-lapse changes for the zero offset or intercept term ΔI are now written as
ΔI12(Δρρ+Δαα)res12(Δρρ+Δαα)ob.
(18)
For the overburden, the vertical P-wave velocity for the preproduction case is given by αi=c330/ρ. In addition, in the overburden, the postproduction vertical P-wave velocity is αf=(c330(12(R1obR2ob)ε33))/ραi{1(R1obR2ob)ε33}. Thus, the relative velocity changes for the overburden are to first order in the time-lapse perturbation (Δα/α)ob=(R1obR2ob)ε33ob. However, unlike the reservoir, the density term is (Δρ/ρ)ob=0. Substituting these and the previous reservoir expression into equation 18 yields
ΔI=12[(R1obR2ob)ε33ob(1+R1res)ε33res].
(19)
For the gradient term, we find
ΔG12{Δαα4(βα)2Δμμ+Δδ}res12{Δαα4(βα)2Δμμ+Δδ}ob,
(20)
where δ for the overburden contribution is now a function of the azimuth. When the TOE strain model is applied and we focus on contributions in the two vertical planes of symmetry, this gives the additional overburden contributions for the [x1, x3] plane
ΔG=12{(2R1res3R2res)ε33res}12{(R1obR2ob)ε33ob[(R1obR2ob)4(R2obR3ob)]ε22ob}
(21)
with a similar result found for the [x2, x3] plane by substituting ε11 by ε22. The anisotropy term Δδ now also brings the overburden horizontal strain term into this equation.
For significant deformation in the overburden shale, the intercept term is dependent on the R-factors of the reservoir and overburden layers, in addition to their vertical strains. Therefore, the overburden term contributes an additional error to the time-lapse intercept. We observe that this error in equation 19 does not depend on the horizontal strain, but only on the vertical strain — as one would expect. Interestingly, because the reservoir is under uniaxial compression when the overburden is in extension, and vice versa, reservoir and overburden strain polarities are in opposite directions, i.e., |ε33res|=|ε33ob|, and thus the reservoir and overburden contributions reinforce each other
|ΔI|=|(R1obR2ob)ε33ob|+|(1+R1res)ε33res|.
(22)
It should be noted that the mean strain in the overburden is typically much smaller by more than an order of magnitude than that of the reservoir; thus, it is likely that the ε33res term will mostly dominate because the R-factor combinations are not significantly different: R1obR2ob5 and R1res2 (Hatchell and Bourne, 2005a). Using these R-factor values, the relative error between the shale effect and the no-shale effect cases amounts to no more than 5/3(ε33ob/ε33res). The more stress sensitive the reservoir is relative to the overburden, the smaller this factor and the smaller the dependence on the strain ratio.

A rough estimate of ε33ob/ε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 ε33ob/ε33res=0.5×(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 ε33ob/ε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 1. For the extremes of this phenomenon, the error in the intercept is now observed to be sizable and up to 300%.

The overburden term that appears in the AVO gradient of equation 21 is also dependent on the vertical strain, but it also depends on the horizontal principal strain in the overburden. Under the assumption of zero volumetric strain and equal principal horizontal strains, equation 21 becomes
ΔG=12(2R1res3R2res)ε33res14{3(R1obR2ob)4(R2obR3ob)}ε33ob.
(23)
If the horizontal strains are unequal, such that ε22/ε11=x, then we may substitute ε22 in equation 21 by xε33/(1+x)
ΔG=12{(2R1res3R2res)ε33res}+12{(R1obR2ob)[(R1obR2ob)4(R2obR3ob)]x1+x}ε33ob.
(24)
The overburden error term is more complicated than for the intercept, and it requires calculation using measured R-values. Computation (see later) shows the error as a function of the ratio of vertical strains ε33ob/ε33res for equal horizontal strain. It is deduced that errors can reach up to 155%.

Discussion of the approximations

Here, we compare the results above, the consequences for 4D AVO, and in particular the computation of pressure. To do this, we first focus on the magnitudes of the three R-factors R1, R2, and R3. Values of these factors can be readily derived from studies in which the third-order coefficients have been measured in the laboratory. This typically requires measuring the P- and S-wave velocities in directions 0o, 45o, and 90o 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||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
Δαα=Rε33.
(25)
This relation is used to evaluate, measure, and model time-lapse time shifts, for which the agreement is generally excellent (MacBeth et al., 2019). We note, as did Herwanger (2008) and Rodriguez-Herrera et al. (2015), that for uniaxial strain conditions as experienced in the reservoir and a strain model based on the TOE coefficients as described above, the relevant R-factor is Runi=R1 (where the superscript uni denotes the uniaxial strain). In addition, in the overburden, the zero volumetric strain assumption (ε11+ε22+ε33=0) holds; thus, the relevant R-factor is Rzvs=R1R2 (where the superscript zvs denotes the zero volumetric strain condition). It is noted that the AVO equations that we derive contain R-factors in combinations of R1, R1-R2, or R2-R3. In other words, two out of the three combinations of field values are already measured. Neglecting R3 as the smallest of the R-factors and assuming lithology dependence allows us to make a simple link with the seismic time-shift analysis: R1resR1ob2 and R1resR2resR1obR2ob5 (using the two most commonly quoted values for the R-factor). It should be noted that there is some variation from these widely quoted values as noted by MacBeth et al. (2018); thus, these approximations must be considered as generic values. Let us now rewrite the four approximations for the 4D AVO to recognize the individual contributing terms
ΔR=12Runiε33res+12{Runi+2Rzvs}ε33ressin2θ,
(26a)
ΔR=12(1+Runi)ε33res+12{Runi+2Rzvs}ε33ressin2θ,
(26b)
ΔR=12(1+Runi)ε33res+12{Runi+3Rzvs}ε33ressin2θ,
(26c)
ΔR=12{Rzvsε33ob(1+Runi)ε33res}+12{(Runi+3Rzvs)ε33res32Rzvsε33ob}sin2θ.
(26d)
Equation 26a is the strain model equivalent of Landrø’s (2001) original equation. Defining the compressive strain to be negative, the 4D AVO is observed to have a positive intercept and a negative gradient. Equation 26b includes the density correction that we can immediately see is not insignificant compared with an R-factor of 2, whereas equation 26c includes the corrections due to density and to anisotropy. Finally, equation 26d is the 4D AVO for all corrections including the overburden term. Substituting the field values proposed by Hatchell and Bourne (2005a), these equations now simplify further
ΔR=ε33res(1+4sin2θ),
(27a)
ΔR=ε33res(32+4sin2θ),
(27b)
ΔR=12ε33res(3+13sin2θ),
(27c)
ΔR=12{5ε33ob3ε33res}+12{13ε33res152ε33ob}sin2θ.
(27d)
Because the overburden strain is likely to be equal to or, most usually, less than the reservoir strain, then the contribution of pressure on the 4D AVO in all cases decreases with offset. As an example, for equation 26a, a typical reservoir strain of 5×103 (corresponding to a pressure change of 10 MPa and a uniaxial compressibility Cuni of 5×104  MPa1) gives an intercept of 0.005 and a gradient of 0.02. Figure 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., ε33ob=±|ε33res|, but we also calculate the intercept and gradient for more typical values of 104. 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.

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.

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 δ 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.

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

UNIAXIAL COMPRESSIBILITY VERSUS HYDROSTATIC COMPRESSIBILITY

According to Zimmerman (2017), in a porous rock, the uniaxial pore volume compressibility Cuni is related to the hydrostatic pore compressibility Cpp (defined as the fractional change in volume with pore pressure change) via the relation
Cuni(12(12ν)α3(1ν))Cpp,
(A-1)
where ν is the Poisson’s ratio and α is the Biot-Willis coefficient. The approximation is valid when the compressibility of the pore space is much larger than the grains constituting the rock matrix. Under the often-made assumption of incompressible grains α = 1, the term in brackets becomes (1 + ν)/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

We have converted hydrostatic measurements in the laboratory as published by MacBeth (2004) into equivalent uniaxial measurements for an identical drop in pressure (or effective stress) ΔP. This can be achieved by first considering equation 1, which describes the change of elastic modulus responsible for vertical P-wave velocity α=c33/ρ
c33=c330+c111ε33+c112(ε11+ε22).
(B-1)
For the case of uniaxial strain ε11=ε22=0,ε330 and for a sample subject to a stress path that creates uniaxial loading, this leads to a velocity change along the axial (three-axis) direction given by
(Δαα)uni=12c111c330ε33
(B-2a)
or
(Δαα)uniR1ε33.
(B-2b)
For uniaxial strain due to pressure depletion, the vertical strain is related to the particular pressure change ΔP via ε33=Cuni.ΔP. From equation B-1, if we now consider hydrostatic loading, for which all three strains are assumed to be equal (ε11=ε22=ε33=ε) such that the total volumetric strain εvol is 3ε, we have the velocity perturbation
(Δαα)iso12c111c330ε12c112c330(ε+ε)
(B-3a)
or
(Δαα)isoR1ε2R2ε,
(B-3b)
where in this case the volumetric strain is given by εvol=3ε=Cpp.ΔP. Now as the ratio of uniaxial strain to volumetric strain ε33/εvol=Cuni/Cpp, we can substitute back into equation B-2b to find the equivalent uniaxial velocity change for any given hydrostatic change
(Δαα)uniR1CuniCpp(3ε)
(B-4a)
or
(Δαα)uni3R1R1+2R2CuniCpp(Δαα)iso,
(B-4b)
where the ratio of compressibilities is a function of the Poisson’s ratio as defined in Appendix  A. Equation B-4b provides the relationship between the hydrostatic measurement for a given drop in pressure and the uniaxial measurement equivalent to the same drop in pressure. Although we have a lack of measurements for R1 and R2 at the field scale, we assume that the ratio of R-values from the laboratory is nevertheless accurate. In Table 1, we calculate the ratio 3R1/(R1+2R2) to be 1.5 for Colton Sandstone. In addition, we observe that for a Poisson’s ratio of 1/4 that Cuni=2/3Cpp. With these values, we thus conclude that
(Δαα)uni(Δαα)iso
(B-5)
for an identical change of pressure ΔP.

CURVATURE TERM IN THE 4D AVO

In the main text, we chose to consider only the intercept and gradient terms of the 4D AVO to remain consistent with previous treatments of pressure-saturation estimation. Here, for completeness, we also provide the curvature term H. It is given by the following:
ΔH12(Δαα¯+Δε)res12(Δαα¯+Δε)ob.
(C-1)
This term may be calculated in the same way as before, by substituting the time lapse into this equation the perturbations that are a function of the R-factors and principal strains. After some rearrangement, we obtain
ΔH12{R2resε33res(R1obR2ob)ε22ob}.
(C-2)
With a slight rewrite to emphasize the empirical R-factors, and for equal principal horizontal strains, this reduces to
ΔH12{(RzvsRuni)ε33res+12Rzvsε33ob}.
(C-3)

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.