Elastic wave equation

Derivation

As in the derivation of the scalar wave equation, we use Newton’s second law.

(1)\[\nabla\cdot\left(\mathbf{\sigma} + \mathbf{M}\right) + \mathbf{f} = \rho \ddot{\mathbf{u}},\]

where \(\mathbf{\sigma}\) is the stress tensor, \(\mathbf{M}\) is the moment tensor, \(\mathbf{f}\) is the body force, \(\rho\) is the density, and \(\mathbf{u}\) is the particle displacement.

The strain can be approximated as

(2)\[\epsilon_{ij} = \frac{1}{2}\left(\partial_iu_j + \partial_ju_i\right),\]

and, for an isotropic medium, the stress simplifies to

(3)\[\sigma_{ij} = \lambda\delta_{ij}\epsilon_{kk} + 2\mu\epsilon_{ij},\]

where \(\lambda\) and \(\mu\) are the Lamé parameters.

Implementation

As moment tensor sources can be incorporated into body forces (but may need to be split over multiple adjacent grid cells), we will only use body sources from now on. For example, an explosion, for which the moment tensor is diagonal, can be approximated by negative body force terms in the x and y cells immediately before the source and positive terms in the x and y cells immediately after the source.

Deepwave uses the velocity-stress formulation of the elastic wave equation on a staggered grid. To obtain the velocity-stress formulation, we rewrite (1) as

(4)\[\nabla\cdot\left(\mathbf{\sigma}\right) + \mathbf{f} = \rho \dot{\mathbf{v}},\]

where \(v\) is the particle velocity. We also rewrite (2) as

(5)\[\dot\epsilon_{ij} = \frac{1}{2}\left(\partial_iv_j + \partial_jv_i\right),\]

and (3) as

(6)\[\dot\sigma_{ij} = \lambda\delta_{ij}\dot\epsilon_{kk} + 2\mu\dot\epsilon_{ij}.\]

The staggered grid has \(\sigma_{yy}\) and \(\sigma_{xx}\) positioned on \((y, x)\) grid cells, \(v_y\) at \((y+\frac{1}{2}, x)\), \(v_x\) at \((y, x+\frac{1}{2})\), and \(\sigma_{xy}\) at \((y+\frac{1}{2}, x+\frac{1}{2})\). The \(v_y\) and \(\sigma_{xy}\) components in the last row of the model, and the \(v_x\) and \(\sigma_{xy}\) components in the last column, are not used. This is depicted in the diagram below, where lowercase means that the component is not used. In this diagram, equals symbols and square brackets are used to denote the edges of the computational domain (so, if the PML width is zero, this will be the surface), while pipes and hyphens correspond to grid cell boundaries. The model parameters (lambda, mu, and buoyancy) are at the same locations as \(\sigma_{yy}\) and \(\sigma_{xx}\):

o--------->x
| SII=VX==|=SII=VX==|=SII vx
| VY  SXY | VY  SXY | VY  sxy
v [-------------------]------
y SII VX  | SII VX  | SII vx
  VY  SXY | VY  SXY | VY  sxy
  [-------------------]------
  SII=VX==|=SII=VX==|=SII vx
  vy  sxy | vy  sxy | vy  sxy

Staggering also occurs in time, with the stress (\(\sigma\)) components defined at time steps \(t\) while velocity components are at times \(t+\frac{1}{2}\). Staggering the components in time and space, we obtain

\[\begin{split}\begin{align} v_y^{t+\frac{1}{2}} &= v_y^{t-\frac{1}{2}} + B^{y+\frac{1}{2}, x}h_t\left(\partial_y\sigma_{yy} + \partial_x\sigma_{xy} + f_y^t\right) \\ v_x^{t+\frac{1}{2}} &= v_x^{t-\frac{1}{2}} + B^{y, x+\frac{1}{2}}h_t\left(\partial_x\sigma_{xx} + \partial_y\sigma_{xy} + f_x^t\right) \\ \sigma_{yy}^{t} &= \sigma_{yy}^{t-1} + h_t\left(\left(\lambda^{y, x}+2\mu^{y, x}\right)\partial_yv_{y} + \lambda^{y, x}\partial_xv_{x}\right)\\ \sigma_{xx}^{t} &= \sigma_{xx}^{t-1} + h_t\left(\left(\lambda^{y, x}+2\mu^{y, x}\right)\partial_xv_{x} + \lambda^{y, x}\partial_yv_{y}\right)\\ \sigma_{xy}^{t} &= \sigma_{xy}^{t-1} + h_t\left(\mu^{y+\frac{1}{2}, x+\frac{1}{2}}\left(\partial_xv_{y} + \partial_yv_{x}\right)\right), \\ \end{align}\end{split}\]

where \(B\) is the buoyancy (reciprocal of the density), and \(h_t\) is the finite difference time step interval.

To use the C-PML method to implement an absorbing boundary, we replace spatial derivatives by

\[\frac{\partial}{\partial \tilde{x}} = \frac{\partial}{\partial x} + \psi.\]

For example,

\[\partial_{\tilde{y}}\sigma_{yy} = \partial_y\sigma_{yy} + \partial_y\sigma_{yy}^m,\]

where \(\partial_y\sigma_{yy,m}\) is a “memory” variable, an auxiliary wavefield that is needed for the calculation. We require one of these memory variables for each spatial derivative, resulting in the following eight auxiliary wavefields that depend on the PML profiles \(a\) and \(b\).

\[\begin{split}\begin{align} \partial_y\sigma_{yy,m} &= a^{y+\frac{1}{2}}\partial_y\sigma_{yy,m} + b^{y+\frac{1}{2}} \partial_y\sigma_{yy} \\ \partial_x\sigma_{xx,m} &= a^{x+\frac{1}{2}}\partial_x\sigma_{xx,m} + b^{x+\frac{1}{2}} \partial_x\sigma_{xx} \\ \partial_y\sigma_{xy,m} &= a^{y}\partial_y\sigma_{xy,m} + b^{y} \partial_y\sigma_{xy} \\ \partial_x\sigma_{xy,m} &= a^{x}\partial_x\sigma_{xy,m} + b^{x} \partial_x\sigma_{xy} \\ \partial_yv_{y,m} &= a^{y}\partial_yv_{y,m} + b^{y} \partial_yv_{y} \\ \partial_xv_{y,m} &= a^{x+\frac{1}{2}}\partial_xv_{y,m} + b^{x+\frac{1}{2}} \partial_xv_{y} \\ \partial_yv_{x,m} &= a^{y+\frac{1}{2}}\partial_yv_{x,m} + b^{y+\frac{1}{2}} \partial_yv_{x} \\ \partial_xv_{x,m} &= a^{x}\partial_xv_{x,m} + b^{x} \partial_xv_{x} \\ \end{align}\end{split}\]

Combining these, we can express a time step in matrix form as

\[\begin{split}\begin{pmatrix} v_y^{t+\frac{1}{2}} \\ v_x^{t+\frac{1}{2}} \\ \partial_y\sigma_{yy, m}^{t} \\ \partial_x\sigma_{xx, m}^{t} \\ \partial_y\sigma_{xy, m}^{t} \\ \partial_x\sigma_{xy, m}^{t} \\ r_{v_y}^{t-\frac{1}{2}} \\ r_{v_x}^{t-\frac{1}{2}} \\ \end{pmatrix} = \begin{pmatrix} B^{y'x}h_t(1+b^{y'})\partial_{y'}&0&B^{y'x}h_t(1+b^{x}) \partial_{x}&1&0&B^{y'x}h_ta^{y'}&0&0&B^{y'x}h_ta^{x} \\ 0&B^{yx'}h_t(1+b^{x'})\partial_{x'}&B^{yx'}h_t(1+b^y)\partial_y&0&1&0&B^{yx'}h_ta^{x'}&B^{yx'}h_ta^y&0 \\ b^{y'}\partial_{y'}&0&0&0&0&a^{y'}&0&0&0 \\ 0&b^{x'}\partial_{x'}&0&0&0&0&a^{x'}&0&0 \\ 0&0&b^y\partial_y&0&0&0&0&a^y&0 \\ 0&0&b^x \partial_x&0&0&0&0&0&a^x \\ 0&0&0&\delta_{r_y}&0&0&0&0&0 \\ 0&0&0&0&\delta_{r_x}&0&0&0&0 \\ \end{pmatrix} \begin{pmatrix} \sigma_{yy}^t \\ \sigma_{xx}^t \\ \sigma_{xy}^t \\ v_y^{t-\frac{1}{2}} \\ v_x^{t-\frac{1}{2}} \\ \partial_y\sigma_{yy, m}^{t-1} \\ \partial_x\sigma_{xx, m}^{t-1} \\ \partial_y\sigma_{xy, m}^{t-1} \\ \partial_x\sigma_{xy, m}^{t-1} \\ \end{pmatrix} \label{elastic_matrix_v}\end{split}\]

and

\[\begin{split}\begin{pmatrix} \sigma_{yy}^{t+1} \\ \sigma_{xx}^{t+1} \\ \sigma_{xy}^{t+1} \\ \partial_yv_{y, m}^{t+\frac{1}{2}} \\ \partial_xv_{y, m}^{t+\frac{1}{2}} \\ \partial_yv_{x, m}^{t+\frac{1}{2}} \\ \partial_xv_{x, m}^{t+\frac{1}{2}} \\ r_p^t \\ \end{pmatrix} = \begin{pmatrix} \left(\lambda^{yx}+2\mu^{yx}\right)h_t(1+b^y)\partial_y&\lambda^{yx}h_t(1+b^x)\partial_x&1&0&0&\left(\lambda^{yx}+2\mu^{yx}\right)h_ta^y&0&0&\lambda^{yx}h_ta^x \\ \lambda^{yx}h_t(1+b^y)\partial_y&\left(\lambda^{yx}+2\mu^{yx}\right)h_t(1+b^x)\partial_x&0&1&0&\lambda^{yx}h_ta^y&0&0&\left(\lambda^{yx}+2\mu^{yx}\right)h_ta^x \\ \mu^{y'x'}h_t(1+b^{x'})\partial_{x'}&\mu^{y'x'}h_t(1+b^{y'})\partial_{y'}&0&0&1&0&\mu^{y'x'}h_ta^{x'}&\mu^{y'x'}h_ta^{y'}&0 \\ b^y\partial_y&0&0&0&0&a^y&0&0&0 \\ b^{x'}\partial_{x'}&0&0&0&0&0&a^{x'}&0&0 \\ 0&b^{y'}\partial_{y'}&0&0&0&0&0&a^{y'}&0 \\ 0&b^x\partial_x&0&0&0&0&0&0&a^x \\ 0&0&-\delta_{r_p}&-\delta_{r_p}&0&0&0&0&0 \\ \end{pmatrix} \begin{pmatrix} v_y^{t+\frac{1}{2}} \\ v_x^{t+\frac{1}{2}} \\ \sigma_{yy}^{t} \\ \sigma_{xx}^{t} \\ \sigma_{xy}^{t} \\ \partial_yv_{y, m}^{t-\frac{1}{2}} \\ \partial_xv_{y, m}^{t-\frac{1}{2}} \\ \partial_yv_{x, m}^{t-\frac{1}{2}} \\ \partial_xv_{x, m}^{t-\frac{1}{2}} \\ \end{pmatrix} \label{elastic_matrix_sigma}\end{split}\]

where, for conciseness, half grid cell shifts are represented by a prime, for example \(B^{y'x'}\) is the buoyancy at locations \((y+\frac{1}{2}, x+\frac{1}{2})\). \(\delta_{r_y}\), \(\delta_{r_x}\), and \(\delta_{r_p}\), are the locations of \(v_y\), \(v_x\), and pressure receivers, respectively. The receiver data cover the same time range as the corresponding input sources, so the velocity receiver recordings are from -half a time step earlier than the pressure recordings.

A “free surface” refers to a surface where the traction is zero. For example, if the \(y\) dimension is depth, then for the top surface to be a free surface we need \(\sigma_{yy}=0\) and \(\sigma_{xy}=0\) there. Different methods have been proposed to implement this. Deepwave uses the improved vacuum method approach, so setting the property models to zero at any point will result in a free surface there.