Stress-dependent elasticity and wave propagation — New insights and connections

To establish a consistent framework for seismic wave propagation that accommodates the effects of stress changes, it is critical to take into account the different definitions of stress and their corresponding effects on seismic quantities (e.g., wave speeds) as dictated by continuum mechanics. Revisiting this fundamental theoretical foundation, we first emphasize the role of stress within various forms of the wave equation resulting from different choices of stress definitions. Subsequently, using this basis, we investigate connections among existing theories that describe the variation of elastic moduli as a function of changes in stress. We find that there is a direct connection between predicting stress-in-duced elastic changes with the well-known third-order elasticity tensor and the recently proposed adiabatic pressure derivatives of elastic moduli. However, each of these approaches has different merits and drawbacks in terms of experimental validation as well as in their use. In addition, we investigate the connection with another general approach that relies on micromechanical structures (e.g., cracks and pores). Although this can be done algebraically, it remains unclear as to which definition of stress and which corresponding constitutive relationship should be considered in practical scenarios. We support our analysis with validations using previously published benchmark experimental data.


INTRODUCTION
The knowledge of stress variations and their effects on seismic wave propagation is crucial for seismic processing in stressed for-mations, reservoir characterization, and 4D monitoring. Therefore, it is critical that the theoretical foundation underlying the calculation of the effects of stress on seismic waves is developed self-consistently, so that seismic-based technologies can be built upon it. To address stress-induced wavefield effects, there are two primary problems that need to be considered: 1) How do seismic waves propagate in media under stress? How do we take into account the effects of the state of stress in the wave equation? 2) How do elastic moduli vary with stress? Can we establish a parameter model for wave propagation experiments that is consistent at the field and laboratory scales?
To our knowledge, the treatment of the first problem appears to be sparse in the seismic exploration literature. This may be partly due to the success of the use of a familiarbut approximatewave equation that is valid when the initial stress of the medium is ignored. This point is noted in, for example, box 8.5 of Aki and Richards (2002). However, a complete treatment of this problem is given in the neighboring fields of global seismology (Dahlen and Tromp, 1998) and ultrasonics (Wallace, 1967;Thurston, 1974). For this reason, we first present a synopsis of how the effects of stress can be incorporated into the wave equation and the importance of various possible definitions of stress that can be used to accomplish this.
With regard to the second problem, there are two general approaches that exploration geophysicists adopt to quantify the effects of stress on elastic moduli, namely, the theory of third-order elasticity (TOE) and theories based on micromechanical structure analysis. The former builds on the original framework from the ultrasonic community for characterizing the linear effects of stress on wave speeds, which are in turn related to the crystal properties under stress. The literature on TOE commonly refers to the original formulation of Brugger (1964) and Thurston and Brugger (1964).
Some additional details can also be found in Thurston (1965bThurston ( , 1974. In seismic exploration, several variants of TOE theory have been proposed (Johnson and Rasolofosaon, 1996;Rasolofosaon, 1998;Winkler et al., 1998;Sinha and Plona, 2001;Sarkar et al., 2003;Prioul and Lebrat, 2004;. Each of these theories involves certain modifications or assumptions. In general, the TOE theory has the benefit of not relying on micromechanical models, which, in principle, makes it more general than other theories that make assumptions, e.g., the crack shape or grain contacts. Moreover, TOE represents a full sixth-ranked tensor with, at most, 56 independent parameters that can be used to compute any wave signatures in arbitrarily anisotropic media ). This number reduces to three when assuming the TOE tensor to be isotropic. Several studies report that these three additional parameters appear to be sufficient to characterize the stress-dependent effects in many common applications (Sinha and Plona, 2001;Sarkar et al., 2003;. Others advocate the use of theories based on analyses of micromechanical structures, such as cracks and pores, and their responses under stress. In this approach, these features are considered to be the most compliant parts of rock pore spaces and most sensitive to deformation upon a change in stress. Therefore, the variation of elastic moduli with stress can be interpreted in terms of parameters related to deformations of the corresponding pore spaces. Several variants of the micromechanical theory exist with different modifications and assumptions (Mavko et al., 1995;Sayers and Kachanov, 1995;Sayers, 2002Sayers, , 2010Tod, 2002;Shapiro, 2003Shapiro, , 2017Shapiro and Kaselow, 2005;Verdon et al., 2008;Angus et al., 2009;Gurevich et al., 2011;Collet et al., 2014), and some comparisons may be found in Price et al. (2017). Most variants make use of micromechanical structure parameters that are described by second-and fourth-rank tensors. Only in one variant (Shapiro, 2017) is a sixthorder tensor involved. In general, this approach allows us to potentially gain more intuitive insights into rock properties than TOE theory. In some variants of the theory, one can characterize a nonlinear stress dependence of wave speeds in a larger range of stress changes because some parameters related to micromechanical structures may be considered exponentially related to stress (Tod, 2002;Shapiro, 2003). The required number of parameters in this approach can be large because each micromechanical structure (e.g., crack) must be described individually. Several studies have reported that, in practice, the total contribution from multiple micromechanical structures can be described and estimated in an effective sense and some stringent requirements, such as crack shape, may be relaxed (Mavko et al., 1995;Sayers, 2002;Gurevich et al., 2011;Shapiro, 2017).
Recently, Tromp and Trampert (2018) propose a theory to characterize the stress dependence of elastic moduli in terms of their adiabatic pressure derivatives. Its preliminary results were validated by ab initio calculations (e.g., Karki et al., 2001;Militzer et al., 2011), in which predicted elastic moduli under changes of stress were compared with those computed directly from molecular dynamics (Tromp et al., 2019). Analogously to the TOE theory, the approach aims at describing the linear effects of stress on elastic moduli but stems from a different rootthe nonuniqueness of the stiffness tensor in continuum mechanics. No micromechanical model is assumed, and the linear effects of stress on elastic parameters are captured through more intuitive mathematical parameters the adiabatic pressure derivatives of the elastic moduli themselves. These derivatives can be written out as a fourth-order tensor in the same manner as for the usual stiffness tensor. Tromp and Trampert (2018) also show that the results from this theory can easily be used with an appropriate wave equation that honors the effects of stress, which, in principle, fully enables its application in seismic modeling and inversion.
To address the two general problems posed above, in this expository paper we first emphasize the importance of the initial stress, definitions of the stress tensor, and their effects on the wave equation following from continuum mechanics. Subsequently, we use this theoretical foundation as a common ground to look at previously proposed theories for predicting stress-dependent elastic moduli and discuss some of their relations. Note that we do not aim to present a new theory for characterizing stress-dependent elastic moduli in this paper. Rather, we present a concise overview of the topic and new insights on connections among existing theories suitable for an audience of exploration geophysicists in the hope that these connections will pave the way for future developments to fully address the aforementioned key problems.

WAVE EQUATIONS IN THE PRESENCE OF STRESS
From continuum mechanics it follows that the wave equation in a stressed medium can be written in three equivalent forms (Dahlen and Tromp, 1998), namely, where ρ 0 is the density, s is the displacement vector, and T 0 ¼ −p 0 I þ τ 0 denotes the state of stress at which waves must propagate, with pressure p 0 ¼ −ð1∕3ÞtrðT 0 Þ and deviatoric stress τ 0 ¼ p 0 I þ T 0 . The three types of stress perturbations used to describe wave propagation are the incremental Lagrangian description of the Cauchy stress (T L1 ), the incremental first (T PK1 ), and second (T SK1 ) Piola-Kirchhoff stresses. A synopsis on the derivation of equation 1 and the three definitions of stress are presented in Appendix A.
If the current state of the medium is stress-free (T 0 ¼ 0), any of the three choices lead to the same wave equation normally used in seismic processing. In the presence of a nonzero prestress, the three definitions differ and involve different (but related) constitutive relations. Here, we express the three possible constitutive relations as with three choices of stiffness tensors, Λ, Υ, and Ξ, which are related to each other as follows: The Christoffel equation corresponding to equation 1 can then be expressed as where a denotes an eigenvector, c denotes the corresponding phase speed, and B denotes the Christoffel matrix with elements Here,k ¼ k∕k denotes a unit wave vector that points in the phase direction. The phase speed c may be found by solving the eigenproblem 5.
To choose among the three stress options for the wave equation and associated constitutive relations, there are at least three points to consider: As a consequence, Λ and Υ only have the following symmetries: whereas Ξ has the "full" symmetry 3) It is either Λ or Υ but not Ξ that appears in the Christoffel equation and thus in expressions of phase speed determined in equation 5.
As a result, previous researchers have noted that stress-induced anisotropy will lead to stiffness tensors with limited symmetry because we have to work with Λ or Υ but not Ξ. We summarize the important aspects associated with three stress definitions in Table 1.
Considering the above three points and Table 1, it appears that working with T L1 is the most convenient choice because (1) the Lagrangian stress perturbation (T L1 ) still maintains its symmetry over a transposition and (2) its associated stiffness tensor elements directly relate to the Christoffel matrix and, in turn, the wave speeds. Moreover, under the static equilibrium condition (∇ · T 0 ¼ 0), the wave equation related to the Lagrangian stress perturbation (T L1 ) becomes which leads us to yet another advantage. According to equation 9, if there are no deviatoric components, τ 0 ¼ 0, the wave equation reduces to the usual one in seismic processing and is valid for hydrostatically stressed media (including the case of unstressed media), namely, This τ 0 ¼ 0 assumption also results in T 0 ¼ −p 0 I, which leads to Λ and Υ becoming and exhibiting the same symmetries as Ξ in equation 8. Consequently, if we choose to consider a hydrostatically stressed (or unstressed) medium and make use of wave equation 10, together with a Υ, whose elements appear naturally in the wave speed expressions, we will have a consistent framework for studies/experiments that properly honor the effects due to the state of stress according to continuum mechanics. At the same time, we will also have a framework that is closely similar to one that is conventionally used in the seismic exploration communitysimilar looking wave equation (equation 10) with stiffness tensor that includes the effects of p 0 terms (equation 12).
We observe that the presence of stress can directly affect the corresponding wave speeds and that the degree of its influence directly depends on its magnitude relative to the magnitudes of the elastic moduli (Ξ ijkl ). In the 40 MPa case, the magnitude of the prestress is much smaller than the magnitude of the elastic moduli, and its influence is less than 0.4%. The error grows by approximately 10 times in the case of 400 MPa given that everything else is kept constant.
In practice, we strive to estimate properties of the subsurface via wave experiments (laboratory or seismic). Without a proper and consistent treatment of stress, its associated constitutive relation, and the wave equation, the obtained parameters will certainly be flawed and subsequent uses of such results will further aggregate errors, especially when ones tries to improve subsurface models via iterative wave-equation-based inversions.

DEFINING DIFFERENT STATES OF STRESS
Before proceeding further, it is advantageous to go over the concept of states of stress. Each state describes a particular configuration of the medium, which provides a basis for us to relate changes from one to another and eventually establish an appropriate wave equation. In the context of the current problem, we identify three states of stress, as shown in Figure 3: 1) The reference configuration is conveniently taken where present data on elastic moduli are available (e.g., seismically inverted models or laboratory measurements at some pressure).
2) The intermediate configuration denoting the same medium but subjected to arbitrary choices of stress setups (e.g., with an additional triaxial stress change on the medium). Several authors also refer to this configuration as a medium with static bias Winkler et al., 1998;Sinha and Plona, 2001;Lei et al., 2012).  Figure 1. A comparison of (a) phase speed, (c) group speed, and (b and d) their differences in the presence of 40 MPa pressure and additional (approximately 15%) deviatoric stress. The blue color denotes P-waves, whereas the orange color denotes SV-waves. Note that the phase and group speeds are almost identical (the dashed and solid lines overlapping) despite the presence of stress due to its relatively small magnitude. 3) The current configuration, which describes the medium in the intermediate state 2 with additional deformations from wave propagation.
Going from reference state 1 to intermediate state 2, one can choose to use any of the three theories to relate the elastic moduli at some state of stress to another. However, we emphasize that there are several options for what one can and should define as a reference. To study wave propagation in the current state 3, we need both knowledge of the medium in 2 and a wave equation that properly accounts for the effects of stress as defined by continuum mechanics. In the previous section, it is shown that the wave equation that defines how waves propagate in a medium under some arbitrary stress field T 0 is given in equation 1 or 9. Given the appropriate linearized constitutive relations (equation 3 or 4) and the knowledge of T 0 , we can propagate wavefields that correctly honor the effects of the state of stress. In other words, this wave equation describes how to relate the intermediate state 2 to the current state 3.
With these considerations in mind, the goal then is to construct a rock-physics model that can connect changes in elastic moduli (Λ, Υ, or Ξ) between the reference 1 and intermediate 2 states so that proper wave studies can be conducted. We emphasize that there are at least two more points that we have to keep in mind: 1) In seismology, we always invert for medium parameters in the presence of some in situ stress in the subsurface. This is also the case in laboratory measurements of elastic moduli, which are normally done under some state of nonzero stress. The fact that these parameter estimation processes commonly rely on the wave equation for hydrostatically stressed media equivalent to equation 10 suggests that the recovered elastic moduli from these approaches represent an approximation of Υ related to some in situ pressure p 0 (equation 12), which in turn relates to Ξ at the same in situ pressure. Therefore, in view of these previous developments, it is convenient to define the reference state to be associated with hydrostatically stressed media (Figure 3), where the elastic moduli can be obtained from the current model building workflow or laboratory measurements at some effective hydrostatic stress. Strictly speaking, in view of equation 10, these moduli provided by the conventional model building workflow are only approximate because the effects of subsurface in situ deviatoric stresses are neglected due to the use of an approximate wave equation. This is generally valid in seismology because the effects of deviatoric stress are assumed to be much smaller than those due to hydrostatic pressure (Dahlen and Tromp, 1998). However, this problem may not be pertinent in laboratory settings because the choice of applied stress in the reference state can be controlled. 2) Another important point must also be made regarding the difference between isothermal and adiabatic quantities. In seismology, we generally assume that adiabatic elastic moduli at constant entropy S (Λ, Υ, or Ξ) are well-approximated by the isothermal moduli at a constant temperature T (Λ T , Υ T , or Ξ T ) and vice versa (e.g., Johnson and Rasolofosaon, 1996;Dahlen and Tromp, 1998). Therefore, given an appropriate frequency range, the model parameters obtained from seismic inversion, which are adiabatic moduli, can directly be related to those isothermal moduli obtained in laboratory settings at a constant temperature. Even though we use a similar assumption here, distinguishing them becomes important when we consider some variants of the TOE theory and try to theoretically relate them to the theory with adiabatic pressure derivatives in the "Third-order elasticity" section.

THEORIES FOR STRESS-DEPENDENT ELASTICITY
Predicting changes in elastic moduli due to changes in stress is a long-standing problem, and many previous studies have proposed theories to characterize this behavior (Murnaghan, 1951;Toupin and Bernstein, 1961;Brugger, 1964;Ghate, 1964;Nur and Simmons, 1969;Wallace, 1970;Nur, 1971;Thurston, 1974). In this section, we go over three general approachesnamely, adiabatic pressure derivatives, TOE, and micromechnical structuresused to characterize changes in elastic moduli due to changes in stress. The first two represent a linear approximation (Figure 4), whereas the third one may be nonlinear. We use continuum mechanics as the common ground of each of the three theories and discuss some of their direct implications. This enables us to provide some connections between the latest theory with adiabatic pressure derivatives and the TOE theory because they stem from the same continuum mechanics basis.
In this section, we use a general expression of the form to describe the variation of elastic moduli with stress. Here, Γ ijkl denotes the elastic moduli in the unstressed reference state 1, f ijkl ðT 0 mn Þ denotes a linear or nonlinear function describing their

Intermediate state
Arbitrarily pre-stressed

Intermediate state
Arbitrarily pre-stressed

Intermediate state
Arbitrarily prestressed Figure 3. Configuration states for studying stress-dependence elasticity with seismic waves. The reference state is defined where the elastic moduli can be measured.

Option 1: Several variants of the third-order elasticity (TOE) theory
Option 2: Theory based on pressure derivatives of elastic moduli Figure 4. Schematic showing the notion of a local linear approximation adopted by the TOE theory and the theory based on adiabatic pressure derivatives (Tromp and Trampert, 2018). Other theories that are based on micromechanical structrues aim to describe this nonlinear variation in a larger range of stress change in terms of parameters related to deformations of micromechanical structures.
change with stress T 0 , and Ξ is the corresponding modulus appropriate for T SK1 , which possesses full symmetry relations 8 and lives in the intermediate state 2. The corresponding expressions for Λ and Υ may be obtained from equations 3 and 4. In principle, additional consideration on the deformation gradient and density variation between the two states must also be honored. Different theories take this point into account variously, and it shall become clear once we specify f ijkl ðT 0 mn Þ in equation 18 depending on the theory.

Adiabatic pressure derivatives
The theory of Tromp and Trampert (2018) aims to capture the linear effects of stress on elastic moduli with the use of adiabatic pressure derivatives of the elastic moduli themselves. They approach this problem on the basis of the nonuniqueness of the stiffness tensor, which was previously used by Dahlen (1972) with the same goal but in a different formulation. In the context of this theory, we can express f ijkl in equation 18 as where we enforce that Ψ ijklmn T 0 mn also has the required symmetry among indices i, j, k, and l (equation 8). Equations 18 and 19 then represent the most general form of elastic tensor that can be used to capture the linear effects of T 0 . In other words, they represent a linear approximation of the elastic tensor with respect to its value at T 0 ¼ 0 given by Γ ijkl . By making different choices of Ψ ijklmn T 0 mn , variants of this theory can be derived.
1) Let us first consider the following choice proposed by Dahlen (1972): where the summation in parentheses ensures the necessary symmetry among indices i, j, k, and l. The scalar parameters a and b can be any arbitrarily chosen quantities, which may vary spatially. Dahlen (1972) shows that choosing a ¼ −b ¼ 1∕2 leads to expressions of T L1 and P-wave phase speeds that are independent of pressure p 0 , which may be desirable in some cases. 2) Consider an isotropic reference Γ given by where κ and μ denote the spatially variable isotropic bulk and shear moduli at p 0 ¼ 0, respectively (or at any other p 0 with Ξ ijkl at that pressure). Building on the choice of an isotropic background and equation 20, Tromp and Trampert (2018) propose an alternative choice for a and b, namely, where κ 0 and μ 0 denote the spatially variable adiabatic pressure derivatives of the bulk and shear moduli, respectively. The derivatives with respect to the deviatoric stress τ 0 are assumed to be negligible in magnitude with respect to κ 0 and μ 0 .

Substituting equation 22 into equation 20 gives
which leads to the Christoffel matrix We observe that the choice of using these adiabatic pressure derivatives (κ 0 and μ 0 ) conforms to the general objective of the theory that aims to capture the first-order (linear) effects of T 0 and quantify them in an intuitive way. Furthermore, it allows us to conveniently isolate pressure-dependent terms in the form of linear approximations κ þ κ 0 p 0 and μ þ μ 0 p 0 from the those that depend on τ 0 . 3) Tromp et al. (2019) generalize the theory for an arbitrary anisotropic background, still in the scope of the linear regime, as follows: The first two terms are similar to those in equation 20 with The last term is dependent on the adiabatic pressure derivative of the background elastic tensor Γ 0 ijkl . If Γ ijkl is assumed to be isotropic, with its pressure derivatives characterized by κ 0 and μ 0 , equation 25 reduces to equation 23 originally proposed by Tromp and Trampert (2018).
Equation 25 maintains the same advantage of conveniently isolating pressure-dependent terms, which becomes apparent when we write out the corresponding Ξ, Λ, and Υ in terms of p 0 and τ 0 as follows: The first term in equations 26-28 is Γ ijkl þ Γ 0 ijkl p 0 , which represents the desired linear approximation in terms of the adiabatic pressure derivative on the background elastic tensor Γ ijkl . The remaining dependence on p 0 in Λ exists in combinations of delta functions but disappears upon contraction withk ikk when constructing the Christoffel equation (equation 5). Substituting equation 27 or 28 into equation 5, we obtain the Christoffel equation: which possesses a noted advantage of having only one pressuredependent term, namely, Γ ijkl þ Γ 0 ijkl p 0 . It is apparent from the generalized theory, captured in equation 25, that when working with adiabatic pressure derivatives, the additional contribution f ijkl ðT 0 mn Þ can be directly expressed in terms of stress and we do not need to involve strain, which is the case in the TOE theory. Moreover, without a loss of generality, because we are looking at a linear function in p 0 and τ 0 , we can also consider Γ as Ξ measured at some other reference state, such as at some convenient hydrostatic stress (equation 12), as previously noted. Equations 26-28 then predict the linear changes in the elastic tensor from that reference state due to additional induced stress T 0 from some external process. In other words, suppose the reference state is chosen at p 1 and we would like to predict elastic moduli at a new p 2 and τ 2 . We can achieve this by substituting into equations 26-28. Similar substitutions can be made in the Christoffel equation 29 for calculations of the wave speeds.

Third-order elasticity
In TOE theory, one makes use of a sixth-ranked tensor and its contraction with a variation in strain to predict changes in the elastic moduli with stress. In general, the TOE tensor can be defined in two ways: • Adiabatic TOE tensor: where U L denotes the Lagrangian internal energy density per unit mass expressed in terms of the local instantaneous Lagrangian strain E L and the local entropy S L . Recall from our synopsis in Appendix A (equation A-41) that the definition of the stiffness tensor Ξ (or Γ in unstressed media) is the second derivative of U L with respect to E L . Therefore, the adiabatic TOE tensor is the derivative of such stiffness tensor while keeping the entropy S L constant.
• Isothermal TOE tensor: where we note that the isothermal TOE is the derivative of Γ while keeping the temperature T constant.
Note that in both definitions, we are considering a deformation from the unstressed reference state 1 to an intermediate state 2. As a result, we consider Γ as the counterpart of Ξ in an unstressed reference, and it is an adiabatic quantity by definition.
The first is the third-order derivative of the Lagrangian internal energy density, which is a purely adiabatic quantity, whereas the latter is the isothermal derivative of the adiabatic Ξ, which is a "mixed" quantity as denoted by the superscript. Both of them are defined with respect to a reference medium that is chosen to be unstressed in the context of ultrasonic studies. The method to solve for the TOE tensor in Thurston and Brugger (1964) was proposed based on equation 31, and its relation to equation 30 was provided by Brugger (1964). Skove and Powell (1967) note that the TOE tensor in equation 31 does not have the same symmetry as those expressed in equation 30 because the pair of indices from isothermal differentiation cannot necessarily be interchanged with the other two pairs. Thurston (1974) points out that the assumption of full symmetry of the TOE tensor in equation 31 leads only to small errors that can be neglected in the presence of other experimental errors. Symmetry properties of the TOE tensor in equation 30 are studied by . Distinguishing between both definitions of TOE is important for linking it with the theory with adiabatic pressure derivatives from the previous section.
To our knowledge, most variants of the TOE theory in seismic exploration appear to adopt the TOE tensor in equation 30. The only exceptions are the studies of Johnson and Rasolofosaon (1996) and Rasolofosaon (1998), which made use of the original workflow of Thurston and Brugger (1964), leading to equation 31. Therefore, in view of the general formula in equation 18, we can consider the TOE tensor in equation 30 and write where ϵ mn ¼ Γ −1 mnpq T 0 pq denotes the corresponding linearized Lagrangian strain that governs the change in the elastic tensor from the reference state upon contraction with the sixth-ranked TOE tensor in equation 30. Strictly speaking, ϵ mn must also include a nonlinear term that involves the inverse of the TOE tensor, but several authors have argued that the magnitude of such a term is small and can be neglected in practice (Sarkar et al., 2003;. Equation 32 serves as the basis for three notable variants of the TOE theory, as listed below: 1)  consider equation 32 together with two additional terms: one stress term that stems from the difference between Ξ and Λ (equation 3) and the other related to a static deformation from the unstressed reference to the intermediate state, expressed in terms of Γ ijkl . Both terms are listed in their equation A-5. An extension of the theory to the elasticplastic regime was proposed by Sinha and Plona (2001). 2)  make use of the result by , but they argue that both additional terms are relatively small in practice and can effectively be neglected in an estimation of the TOE tensor. This results in their equation 1, which is an approximation to Λ but has the full symmetry expressed in equation 8. 3) Sarkar et al. (2003) reach the same final equation 32, but their derivation stems from a different approach. They assume that the range of strain (and stress) change is small, such that the deformation gradient relating the unstressed reference and the intermediate states is well-approximated by a delta function. The difference between Ξ and Λ (equation 3) is kept in this scheme.
Analogous to the theory with pressure derivatives, the TOE theory is linear in stress changes. Therefore, we can adopt a similar strategy to change the unstressed reference medium to a hydrostatically stressed medium. Subsequently, estimation of the TOE tensor in equation 30 is accomplished by a regression analysis with appropriate restrictions on the range of stress variation to ensure consistency with the underlying linear approximation.

Connections to the adiabatic pressure-derivative theory
At the most fundamental level, the TOE theory and the theory based on pressure derivatives (equation 25) are after the same goal of predicting linear effects of changes in stress on seismic wave propagation. Even though the latter was not established until recently, the adiabatic pressure derivative as a quantity is by no means foreign. In principle, these derivatives can be used to calculate elastic moduli at elevated pressure without relying on the TOE tensor (Thurston, 1965b(Thurston, , 1974. The theory of Tromp and Trampert (2018) is a significant step forward in this regard because we can now use the same information on the pressure derivatives to predict elastic moduli under an arbitrary state of stress nearbya task that was previously only possible by using the TOE tensor or an entirely different formulation based on microstructural parameters.
Following Thurston (1974), we show that there exists a connection between the TOE tensor (equation 30) and the adiabatic pressure derivatives (Γ 0 ijkl in the unstressed reference or Ξ 0 ijkl for the hydrostatically stressed reference). The starting point of our derivation is the assumption that the reference medium under consideration is in equilibrium under hydrostatic pressure, which we first assume to be unstressed. As a result, the Lagrangian strain (E L ) that corresponds to the change from the reference at zero pressure to the intermediate state at some elevated pressure can be described as a function of two variables: pressure change p and entropy S L . Consequently, in the unstressed reference, we can write Γ ijkl ðp; S L Þ and The latter derivative of strain with respect to pressure may be found from Finally, we may write which provides a direct bridge between the adiabatic pressure derivatives and the adiabatic TOE tensor (equation 30). We use s mnpq ¼ ðΓ −1 Þ mnpq to denote elements of the corresponding compliance tensor. It is apparent from equation 35 that given the stiffness/compliance tensor and the TOE tensor of any known medium, the adiabatic pressure derivatives essential for the theory of Tromp and Trampert (2018) can be computed. Similarly, equation 35 allows us to compute the TOE tensor needed by various formulations of TOE theory from the adiabatic pressure derivatives and the stiffness/compliance tensor. We further investigate this connection in more detail in the "Numerical comparison between TOE and pressure-derivative theories" section. An analogous bridging expression for isothermal pressure derivatives has a similar form as equation 35, but it involves isothermally measured quantities. A different expression that involves the mixed TOE tensor (equation 31) can also be established. Some details of the connections among these expressions may be found in Barsch (1967) and sections 11.11 and 29.2 of Thurston (1974). Without a loss of generality, it is also possible to define a reference medium to be hydrostatically stressed. The pressure change is then expressed W54 Sripanich et al.
Downloaded 08/08/21 to 83.82.230.117. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms with respect to this initial value, but the general form of equation 35 remains the same with Ξ at the reference pressure instead of Γ. The TOE tensor is then defined at this reference pressure instead of at zero pressure, as is the case for Γ.
It should be emphasized that in the derivation of equation 35, elastic tensors are assumed to be dependent upon only pressure change p and entropy S Lthe dependence on deviatoric stress is neglected. The same assumption is central to the development of equation 25 by Tromp et al. (2019). This has led to a decrease in the number of dependent parameters from a total of 56 in the TOE tensor to a maximum of 21 in the case of Tromp et al. (2019), i.e., one derivative for each of the stiffness coefficients. A discussion on the potential importance of the absent derivatives with respect to the deviatoric stress and rotational properties of the deformation gradient can be found in Maitra and Al-Attar (2021).
Based on this connection 35 and satisfactory preliminary results from the pressure-derivative theory (Tromp et al., 2019), using equation 25 to accomplish the task of predicting linear variations in elastic moduli with stress instead of TOE theory may represent a reasonable alternative in practice.

Microscale inclusions
Another approach to characterize the change in elastic moduli with stress is based on additional compliances related to deformations of microstructuresmicroscale inclusionswhich are often thought to be the most compliant parts of rock pore spaces and are most sensitive to deformations upon a change in stress. (Though often used in the context of "soft" inclusions, these approaches can also be used for inclusions that are stiffer than the background medium [Kachanov and Sevostianov, 2005]). Therefore, the common ground of this general approach is that one should think of the change in elastic moduli in terms of additional compliances, not stiffnesses. We refer the reader to the cited works for references and specific details on the developments based on this approach: Kachanov (1980); Sayers and Kachanov (1995); Schoenberg and Sayers (1995); Sayers (1999Sayers ( , 2002Sayers ( , 2009Sayers ( , 2010; Tod (2002); Shapiro (2003Shapiro ( , 2017; Kachanov and Sevostianov (2005); Shapiro and Kaselow (2005); ; ; Hall et al. (2008);Verdon et al. (2008); Angus et al. (2009);and Kachanov et al. (2010).
According to the derivation in Sayers and Kachanov (1995), which serves as the basis for most variants under this general approach, the linearized strain is used as opposed to the exact Lagrangian strain needed to derive equations 2-4. Therefore, in view of the above theoretical foundation, it is unclear to which of the elastic moduli (Λ, Υ, or Ξ) that this excess compliance should be added. This is notably different from the connection 35 that strictly follows from continuum mechanics. More importantly, the inability to connect the inclusion-based approaches above to one of the three elastic moduli has an important consequence: It means that a strict theoretical consistency between the Lagrangian-and inclusion-based approaches cannot be established at this time.
Because of this fact, we point out some connections that follow from algebraically equating the estimated Ξ from the theory captured in equation 25 to that from the microstructural approach under the assumption that the excess compliance should be added to Ξ with full symmetry (equation 8). As a result, we have where the additional compliance can be written as with α ij and β ijkl denoting the crack-related tensors. For the case of oblate-spheroidal cracks with a particular orientation n with normal B N and shear B T compliances, we may write where P r A ðrÞ ∕V denotes the crack-specific area (Gurevich et al., 2011).
Another notable alternative that relies on a different micromechanical structure (stiff and compliant pore spaces) is the porosity deformation approach (PDA) (Shapiro and Kaselow, 2005). Shapiro (2017) shows that analytical relationships between the TOE tensor and parameters related to such microstructure can be derived by equating the predicted changes in the elastic moduli from the PDA and the TOE theory. This may be accomplished by using the following relationship (Thurston, 1974): s pqrsuv ¼ −s uvmn s pqij s rskl c ijklmn : In the context of the PDA approach, the term s pqrsuv can be related to parameters controlling the properties of stiff and compliant pore spaces, and the total additional compliance can then be found from g ijkl ðT 0 mn Þ ¼ s pqrsuv T 0 uv . As we have discussed, a strict theoretical link between inclusionand Lagrangian-based approaches to predicting compliance changes is not yet available. As such, in view of continuum mechanics, the relationships between the additional compliance and the microstructure parameters for both approaches in this section are purely algebraic and are established by expressing the additional compliance in terms of the microstructure parameters of choice. In other words, the identities in this section would hold only if the different theories were indeed consistent with each otherwhich would have to be assumed because it cannot be theoretically proven. This is not the case for the identity in equation 35, which follows directly from the same theoretical foundation.

A TEST CASE
Equipped with our findings in the previous sections, we investigate the properties of the newly proposed theory based on the adiabatic pressure derivative in equation 19. We consider the specific case of an orthorhombic (ORT) Ξ subjected to triaxial stress changes that are aligned with its principal directions (T 0 contains no off-diagonal components), the change in stiffness moduli δΞ ijkl in Voigt notation can then be expressed according to equation 26 as Stress-dependent elasticity W55 Downloaded 08/08/21 to 83.82.230.117. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/page/policies/terms One can observe that the change of each elastic modulus is characterized by a combination of its derivative with respect to pressure and the magnitude of the pertinent stress components. For other more general cases, the relationships for δΞ ijkl can be straightforwardly obtained from equation 26 and expressed in Voigt notation if preferred.

Connection with TOE
To relate between the pressure derivatives of elastic moduli and the elements of the TOE tensor, we can expand equation 35 to obtain the desired derivatives. With the assumption of an isotropic TOE tensor with three independent parameters (c 111 , c 112 , and c 123 ) , the pressure derivatives for the ORT case in equation 40 in Voigt notation can be expressed as

Connection with microstructural theory
Unlike the connection with the TOE tensor in equation 35, the connection between the pressure derivatives used by Tromp and Trampert (2018) and the microstructure parameters can only be algebraic as discussed previously. To achieve this, we can equate the change in stiffness moduli as predicted by Tromp and Trampert (2018) in equation 25 to that from the microstructual theory (inverse of equation 37). Choosing n ¼ ½1; 0; 0, together with equation 38 and the previous ORT case, we can then express the pressure derivatives in terms of the crack parameters as follows: where p 0 and τ 0 ij are the pressure and deviatoric stress changes to the background as a result of adding microstructural inclusions, respectively (i.e., cracks, in this case). Because equation 36 represents the overall effective elasticity as a result of microinclusions, the resulting p 0 and τ 0 ij in equation 42 represent volume-averaged stress changes.

NUMERICAL COMPARISON BETWEEN TOE AND PRESSURE-DERIVATIVE THEORIES
Apart from the theoretical connection between the adiabatic pressure derivatives of elastic moduli and the TOE tensor as discussed in previous sections, we further illustrate in this section the similarities/differences between the two theories (that of Tromp and Trampert [2018] and the TOE theory) that rely on these two quantities. For this comparison, we rely on two VTI examples described by : a hydrostatic experiment on Jurassic North Sea shale and a hydrostatic/biaxial experiment on Colton sandstone.

Jurassic North Sea shale (hydrostatic)
For the shale data set, we use equation 41 and a similar reference background as  to convert the reported c 111 , c 112 , and c 123 to the following pressure derivatives that relate the elastic moduli (GPa) to the change in stress (MPa): We then use these results and equation 40 to predict the change in elastic moduli with stress according to equation 25. Figure 5 shows the results of this experiment with the solid lines indicating the predicted values of Υ. Because the ultrasonic measurements used in this example are dynamic moduli (obtained from measured wave speeds) and are related to either Λ or Υ, we must use equations 3 and 4 to relate them to Ξ. The dashed lines, which are almost completely overlain by the solid lines in the same figure denote the predicted values from the TOE theory. Even though no conversion is applied to the TOE results (estimated Ξ) to properly obtain Υ, the predictions fit the measurements well because the pressure range is 0-100 MPa, which is small in comparison with the magnitude (GPa) of the moduli Υ themselves.

Colton sandstone (hydrostatic and biaxial)
For the second example, we consider a data set from hydrostatic and biaxial experiments on Colton sandstone (Figure 6a). We follow the same process as before and convert the reported third-order elastic constants between the first and last states to Γ  Figure 6 with estimates from the TOE theory in dashed lines similar to those in  and those from the theory with pressure derivatives (equation 25) in solid lines.
Several observations can be made after examining the results from both experiments in this study: 1) Predictions from the pressure-derivative and TOE theories generally agree with each other, although a slightly better performance from the latter is observed in the states in which the applied stress is further away from being hydrostatic. Table 2 shows a comparison between the mean difference of the predicted values from either the TOE or the pressure-derivative theory and the values from laboratory measurements. Note that the TOE appears to perform better the further away the state of stress is from hydrostatic, but it requires third-order constants that are difficult to measure. 2) According to the pressure-derivative theory, the predicted values of δΞ ijkl and its corresponding Υ will follow the associated stress path. For example, we can observe in Figure 6b that the predicted Υ 33 follows the same path as the T 0 33 shown in Figure 6a.
3) To predict changes in the elastic moduli with stress, the pressure-derivative theory only needs the derivative of the wanted element, whereas the TOE theory requires all elastic moduli of the reference medium for strain computation.

DISCUSSION
We have shown in this study that there exists a consistent theoretical foundation that takes into account the effects of different stress definitions and states of stress, which directly influences the choice of the wave equation, constitutive relation, and speed of wave propagation. This consistent foundation follows directly from continuum mechanics and represents the answer to the first primary question that we raised in the "Introduction" section. In particular, we highlight the importance of the additional wave-equation stress terms connected to the absolute magnitudes of principal stressesthis implies that the additional terms may be neglected at relatively small stress variations (e.g., reservoir-scale compaction) and becomes important for larger magnitude stress perturbations (e.g., after larger earthquakes or at continental scales). Note that this is independent of the various choices one can make to describe changes in elastic moduli with respect to stress changes, e.g., TOE, microstructure-based theory, or other. Nonetheless, two points are worth keeping in mind: 1) Although the state of stress may be small in the context of exploration seismology (in the range of tens or hundreds of MPa), its effects on the wavefield can be quantified via proper choices of the wave equation and constitutive relation (equations 1 and 6). Errors from any form of approximation (e.g., ignoring the effects of the state of stress) can then be analyzed. In light of the growing use of waveform inversion for subsurface model estimation, choosing an appropriate wave equation becomes even more important.
2) Concerning the speed of wave propagation, equation 6 suggests that it is either Λ or Υ but not Ξ that is involved. Therefore, one must be mindful of what stiffness coefficients (Λ, Υ, or Ξ) are obtained from wave experiments at the field and laboratory (e.g., dynamic moduli measurement) scales. This has to be carried out regardless of the theory that one wishes to use (TOE, microstructure-based, or pressure derivative-based). Without a proper treatment, the estimated stiffness coefficients (or microstructure parameters) will surely be biased.
In view of the second primary problem identified in the "Introduction" section, we have discussed three different theories that can be used to describe how elastic moduli vary with stress changes: TOE, inclusion-based, and pressure-derivative theories. Each one of them stems from a different basis and has its own advantages and disadvantages. Building on the foundation of continuum mechanics, we have shown that there is a theoretical connection (equation 35) between the TOE tensor and the pressure derivatives; however, such a connection relating the TOE tensor and the pressure derivatives to inclusion-based parameters is not yet availableonly an algebraic connection is possible. Furthermore, the particular connection between the TOE tensor and the pressure derivatives also allows us to unravel some similarities and differences between the characteristics of the TOE theory and the pressure-derivative theory when it comes to their use in practice. However, we emphasize that further study is needed to decisively argue for any single theory in any particular situations and/or applications.
Although it is not the objective of this study, one important application of the models herein is the inverse problem of inferring the medium's stress state from observed wave data. The more likely scenario being a time-lapse change in the medium, e.g., due to reservoir compaction, preseismic stress build-up, or postseismic stressrelease accommodationin which case, the observations would be time-lapse waveforms from either active or passive data. To solve the inverse problem for stress changes, in addition to the wave data one would also need a priori knowledge of the baseline elasticity and the stress sensitivity of the medium. The stress sensitivity comes in the form of the sixth-rank tensor for TOE, or the fourthrank tensor for the pressure-derivative approach. Presumably, these would come from supporting experimental data. Here, we note that the pressure-derivative approach may have the practical advantage of not only requiring the observation of fewer parameters but also simpler experiments because only (isotropic) pressure derivatives are needed. A more thorough study on this inverse problem is the topic of further research.
With regard to the inverse problem of subsurface model building, we note that some previous studies such as  and Keys et al. (2017) have investigated the theoretical bounds of various elements of the TOE tensor in the hope of constraining the inversion results and helping with the convergence during the process. However, we emphasize that regardless of the value of various elements of the TOE tensor and their associated bounds, our connection in equation 35 between the TOE tensor and the pressure derivatives remains valid. Although our work here is concerned with providing physical connections and insights into the effects of stress on wave propagation, further understanding into the role and behavior of the different formulation choicesin the context of different stress regimes (e.g., reservoir or lithospheric scales)is the subject of future study. Note that the TOE appears to perform better the further away the state of stress is from hydrostatic, but it requires third-order constants that are difficult to measure.
In addition, it is important to point out that the elasticity models that we discuss in this work have limitations tied to the intrinsic assumptions of theories behind themwe see our work as a starting point for a comprehensive analysis of the effects of stress changes on wave properties. First, here we envisage scenarios in which the dominant wavelengths are large compared with heterogeneity scales, the wave-strain magnitudes are small, and the medium is in static stress equilibrium during the wave propagation experimentas such, nonlinear wave-medium interactions (e.g., clap cracks) or strong resonance effects are assumed to be not present/observed. Second, we here treat the purely elastic case, thus neglecting effects such as wave attenuation due to intrinsic and/or effective absorption or fluid-related effects (e.g., poroelasticity). These effects are of course important in their own contexts, and their relation to stress-induced wave elasticity deserves investigation in future studies. Third, we ignore effects due to thermal expansion and strictly consider the effects of stress changes on moduli. In principle, when there are fluid-content changes in a reservoir, the temperature and the state of stress change. Further investigations are required to incorporate the stress and temperature effects (Wang and Nur, 1988;Eastwood, 1993) on elastic moduli.
Finally, we note that here, as in general exploration seismology, rocks are generally assumed to be hyperelastic and hysteresis (stress-strain curves being different for loading and unloading) is often ignored (Mavko et al., 2009)implying that plastic deformation effects are not taken into account. Together with the fact that the Lagrangian and Eulerian descriptions of motions are rarely differentiated (in contrast to what we present here), this has led to challenges in describing the shapes of stress-strain curves over large ranges of stress and strain. However, we note that research in the ultrasonic community has led to, for example, the usage of bimodulus (bilinear) materials, in which material stiffnesses can be modeled in accordance with, e.g., crack characteristicsopen or closed, and other modifications to describe hysteresis (Dyskin et al., 2012;Broda et al., 2014). Therefore, in principle, it is possible to bring in these developments and further improve our description of stress-dependent elasticity for exploration seismology, e.g., by accommodating plastic deformation mechanisms in the future.

CONCLUSION
As dictated by continuum mechanics, we emphasize that in the presence of nonzero stress it is important to distinguish among different stress definitions, in which each choice pertains to a different wave equation and a different set of elastic moduli. Their differences become more prominent with the larger magnitudes of stress changes with respect to the reference state. Building on this foundation, we show that there is an explicit connection between the TOE tensor and adiabatic pressure derivatives of elastic moduli, which allows the recently developed theory with pressure derivatives to benefit from existing research on TOE. Finally, we point out that with regard to the definition of strain there are fundamental differences in the underlying framework and assumptions between the pressure-derivative and TOE theories on the one hand and those based on micromechanical structures on the other hand. Therefore, further work is needed to connect microstructure-based approaches with those stemming from continuum mechanics in a consistent and verifiable manner.

ACKNOWLEDGMENTS
We thank L. Adam for helpful discussions and R. Prioul for sharing the experimental data. We are grateful to J. Shragge, M. Sen, D. Foster, and the anonymous reviewers for their constructive comments during the reviewing process. This work is supported by the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP/2007(FP/ -2013 grant agreement number 320639 (iGEO).

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

DEFINITIONS OF STRESS AND CORRESPONDING WAVE EQUATION
The background for this review is the treatise of Dahlen and Tromp (1998), from which we adopt the same notation convention for convenience. The bold font is used to denote tensors (including vectors and matrices), but we will switch to index notation when indices become important, especially when working with fourthranked tensors. We note that indices involving derivatives are arranged as ð∇ x rÞ ij ¼ ∇ i r j ; consequently, ð∇ x rÞ T ij ¼ ∇ j r i . The divergence of a tensor is obtained by contracting on its first index; i.e., ð∇ · TÞ j ¼ ∇ i T ij .

Lagrangian versus Eulerian descriptions
In seismology, we treat earth as a continuuma continuous distribution of matter that can interact with any subjected forces. The associated motion of a continuum can be described in two ways: Lagrangian or Eulerian. The former involves labeling particles in the material by their position and following them through time to provide a kinematic description of the motion. We shall refer to the particle, whose position at time t ¼ 0 is x, as particle x and denote its position at t ≥ 0 by rðx; tÞ, i.e., rðx; 0Þ ¼ x. The corresponding velocity of particle x can be written as u L ðx; tÞ ¼ ∂ t r, where we use a superscript L to denote Lagrangian quantities.
On the other hand, the Eulerian description does not focus on the motion of individual particles through time, but instead it focuses on fixed points in the considered coordinate space. We label those fixed points by position r and the velocity of the particle that is at point r can be written as u E ðr; tÞ with the superscript E denoting Eulerian quantities. The knowledge of u E ðr; tÞ at all points r at all time t ≥ 0 is an alternative yet equally valid description of the motion of a continuum.
For every scalar, vector, or tensor quantity q in a continuum, there exists a Lagrangain q L ðx; tÞ and an Eulerian q E ðr; tÞ description. These descriptions can be related by q L ðx; tÞ ¼ q E ðrðx; tÞ; tÞ; (A-1) which represents the quantity q for particle x at position r at time t.
Differentiating equation A-1 with respect to time t leads to Equation A-2 describes two types of changes in variable q experienced by an observer riding on particle x. The term ∂ t q E represents a change experienced by a stationary observer, whereas u E · ∇ r q E represents an additional change due to the motion of the particle itself. The combination of both changes, D t ¼ ∂ t þ u E · ∇ r , is called the material derivative, which only acts on Eulerian quantities while holding x fixed. Receivers in seismic experiments directly record the motion rðx; tÞ of particle x to which they are attached. Hence, adopting the Lagrangian viewpoint is a natural and appropriate choice in seismological studies. It will become apparent in the subsequent sections why distinguishing between the Lagrangian and Eulerian viewpoints is crucial in the studies of stress-dependent effects on the material and the development of an appropriate wave equation in stressed media. This is primarily due to the fact that exploration geophysicists generally use the Eulerian Cauchy stress tensor to approximate its Lagrangian counterpart under the assumption of linear elasticity.

Lagrangian strain
We can consider relative motion at two adjacent particles (Lagrangian) or two adjacent points in the coordinate space (Eulerian) to characterize the deformation. Here, we focus on the Lagrangian case and consider two particles that were initially at x and x þ dx but are now located at r and r þ dr, where and T denotes the transpose. The term F is referred to as the deformation tensor, which measures the total deformation experienced by a moving material particle and its vicinity. Its determinant (det F ¼ J) governs the volumetric change between the volume at the initial undeformed state (dV 0 ) and that after the deformation at later time (dV t ): The Lagrangian strain E L is defined in terms of the squared length of kdxk 2 and kdrk 2 as follows: -5) or, alternatively, Here, E L measures the total finite strain accumulated by the particle x and its vicinity from time t ¼ 0. We emphasize that the definition of Lagrangian strain (equation A-6) that we should use in seismological studies is directly related to the deformation tensor F.

Stress definitions
Forces in continuum mechanics are generally represented in terms of a stress tensor acting on a surface element. Let us first consider an infinitesimal-oriented surface elementn 0 dΣ 0 centered at a particle x with normaln 0 and size dΣ 0 in the undeformed state. At a later time t, the surface element is deformed ton t dΣ t at position r ( Figure A-1), where the relationship between the two states can be expressed asn Equation A-7 is the areal counterpart of equation A-4. Using this information, the Eulerian Cauchy stress tensor T E is defined as where df E denotes the instantaneous surface force exerted by all particles immediately in the front of the surface element on the particles immediately in the back. The corresponding Lagrangian Cauchy stress can then be defined with the help of equation A-1 as T L ðx; tÞ ¼ T E ðrðx; tÞ; tÞ. We note that T E naturally appears in the derivation of Eulerian momentum conservation law. The second definition of stress that offers convenience in the derivation of Lagrangian momentum conservation law is the first Piola-Kirchhoff stress T PK given by df E ¼n 0 dΣ 0 · T PK : (A-9) We emphasize that T PK gives df E acting across the deformed surface elementn t dΣ t at r in terms of the undeformed elementn 0 dΣ 0 . Therefore, T PK can be interpreted as the force per unit of undeformed area. Using equation A-7 together with equations A-8 and A-9, we can thus write The third definition of stress that is convenient for the derivation of the constitutive relation for perfectly elastic media is the second Piola-Kirchhoff stress T SK given by Here, T SK acts on the undeformed surface elementn 0 dΣ 0 and gives the transformed force df L acting upon the initial position x rather than the deformed position r. Using equations A-10 and A-11, we can also write

Lagrangian conservation laws
The preceding definitions of stress and strain can be used in the derivation of various conservation laws from either Lagrangian or Eulerian viewpoints. However, the details of such derivation are unimportant for the objective of this paper, and we refer readers to Chapter 2 of Dahlen and Tromp (1998). We only summarize the important final results from Lagrangian conservation laws: 1) Conservation of mass: which relates the instantaneous density of a particle ρ L to its initial value ρ 0 . 2) Conservation of momentum: where we neglect the effects from rotational reference frame and the gravitational forces. 3) Conservation of angular momentum: which indicates that the Lagrangian Cauchy stress is symmetric. It also follows from equations A-10 and A-12 that T PK is not symmetric because ðT PK Þ T ¼ F · T PK · F −T , but T SK is symmetric. 4) Conservation of energy: where U L denotes the Lagrangian internal energy density per unit mass that will serve as the basis for constructing constitutive relations for perfectly elastic material in the next section.

Constitutive relations
Let us now consider the Lagrangian formulation of the constitutive relations for perfectly elastic media, which governs how the material particles deform and return to their natural reference configuration by any applied stress. The most general perfect elastic material is the one whose Lagrangian internal energy density U L can be expressed in terms of the local instantaneous strain E L and the local entropy S L . The first dependency on E L stipulates that U L only changes when there is deformation (strain) but remains invariant due to a rigid rotation. The second dependency is generally neglected (∂ t S L ¼ D t S E ¼ 0) because seismic deformation is assumed to be isentropicadiabatic and reversible, where S L is unchanged. Therefore, we can write the derivative of the Lagrangian internal energy density with respect to time as where we multiply both sides by ρ 0 . Equating the conservation of energy law (equation A-16) and equation A-17 leads to which is expressed in terms of E L . Using the relationship between different measures of stress (equations A-10 and A-12), we can alternatively write Equations A-18-A-20 represent the most general form of constitutive relations for perfectly elastic media. We emphasize that from equation A-18, the Lagrangian strain E L is directly related to the second Piola-Kirchhoff stress T SK not the Lagrangian Cauchy stress T L nor the Eulerian Cauchy stress T E .

Linearized seismic deformation
In practice, we generally assume that seismic deformations are small in magnitude, which permits linearization of the constitutive relations, the conservation laws, and eventually the wave equation. Prior to such deformations the material is presumed to be in a state of mechanical equilibrium (∇ · T 0 ¼ 0), where there is zero strain and potentially nonzero initial static stress T 0 (intermediate configuration). In such equilibrium, there is no need to distinguish between Lagrangian and Eulerian quantities because that r ¼ x and F ¼ I. All the measures of stress, therefore, coincide and are equal to T 0 , which can be decomposed into the initial pressure p 0 and the initial deviatoric stress τ 0 as follows: where tr denotes the trace operator. Due to the static equilibrium condition, we can also deduce that ∇p 0 ¼ ∇ · τ 0 , and we can always express T 0 in terms of τ 0 alone. To characterize additional seismic deformation superposed on this initially stressed medium, we consider the Lagrangian description of motion expressed as rðx; tÞ ¼ x þ sðx; tÞ; (A-22) where s is the displacement away from equilibrium of particle x at time t, with particle velocity For any other physical quantity q, we can define first-order perturbations q E1 and q L1 as follows: q E ðr; tÞ ¼ q 0 ðrÞ þ q E1 ðr; tÞ; q L ðx; tÞ ¼ q 0 ðxÞ þ q L1 ðx; tÞ; (A-24) where q 0 ðrÞ ¼ q 0 ðxÞ denotes the zeroth-order initial value in the undeformed (intermediate) state. Because the magnitude of s is assumed to be small, we systematically neglect all terms of order jjsjj 2 or higher. This assumption allows us to write q E1 ðr; tÞ ¼ q E1 ðx; tÞ and q L1 ðr; tÞ ¼ q L1 ðx; tÞ; (A-25) which is correct to first order in s. Therefore, we may henceforth regard all zeroth-and first-order quantities as functions of the position x and the domain of all pertaining linearized equations will be associated with the volume of the undeformed equilibrium state. We will also henceforth omit the subscript x for conciseness. With regard to equation A-2, we can then rewrite it as which represents a linearized and integrated version of the material derivative. Its physical interpretation can be summarized as follows.
The first-order change q L1 experienced by an observer riding on a moving particle x is the summation of the first-order change q E1 at a fixed initial position x in the coordinate space, plus the change due to the particle displacement s with respect to the initial spatial gradient ∇q 0 . In view of the deformation gradient F, substituting equation A-22 into equation A-3 gives where the symmetric contribution ϵ is the infinitesimal strain tensor generally used in geophysical studies and the antisymmetric part ω is the infinitesimal rotation tensor given by A similar linearization result can also be obtained via the Eulerian framework, which leads to a common conclusion that there is no need to distinguish between the Lagrangian and Eulerian viewpoints in the case of infinitesimal (seismic) deformation. However, this is not true for variables that have zeroth-order values q 0 , e.g., density and stress, as we shall show subsequently.

Stress perturbations
Using the general definition in equation A-24, we can write the first-order perturbations to the Lagrangian and Eulerian Cauchy stress as follows: (A-33) where the perturbations themselves are related by (equation A-26) We can clearly observe that the perturbations of the Eulerian Cauchy stress from seismic deformation differ from its Lagrangian counterpart by a term that depends on the initial stress T 0 . Only when ignoring the effects from T 0 do the two become similar. Following the same procedure for T PK and T SK and using equations A-10, A-12, and A-27-A-30, we have Equations A-34, A-36, and A-37 represent the first-order perturbations in stress due to seismic deformations that we will use to establish the wave equation in the next section.

Wave equation under stress
To derive the equation of motion (wave equation) for seismic deformation, we consider the linearized version of the momentum conservation law (equation A-14). Substituting the expression for r (equation A-22) and T PK (equations A-35-A-37), we can write correctly to the first order in ksk wave equations under stress as shown in equation 1 in the main text.
We note that equation 1 consists of three equivalent wave equations with respect to the different choices of stress measure. They generally involve terms related to the initial static stress T 0 except in the case of the nonsymmetric T PK1 . Even though the wave equation with T PK1 appears the simplest at a glance, T PK1 ≠ ðT PK1 Þ T as follows from the conservation of angular momentum in equation A-15 and, thus, T PK1 may not be the most desirable measure of stress with which to work. Depending on the choice of stress measure, one also has to work with different constitutive relationships (equations A-18-A-20). Therefore, what we commonly refer to as "elastic moduli" or "stiffness coefficients" must conform to the type of considered stress. In the next section, we show how each choice of stress perturbation can be related to its corresponding linearized constitutive relation, which, in turn, governs how the material responds to small (seismic) deformation under the influence of the initial state of stress T 0 .

Linearized constitutive relations
In view of the linearized constitutive relations, it is extremely crucial that we calculate the Lagrangian internal energy density U L correct to the second order of ksk. This is essential because the general constitutive relations (equations A-18-A-20) are associated with the derivative of U L with respect to E L or F. In other words, because we want the derivative of U L to be of order one, U L must be of order two. Therefore, this requirement does not allow us to use the linearized approximation of E L ¼ ϵ (equation A-32) in the derivation of linearized constitutive relations from U L .
Because we will also be working with higher order tensors in constitutive relations, it is convenient to consider index notation for such quantities. We will be, from this point onward, switching between the previous invariant notation and index notation as appropriate. Equation A-31 can be expressed in index notation as E L ij ¼ 1 2 ð∂ i s j þ ∂ j s i Þ þ 1 2 ∂ i s k ∂ j s k : (A-38) Expanding the Lagrangian internal energy density U L in terms of E L to second order in ksk gives where ρ 0 U 0 is the zeroth-order term that represents the energy in the underformed reference state. The first-order term T 0 ∶E L ensures that all three stresses (T L , T PK , and T SK ) reduce to T 0 (equations A-18-A-20) correct to zeroth order in jjsjj. Differentiating equation A-39 with respect to E L , we can express the second Piola-Kirchhoff stress T SK from equations A-18 and A-35 as where it follows that -41) and Ξ ijkl has the symmetry as shown in equation 8 in the main text. The symmetry of the first two indices is due to the symmetry of the stress tensor (equation A-15). The symmetry of the last two indices is due to the symmetry of the strain tensor by definition (equation A-6). The final symmetry between the front and back pairs of indices is due to equation A-41, where the mixed second-order derivative remains the same regardless of the order of differentiation. At this point, we conclude our consideration on U L (equation A-39), where the exact definition of E L (equations A-31 and A-32) correct to the second order in ksk must be used. Now that the linear relation (equation A-40) is obtained with the correct definition of Ξ, it is permissible to use E L ≈ ϵ, which leads to equations 2-4 in the main text.
At this point, we can see that there are several possible measures of stress and each one is associated with a different fourth-rank tensor that governs the material's responses to small seismic deformation. Only in the absence of any initial static stress T 0 ¼ 0 do all measures of stress perturbations become equal (equations A-34 and A-36-A-37) and lead to the classic stress-strain (Hooke's law) relation written as T ¼ Γ∶ϵ; (A-42) where T can be any measure of stress perturbations and the elastic (stiffness) tensor has 21 independent coefficients and the property In the more general case of nonzero initial static stress T 0 ≠ 0 with superimposed incremental stress perturbations associated with seis-mic deformation, it is crucial to distinguish among different measures of stress and the corresponding elastic tensor because it is not perturbations in the Eulerian Cauchy stress T E1 that are related to ∇s by the elastic parameters but rather perturbations in the Lagrangian stress T L1 , or alternatively T PK1 or T SK1 . However, this requirement is actually not as troublesome as it may seem. In fact, we have great flexibility in the construction of Ξ as long as its symmetry properties (equation 8) are satisfied. In other words, there is a case of nonuniqueness for the stiffness tensor Ξ ( Barron and Klein, 1965;Thurston, 1965aThurston, , 1974Wallace, 1967;Dahlen, 1972). The corresponding expressions for Λ and Υ can then be obtained according to equations 3 and 4.