A free boundary condition for the elastic wave equation in the pseudospectral method

Abstract

Free surface boundary conditions constitute an usefu1 idealization of the behavior of the Earth’s surface when stimulated by seismic waves, either in its solid or liquid portions. Numerical methods that rely on direct discretization of the wave equation on a grid, such as finite differences (Virieux, 1986; Levander, 1988) and pseudo-spectral (Gazdag., 1982; Kosloff and Baysa1, 1982; Kosloff et al. 1984), seem to show a remarkable difficulty in addressing the apparently simple question of zeroing out some stress components on a free boundary. Some authors report that these difficulties increase with the size of the differential operators (Kosloff et al. 1990), becoming almost intractable in the pseudo-spectral method, when operators become as big as the whole grid and wrap the wavefield around it. Soon it was realized that the mere annihilation of stresses in a boundary lead to instability in propagation (Bayliss et al. 1986). Ingenious, but sometimes very involved or machine-costly, attempts to simulate a free boundary in the context of grid methods have then included special propagator stencils for these bounties (Vidale and Clayton 1986), the padding of the grid with zeros beyond the free surface (a “heavy vacuum”, or a material with zero stiffiess but finite density) (Kosloff et al. 1984) and using non-uniform grid spacing close to the reflecting boundary (Kosloff et al. 1990). We present a simple and cheap scheme to simulate a free boundary for the elastic wavefield that performs a decent job in matching traveltimes, polarities, phases and amplitudes and is also stable. Our quite standard, somewhat dispersive, but cost-effective propagation algorithm is pseudo-spectral in space and second order finite difference in time (Witte, 1989). In order to comply with public domain, generally available, simple and portable radix 2 FFT algorithms, but keep the differential operators as short and local as possible, it uses a 2D staggered grid and staggered derivatives (Fomberg, 1989; Witte, 1989). Particle displacements are used as primary variables (Witte, 1989), as opposed to velocities and stresses (Kosloff et al. 1984, 1990). We suggest that the pseudo-spectral method is particularly sensitive to boundary conditions that interfere directly with the primary variables being propagated, but more compliant The 2D elastic wave equation in an homogeneous isotropic medium is written in terms of displacement gradient and time derivative, densities and Lame’ parameters of the medium (or alternatively P and S wave velocities). This process can be performed in two stages, that involve the combination of Newton’s Second Law, where displacement second time derivative and stress gradient are present and the constitutive equation of linear elastic materials, where stress and displacement gradient are present. Gradients are calculated by Fourier transforms (pseudo-spectral method) in the staggered grid of Fig. 1 using staggered derivatives (Witte, 1989). In Fig. 1 “s” stands for stresses, m for moment tensors, u for displacements and “d” for densities. As already mentioned, staggering is needed only to keep derivative operators short in the world of radix 2 FFTs and power of two sized grids. Second order time derivatives are calculated by the ordinary centered scheme. As is well known, the stability of this algorithm is controlled by the maximum velocity of the medium and its ability to avoid aliasing by the minimum velocity of the medium (Kosloff et al., 1984). Hence, wide velocity dynamic ranges reduce the efficiency of the algorithm. Sources are introduced by means of the components of a moment tensor, with a given time signature (Witte: 1989). The sharp global spatial derivative operators allow us to use point sources (Witte, 1989). Algorithmically this procedure gives full access to both displacements as well as stresses in each time step. If the top of our grid is supposed to represent a free boundary. Z is the vertical downward and X the horizontal rightward axes, the stress components and are zeroed out in each time step along the grid rows shown in Fig. 2, where a 4x4 mesh is represented. Optionally, the telegraphy equation absorbing boundary conditions (Cerjan et al. 1985; Levander, 1985) can be applied along the sides and bottom of the grid, significantly reducing the wrapping around of the wavefield there and improving the results. We point out that the free boundary condition at the top of the grid works independently of the absorbing boundary conditions at the bottom and sides of the grid, as some of the following examples show.