Abstract Modeling of seismic wave propagation plays a key role in almost every aspect of exploration seismology. Fundamentally, it provides us with a means of understanding the character of recorded seismic data. Although analytical or semianalytical solutions exist for a number of canonical models, these are often insufficient to explain the full range of phenomena that arise in complex heterogeneous earth models governed by, for example, anisotropic, viscoelastic, or poroelastic rheologies. Typical manifestations of such complex phenomena that are observed widely in surface seismic data include scattering, generation of multiples, or interface waves. Another application area of seismic modeling is in survey evaluation and design, where different acquisition geometries and subsurface model hypotheses are assessed to choose an optimal acquisition and processing strategy. Furthermore, modeling is at the core of many seismic data processing and imaging/inversion algorithms, including noise attenuation (removal of coherent noise by first modeling it), multiple suppression, migration, inversion, etc. Further up the seismic value chain, modeling also is very useful to guide final interpretation of an image, both to evaluate whether key features in the data are real or the result of (for instance) multiples, and to assess lithology or fluid content. Finally, generating synthetic data sets for research has proven extremely valuable for the exploration-seismics community in benchmarking and testing new processing algorithms, e.g., Marmousi (Versteeg, 1994), SMART JV Pluto model (Stoughton et al., 2001), the SEG salt model (Aminzadeh et al., 1995), and the SEAM initiative (Fehler, 2009; Symes et al., 2009). Thus, it should

Abstract This work presents a new constitutive law for linear viscoelastic and anisotropic media, to model rock behaviour and its effects on wave propagation. In areas with high dissipation properties (e.g. hydrocarbon reservoirs), the interpretation of seismic data based on the isotropic and purely elastic assumption might lead to misinterpretations or, even worse, to overlooking useful information. Thus, a proper description of wave propagation requires a rheology which accounts for the anisotropic and anelastic behaviour of rocks. The present model is based on the following mechanical interpretation; each eigenvector (eigenstrain) of the stiffness tensor of an anisotropic solid defines a fundamental deformation state of the medium. The six eigenvalues (eigenstiffnesses) represent the genuine elastic parameters. Since they are independent of the reference system, they have an intrinsic physical content. From this fact and the correspondence principle we infer that in a real medium the rheological properties depend essentially on six relaxation functions, which are the generalization of the eigenstiffnesses to the viscoelastic case. The existence of six or less complex moduli depends on the symmetry class of the medium. We probe the new stress-strain relation with homogeneous viscoelastic plane waves, and give expressions for the slowness, attenuation, phase velocity, energy velocity (wavefront) and quality factor of the different wave modes.

Abstract Real earth media disperse and attenuate propagating mechanical waves. This anelastic behavior can be described well by a viscoelastic model. We have developed a finite-difference simulator to model wave propagation in viscoelastic media. The finite-difference method was chosen in favor of other methods for several reasons. Finite-difference codes are more portable than, for example, pseudospectral codes. Moreover, finite-difference schemes provide a convenient environment in which to define complicated boundaries. A staggered scheme of second-order accuracy in time and fourth-order accuracy in space appears to be optimally efficient. Because of intrinsic dispersion, no fixed grid points per wavelength rule can be given; instead, we present tables, which enable a choice of grid parameters for a given level of accuracy. Since the scheme models energy absorption, natural and efficient absorbing boundaries may be implemented merely by changing the parameters near the grid boundary. The viscoelastic scheme is only marginally more expensive than analogous elastic schemes. The efficient implementation of absorbing boundaries may therefore be a good reason for also using the viscoelastic scheme in purely elastic simulations. We illustrate our method and the importance of accurately modeling anelastic media through 2-D and 3-D examples from shallow marine environments.

Abstract We address the basic theoretical and algorithmic aspects of memory-efficient implementation of realistic attenuation in the staggered-grid finite-difference modeling of seismic-wave propagation in media with material discontinuities. We show that if averaging is applied to viscoelastic moduli in the frequency domain, it is possible to determine anelastic coefficients of the averaged medium representing a material discontinuity. We define (1) the anelastic functions in a new way, as being independent of anelastic coefficients, and (2) a new coarse spatial distribution of the anelastic functions in order to properly account for material discontinuities and, at the same time, to have it memory efficient. Numerical tests demonstrate that our approach enables more accurate viscoelastic modeling than other approaches.

Abstract A numerical modelling method for wave propagation in a linear viscoelastic medium with singular memory is developed in this paper. For a demonstration of the method, the Cole–Cole model of viscoelastic relaxation is adopted here. A formulation of the Cole–Cole model based on internal variables satisfying fractional relaxation equations is applied. In order to avoid integrating and storing of the entire history of the variables, a new method for solving fractional differential equations of arbitrary order based on a set of secondary internal variables is developed. Using the new method, the velocity–stress equations and the fractional relaxation equations are reduced to a system of first-order ordinary differential equations for the velocities, stresses, primary internal variables as well as the secondary internal variables. The horizontal spatial derivatives involved in the governing equations are calculated by the Fourier pseudospectral (PS) method, while the vertical ones are calculated by the Chebychev PS method. The physical boundary conditions and the non-reflecting conditions for the Chebychev PS method are also discussed. The global solution of the first-order system of ordinary differential equations is advanced in time by the Euler predictor–corrector methods. For the demonstration of our method, some numerical results are presented.

Abstract Until recently, the term “elastic” usually implied two-dimensional (2-D) and isotropic. In this limited context, the divergence and curl operators have found wide use as wave separation operators. For example, Mora (1987) used them in his inversion method to allow separate correlation of P and S arrivals, although the separation is buried in the math and not obvious. Clayton (1981) used them explicitly in several modeling and inversion methods. Devaney and Oristaglio (1986) used closely related operators to separate P and S arrivals in elastic VSP data. With the current widespread interest in anisotropy, it seems useful to extend the wave-type separation concept to anisotropic media. We give a simple geometrical explanation of why divergence and curl are wave-type separation operators in the isotropic case and then show how to construct wave-type separation operators for general 2-D anisotropic media. We demonstrate the method on a heterogeneous strongly anisotropic finite-difference example. Extending existing isotropic 2-D algorithms based on wave-type separation to include anisotropy seems to be straightforward.

Abstract An algorithm is presented to solve the elastic-wave equation by replacing the partial differentials with finite differences. It enables wave propagation to be simulated in three dimensions through generally anisotropic and heterogeneous models. The space derivatives are calculated using discrete convolution sums, while the time derivatives are replaced by a truncated Taylor expansion. A centered finite difference scheme in cartesian coordinates is used for the space derivatives leading to staggered grids. The use of finite difference approximations to the partial derivatives results in a frequency-dependent error in the group and phase velocities of waves. For anisotropic media, the use of staggered grids implies that some of the elements of the stress and strain tensors must be interpolated to calculate the Hook sum. This interpolation induces an additional error in the wave properties. The overall error depends on the precision of the derivative and interpolation operators, the anisotropic symmetry system, its orientation and the degree of anisotropy. The dispersion relation for the homogeneous case was derived for the proposed scheme. Since we use a general description of convolution sums to describe the finite difference operators, the numerical wave properties can be calculated for any space operator and an arbitrary homogeneous elastic model. In particular, phase and group velocities of the three wave types can be determined in any direction. We demonstrate that waves can be modeled accurately even through models with strong anisotropy when the operators are properly designed.

Abstract We describe the application of the rotated staggered grid (RSG) finite-difference technique to the wave equations for anisotropic and viscoelastic media. The RSG uses rotated finite-difference operators, leading to a distribution of modeling parameters in an elementary cell where all components of one physical property are located only at one single position. This can be advantageous for modeling wave propagation in anisotropic media or complex media, including high-contrast discontinuities, because no averaging of elastic moduli is needed. The RSG can be applied both to displacement-stress and to velocity-stress finite-difference (FD) schemes, whereby the latter are commonly used to model viscoelastic wave propagation. With a von Neumann-style analysis, we estimate the dispersion error of the RSG scheme in general anisotropic media. In three different simulation examples, all based on previously published problems, we demonstrate the application and the accuracy of the proposed numerical approach.

Abstract Several staggered grid schemes have been suggested for performing finite-difference calculations for the elastic wave equations. In this paper, thedispersion relationships and related computational requirements for the Lebedev and rotated staggered grids for anisotropic, elastic, finite-differencecalculations in smooth models are analyzed and compared. These grids are related to a popular staggered grid for the isotropic problem, the Virieux grid. The Lebedev grid decomposes into Virieux grids, two in two dimensions and fourin three dimensions, which decouple in isotropicmedia. Therefore the Lebedev scheme will havetwice or four times the computational requirements, memory, and CPU as the Virieux gridbut can be used with general anisotropy. In twodimensions, the rotated staggered grid is exactlyequivalent to the Lebedev grid, but in three dimensions it is fundamentally different. The numericaldispersion in finite-difference grids depends on the direction of propagation and the grid typeand parameters. A joint numerical dispersion relation for the two grids types in the isotropic case is derived. In order to compare the computationalrequirements for the two grid types, the dispersion, averaged over propagation direction and mediumvelocity are calculated. Setting the parameters sothe average dispersion is equal for the two grids, the computational requirements of the two grid types are compared. In three dimensions, the rotated staggered grid requires at least 20% more memory for the field dataand at least twice as many number of floating point operations and memory accesses, so the Lebedev grid is more efficient and is to be preferred.

Abstract We introduce in this work a new spectral collocation scheme for the elastic wave equation transformed from Cartesian to generalized coordinates. Both the spatial derivatives of field variables and the metrics of the transformation are calculated by the Chebyshev pseudospectral method. The technique requires a special treatment of the boundary conditions based on one-dimensional characteristics perpendicular to the boundaries. An explicit Runge-Kutta time integration scheme is used for time marching. The numerical solution of Lamb’s problem (wave propagation over the surface of an elastic solid) requires two one-dimensional stretching transformations for each Cartesian direction of the 2-D Chebyshev grid. The results show excellent agreement between the numerical and analytical solutions, demonstrating the effectiveness of the 2-D differential operator and boundary treatment. The second example uses a 2-D transformation to simulate wave propagation over a smooth step discontinuity at the surface. The snapshots yield the wave pattern expected from such a structure.

Abstract Few problems in elastodynamics have a closed-form analytical solution. The others can be investigated with semianalytical methods, but in general one is not sure whether these methods give reliable solutions. The same happens with numerical techniques: for instance, finite difference methods solve, in principle, any complex problem, including those with arbitrary inhomogeneities and boundary conditions. However, there is no way to verify the quantitative correctness of the solutions. The major problems are stability with respect to material properties, numerical dispersion, and the treatment of boundary conditions. In practice, these problems may produce inaccurate solutions. In this paper, the study of complex problems with two different numerical grid techniques in order to cross-check the solutions is proposed. Interface waves, in particular, are emphasized, since they pose the major difficulties due to the need to implement boundary conditions. The first method is based on global differential operators where the solution is expanded in terms of the Fourier basis and Chebyshev polynomials, while the second is the spectral element method, an extension of the finite element method that uses Chebyshev polynomials as interpolating functions. Both methods have spectral accuracy up to approximately the Nyquist wave number of the grid. Moreover, both methods implement the boundary conditions in a natural way, particularly the spectra element algorithm. We first solve Lamb’s problem and compare numerical and analytical solutions; then, the problem of dispersed Rayleigh waves, and finally, the two-quarter space problem. We show that the modeling algorithms correctly reproduce the analytical solutions and yield a perfect matching when these solutions do not exist. The combined modeling techniques provide a powerful tool for solving complex problems in elastodynamics.

Abstract For the propagation of elastic waves in unbounded domains, absorbing boundary conditions (ABCs) at the fictitious numerical boundaries have been proposed. In this paper we focus on both first- and second-order ABCs in the framework of variational (weak) approximations, like those stemming from Galerkin method (or its variants) for finite element or spectral approximations (1). In particular, we recover first order conditions as natural (or Neumann) conditions, whereas we propose a penalty residual method for the treatment of second order ABCs. The time discretization is based on implicit backward finite differences, whereas we use spectral Legendre collocation methods set in a variational form for the spatial discretization (treatment of finite element or spectral element approximations is completely similar). Numerical experiments exhibit that the present formulation of second-order ABCs improves the one based on first-order ABCs with regard to both the reduction of the total energy in the computational domain, and the Fourier spectrum of the displacement held at selected points of the elastic medium. A stability analysis is developed for the variational problem in the continuous case both for first- and second-order ABCs. A suitable treatment of ABCs at comers is also proposed.

Abstract We develop a two-dimensional solver for the acoustic wave equation with spatially varying coefficients. In what is a new approach, we use a basis of approximate prolate spheroidal wavefunctions and construct derivative operators that incorporate boundary and interface conditions. Writing the wave equation as a first-order system, we evolve the equation in time using the matrix exponential. Computation of the matrix exponential requires efficient representation of operators in two dimensions and for this purpose we use short sums of one-dimensional operators. We also use a partitioned low-rank representation in one dimension to further speed up the algorithm. We demonstrate that the method significantly reduces numerical dispersion and computational time when compared with a fourth-order finite difference scheme in space and an explicit fourth-order Runge–Kutta solver in time.

Abstract We have developed a new hybrid numerical method for 3-D viscoelastic modelling of seismic wave propagation and earthquake motion in heterogeneous media. The method is based on a combination of the fourth-order velocity–stress staggered-grid finite-difference (FD) scheme, that covers a major part of a computational domain, with the second-order finite-element (FE) method which can be applied to one or several relatively small subdomains. The FD and FE parts causally communicate at each time level in the FD–FE transition zone consisting of the FE Dirichlet boundary, FD–FE averaging zone and FD Dirichlet zone. The implemented FE formulation makes use of the concept of the global restoring-force vector which significantly reduces memory requirements compared to the standard formulation based on the global stiffness matrix. The realistic attenuation in the whole medium is incorporated using the rheology of the generalized Maxwell body in a definition equivalent to the generalized Zener body. The FE subdomains can comprise extended kinematic or dynamic models of the earthquake source or the free-surface topography. The kinematic source can be simulated using the body-force term in the equation of motion. The traction-at-split-node method is implemented in the FE method for simulation of the spontaneous rupture propagation. The hybrid method can be applied to a variety of problems related to the numerical modelling of earthquake ground motion in structurally complex media and source dynamics.

Abstract This article provides an overview of the application of the staggered-grid finite-difference technique to model wave propagation problems in 3D elastic media. In addition to presenting generalized, discrete representations of the differential equations of motion using the staggered-grid approach, we also provide detailed formulations that describe the incorporation of moment-tensor sources, the implementation of a stable and accurate representation of a planar free-surface boundary for 3D models, and the derivation and implementation of an approximate technique to model spatially variable anelastic attenuation within time-domain finite-difference computations. The comparison of results obtained using the staggered-grid technique with those obtained using a frequency-wavenumber algorithm shows excellent agreement between the two methods for a variety of models. In addition, this article also introduces a memory optimization procedure that allows large-scale 3D finite-difference problems to be computed on a conventional, single-processor desktop workstation. With this technique, model storage is accommodated using both external (hard-disk) and internal (core) memory. To reduce system overhead, a cascaded time update procedure is utilized to maximize the number of computations performed between I/O operations. This formulation greatly expands the applicability of the 3D finite-difference technique by providing an efficient and practical algorithm for implementation on commonly available workstation platforms.

Abstract Evaluating the performance of finite-difference algorithms typically uses a technique known as von Neumann analysis. For a given algorithm, application of the technique yields both a dispersion relation valid for the discrete timespace grid and a mathematical condition for stability. In practice, a major shortcoming of conventional von Neumann analysis is that it can be applied only to an idealized numerical model — that of an infinite, homogeneous whole space. Experience has shown that numerical instabilities often arise in finite-difference simulations of wave propagation at interfaces address this issue, I generalize von Neumann analysis for a model of two half-spaces. I perform the analysis for the case of acoustic wave propagation using a standard staggered-grid finite-difference numerical scheme. By deriving expressions for the discrete reflection and transmission coefficients, I study under what conditions the discrete reflection and transmission coefficients become unbounded. I find that instabilities encountered in numerical modeling near interfaces with strong material contrasts are linked to these cases and develop a modified stability criterion that takes into account the resulting instabilities. I test and verify the stability criterion by executing a finite-difference algorithm under conditions predicted to be stable and unstable.