Derivation
As in the derivation of the scalar wave equation, we use Newton’s second law.
(1)
where is the stress tensor, is the moment tensor, is the body force, is the density, and is the particle displacement.
The strain can be approximated as
(2)
and, for an isotropic medium, the stress simplifies to
(3)
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)
where is the particle velocity. We also rewrite (2) as
(5)
and (3) as
(6)
The staggered grid has positioned on grid cells, at , and at , and at . The and (where refers to or ) components in the first row of the model, and the and 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 :
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 while velocity components are at times . Staggering the components in time and space, we obtain
where is the buoyancy (reciprocal of the density), and is the finite difference time step interval.
To use the C-PML method to implement an absorbing boundary, we replace spatial derivatives by
For example,
where 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 and .
Combining these, we can express a time step in matrix form as
and
where, for conciseness, half grid cell shifts are represented by a prime, for example is the buoyancy at locations . , , and , are the locations of , , 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 dimension is depth, then for the top surface to be a free surface we need and 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.