Tensorial Elastodynamics for Isotropic Media on Vertically Deformed Meshes

Jeffrey Shragge⇤ , University of Western Australia SUMMARY

Solutions of the elastodynamic wave equation are sometimes required in reverse-time migration and full waveform inversion applications involving vertically deformed meshes that incorporate irregular free-surface topography and internal mesh boundaries (e.g., water bottom). By following a tensorial approach, I develop a system of analytic equations governing elastodynamics and the free-surface boundary condition. These equations are both straightforward and efficient to solve using a velocity-stress formulation with mimetic finite-difference (MFD) operators implemented on a fully staggered grid. Examples show that high-quality elastic wavefields can be generated for free surfaces exhibiting significant topographic relief.

INTRODUCTION Numerical solutions of elastodynamics increasingly are being used for a wide variety of 3D earth imaging problems within the earthquake seismology, microseismic and exploration communities. One of the main drivers is a need to accurately model and reproduce the full-wavefield elastic phenomena that are required to undertake a number of data-domain imaging and inverse problems such as elastic reverse-time migration and full waveform inversion, respectively. While the majority of elastodynamic modeling efforts are undertaken on Cartesian solution domains, there have been numerous examples of computing such responses on more generalized non-Cartesian meshes (Ohminato and Chouet, 1997; Hestholm, 1999; Hestholm and Ruud, 2002; Hestholm et al., 2006; Appel¨o and Petersson, 2009). Of particular interest here are solutions computed on vertically deformed meshes that straightforwardly adapt to irregular computational geometries encountered in real-world applications. Such scenarios include handling free-surface topography and/or irregular internal boundaries that arise in land and marine applications. Finite-difference (FD) methods are straightforward to implement on simple, rotated or fully staggered grids (FSGs), and have compact stencils for derivative operators for fairly high numerical orders of accuracy that are well suited for use in highly parallelized and efficient accelerator-based numerical solvers. However, a key challenge associated with applying this approach on irregular grids is devising FD methods with stable and accurate numerical solutions on computational meshes characterized by moderate to strong roughness and skewness. These challenges become increasingly acute in boundary regions where applying higher-order stencils is problematic. The development of stable and high-order FD approaches able to address these challenges in the past decades [e.g., mimetic finite difference (MFD) operators (Castillo and Miranda, 2013; de la Puente et al., 2014, and references therein) applied on FSGs (Lebedev, 1964; Lisitsa and Vishnevskiy, 2010)] arguably

has lead to renewed interest in applying FD operators for modeling elastodynamics in scenarios involving vertically deformed meshes. To generate these solutions practitioners generally follow one of two mathematical strategies: the tensorial or the chain-rule approach. Tensorial approaches pose elastodynamics problems and subsequent numerical solutions directly on a generalized coordinate mesh. Though this requires only the same number of wavefield derivatives as Cartesian solutions (Komatitsch et al., 1996), additional memory variables are required to account for the non-Cartesian metric tensor and corresponding Christoffel symbols. On the other hand, chain-rule approaches first compute the spatial derivatives in the generalized coordinate system before using the chain-rule to transform the result back onto the Cartesian mesh. These operations have an associated additional computational and memory overhead relative to solutions computed in a Cartesian solution. Relative to the tensorial approach, chain-rule applications have higher computational complexity, but lower memory requirements. Thus, during an era of memory-bound computing, the chain-rule approach tended to be more widely applied. This work re-examines tensorial formulation of the elastodynamic wave equation using a new mathematical strategy. Instead of specifying a numerical mesh and directly computing the corresponding geometric memory variables, I develop analytic expressions for the equations governing isotropic elastodynamics for all admissible (i.e., C2 differentiable) vertically deformed coordinate systems. Importantly, the development of analytic expressions reduces the memory requirements associated with geometric variables to simple expressions of one dimension. An additional benefit is the specification of an analytic expressions for the free-surface boundary condition for all admissible vertically deformed grids, which is generally straightforward to implement in practice. The paper begins with the tensorial theory of elastodynamics, and is followed by a development of analytic equations governing elastodynamics on a family of vertically deformed meshes using the example of a 2D topographic coordinate system. I then describe the mimetic FD and fully staggered grid (MFD+FSG) approach used to solve the analytic expressions using a finite-difference time-domain (FDTD) technique. I then present two topographic numerical 2D examples for (1) a dipping planar surface and (2) an undulating surface with 585 m of overall relief. Finally, while only 2D examples are presented herein, the approach is generally valid in 3D, which will be reported in a follow-on manuscript. TENSORIAL ELASTODYNAMICS WAVE EQUATION This section presents the tensorial formulation of elastodynamics for a generalized coordinate system, x = [x 1 , x 2 , x 3 ], which represents a computational domain that is linked to an underlying Cartesian mesh x = [x1 , x2 , x3 ] through a set of unique mappings. The goal is to develop generalized expressions for

Tensorial Elastodynamics the two key equations linking the stress and particle velocity field variables: the linear conservation of momentum and the material constitutive relationship. Appendix A presents an overview of some of the differential geometry and tensorial calculus used in the following treatment. The tensorial formulation of the linear conservation of momentum (LCM) (McConnell, 1957; Brillouin, 1964; Komatitsch et al., 1996) is given by: j

r v˙i = — j si + fi ,

(1)

where r = r(x ) is density; v˙i = v˙i (x ) is the particle acceleration (i.e., temporal derivative of particle velocity) in the ith direction; a dot over a variable indicates temporal derivative; — j is the covariant derivative defined in x -coordinates; fi = fi (x ) j j is a covariant component of force vector; and si = si (x ) is the mixed (1,1) stress tensor. Inserting the expression for the covariant derivative of the mixed tensor field (see equation A-8) into the 3D elastodynamics wave equation yields j

∂ si j + G jk sik ∂x j

2D TOPOGRAPHIC COORDINATES In this section I derive an analytic 2D elastodynamic wave equation for a mesh with a topographic upper boundary using the following straightforward coordinate mapping  1  1 x x +T = , (7) x2 x2 where T = T (x 2 ) = T (x2 ) is a generalised topographic function. Using a notation where T2 ⌘ ∂∂xT2 , the inverse metric tensor, gi j , for mapping in equation 7 may be written h i  1 + T22 T2 gi j = , T2 1

and the only non-zero Christoffel symbol is G122 = ∂∂ xT22 . After introducing these geometric variables in the above equations, and through some algebraic manipulations, one can show that the 2D linear conservation of momentum equation under this coordinate transformation is

(2)

r v˙1 =

∂ s11 ∂ s12 + + f1 , ∂ x1 ∂ x2

where Gkij is a Christoffel symbol of the second kind (see equation A-7). The constitutive relationship between the stress and strain tensors for an isotropic elastic medium is

r v˙2 =

∂ s21 ∂ s22 ∂s2 + + T2 1 ∂ x1 ∂ x2 ∂ x2

r v˙i =

j

s˙ i = l

j

Gkji sk + fi ,

j j di e˙kk + 2µ e˙i ,

(3)

j where l and µ are Lam´e parameters; and e˙i is the temporal derivative of the mixed (1,1) infinitesimal strain tensor (IST). An expression for the latter can be derived by raising the index of the IST in its covariant basis e˙in : j e˙i = gn j e˙in ,

(4)

resulting in the following isotropic stress tensor j j s˙ i = l di gkn e˙nk + 2µgk j e˙ik .

The covariant representation of the IST is given by  1 1 ∂ vk ∂ vi e˙ik = [—i vk + —k vi ] = + Glik vl . 2 2 ∂xi ∂xk

(5)

(6)

Combining equations 3 and 6 yields an expression for the time derivative of the generalized coordinate stress tensor. Equations 2 and 5-6 represent the two key relationships used to simulate numerical elastodynamic solutions using the velocitystress numerical approaches. The RHS of equations 2 and 6 contain spatial derivatives of the stress and velocity variables, respectively; however, as written they also both dependent on j si and vi , which presents a challenge for staggered-grid implementations. In practice, by defining a specific x -coordinate transformation along with its associated Christoffel symbols, one can apply the chain rule to shift the derivative operators from Gijk back onto v1 , leaving a system of equations with the RHS expressed as function of spatial derivatives and not explicitly on the stresses and particle velocities.

(8)

(9a) ∂ T2 s12 + f2 , ∂ x2

and that the 2D isotropic stress-strain relationship is  (l + 2µ) 1 + T22 ∂ v1 s˙ 11 µ ∂ v1 l + µ = T + 2 l l l ∂x2 l ∂x1 ∂ (v2 T2 u1 ) + , ∂ x2 ⌘ ∂v s˙ 21 ⇣ ∂ v2 ∂ 1 + T22 v1 2 = 1 + T22 2T2 2 + , 1 µ ∂ x2 ∂x ∂x

(9b)

∂ v2 ∂x1 (10a) (10b)

s˙ 12 ∂ v1 ∂ v2 ∂ v1 = + 2T2 1 , (10c) µ ∂x2 ∂x1 ∂x  s˙ 22 µ ∂ v1 l + µ ∂ v2 l + 2µ ∂ (v2 T2 v1 ) = T2 + l l ∂x2 l ∂x1 l ∂ x2 " # 2 ⇣ ⌘ µ ∂ T2 v1 ∂ v1 ∂ v1 + T22 1 + 1 + T22 . (10d) l ∂ x1 ∂x ∂x1

These equations have been verified using the Mathematica software package. Equations 9 and 10 represent analytic expressions of the 2D elastodynamics on a simple topographic mesh. Both equations have a temporal derivative on the LHS and are solely formed by spatial derivatives on the RHS. Thus, these expressions in an appropriate format for generating the MFD+FSG numerical scheme discussed below. The free-surface boundary condition (FSBC) for this topographic coordinate system may be obtained by directly solving j

si n j

x 1 =0

= 0,

(11)

where n j is the normal vector at the free surface (i.e., x 1 = 0). From the relationships in equations 10a-d, one may solve these equations in terms of spatial derivatives of particle velocity: ∂ v1 = zT ∂x1

and

∂ v2 = ∂x1

∂ v1 + 2T2 zT , ∂x2

(12)

Tensorial Elastodynamics where zT =

l 1



l ∂ v1 ∂ (v2 T2 u1 ) T2 2 + . ∂x ∂x2 T22 + 2µ

(13)

The zT field represents a correction filter that accounts for n j not necessarily being orthogonal to the computational grid. Finally, injecting an effective Cartesian force source fi into the x -coordinate mesh requires applying the coordinate transformation of a vector field (equation A-4)   v1 v1 = , (14) v2 x v2 + T2 v1 x which is the wavefield source applied in following examples. Similarly, extracting modeled particle velocities from the computational grids requires undergoing the inverse transformation back to Cartesian variables for visualization purposes. NUMERICAL IMPLEMENTATION I implement the equations of elastodynamics for the 2D simple topographic coordinate system presented above using a MFD strategy (Castillo and Miranda, 2013) on a fully staggered grid (FSG). The O(Dx4 ) gradient operators used herein are equivalent to those presented in de la Puente et al. (2014), which provides an excellent overview of the use of MFD in 3D elastodynamics simulations. MFD operators increasingly are being used to develop (more) uniform FD approximations of higherorder accuracy throughout the entire computational domain not just within the domain interior away from boundary regions where using FD-based stencil introduces numerical challenges. In the boundary region, MFD operators satisfy the discrete equivalent of underlying conservation laws, and are designed to prevent erroneous numerical wavefield flux from entering the computational domain. When progressing from the boundary into the interior, the MFD coefficients become equivalent to those from standard Taylor-series expansions. The MFD operators are applied within a fully-staggered grid (FSG) (Lebedev, 1964; Lisitsa and Vishnevskiy, 2010). In 2D (3D) applications a FSG system requires computing numerical solutions on two (four) sub-grids that are coupled through crossderivative terms. In this 2D velocity-stress implementation, I locate the particle velocity and stress variables on the regular and spatially staggered mesh of each sub-grid, respectively. Overall, evaluating the 2D (3D) linear conservation of momentum requires computing 2x5 (4x13) partial derivatives, while the stress-strain relationship requires 2x5 (4x13) partial derivative computations. The 2D (3D) MFD+FSG implementation requires holding 2x6 (4x12) velocity/stress variables in memory. The memory requirement for 2D (3D) geometric variables is 1D (2D) and comparably negligible. Assuming initial quiescence, I first update the stress variables from equation 10 at the half-time step using a standard secondorder temporal FD approximation. I then inject the point force source and update of the velocity variables at the full time step again using the same second-order temporal FD approximation. The velocity-stress updating process is subsequently

continued until reaching the maximum propagation time. Multicomponent wavefield snapshots are projected back to an orthogonal coordinate system, output at regular time intervals, and interpolated back to a Cartesian domain for visualization purposes using sinc-based operators (Shragge, 2014) 2D TOPOGRAPHIC EXAMPLES The first example is a flat interface dipping at a 30 angle, which represents a tilted Garvin test commonly used for evaluating the accuracy of implemented elastodynamic schemes. I use constant P- and S-wave velocities of Vp = 3.0 km/s and Vs = 1.5 km/s, and set the density to r = 1.0 g/cc using a uniform 6.25 m mesh sample interval. The 25 Hz Ricker wavelet is injected as a force source oriented in the x1 (vertical) direction two grid points in from the surface. The FSBC has been applied according to the analytic expression above. The upper (lower) panel shows the v1 (v2 ) wavefield component after interpolating the solution back to Cartesian. Note that because the point source is not injected perpendicularly to the topographic surface, the resulting free-surface response is asymmetric.

Figure 1: Wavefield snapshot for the v1 (top) and v2 (bottom) Cartesian components for a 30 dipping interface using a point force source oriented in the x1 direction. The second example shows the 2D elastodynamic responses computed on simple topographic mesh defined by a randomly generated topography with a 585 m of total relief. The elastic model components and the source force are the same as in the previous example. Figure 2 shows the v1 and v2 wavefield components extracted at t = 0.43 s. The undulating free surface introduces a significant amount of topographic scattering and a more complicated free-surface response. Overall,

Tensorial Elastodynamics these examples demonstrate that accurate and stable impulse responses can be achieved with the equations of elastodynamics and MFD+FSG FDTD implementation described above.

APPENDIX A TENSOR CALCULUS REVIEW ⇥ ⇤ A generalized 3D coordinate system x = x 1 , x 2 , x 3 can be ⇥ 1 2 3⇤ related to an underlying 3D Cartesian field x = x , x , x according to the forward and inverse mappings: x i = x i (x) and

xi = xi (x ),

(A-1)

where a raised (lower) index indicates a contravariant (covariant) basis. The geometric connection between the two coordinate systems is the symmetric covariant metric tensor gi j : gi j =

∂ xk ∂ xk , ∂xi ∂x j

(A-2)

where summation notation is used over repeated indices. The inverse metric tensor, gi j , is found by raising the indices of gi j gi j g jk = dki = dik

(A-3)

dki

where is Kroenker delta function. When writing the physical equation in a generalized sense, one must appropriately define a vector and higher-order tensors in x -coordinates by applying a change-of-variables (COV) transformation. The COV of covariant vector vi is: vi |x = Figure 2: Two sets of wavefield snapshots for a simple topographic coordinate system that conforms to the topographic surface (black line). The upper and lower panels (again v1 and v2 ) show wavefield snapshots at t = 0.63 s.

CONCLUSIONS By deriving analytic equations of elastodynamics for all permissible vertically deformed coordinate meshes, it is possible to specify a set of expressions that directly incorporate free-surface topography and allow for straightforward implementations of the free-surface boundary condition. By using a MFD+FSG FDTD strategy, it is possible to generate highquality wavefield responses that exhibit the expected complex elastic full-wavefield behavior. Thus, these modeling operators would be appropriate for elastic RTM and FWI applications involving vertically deformed meshes. Finally, while not illustrated here, additional analysis has demonstrated that it is straightforward to develop expressions for elastodynamics in 3D vertically deformed coordinates; however, this extension is left for a more detailed development in a future manuscript. ACKNOWLEDGMENTS I acknowledge the support of Woodside Pty Ltd. and the UWA:RM consortium sponsors. Analytic results were verified using Mathematica. Reproducible numerical examples were generated using the Madagascar package.

∂xj vj , ∂xi x

(A-4)

where the evaluation domains are explicitly identified. The COV of the mixed (1,1) stress tensor is: s ij

x

= gip

∂ xk ∂ xl k s ∂x p ∂x j l

x

(A-5)

.

Applying differential operators to the tensor fields requires accounting for spatially varying coordinate geometry. Of importance to this study is the covariant derivative — j of a covariant vector field vi . This operator is expressed through “semi-colon” notation (i.e., — j vi ⌘ vi; j ) that is related to the partial derivative operator through vi; j =

∂ vi ∂x j

Gkij vk ,

(A-6)

where Gkij is a Christoffel symbol of the second kind: 1 Gkij = gkm 2

✓

∂ gmi ∂ gm j + ∂x j ∂xi

∂ gi j ∂xm

◆

.

(A-7)

Christoffel symbols are symmetric about the lower indices (i.e., Gkij = Gkji ) and are not tensors since they are not invariant under coordinate transformation. The covariant derivative of a mixed a , given by (1,1) tensor, sb;c a sb;c =

∂ sba + Gacd sbd ∂xc

Gdcb sda .

(A-8)