Modeling Of Seismic Wave Propagation Research Articles

Analysis of transient wave scattering and its applications to site response analysis using the scaled boundary finite-element method

In this paper, the problem of scattering and amplification of seismic waves by topographical and geological irregularities is addressed directly in the time domain through the scaled boundary finite-element method (SBFEM). The quadtree domain decomposition technique is utilized for the SBFEM discretization of the near field. The far field is rigorously modeled by the displacement unit-impulse response matrix. The computational cost of analyses is reduced through using local formulations in space and time. Considering incident fields of obliquely plane waves coming from far field, the seismic wave inputs are formulated as boundary tractions applied to the near field. Implementing all these aspects in the SBFEM, leads to an elegant technique for time-domain modeling of seismic wave propagation in heterogeneous media with topographical irregularities. Four numerical examples considering various site effects and wave patterns demonstrate the accuracy, versatility and applicability of the approach.The approach is straightforwardly applicable to nonlinear and 3D wave scattering problems with complex site effects.

An orthorhombic representation of a heterogeneous medium for the finite-difference modelling of seismic wave propagation

The possibility of applying one explicit finite-difference (FD) scheme to all interior grid points (points not lying on a grid border) no matter what their positions are with respect to the material interface is one of the key factors of the computational efficiency of the FD modelling. Smooth or discontinuous heterogeneity of the medium is accounted for only by values of the effective grid moduli and densities. Accuracy of modelling thus very much depends on how these effective grid parameters are evaluated. We present an orthorhombic representation of a heterogeneous medium for the FD modelling. We numerically demonstrate its superior accuracy. Compared to the harmonic-averaging representation the orthorhombic representation is more accurate mainly in the case of strong surface waves that are especially important in local surface sedimentary basins. The orthorhombic representation is applicable to modelling seismic wave propagation and earthquake motion in isotropic models with material interfaces and smooth heterogeneities using velocity–stress, displacement–stress and displacement FD schemes on staggered, partly staggered, Lebedev and collocated grids.

An improved immersed boundary finite-difference method for seismic wave propagation modeling with arbitrary surface topography

To accurately simulate seismic wave propagation for the purpose of developing modern land data processing tools, especially full-waveform inversion (FWI), we have developed an efficient high-order finite-difference forward-modeling algorithm with the capability of handling arbitrarily shaped free-surface topography. Unlike most existing forward-modeling algorithms using curvilinear grids to fit irregular surface topography, this finite-difference algorithm, based on an improved immersed boundary method, uses a regular Cartesian grid system without suffering from staircasing error, which is inevitable in a conventional finite-difference method. In this improved immersed boundary finite-difference (IBFD) algorithm, arbitrarily curved surface topography is accounted for by imposing the free-surface boundary conditions at exact boundary locations instead of using body-conforming grids or refined grids near the boundaries, thus greatly reducing the complexity of its preprocessing procedures and the computational cost. Furthermore, local continuity, large curvatures, and subgrid curvatures are represented precisely through the employment of the so-called dual-coordinate system — a local cylindrical and a global Cartesian coordinate. To properly describe the wave behaviors near complex free-surface boundaries (e.g., overhanging structures and thin plates, or other fine geometry features), the wavefields in a ghost zone required for the boundary condition enforcement are reconstructed accurately by introducing a special recursive interpolation technique into the algorithm, which substantially simplifies the boundary treatment procedures and further improves the numerical performance of the algorithm, as demonstrated by the numerical experiments. Numerical examples revealed the performance of the IBFD method in comparison with a conventional finite-difference method.

2D spectral element modeling of GPR wave propagation in inhomogeneous media

We present a spectral element method, for simulation of ground-penetrating radar (GPR) in two dimensions. The technique is based upon a weak formulation of the equations of Maxwell and combines the flexibility of the elemental-based methods with the accuracy of the spectral based methods. The wave field on the elements is discretized using high-degree Lagrange interpolation and integration over an element is accomplished based upon the Gauss–Lobatto–Legendre integration rule. As a result, the mass matrix and the damping matrix are always diagonal, which drastically reduces the computational cost. We first develop the formulation of 2D spectral element method (SEM) in the time-domain based on Maxwell's equations. The presented formulation is with matrix notation that simplifies the implementation of the relations in computer programs, especially in MATLAB application. We discuss the differences between spectral element method and finite-element method in the time-domain. Also, we show that the SEM numerical dispersion is much lower than FEM. To absorb waves at the edges of the modeling domain, we implement first order Clayton and Engquist absorbing boundary conditions (CE-ABC) introduced in numerical finite-difference modeling of seismic wave propagation. We used the SEM to simulate a complex model to show its abilities and limitations. As well as, one distinct advantage of SEM is that we can easily define our model features in nodal points, because the integration points and the interpolation points are similar that makes it very flexible in simulation of complex models.

Numerical Investigation on the Feasibility of Estimating the Thickness of Europa’s Ice Shell by a Planned Impact

AbstractThe existence of a liquid water ocean under the icy surface of Jupiter’s moon Europa has been considered highly probable for decades. It has been speculated that such an environment could host extraterrestrial life. However, the thickness of Europa’s top layer of ice and the depth of the ocean underneath are currently unknown. The authors investigate the feasibility of using seismic-inversion analysis to estimate the thickness of the ice shell and the depth of the subsurface ocean in a hypothetical event, in which a spacecraft-released impactor makes a planned collision with Europa’s surface. The wave motions on Europa’s surface that would result from a planned impact could be measured by impactor-probes penetrated into the ice shell. This work uses finite element method (FEM) modeling of seismic wave propagation to compute wave motions with a guessed ice thickness and ocean depth as input. Inverse modeling is used to determine the input that produces the best results when comparing the computed w...

Finite-difference strategy for elastic wave modelling on curved staggered grids

Waveform modelling is essential for seismic imaging and inversion. Because including more physical characteristics can potentially yield more accurate Earth models, we analyse strategies for elastic seismic wave propagation modelling including topography. We focus on using finite differences on modified staggered grids. Computational grids can be curved to fit the topography using distribution functions. With the chain rule, the elasto-dynamic formulation is adapted to be solved directly on curved staggered grids. The chain-rule approach is computationally less expensive than the tensorial approach for finite differences below the 6th order, but more expensive than the classical approach for flat topography (i.e. rectangular staggered grids). Free-surface conditions are evaluated and implemented according to the stress image method. Non-reflective boundary conditions are simulated via a Convolutional Perfect Matching Layer. This implementation does not generate spurious diffractions when the free-surface topography is not horizontal, as long as the topography is smoothly curved. Optimal results are obtained when the angle between grid lines at the free surface is orthogonal. The chain-rule implementation shows high accuracy when compared to the analytical solution in the case of the Lamb’s problem, Garvin’s problem and elastic interface.

Технология суперкомпьютерного 3D моделирования сейсмических волновых полей в сложно построенных средах

The paper considered computing technology solving problems related to the modeling of seismic wave propagation in inhomogeneous media typical of volcanic structures using supercom-puter simulations in order to create systems of vibroseis monitoring for quake-prone objects. The physico-mathematical model of the magmatic volcano is constructed and software implementation on the basis of the known numerical method that effectively using the architecture of modern su-percomputers equipped with GPU is developed. The parallel 2D and 3D algorithms and software for simulation of elastic wave propagation in a complicated medium (2D model is separation of original 3D model using various angles and planes) on basis of the explicit finite-difference scheme for the shifted grids and CFS-PML method of absorbing boundaries is developed. Scalability of algorithms is investigated. The application of the developed technology allows for much more effi-cient to carry out studies of the structure of the wave field due to the geometry of the internal boundaries and refinement of its kinematic and dynamic characteristics.

Wave propagation in soils problems using the Generalized Finite Difference Method

The consideration of highly irregular topography or heterogeneities and singularities in the domain are necessary for solving a wide variety of seismic problems, e. g. earthquake site studies. Therefore, the use of a meshless method (MM) with the possibility of employing an irregular grid point distribution can be of interest for modeling this kind of problems. The tradition of using Finite Different Methods for modeling seismic wave propagation problems have allowed researches to solve it and, in fact, a numerical model using the geotechnical calculation software FLAC has been used here to validate the results. It is clear the improvement of using meshless method techniques to incorporate irregularities in a natural way.

Multicomponent reduced scale seismic modelling: upgrade of the MUSC laboratory with application to polarization observations

Reduced scale physical modelling is a useful intermediate step to fill the gap between field measurements and numerical simulations. In many different areas of seismology, multicomponent measurements are now commonly recorded and processed. However, up to now, laboratory facilities providing flexible and accurate reproduction of large multicomponent seismic experiments have not reached maturity. Within this context, we present an improvement of a measurement facility (MUSC) developed previously to enable the flexible and versatile reproduction of seismic experiments at reduced scale. This new measurement system provides simultaneous measurements of the vertical and the horizontal components of seismic wavefields by taking advantage of recent developments in laser ultrasonic sensing technologies. Reference measurements on an aluminum block are carried out in order to quantitatively evaluate multicomponent measurements. Noise characterization, surface wave ellipticity analysis and waveform comparisons with a semi-analytical model are presented. To enable the simultaneous comparison of phase and polarization attributes, we also introduce a regularized version of the polarization attributes and apply it to synthetic and experimental data. These comparisons demonstrate that high fidelity multicomponent reduced scale physical modelling is possible with this new facility. An illustrative example with a model containing a shallow cavity is presented. In this case, multicomponent measurements help to identify the cavity signature but multiples have a detrimental effect on regularized instantaneous polarization attributes. By assessing the reliability of the MUSC measurement facility in the performance of massive high-quality multicomponent data acquisition at the laboratory scale, this study highlights the interesting potentials of reduced scale modelling of seismic wave propagation for a large variety of applications

Finite difference elastic wave modeling with an irregular free surface using ADER scheme

In numerical modeling of seismic wave propagation in the earth, we encounter two important issues: the free surface and the topography of the surface (i.e. irregularities). In this study, we develop a 2D finite difference solver for the elastic wave equation that combines a 4th- order ADER scheme (Arbitrary high-order accuracy using DERivatives), which is widely used in aeroacoustics, with the characteristic variable method at the free surface boundary. The idea is to treat the free surface boundary explicitly by using ghost values of the solution for points beyond the free surface to impose the physical boundary condition. The method is based on the velocity-stress formulation. The ultimate goal is to develop a numerical solver for the elastic wave equation that is stable, accurate and computationally efficient. The solver treats smooth arbitrary-shaped boundaries as simple plane boundaries. The computational cost added by treating the topography is negligible compared to flat free surface because only a small number of grid points near the boundary need to be computed. In the presence of topography, using 10 grid points per shortest shear-wavelength, the solver yields accurate results. Benchmark numerical tests using several complex models that are solved by our method and other independent accurate methods show an excellent agreement, confirming the validity of the method for modeling elastic waves with an irregular free surface.

Seismic attenuation in fractured porous media: insights from a hybrid numerical and analytical model

Seismic attenuation in fluid-saturated porous rocks can occur by geometric spreading, wave scattering or the internal dissipation of energy, most likely due to the squirt-flow mechanism. In principle, the pattern of seismic attenuation recorded on an array of sensors contains information about the medium, in terms of material heterogeneity and anisotropy, as well as material properties such as porosity, crack density, and pore-fluid composition and mobility. In practice, this inverse problem is challenging. Here we provide some insights into the effects of internal dissipation by analysing synthetic data produced by a hybrid numerical and analytical model for seismic wave propagation in a fractured medium embedded within a layered geological structure. The model is made up of one anisotropic and three isotropic horizontal layers. The anisotropic layer consists of a porous, fluid-saturated material containing vertically aligned inclusions representing a set of fractures. This combination allows squirt-flow to occur between the pores in the matrix and the model fractures. Our results show that the fluid mobility and the associated relaxation time of the fluid-pressure gradient control the frequency range over which attenuation occurs. This induced attenuation increases with incidence angle and azimuth away from the fracture strike-direction. Azimuthal variations in the induced attenuation are elliptical allowing the fracture orientations to be obtained from the axes of the ellipse. These observations hold out the potential of using seismic attenuation as an additional diagnostic in the characterisation of rock formations for a variety of applications including hydrocarbon exploration and production, subsurface storage of CO2, and geothermal energy extraction.

Acoustic and elastic modeling by optimal time-space-domain staggered-grid finite-difference schemes

Staggered-grid finite-difference (SFD) methods are widely used in modeling seismic-wave propagation, and the coefficients of finite-difference (FD) operators can be estimated by minimizing dispersion errors using Taylor-series expansion (TE) or optimization. We developed novel optimal time-space-domain SFD schemes for acoustic- and elastic-wave-equation modeling. In our schemes, a fourth-order multiextreme value objective function with respect to FD coefficients was involved. To yield the globally optimal solution with low computational cost, we first used variable substitution to turn our optimization problem into a quadratic convex one and then used least-squares (LS) to derive the optimal SFD coefficients by minimizing the relative error of time-space-domain dispersion relations over a given frequency range. To ensure the robustness of our schemes, a constraint condition was imposed that the dispersion error at each frequency point did not exceed a given threshold. Moreover, the hybrid absorbing boundary condition was applied to remove artificial boundary reflections. We compared our optimal SFD with the conventional, TE-based time-space-domain, and LS-based SFD schemes. Dispersion analysis and numerical simulation results suggested that the new SFD schemes had a smaller numerical dispersion than the other three schemes when the same operator lengths were adopted. In addition, our LS-based time-space-domain SFD can obtain the same modeling accuracy with shorter spatial operator lengths. We also derived the stability condition of our schemes. The experiment results revealed that our new LS-based SFD schemes needed a slightly stricter stability condition.

Dreamlet-based interpolation using POCS method

Due to incomplete and non-uniform coverage of the acquisition system and dead traces, real seismic data always has some missing traces which affect the performance of a multi-channel algorithm, such as Surface-Related Multiple Elimination (SRME), imaging and inversion. Therefore, it is necessary to interpolate seismic data. Dreamlet transform has been successfully used in the modeling of seismic wave propagation and imaging, and this paper explains the application of dreamlet transform to seismic data interpolation. In order to avoid spatial aliasing in transform domain thus getting arbitrary under-sampling rate, improved Jittered under-sampling strategy is proposed to better control the dataset. With L0 constraint and Projection Onto Convex Sets (POCS) method, performances of dreamlet-based and curvelet-based interpolation are compared in terms of recovered signal to noise ratio (SNR) and convergence rate. Tests on synthetic and real cases demonstrate that dreamlet transform has superior performance to curvelet transform.

Finite-difference time-domain method for modelling of seismic wave propagation in viscoelastic media

A finite difference formulation for the equations of SH waves in viscoelastic media has been developed and applied to the problem of an infinite half space with a buried source emitting a compressional pulse. A dispersion relation has been developed by employing the finite difference (FD) method, which has been used to discuss the stability criteria and also for obtaining the relation for phase and group velocities in viscoelastic media. Furthermore, the effect of variation of stability parameter on phase and group velocities; effect on the wave propagation with increase in dispersion parameter; and change in the wave propagation at different locations from the source have also been studied. The SH waves are found to be dispersive in viscoelastic media. The phase velocity increases with decrease in stability parameter. Furthermore, the location of the source has an important effect on the strength of the wave propagation in that the intensity of wave and therefore its longevity is inversely proportional to the distance from the source, implying that the wave will progressively die out exponentially.

Time-Lapse Seismic Data Integration Incorporates Pressure and Saturation Effects

This article, written by Special Publications Editor Adam Wilson, contains highlights of paper SPE 166395, ’Streamline-Based Time-Lapse Seismic-Data Integration Incorporating Pressure and Saturation Effects,’ by Shingo Watanabe, SPE, Jichao Han, SPE, Akhil Datta-Gupta, SPE, and Michael J. King, SPE, Texas A&M University, prepared for the 2013 SPE Annual Technical Conference and Exhibition, New Orleans, 30 September-2 October. The paper has not been peer reviewed. This paper presents an efficient history matching approach that simultaneously integrates 4D repeat seismic surveys with well production data. This approach is particularly well-suited for the calibration of the reservoir properties of high-resolution geologic models because the seismic data are areally dense but sparse in time, while the production data are continuous in time but averaged over interwell spacing. The joint history matching is performed using streamline-based sensitivities derived from either finite-difference or streamline-based flow simulation. Introduction Previous studies have made attempts at quantitative 4D seismic history matching. However, the reconciliation of reservoir model heterogeneity with temporal changes in seismic attributes remains a particularly complex task. Several dynamic data integration algorithms have been proposed in the literature, which can be categorized broadly into three data integration levels: (1) reservoir-simulation- based integration between the pressure and saturation estimates inverted from seismic observation data and the direct simulated saturation and pressure responses, (2) petroelastic-based integration between the seismic inverted rock elastic properties derived from a geophysical inversion and the simulated rock elastic responses from the simulated saturation and pressure responses by means of petroelastic models, and (3) seismic-forward-model-based integration between the direct seismic traces attributes and the simulated seismic responses by means of seismic wave propagation modeling. The seismic-forward-model-based seismic data integration approach uses direct seismic trace data and circumvents the uncertainties associated with the seismic-inversion process. The general workflow for generating synthetic seismic data from a reservoir simulation model generally involves (1) static reservoir properties on the simulation grid, (2) simulation of dynamic pressure and fluid saturations at each cell, (3) computation of seismic elastic properties, and (4) simulation of the seismic attribute by applying a seismic wave propagation model over the reservoir interval and the overburden rock. This study examines both reservoir-simulation- based calibration and rock-physics- based calibration. In the first instance, the authors calibrate against time-lapse changes in both pressure and saturation. In the second instance, they calibrate against the time-lapse change in acoustic impedance. In both cases, the inversion will provide local updates to the reservoir model using analytic streamline- based sensitivity calculations.

Numerical modeling of seismic and acoustic-gravity waves propagation in an “Earth-Atmosphere” model in the presence of wind in the air

In this paper, an effective numerical algorithm for 2.5D seismic and acoustic-gravitational wave propagation is applied to a combined “Earth-Atmosphere” model in the presence of wind in the air. Seismic wave propagation in an elastic half-space is described by a system of first-order dynamic equations of elasticity theory. The propagation of acoustic-gravitational waves in the atmosphere in the presence of wind is described by the linearized Navier-Stokes equations. The algorithm is based on the integral Laguerre transform with respect to time, the finite integral Fourier transform with respect to a spatial coordinate combined with a finite difference method for the reduced problem.

Perfectly matched layer-absorbing boundary condition for finite-element time-domain modeling of elastic wave equations

The perfectly matched layer (PML) is a highly efficient absorbing boundary condition used for the numerical modeling of seismic wave equation. The article focuses on the application of this technique to finite-element time-domain numerical modeling of elastic wave equation. However, the finite-element time-domain scheme is based on the second-order wave equation in displacement formulation. Thus, the first-order PML in velocity-stress formulation cannot be directly applied to this scheme. In this article, we derive the finite-element matrix equations of second-order PML in displacement formulation, and accomplish the implementation of PML in finite-element time-domain modeling of elastic wave equation. The PML has an approximate zero reflection coefficients for bulk and surface waves in the finite-element modeling of P-SV and SH wave propagation in the 2D homogeneous elastic media. The numerical experiments using a two-layer model with irregular topography validate the efficiency of PML in the modeling of seismic wave propagation in geological models with complex structures and heterogeneous media.

Finite-frequency sensitivity kernels for two-station surface wave measurements

SUMMARY We present a method for the computation offinite-frequency sensitivity kernels for two-station surface wave measurements. It is based on the combination of spectral-element modelling of seismic wave propagation, adjoint techniques and traveltime estimates from two-station crosscorrelations of seismograms. The analysis of sensitivity kernels for a 1-D earth model resulted in two major conclusions: (1) The finite-frequency sensitivity is zero along the interstation ray path for group velocity measurements obtained by cross-correlations. It follows that interstation group velocity measurements should not be made by cross-correlation if a raybasedinterpretationisused.(2)Althoughsensitivityalongtheinterstationraypathisdominant for phase velocity measurements, sensitivity far from the ray path can be large, depending on the details of the source–receiver geometry. The complexities of finite-frequency two-station sensitivity kernels should be taken into account to avoid misinterpretations and to improve the quality of tomographic inversions.

A Coupled Method of Spectral Element and Normal Mode: A Powerful Tool for Modeling the Seismic Waves Propagation in the Anisotropic D″ Region

Abstract:We extended the coupled mode/spectral element method (CSEM hereafter) to incorporate three‐dimensional, anisotropic media in the lowest 300 km above the core‐mantle boundary, the D″ area. This extended method was applied to modeling of seismic waves propagation in earth models under the so‐called “sandwich” structure, which uses a spherical shell of spectral elements between two modal solutions to represent the D″ area. Using a spherically symmetric reference model with modified anisotropic parameters in the spherical shell, the CSEM yields waveforms in excellent agreement with that of normal‐mode. The two methods match very well for all body wave arrivals with periods greater than 8 s. Besides, the S, ScS, SKS, and SKKS phases in the synthetic seismograms of the transverse isotropic model show obvious arrival time differences between SH and SV components. This proves the capability of the extended CSEM code for modeling S wave splitting caused by the radial anisotropic media in the D″ area.

Geomechanical Modeling of Fault Responses and the Potential for Notable Seismic Events During Underground CO2 Injection

We summarize a number of recent modeling studies related to the potential for fault reactivations and induced seismicity during underground CO2 injection. The model simulations were conducted using coupled multiphase fluid flow and geomechanics, including fault-frictional weakening enabling analysis of sudden (seismic) fault rupture, with some of the numerical analyses extended to dynamic modeling of seismic source, wave propagation, and ground motion. These model simulations show that the critical factors in determining the likelihood and magnitude of such an event are the local in situ stress field, fault orientation and size, fault strength, and injection pressure. We analyzed the case of activation of a 1km long minor fault that might have gone undetected during the site investigation and show that the maximum seismic magnitudes would likely be less than about 3.6, even if the entire 1km fault would to be activated. We then include seismic wave propagation generated by the rupture and show how the acceleration and deceleration of the rupture generate waves and result in a peak ground acceleration of about 0.1g, except for a localized –0.6g of horizontal peak acceleration at the faults intersection with the ground surface. The modeling shows that these are high frequency events that would not cause any substantial damage but could certainly be felt by the local population. We may also considered that fault reactivation, even associated with relatively small seismic or aseismic events, could potentially increase CO2 seepage out of the intended storage complex and therefore reduce the effectiveness of a CO2 storage operation. Under these circumstances, we recommend a staged, learn-as-you-go approach, involving a gradual increase of injection rates combined with continuous monitoring of geomechanical changes, as well as siting beneath a multiple layered overburden for multiple flow barrier protection, should an unexpected deep fault activation occur.