Elastic wave equation

Derivation

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

(1)(σ+M)+f=ρu¨,

where σ is the stress tensor, M is the moment tensor, f is the body force, ρ is the density, and u is the particle displacement.

The strain can be approximated as

(2)ϵij=12(iuj+jui),

and, for an isotropic medium, the stress simplifies to

(3)σij=λδijϵkk+2μϵij,

where λ and μ 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)(σ)+f=ρv˙,

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

(5)ϵ˙ij=12(ivj+jvi),

and (3) as

(6)σ˙ij=λδijϵ˙kk+2μϵ˙ij.

The staggered grid has vx positioned on (y,x) grid cells, vy at (y+12,x+12), σxx and σyy at (y,x+12), and σxy at (y+12,x). The vx and σii (where ii refers to xx or yy) components in the first row of the model, and the vy and σii 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 (lamba, mu, and buoyancy) are at the same locations as vx:

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

Staggering also occurs in time, with the stress (σ) components defined at time steps t while velocity components are at times t+12. Staggering the components in time and space, we obtain

vyt+12=vyt12+By+12,x+12ht(yσyy+xσxy+fyt)vxt+12=vxt12+By,xht(xσxx+yσxy+fxt)σyyt=σyyt1+ht((λy,x+12+2μy,x+12)yvy+λy,x+12xvx)σxxt=σxxt1+ht((λy,x+12+2μy,x+12)xvx+λy,x+12yvy)σxyt=σxyt1+ht(μy+12,x(xvy+yvx)),

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

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

x~=x+ψ.

For example,

y~σyy=yσyy+yσyym,

where yσ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.

yσyy,m=ay+12yσyy,m+by+12yσyyxσxx,m=axxσxx,m+bxxσxxyσxy,m=ayyσxy,m+byyσxyxσxy,m=ax+12xσxy,m+bx+12xσxyyvy,m=ayyvy,m+byyvyxvy,m=axxvy,m+bxxvyyvx,m=ay+12yvx,m+by+12yvxxvx,m=ax+12xvx,m+bx+12xvx

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

(vyt+12vxt+12yσyy,mtxσxx,mtyσxy,mtxσxy,mtrvyt12rvxt12)=(Byxht(1+by)y0Byxht(1+bx)x10Byxhtay00Byxhtax0Byxht(1+bx)xByxht(1+by)y010ByxhtaxByxhtay0byy0000ay0000bxx0000ax0000byy0000ay000bxx00000ax000δry000000000δrx0000)(σyytσxxtσxytvyt12vxt12yσyy,mt1xσxx,mt1yσxy,mt1xσxy,mt1)

and

(σyyt+1σxxt+1σxyt+1yvy,mt+12xvy,mt+12yvx,mt+12xvx,mt+12rpt)=((λyx+2μyx)ht(1+by)yλyxht(1+bx)x100(λyx+2μyx)htay00λyxhtaxλyxht(1+by)y(λyx+2μyx)ht(1+bx)x010λyxhtay00(λyx+2μyx)htaxμyxht(1+bx)xμyxht(1+by)y0010μyxhtaxμyxhtay0byy0000ay000bxx00000ax000byy00000ay00bxx000000ax00δrpδrp00000)(vyt+12vxt+12σyytσxxtσxytyvy,mt12xvy,mt12yvx,mt12xvx,mt12)

where, for conciseness, half grid cell shifts are represented by a prime, for example Byx is the buoyancy at locations (y+12,x+12). δry, δrx, and δrp, are the locations of vy, vx, and pressure receivers, respectively. In order for the velocity receiver data to cover the same time range as the input sources, the recordings are shifted by half a time step before being returned to the user.

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 σyy=0 and σxy=0 there. Different methods have been proposed to implement this. Deepwave currently uses the W-AFDA approach, and applies it to all four edges so that setting the PML width to zero on any of them will result in a free surface. W-AFDA uses non-symmetric finite difference stencils near the edges so that no values beyond the free surface are needed for calculations, and imposes the constraint that the traction is zero by setting the relevant stresses to zero in the calculations.