Modeling Of Seismic Wave Propagation Research Articles

Shallow structure of the Chalco and Xochimilco sub-basins (southern Mexico basin) based on wave propagation modelling and seismic data

The shallow structures of the Chalco and Xochimilco sub-basins (southeastern Mexico Basin) are established based on a seismic line that stretches trough both depressions. Interpretation of the seismic line was constraint by the stratigraphic column cut by the approximately 3000 m deep Tulyehualco well. The seismic line was converted to depth based on a vertical transverse isotropy velocity model determined from the interval, and root mean square (rms) velocities of the geologic units as established at the Tulyehualco deep well. P-wave modelling confirmed the results. The volcano-sedimentary infill established is consistent with the stratigraphic column of the 502 m deep MexiDrill borehole. It was possible to establish 1) the base of the fine stratified lacustrine sediments, and 2) the top of the basalts along the Xochimilco and Chalco sub-basins. The Xochimilco and Chalco sub-basins are separated by a shallow structural high (interpreted as due volcanic flows) against which the lacustrine sediments abut. A conspicuous anticline shaped structural high (adjacent to Xico tuff ring) divides the Chalco sub-basin into two depressions. The top of the basaltic basement matches with the upper volcanic unit horizon. In the western trough, the basement gets deeper from east to west (from 336 m at the structural high, to 660 m to the west). Immediately east of the anticline-shaped structural high, the basement is located at 510 m. Beneath the Xochimilco trough (at 180 m approximately) the volcanic basement shallows up smoothly northwestwards. In the eastern part of the Xochimilco sub-basin, the shallow lacustrian sediments are relatively thin. Northwestwards of the Tulyehualco well, their thickness diminishes. This apparent thinning is accounted by the fact that there the seismic line occupies the northern rim of the depression. Reflections in the upper 300 m of the section are not continuous. Beneath the central part of Chalco lake the reflectors are horizontal with a slight westward dip. At both sides of this central zone, reflectors are inclined or have dome shapes. This indicates that deposition along the Chalco sub-basin has been deformed. The thickness of the reflections associated with the fine laminated to volcanoclastic material transition zone is approximately constant, but thins and finally abuts against the structural high. In the eastern trough, seismic reflections at the shallow first 300 m are not as fine as in the western trough. These low magnitude seismic reflections might be associated with coarse fluvial volcanoclastic material from nearby Popocatéptl volcano. Seismic amplification in these two tectonic depressions by numerical modelling of seismic wave propagation can be based on the structure here established. Between 8 and 9 high amplitude reflectors in the first 300 m of the infill sequence are interpreted as due to impedance contrasts between geological units. Density and velocity logs along the MexiDrill borehole support this interpretation. Shallow strata are approximately 30 m thick. Towards the basin base, the strata are thicker (70 m approximately). These thicknesses are consistent with the nature of the volcano-sedimentary infill cut by the MexiDrill borehole. In particular, each seismic reflection would represent a sedimentation cycle of approximately 30,000 years (corresponding to a mean thickness of 30 m). If a sedimentation rate of 0.8 mm per year holds along the first 200–300 m of the basin infill (i.e., neglecting compaction effects, etc.) then this upper sequence comprises a record of paleo-ambiental changes that affected the southern Mexico Basin during the last 250,000 to 375,000 years.

Optimal Third-Order Symplectic Integration Modeling of Seismic Acoustic Wave Propagation

ABSTRACTSeismic wavefield modeling based on the wave equation is widely used in understanding and predicting the dynamic and kinematic characteristics of seismic wave propagation through media. This article presents an optimal numerical solution for the seismic acoustic wave equation in a Hamiltonian system based on the third-order symplectic integrator method. The least absolute truncation error analysis method is used to determine the optimal coefficients. The analysis of the third-order symplectic integrator shows that the proposed scheme exhibits high stability and minimal truncation error. To illustrate the accuracy of the algorithm, we compare the numerical solutions generated by the proposed method with the theoretical analysis solution for 2D and 3D seismic wave propagation tests. The results show that the proposed method reduced the phase error to the eighth-order magnitude accuracy relative to the exact solution. These simulations also demonstrated that the proposed third-order symplectic method can minimize numerical dispersion and preserve the waveforms during the simulation. In addition, comparing different central frequencies of the source and grid spaces (90, 60, and 20 m) for simulation of seismic wave propagation in 2D and 3D models using symplectic and nearly analytic discretization methods, we deduce that the suitable grid spaces are roughly equivalent to between one-fourth and one-fifth of the wavelength, which can provide a good compromise between accuracy and computational cost.

Optimizing the finite-difference implementation of three-dimensional free-surface boundary in frequency-domain modeling of elastic waves

The problem of modeling seismic wave propagation for multiple sources, such as in the solution of gradient-based elastic full-waveform inversion, is an important topic in seismic exploration. The frequency-domain finite-difference (FD) method is a good choice for this purpose, mainly because of its simple discretization and high computational efficiency. However, when it comes to modeling the complete elastic wavefields, this approach has limited surface-wave accuracy because, when modeling with the strong form of the wave equation, it is not always easy to implement an accurate stress-free boundary condition. Although a denser spatial sampling is helpful for overcoming this problem, the additional discrete points will significantly increase the computational cost in the resolution of its resulting discrete system, especially in 3D problems. Furthermore, sometimes, when modeling with optimized schemes, an inconsistency in the computation precision between the regions at the free surface and inside the model volume would happen and introduce numerical artifacts. To overcome these issues, we have considered optimizing the FD implementation of the free-surface boundary. In our method, the problem was formulated in terms of a novel system of partial differential equations satisfied at the free surface, and the weighted-averaging strategy was introduced to optimize its discretization. With this approach, we can impose FD schemes for the free surface and internal region consistently and improve their discretization precision simultaneously. Benchmark tests for Lamb’s problem indicate that the proposed free-surface implementation contributes to improving the simulation accuracy on surface waves, without increasing the number of grid points per wavelength. This reveals the potential of developing optimized schemes in the free-surface implementation. In particular, through the successful introduction of weighting coefficients, this free-surface FD implementation enables adaptation to the variation of Poisson’s ratio, which is very useful for modeling in heterogeneous near-surface weathered zones.

Multilevel Monte Carlo acceleration of seismic wave propagation under uncertainty

We interpret uncertainty in a model for seismic wave propagation by treating the model parameters as random variables, and apply the Multilevel Monte Carlo method to reduce the cost of approximating expected values of selected, physically relevant, quantities of interest (QoI) with respect to the random variables. Targeting source inversion problems, where the source of an earthquake is inferred from ground motion recordings on the Earth’s surface, we consider two QoIs that measure the discrepancies between computed seismic signals and given reference signals: one QoI, $${\mathscr {Q}}_E$$, is defined in terms of the $$L^2$$-misfit, which is directly related to maximum likelihood estimates of the source parameters; the other, $${\mathscr {Q}}_W$$, is based on the quadratic Wasserstein distance between probability distributions, and represents one possible choice in a class of such misfit functions that have become increasingly popular to solve seismic inversion in recent years. We simulate seismic wave propagation, including seismic attenuation, using a publicly available code in widespread use, based on the spectral element method. Using random coefficients and deterministic initial and boundary data, we present benchmark numerical experiments with synthetic data in a two-dimensional physical domain and a one-dimensional velocity model where the assumed parameter uncertainty is motivated by realistic Earth models. Here, the computational cost of the standard Monte Carlo method was reduced by up to $$97\%$$ for $${\mathscr {Q}}_E$$, and up to $$78\%$$ for $${\mathscr {Q}}_W$$, using a relevant range of tolerances. Shifting to three-dimensional domains is straight-forward and will further increase the relative computational work reduction.

A compact high order Alternating Direction Implicit method for three-dimensional acoustic wave equation with variable coefficient

Efficient and accurate numerical simulation of seismic wave propagation is important in various Geophysical applications such as seismic full waveform inversion (FWI) problem. However, due to the large size of the physical domain and requirement on low numerical dispersion, many existing numerical methods are inefficient for numerical modelling of seismic wave propagation in an heterogeneous media. Despite the great efforts that have been devoted during the past decades, it still remains a challenging task in the development of efficient and accurate finite difference method for multi-dimensional acoustic wave equation with variable velocity. In this paper we proposed a Padé approximation based finite difference scheme for solving the acoustic wave equation in three-dimensional heterogeneous media. The new method is obtained by combining the Padé approximation and a novel algebraic manipulation. The efficiency of the new algorithm is further improved through the Alternative Directional Implicit (ADI) method. The stability of the new algorithm has been theoretically proved by energy method. The new method is conditionally stable with a better Courant–Friedrichs–Lewy condition (CFL) condition, which has been verified numerically. Extensive numerical examples have been solved, which demonstrated that the new method is accurate, efficient and stable.

A discrete representation of a heterogeneous viscoelastic medium for the finite-difference modelling of seismic wave propagation

A discrete representation of a heterogeneous viscoelastic medium for the finite-difference modelling of seismic wave propagation

Numerical simulation of the seismic wave propagation and fluid pressure in complex porous media at the mesoscopic scale

In this study, we apply an improved high-order staggered-grid finite difference method to model seismic wave propagation in complex heterogeneous porous media with cracks. Upon changing the mesh step and applying a broadband wavelet, we obtain numerical simulation results at the mesoscopic scale. As a fundamental result of applying Biot’s theory, three types of waves are also found. We discuss the relationship between the shear wave splitting parameter and the crack density of an equivalent medium in detail. Based on Darcy’s law, the finite element method is used to investigate the distribution of the fluid pressure field in porous media at the mesoscopic scale. A flat ellipse is used as an approximation of a crack. According to the simulation results, we further study the influence of variations in the crack angle on the distribution of the fluid pressure. Moreover, the formulas derived from the Fourier transform are used to calculate the inverse attenuation quality factors and phase velocities of the fast compressional wave, slow compressional wave and shear wave. The effects of the porosity and elastic modulus of an equivalent medium on wave attenuation are also taken into account in this study.

2-D P-SV and SH spectral element modelling of seismic wave propagation in non-linear media with pore-pressure effects

SUMMARYIt has long been recognized that the effects of superficial geological layers, or site effects, can play a major role on the seismic ground motion at the free surface. In this study, we compute wave propagation in a 2-D asymmetrical basin considering both soil non-linearity and pore-pressure effects. Equations of elastodynamics of wave propagation are solved using the spectral element method (SEM). The geometry of the basin gives rise to basin-edge generated waves, that are different for in-plane (P-SV) and out-of-plane (SH) wave propagation and resulting in different non-linear response. Moreover, the excess-pore pressure development in superficial liquefiable layers (effective stress analysis) brings larger deformation and loss of strength than the analysis without pore-pressure effects (total stress analysis). The coupling of vertically propagating waves and the waves specifically generated in 1-D model leads to waves whose amplitude and duration are higher than the 1-D case. This multidimensional effect increases material non-linearity. Such complex wavefield provokes larger deformation and higher pore-pressure rise that cannot be predicted by 1-D modelling. Therefore, our paper suggests the use of multidimensional modelling while studying seismic wave propagation in both linear and non-linear complex media.

Pre-mission InSights on the Interior of Mars

Abstract The Interior exploration using Seismic Investigations, Geodesy, and Heat Trans- port (InSight) Mission will focus on Mars’ interior structure and evolution. The basic structure of crust, mantle, and core form soon after accretion. Understanding the early differentiation process on Mars and how it relates to bulk composition is key to improving our understanding of this process on rocky bodies in our solar system, as well as in other solar systems. Current knowledge of differentiation derives largely from the layers observed via seismology on the Moon. However, the Moon’s much smaller diameter make it a poor analog with respect to interior pressure and phase changes. In this paper we review the current knowledge of the thickness of the crust, the diameter and state of the core, seismic attenuation, heat flow, and interior composition. InSight will conduct the first seismic and heat flow measurements of Mars, as well as more precise geodesy. These data reduce uncertainty in crustal thickness, core size and state, heat flow, seismic activity and meteorite impact rates by a factor of 3–10× relative to previous estimates. Based on modeling of seismic wave propagation, we can further constrain interior temperature, composition, and the location of phase changes. By combining heat flow and a well constrained value of crustal thickness, we can estimate the distribution of heat producing elements between the crust and mantle. All of these quantities are key inputs to models of interior convection and thermal evolution that predict the processes that control subsurface temperature, rates of volcanism, plume distribution and stability, and convective state. Collectively these factors offer strong controls on the overall evolution of the geology and habitability of Mars.

Seismic wave propagation in nonlinear viscoelastic media using the auxiliary differential equation method

In previous studies, the auxiliary differential equation (ADE) method has been applied to 2-D seismic-wave propagation modelling in viscoelastic media. This method is based on the separation of the wave propagation equations derived from the constitutive law defining the stress–strain relation. We make here a 3-D extension of a finite-difference (FD) scheme to solve a system of separated equations consisting in the stress–strain rheological relation, the strain–velocity and the velocity–stress equations. The current 3-D FD scheme consists in the discretization of the second order formulation of a non-linear viscoelastic wave equation with a time actualization of the velocity and displacement fields. Compared to the usual memory variable formalism, the ADE method allows flexible implementation of complex expressions of the desired rheological model such as attenuation/viscoelastic models or even non-linear behaviours, with physical parameters that can be provided from dispersion analysis. The method can also be associated with optimized perfectly matched layers-based boundary conditions that can be seen as additional attenuation (viscoelastic) terms. We present the results obtained for a non-linear viscoelastic model made of a Zener viscoelastic body associated with a non-linear quadratic strain term. Such non-linearity is relevant to define unconsolidated granular model behaviour. Thanks to a simple model, but without loss of generality, we demonstrate the accuracy of the proposed numerical approach.

A discrete representation of material heterogeneity for the finite-difference modelling of seismic wave propagation in a poroelastic medium

A discrete representation of material heterogeneity for the finite-difference modelling of seismic wave propagation in a poroelastic medium

Enriched Galerkin finite element approximation for elastic wave propagation in fractured media

Finite Element Methods (FEM) are becoming increasingly popular in modeling seismic wave propagation. These methods provide higher order accuracy, geometrical flexibility and adaptive gridding capabilities that are not easy to incorporate in traditional finite difference methods employed for generation of synthetic seismograms. Moreover, several studies have shown that Discontinuous Galerkin FEM (DGM) is a promising approach for modeling wave propagation in fractured media. Here we propose an Enriched Galerkin FEM (EGM) for elastic wave propagation. EGM uses the same bilinear form as DGM and the continuous Galerkin finite element spaces enriched by discontinuous piecewise constants or bilinear functions. EGM satisfies local equilibrium while reducing the degrees of freedom in DGM formulations. In this paper, we consider elastic wave propagation and derive optimal a priori error estimates for DGM and EGM. We present numerical examples in two spatial dimensions that confirm these theoretical results. In addition, we provide numerical comparisons with the Spectral element method. In previous work, DGM has been shown to be effective in modeling elastic wave propagation in fractured media using the linear slip model. We now extend these results to EGM with reduced computational costs over DGM.

Numerical modeling of Gaussian beam propagation and diffraction in inhomogeneous media based on the complex eikonal equation

Gaussian beam is an important complex geometrical optical technology for modeling seismic wave propagation and diffraction in the subsurface with complex geological structure. Current methods for Gaussian beam modeling rely on the dynamic ray tracing and the evanescent wave tracking. However, the dynamic ray tracing method is based on the paraxial ray approximation and the evanescent wave tracking method cannot describe strongly evanescent fields. This leads to inaccuracy of the computed wave fields in the region with a strong inhomogeneous medium. To address this problem, we compute Gaussian beam wave fields using the complex phase by directly solving the complex eikonal equation. In this method, the fast marching method, which is widely used for phase calculation, is combined with Gauss–Newton optimization algorithm to obtain the complex phase at the regular grid points. The main theoretical challenge in combination of this method with Gaussian beam modeling is to address the irregular boundary near the curved central ray. To cope with this challenge, we present the non-uniform finite difference operator and a modified fast marching method. The numerical results confirm the proposed approach.

An incrementation-adaptive multi-transmitting boundary for seismic fracture analysis of concrete gravity dams

This study presents a numerical framework and methodology that features a nonlinear extended finite element formulation for earthquake fracture analysis of dam-reservoir-sediment-foundation coupled systems and an adaptive time incrementation multi-transmitting boundary for modeling seismic wave propagation in both infinite reservoir and foundation domains. To improve the convergence rate during the iterative solution of the highly nonlinear dynamic equilibrium equation, the original Multi-Transmitting Formula (MTF) is modified and coupled to the automatic incrementation time-stepping scheme through the introduction of a fictitious time scale in conjunction with an interpolation scheme that links the fictitious time scale to the real time scale. Meanwhile, the dynamic crack initiation and propagation in the body of the dam is represented by the extended finite element method including the effects of branching and intersecting cracks. As an application of the proposed formulation, the seismic cracking behavior of a representative concrete gravity dam system is analyzed.

Numerical Modeling of 3D Seismic Wave Propagation around Yogyakarta, the Southern Part of Central Java, Indonesia, Using Spectral-Element Method on MPI-GPU Cluster

A strong tectonic earthquake with a magnitude of 5.9 Richter scale has been occurred in Yogyakarta and Central Java on May 26, 2006. The earthquake has caused severe damage in Yogyakarta and the southern part of Central Java, Indonesia. The understanding of seismic response of earthquake among ground shaking and the level of building damage is important. We present numerical modeling of 3D seismic wave propagation around Yogyakarta and the southern part of Central Java using spectral-element method on MPI-GPU (Graphics Processing Unit) computer cluster to observe its seismic response due to the earthquake. The homogeneous 3D realistic model is generated with detailed topography surface. The influences of free surface topography and layer discontinuity of the 3D model among the seismic response are observed. The seismic wave field is discretized using spectral-element method. The spectral-element method is solved on a mesh of hexahedral elements that is adapted to the free surface topography and the internal discontinuity of the model. To increase the data processing capabilities, the simulation is performed on a GPU cluster with implementation of MPI (Message Passing Interface).

An adaptive free‐surface expression for three‐dimensional finite‐difference frequency‐domain modelling of elastic wave

ABSTRACTFinite‐difference frequency‐domain modelling of seismic wave propagation is attractive for its efficient solution of multisource problems, and this is crucial for full‐waveform inversion and seismic imaging, especially in the three‐dimensional seismic problem. However, implementing the free surface in the finite‐difference method is nontrivial. Based on an average medium method and the limit theorem, we present an adaptive free‐surface expression to describe the behaviour of wavefields at the free surface, and no extra work for the free‐surface boundary condition is needed. Essentially, the proposed free‐surface expression is a modification of density and constitutive relation at the free surface. In comparison with a direct difference approximate method of the free‐surface boundary condition, this adaptive free‐surface expression can produce more accurate and stable results for a broad range of Poisson's ratio. In addition, this expression has a good performance in handling the lateral variation of Poisson's ratio adaptively and without instability.

Modeling of frequency-domain scalar wave equation with the average-derivative optimal scheme based on a multigrid-preconditioned iterative solver

An efficient finite-difference frequency-domain modeling of seismic wave propagation relies on the discrete schemes and appropriate solving methods. The average-derivative optimal scheme for the scalar wave modeling is advantageous in terms of the storage saving for the system of linear equations and the flexibility for arbitrary directional sampling intervals. However, using a LU-decomposition-based direct solver to solve its resulting system of linear equations is very costly for both memory and computational requirements. To address this issue, we consider establishing a multigrid-preconditioned BI-CGSTAB iterative solver fit for the average-derivative optimal scheme. The choice of preconditioning matrix and its corresponding multigrid components is made with the help of Fourier spectral analysis and local mode analysis, respectively, which is important for the convergence. Furthermore, we find that for the computation with unequal directional sampling interval, the anisotropic smoothing in the multigrid precondition may affect the convergence rate of this iterative solver. Successful numerical applications of this iterative solver for the homogenous and heterogeneous models in 2D and 3D are presented where the significant reduction of computer memory and the improvement of computational efficiency are demonstrated by comparison with the direct solver. In the numerical experiments, we also show that the unequal directional sampling interval will weaken the advantage of this multigrid-preconditioned iterative solver in the computing speed or, even worse, could reduce its accuracy in some cases, which implies the need for a reasonable control of directional sampling interval in the discretization.

Spectral element modelling of seismic wave propagation in visco-elastoplastic media including excess-pore pressure development

Numerical modelling of seismic wave propagation, considering soil nonlinearity, has become a major topic in seismic hazard studies when strong shaking is involved under particular soil conditions. Indeed, when strong ground motion propagates in saturated soils, pore pressure is another important parameter to take into account when successive phases of contractive and dilatant soil behaviour are expected. Here, we model 1-D seismic wave propagation in linear and nonlinear media using the spectral element numerical method. The study uses a three-component (3C) nonlinear rheology and includes pore-pressure excess. The 1-D– 3C model is used to study the 1987 Superstition Hills earthquake (ML 6.6), which was recorded at the Wildlife Refuge Liquefaction Array, USA. The data of this event present strong soil nonlinearity involving pore-pressure effects. The ground motion is numerically modelled for different assumptions on soil rheology and input motion (1C versus 3C), using the recorded borehole signals as input motion. The computed acceleration–time histories show low-frequency amplification and strong high-frequency damping due to the development of pore pressure in one of the soil layers. Furthermore, the soil is found to be more nonlinear and more dilatant under triaxial loading compared to the classical 1C analysis, and significant differences in surface displacements are observed between the 1C and 3C approaches. This study contributes to identify and understand the dominant phenomena occurring in superficial layers, depending on local soil properties and input motions, conditions relevant for site- specific studies.

Numerical modeling of seismic wave propagation and ground motion in underground mines

Dynamic rock support plays a critical role in ensuring mine safety and its design depends on several key design inputs. For rock support design and forensic analysis of rockburst damage, it is important to understand seismic wave propagation and ground motion distribution near excavation boundaries for seismic waves generated from a remote fault slip seismic event. This study aims at achieving a better understanding of wave propagation patterns around mine tunnels and capturing ground motion (Peak Particle Velocity - PPV) accurately for dynamic support design and forensic analysis. An advanced seismic wave propagation modeling tool, SPECFEM2D, is used to study complex wave propagation in underground mines. Attentions are paid to studying the influence of different mine excavations and geological structures on wavefields and ground motions near excavation boundaries. The simulation results show that the wavefields and the ground motion distributions become more complicated as more mine excavations and geological structures are involved. Moreover, the PPV distribution around a tunnel can be altered, leading to high and low PPV zones around the tunnel. In addition, modulation of travel time and long S-coda waves can be observed in the complex seismograms. Using the modeling approach, areas in an underground mine that may experience high potentials of rockburst damage could be identified and mine safety could be improved by implementing dynamic rock ground support in these areas.

Immersive boundary conditions: Theory, implementation, and examples

Many applications in computational geophysics involve the modeling of seismic wave propagation on a set of closely related subsurface models. In such scenarios, it is of interest to recompute the seismic wavefields locally (only in the regions of change), instead of in the full subsurface model. We have developed a method for local acoustic wavefield recomputation that makes it possible to fully immerse a local modeling domain within a larger domain of arbitrary extent and complexity, such that the wave propagation in the full domain is completely accounted for. The method enables wavefield modeling on much smaller local domains, while relying on the up-front generation of a large number of Green’s functions and a wavefield extrapolation step at each time step of the simulation. A Kirchhoff-Helmholtz extrapolation integral is used to predict the interaction of the wavefield leaving the local domain with the exterior domain. The outward propagating wavefield and the wavefield reentering the local domain are applied as a boundary condition along the edges. Thanks to these dynamically calculated boundary conditions, all higher order long-range interactions between the two domains are fully accounted for. We have implemented the method in a conventional finite-difference time-domain scheme and determined that the locally calculated wavefields are equal to wavefields generated on the full domain to within numerical precision. The efficiency of the local modeling algorithm will greatly depend on the nature and size of the problem.