Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves

A method is proposed for accurately describing arbitrary-shaped free boundaries in single-grid finite-difference schemes for elastodynamics, in a time-domain velocity-stress framework. The basic idea is as follows: fictitious values of the solution are built in vacuum, and injected into the numerical integration scheme near boundaries. The most original feature of this method is the way in which these fictitious values are calculated. They are based on boundary conditions and compatibility conditions satisfied by the successive spatial derivatives of the solution, up to a given order that depends on the spatial accuracy of the integration scheme adopted. Since the work is mostly done during the preprocessing step, the extra computational cost is negligible. Stress-free conditions can be designed at any arbitrary order without any numerical instability, as numerically checked. Using 10 grid nodes per minimal S-wavelength with a propagation distance of 50 wavelengths yields highly accurate results. With 5 grid nodes per minimal S-wavelength, the solution is less accurate but still acceptable. A subcell resolution of the boundary inside the Cartesian meshing is obtained, and the spurious diffractions induced by staircase descriptions of boundaries are avoided. Contrary to what occurs with the vacuum method, the quality of the numerical solution obtained with this method is almost independent of the angle between the free boundary and the Cartesian meshing.Comment: accepted and to be published in Geophys. J. In

Generalized multiscale finite element methods for wave propagation in heterogeneous media

Numerical modeling of wave propagation in heterogeneous media is important in many applications. Due to the complex nature, direct numerical simulations on the fine grid are prohibitively expensive. It is therefore important to develop efficient and accurate methods that allow the use of coarse grids. In this paper, we present a multiscale finite element method for wave propagation on a coarse grid. The proposed method is based on the Generalized Multiscale Finite Element Method (GMsFEM). To construct multiscale basis functions, we start with two snapshot spaces in each coarse-grid block where one represents the degrees of freedom on the boundary and the other represents the degrees of freedom in the interior. We use local spectral problems to identify important modes in each snapshot space. These local spectral problems are different from each other and their formulations are based on the analysis. To our best knowledge, this is the first time where multiple snapshot spaces and multiple spectral problems are used and necessary for efficient computations. Using the dominant modes from local spectral problems, multiscale basis functions are constructed to represent the solution space locally within each coarse block. These multiscale basis functions are coupled via the symmetric interior penalty discontinuous Galerkin method which provides a block diagonal mass matrix, and, consequently, results in fast computations in an explicit time discretiza- tion. Our methods' stability and spectral convergence are rigorously analyzed. Numerical examples are presented to show our methods' performance. We also test oversampling strategies. In particular, we discuss how the modes from different snapshot spaces can affect the proposed methods' accuracy

Modeling seismic wave propagation and amplification in 1D/2D/3D linear and nonlinear unbounded media

To analyze seismic wave propagation in geological structures, it is possible to consider various numerical approaches: the finite difference method, the spectral element method, the boundary element method, the finite element method, the finite volume method, etc. All these methods have various advantages and drawbacks. The amplification of seismic waves in surface soil layers is mainly due to the velocity contrast between these layers and, possibly, to topographic effects around crests and hills. The influence of the geometry of alluvial basins on the amplification process is also know to be large. Nevertheless, strong heterogeneities and complex geometries are not easy to take into account with all numerical methods. 2D/3D models are needed in many situations and the efficiency/accuracy of the numerical methods in such cases is in question. Furthermore, the radiation conditions at infinity are not easy to handle with finite differences or finite/spectral elements whereas it is explicitely accounted in the Boundary Element Method. Various absorbing layer methods (e.g. F-PML, M-PML) were recently proposed to attenuate the spurious wave reflections especially in some difficult cases such as shallow numerical models or grazing incidences. Finally, strong earthquakes involve nonlinear effects in surficial soil layers. To model strong ground motion, it is thus necessary to consider the nonlinear dynamic behaviour of soils and simultaneously investigate seismic wave propagation in complex 2D/3D geological structures! Recent advances in numerical formulations and constitutive models in such complex situations are presented and discussed in this paper. A crucial issue is the availability of the field/laboratory data to feed and validate such models.Comment: of International Journal Geomechanics (2010) 1-1

Solving seismic wave propagation in elastic media using the matrix exponential approach

Three numerical algorithms are proposed to solve the time-dependent elastodynamic equations in elastic solids. All algorithms are based on approximating the solution of the equations, which can be written as a matrix exponential. By approximating the matrix exponential with a product formula, an unconditionally stable algorithm is derived that conserves the total elastic energy density. By expanding the matrix exponential in Chebyshev polynomials for a specific time instance, a so-called ``one-step'' algorithm is constructed that is very accurate with respect to the time integration. By formulating the conventional velocity-stress finite-difference time-domain algorithm (VS-FDTD) in matrix exponential form, the staggered-in-time nature can be removed by a small modification, and higher order in time algorithms can be easily derived. For two different seismic events the accuracy of the algorithms is studied and compared with the result obtained by using the conventional VS-FDTD algorithm.Comment: 13 pages revtex, 6 figures, 2 table

Three-dimensional dynamic rupture simulation with a high-order discontinuous Galerkin method on unstructured tetrahedral meshes

Accurate and efficient numerical methods to simulate dynamic earthquake rupture and wave propagation in complex media and complex fault geometries are needed to address fundamental questions in earthquake dynamics, to integrate seismic and geodetic data into emerging approaches for dynamic source inversion, and to generate realistic physics-based earthquake scenarios for hazard assessment. Modeling of spontaneous earthquake rupture and seismic wave propagation by a high-order discontinuous Galerkin (DG) method combined with an arbitrarily high-order derivatives (ADER) time integration method was introduced in two dimensions by de la Puente et al. (2009). The ADER-DG method enables high accuracy in space and time and discretization by unstructured meshes. Here we extend this method to three-dimensional dynamic rupture problems. The high geometrical flexibility provided by the usage of tetrahedral elements and the lack of spurious mesh reflections in the ADER-DG method allows the refinement of the mesh close to the fault to model the rupture dynamics adequately while concentrating computational resources only where needed. Moreover, ADER-DG does not generate spurious high-frequency perturbations on the fault and hence does not require artificial Kelvin-Voigt damping. We verify our three-dimensional implementation by comparing results of the SCEC TPV3 test problem with two well-established numerical methods, finite differences, and spectral boundary integral. Furthermore, a convergence study is presented to demonstrate the systematic consistency of the method. To illustrate the capabilities of the high-order accurate ADER-DG scheme on unstructured meshes, we simulate an earthquake scenario, inspired by the 1992 Landers earthquake, that includes curved faults, fault branches, and surface topography

Energy-conserving 3D elastic wave simulation with finite difference discretization on staggered grids with nonconforming interfaces

In this work, we describe an approach to stably simulate the 3D isotropic elastic wave propagation using finite difference discretization on staggered grids with nonconforming interfaces. Specifically, we consider simulation domains composed of layers of uniform grids with different grid spacings, separated by planar interfaces. This discretization setting is motivated by the observation that wave speeds of earth media tend to increase with depth due to sedimentation and consolidation processes. We demonstrate that the layer-wise finite difference discretization approach has the potential to significantly reduce the simulation cost, compared to its counterpart that uses holistically uniform grids. Such discretizations are enabled by summation-by-parts finite difference operators, which are standard finite difference operators with special adaptations near boundaries or interfaces, and simultaneous approximation terms, which are penalty terms appended to the discretized system to weakly impose boundary or interface conditions. Combined with specially designed interpolation operators, the discretized system is shown to preserve the energy-conserving property of the continuous elastic wave equation, and a fortiori ensure the stability of the simulation. Numerical examples are presented to corroborate these analytical developments

Stability analysis of second- and fourth-order finite-difference modelling of wave propagation in orthotropic media

The stability of the finite-difference approximation of elastic wave propagation in orthotropic homogeneous media in the three-dimensional case is discussed. The model applies second- and fourth-order finite-difference approaches with staggered grid and stress-free boundary conditions in the space domain and second-order finite-difference approach in the time domain. The numerical integration of the wave equation by central differences is conditionally stable and the corresponding stability criterion for the time domain discretisation has been deduced as a function of the material properties and the geometrical discretization. The problem is discussed by applying the method of VonNeumann. Solutions and the calculation of the critical time steps is presented for orthotropic material in both the second- and fourth-order case. The criterion is verified for the special case of isotropy and results in the well-known formula from the literature. In the case of orthotropy the method was verified by long time simulations and by calculating the total energy of the system

Modelling Seismic Wave Propagation for Geophysical Imaging

International audienceThe Earth is an heterogeneous complex media from the mineral composition scale (10−6m) to the global scale ( 106m). The reconstruction of its structure is a quite challenging problem because sampling methodologies are mainly indirect as potential methods (Günther et al., 2006; Rücker et al., 2006), diffusive methods (Cognon, 1971; Druskin & Knizhnerman, 1988; Goldman & Stover, 1983; Hohmann, 1988; Kuo & Cho, 1980; Oristaglio & Hohmann, 1984) or propagation methods (Alterman & Karal, 1968; Bolt & Smith, 1976; Dablain, 1986; Kelly et al., 1976; Levander, 1988; Marfurt, 1984; Virieux, 1986). Seismic waves belong to the last category. We shall concentrate in this chapter on the forward problem which will be at the heart of any inverse problem for imaging the Earth. The forward problem is dedicated to the estimation of seismic wavefields when one knows the medium properties while the inverse problem is devoted to the estimation of medium properties from recorded seismic wavefields