Modeling Of Seismic Wave Propagation Research Articles

Scaling behavior and the effects of heterogeneity on shallow seismic imaging of mineral deposits: A case study from Brunswick No. 6 mining area, Canada

Abstract We have studied the scaling behavior of compressional-wave velocity and density logs from an exploration borehole that extends down to about 700 m depth in the Brunswick No. 6 mining area, Bathurst Mining Camp, Canada. Using statistical methods, vertical and horizontal scale lengths of heterogeneity were estimated. Vertical scale length estimates from the velocity, density and calculated acoustic impedance are 14 m, 33 m, and about 20 m, respectively. Although the estimated scale length for the acoustic impedance implies a weak scattering environment, elastic finite difference modeling of seismic wave propagation in 2D heterogeneous media demonstrates that even this weak scattering medium can mask seismic signals from small, but yet economically feasible, massive sulfide deposits. Further analysis of the synthetic seismic data suggests that in the presence of heterogeneity, lenticular-shaped targets may only exhibit incomplete diffraction signals whereby the down-dip tails of these diffractions are mainly visible on the stacked sections. Therefore, identification of orebody generated diffractions is much easier on the unmigrated stacked sections than on migrated stacked sections. The numerical seismic modeling in 2D heterogeneous media indicates that in the presence of large horizontal, but small vertical scale lengths (structural anisotropy), identification of massive sulfide deposits is possible, but their delineation at depth requires detailed velocity modeling and processing algorithms which can handle the anisotropy.

Energy efficiency vs. performance of the numerical solution of PDEs: An application study on a low-power ARM-based cluster

Power consumption and energy efficiency are becoming critical aspects in the design and operation of large scale HPC facilities, and it is unanimously recognised that future exascale supercomputers will be strongly constrained by their power requirements. At current electricity costs, operating an HPC system over its lifetime can already be on par with the initial deployment cost. These power consumption constraints, and the benefits a more energy-efficient HPC platform may have on other societal areas, have motivated the HPC research community to investigate the use of energy-efficient technologies originally developed for the embedded and especially mobile markets. However, lower power does not always mean lower energy consumption, since execution time often also increases. In order to achieve competitive performance, applications then need to efficiently exploit a larger number of processors. In this article, we discuss how applications can efficiently exploit this new class of low-power architectures to achieve competitive performance. We evaluate if they can benefit from the increased energy efficiency that the architecture is supposed to achieve. The applications that we consider cover three different classes of numerical solution methods for partial differential equations, namely a low-order finite element multigrid solver for huge sparse linear systems of equations, a Lattice-Boltzmann code for fluid simulation, and a high-order spectral element method for acoustic or seismic wave propagation modelling. We evaluate weak and strong scalability on a cluster of 96 ARM Cortex-A9 dual-core processors and demonstrate that the ARM-based cluster can be more efficient in terms of energy to solution when executing the three applications compared to an x86-based reference machine.

Elastic modeling of seismic wave propagation in partially saturated rocks

SUMMARY The effect of single-phase fluid saturation on a rock’s bulk modulus is well understood using Gassmann’s equation. However when multiple fluids are involved the behaviour is not as well understood. Several fluid mixing averages have been suggested (Voigt, Reuss, Hill), and each apply in certain situations, however it is often not clear which model to select in a specific scenario and in some scenarios none of the models are accurate. The critical factor in deciding which average to use depends on the way the fluids are spatially distributed within the rock. We have applied elastic finite difference computational modelling to many different fluid distribution scenarios and have replicated behaviour described by various theoretical, empirical and lab data results, as well as generating results that span the space between these models. Importantly, our results compare well with observations in lab experiments, without relying on poroelastic or squirt-flow models which require parameters that are difficult to estimate for real reservoirs. Our elastic scattering approach is less computationally expensive than poroelastic modelling and can be more easily applied to actual reservoir rock and fluid distributions. Our results provide us with a powerful new method to analyse and predict the effects of multiple fluids and ‘patchy’ saturation on saturated rock bulk moduli and velocity. They also challenge traditional assumptions about the controlling factors on saturated bulk moduli suggesting it is more dominantly affected by the spatial fluid distribution properties rather than pore-scale fluid flow effects.

Australian Seismological Reference Model (AuSREM): crustal component

Although Australia has been the subject of a wide range of seismological studies, these have concentrated on specific features of the continent at crustal scales and on the broad scale features in the mantle. The Australian Seismological Reference Model (AuSREM) is designed to bring together the existing information, and provide a synthesis in the form of a 3-D model that can provide the basis for future refinement from more detailed studies. Extensive studies in the last few decades provide good coverage for much of the continent, and the crustal model builds on the various data sources to produce a representative model that captures the major features of the continental structure and provides a basis for a broad range of further studies. The model is grid based with a 0.5° sampling in latitude and longitude, and is designed to be fully interpolable, so that properties can be extracted at any point. The crustal structure is built from five-layer representations of refraction and receiver function studies and tomographic information. The AuSREM crustal model is available at 1 km intervals. The crustal component makes use of prior compilations of sediment thicknesses, with cross checks against recent reflection profiling, and provides P and S wavespeed distributions through the crust. The primary information for P wavespeed comes from refraction profiles, for S wavespeed from receiver function studies. We are also able to use the results of ambient noise tomography to link the point observations into national coverage. Density values are derived using results from gravity interpretations with an empirical relation between P wavespeed and density. AuSREM is able to build on a new map of depth to Moho, which has been created using all available information including Moho picks from over 12 000 km of full crustal profiling across the continent. The crustal component of AuSREM provides a representative model that should be useful for modelling of seismic wave propagation and calculation of crustal corrections for tomography. Other applications include gravity studies and dynamic topography at the continental scale.

Numerical modelling of the upper-mantle anisotropy beneath a migrating strike-slip plate boundary: the San Andreas Fault system

We performed forward modelling of seismic anisotropy beneath a migrating strike-slip plate boundary to: (1) test if such a geodynamic context might explain teleseismic shear wave splitting data in the vicinity of the central part of the San Andreas Fault and (2) constrain the power of such data to unravel vertical and lateral variations in deformation patterns in the upper mantle. The modelling involves five steps: (1) thermomechanical modelling, using a finite-element code, of the deformation field, (2) viscoplastic self-consistent modelling of the resulting olivine and pyroxene crystal preferred orientations, (3) calculation of the elasticity tensors for different domains of the finite elements (FEs) model, (4) forward modelling of seismic wave propagation through the model using ray theory, finite-frequency theory and a full-wave approach and (5) performing splitting measurements on the synthetic seismograms. SKS splitting data in central California are best fitted by a model with a hotter geotherm within 60 km of the plate boundary accounting for the opening of an asthenospheric window due to the northward migration of the Mendocino Triple Junction. The westward motion of the plate boundary cannot however explain the rotation of fast polarization directions east of the San Andreas Fault in central California. Comparison between modelled and measured individual shear wave splitting also implies that the homogeneity of the two-layer models accounting for the observations in the vicinity of the San Andreas Fault indicates a sharp transition between lithospheric and asthenospheric deformations beneath this plate boundary. The ability of different synthetic approaches to localize horizontally and vertically the plate boundary-related deformation differs significantly. Splitting in ray theory synthetics closely follow variations in olivine crystallographic preferred orientations in the model. In contrast, splitting analysis on full-wave synthetics, which should be more representative of actual long-period SKS waves, results in smooth apparent lateral variations of anisotropy; the exact location and width of the plate boundary may only be retrieved by comparing fast-polarization profiles obtained using a multichannel analysis on waves with different periods.

3D finite-difference modeling of elastic wave propagation in the Laplace-Fourier domain

With the recent interest in the Laplace-Fourier domain full waveform inversion, we have developed new heterogeneous 3D fourth- and second-order staggered-grid finite-difference schemes for modeling seismic wave propagation in the Laplace-Fourier domain. Our approach is based on the integro-interpolation technique for the velocity-stress formulation in the Cartesian coordinate system. Five averaging elastic coefficients and three averaging densities are necessary to describe the heterogeneous medium, with harmonic averaging of the bulk and shear moduli, and arithmetic averaging of density. In the fourth-order approximation, we improved the accuracy of the scheme using a combination of integral identities for two elementary volumes — “small” and “large” around spatial gridpoints where the wave variables are defined. Two solution approaches are provided, both of which are solved with iterative Krylov methods. In the first approach, the stress variables are eliminated and a linear system for the velocity components is solved. In the second approach, we worked directly with the first-order system of velocity and stress variables. This reduced the computer memory required to store the complex matrix, along with reducing (by 30%) the number of arithmetic operations needed for the solution compared to the fourth-order scheme for velocity only. Numerical examples show that our finite-difference formulations for elastic wavefield simulations can achieve more accurate solutions with fewer grid points than those from previously published second and fourth-order frequency-domain schemes. We applied our simulator to the investigation of wavefields from the SEG/EAGE model in the Laplace-Fourier domain. The calculation is sensitive to the heterogeneity of the medium and capable of describing the structures of complex objects. Our technique can also be extended to 3D elastic modeling within the time domain.

Accurate and efficient modeling of global seismic wave propagation for an attenuative Earth model including the center

We propose a method for modeling global seismic wave propagation through an attenuative Earth model including the center. This method enables accurate and efficient computations since it is based on the 2.5-D approach, which solves wave equations only on a 2-D cross section of the whole Earth and can correctly model 3-D geometrical spreading. We extend a numerical scheme for the elastic waves in spherical coordinates using the finite-difference method (FDM), to solve the viscoelastodynamic equation. For computation of realistic seismic wave propagation, incorporation of anelastic attenuation is crucial. Since the nature of Earth material is both elastic solid and viscous fluid, we should solve stress–strain relations of viscoelastic material, including attenuative structures. These relations represent the stress as a convolution integral in time, which has had difficulty treating viscoelasticity in time-domain computation such as the FDM. However, we now have a method using so-called memory variables, invented in the 1980s, followed by improvements in Cartesian coordinates. Arbitrary values of the quality factor (Q) can be incorporated into the wave equation via an array of Zener bodies. We also introduce the multi-domain, an FD grid of several layers with different grid spacings, into our FDM scheme. This allows wider lateral grid spacings with depth, so as not to perturb the FD stability criterion around the Earth center. In addition, we propose a technique to avoid the singularity problem of the wave equation in spherical coordinates at the Earth center. We develop a scheme to calculate wavefield variables on this point, based on linear interpolation for the velocity–stress, staggered-grid FDM. This scheme is validated through a comparison of synthetic seismograms with those obtained by the Direct Solution Method for a spherically symmetric Earth model, showing excellent accuracy for our FDM scheme. As a numerical example, we apply the method to simulate seismic waves affected by hemispherical variations of P-wavespeed and attenuation in the top 300km of the inner core.

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.

Application of the Nearly Perfectly Matched Layer to Seismic-Wave Propagation Modeling in Elastic Anisotropic Media

Seismic anisotropy has become a main research topic in both seismic data processing and interpretation (e.g., migration and impedance inversion). As an important part of seismic anisotropy, numerical modeling can help us better understand seismic responses when seismic waves propagate through the interior of the Earth. For the numerical implementation of finite-difference time-domain algorithms, boundary conditions are needed to suppress the spurious reflections from the artificially truncated edges of the discrete model. Here, the nearly perfectly matched layer (NPML) technique is applied to a numerical simulation of seismic-wave propagation in a 2D elastic anisotropic medium, which has been proven to be a very effective absorbing boundary condition in electromagnetic, acoustic, and elastic isotropic media. The equivalence of the wave absorbing performance between the NPML and the standard perfectly matched layer (PML) for elastic wave modeling is theoretically proven in this paper. To check the absorbing ability of the NPML in the cases of grazing incidence and longer simulations (up to 200 s), we choose an elongated anisotropic model in our test. Test results indicate that the NPML has essentially the same computational accuracy as the standard PML but with higher computational savings.

FDM Simulation of an Anomalous Later Phase from the Japan Trench Subduction Zone Earthquakes

We investigated the development of a distinct later phase observed at stations near the Japan Trench associated with shallow, outer-rise earthquakes off the coast of Sanriku, northern Japan based on the analysis of three-component broadband seismograms and FDM simulations of seismic wave propagation using a heterogeneous structural model of the Japan Trench subduction zone. Snapshots of seismic wave propagation obtained through these simulations clearly demonstrate the complicated seismic wavefield constructed by a coupling of the ocean acoustic waves and the Rayleigh waves propagating within seawater and below the sea bottom by multiple reflections associated with shallow subduction zone earthquakes. We demonstrated that the conversion to the Rayleigh wave from the coupled ocean acoustic waves and the Rayleigh wave as they propagate upward along the slope of seafloor near the coast is the primary cause of the arrival of the distinct later phase at the station near the coast. Through a sequence of simulations using different structural models of the Japan Trench subduction zone, we determined that the thick layer of seawater along the trench and the suddenly rising sea bottom onshore of the Japanese island are the major causes of the distinct later phase. The results of the present study indicate that for realistic modeling of seismic wave propagation from the subduction zone earthquakes, a high-resolution bathymetry model is very crucial, although most current simulations do not include a water column in their simulation models.

On 3D modeling of seismic wave propagation via a structured parallel multifrontal direct Helmholtz solver

ABSTRACTWe consider the modeling of (polarized) seismic wave propagation on a rectangular domain via the discretization and solution of the inhomogeneous Helmholtz equation in 3D, by exploiting a parallel multifrontal sparse direct solver equipped with Hierarchically Semi‐Separable (HSS) structure to reduce the computational complexity and storage. In particular, we are concerned with solving this equation on a large domain, for a large number of different forcing terms in the context of seismic problems in general, and modeling in particular. We resort to a parsimonious mixed grid finite differences scheme for discretizing the Helmholtz operator and Perfect Matched Layer boundaries, resulting in a non‐Hermitian matrix. We make use of a nested dissection based domain decomposition, and introduce an approximate direct solver by developing a parallel HSS matrix compression, factorization, and solution approach. We cast our massive parallelization in the framework of the multifrontal method. The assembly tree is partitioned into local trees and a global tree. The local trees are eliminated independently in each processor, while the global tree is eliminated through massive communication. The solver for the inhomogeneous equation is a parallel hybrid between multifrontal and HSS structure. The computational complexity associated with the factorization is almost linear in the size, n say, of the matrix, viz. between O(n log n) and O(n4/3 log n), while the storage is almost linear as well, between O(n) and O(n log n). We exploit the use of a regular (Cartesian) mesh common in many seismic applications.

Three-Dimensional Wave-Field Simulation in Heterogeneous Transversely Isotropic Medium with Irregular Free Surface

Modeling of seismic-wave propagation in anisotropic medium with irregular topography is beneficial to interpret seismic data acquired by active and pas- sive source seismology conducted in areas of interest such as mountain ranges and basins. The major challenge in this context is the difficulty in tackling the irregular free-surface boundary condition in a Cartesian coordinate system. To implement sur- facetopography,weusetheboundary-conforminggridandmaparectangulargridonto a curved grid. We use a stable and explicit second-order accurate finite-difference scheme to discretize the elastic wave equations (in a curvilinear coordinate system) in a 3D heterogeneous transversely isotropic medium. The free-surface boundary con- ditions are accurately applied by introducing a discretization that uses boundary- modified difference operators for the mixed derivatives in the governing equations. The accuracy of the proposed method is checked by comparing the numerical results obtained by the trial algorithm with the analytical solutions of the Lamb's problem, for anisotropicmediumandatransverselyisotropicmediumwithaverticalsymmetryaxis, respectively. Efficiency tests performed by different numerical experiments illustrate clearlytheinfluenceofanirregular(nonflat)freesurfaceonseismic-wavepropagation.

EPcrust: a reference crustal model for the European Plate

We present a new crustal model for the European Plate, derived from collection and critical integration of information selected from the literature. The model covers the whole European Plate from North Africa to the North Pole (20°N–90°N) and from the Mid-Atlantic Ridge to the Urals (40°W–70°E). The chosen parametrization represents the crust in three layers (sediments, upper crust and lower crust), and describes the 3-D geometry of the interfaces and seismologically relevant parameters—isotropic P- and S-wave velocity, plus density—with a resolution of 0.5°× 0.5° on a geographical latitude–longitude grid. We selected global and local models, derived from geological assumptions, active seismic experiments, surface wave studies, noise correlation, receiver functions. Model EPcrust presents significant advantages with respect to previous models: it covers the whole European Plate; it is a complete and internally-consistent model (with all the parameters provided, also for the sedimentary layer); it is reproducible; it is easy to update in the future by adding new contributions; and it is available in a convenient digital format. EPcrust could be used to account for crustal structure in seismic wave propagation modelling at continental scale or to compute linearized crustal corrections in continental scale seismic tomography, gravity studies, dynamic topography and other applications that require a reliable crustal structure. Because of its resolution, our model is not suited for local-scale studies, such as the computation of earthquake scenarios, where more detailed knowledge of the structure is required. We plan to update the model as new data will become available, and possibly improve its resolution for selected areas in the future.

Seismic-Wave Propagation Modeling in Viscoelastic Media Using the Auxiliary Differential Equation Method

In many seismic applications, a constant quality factor is used to describe the constitutive laws of viscoelastic materials, characterized by frequency-independent attenuation characteristics. In such cases, the frequency dependence of the medium’s properties is not taken into account. To overcome this drawback, we proposed an elegant finite difference time domain implementation, with an auxiliary differential equation technique to explicitly solve any stress-strain relation. This scheme is inherited from the formalism of electromagnetism and is based on the separation of the propagation equations from the constitutive law defined in the stress-strain equation. The conventional assumption of a constant quality factor assumption can then be easily avoided in the modeling of seismic-wave propagation in viscoelastic media. We developed such a method and simulated synthetic traces over a simple, 2D viscoelastic homogeneous medium using a Zener model. Wave propagation phase velocities were estimated by means of dispersion analysis and appeared to match theoretical values over a reasonable frequency range. We also measured the material’s attenuation behavior by studying the quality factor, thereby validating our approach.

Fluid–solid coupling on a cluster of GPU graphics cards for seismic wave propagation

We develop a hybrid multiGPUs and CPUs version of an algorithm to model seismic wave propagation based on the spectral-element method in the case of models of the Earth containing both fluid and solid layers. Thanks to the overlapping of communications between processing nodes on the computer with calculation by means of non-blocking message passing, we obtain excellent weak scalability of this finite-element code on a cluster of 192 GPUs and speedup factors of more than one order of magnitude compared to the same code run on a cluster of traditional CPUs. This enables us to show a new geophysical phenomenon concerning wave propagation of diffracted shear waves in a layer called D″ located at the base of the Earth's mantle, namely that in this layer the transverse and radial components of these waves can undergo a relative shift even in an isotropic Earth model, whereas this observation in real seismological data was interpreted until now as an indication of the presence of anisotropy in this layer.

Non-conforming hybrid meshes for efficient 2-D wave propagation using the Discontinuous Galerkin Method

We present a Discontinuous Galerkin finite element method using a high-order time integration technique for seismic wave propagation modelling on non-conforming hybrid meshes in two space dimensions. The scheme can be formulated to achieve the same approximation order in space and time and avoids numerical artefacts due to non-conforming mesh transitions or the change of the element type. A point-wise Gaussian integration along partially overlapping edges of adjacent elements is used to preserve the schemes accuracy while providing a higher flexibility in the problem-adapted mesh generation process. We describe the domain decomposition strategy of the parallel implementation and validate the performance of the new scheme by numerical convergence test and experiments with comparisons to independent reference solutions. The advantage of non-conforming hybrid meshes is the possibility to choose the mesh spacing proportional to the seismic velocity structure, which allows for simple refinement or coarsening methods even for regular quadrilateral meshes. For particular problems of strong material contrasts and geometrically thin structures, the scheme reduces the computational cost in the sense of memory and run-time requirements. The presented results promise to achieve a similar behaviour for an extension to three space dimensions where the coupling of tetrahedral and hexahedral elements necessitates non-conforming mesh transitions to avoid linking elements in form of pyramids.

Stable discontinuous staggered grid in the finite-difference modelling of seismic motion

SUMMARY We present an algorithm of the spatial discontinuous grid for the 3-D fourth-order velocity–stress staggered-grid finite-difference modelling of seismic wave propagation and earthquake motion. The ratio between the grid spacing of the coarser and finer grids can be an arbitrary odd number. The algorithm allows for large numbers of time levels without inaccuracy and eventual instability due to numerical noise inevitably generated at the contact of two grids with different spatial grid spacings. The key feature of the algorithm is the application of the Lanczos downsampling filter. The algorithm of the discontinuous grid is directly applicable also to the displacement-stress staggered-grid finite-difference scheme.

High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster

We implement a high-order finite-element application, which performs the numerical simulation of seismic wave propagation resulting for instance from earthquakes at the scale of a continent or from active seismic acquisition experiments in the oil industry, on a large cluster of NVIDIA Tesla graphics cards using the CUDA programming environment and non-blocking message passing based on MPI. Contrary to many finite-element implementations, ours is implemented successfully in single precision, maximizing the performance of current generation GPUs. We discuss the implementation and optimization of the code and compare it to an existing very optimized implementation in C language and MPI on a classical cluster of CPU nodes. We use mesh coloring to efficiently handle summation operations over degrees of freedom on an unstructured mesh, and non-blocking MPI messages in order to overlap the communications across the network and the data transfer to and from the device via PCIe with calculations on the GPU. We perform a number of numerical tests to validate the single-precision CUDA and MPI implementation and assess its accuracy. We then analyze performance measurements and depending on how the problem is mapped to the reference CPU cluster, we obtain a speedup of 20x or 12x.

Impact of time-lapse seismic data for permeability estimation

We consider the impact of using time-lapse seismic data in addition to production data for permeability estimation in a porous medium with multiphase fluid flows, such as a petroleum reservoir under water-assisted production. Since modeling seismic wave propagation in addition to modeling fluid flows in the reservoir is quite involved, it is assumed that the time-lapse seismic data have already been inverted into fluid saturation differences (pseudoseismic data). Because an inversion process often leads to considerable error growth, we will consider pseudoseismic data with large uncertainties. The impact of pseudoseismic data is assessed through permeability estimation with and without such data and through application of some uncertainty measures for the estimated parameters. A multiscale algorithm is used for the parameter estimations, so that potential differences in attainable permeability resolution will be easily revealed. The numerical examples clearly indicate that the permeability estimation problem is stabilized at a higher level of resolution when pseudoseismic data are applied in addition to production data, even if the pseudoseismic data have large associated uncertainties. Use of the parameter uncertainty measures confirm these results.

Seismic vibrations in bulk blasting with high-precise electronic and nonelectric blasting systems at quarries

The authors report measurement data on seismic waves in bulk blasting at quarries by using new high-precise electronic and pyrotechnic blasting systems. It is proved that both systems are efficient, intensity of seismic waves is much lower in large-scale bulk blasting. The authors implemented numerical modeling of seismic wave propagation under a short-delay bulk blast at a quarry. Influence of the blast delay parameters and their precision on the maximum level of seismic vibrations has been studied, and their optimal ranges have been established.