Effects of induced stress on seismic forward modelling and inversion

Effects of induced stress on seismic forward modelling and inversion Summary We demonstrate how effects of induced stress may be incorporated in seismic modelling and inversion. Our approach is motivated by the accommodation of pre-stress in global seismology. Induced stress modifies both the equation of motion and the constitutive relationship. The theory predicts that induced pressure linearly affects the unstressed isotropic moduli with a slope determined by their adiabatic pressure derivatives. The induced deviatoric stress produces anisotropic compressional and shear wave speeds; the latter result in shear wave splitting. For forward modelling purposes, we determine the weak form of the equation of motion under induced stress. In the context of the inverse problem, we determine induced stress sensitivity kernels, which may be used for adjoint tomography. The theory is illustrated by considering 2-D propagation of SH waves and related Fréchet derivatives based on a spectral-element method. Elasticity and anelasticity, Equations of state, High-pressure behaviour, Computational seismology, Theoretical seismology, Wave propagation 1 INTRODUCTION In many engineering applications, it is important to be able to monitor stress changes during operations. Obvious examples are drilling, oil and gas exploration, geothermal activity, and mining. In geomechanical modelling, which describes the response of the underground to mechanical perturbations, knowledge of the state of stress is essential. Most often estimates of stress are obtained from seismic data using simplifying assumptions or empirical models. These empirical models and assumptions are derived from laboratory rock physics experiments. We propose a method which quantitatively links changes in seismic wavefields directly to changes in stress using first principles. Basic effects of changes in stress on seismic wave speeds have been known for a long time (e.g. Birch 1961; Nur & Simmons 1969), and such effects have been observed in laboratory (e.g. Eberhart-Phillips et al.1989; Verdon et al.2008) and field studies (e.g. Fazio et al.1973; Silver et al.2007). Exploration for oil & gas induces stress changes inside and outside reservoirs (e.g. Hatchell & Bourne 2005). As a consequence, seismic waves are observed to have different traveltimes before and after depletion of the reservoir or after waste water injection. Such time-lapse time shifts are caused by two competing factors, namely, changes in layer geometry and stress-induced changes in seismic wave speeds. Two descriptions of the effects of an induced stress on seismic wave speeds have been commonly used. The first approach focuses on stress effects on pre-existing or induced cracks, which manifest themselves in the form of seismic anisotropy (e.g. Nur 1971; O’Connell & Budiansky 1974; Zheng 2000). O’Connell & Budiansky (1974) incorrectly accounted for the crack energy as the crack density increases, an issue addressed by Bruner (1976) and Henyey & Pomphrey (1982). The second approach is based on third-order elasticity theory (Murnaghan 1951; Hughes & Kelly 1953; Egle & Bray 1976), which requires the knowledge of higher-order elastic constants. The latter are not easily measured in the laboratory, with many discrepancies remaining, although new methods are continuously proposed (e.g. Renaud et al.2012). In this paper, we take an approach motivated by the accommodation of pre-stress in global seismology. As first discussed in Dahlen (1972a,b) and also in Dahlen & Tromp (1998, sections 3.3.2 and 3.6.2), pre-stress affects both the equation of motion and the constitutive relationship. Here we explore these effects both from a forward modelling perspective and from the point of view of the inverse problem. 2 INITIAL EQUATION OF MOTION Prior to inducing stress, the equation of motion is   \begin{equation} \rho \,\partial _t^2{{\bf s}}-\boldsymbol \nabla \cdot {{\bf T}}={{\bf 0}}, \end{equation} (1)where ρ denotes mass density, s the displacement vector, and T the symmetric stress tensor. At internal discontinuities and the free surface, the traction needs to be continuous, that is,   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}]_-^+={{\bf 0}}, \end{equation} (2)where $$\hat{{{\bf n}}}$$ denotes the unit normal pointing from the − side of a discontinuity to the + side, and where the notation $$[\,\cdot \,]_-^+$$ denotes the jump in the enclosed quantity when going from the − to the + side. On a traction-free surface, such as Earth’s surface, traction vanishes: $$\hat{{{\bf n}}}\cdot {{\bf T}}={{\bf 0}}$$. The stress tensor, T, is related to the symmetric infinitesimal strain tensor,   \begin{eqnarray} \boldsymbol \epsilon =&\frac{1}{2}[\boldsymbol \nabla {{\bf s}}+(\boldsymbol \nabla {{\bf s}})^T)] , \end{eqnarray} (3)(a superscript T denotes the transpose) via Hooke’s law:   \begin{equation} {{\bf T}}= \boldsymbol {\Gamma} :\boldsymbol {\epsilon} . \end{equation} (4)The fourth-order elastic tensor, $$\boldsymbol \Gamma$$, linearly relates the stress and strain tensors. Let U denote the elastic energy per unit mass. Then the stress may be related to the elastic energy density, defined to second order in strain, as   \begin{eqnarray} \rho \,U=&\frac{1}{2}\,\boldsymbol \epsilon :\boldsymbol \Gamma :\boldsymbol \epsilon\!, \end{eqnarray} (5)via   \begin{equation} {{\bf T}}=\rho \,\frac{\partial U}{\partial {\boldsymbol \epsilon} } . \end{equation} (6)The symmetries of the stress and strain tensors and the quadratic nature of the elastic energy density dictate that the fourth-order elastic tensor, with elements   \begin{equation} \Gamma _{ijk\ell}=\frac{\partial ^2 U}{\partial \epsilon _{ij}\,\partial \epsilon _{k\ell }}, \end{equation} (7)exhibits the symmetries Γijkℓ = Γjikℓ, Γijkℓ = Γijℓk, and Γijkℓ = Γkℓij, which, in the most general case, reduces the number of independent parameters from 81 to 21. It is often convenient to express the elastic tensor in terms of its isotropic and purely anisotropic parts as   \begin{eqnarray} \Gamma _{ijk\ell }=&\left(\kappa -\frac{2}{3}\,\mu \right)\,\delta _{ij}\,\delta _{k\ell }+\mu \,(\delta _{ik}\,\delta _{j\ell }+\delta _{i\ell }\,\delta _{jk})+\gamma _{ijk\ell } , \end{eqnarray} (8)where κ and μ denote the isotropic bulk and shear moduli, respectively, and γijkℓ the purely anisotropic contribution. The elements γijkℓ exhibit the same symmetries as the elements Γijkℓ, and for purely isotropic media γijkℓ = 0. 3 INDUCED STRESS We induce a symmetric static stress field, T0, subject to the equilibrium condition   \begin{equation} {\boldsymbol \nabla} \cdot { {\bf T}}^0 ={{\bf 0}}. \end{equation} (9)The associated boundary condition is   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}^0]_-^+={{\bf 0}}. \end{equation} (10)The induced stress may be written in terms of an induced pressure,   \begin{eqnarray} p^0=&-\frac{1}{3}\,{\rm tr}({{\bf T}}^0) , \end{eqnarray} (11)and a symmetric trace-free induced deviatoric stress   \begin{eqnarray} {\boldsymbol \tau}^{0}=&{{\bf T}}^{0}-\frac{1}{3}\,{\rm tr}({{\bf T}}^{0})\,{{\bf I}}, \end{eqnarray} (12)such that,   \begin{equation} {{\bf T}}^0=-p^0\,{{\bf I}}+{\boldsymbol \tau} ^0 , \end{equation} (13)where I denotes the identity tensor. Taken together, the equilibrium condition and decomposition imply that   \begin{equation} {\boldsymbol \nabla} p^0={\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0 , \end{equation} (14)which enables us to eliminate $${\boldsymbol \nabla} p^0$$ in terms of $${\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0$$. For example, we have   \begin{equation} {\boldsymbol \nabla} {\bf T}^0=-{\boldsymbol \nabla} p^0\,{{\bf I}}+{\boldsymbol \nabla} {\boldsymbol \tau} ^0=-{\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0\,{{\bf I}}+{\boldsymbol \nabla} {\boldsymbol \tau} ^0 . \end{equation} (15)Note that the equilibrium condition (9) constrains only three of the six elements of the induced stress tensor. 4 EQUATION OF MOTION UNDER INDUCED STRESS Seismic waves propagating in a stress-induced medium are governed by the wave equation (Dahlen & Tromp 1998, eq. 3.56 without density and gravity perturbations or rotational terms)   \begin{equation} \rho \,\partial _t^2{{\bf s}}-{\boldsymbol \nabla} \cdot {{\bf T}}^{\rm L1}+{\boldsymbol \nabla} \cdot ({{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0)={{\bf 0}}, \end{equation} (16)where TL1 denotes the symmetric incremental Lagrangian Cauchy stress, and where the term $$-{{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0$$ captures advection of the induced stress. Using decomposition (13) and equilibrium eq. (9), the equation of motion (16) may be rewritten in the form (Dahlen & Tromp 1998, eq. 3.58 without density and gravity perturbations or rotational terms)   \begin{equation} \rho \,\partial _t^2{{\bf s}}-{\boldsymbol \nabla} \cdot {{\bf T}}^{\rm L1}-{\boldsymbol \nabla} [{{\bf s}}\cdot ({\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0)]+{\boldsymbol \nabla} \cdot ({{\bf s}}\cdot {\boldsymbol \nabla} {\boldsymbol \tau} ^0) ={{\bf 0}}. \end{equation} (17)The associated boundary condition is most conveniently expressed in terms of the asymmetric incremental first Piola–Kirchhoff stress, TPK1, namely (Dahlen & Tromp 1998, eq. 3.65),   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}^{\rm PK1}]_-^+={{\bf 0}}. \end{equation} (18)The incremental first Piola–Kirchhoff stress is related to the incremental Lagrangian Cauchy stress via (Dahlen & Tromp 1998, eq. 3.36)   \begin{equation} {{\bf T}}^{\rm PK1}={{\bf T}}^{\rm L1}+{{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0 . \end{equation} (19)For numerical simulations, in particular spectral-element simulations, the weak form of the equation of motion is discussed in Appendix. What we need next is a constitutive relationship for seismic wave propagation in a medium under induced stress. 5 ELASTIC STRAIN ENERGY DENSITY UNDER INDUCED STRESS A discussion of elastic strain energy requires careful definition of ‘the strain tensor’, and consideration of second-order quantities in the displacement gradient. As discussed in Dahlen & Tromp (1998, section 3.6.1), in an adiabatic, perfectly elastic material, the Lagrangian internal energy density depends only on the Lagrangian strain tensor   \begin{eqnarray} {{\bf E}}^{\rm L} =& \frac{1}{2}[{\boldsymbol \nabla} {{\bf s}}+({\boldsymbol \nabla} {{\bf s}})^T)]+\frac{1}{2}({\boldsymbol \nabla} {{\bf s}})\cdot ({\boldsymbol \nabla} {{\bf s}})^T , \end{eqnarray} (20)which may, alternatively, be expressed in terms of the infinitesimal strain tensor $${\boldsymbol \epsilon}$$ defined in eq. (3) and the antisymmetric infinitesimal vorticity tensor   \begin{eqnarray} {\boldsymbol \omega} &=-\frac{1}{2}\,[{\boldsymbol \nabla} {{\bf s}}-({\boldsymbol \nabla} {{\bf s}})^T] , \end{eqnarray} (21)as   \begin{eqnarray} {{\bf E}}^{\rm L} &= {\boldsymbol \epsilon} +\frac{1}{2}\,{\boldsymbol \epsilon} \cdot {\boldsymbol \omega} -\frac{1}{2}\,{\boldsymbol \omega} \cdot {\boldsymbol \epsilon} -\frac{1}{2}\,{\boldsymbol \omega} \cdot {\boldsymbol \omega} +\frac{1}{2}\,{\boldsymbol \epsilon} \cdot {\boldsymbol \epsilon} . \end{eqnarray} (22)The internal energy density under induced stress must be calculated correct to second order in the displacement gradient, $${\boldsymbol \nabla} {{\bf s}}$$, which prohibits the linearization $${{\bf E}}^{\rm L}\approx {\boldsymbol \epsilon}$$ in energetic considerations. Keeping this in mind, to second order in $${\boldsymbol \nabla} {{\bf s}}$$ the Lagrangian internal energy density may be expressed in the form   \begin{eqnarray} \rho ^0\,U^{\rm L}&={{\bf T}}^0:{{\bf E}}^{\rm L}+\frac{1}{2}\,{{\bf E}}^{\rm L}:{\boldsymbol \Xi} :{{\bf E}}^{\rm L} , \end{eqnarray} (23)where ρ0 denotes the density before straining the material, making UL the Lagrangian internal energy per unit mass. We have assumed for convenience that the Lagrangian internal energy density vanishes in the absence of strain. The symmetric second Piola-Kirchhoff stress is defined in terms of the Lagrangian internal energy via   \begin{eqnarray} {{\bf T}}^{\rm SK} \ &\equiv & \ \rho ^0\,\frac{\partial U^{\rm L}}{\partial {{\bf E}}^{\rm L}} \nonumber\\ & = & \ {{\bf T}}^0+{{\bf T}}^{\rm SK1} , \end{eqnarray} (24)where TSK1 denotes the symmetric incremental second Piola–Kirchhoff stress, namely,   \begin{eqnarray} {{\bf T}}^{\rm SK1} \ & = & \ {\boldsymbol \Xi} :{{\bf E}}^{\rm L} \nonumber\\ &\approx & \ {\boldsymbol \Xi} :{\boldsymbol \epsilon} , \end{eqnarray} (25)where the last equality holds to first order in the displacement gradient. The components of the fourth-order tensor $${\boldsymbol \Xi}$$ are given by   \begin{equation} \Xi _{ijk\ell }=\rho ^0\,\frac{\partial ^2 U^{\rm L}}{\partial E_{ij}^{\rm L}\,\partial E_{k\ell }^{\rm L}} , \end{equation} (26)implying that $${\boldsymbol \Xi}$$ exhibits the symmetries Ξijkℓ = Ξjikℓ, Ξijkℓ = Ξijℓk, and Ξijkℓ = Ξkℓij, thereby reducing its number of independent elements from 81 to 21. As discussed in Dahlen & Tromp (1998, section 3.6.2) and Dahlen (1972a), the tensor $${\boldsymbol \Xi}$$ is one of a family of permissible fourth-order tensors which may be used to define the Lagrangian internal energy density, and, therefore, the incremental second Piola–Kirchhoff stress. Specifically, we may write   \begin{equation} \Xi _{ijk\ell }=\Gamma _{ijk\ell } +a\,\left(T^0_{ij}\,\delta _{k\ell }+T^0_{k\ell }\,\delta _{ij}\right) +b\,\left(T^0_{ik}\,\delta _{j\ell }+T^0_{jk}\,\delta _{i\ell }+T^0_{i\ell }\,\delta _{jk}+T^0_{j\ell }\,\delta _{ik}\right) , \end{equation} (27)where the elements Γijkℓ are given by (8), and for arbitrary values of the quantities a and b, which are allowed to vary spatially. We now turn to the task of determining suitable values for these quantities. 5.1 Plane-wave solutions In this section we consider plane-wave solutions in a homogeneous medium with induced stress. We consider a plane wave of the form   \begin{equation} {{\bf s}}={{\bf a}}\,\exp ({{\bf k}}\cdot {{\bf x}}-{\omega}\,t) , \end{equation} (28)where a denotes the amplitude, k the wavevector, and ω the angular frequency. Substitution of plane wave (28) in the equation of motion (17) leads to the eigenvalue problem   \begin{equation} {{\bf B}}\cdot {{\bf a}}=c^2\,{{\bf a}}, \end{equation} (29)where c = ω/k denotes the phase speed, and where B denotes the symmetric Christoffel tensor with elements (Dahlen 1972b; Dahlen & Tromp 1998, eq. 3.156)   \begin{equation} \rho B_{j\ell }=\Lambda _{ijk\ell }\,\hat{k}_i\,\hat{k}_k . \end{equation} (30)The elements Λijkℓ are related to the elements Ξijkℓ via (Dahlen & Tromp 1998, eq. 3.122)   \begin{equation} \Lambda _{ijk\ell }=\Xi _{ijk\ell }+T^0_{ik}\,\delta _{j\ell } , \end{equation} (31)and are used to relate the first Piola-Kirchhoff stress defined in eq. (19) to the displacement gradient, that is,   \begin{equation} {{\bf T}}^{\rm PK1} = {\boldsymbol \Lambda} :{\boldsymbol \nabla} {{\bf s}}. \end{equation} (32)More specifically, using eqs (8), (13), (27), and (31), we have   \begin{eqnarray} \rho B_{j\ell } &=& \left(\kappa +{\textstyle\frac{1}{3}}\,\mu \right)\,\hat{k}_j\,\hat{k}_\ell +\mu \,\delta _{j\ell } +\gamma _{ijk\ell }\,\hat{k}_i\,\hat{k}_k -2\,(a+b)\,p^0\,\hat{k}_j\,\hat{k}_\ell -(1+2\,b)\,p^0\,\delta _{j\ell } + a\,\tau ^0_{ij}\,\hat{k}_i\,\hat{k}_\ell +a\,\tau ^0_{k\ell }\,\hat{k}_j\,\hat{k}_k +(1+b)\,\tau ^0_{ik}\,\hat{k}_i\,\hat{k}_k\,\delta _{j\ell }\nonumber\\ && +\,b\,\tau ^0_{jk}\,\hat{k}_\ell \,\hat{k}_k +b\,\tau ^0_{i\ell }\,\hat{k}_i\,\hat{k}_j +b\,\tau ^0_{j\ell } . \end{eqnarray} (33)Here $$\hat{k}_i$$ denotes an element of the unit wavevector, $$\hat{{{\bf k}}}={{\bf k}}/k$$. Let us define   \begin{equation} \mu ^{\prime }\equiv -1-2\,b , \end{equation} (34)  \begin{equation} \kappa ^{\prime }\equiv \frac{1}{3}-2\,a-\frac{4}{3}\,b . \end{equation} (35)Then we find that the Christoffel tensor has elements determined by   \begin{eqnarray} \rho B_{j\ell } \ & = & \ \left[(\kappa +\kappa ^{\prime }\,p^0)+{\textstyle\frac{1}{3}}\,(\mu +\mu ^{\prime }\,p^0)\right]\,\hat{k}_j\,\hat{k}_\ell +(\mu +\mu ^{\prime }\,p^0)\,\delta _{j\ell } +\gamma _{ijk\ell }\,\hat{k}_i\,\hat{k}_k - {\textstyle\frac{1}{2}}\,\left(\kappa ^{\prime }+{\textstyle\frac{1}{3}}\,\mu ^{\prime }\right)\,\left(\tau ^0_{ij}\,\hat{k}_i\,\hat{k}_\ell +\tau ^0_{k\ell }\,\hat{k}_j\,\hat{k}_k\right)\nonumber\\ && +\,{\textstyle\frac{1}{2}}(1-\mu ^{\prime })\,\tau ^0_{ik}\,\hat{k}_i\,\hat{k}_k\,\delta _{j\ell } -{\textstyle\frac{1}{2}}(1+\mu ^{\prime })\,\tau ^0_{j\ell } . \end{eqnarray} (36)The terms involving κ + κ΄ p0 and μ + μ΄ p0 enable us to identify κ΄ and μ΄ as the adiabatic pressure derivatives of the bulk and shear moduli, respectively (Nur & Simmons 1969). In fact, it is this desirable result that motivated definitions (34) & (35). We may regard κ and μ as spatially variable elastic moduli at zero pressure, whereas κ΄ and μ΄ represent their spatially variable adiabatic pressure derivatives. Because B is a symmetric positive-definite tensor, the eigenvalue problem (29) has three positive eigenvalues, c2, and associated orthogonal eigenvectors, a. To first-order in the induced stress, assuming the medium before inducing stress is isotropic, that is, γijkℓ = 0, the compressional wave speed is determined by   \begin{eqnarray} \rho \,c_P^2=&(\kappa +\kappa ^{\prime }\,p^0)+\frac{4}{3}\,(\mu +\mu ^{\prime }\,p^0)-(\kappa ^{\prime }+\frac{4}{3}\,\mu ^{\prime })\,\hat{{{\bf k}}}^0\cdot \,{\boldsymbol \tau} ^0\cdot \hat{{{\bf k}}}^0 , \end{eqnarray} (37)and the two shear wave speeds are determined by   \begin{eqnarray} \rho \,c^2_{{\rm S}_{1,2}}=&(\mu +\mu ^{\prime }\,p^0)+\frac{1}{2}\,(1-\mu ^{\prime })\,\hat{{{\bf k}}}^0\cdot \,{\boldsymbol \tau} ^0\cdot \hat{{{\bf k}}}^0-\frac{1}{2}\,(1+\mu ^{\prime })\,\hat{{{\bf a}}}^0_{1,2}\cdot {\boldsymbol \tau} ^0\cdot \hat{{{\bf a}}}^0_{1,2} . \end{eqnarray} (38)Here $$\hat{{{\bf k}}}^0$$ denotes the unit wavevector prior to inducing stress, and $$\hat{{{\bf a}}}^0_{1,2}$$ the unit shear wave polarization directions prior to inducing stress. Note that $$\hat{{{\bf k}}}^0\cdot \hat{{{\bf a}}}^0_{1,2}=0$$. We conclude that definitions (34) and (35) lead to a desirable dependence of the compressional and shear wave speeds (37) & (38) on induced pressure, p0. We observe that the induced deviatoric stress, $${\boldsymbol \tau} ^0$$, manifests itself as an anisotropy in seismic wave speeds. In global seismology, the quantities a and b are selected such that the seismic wave speeds (37) and (38) are independent of the hydrostatic pre-stress, p0. This amounts to setting κ΄ = 0 and μ΄ = 0, and corresponds to the choice a = −b = 1/2. With this choice, the wave speeds (37) and (38) are in agreement with those of Dahlen (1972b, eqs 22 and 25). Note, in particular, that with this choice the compressional wave speed is independent of the induced stress. In this paper, we argue that, rather than choosing the parameters a and b, the parameters μ΄ and κ΄ should be recognized as the derivatives of the moduli with regards to hydrostatic pressure. With this interpretation, these parameters are not chosen, but rather determined by the material properties. For an S wave travelling in the z-direction there are two shear wave speeds, namely, for the wave polarized in the $$\hat{{{\bf x}}}$$ direction   \begin{equation} c^2_{{\rm S}_1}=\beta ^2+{\textstyle\frac{1}{2}}\,[(1-\mu ^{\prime })\,\tau ^0_{zz}-(1+\mu ^{\prime })\tau ^0_{xx}]/\rho , \end{equation} (39)and for the wave polarized in the $$\hat{{{\bf y}}}$$ direction   \begin{equation} c^2_{{\rm S}_2}=\beta ^2+{\textstyle\frac{1}{2}}\,[(1-\mu ^{\prime })\,\tau ^0_{zz}-(1+\mu ^{\prime })\tau ^0_{yy}]/\rho , \end{equation} (40)where $$\beta =\sqrt{\mu /\rho }$$ denotes the unperturbed shear wave speed. These polarization-dependent differences in shear wave speed lead to shear wave splitting, such that the split time between the arrival polarized in the x direction and the arrival polarized in the y direction is determined by   \begin{eqnarray} \delta T \ & = & \ \frac{1}{4}\int (1+\mu ^{\prime })(\tau ^0_{xx}-\tau ^0_{yy})/(\rho \,\beta ^3)\,{\mathrm{d}}z \nonumber\\ & = & \ \frac{1}{4}\int (1+\mu ^{\prime })(\tau ^0_{xx}/\mu -\tau ^0_{yy}/\mu )\,{\mathrm{d}}z/\beta . \end{eqnarray} (41)This expression provides a direct relationship between the split time and the induced deviatoric stress. 5.2 Constitutive relationship Definitions (34) and (35), which lead to a desirable dependence of the compressional and shear wave speeds on induced pressure via (37) and (38), provide the sought values of the quantities a and b, namely,   \begin{equation} b=-\frac{1}{2}\,(1+\mu ^{\prime }) , \end{equation} (42)  \begin{equation} a=\frac{1}{2}\,\left(1-\kappa ^{\prime }+\frac{2}{3}\,\mu ^{\prime }\right) . \end{equation} (43)This implies that eq. (27) becomes   \begin{equation} \Xi _{ijk\ell }=\Gamma _{ijk\ell } +{\textstyle\frac{1}{2}}\,\left(1-\kappa ^{\prime }+{\textstyle\frac{2}{3}}\,\mu ^{\prime }\right)\,\left(T^0_{ij}\,\delta _{k\ell }+T^0_{k\ell }\,\delta _{ij}\right) -{\textstyle\frac{1}{2}}\,(1+\mu ^{\prime })\,\left(T^0_{ik}\,\delta _{j\ell }+T^0_{jk}\,\delta _{i\ell }+T^0_{i\ell }\,\delta _{jk}+T^0_{j\ell }\,\delta _{ik}\right) . \end{equation} (44) As discussed in Section 4, seismic wave propagation under induced stress is described in terms of the incremental Lagrangian Cauchy stress, TL1, which is related to the incremental second Piola–Kirchhoff stress, TSK1, via (Dahlen & Tromp 1998, eqs 3.37)   \begin{equation} {{\bf T}}^{\rm L1}={{\bf T}}^{\rm SK1}-{{\bf T}}^0\,({\boldsymbol \nabla} \cdot {{\bf s}})+({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0+{{\bf T}}^0\cdot {\boldsymbol \nabla} {{\bf s}}. \end{equation} (45)Upon using eqs (25) and (44) in eq. (45), we find that the incremental Lagrangian Cauchy stress has components   \begin{eqnarray} T^{\rm L1}_{ij} & = & \left[(\kappa +\kappa ^{\prime }\,p^0)-{\textstyle\frac{2}{3}}\,(\mu +\mu ^{\prime }\,p^0)\right]\,\delta _{ij}\,\epsilon _{kk} +2\,(\mu +\mu ^{\prime }\,p^0)\,\epsilon _{ij} +\gamma _{ijk\ell }\,\epsilon _{k\ell } + {\textstyle\frac{1}{2}}\,\left(1-\kappa ^{\prime }+{\textstyle\frac{2}{3}}\,\mu ^{\prime }\right)\,\tau ^0_{k\ell }\,\epsilon _{k\ell }\,\delta _{ij} -{\textstyle\frac{1}{2}}\,\left(1+\kappa ^{\prime }-{\textstyle\frac{2}{3}}\,\mu ^{\prime }\right)\,\tau ^0_{ij}\,\epsilon _{kk} \nonumber\\ &&-\,\mu ^{\prime }\,\left(\tau ^0_{ik}\,\epsilon _{kj} +\tau ^0_{jk}\,\epsilon _{ki}\right) +\omega _{ik}\,\tau ^0_{kj}-\tau ^0_{ik}\,\omega _{kj} . \end{eqnarray} (46)In the absence of induced deviatoric stress, $$\tau ^0_{ij}=0$$, we have   \begin{equation} T^{\rm L1}_{ij} = \left[(\kappa +\kappa ^{\prime }\,p^0)-{\textstyle\frac{2}{3}}\,(\mu +\mu ^{\prime }\,p^0)\right]\,\delta _{ij}\,\epsilon _{kk} +2\,(\mu +\mu ^{\prime }\,p^0)\,\epsilon _{ij} +\gamma _{ijk\ell }\,\epsilon _{k\ell } , \end{equation} (47)and here again we clearly recognize the contributions κ + κ΄ p0 and μ + μ΄ p0 as capturing the pressure dependence of the isotropic moduli. In terms of the stress tensor prior to inducing stress, given by eq. (4), we may write   \begin{equation} {{\bf T}}^{\rm L1} = {{\bf T}}+\Delta {{\bf T}}, \end{equation} (48)where   \begin{eqnarray} \Delta {{\bf T}} &= & \overbrace{ \left[\left(\kappa ^{\prime }-\textstyle{\frac{2}{3}}\,\mu ^{\prime }\right)\,{\rm tr}({\boldsymbol \epsilon} )\,{{\bf I}}+2\,\mu ^{\prime }\,{\boldsymbol \epsilon}\right ]\,p^0 }^{\rm {{strain-dependent\, stress\, (due\, to\, induced\, pressure)}}} + \underbrace{ \textstyle{\frac{1}{2}}\,[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}-{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0] -\textstyle{\frac{1}{2}}\,\left(\kappa ^{\prime }-\textstyle{\frac{2}{3}}\,\mu ^{\prime }\right)\,[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}+{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0] -\mu ^{\prime }\,({\boldsymbol \tau} ^0\cdot {\boldsymbol \epsilon} +{\boldsymbol \epsilon} \cdot {\boldsymbol \tau} ^0) }_{{\rm {strain-dependent\, stress\, (due\, to\, induced\, deviatoric\, stress)}}} \nonumber\\ &&+ \,\underbrace{ \boldsymbol \omega \cdot {\boldsymbol \tau} ^0 -{\boldsymbol \tau} ^0\cdot \boldsymbol \omega }_{\rm{rotated \,induced\, stress}}. \end{eqnarray} (49)We conclude that, even in an initially isotropic medium, induced stress manifests itself as an anisotropy in seismic wave speeds. Part of this involves a strain-dependent stress—with distinct contributions from the induced pressure and from the induced deviatoric stress—and another part involves a rotation of the induced stress. It is interesting to note that eq. (49) may be expressed also in terms of the full induced stress tensor, T0, as   \begin{eqnarray} \Delta {{\bf T}} &=& -\textstyle{\frac{1}{2}}\,\left(\kappa ^{\prime }-\textstyle{\frac{2}{3}}\,\mu ^{\prime }\right)\,[({{\bf T}}^0:{\boldsymbol \epsilon} )\,{{\bf I}}+{\rm tr}({\boldsymbol \epsilon} )\,{{\bf T}}^0] -\mu ^{\prime }\,({{\bf T}}^0\cdot {\boldsymbol \epsilon} +{\boldsymbol \epsilon} \cdot {{\bf T}}^0) + \textstyle{\frac{1}{2}}\,[({{\bf T}}^0:{\boldsymbol \epsilon} )\,{{\bf I}}-{\rm tr}({\boldsymbol \epsilon} )\,{{\bf T}}^0] + {\boldsymbol \omega} \cdot {{\bf T}}^0 -{{\bf T}}^0\cdot {\boldsymbol \omega} . \end{eqnarray} (50) We end this section by noting that in global seismology the quantities a and b are selected such that the incremental Lagrangian Cauchy stress (48) is independent of the hydrostatic pre-stress, p0, thereby also rendering the seismic wave speeds (37) and (38) independent of the hydrostatic pre-stress, as discussed previously. This amounts to setting κ΄ = 0 and μ΄ = 0 in eq. (49), implying via eqs (42) and (43) that a = −b = 1/2, and resulting in (Dahlen & Tromp 1998, eq. 3.145)   \begin{eqnarray} \Delta {{\bf T}}\ = & \ \frac{1}{2}\,[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}-{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0] +{\boldsymbol \omega} \cdot {\boldsymbol \tau} ^0 -{\boldsymbol \tau} ^0\cdot {\boldsymbol \omega} . \end{eqnarray} (51)Motivated by expressions (37) and (38) for the compressional and shear wave speeds, perhaps a better strategy is to select κ + κ΄ p0 and μ + μ΄ p0 as ‘the’ bulk and shear moduli in a pre-stressed Earth model, thereby building the pressure dependence into the moduli and capturing the effects of deviatoric pre-stress in terms of   \begin{eqnarray} \Delta {{\bf T}} = & \frac{1}{2}\,[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}-{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0] + {\boldsymbol \omega} \cdot {\boldsymbol \tau} ^0 -{\boldsymbol \tau} ^0\cdot {\boldsymbol \omega} -\frac{1}{2}\,\left(\kappa ^{\prime }-\frac{2}{3}\,\mu ^{\prime }\right)\,\left[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}+{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0\right] -\mu ^{\prime }\,({\boldsymbol \tau} ^0\cdot {\boldsymbol \epsilon} +{\boldsymbol \epsilon} \cdot {\boldsymbol \tau} ^0) . \end{eqnarray} (52)Compared to eq. (51), the additional terms are generally not negligible, because typical values of κ΄ and μ΄ are of order unity (e.g. Stacey 1992). Arguably, such terms could be captured by the purely anisotropic elements γijkℓ, but such an approach fails to reflect the full influence of the deviatoric pre-stress. 6 BORN APPROXIMATION To compare and contrast the equation of motion and boundary condition before inducing stress, eqs (1) and (2), with those after inducing stress, eqs (17) and (18), we use eq. (48) to rewrite the latter as   \begin{equation} \rho \,\partial _t^2{{\bf s}}-{\boldsymbol \nabla} \cdot {{\bf T}}={\boldsymbol \nabla} \cdot \Delta {{\bf T}}+{\boldsymbol \nabla} [{{\bf s}}\cdot ({\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0)]-{\boldsymbol \nabla} \cdot ({{\bf s}}\cdot {\boldsymbol \nabla} {\boldsymbol \tau} ^0) , \end{equation} (53)and   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}]_-^+=-[ \hat{{{\bf n}}}\cdot {{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-\hat{{{\bf n}}}\cdot ({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0 +\hat{{{\bf n}}}\cdot \Delta {{\bf T}}]_-^+ . \end{equation} (54)The terms on the right-hand sides of these equations capture effects due to induced stress. If the induced stress is accompanied by changes in density, δρ, and the elastic parameters, $$\delta {\boldsymbol \Gamma}$$, then the modified equations take the form   \begin{equation} \rho \,\partial _t^2{{\bf s}}-{\boldsymbol \nabla} \cdot {{\bf T}}= -\delta \rho \,\partial _t^2{{\bf s}}+{\boldsymbol \nabla} \cdot \delta {{\bf T}}+{\boldsymbol \nabla} \cdot \Delta {{\bf T}}+{\boldsymbol \nabla} [{{\bf s}}\cdot ({\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0)]-{\boldsymbol \nabla} \cdot ({{\bf s}}\cdot {\boldsymbol \nabla} {\boldsymbol \tau} ^0) , \end{equation} (55)and   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}]_-^+=-[ \hat{{{\bf n}}}\cdot \Delta {{\bf T}}+\hat{{{\bf n}}}\cdot {{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-\hat{{{\bf n}}}\cdot ({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0 +\hat{{{\bf n}}}\cdot \delta {{\bf T}}]_-^+ , \end{equation} (56)where   \begin{equation} \delta {{\bf T}}=\delta {\boldsymbol \Gamma} :{\boldsymbol \epsilon} . \end{equation} (57)Additional effects associated with changes in layer geometry are captured by perturbations in boundary locations. Such topographic perturbations on internal discontinuities and the free surface modify the boundary conditions. The theoretical implications are well known (see e.g. Dahlen 2005; Tromp et al.2005) and readily accommodated. For the sake of brevity, in this paper we focus mainly on the direct effects of induced stress, but in Section 7 we also summarize effects associated with perturbations in density, the elastic tensor, and boundary topography. If the wavefield prior to inducing stress, s, is governed by the equation of motion (1) subject to boundary conditions (2), and the wavefield after inducing stress, s + δs, is governed by the equation of motion (53) subject to boundary conditions (54), then, invoking the Born approximation, the perturbed wavefield, δs, due to the induced stress is determined by   \begin{eqnarray} \delta s_i({{\bf x}},t) &=& \ \int _{-\infty }^t\int _VG_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime })\,\partial ^{\prime }_k[ \Delta T_{kj}({{\bf x}}^{\prime },t^{\prime }) +s_\ell ({{\bf x}}^{\prime },t^{\prime })\,\partial ^{\prime }_m \tau _{m\ell }^0({{\bf x}}^{\prime },t^{\prime })\,\delta _{kj} -s_\ell ({{\bf x}}^{\prime },t^{\prime })\,\partial ^{\prime }_\ell \tau _{kj}^0({{\bf x}}^{\prime },t^{\prime }) ]\,{\mathrm{d}}^3{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime}\nonumber\\ && + \int _{-\infty }^t\int _\Sigma G_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime })\, \hat{n}_k({{\bf x}}^{\prime })\,[ \Delta T_{kj}({{\bf x}}^{\prime },t^{\prime }) +T_{kj}^0({{\bf x}}^{\prime },t^{\prime })\,\partial ^{\prime }_\ell s_\ell ({{\bf x}}^{\prime },t^{\prime }) -(\partial ^{\prime }_\ell s_k({{\bf x}}^{\prime },t^{\prime }))\,T_{\ell j}^0({{\bf x}}^{\prime },t^{\prime }) ]_-^+\,{\mathrm{d}}^2{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime}. \end{eqnarray} (58)Here V denotes the model volume, Σ the collection of all discontinuities, including internal discontinuities and the free surface, and Gij(x, x΄; t − t΄) the Green function for the medium prior to inducing stress. A partial derivative with respect to the primed coordinates, x΄, is denoted by $$\partial ^{\prime }_i$$. Additional terms associated with perturbations in density, δρ, elements of the elastic tensor, δΓijkℓ, and boundary topography, δh, may be found in Tromp et al. (2005, eqs 3 and 21). After a lengthy, tedious calculation, this result may be rewritten in the form   \begin{eqnarray} \delta s_i({{\bf x}},t) & = & -\int _{-\infty }^t\int _V\partial ^{\prime }_kG_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime })\, [\Delta T_{kj}({{\bf x}}^{\prime },t^{\prime }) -(\partial ^{\prime }_ms_\ell ({{\bf x}}^{\prime },t^{\prime }))\,\tau _{m\ell }^0({{\bf x}}^{\prime },t^{\prime })\,\delta _{kj} +(\partial ^{\prime }_\ell s_\ell ({{\bf x}}^{\prime },t^{\prime }))\,\tau _{kj}^0({{\bf x}}^{\prime },t^{\prime })]\,{\mathrm{d}}^3{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime } \nonumber\\ &&+ \,\int _{-\infty }^t\int _V(\partial ^{\prime }_m\partial ^{\prime }_jG_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime }))\, s_\ell ({{\bf x}}^{\prime },t^{\prime })\,\tau _{m\ell }^0({{\bf x}}^{\prime },t^{\prime }) \,{\mathrm{d}}^3{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime } -\int _{-\infty }^t\int _V(\partial ^{\prime }_\ell \partial ^{\prime }_kG_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime }))\, s_\ell ({{\bf x}}^{\prime },t^{\prime })\,\tau _{kj}^0({{\bf x}}^{\prime },t^{\prime }) \,{\mathrm{d}}^3{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime } . \end{eqnarray} (59)Note that the boundary integral has been eliminated, and, in particular, that eq. (59) contains no derivatives of the induced deviatoric stress. Result (59) forms the basis of the calculation of Fréchet derivatives in the next section. 7 KERNELS In this section we determine induced deviatoric stress sensitivity kernels or Fréchet derivatives. Suppose we record a set of three-component data with an array of Nr receivers. We denote the three component time series in this data set by di(xr, t), i = 1, 2, 3, r = 1, …, Nr. After a time lapse, ΔT, we record a second data set $$d_i^{\Delta T}({{\bf x}}_r,t)$$, i = 1, 2, 3, r = 1, …, Nr. We seek to find model perturbations such that the related simulated wavefield perturbations, δsi, capture the change in the data. This is accomplished by defining the misfit function   \begin{equation} \chi =\frac{1}{2}\,\sum _{r=1}^{N_r}\sum _{i=1}^3\int ||d_i({{\bf x}}_r,t)+\delta s_i({{\bf x}}_r,t)-d_i^{\Delta T}({{\bf x}}_r,t)||^2\,{\mathrm{d}}t , \end{equation} (60)such that   \begin{equation} \delta \chi =\sum _{r=1}^{N_r}\sum _{i=1}^3\int [d_i({{\bf x}}_r,t)-d_i^{\Delta T}({{\bf x}}_r,t)]\,\delta s_i({{\bf x}}_r,t)\,{\mathrm{d}}t . \end{equation} (61)Using the Born expression determined in Section 6, the derivative of the misfit function with respect to the induced deviatoric stress may be expressed as   \begin{equation} \delta \chi =\int _V p^0\,K_{p^0}+{\boldsymbol \tau} ^0:{{\bf K}}_{{\bf \tau }^0}\,{\mathrm{d}}V , \end{equation} (62)where   \begin{equation} K_{p^0}=-\int _0^T[\kappa ^{\prime }({{\bf x}})\,{\boldsymbol \nabla} \cdot {{\bf s}}({{\bf x}},t)\,{\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t) +2\,\mu ^{\prime }({{\bf x}})\,{{\bf D}}({{\bf x}},t):{{\bf D}}^\dagger ({{\bf x}},T-t)]\,{\mathrm{d}}t , \end{equation} (63)  \begin{eqnarray} {{\bf K}}_{{\bf \tau }^0}({{\bf x}}) &= & -\int _0^T\left\lbrace \frac{1}{2}\,{\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t)\,{{\bf D}}({{\bf x}},t) -\frac{1}{2}\,{\boldsymbol \nabla} \cdot {{\bf s}}({{\bf x}},t)\,{{\bf D}}^\dagger ({{\bf x}},T-t)\right. \nonumber\\ &&\left. -\,\frac{1}{2}\left[\kappa ^{\prime }({{\bf x}})+\frac{2}{3}\,\mu ^{\prime }({{\bf x}})\right][{\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t)\,{{\bf D}}({{\bf x}},t)+{\boldsymbol \nabla} \cdot {{\bf s}}({{\bf x}},t)\,{{\bf D}}^\dagger ({{\bf x}},T-t)]\right. \nonumber\\ && \left. -\,\mu ^{\prime }({{\bf x}})\,[{{\bf D}}({{\bf x}},t)\cdot {{\bf D}}^\dagger ({{\bf x}},T-t)+{{\bf D}}^\dagger ({{\bf x}},T-t)\cdot {{\bf D}}({{\bf x}},t)] -{\boldsymbol \omega} ({{\bf x}},t)\cdot {{\bf D}}^\dagger ({{\bf x}},T-t) +{{\bf D}}^\dagger ({{\bf x}},T-t)\cdot {\boldsymbol \omega} ({{\bf x}},t)\right. \nonumber\\ && \left. +\,{{\bf s}}({{\bf x}},t)\cdot {\boldsymbol \nabla} {\boldsymbol \epsilon} ^\dagger ({{\bf x}},T-t) -\frac{1}{2}\,{{\bf s}}({{\bf x}},t)\,{\boldsymbol \nabla} {\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t) -\frac{1}{2}\,[{\boldsymbol \nabla} {\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t)]\,{{\bf s}}({{\bf x}},t) \right\rbrace \,{\mathrm{d}}t . \end{eqnarray} (64)Here T denotes the record length. The adjoint wavefield, s†, satisfies the same wave equation as the regular wavefield, s, but is generated by the adjoint source   \begin{equation} f_i^\dagger ({{\bf x}},t)=\sum _{r=1}^{N_r}\left[d_i({{\bf x}}_r,T-t)-d_i^{\Delta T}({{\bf x}}_r,T-t)\right]\,\delta ({{\bf x}}-{{\bf x}}_r) . \end{equation} (65)This equation illustrates how time-lapse data, which are simultaneously insert at all receiver locations in reverse time, give rise to the adjoint wavefield. The strain deviator and its adjoint are denoted by D and D†, respectively. Note that the kernel $${{\bf K}}_{{\bf \tau }^0}$$ has zero trace, as expected, and that the second derivative of the adjoint wavefield, $${\boldsymbol \nabla} {\boldsymbol \nabla} {{\bf s}}^\dagger$$, is required for its calculation. Effects due to changes in density, δρ, elastic parameters, $$\delta \boldsymbol \Gamma$$, and topography on solid-solid boundaries, δh, are readily incorporated (Tromp et al.2005). In this case the derivative of the misfit function is modified as   \begin{equation} \delta \chi =\int _V (\delta \rho \,K_\rho +\delta {\boldsymbol \Gamma} :: {{\bf K}}_{{{\boldsymbol \Gamma} }}+p^0\,K_{p^0}+{\boldsymbol \tau} ^0:{{\bf K}}_{{\bf \tau }^0})\,{\mathrm{d}}V+\int _\Sigma \delta h\,K_h\,{\mathrm{d}}\,\Sigma , \end{equation} (66)where   \begin{equation} K_\rho =\int _0^T\partial _t{{\bf s}}^\dagger \cdot \partial _t{{\bf s}}\,{\mathrm{d}}t , \end{equation} (67)   \begin{equation} {{\bf K}}_{{{\boldsymbol \Gamma} }}=-\frac{1}{2}\int _0^T({\boldsymbol \epsilon} ^\dagger \,{\boldsymbol \epsilon} +{\boldsymbol \epsilon} \,{\boldsymbol \epsilon} ^\dagger )\,{\mathrm{d}}t , \end{equation} (68)  \begin{equation} K_h=-\int _0^T\left[\rho \,\partial _t{{\bf s}}^\dagger \cdot \partial _t{{\bf s}}-{\boldsymbol \epsilon} ^\dagger :{\boldsymbol \Gamma} :{\boldsymbol \epsilon} +\hat{{{\bf n}}}\,\partial _n{{\bf s}}^\dagger :{\boldsymbol \Gamma} :{\boldsymbol \epsilon} +\hat{{{\bf n}}}\,\partial _n{{\bf s}}:{\boldsymbol \Gamma} :{\boldsymbol \epsilon} ^\dagger \right]_-^+\,{\mathrm{d}}t . \end{equation} (69)Here $$\partial _n=\hat{{{\bf n}}}\cdot {\boldsymbol \nabla}$$, and the dependence of the adjoint wavefield on time T − t and the forward wavefield on time t has been omitted for brevity. Topographic perturbations on fluid-solid boundaries may also be accommodated (Dahlen 2005; Tromp et al.2005), but are omitted for brevity. For isotropic models, kernel KΓ defined in eq. (68) is replaced by kernels for the bulk and shear moduli (Tromp et al.2005), namely   \begin{equation} K_\kappa =-\int _0^T{\boldsymbol \nabla} \cdot {{\bf s}}({{\bf x}},t)\,{\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t) \,{\mathrm{d}}t , \end{equation} (70)  \begin{equation} K_\mu =-2\int _0^T{{\bf D}}({{\bf x}},t):{{\bf D}}^\dagger ({{\bf x}},T-t)]\,{\mathrm{d}}t . \end{equation} (71)Upon comparing the isotropic kernels (70) and (71) with kernel (63) for the induced pressure, we conclude that induced pressure trades off with the moduli, an unfortunate fact to be expected from expressions (37) and (38) for the wave speeds. As mentioned previously, for this reason it makes sense in both active and passive seismology to build the induced pressure dependence into the moduli by selecting κ + κ΄ p0 and μ + μ΄ p0 as model parameters, thereby eliminating induced pressure as an unknown parameter. However, one can imagine engineering applications in which the elastic moduli, κ & μ, and their pressure derivatives, κ΄ & μ΄, are known, and one seeks to determine induced pressure variations, for example, based on time-lapse ultrasonics. In that case kernel (63) may prove to be a useful tool. We note that the formulation of an iterative inverse problem based on the Fréchet derivatives presented in this section is non-trivial, because updates in induced stress need to continue to satisfy the equilibrium eq. (9), which, as previously mentioned, constrains three out of six elements of the induced stress. This might be accomplished by making use of the ‘Beltrami stress function’ (Gurtin 1963), or based on the orthogonal decompositions of stress tensors introduced by Al-Attar (2010). This important topic is beyond the scope of this paper. 8 2-D SH WAVES 8.1 Equation of motion and kernels Consider a 2-D SH wavefield polarized in the y-direction with displacement field s(x, z, t). The associated scalar wave equation is   \begin{equation} \rho \,\partial _t^2s=\partial _x(\mu \,\partial _x s)+\partial _z(\mu \,\partial _z s) , \end{equation} (72)where ρ = ρ(x, z) denotes mass density and μ = μ(x, z) the shear modulus. Upon inducing a purely deviatoric stress of the form $$\tau ^0_{xz}=\tau$$, the 2-D SH equation of motion becomes   \begin{eqnarray} \rho \,\partial _t^2s &=\partial _x\left(\mu \,\partial _x s+\frac{1}{4}\,\tau \,\partial _zs\right)+\partial _z\left(\mu \,\partial _z s+\frac{1}{4}\,\tau \,\partial _xs\right) . \end{eqnarray} (73)The misfit gradient may be written in the form   \begin{equation} \delta \chi =\int _V (\delta \ln \rho \,K_\rho +\delta \ln \mu \,K_\mu +\tau \,K_{\tau })\,{\mathrm{d}}^2{{\bf x}}+\int _S \delta h\,K_h\,{\mathrm{d}}^1{{\bf x}}, \end{equation} (74)where   \begin{equation} K_{\rho }= \int _0^T\rho \,\partial _ts(t)\,\partial _ts^\dagger (T-t)\,{\mathrm{d}}t , \end{equation} (75)  \begin{eqnarray} K_{\mu } =& -\int _0^T\mu \,[ \partial _xs(t)\,\partial _xs^\dagger (T-t)+\partial _zs(t)\,\partial _zs^\dagger (T-t)]\,{\mathrm{d}}t , \end{eqnarray} (76)  \begin{eqnarray} K_{\tau } =& -\frac{1}{4}\,\int _0^T [\partial _zs(t)\,\partial _xs^\dagger (T-t)+\partial _xs(t)\,\partial _zs^\dagger (T-t)]\,{\mathrm{d}}t , \end{eqnarray} (77)  \begin{eqnarray} K_h & = & -\int _0^T[\rho \,\partial _ts(t)\,\partial _ts^\dagger (T-t) -\mu \,\partial _xs(t)\,\partial _xs^\dagger (T-t) -\mu \,\partial _zs(t)\,\partial _zs^\dagger (T-t) \nonumber\\ && +\,\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _ns^\dagger (t)\,\partial _xs(T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial _ns^\dagger (t)\,\partial _xs(T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _ns(t)\,\partial _xs^\dagger (T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial _ns(t)\,\partial _xs^\dagger (T-t) ]_-^+ \,{\mathrm{d}}t . \end{eqnarray} (78)Here $$\partial _n=\hat{{{\bf n}}}\cdot {\boldsymbol \nabla} =\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _x+\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial_z$$. In terms of $$\beta =\sqrt{\mu /\rho }$$ and τ΄ = τ/μ, we have the preferred alternative expression   \begin{equation} \delta \chi =\int _V (\delta \ln \rho \,K^{\prime }_\rho +\delta \ln \beta \,K_\beta +\tau ^{\prime }\,K_{\tau ^{\prime }})\,{\mathrm{d}}^2{{\bf x}}+\int _S \delta \ln h\,K^{\prime }_h\,{\mathrm{d}}^1{{\bf x}}, \end{equation} (79)where   \begin{equation} K^{\prime }_\rho =K_\rho +K_\mu , \end{equation} (80)  \begin{equation} K_\beta =2\,K_\mu , \end{equation} (81)  \begin{eqnarray} K_{\tau ^{\prime }} =& -\frac{1}{4}\,\int _0^T\mu \, [\partial _zs(t)\,\partial _xs^\dagger (T-t)+\partial _xs(t)\,\partial _zs^\dagger (T-t)]\,{\mathrm{d}}t , \end{eqnarray} (82)  \begin{eqnarray} K^{\prime }_h & = & -\int _0^Th[\rho \,\partial _ts(t)\,\partial _ts^\dagger (T-t) -\mu \,\partial _xs(t)\,\partial _xs^\dagger (T-t) -\mu \,\partial _zs(t)\,\partial _zs^\dagger (T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _ns^\dagger (t)\,\partial _xs(T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial _ns^\dagger (t)\,\partial _xs(T-t) \nonumber\\ &&+\,\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _ns(t)\,\partial _xs^\dagger (T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial _ns(t)\,\partial _xs^\dagger (T-t) ]_-^+ \,{\mathrm{d}}t . \end{eqnarray} (83) 8.2 Examples Illustrating the effects of induced stress on elastic wave propagation for 2-D SH waves is both straightforward and representative of more realistic cases. We use the spectral-element method (Komatitsch & Vilotte 1998; Komatitsch & Tromp 1999; Peter et al.2011) to simulate wave propagation with and without induced stress, and for the calculation of the kernels. The first example is taken from Tromp et al. (2005) and is based on a 2-D homogeneous medium in a 80 km × 200 km box with values of ρ = 2600 kg m−3 and β = 3200 m s−1, corresponding to a shear modulus of μ = 26.6 GPa. The box has a free surface at the top and absorbing boundaries at the bottom and at the sides. The source is a point force in the direction perpendicular to the medium plane with a Gaussian source time function. The wavefield is perfectly circular, before touching any boundaries. If we add a 20 per cent induced stress to this medium, we see (Fig. 1) that the wavefield now becomes slightly elliptical, as predicted in Sections 5 and 5.1. Effectively, the wavefield sees an apparent anisotropic medium. Figure 1. View largeDownload slide Snapshot after 16 s of SH wave propagation in a homogeneous background medium. Without induced stress, the wavefield is perfectly circular (bottom panel). In the presence of induced stress, the wave sees an apparent anisotropic material (middle panel). The difference of the two wavefields (top panel) gives an indication of the magnitude of the stress effect. The wavefield difference is about 10 times smaller than the wavefield itself in the presence of an induced stress of 20 per cent of the shear modulus. The units of the gird are in m and of the colour scale in 10 − 2 μm. Figure 1. View largeDownload slide Snapshot after 16 s of SH wave propagation in a homogeneous background medium. Without induced stress, the wavefield is perfectly circular (bottom panel). In the presence of induced stress, the wave sees an apparent anisotropic material (middle panel). The difference of the two wavefields (top panel) gives an indication of the magnitude of the stress effect. The wavefield difference is about 10 times smaller than the wavefield itself in the presence of an induced stress of 20 per cent of the shear modulus. The units of the gird are in m and of the colour scale in 10 − 2 μm. Keeping the same background model, we now calculate the kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ as specified by eqs (80)–(82). In terms of shape, we see (Fig. 2) that the kernels have the structure we would expect. In particular the diagonal symmetry seen in $$K_{\tau ^{\prime }}$$ is the result of the cross-product of partial derivatives (compare equations 76 and 82). In terms of amplitude, although we plot the kernels on a saturated scale to visualize the details, the β kernel is roughly an order of magnitude larger than both $$K^{\prime }_\rho$$ and $$K_{\tau ^{\prime }}$$. Figure 2. View largeDownload slide Banana-doughnut kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in a homogeneous background medium. The units of the grid are in m and of the colour scale in 10−9 m−2. Figure 2. View largeDownload slide Banana-doughnut kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in a homogeneous background medium. The units of the grid are in m and of the colour scale in 10−9 m−2. To investigate a more representative case for oil & gas exploration, we implement a simplified model for the Groningen gas field in the Netherlands (Kraaijpool & Dost 2013) in a 4 km × 4 km box. The top boundary acts as a free surface and the other 3 boundaries are absorbing. The model consist of seven layers (Table 1), where the gas is contained in the Rotliegend group overlain by a fast salt layer of the Zechstein group. We implement an array of 10 seismometers at the surface and another array of 10 instruments in a borehole. We assume that gas extraction induces an earthquake in the gas layer. The salt layer acts like a strong barrier to wave propagation and only little energy escapes to the surface (Fig. 3). Figure 3. View largeDownload slide Snapshot after 1.5 s of SH wave propagation in the isotropic Groningen model (Table 1). The seismometers are marked by green squares and the event by an orange cross. The model is represented by the grey background, where the darkest colour represents the fastest layer. Most energy of the waveform remains in and below the salt layer, only the first wave being transmitted. Figure 3. View largeDownload slide Snapshot after 1.5 s of SH wave propagation in the isotropic Groningen model (Table 1). The seismometers are marked by green squares and the event by an orange cross. The model is represented by the grey background, where the darkest colour represents the fastest layer. Most energy of the waveform remains in and below the salt layer, only the first wave being transmitted. Table 1. Model parameters for the simplified Groningen gas field. Unit  Depth (m)  β (m s−1)  ρ (kg m−3)  North Sea supergroup  0  372  1960  Chalk group  828  756  2300  Rijnland group  1713  663  2650  Upper and lower germanic Trias groups  1820  1427  2600  Zechstein group  1925  2725  2100  Upper Rotliegend group  2866  1184  2550  Limburg group  3126  1695  2700  Unit  Depth (m)  β (m s−1)  ρ (kg m−3)  North Sea supergroup  0  372  1960  Chalk group  828  756  2300  Rijnland group  1713  663  2650  Upper and lower germanic Trias groups  1820  1427  2600  Zechstein group  1925  2725  2100  Upper Rotliegend group  2866  1184  2550  Limburg group  3126  1695  2700  View Large An important question is whether induced stress is measurable by seismic time-lapse experiments. Making simplifying assumptions, Hatchell & Bourne (2005) related measured time-lapse time shifts to strain changes induced by extraction, showing values of 10–20 milliseconds for several off-shore reservoirs. In our simulations, we keep the elastic constants fixed and only add induced stresses due to gas extraction, which is then equivalent to induced strains. Dempsey & Suckale (2017) estimated induced stresses in the Groningen gas field and found values of the order of 75 MPa in the gas layer, corresponding to 2 per cent of the shear modulus in that layer. We systematically vary the induced stress up to 20 per cent of the shear modulus in the Rotliegend and measure cross-correlation time shifts compared to the case with no induced stress for a station at the surface vertically above the event (S1), the most extreme right station at the surface (S10), and the borehole station just below the gas layer (B18), which receives the bulk of the seismic energy. No time shifts are detected at S1, which receives a direct wave that travelled the least in the stressed medium. S10 and B18 both see maximum time shifts of 1.5 and 3.0 ms, respectively, for stresses at 20 per cent of the shear modulus. These time shifts are an order of magnitude smaller than those reported by Hatchell & Bourne (2005). The time shifts of a few tens of milliseconds seen by Hatchell & Bourne (2005) are thus more likely due to boundary or velocity changes induced by the stress changes, rather than the stress changes themselves, as the latter would require substantial stresses to achieve such large time shifts. Looking at the waveforms however, we notice that the seismic coda shows more significant differences while the first arriving waves present tiny time shifts (Fig. 4). To quantify this, we evaluate the average relative waveform difference δΦ defined as   \begin{equation} \delta \Phi = \frac{1}{2}\,\log \left(\frac{\int d_r^2\,{\mathrm{d}}t}{\int d_0^2\,{\mathrm{d}}t}\right) , \end{equation} (84)where d0 is the seismogram in the unstressed medium and dτ the seismogram in the medium with induced stress. Measuring δΦ systematically as a function of induced stress, we see that significant changes in waveforms are noticeable even for modest stresses (Fig. 5). While in this particular setting, some induced stress kernels are small, others, especially in deep boreholes, clearly show that induced stress has a large sensitivity on the waveforms (Figs 6–8). The conclusion of these observations is that it is possible to detect induced stresses of magnitudes we expect in reservoir settings in time lapse experiments by looking at full waveform differences rather than time shifts of first arriving waves. Our 2-D examples are designed to show the order of magnitude of the waveform effects to be expected, and hence our simple waveform difference measure δΦ is sufficient. In real applications, a more rigorous quantification of the waveform difference should be used (e.g. Kristekova et al.2009). We anticipate that in three dimensions, shear wave splitting caused by induced deviatoric stress may also be a detectable. Figure 4. View largeDownload slide Sample seismogram at station B18 (location shown on Fig. 3) in the Groningen medium with a 20 per cent induced stress relative to the shear modulus (black) compared to the seismogram in the unstressed Groningen medium (red). The bottom panel shows the complete seismograms. The top right panel is a zoom into the first arriving wave. The top left panel is a zoom into the coda. All three panels used the same units specified in the bottom panel. Figure 4. View largeDownload slide Sample seismogram at station B18 (location shown on Fig. 3) in the Groningen medium with a 20 per cent induced stress relative to the shear modulus (black) compared to the seismogram in the unstressed Groningen medium (red). The bottom panel shows the complete seismograms. The top right panel is a zoom into the first arriving wave. The top left panel is a zoom into the coda. All three panels used the same units specified in the bottom panel. Figure 5. View largeDownload slide Average relative waveform differences δΦ measured in stressed models with respect to the unstressed Groningen medium. The different curves correspond to waveform changes at different stations located in Fig. 3. Figure 5. View largeDownload slide Average relative waveform differences δΦ measured in stressed models with respect to the unstressed Groningen medium. The different curves correspond to waveform changes at different stations located in Fig. 3. Figure 6. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station S1 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 6. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station S1 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 7. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station S10 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 7. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station S10 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 8. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station B18 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 8. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station B18 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. 9 CONCLUDING REMARKS It is well known that seismic wave speeds in crustal rocks increase with hydrostatic pressure (Birch 1961; Nur & Simmons 1969). Under non-hydrostatic stress, the same wave speeds develop clear characteristics of anisotropy (Nur & Simmons 1969). Mavko et al. (1995), for instance, explained these observations by introducing cracks with prescribed compliances and further assuming a porosity for each crack. Another approach is to use higher-order elasticity (e.g. Johnson & Rasolofosaon 1996; Prioul et al.2004) to explain effective anisotropy induced by deviatoric stress. Johnson & Rasolofosaon (1996) plot squared compressional wave speed $$c_p^2$$ as a function of applied uniaxial stress. Using our expression (37), we find that $$\rho \,c_p^2$$ changes linearly with uniaxial stress P. Recalling Birch’s law (Birch 1961), stating that density is proportional to P-wave speed, we find that $$c_p^2 \sim P^{2/3}$$, similar to the data shown in Johnson & Rasolofosaon (1996). These authors further plot S-wave birefringence data, which present a trend in sin 2θ, where θ is the angle between the propagation direction and the stress direction. Our eq. (38) similarly exhibits such a sin 2θ dependence without the need for third-order elasticity. Prioul et al. (2004) show data where elastic constants vary linearly with applied hydrostatic stress. In our formulation we clearly identify the linear hydrostatic pressure dependence of the elastic constants (eqs 37 and 38). Motivated by the formulation common in global seismology, we show the effects of induced stress on the elastic wave equation and constitutive relation. Without employing higher-order elasticity, our formulation leads to trends observed in measurements based on laboratory data. We also, for the first time, give expressions for induced stress sensitivity kernels, which can be used in adjoint tomography. From an inverse theory perspective, induced pressure trades off with the isotropic moduli, and for this reason it should be built into the definition of ‘the moduli’, thereby eliminating it as an unknown parameter. Induced deviatoric stress results in anisotropic seismic wave speeds, and thus it is possible to image deviatoric stress using time-lapse experiments. Simple 2-D SH experiments suggest waveforms are sufficiently sensitive to stress changes expected in reservoir settings, but we will have to deal with trade-offs between changes in elastic constants, boundaries, and induced deviatoric stress, which all affect the wavefield simultaneously. The theory contains two quantities a and b (see eq. 27), which may be chosen to obtain a particular equation of state. Global seismologists (e.g. Dahlen 1972a; Dahlen & Tromp 1998) prefer to choose a = −b = 1/2, thereby rendering the formulation independent of the hydrostatic pre-stress. An alternative second-order equation of state is based on logarithmic or Hencky strain (Poirier & Tarantola 1998) rather than the Lagrangian strain, and may be recovered from our theory by setting a = 0 and b = −1/2. Our analysis suggests a physical interpretation of the quantities a and b that gives the isotropic bulk and shear moduli an induced pressure dependence of the form κ + κ΄ p0 and μ + μ΄ p0, respectively, where κ and μ denote isotropic elastic moduli at zero pressure, and κ΄ and μ΄ represent their adiabatic pressure derivatives. Note that—because the derivatives κ΄ and μ΄ may be depth-dependent—the quantities κ + κ΄ p0 and μ + μ΄ p0 do not necessarily vary linearly with pressure. For moderate stress changes, Birch-Murnaghan’s second-order equation of state in Eulerian strain is often used (e.g. Stacey 1992). This formulations leads to κ΄ = 4, and together with the empirical relation μ = 0.631 κ − 0.899 P (see Stacey 1992), results in a = −0.95833 and b = −1.3125, values that are quite different from those used in global seismology. It is important to recognize that whereas Dahlen (1972a,b) and Dahlen & Tromp (1998) choose to set a = −b = 1/2, we do not select specific values for a and b. By rewriting a and b in terms of two new parameters, μ΄ and κ΄, we show that the predicted seismic wave speeds, defined by eqs (37) and (38), take on experimentally expected forms when one interprets these new parameters as the adiabatic pressure derivatives of the shear and bulk moduli with respect to pressure. But note also that based on this interpretation one no longer simply ‘chooses’ a and b. Instead, the values of a and b are determined by the pressure derivatives of the shear and bulk moduli. Our theory makes testable predictions given values of the shear and bulk moduli and their pressure derivatives and an applied pre-stress, for example, a directional dependence of the compressional wave speed (37) and the shear wave split time (41). We predict that the compressional wave speed acquires a directional dependence with an amplitude linearly proportional to κ΄ and μ΄, whereas setting a = −b = 1/2—equivalent to setting μ΄ = κ΄ = 0 in our formulation—results in an isotropic P wave. Similarly, we predict that for two materials with the same shear modulus μ but different pressure derivatives μ΄, the split time will be different, whereas setting a = −b = 1/2 would not show any differences in split time. One could imagine lab experiments that could attempt to verify these predictions given experimentally determined values of the moduli μ and κ and their pressure derivatives μ΄ and κ΄ for various materials. Acknowledgements We thank an anonymous reviewer, David Al-Attar, and Editor Martin Mai for helpful comments which helped to improve the manuscript. The open source spectral-element package SPECFEM2D used in this study is freely available via the Computational Infrastructure for Geodynamics (CIG). The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP/2007-2013) grant agreement number 320639 (iGEO). REFERENCES Al-Attar D., 2010. On the parametrisation of equilibrium stress fields in the Earth, Geophys. J. Int. , 181, 567– 576. https://doi.org/10.1111/j.1365-246X.2010.04527.x Birch F., 1961. The velocity of compressional waves in rocks to 10 kilobars: 2, J. geophys. Res. , 66, 2199– 2224. Google Scholar CrossRef Search ADS   Bruner W.M., 1976. Comment on ‘Seismic velocities in dry and saturated cracked solids’ by Richard J. O’Connell and Bernard Budiansky, J. geophys. Res. , 81, 2573– 2576. Google Scholar CrossRef Search ADS   Dahlen F.A., 1972a. Elastic dislocation theory for a self-gravitating elastic configuration with an initial static stress field, Geophys. J. R. astr. Soc. , 28, 357– 383. https://doi.org/10.1111/j.1365-246X.1972.tb06798.x Google Scholar CrossRef Search ADS   Dahlen F.A., 1972b. Elastic velocity anisotropy in the presence of an anisotropic initial stress, Bull. seism. Soc. Am. , 62, 1183– 1193. Dahlen F.A., 2005. Finite-frequency sensitivity kernels for boundary topography perturbations, Geophys. J. Int. , 162, 525– 540. https://doi.org/10.1111/j.1365-246X.2005.02682.x Dahlen F.A., Tromp J., 1998. Theoretical Global Seismology , Princeton University Press. Dempsey D., Suckale J., 2017. Physics-based forecasting of induced seismicity at the Groningen gas filed, The Netherlands, Geophys. Res. Lett. , 44, doi:10.10022017GL073878. https://doi.org/10.1002/2017GL073878 Eberhart-Phillips D., Han D.-H., Zoback M., 1989. Empirical relationships among seismic velocity, effective pressure, porosity, and clay content in sandstone, Geophysics , 54, 82– 89. https://doi.org/10.1190/1.1442580 Google Scholar CrossRef Search ADS   Egle D., Bray D.E., 1976. Measurement of acoustoelastic and third-order elastic constants for rail steel, J. acoust. Soc. Am. , 60, 741– 744. https://doi.org/10.1121/1.381146 Google Scholar CrossRef Search ADS   Fazio T., Aki L., Alba K., 1973. Solid earth tide and observed change in the in situ seismic velocity, J. geophys. Res. , 78, 1319– 1322. Google Scholar CrossRef Search ADS   Gurtin M., 1963. A generalization of the Beltrami stress functions in continuum mechanics, Arch. Ration. Mech. Anal. , 13, 321– 329. https://doi.org/10.1007/BF01262700 Google Scholar CrossRef Search ADS   Hatchell P., Bourne S., 2005. Rocks under strain: strain-induced time-lapse time shifts are observed for depleting reservoirs, Leading Edge , December, 1222– 1225. Henyey F.S., Pomphrey N., 1982. Self-consistent elastic moduli of a cracked solid, Geophys. Res. Lett. , 81, 903– 906. https://doi.org/10.1029/GL009i008p00903 Google Scholar CrossRef Search ADS   Hughes D., Kelly J.L., 1953. Second-order elastic deformation of solids, Phys. Rev. , 92, 1145– 1149. https://doi.org/10.1103/PhysRev.92.1145 Google Scholar CrossRef Search ADS   Johnson P., Rasolofosaon P., 1996. Nonlinear elasticity and stress-induced anisotropy in rock, J. geophys. Res. , 101, 3113– 3124. Google Scholar CrossRef Search ADS   Komatitsch D., Tromp J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation, Geophys. J. Int. , 139, 806– 822. https://doi.org/10.1046/j.1365-246x.1999.00967.x Komatitsch D., Vilotte J.-P., 1998. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. seism. Soc. Am. , 88, 368– 392. Kraaijpool D., Dost B., 2013. Implications of salt-related propagation and mode conversion effects on the analysis of induced seismicity, J. Seismol. , 17, 95– 107. https://doi.org/10.1007/s10950-012-9309-4 Google Scholar CrossRef Search ADS   Kristekova M., Kristek J., Moczo P., 2009. Time frequency misfit and goodness-of-fit criteria, Geophys. J. Int. , 178, 813– 825. https://doi.org/10.1111/j.1365-246X.2009.04177.x Mavko G., Mukerji T., Godfrey N., 1995. Predicting stress-induced velocity anisotropy in rocks, Geophysics , 60, 1081– 1087. https://doi.org/10.1190/1.1443836 Google Scholar CrossRef Search ADS   Murnaghan F., 1951. Finite Deformation of an Elastic Solid , Dover Publications Inc. Nur A., 1971. Effects of stress on velocity anisotropy in rocks with cracks, J. geophys. Res. , 76, 2022– 2024. Google Scholar CrossRef Search ADS   Nur A., Simmons G., 1969. Stress-induced velocity anisotropy in rock: an experimental study, J. geophys. Res. , 74, 6667– 6674. Google Scholar CrossRef Search ADS   O’Connell R., Budiansky B., 1974. Seismic velocities in dry and saturated cracked solids, J. geophys. Res. , 79, 5412– 5426. Google Scholar CrossRef Search ADS   Peter D. et al.  , 2011. Forward and adjoint simulations of seismic wave propagation on unstructured hexahedral meshes, Geophys. J. Int. , 186, 721– 739. https://doi.org/10.1111/j.1365-246X.2011.05044.x Poirier J., Tarantola A., 1998. A logarithmic equation of state, Phys. Earth planet. Inter. , 109, 1– 8. https://doi.org/10.1016/S0031-9201(98)00112-5 Google Scholar CrossRef Search ADS   Prioul R., Bakulin A., Bakulin V., 2004. Nonlinear rock physics model for estimation of 3D subsurface stress in anisotropic formations: theory and laboratory verification, Geophysics , 69, 415– 425. https://doi.org/10.1190/1.1707061 Google Scholar CrossRef Search ADS   Renaud G., Bas P.L., Johnson P.A., 2012. Revealing highly complex elastic nonlinerar (anelastic) behaviour of earth materials applying a newprobe: dynamic accoustic testing, J. geophys. Res. , 117, doi:10.10292011JB009127. Silver P., Daley T., Niu F., Majer E., 2007. Active source monitoring of cross-well seismic travel time for stress-induced changes, Bull. seism. Soc. Am. , 97, 281– 293. https://doi.org/10.1785/0120060120 Stacey F., 1992. Physics of the Earth , Brookfield Press. Tromp J., Tape C., Liu Q.Y., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, Geophys. J. Int. , 160, 195– 216. https://doi.org/10.1111/j.1365-246X.2004.02453.x Verdon J., Angus D., Kendall M., Hall S., 2008. The effect of microstructure and nonlinear stress on anisotropic seismic velocities, Geophysics , 78, D41– D51. https://doi.org/10.1190/1.2931680 Google Scholar CrossRef Search ADS   Zheng Z., 2000. Seismic anisotropy due to stress induced cracks, Int. J. Rock Mech. Min. Sci. , 37, 39– 49. https://doi.org/10.1016/S1365-1609(99)00090-8 Google Scholar CrossRef Search ADS   APPENDIX: WEAK FORM OF THE EQUATION OF MOTION UNDER INDUCED STRESS We have the equation of motion eq. (16), that is,   \begin{equation} \rho \,\partial _t^2{{\bf s}}={\boldsymbol \nabla} \cdot ({{\bf T}}^{\rm L1} -{{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0) , \end{equation} (A1)subject to the boundary condition (18), which is   \begin{equation} \hat{{{\bf n}}}\cdot [{{\bf T}}^{\rm L1}+{{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0]_-^+={{\bf 0}}. \end{equation} (A2)Upon dotting the equation of motion (A1) with a test vector w, followed by an integration by parts in which we use the boundary condition (A2), the weak form of the equation of motion is   \begin{eqnarray} \int _V\rho \,{{\bf w}}\cdot \partial _t^2{{\bf s}}\,{\mathrm{d}}V \ & = & \ -\int _V({\boldsymbol \nabla} {{\bf w}}):({{\bf T}}^{\rm L1}-{{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0):({\boldsymbol \nabla} {{\bf w}})\,{\mathrm{d}}\,V - \int _\Sigma \hat{{{\bf n}}}\cdot [({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0 -{{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-{{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0]_-^+\cdot {{\bf w}}\,{\mathrm{d}}\,\Sigma . \end{eqnarray} (A3)Another lengthy, tedious calculation enables us to eliminate the surface integral and derivatives of the induced deviatoric stress, resulting in the weak form   \begin{eqnarray} \int _V\rho \,{{\bf w}}\cdot \partial _t^2{{\bf s}}\,{\mathrm{d}}V & = & -\int _V({\boldsymbol \nabla} {{\bf w}}):({{\bf T}}+\Delta {{\bf T}})\,{\mathrm{d}}\,V + \int _V[({\boldsymbol \tau} ^0:{\boldsymbol \nabla} {{\bf s}})\,{\boldsymbol \nabla} \cdot {{\bf w}}+{{\bf s}}\cdot {\boldsymbol \tau} ^0\cdot ({\boldsymbol \nabla} {\boldsymbol \nabla} \cdot {{\bf w}}) -({\boldsymbol \nabla} \cdot {{\bf s}})\,{\boldsymbol \tau} ^0:{\boldsymbol \nabla} {{\bf w}}-{\boldsymbol \tau} ^0:({{\bf s}}\cdot {\boldsymbol \nabla} {\boldsymbol \nabla} {{\bf w}}) ]\,{\mathrm{d}}\,V . \end{eqnarray} (A4)Unfortunately, the second derivative of the test function, $${\boldsymbol \nabla} {\boldsymbol \nabla} {{\bf w}}$$, is required. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Effects of induced stress on seismic forward modelling and inversion

Loading next page...
 
/lp/ou_press/effects-of-induced-stress-on-seismic-forward-modelling-and-inversion-UDQW4fWTIU
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy020
Publisher site
See Article on Publisher Site

Abstract

Summary We demonstrate how effects of induced stress may be incorporated in seismic modelling and inversion. Our approach is motivated by the accommodation of pre-stress in global seismology. Induced stress modifies both the equation of motion and the constitutive relationship. The theory predicts that induced pressure linearly affects the unstressed isotropic moduli with a slope determined by their adiabatic pressure derivatives. The induced deviatoric stress produces anisotropic compressional and shear wave speeds; the latter result in shear wave splitting. For forward modelling purposes, we determine the weak form of the equation of motion under induced stress. In the context of the inverse problem, we determine induced stress sensitivity kernels, which may be used for adjoint tomography. The theory is illustrated by considering 2-D propagation of SH waves and related Fréchet derivatives based on a spectral-element method. Elasticity and anelasticity, Equations of state, High-pressure behaviour, Computational seismology, Theoretical seismology, Wave propagation 1 INTRODUCTION In many engineering applications, it is important to be able to monitor stress changes during operations. Obvious examples are drilling, oil and gas exploration, geothermal activity, and mining. In geomechanical modelling, which describes the response of the underground to mechanical perturbations, knowledge of the state of stress is essential. Most often estimates of stress are obtained from seismic data using simplifying assumptions or empirical models. These empirical models and assumptions are derived from laboratory rock physics experiments. We propose a method which quantitatively links changes in seismic wavefields directly to changes in stress using first principles. Basic effects of changes in stress on seismic wave speeds have been known for a long time (e.g. Birch 1961; Nur & Simmons 1969), and such effects have been observed in laboratory (e.g. Eberhart-Phillips et al.1989; Verdon et al.2008) and field studies (e.g. Fazio et al.1973; Silver et al.2007). Exploration for oil & gas induces stress changes inside and outside reservoirs (e.g. Hatchell & Bourne 2005). As a consequence, seismic waves are observed to have different traveltimes before and after depletion of the reservoir or after waste water injection. Such time-lapse time shifts are caused by two competing factors, namely, changes in layer geometry and stress-induced changes in seismic wave speeds. Two descriptions of the effects of an induced stress on seismic wave speeds have been commonly used. The first approach focuses on stress effects on pre-existing or induced cracks, which manifest themselves in the form of seismic anisotropy (e.g. Nur 1971; O’Connell & Budiansky 1974; Zheng 2000). O’Connell & Budiansky (1974) incorrectly accounted for the crack energy as the crack density increases, an issue addressed by Bruner (1976) and Henyey & Pomphrey (1982). The second approach is based on third-order elasticity theory (Murnaghan 1951; Hughes & Kelly 1953; Egle & Bray 1976), which requires the knowledge of higher-order elastic constants. The latter are not easily measured in the laboratory, with many discrepancies remaining, although new methods are continuously proposed (e.g. Renaud et al.2012). In this paper, we take an approach motivated by the accommodation of pre-stress in global seismology. As first discussed in Dahlen (1972a,b) and also in Dahlen & Tromp (1998, sections 3.3.2 and 3.6.2), pre-stress affects both the equation of motion and the constitutive relationship. Here we explore these effects both from a forward modelling perspective and from the point of view of the inverse problem. 2 INITIAL EQUATION OF MOTION Prior to inducing stress, the equation of motion is   \begin{equation} \rho \,\partial _t^2{{\bf s}}-\boldsymbol \nabla \cdot {{\bf T}}={{\bf 0}}, \end{equation} (1)where ρ denotes mass density, s the displacement vector, and T the symmetric stress tensor. At internal discontinuities and the free surface, the traction needs to be continuous, that is,   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}]_-^+={{\bf 0}}, \end{equation} (2)where $$\hat{{{\bf n}}}$$ denotes the unit normal pointing from the − side of a discontinuity to the + side, and where the notation $$[\,\cdot \,]_-^+$$ denotes the jump in the enclosed quantity when going from the − to the + side. On a traction-free surface, such as Earth’s surface, traction vanishes: $$\hat{{{\bf n}}}\cdot {{\bf T}}={{\bf 0}}$$. The stress tensor, T, is related to the symmetric infinitesimal strain tensor,   \begin{eqnarray} \boldsymbol \epsilon =&\frac{1}{2}[\boldsymbol \nabla {{\bf s}}+(\boldsymbol \nabla {{\bf s}})^T)] , \end{eqnarray} (3)(a superscript T denotes the transpose) via Hooke’s law:   \begin{equation} {{\bf T}}= \boldsymbol {\Gamma} :\boldsymbol {\epsilon} . \end{equation} (4)The fourth-order elastic tensor, $$\boldsymbol \Gamma$$, linearly relates the stress and strain tensors. Let U denote the elastic energy per unit mass. Then the stress may be related to the elastic energy density, defined to second order in strain, as   \begin{eqnarray} \rho \,U=&\frac{1}{2}\,\boldsymbol \epsilon :\boldsymbol \Gamma :\boldsymbol \epsilon\!, \end{eqnarray} (5)via   \begin{equation} {{\bf T}}=\rho \,\frac{\partial U}{\partial {\boldsymbol \epsilon} } . \end{equation} (6)The symmetries of the stress and strain tensors and the quadratic nature of the elastic energy density dictate that the fourth-order elastic tensor, with elements   \begin{equation} \Gamma _{ijk\ell}=\frac{\partial ^2 U}{\partial \epsilon _{ij}\,\partial \epsilon _{k\ell }}, \end{equation} (7)exhibits the symmetries Γijkℓ = Γjikℓ, Γijkℓ = Γijℓk, and Γijkℓ = Γkℓij, which, in the most general case, reduces the number of independent parameters from 81 to 21. It is often convenient to express the elastic tensor in terms of its isotropic and purely anisotropic parts as   \begin{eqnarray} \Gamma _{ijk\ell }=&\left(\kappa -\frac{2}{3}\,\mu \right)\,\delta _{ij}\,\delta _{k\ell }+\mu \,(\delta _{ik}\,\delta _{j\ell }+\delta _{i\ell }\,\delta _{jk})+\gamma _{ijk\ell } , \end{eqnarray} (8)where κ and μ denote the isotropic bulk and shear moduli, respectively, and γijkℓ the purely anisotropic contribution. The elements γijkℓ exhibit the same symmetries as the elements Γijkℓ, and for purely isotropic media γijkℓ = 0. 3 INDUCED STRESS We induce a symmetric static stress field, T0, subject to the equilibrium condition   \begin{equation} {\boldsymbol \nabla} \cdot { {\bf T}}^0 ={{\bf 0}}. \end{equation} (9)The associated boundary condition is   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}^0]_-^+={{\bf 0}}. \end{equation} (10)The induced stress may be written in terms of an induced pressure,   \begin{eqnarray} p^0=&-\frac{1}{3}\,{\rm tr}({{\bf T}}^0) , \end{eqnarray} (11)and a symmetric trace-free induced deviatoric stress   \begin{eqnarray} {\boldsymbol \tau}^{0}=&{{\bf T}}^{0}-\frac{1}{3}\,{\rm tr}({{\bf T}}^{0})\,{{\bf I}}, \end{eqnarray} (12)such that,   \begin{equation} {{\bf T}}^0=-p^0\,{{\bf I}}+{\boldsymbol \tau} ^0 , \end{equation} (13)where I denotes the identity tensor. Taken together, the equilibrium condition and decomposition imply that   \begin{equation} {\boldsymbol \nabla} p^0={\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0 , \end{equation} (14)which enables us to eliminate $${\boldsymbol \nabla} p^0$$ in terms of $${\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0$$. For example, we have   \begin{equation} {\boldsymbol \nabla} {\bf T}^0=-{\boldsymbol \nabla} p^0\,{{\bf I}}+{\boldsymbol \nabla} {\boldsymbol \tau} ^0=-{\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0\,{{\bf I}}+{\boldsymbol \nabla} {\boldsymbol \tau} ^0 . \end{equation} (15)Note that the equilibrium condition (9) constrains only three of the six elements of the induced stress tensor. 4 EQUATION OF MOTION UNDER INDUCED STRESS Seismic waves propagating in a stress-induced medium are governed by the wave equation (Dahlen & Tromp 1998, eq. 3.56 without density and gravity perturbations or rotational terms)   \begin{equation} \rho \,\partial _t^2{{\bf s}}-{\boldsymbol \nabla} \cdot {{\bf T}}^{\rm L1}+{\boldsymbol \nabla} \cdot ({{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0)={{\bf 0}}, \end{equation} (16)where TL1 denotes the symmetric incremental Lagrangian Cauchy stress, and where the term $$-{{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0$$ captures advection of the induced stress. Using decomposition (13) and equilibrium eq. (9), the equation of motion (16) may be rewritten in the form (Dahlen & Tromp 1998, eq. 3.58 without density and gravity perturbations or rotational terms)   \begin{equation} \rho \,\partial _t^2{{\bf s}}-{\boldsymbol \nabla} \cdot {{\bf T}}^{\rm L1}-{\boldsymbol \nabla} [{{\bf s}}\cdot ({\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0)]+{\boldsymbol \nabla} \cdot ({{\bf s}}\cdot {\boldsymbol \nabla} {\boldsymbol \tau} ^0) ={{\bf 0}}. \end{equation} (17)The associated boundary condition is most conveniently expressed in terms of the asymmetric incremental first Piola–Kirchhoff stress, TPK1, namely (Dahlen & Tromp 1998, eq. 3.65),   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}^{\rm PK1}]_-^+={{\bf 0}}. \end{equation} (18)The incremental first Piola–Kirchhoff stress is related to the incremental Lagrangian Cauchy stress via (Dahlen & Tromp 1998, eq. 3.36)   \begin{equation} {{\bf T}}^{\rm PK1}={{\bf T}}^{\rm L1}+{{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0 . \end{equation} (19)For numerical simulations, in particular spectral-element simulations, the weak form of the equation of motion is discussed in Appendix. What we need next is a constitutive relationship for seismic wave propagation in a medium under induced stress. 5 ELASTIC STRAIN ENERGY DENSITY UNDER INDUCED STRESS A discussion of elastic strain energy requires careful definition of ‘the strain tensor’, and consideration of second-order quantities in the displacement gradient. As discussed in Dahlen & Tromp (1998, section 3.6.1), in an adiabatic, perfectly elastic material, the Lagrangian internal energy density depends only on the Lagrangian strain tensor   \begin{eqnarray} {{\bf E}}^{\rm L} =& \frac{1}{2}[{\boldsymbol \nabla} {{\bf s}}+({\boldsymbol \nabla} {{\bf s}})^T)]+\frac{1}{2}({\boldsymbol \nabla} {{\bf s}})\cdot ({\boldsymbol \nabla} {{\bf s}})^T , \end{eqnarray} (20)which may, alternatively, be expressed in terms of the infinitesimal strain tensor $${\boldsymbol \epsilon}$$ defined in eq. (3) and the antisymmetric infinitesimal vorticity tensor   \begin{eqnarray} {\boldsymbol \omega} &=-\frac{1}{2}\,[{\boldsymbol \nabla} {{\bf s}}-({\boldsymbol \nabla} {{\bf s}})^T] , \end{eqnarray} (21)as   \begin{eqnarray} {{\bf E}}^{\rm L} &= {\boldsymbol \epsilon} +\frac{1}{2}\,{\boldsymbol \epsilon} \cdot {\boldsymbol \omega} -\frac{1}{2}\,{\boldsymbol \omega} \cdot {\boldsymbol \epsilon} -\frac{1}{2}\,{\boldsymbol \omega} \cdot {\boldsymbol \omega} +\frac{1}{2}\,{\boldsymbol \epsilon} \cdot {\boldsymbol \epsilon} . \end{eqnarray} (22)The internal energy density under induced stress must be calculated correct to second order in the displacement gradient, $${\boldsymbol \nabla} {{\bf s}}$$, which prohibits the linearization $${{\bf E}}^{\rm L}\approx {\boldsymbol \epsilon}$$ in energetic considerations. Keeping this in mind, to second order in $${\boldsymbol \nabla} {{\bf s}}$$ the Lagrangian internal energy density may be expressed in the form   \begin{eqnarray} \rho ^0\,U^{\rm L}&={{\bf T}}^0:{{\bf E}}^{\rm L}+\frac{1}{2}\,{{\bf E}}^{\rm L}:{\boldsymbol \Xi} :{{\bf E}}^{\rm L} , \end{eqnarray} (23)where ρ0 denotes the density before straining the material, making UL the Lagrangian internal energy per unit mass. We have assumed for convenience that the Lagrangian internal energy density vanishes in the absence of strain. The symmetric second Piola-Kirchhoff stress is defined in terms of the Lagrangian internal energy via   \begin{eqnarray} {{\bf T}}^{\rm SK} \ &\equiv & \ \rho ^0\,\frac{\partial U^{\rm L}}{\partial {{\bf E}}^{\rm L}} \nonumber\\ & = & \ {{\bf T}}^0+{{\bf T}}^{\rm SK1} , \end{eqnarray} (24)where TSK1 denotes the symmetric incremental second Piola–Kirchhoff stress, namely,   \begin{eqnarray} {{\bf T}}^{\rm SK1} \ & = & \ {\boldsymbol \Xi} :{{\bf E}}^{\rm L} \nonumber\\ &\approx & \ {\boldsymbol \Xi} :{\boldsymbol \epsilon} , \end{eqnarray} (25)where the last equality holds to first order in the displacement gradient. The components of the fourth-order tensor $${\boldsymbol \Xi}$$ are given by   \begin{equation} \Xi _{ijk\ell }=\rho ^0\,\frac{\partial ^2 U^{\rm L}}{\partial E_{ij}^{\rm L}\,\partial E_{k\ell }^{\rm L}} , \end{equation} (26)implying that $${\boldsymbol \Xi}$$ exhibits the symmetries Ξijkℓ = Ξjikℓ, Ξijkℓ = Ξijℓk, and Ξijkℓ = Ξkℓij, thereby reducing its number of independent elements from 81 to 21. As discussed in Dahlen & Tromp (1998, section 3.6.2) and Dahlen (1972a), the tensor $${\boldsymbol \Xi}$$ is one of a family of permissible fourth-order tensors which may be used to define the Lagrangian internal energy density, and, therefore, the incremental second Piola–Kirchhoff stress. Specifically, we may write   \begin{equation} \Xi _{ijk\ell }=\Gamma _{ijk\ell } +a\,\left(T^0_{ij}\,\delta _{k\ell }+T^0_{k\ell }\,\delta _{ij}\right) +b\,\left(T^0_{ik}\,\delta _{j\ell }+T^0_{jk}\,\delta _{i\ell }+T^0_{i\ell }\,\delta _{jk}+T^0_{j\ell }\,\delta _{ik}\right) , \end{equation} (27)where the elements Γijkℓ are given by (8), and for arbitrary values of the quantities a and b, which are allowed to vary spatially. We now turn to the task of determining suitable values for these quantities. 5.1 Plane-wave solutions In this section we consider plane-wave solutions in a homogeneous medium with induced stress. We consider a plane wave of the form   \begin{equation} {{\bf s}}={{\bf a}}\,\exp ({{\bf k}}\cdot {{\bf x}}-{\omega}\,t) , \end{equation} (28)where a denotes the amplitude, k the wavevector, and ω the angular frequency. Substitution of plane wave (28) in the equation of motion (17) leads to the eigenvalue problem   \begin{equation} {{\bf B}}\cdot {{\bf a}}=c^2\,{{\bf a}}, \end{equation} (29)where c = ω/k denotes the phase speed, and where B denotes the symmetric Christoffel tensor with elements (Dahlen 1972b; Dahlen & Tromp 1998, eq. 3.156)   \begin{equation} \rho B_{j\ell }=\Lambda _{ijk\ell }\,\hat{k}_i\,\hat{k}_k . \end{equation} (30)The elements Λijkℓ are related to the elements Ξijkℓ via (Dahlen & Tromp 1998, eq. 3.122)   \begin{equation} \Lambda _{ijk\ell }=\Xi _{ijk\ell }+T^0_{ik}\,\delta _{j\ell } , \end{equation} (31)and are used to relate the first Piola-Kirchhoff stress defined in eq. (19) to the displacement gradient, that is,   \begin{equation} {{\bf T}}^{\rm PK1} = {\boldsymbol \Lambda} :{\boldsymbol \nabla} {{\bf s}}. \end{equation} (32)More specifically, using eqs (8), (13), (27), and (31), we have   \begin{eqnarray} \rho B_{j\ell } &=& \left(\kappa +{\textstyle\frac{1}{3}}\,\mu \right)\,\hat{k}_j\,\hat{k}_\ell +\mu \,\delta _{j\ell } +\gamma _{ijk\ell }\,\hat{k}_i\,\hat{k}_k -2\,(a+b)\,p^0\,\hat{k}_j\,\hat{k}_\ell -(1+2\,b)\,p^0\,\delta _{j\ell } + a\,\tau ^0_{ij}\,\hat{k}_i\,\hat{k}_\ell +a\,\tau ^0_{k\ell }\,\hat{k}_j\,\hat{k}_k +(1+b)\,\tau ^0_{ik}\,\hat{k}_i\,\hat{k}_k\,\delta _{j\ell }\nonumber\\ && +\,b\,\tau ^0_{jk}\,\hat{k}_\ell \,\hat{k}_k +b\,\tau ^0_{i\ell }\,\hat{k}_i\,\hat{k}_j +b\,\tau ^0_{j\ell } . \end{eqnarray} (33)Here $$\hat{k}_i$$ denotes an element of the unit wavevector, $$\hat{{{\bf k}}}={{\bf k}}/k$$. Let us define   \begin{equation} \mu ^{\prime }\equiv -1-2\,b , \end{equation} (34)  \begin{equation} \kappa ^{\prime }\equiv \frac{1}{3}-2\,a-\frac{4}{3}\,b . \end{equation} (35)Then we find that the Christoffel tensor has elements determined by   \begin{eqnarray} \rho B_{j\ell } \ & = & \ \left[(\kappa +\kappa ^{\prime }\,p^0)+{\textstyle\frac{1}{3}}\,(\mu +\mu ^{\prime }\,p^0)\right]\,\hat{k}_j\,\hat{k}_\ell +(\mu +\mu ^{\prime }\,p^0)\,\delta _{j\ell } +\gamma _{ijk\ell }\,\hat{k}_i\,\hat{k}_k - {\textstyle\frac{1}{2}}\,\left(\kappa ^{\prime }+{\textstyle\frac{1}{3}}\,\mu ^{\prime }\right)\,\left(\tau ^0_{ij}\,\hat{k}_i\,\hat{k}_\ell +\tau ^0_{k\ell }\,\hat{k}_j\,\hat{k}_k\right)\nonumber\\ && +\,{\textstyle\frac{1}{2}}(1-\mu ^{\prime })\,\tau ^0_{ik}\,\hat{k}_i\,\hat{k}_k\,\delta _{j\ell } -{\textstyle\frac{1}{2}}(1+\mu ^{\prime })\,\tau ^0_{j\ell } . \end{eqnarray} (36)The terms involving κ + κ΄ p0 and μ + μ΄ p0 enable us to identify κ΄ and μ΄ as the adiabatic pressure derivatives of the bulk and shear moduli, respectively (Nur & Simmons 1969). In fact, it is this desirable result that motivated definitions (34) & (35). We may regard κ and μ as spatially variable elastic moduli at zero pressure, whereas κ΄ and μ΄ represent their spatially variable adiabatic pressure derivatives. Because B is a symmetric positive-definite tensor, the eigenvalue problem (29) has three positive eigenvalues, c2, and associated orthogonal eigenvectors, a. To first-order in the induced stress, assuming the medium before inducing stress is isotropic, that is, γijkℓ = 0, the compressional wave speed is determined by   \begin{eqnarray} \rho \,c_P^2=&(\kappa +\kappa ^{\prime }\,p^0)+\frac{4}{3}\,(\mu +\mu ^{\prime }\,p^0)-(\kappa ^{\prime }+\frac{4}{3}\,\mu ^{\prime })\,\hat{{{\bf k}}}^0\cdot \,{\boldsymbol \tau} ^0\cdot \hat{{{\bf k}}}^0 , \end{eqnarray} (37)and the two shear wave speeds are determined by   \begin{eqnarray} \rho \,c^2_{{\rm S}_{1,2}}=&(\mu +\mu ^{\prime }\,p^0)+\frac{1}{2}\,(1-\mu ^{\prime })\,\hat{{{\bf k}}}^0\cdot \,{\boldsymbol \tau} ^0\cdot \hat{{{\bf k}}}^0-\frac{1}{2}\,(1+\mu ^{\prime })\,\hat{{{\bf a}}}^0_{1,2}\cdot {\boldsymbol \tau} ^0\cdot \hat{{{\bf a}}}^0_{1,2} . \end{eqnarray} (38)Here $$\hat{{{\bf k}}}^0$$ denotes the unit wavevector prior to inducing stress, and $$\hat{{{\bf a}}}^0_{1,2}$$ the unit shear wave polarization directions prior to inducing stress. Note that $$\hat{{{\bf k}}}^0\cdot \hat{{{\bf a}}}^0_{1,2}=0$$. We conclude that definitions (34) and (35) lead to a desirable dependence of the compressional and shear wave speeds (37) & (38) on induced pressure, p0. We observe that the induced deviatoric stress, $${\boldsymbol \tau} ^0$$, manifests itself as an anisotropy in seismic wave speeds. In global seismology, the quantities a and b are selected such that the seismic wave speeds (37) and (38) are independent of the hydrostatic pre-stress, p0. This amounts to setting κ΄ = 0 and μ΄ = 0, and corresponds to the choice a = −b = 1/2. With this choice, the wave speeds (37) and (38) are in agreement with those of Dahlen (1972b, eqs 22 and 25). Note, in particular, that with this choice the compressional wave speed is independent of the induced stress. In this paper, we argue that, rather than choosing the parameters a and b, the parameters μ΄ and κ΄ should be recognized as the derivatives of the moduli with regards to hydrostatic pressure. With this interpretation, these parameters are not chosen, but rather determined by the material properties. For an S wave travelling in the z-direction there are two shear wave speeds, namely, for the wave polarized in the $$\hat{{{\bf x}}}$$ direction   \begin{equation} c^2_{{\rm S}_1}=\beta ^2+{\textstyle\frac{1}{2}}\,[(1-\mu ^{\prime })\,\tau ^0_{zz}-(1+\mu ^{\prime })\tau ^0_{xx}]/\rho , \end{equation} (39)and for the wave polarized in the $$\hat{{{\bf y}}}$$ direction   \begin{equation} c^2_{{\rm S}_2}=\beta ^2+{\textstyle\frac{1}{2}}\,[(1-\mu ^{\prime })\,\tau ^0_{zz}-(1+\mu ^{\prime })\tau ^0_{yy}]/\rho , \end{equation} (40)where $$\beta =\sqrt{\mu /\rho }$$ denotes the unperturbed shear wave speed. These polarization-dependent differences in shear wave speed lead to shear wave splitting, such that the split time between the arrival polarized in the x direction and the arrival polarized in the y direction is determined by   \begin{eqnarray} \delta T \ & = & \ \frac{1}{4}\int (1+\mu ^{\prime })(\tau ^0_{xx}-\tau ^0_{yy})/(\rho \,\beta ^3)\,{\mathrm{d}}z \nonumber\\ & = & \ \frac{1}{4}\int (1+\mu ^{\prime })(\tau ^0_{xx}/\mu -\tau ^0_{yy}/\mu )\,{\mathrm{d}}z/\beta . \end{eqnarray} (41)This expression provides a direct relationship between the split time and the induced deviatoric stress. 5.2 Constitutive relationship Definitions (34) and (35), which lead to a desirable dependence of the compressional and shear wave speeds on induced pressure via (37) and (38), provide the sought values of the quantities a and b, namely,   \begin{equation} b=-\frac{1}{2}\,(1+\mu ^{\prime }) , \end{equation} (42)  \begin{equation} a=\frac{1}{2}\,\left(1-\kappa ^{\prime }+\frac{2}{3}\,\mu ^{\prime }\right) . \end{equation} (43)This implies that eq. (27) becomes   \begin{equation} \Xi _{ijk\ell }=\Gamma _{ijk\ell } +{\textstyle\frac{1}{2}}\,\left(1-\kappa ^{\prime }+{\textstyle\frac{2}{3}}\,\mu ^{\prime }\right)\,\left(T^0_{ij}\,\delta _{k\ell }+T^0_{k\ell }\,\delta _{ij}\right) -{\textstyle\frac{1}{2}}\,(1+\mu ^{\prime })\,\left(T^0_{ik}\,\delta _{j\ell }+T^0_{jk}\,\delta _{i\ell }+T^0_{i\ell }\,\delta _{jk}+T^0_{j\ell }\,\delta _{ik}\right) . \end{equation} (44) As discussed in Section 4, seismic wave propagation under induced stress is described in terms of the incremental Lagrangian Cauchy stress, TL1, which is related to the incremental second Piola–Kirchhoff stress, TSK1, via (Dahlen & Tromp 1998, eqs 3.37)   \begin{equation} {{\bf T}}^{\rm L1}={{\bf T}}^{\rm SK1}-{{\bf T}}^0\,({\boldsymbol \nabla} \cdot {{\bf s}})+({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0+{{\bf T}}^0\cdot {\boldsymbol \nabla} {{\bf s}}. \end{equation} (45)Upon using eqs (25) and (44) in eq. (45), we find that the incremental Lagrangian Cauchy stress has components   \begin{eqnarray} T^{\rm L1}_{ij} & = & \left[(\kappa +\kappa ^{\prime }\,p^0)-{\textstyle\frac{2}{3}}\,(\mu +\mu ^{\prime }\,p^0)\right]\,\delta _{ij}\,\epsilon _{kk} +2\,(\mu +\mu ^{\prime }\,p^0)\,\epsilon _{ij} +\gamma _{ijk\ell }\,\epsilon _{k\ell } + {\textstyle\frac{1}{2}}\,\left(1-\kappa ^{\prime }+{\textstyle\frac{2}{3}}\,\mu ^{\prime }\right)\,\tau ^0_{k\ell }\,\epsilon _{k\ell }\,\delta _{ij} -{\textstyle\frac{1}{2}}\,\left(1+\kappa ^{\prime }-{\textstyle\frac{2}{3}}\,\mu ^{\prime }\right)\,\tau ^0_{ij}\,\epsilon _{kk} \nonumber\\ &&-\,\mu ^{\prime }\,\left(\tau ^0_{ik}\,\epsilon _{kj} +\tau ^0_{jk}\,\epsilon _{ki}\right) +\omega _{ik}\,\tau ^0_{kj}-\tau ^0_{ik}\,\omega _{kj} . \end{eqnarray} (46)In the absence of induced deviatoric stress, $$\tau ^0_{ij}=0$$, we have   \begin{equation} T^{\rm L1}_{ij} = \left[(\kappa +\kappa ^{\prime }\,p^0)-{\textstyle\frac{2}{3}}\,(\mu +\mu ^{\prime }\,p^0)\right]\,\delta _{ij}\,\epsilon _{kk} +2\,(\mu +\mu ^{\prime }\,p^0)\,\epsilon _{ij} +\gamma _{ijk\ell }\,\epsilon _{k\ell } , \end{equation} (47)and here again we clearly recognize the contributions κ + κ΄ p0 and μ + μ΄ p0 as capturing the pressure dependence of the isotropic moduli. In terms of the stress tensor prior to inducing stress, given by eq. (4), we may write   \begin{equation} {{\bf T}}^{\rm L1} = {{\bf T}}+\Delta {{\bf T}}, \end{equation} (48)where   \begin{eqnarray} \Delta {{\bf T}} &= & \overbrace{ \left[\left(\kappa ^{\prime }-\textstyle{\frac{2}{3}}\,\mu ^{\prime }\right)\,{\rm tr}({\boldsymbol \epsilon} )\,{{\bf I}}+2\,\mu ^{\prime }\,{\boldsymbol \epsilon}\right ]\,p^0 }^{\rm {{strain-dependent\, stress\, (due\, to\, induced\, pressure)}}} + \underbrace{ \textstyle{\frac{1}{2}}\,[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}-{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0] -\textstyle{\frac{1}{2}}\,\left(\kappa ^{\prime }-\textstyle{\frac{2}{3}}\,\mu ^{\prime }\right)\,[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}+{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0] -\mu ^{\prime }\,({\boldsymbol \tau} ^0\cdot {\boldsymbol \epsilon} +{\boldsymbol \epsilon} \cdot {\boldsymbol \tau} ^0) }_{{\rm {strain-dependent\, stress\, (due\, to\, induced\, deviatoric\, stress)}}} \nonumber\\ &&+ \,\underbrace{ \boldsymbol \omega \cdot {\boldsymbol \tau} ^0 -{\boldsymbol \tau} ^0\cdot \boldsymbol \omega }_{\rm{rotated \,induced\, stress}}. \end{eqnarray} (49)We conclude that, even in an initially isotropic medium, induced stress manifests itself as an anisotropy in seismic wave speeds. Part of this involves a strain-dependent stress—with distinct contributions from the induced pressure and from the induced deviatoric stress—and another part involves a rotation of the induced stress. It is interesting to note that eq. (49) may be expressed also in terms of the full induced stress tensor, T0, as   \begin{eqnarray} \Delta {{\bf T}} &=& -\textstyle{\frac{1}{2}}\,\left(\kappa ^{\prime }-\textstyle{\frac{2}{3}}\,\mu ^{\prime }\right)\,[({{\bf T}}^0:{\boldsymbol \epsilon} )\,{{\bf I}}+{\rm tr}({\boldsymbol \epsilon} )\,{{\bf T}}^0] -\mu ^{\prime }\,({{\bf T}}^0\cdot {\boldsymbol \epsilon} +{\boldsymbol \epsilon} \cdot {{\bf T}}^0) + \textstyle{\frac{1}{2}}\,[({{\bf T}}^0:{\boldsymbol \epsilon} )\,{{\bf I}}-{\rm tr}({\boldsymbol \epsilon} )\,{{\bf T}}^0] + {\boldsymbol \omega} \cdot {{\bf T}}^0 -{{\bf T}}^0\cdot {\boldsymbol \omega} . \end{eqnarray} (50) We end this section by noting that in global seismology the quantities a and b are selected such that the incremental Lagrangian Cauchy stress (48) is independent of the hydrostatic pre-stress, p0, thereby also rendering the seismic wave speeds (37) and (38) independent of the hydrostatic pre-stress, as discussed previously. This amounts to setting κ΄ = 0 and μ΄ = 0 in eq. (49), implying via eqs (42) and (43) that a = −b = 1/2, and resulting in (Dahlen & Tromp 1998, eq. 3.145)   \begin{eqnarray} \Delta {{\bf T}}\ = & \ \frac{1}{2}\,[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}-{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0] +{\boldsymbol \omega} \cdot {\boldsymbol \tau} ^0 -{\boldsymbol \tau} ^0\cdot {\boldsymbol \omega} . \end{eqnarray} (51)Motivated by expressions (37) and (38) for the compressional and shear wave speeds, perhaps a better strategy is to select κ + κ΄ p0 and μ + μ΄ p0 as ‘the’ bulk and shear moduli in a pre-stressed Earth model, thereby building the pressure dependence into the moduli and capturing the effects of deviatoric pre-stress in terms of   \begin{eqnarray} \Delta {{\bf T}} = & \frac{1}{2}\,[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}-{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0] + {\boldsymbol \omega} \cdot {\boldsymbol \tau} ^0 -{\boldsymbol \tau} ^0\cdot {\boldsymbol \omega} -\frac{1}{2}\,\left(\kappa ^{\prime }-\frac{2}{3}\,\mu ^{\prime }\right)\,\left[({\boldsymbol \tau} ^0:{\boldsymbol \epsilon} )\,{{\bf I}}+{\rm tr}({\boldsymbol \epsilon} )\,{\boldsymbol \tau} ^0\right] -\mu ^{\prime }\,({\boldsymbol \tau} ^0\cdot {\boldsymbol \epsilon} +{\boldsymbol \epsilon} \cdot {\boldsymbol \tau} ^0) . \end{eqnarray} (52)Compared to eq. (51), the additional terms are generally not negligible, because typical values of κ΄ and μ΄ are of order unity (e.g. Stacey 1992). Arguably, such terms could be captured by the purely anisotropic elements γijkℓ, but such an approach fails to reflect the full influence of the deviatoric pre-stress. 6 BORN APPROXIMATION To compare and contrast the equation of motion and boundary condition before inducing stress, eqs (1) and (2), with those after inducing stress, eqs (17) and (18), we use eq. (48) to rewrite the latter as   \begin{equation} \rho \,\partial _t^2{{\bf s}}-{\boldsymbol \nabla} \cdot {{\bf T}}={\boldsymbol \nabla} \cdot \Delta {{\bf T}}+{\boldsymbol \nabla} [{{\bf s}}\cdot ({\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0)]-{\boldsymbol \nabla} \cdot ({{\bf s}}\cdot {\boldsymbol \nabla} {\boldsymbol \tau} ^0) , \end{equation} (53)and   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}]_-^+=-[ \hat{{{\bf n}}}\cdot {{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-\hat{{{\bf n}}}\cdot ({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0 +\hat{{{\bf n}}}\cdot \Delta {{\bf T}}]_-^+ . \end{equation} (54)The terms on the right-hand sides of these equations capture effects due to induced stress. If the induced stress is accompanied by changes in density, δρ, and the elastic parameters, $$\delta {\boldsymbol \Gamma}$$, then the modified equations take the form   \begin{equation} \rho \,\partial _t^2{{\bf s}}-{\boldsymbol \nabla} \cdot {{\bf T}}= -\delta \rho \,\partial _t^2{{\bf s}}+{\boldsymbol \nabla} \cdot \delta {{\bf T}}+{\boldsymbol \nabla} \cdot \Delta {{\bf T}}+{\boldsymbol \nabla} [{{\bf s}}\cdot ({\boldsymbol \nabla} \cdot {\boldsymbol \tau} ^0)]-{\boldsymbol \nabla} \cdot ({{\bf s}}\cdot {\boldsymbol \nabla} {\boldsymbol \tau} ^0) , \end{equation} (55)and   \begin{equation} [\hat{{{\bf n}}}\cdot {{\bf T}}]_-^+=-[ \hat{{{\bf n}}}\cdot \Delta {{\bf T}}+\hat{{{\bf n}}}\cdot {{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-\hat{{{\bf n}}}\cdot ({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0 +\hat{{{\bf n}}}\cdot \delta {{\bf T}}]_-^+ , \end{equation} (56)where   \begin{equation} \delta {{\bf T}}=\delta {\boldsymbol \Gamma} :{\boldsymbol \epsilon} . \end{equation} (57)Additional effects associated with changes in layer geometry are captured by perturbations in boundary locations. Such topographic perturbations on internal discontinuities and the free surface modify the boundary conditions. The theoretical implications are well known (see e.g. Dahlen 2005; Tromp et al.2005) and readily accommodated. For the sake of brevity, in this paper we focus mainly on the direct effects of induced stress, but in Section 7 we also summarize effects associated with perturbations in density, the elastic tensor, and boundary topography. If the wavefield prior to inducing stress, s, is governed by the equation of motion (1) subject to boundary conditions (2), and the wavefield after inducing stress, s + δs, is governed by the equation of motion (53) subject to boundary conditions (54), then, invoking the Born approximation, the perturbed wavefield, δs, due to the induced stress is determined by   \begin{eqnarray} \delta s_i({{\bf x}},t) &=& \ \int _{-\infty }^t\int _VG_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime })\,\partial ^{\prime }_k[ \Delta T_{kj}({{\bf x}}^{\prime },t^{\prime }) +s_\ell ({{\bf x}}^{\prime },t^{\prime })\,\partial ^{\prime }_m \tau _{m\ell }^0({{\bf x}}^{\prime },t^{\prime })\,\delta _{kj} -s_\ell ({{\bf x}}^{\prime },t^{\prime })\,\partial ^{\prime }_\ell \tau _{kj}^0({{\bf x}}^{\prime },t^{\prime }) ]\,{\mathrm{d}}^3{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime}\nonumber\\ && + \int _{-\infty }^t\int _\Sigma G_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime })\, \hat{n}_k({{\bf x}}^{\prime })\,[ \Delta T_{kj}({{\bf x}}^{\prime },t^{\prime }) +T_{kj}^0({{\bf x}}^{\prime },t^{\prime })\,\partial ^{\prime }_\ell s_\ell ({{\bf x}}^{\prime },t^{\prime }) -(\partial ^{\prime }_\ell s_k({{\bf x}}^{\prime },t^{\prime }))\,T_{\ell j}^0({{\bf x}}^{\prime },t^{\prime }) ]_-^+\,{\mathrm{d}}^2{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime}. \end{eqnarray} (58)Here V denotes the model volume, Σ the collection of all discontinuities, including internal discontinuities and the free surface, and Gij(x, x΄; t − t΄) the Green function for the medium prior to inducing stress. A partial derivative with respect to the primed coordinates, x΄, is denoted by $$\partial ^{\prime }_i$$. Additional terms associated with perturbations in density, δρ, elements of the elastic tensor, δΓijkℓ, and boundary topography, δh, may be found in Tromp et al. (2005, eqs 3 and 21). After a lengthy, tedious calculation, this result may be rewritten in the form   \begin{eqnarray} \delta s_i({{\bf x}},t) & = & -\int _{-\infty }^t\int _V\partial ^{\prime }_kG_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime })\, [\Delta T_{kj}({{\bf x}}^{\prime },t^{\prime }) -(\partial ^{\prime }_ms_\ell ({{\bf x}}^{\prime },t^{\prime }))\,\tau _{m\ell }^0({{\bf x}}^{\prime },t^{\prime })\,\delta _{kj} +(\partial ^{\prime }_\ell s_\ell ({{\bf x}}^{\prime },t^{\prime }))\,\tau _{kj}^0({{\bf x}}^{\prime },t^{\prime })]\,{\mathrm{d}}^3{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime } \nonumber\\ &&+ \,\int _{-\infty }^t\int _V(\partial ^{\prime }_m\partial ^{\prime }_jG_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime }))\, s_\ell ({{\bf x}}^{\prime },t^{\prime })\,\tau _{m\ell }^0({{\bf x}}^{\prime },t^{\prime }) \,{\mathrm{d}}^3{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime } -\int _{-\infty }^t\int _V(\partial ^{\prime }_\ell \partial ^{\prime }_kG_{ij}({{\bf x}},{{\bf x}}^{\prime };t-t^{\prime }))\, s_\ell ({{\bf x}}^{\prime },t^{\prime })\,\tau _{kj}^0({{\bf x}}^{\prime },t^{\prime }) \,{\mathrm{d}}^3{{\bf x}}^{\prime }\,{\mathrm{d}}t^{\prime } . \end{eqnarray} (59)Note that the boundary integral has been eliminated, and, in particular, that eq. (59) contains no derivatives of the induced deviatoric stress. Result (59) forms the basis of the calculation of Fréchet derivatives in the next section. 7 KERNELS In this section we determine induced deviatoric stress sensitivity kernels or Fréchet derivatives. Suppose we record a set of three-component data with an array of Nr receivers. We denote the three component time series in this data set by di(xr, t), i = 1, 2, 3, r = 1, …, Nr. After a time lapse, ΔT, we record a second data set $$d_i^{\Delta T}({{\bf x}}_r,t)$$, i = 1, 2, 3, r = 1, …, Nr. We seek to find model perturbations such that the related simulated wavefield perturbations, δsi, capture the change in the data. This is accomplished by defining the misfit function   \begin{equation} \chi =\frac{1}{2}\,\sum _{r=1}^{N_r}\sum _{i=1}^3\int ||d_i({{\bf x}}_r,t)+\delta s_i({{\bf x}}_r,t)-d_i^{\Delta T}({{\bf x}}_r,t)||^2\,{\mathrm{d}}t , \end{equation} (60)such that   \begin{equation} \delta \chi =\sum _{r=1}^{N_r}\sum _{i=1}^3\int [d_i({{\bf x}}_r,t)-d_i^{\Delta T}({{\bf x}}_r,t)]\,\delta s_i({{\bf x}}_r,t)\,{\mathrm{d}}t . \end{equation} (61)Using the Born expression determined in Section 6, the derivative of the misfit function with respect to the induced deviatoric stress may be expressed as   \begin{equation} \delta \chi =\int _V p^0\,K_{p^0}+{\boldsymbol \tau} ^0:{{\bf K}}_{{\bf \tau }^0}\,{\mathrm{d}}V , \end{equation} (62)where   \begin{equation} K_{p^0}=-\int _0^T[\kappa ^{\prime }({{\bf x}})\,{\boldsymbol \nabla} \cdot {{\bf s}}({{\bf x}},t)\,{\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t) +2\,\mu ^{\prime }({{\bf x}})\,{{\bf D}}({{\bf x}},t):{{\bf D}}^\dagger ({{\bf x}},T-t)]\,{\mathrm{d}}t , \end{equation} (63)  \begin{eqnarray} {{\bf K}}_{{\bf \tau }^0}({{\bf x}}) &= & -\int _0^T\left\lbrace \frac{1}{2}\,{\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t)\,{{\bf D}}({{\bf x}},t) -\frac{1}{2}\,{\boldsymbol \nabla} \cdot {{\bf s}}({{\bf x}},t)\,{{\bf D}}^\dagger ({{\bf x}},T-t)\right. \nonumber\\ &&\left. -\,\frac{1}{2}\left[\kappa ^{\prime }({{\bf x}})+\frac{2}{3}\,\mu ^{\prime }({{\bf x}})\right][{\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t)\,{{\bf D}}({{\bf x}},t)+{\boldsymbol \nabla} \cdot {{\bf s}}({{\bf x}},t)\,{{\bf D}}^\dagger ({{\bf x}},T-t)]\right. \nonumber\\ && \left. -\,\mu ^{\prime }({{\bf x}})\,[{{\bf D}}({{\bf x}},t)\cdot {{\bf D}}^\dagger ({{\bf x}},T-t)+{{\bf D}}^\dagger ({{\bf x}},T-t)\cdot {{\bf D}}({{\bf x}},t)] -{\boldsymbol \omega} ({{\bf x}},t)\cdot {{\bf D}}^\dagger ({{\bf x}},T-t) +{{\bf D}}^\dagger ({{\bf x}},T-t)\cdot {\boldsymbol \omega} ({{\bf x}},t)\right. \nonumber\\ && \left. +\,{{\bf s}}({{\bf x}},t)\cdot {\boldsymbol \nabla} {\boldsymbol \epsilon} ^\dagger ({{\bf x}},T-t) -\frac{1}{2}\,{{\bf s}}({{\bf x}},t)\,{\boldsymbol \nabla} {\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t) -\frac{1}{2}\,[{\boldsymbol \nabla} {\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t)]\,{{\bf s}}({{\bf x}},t) \right\rbrace \,{\mathrm{d}}t . \end{eqnarray} (64)Here T denotes the record length. The adjoint wavefield, s†, satisfies the same wave equation as the regular wavefield, s, but is generated by the adjoint source   \begin{equation} f_i^\dagger ({{\bf x}},t)=\sum _{r=1}^{N_r}\left[d_i({{\bf x}}_r,T-t)-d_i^{\Delta T}({{\bf x}}_r,T-t)\right]\,\delta ({{\bf x}}-{{\bf x}}_r) . \end{equation} (65)This equation illustrates how time-lapse data, which are simultaneously insert at all receiver locations in reverse time, give rise to the adjoint wavefield. The strain deviator and its adjoint are denoted by D and D†, respectively. Note that the kernel $${{\bf K}}_{{\bf \tau }^0}$$ has zero trace, as expected, and that the second derivative of the adjoint wavefield, $${\boldsymbol \nabla} {\boldsymbol \nabla} {{\bf s}}^\dagger$$, is required for its calculation. Effects due to changes in density, δρ, elastic parameters, $$\delta \boldsymbol \Gamma$$, and topography on solid-solid boundaries, δh, are readily incorporated (Tromp et al.2005). In this case the derivative of the misfit function is modified as   \begin{equation} \delta \chi =\int _V (\delta \rho \,K_\rho +\delta {\boldsymbol \Gamma} :: {{\bf K}}_{{{\boldsymbol \Gamma} }}+p^0\,K_{p^0}+{\boldsymbol \tau} ^0:{{\bf K}}_{{\bf \tau }^0})\,{\mathrm{d}}V+\int _\Sigma \delta h\,K_h\,{\mathrm{d}}\,\Sigma , \end{equation} (66)where   \begin{equation} K_\rho =\int _0^T\partial _t{{\bf s}}^\dagger \cdot \partial _t{{\bf s}}\,{\mathrm{d}}t , \end{equation} (67)   \begin{equation} {{\bf K}}_{{{\boldsymbol \Gamma} }}=-\frac{1}{2}\int _0^T({\boldsymbol \epsilon} ^\dagger \,{\boldsymbol \epsilon} +{\boldsymbol \epsilon} \,{\boldsymbol \epsilon} ^\dagger )\,{\mathrm{d}}t , \end{equation} (68)  \begin{equation} K_h=-\int _0^T\left[\rho \,\partial _t{{\bf s}}^\dagger \cdot \partial _t{{\bf s}}-{\boldsymbol \epsilon} ^\dagger :{\boldsymbol \Gamma} :{\boldsymbol \epsilon} +\hat{{{\bf n}}}\,\partial _n{{\bf s}}^\dagger :{\boldsymbol \Gamma} :{\boldsymbol \epsilon} +\hat{{{\bf n}}}\,\partial _n{{\bf s}}:{\boldsymbol \Gamma} :{\boldsymbol \epsilon} ^\dagger \right]_-^+\,{\mathrm{d}}t . \end{equation} (69)Here $$\partial _n=\hat{{{\bf n}}}\cdot {\boldsymbol \nabla}$$, and the dependence of the adjoint wavefield on time T − t and the forward wavefield on time t has been omitted for brevity. Topographic perturbations on fluid-solid boundaries may also be accommodated (Dahlen 2005; Tromp et al.2005), but are omitted for brevity. For isotropic models, kernel KΓ defined in eq. (68) is replaced by kernels for the bulk and shear moduli (Tromp et al.2005), namely   \begin{equation} K_\kappa =-\int _0^T{\boldsymbol \nabla} \cdot {{\bf s}}({{\bf x}},t)\,{\boldsymbol \nabla} \cdot {{\bf s}}^\dagger ({{\bf x}},T-t) \,{\mathrm{d}}t , \end{equation} (70)  \begin{equation} K_\mu =-2\int _0^T{{\bf D}}({{\bf x}},t):{{\bf D}}^\dagger ({{\bf x}},T-t)]\,{\mathrm{d}}t . \end{equation} (71)Upon comparing the isotropic kernels (70) and (71) with kernel (63) for the induced pressure, we conclude that induced pressure trades off with the moduli, an unfortunate fact to be expected from expressions (37) and (38) for the wave speeds. As mentioned previously, for this reason it makes sense in both active and passive seismology to build the induced pressure dependence into the moduli by selecting κ + κ΄ p0 and μ + μ΄ p0 as model parameters, thereby eliminating induced pressure as an unknown parameter. However, one can imagine engineering applications in which the elastic moduli, κ & μ, and their pressure derivatives, κ΄ & μ΄, are known, and one seeks to determine induced pressure variations, for example, based on time-lapse ultrasonics. In that case kernel (63) may prove to be a useful tool. We note that the formulation of an iterative inverse problem based on the Fréchet derivatives presented in this section is non-trivial, because updates in induced stress need to continue to satisfy the equilibrium eq. (9), which, as previously mentioned, constrains three out of six elements of the induced stress. This might be accomplished by making use of the ‘Beltrami stress function’ (Gurtin 1963), or based on the orthogonal decompositions of stress tensors introduced by Al-Attar (2010). This important topic is beyond the scope of this paper. 8 2-D SH WAVES 8.1 Equation of motion and kernels Consider a 2-D SH wavefield polarized in the y-direction with displacement field s(x, z, t). The associated scalar wave equation is   \begin{equation} \rho \,\partial _t^2s=\partial _x(\mu \,\partial _x s)+\partial _z(\mu \,\partial _z s) , \end{equation} (72)where ρ = ρ(x, z) denotes mass density and μ = μ(x, z) the shear modulus. Upon inducing a purely deviatoric stress of the form $$\tau ^0_{xz}=\tau$$, the 2-D SH equation of motion becomes   \begin{eqnarray} \rho \,\partial _t^2s &=\partial _x\left(\mu \,\partial _x s+\frac{1}{4}\,\tau \,\partial _zs\right)+\partial _z\left(\mu \,\partial _z s+\frac{1}{4}\,\tau \,\partial _xs\right) . \end{eqnarray} (73)The misfit gradient may be written in the form   \begin{equation} \delta \chi =\int _V (\delta \ln \rho \,K_\rho +\delta \ln \mu \,K_\mu +\tau \,K_{\tau })\,{\mathrm{d}}^2{{\bf x}}+\int _S \delta h\,K_h\,{\mathrm{d}}^1{{\bf x}}, \end{equation} (74)where   \begin{equation} K_{\rho }= \int _0^T\rho \,\partial _ts(t)\,\partial _ts^\dagger (T-t)\,{\mathrm{d}}t , \end{equation} (75)  \begin{eqnarray} K_{\mu } =& -\int _0^T\mu \,[ \partial _xs(t)\,\partial _xs^\dagger (T-t)+\partial _zs(t)\,\partial _zs^\dagger (T-t)]\,{\mathrm{d}}t , \end{eqnarray} (76)  \begin{eqnarray} K_{\tau } =& -\frac{1}{4}\,\int _0^T [\partial _zs(t)\,\partial _xs^\dagger (T-t)+\partial _xs(t)\,\partial _zs^\dagger (T-t)]\,{\mathrm{d}}t , \end{eqnarray} (77)  \begin{eqnarray} K_h & = & -\int _0^T[\rho \,\partial _ts(t)\,\partial _ts^\dagger (T-t) -\mu \,\partial _xs(t)\,\partial _xs^\dagger (T-t) -\mu \,\partial _zs(t)\,\partial _zs^\dagger (T-t) \nonumber\\ && +\,\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _ns^\dagger (t)\,\partial _xs(T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial _ns^\dagger (t)\,\partial _xs(T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _ns(t)\,\partial _xs^\dagger (T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial _ns(t)\,\partial _xs^\dagger (T-t) ]_-^+ \,{\mathrm{d}}t . \end{eqnarray} (78)Here $$\partial _n=\hat{{{\bf n}}}\cdot {\boldsymbol \nabla} =\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _x+\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial_z$$. In terms of $$\beta =\sqrt{\mu /\rho }$$ and τ΄ = τ/μ, we have the preferred alternative expression   \begin{equation} \delta \chi =\int _V (\delta \ln \rho \,K^{\prime }_\rho +\delta \ln \beta \,K_\beta +\tau ^{\prime }\,K_{\tau ^{\prime }})\,{\mathrm{d}}^2{{\bf x}}+\int _S \delta \ln h\,K^{\prime }_h\,{\mathrm{d}}^1{{\bf x}}, \end{equation} (79)where   \begin{equation} K^{\prime }_\rho =K_\rho +K_\mu , \end{equation} (80)  \begin{equation} K_\beta =2\,K_\mu , \end{equation} (81)  \begin{eqnarray} K_{\tau ^{\prime }} =& -\frac{1}{4}\,\int _0^T\mu \, [\partial _zs(t)\,\partial _xs^\dagger (T-t)+\partial _xs(t)\,\partial _zs^\dagger (T-t)]\,{\mathrm{d}}t , \end{eqnarray} (82)  \begin{eqnarray} K^{\prime }_h & = & -\int _0^Th[\rho \,\partial _ts(t)\,\partial _ts^\dagger (T-t) -\mu \,\partial _xs(t)\,\partial _xs^\dagger (T-t) -\mu \,\partial _zs(t)\,\partial _zs^\dagger (T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _ns^\dagger (t)\,\partial _xs(T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial _ns^\dagger (t)\,\partial _xs(T-t) \nonumber\\ &&+\,\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf x}}}\,\partial _ns(t)\,\partial _xs^\dagger (T-t) +\mu \,\hat{{{\bf n}}}\cdot \hat{{{\bf z}}}\,\partial _ns(t)\,\partial _xs^\dagger (T-t) ]_-^+ \,{\mathrm{d}}t . \end{eqnarray} (83) 8.2 Examples Illustrating the effects of induced stress on elastic wave propagation for 2-D SH waves is both straightforward and representative of more realistic cases. We use the spectral-element method (Komatitsch & Vilotte 1998; Komatitsch & Tromp 1999; Peter et al.2011) to simulate wave propagation with and without induced stress, and for the calculation of the kernels. The first example is taken from Tromp et al. (2005) and is based on a 2-D homogeneous medium in a 80 km × 200 km box with values of ρ = 2600 kg m−3 and β = 3200 m s−1, corresponding to a shear modulus of μ = 26.6 GPa. The box has a free surface at the top and absorbing boundaries at the bottom and at the sides. The source is a point force in the direction perpendicular to the medium plane with a Gaussian source time function. The wavefield is perfectly circular, before touching any boundaries. If we add a 20 per cent induced stress to this medium, we see (Fig. 1) that the wavefield now becomes slightly elliptical, as predicted in Sections 5 and 5.1. Effectively, the wavefield sees an apparent anisotropic medium. Figure 1. View largeDownload slide Snapshot after 16 s of SH wave propagation in a homogeneous background medium. Without induced stress, the wavefield is perfectly circular (bottom panel). In the presence of induced stress, the wave sees an apparent anisotropic material (middle panel). The difference of the two wavefields (top panel) gives an indication of the magnitude of the stress effect. The wavefield difference is about 10 times smaller than the wavefield itself in the presence of an induced stress of 20 per cent of the shear modulus. The units of the gird are in m and of the colour scale in 10 − 2 μm. Figure 1. View largeDownload slide Snapshot after 16 s of SH wave propagation in a homogeneous background medium. Without induced stress, the wavefield is perfectly circular (bottom panel). In the presence of induced stress, the wave sees an apparent anisotropic material (middle panel). The difference of the two wavefields (top panel) gives an indication of the magnitude of the stress effect. The wavefield difference is about 10 times smaller than the wavefield itself in the presence of an induced stress of 20 per cent of the shear modulus. The units of the gird are in m and of the colour scale in 10 − 2 μm. Keeping the same background model, we now calculate the kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ as specified by eqs (80)–(82). In terms of shape, we see (Fig. 2) that the kernels have the structure we would expect. In particular the diagonal symmetry seen in $$K_{\tau ^{\prime }}$$ is the result of the cross-product of partial derivatives (compare equations 76 and 82). In terms of amplitude, although we plot the kernels on a saturated scale to visualize the details, the β kernel is roughly an order of magnitude larger than both $$K^{\prime }_\rho$$ and $$K_{\tau ^{\prime }}$$. Figure 2. View largeDownload slide Banana-doughnut kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in a homogeneous background medium. The units of the grid are in m and of the colour scale in 10−9 m−2. Figure 2. View largeDownload slide Banana-doughnut kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in a homogeneous background medium. The units of the grid are in m and of the colour scale in 10−9 m−2. To investigate a more representative case for oil & gas exploration, we implement a simplified model for the Groningen gas field in the Netherlands (Kraaijpool & Dost 2013) in a 4 km × 4 km box. The top boundary acts as a free surface and the other 3 boundaries are absorbing. The model consist of seven layers (Table 1), where the gas is contained in the Rotliegend group overlain by a fast salt layer of the Zechstein group. We implement an array of 10 seismometers at the surface and another array of 10 instruments in a borehole. We assume that gas extraction induces an earthquake in the gas layer. The salt layer acts like a strong barrier to wave propagation and only little energy escapes to the surface (Fig. 3). Figure 3. View largeDownload slide Snapshot after 1.5 s of SH wave propagation in the isotropic Groningen model (Table 1). The seismometers are marked by green squares and the event by an orange cross. The model is represented by the grey background, where the darkest colour represents the fastest layer. Most energy of the waveform remains in and below the salt layer, only the first wave being transmitted. Figure 3. View largeDownload slide Snapshot after 1.5 s of SH wave propagation in the isotropic Groningen model (Table 1). The seismometers are marked by green squares and the event by an orange cross. The model is represented by the grey background, where the darkest colour represents the fastest layer. Most energy of the waveform remains in and below the salt layer, only the first wave being transmitted. Table 1. Model parameters for the simplified Groningen gas field. Unit  Depth (m)  β (m s−1)  ρ (kg m−3)  North Sea supergroup  0  372  1960  Chalk group  828  756  2300  Rijnland group  1713  663  2650  Upper and lower germanic Trias groups  1820  1427  2600  Zechstein group  1925  2725  2100  Upper Rotliegend group  2866  1184  2550  Limburg group  3126  1695  2700  Unit  Depth (m)  β (m s−1)  ρ (kg m−3)  North Sea supergroup  0  372  1960  Chalk group  828  756  2300  Rijnland group  1713  663  2650  Upper and lower germanic Trias groups  1820  1427  2600  Zechstein group  1925  2725  2100  Upper Rotliegend group  2866  1184  2550  Limburg group  3126  1695  2700  View Large An important question is whether induced stress is measurable by seismic time-lapse experiments. Making simplifying assumptions, Hatchell & Bourne (2005) related measured time-lapse time shifts to strain changes induced by extraction, showing values of 10–20 milliseconds for several off-shore reservoirs. In our simulations, we keep the elastic constants fixed and only add induced stresses due to gas extraction, which is then equivalent to induced strains. Dempsey & Suckale (2017) estimated induced stresses in the Groningen gas field and found values of the order of 75 MPa in the gas layer, corresponding to 2 per cent of the shear modulus in that layer. We systematically vary the induced stress up to 20 per cent of the shear modulus in the Rotliegend and measure cross-correlation time shifts compared to the case with no induced stress for a station at the surface vertically above the event (S1), the most extreme right station at the surface (S10), and the borehole station just below the gas layer (B18), which receives the bulk of the seismic energy. No time shifts are detected at S1, which receives a direct wave that travelled the least in the stressed medium. S10 and B18 both see maximum time shifts of 1.5 and 3.0 ms, respectively, for stresses at 20 per cent of the shear modulus. These time shifts are an order of magnitude smaller than those reported by Hatchell & Bourne (2005). The time shifts of a few tens of milliseconds seen by Hatchell & Bourne (2005) are thus more likely due to boundary or velocity changes induced by the stress changes, rather than the stress changes themselves, as the latter would require substantial stresses to achieve such large time shifts. Looking at the waveforms however, we notice that the seismic coda shows more significant differences while the first arriving waves present tiny time shifts (Fig. 4). To quantify this, we evaluate the average relative waveform difference δΦ defined as   \begin{equation} \delta \Phi = \frac{1}{2}\,\log \left(\frac{\int d_r^2\,{\mathrm{d}}t}{\int d_0^2\,{\mathrm{d}}t}\right) , \end{equation} (84)where d0 is the seismogram in the unstressed medium and dτ the seismogram in the medium with induced stress. Measuring δΦ systematically as a function of induced stress, we see that significant changes in waveforms are noticeable even for modest stresses (Fig. 5). While in this particular setting, some induced stress kernels are small, others, especially in deep boreholes, clearly show that induced stress has a large sensitivity on the waveforms (Figs 6–8). The conclusion of these observations is that it is possible to detect induced stresses of magnitudes we expect in reservoir settings in time lapse experiments by looking at full waveform differences rather than time shifts of first arriving waves. Our 2-D examples are designed to show the order of magnitude of the waveform effects to be expected, and hence our simple waveform difference measure δΦ is sufficient. In real applications, a more rigorous quantification of the waveform difference should be used (e.g. Kristekova et al.2009). We anticipate that in three dimensions, shear wave splitting caused by induced deviatoric stress may also be a detectable. Figure 4. View largeDownload slide Sample seismogram at station B18 (location shown on Fig. 3) in the Groningen medium with a 20 per cent induced stress relative to the shear modulus (black) compared to the seismogram in the unstressed Groningen medium (red). The bottom panel shows the complete seismograms. The top right panel is a zoom into the first arriving wave. The top left panel is a zoom into the coda. All three panels used the same units specified in the bottom panel. Figure 4. View largeDownload slide Sample seismogram at station B18 (location shown on Fig. 3) in the Groningen medium with a 20 per cent induced stress relative to the shear modulus (black) compared to the seismogram in the unstressed Groningen medium (red). The bottom panel shows the complete seismograms. The top right panel is a zoom into the first arriving wave. The top left panel is a zoom into the coda. All three panels used the same units specified in the bottom panel. Figure 5. View largeDownload slide Average relative waveform differences δΦ measured in stressed models with respect to the unstressed Groningen medium. The different curves correspond to waveform changes at different stations located in Fig. 3. Figure 5. View largeDownload slide Average relative waveform differences δΦ measured in stressed models with respect to the unstressed Groningen medium. The different curves correspond to waveform changes at different stations located in Fig. 3. Figure 6. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station S1 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 6. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station S1 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 7. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station S10 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 7. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station S10 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 8. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station B18 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. Figure 8. View largeDownload slide Waveform difference kernels $$K^{\prime }_\rho$$, Kβ and $$K_{\tau ^{\prime }}$$ in the Groningen medium for station B18 (location shown in Fig. 3). The units of the grid are in m and of the colour scale in 10−6 m−2. 9 CONCLUDING REMARKS It is well known that seismic wave speeds in crustal rocks increase with hydrostatic pressure (Birch 1961; Nur & Simmons 1969). Under non-hydrostatic stress, the same wave speeds develop clear characteristics of anisotropy (Nur & Simmons 1969). Mavko et al. (1995), for instance, explained these observations by introducing cracks with prescribed compliances and further assuming a porosity for each crack. Another approach is to use higher-order elasticity (e.g. Johnson & Rasolofosaon 1996; Prioul et al.2004) to explain effective anisotropy induced by deviatoric stress. Johnson & Rasolofosaon (1996) plot squared compressional wave speed $$c_p^2$$ as a function of applied uniaxial stress. Using our expression (37), we find that $$\rho \,c_p^2$$ changes linearly with uniaxial stress P. Recalling Birch’s law (Birch 1961), stating that density is proportional to P-wave speed, we find that $$c_p^2 \sim P^{2/3}$$, similar to the data shown in Johnson & Rasolofosaon (1996). These authors further plot S-wave birefringence data, which present a trend in sin 2θ, where θ is the angle between the propagation direction and the stress direction. Our eq. (38) similarly exhibits such a sin 2θ dependence without the need for third-order elasticity. Prioul et al. (2004) show data where elastic constants vary linearly with applied hydrostatic stress. In our formulation we clearly identify the linear hydrostatic pressure dependence of the elastic constants (eqs 37 and 38). Motivated by the formulation common in global seismology, we show the effects of induced stress on the elastic wave equation and constitutive relation. Without employing higher-order elasticity, our formulation leads to trends observed in measurements based on laboratory data. We also, for the first time, give expressions for induced stress sensitivity kernels, which can be used in adjoint tomography. From an inverse theory perspective, induced pressure trades off with the isotropic moduli, and for this reason it should be built into the definition of ‘the moduli’, thereby eliminating it as an unknown parameter. Induced deviatoric stress results in anisotropic seismic wave speeds, and thus it is possible to image deviatoric stress using time-lapse experiments. Simple 2-D SH experiments suggest waveforms are sufficiently sensitive to stress changes expected in reservoir settings, but we will have to deal with trade-offs between changes in elastic constants, boundaries, and induced deviatoric stress, which all affect the wavefield simultaneously. The theory contains two quantities a and b (see eq. 27), which may be chosen to obtain a particular equation of state. Global seismologists (e.g. Dahlen 1972a; Dahlen & Tromp 1998) prefer to choose a = −b = 1/2, thereby rendering the formulation independent of the hydrostatic pre-stress. An alternative second-order equation of state is based on logarithmic or Hencky strain (Poirier & Tarantola 1998) rather than the Lagrangian strain, and may be recovered from our theory by setting a = 0 and b = −1/2. Our analysis suggests a physical interpretation of the quantities a and b that gives the isotropic bulk and shear moduli an induced pressure dependence of the form κ + κ΄ p0 and μ + μ΄ p0, respectively, where κ and μ denote isotropic elastic moduli at zero pressure, and κ΄ and μ΄ represent their adiabatic pressure derivatives. Note that—because the derivatives κ΄ and μ΄ may be depth-dependent—the quantities κ + κ΄ p0 and μ + μ΄ p0 do not necessarily vary linearly with pressure. For moderate stress changes, Birch-Murnaghan’s second-order equation of state in Eulerian strain is often used (e.g. Stacey 1992). This formulations leads to κ΄ = 4, and together with the empirical relation μ = 0.631 κ − 0.899 P (see Stacey 1992), results in a = −0.95833 and b = −1.3125, values that are quite different from those used in global seismology. It is important to recognize that whereas Dahlen (1972a,b) and Dahlen & Tromp (1998) choose to set a = −b = 1/2, we do not select specific values for a and b. By rewriting a and b in terms of two new parameters, μ΄ and κ΄, we show that the predicted seismic wave speeds, defined by eqs (37) and (38), take on experimentally expected forms when one interprets these new parameters as the adiabatic pressure derivatives of the shear and bulk moduli with respect to pressure. But note also that based on this interpretation one no longer simply ‘chooses’ a and b. Instead, the values of a and b are determined by the pressure derivatives of the shear and bulk moduli. Our theory makes testable predictions given values of the shear and bulk moduli and their pressure derivatives and an applied pre-stress, for example, a directional dependence of the compressional wave speed (37) and the shear wave split time (41). We predict that the compressional wave speed acquires a directional dependence with an amplitude linearly proportional to κ΄ and μ΄, whereas setting a = −b = 1/2—equivalent to setting μ΄ = κ΄ = 0 in our formulation—results in an isotropic P wave. Similarly, we predict that for two materials with the same shear modulus μ but different pressure derivatives μ΄, the split time will be different, whereas setting a = −b = 1/2 would not show any differences in split time. One could imagine lab experiments that could attempt to verify these predictions given experimentally determined values of the moduli μ and κ and their pressure derivatives μ΄ and κ΄ for various materials. Acknowledgements We thank an anonymous reviewer, David Al-Attar, and Editor Martin Mai for helpful comments which helped to improve the manuscript. The open source spectral-element package SPECFEM2D used in this study is freely available via the Computational Infrastructure for Geodynamics (CIG). The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP/2007-2013) grant agreement number 320639 (iGEO). REFERENCES Al-Attar D., 2010. On the parametrisation of equilibrium stress fields in the Earth, Geophys. J. Int. , 181, 567– 576. https://doi.org/10.1111/j.1365-246X.2010.04527.x Birch F., 1961. The velocity of compressional waves in rocks to 10 kilobars: 2, J. geophys. Res. , 66, 2199– 2224. Google Scholar CrossRef Search ADS   Bruner W.M., 1976. Comment on ‘Seismic velocities in dry and saturated cracked solids’ by Richard J. O’Connell and Bernard Budiansky, J. geophys. Res. , 81, 2573– 2576. Google Scholar CrossRef Search ADS   Dahlen F.A., 1972a. Elastic dislocation theory for a self-gravitating elastic configuration with an initial static stress field, Geophys. J. R. astr. Soc. , 28, 357– 383. https://doi.org/10.1111/j.1365-246X.1972.tb06798.x Google Scholar CrossRef Search ADS   Dahlen F.A., 1972b. Elastic velocity anisotropy in the presence of an anisotropic initial stress, Bull. seism. Soc. Am. , 62, 1183– 1193. Dahlen F.A., 2005. Finite-frequency sensitivity kernels for boundary topography perturbations, Geophys. J. Int. , 162, 525– 540. https://doi.org/10.1111/j.1365-246X.2005.02682.x Dahlen F.A., Tromp J., 1998. Theoretical Global Seismology , Princeton University Press. Dempsey D., Suckale J., 2017. Physics-based forecasting of induced seismicity at the Groningen gas filed, The Netherlands, Geophys. Res. Lett. , 44, doi:10.10022017GL073878. https://doi.org/10.1002/2017GL073878 Eberhart-Phillips D., Han D.-H., Zoback M., 1989. Empirical relationships among seismic velocity, effective pressure, porosity, and clay content in sandstone, Geophysics , 54, 82– 89. https://doi.org/10.1190/1.1442580 Google Scholar CrossRef Search ADS   Egle D., Bray D.E., 1976. Measurement of acoustoelastic and third-order elastic constants for rail steel, J. acoust. Soc. Am. , 60, 741– 744. https://doi.org/10.1121/1.381146 Google Scholar CrossRef Search ADS   Fazio T., Aki L., Alba K., 1973. Solid earth tide and observed change in the in situ seismic velocity, J. geophys. Res. , 78, 1319– 1322. Google Scholar CrossRef Search ADS   Gurtin M., 1963. A generalization of the Beltrami stress functions in continuum mechanics, Arch. Ration. Mech. Anal. , 13, 321– 329. https://doi.org/10.1007/BF01262700 Google Scholar CrossRef Search ADS   Hatchell P., Bourne S., 2005. Rocks under strain: strain-induced time-lapse time shifts are observed for depleting reservoirs, Leading Edge , December, 1222– 1225. Henyey F.S., Pomphrey N., 1982. Self-consistent elastic moduli of a cracked solid, Geophys. Res. Lett. , 81, 903– 906. https://doi.org/10.1029/GL009i008p00903 Google Scholar CrossRef Search ADS   Hughes D., Kelly J.L., 1953. Second-order elastic deformation of solids, Phys. Rev. , 92, 1145– 1149. https://doi.org/10.1103/PhysRev.92.1145 Google Scholar CrossRef Search ADS   Johnson P., Rasolofosaon P., 1996. Nonlinear elasticity and stress-induced anisotropy in rock, J. geophys. Res. , 101, 3113– 3124. Google Scholar CrossRef Search ADS   Komatitsch D., Tromp J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation, Geophys. J. Int. , 139, 806– 822. https://doi.org/10.1046/j.1365-246x.1999.00967.x Komatitsch D., Vilotte J.-P., 1998. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. seism. Soc. Am. , 88, 368– 392. Kraaijpool D., Dost B., 2013. Implications of salt-related propagation and mode conversion effects on the analysis of induced seismicity, J. Seismol. , 17, 95– 107. https://doi.org/10.1007/s10950-012-9309-4 Google Scholar CrossRef Search ADS   Kristekova M., Kristek J., Moczo P., 2009. Time frequency misfit and goodness-of-fit criteria, Geophys. J. Int. , 178, 813– 825. https://doi.org/10.1111/j.1365-246X.2009.04177.x Mavko G., Mukerji T., Godfrey N., 1995. Predicting stress-induced velocity anisotropy in rocks, Geophysics , 60, 1081– 1087. https://doi.org/10.1190/1.1443836 Google Scholar CrossRef Search ADS   Murnaghan F., 1951. Finite Deformation of an Elastic Solid , Dover Publications Inc. Nur A., 1971. Effects of stress on velocity anisotropy in rocks with cracks, J. geophys. Res. , 76, 2022– 2024. Google Scholar CrossRef Search ADS   Nur A., Simmons G., 1969. Stress-induced velocity anisotropy in rock: an experimental study, J. geophys. Res. , 74, 6667– 6674. Google Scholar CrossRef Search ADS   O’Connell R., Budiansky B., 1974. Seismic velocities in dry and saturated cracked solids, J. geophys. Res. , 79, 5412– 5426. Google Scholar CrossRef Search ADS   Peter D. et al.  , 2011. Forward and adjoint simulations of seismic wave propagation on unstructured hexahedral meshes, Geophys. J. Int. , 186, 721– 739. https://doi.org/10.1111/j.1365-246X.2011.05044.x Poirier J., Tarantola A., 1998. A logarithmic equation of state, Phys. Earth planet. Inter. , 109, 1– 8. https://doi.org/10.1016/S0031-9201(98)00112-5 Google Scholar CrossRef Search ADS   Prioul R., Bakulin A., Bakulin V., 2004. Nonlinear rock physics model for estimation of 3D subsurface stress in anisotropic formations: theory and laboratory verification, Geophysics , 69, 415– 425. https://doi.org/10.1190/1.1707061 Google Scholar CrossRef Search ADS   Renaud G., Bas P.L., Johnson P.A., 2012. Revealing highly complex elastic nonlinerar (anelastic) behaviour of earth materials applying a newprobe: dynamic accoustic testing, J. geophys. Res. , 117, doi:10.10292011JB009127. Silver P., Daley T., Niu F., Majer E., 2007. Active source monitoring of cross-well seismic travel time for stress-induced changes, Bull. seism. Soc. Am. , 97, 281– 293. https://doi.org/10.1785/0120060120 Stacey F., 1992. Physics of the Earth , Brookfield Press. Tromp J., Tape C., Liu Q.Y., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, Geophys. J. Int. , 160, 195– 216. https://doi.org/10.1111/j.1365-246X.2004.02453.x Verdon J., Angus D., Kendall M., Hall S., 2008. The effect of microstructure and nonlinear stress on anisotropic seismic velocities, Geophysics , 78, D41– D51. https://doi.org/10.1190/1.2931680 Google Scholar CrossRef Search ADS   Zheng Z., 2000. Seismic anisotropy due to stress induced cracks, Int. J. Rock Mech. Min. Sci. , 37, 39– 49. https://doi.org/10.1016/S1365-1609(99)00090-8 Google Scholar CrossRef Search ADS   APPENDIX: WEAK FORM OF THE EQUATION OF MOTION UNDER INDUCED STRESS We have the equation of motion eq. (16), that is,   \begin{equation} \rho \,\partial _t^2{{\bf s}}={\boldsymbol \nabla} \cdot ({{\bf T}}^{\rm L1} -{{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0) , \end{equation} (A1)subject to the boundary condition (18), which is   \begin{equation} \hat{{{\bf n}}}\cdot [{{\bf T}}^{\rm L1}+{{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0]_-^+={{\bf 0}}. \end{equation} (A2)Upon dotting the equation of motion (A1) with a test vector w, followed by an integration by parts in which we use the boundary condition (A2), the weak form of the equation of motion is   \begin{eqnarray} \int _V\rho \,{{\bf w}}\cdot \partial _t^2{{\bf s}}\,{\mathrm{d}}V \ & = & \ -\int _V({\boldsymbol \nabla} {{\bf w}}):({{\bf T}}^{\rm L1}-{{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0):({\boldsymbol \nabla} {{\bf w}})\,{\mathrm{d}}\,V - \int _\Sigma \hat{{{\bf n}}}\cdot [({\boldsymbol \nabla} {{\bf s}})^T\cdot {{\bf T}}^0 -{{\bf T}}^0\,{\boldsymbol \nabla} \cdot {{\bf s}}-{{\bf s}}\cdot {\boldsymbol \nabla} {{\bf T}}^0]_-^+\cdot {{\bf w}}\,{\mathrm{d}}\,\Sigma . \end{eqnarray} (A3)Another lengthy, tedious calculation enables us to eliminate the surface integral and derivatives of the induced deviatoric stress, resulting in the weak form   \begin{eqnarray} \int _V\rho \,{{\bf w}}\cdot \partial _t^2{{\bf s}}\,{\mathrm{d}}V & = & -\int _V({\boldsymbol \nabla} {{\bf w}}):({{\bf T}}+\Delta {{\bf T}})\,{\mathrm{d}}\,V + \int _V[({\boldsymbol \tau} ^0:{\boldsymbol \nabla} {{\bf s}})\,{\boldsymbol \nabla} \cdot {{\bf w}}+{{\bf s}}\cdot {\boldsymbol \tau} ^0\cdot ({\boldsymbol \nabla} {\boldsymbol \nabla} \cdot {{\bf w}}) -({\boldsymbol \nabla} \cdot {{\bf s}})\,{\boldsymbol \tau} ^0:{\boldsymbol \nabla} {{\bf w}}-{\boldsymbol \tau} ^0:({{\bf s}}\cdot {\boldsymbol \nabla} {\boldsymbol \nabla} {{\bf w}}) ]\,{\mathrm{d}}\,V . \end{eqnarray} (A4)Unfortunately, the second derivative of the test function, $${\boldsymbol \nabla} {\boldsymbol \nabla} {{\bf w}}$$, is required. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial