Modeling Of Seismic Wave Propagation Research Articles

Modeling a pure qP wave in attenuative transverse isotropic media using a complex-valued fractional viscoacoustic wave equation

Seismic anisotropy, characterized by direction-dependent seismic velocity and attenuation, is a fundamental property of the earth, widely observed in field and laboratory measurements. This anisotropy significantly influences seismic wavefield attributes, such as amplitude, phase, and traveltime. Accurate quantification of anisotropy-related seismic wave propagation is essential for global- and regional-scale studies in seismology, such as imaging and inversion. Because the elastic assumption of the earth requires contending with S waves, which are often weak in our seismic recording, we derive a new general complex-valued spatial fractional Laplacian (SFL) viscoacoustic anisotropic pure quasi-P (qP) wave equation (WE) in attenuative transverse isotropic (A-TI) media under the acoustic assumption and use the fixed-order SFL operator to quantify seismic wave propagation in heterogeneous media. Our WE accurately captures anisotropic effects in seismic velocity and attenuation, even in strongly attenuative and anisotropic environments. It also describes the decoupled energy attenuation and velocity dispersion behaviors of the frequency-independent quality factor Q, beneficial for seismic imaging and inversion. Compared with existing viscoacoustic anisotropic pure qP WEs in A-TI media, our WE simplifies computational complexity and is easy to implement, as it decomposes into a complex-valued scalar operator and two differential operators. The evaluated dispersion curves further confirm the accuracy of our approach. Numerical examples in complex geologic settings demonstrate the new equation’s precision, stability, and efficiency in modeling seismic wave propagation.

Efficient pure qP-wave modeling and reverse time migration in tilted transversely isotropic media calculated by a finite-difference approach

Subsurface anisotropy is commonly induced by shale layers, aligned cracks, and fine bedding and has a significant impact on seismic wave propagation. Ignoring anisotropic effects during seismic migration degrades image quality. Therefore, we derive a pure qP-wave equation with high accuracy for modeling seismic wave propagation in tilted transversely isotropic (TTI) media. However, the derived pure qP-wave equation requires a computationally expensive spectral-based method for performing numerical simulations. This is unsuitable for large-scale industrial applications, particularly 3D applications. For numerical efficiency, we first decompose the newly derived wave equation into some fractional differential operators and spatial derivatives. The spatial derivatives can be directly solved by conventional finite-difference (FD) approaches. Then, we use an asymptotic approximation to find an equivalent form of fractional differential operators, obtaining scalar operators that we can discretize with the FD method. Numerical tests show that our TTI pure qP-wave equation with an FD discretization can produce accurate and highly efficient wavefield simulations in TTI media. We also use our TTI pure qP-wave equation with an FD discretization as a forward engine to implement TTI reverse time migration (RTM). Synthetic examples and a field data test demonstrate that our TTI RTM can effectively correct the anisotropic effects, providing high-quality imaging results while maintaining good computational efficiency.

SeismicNet: Physics-informed neural networks for seismic wave modeling in semi-infinite domain

Recently, there has been an increasing interest in leveraging physics-informed neural networks (PINNs) for modeling dynamical systems. However, limited studies have been conducted along this horizon on seismic wave modeling tasks. A critical challenge is that these geophysical problems are typically defined in large domains (i.e., semi-infinite), which leads to high computational costs. We present a new PINN model for seismic wave modeling in semi-infinite domain without the need for labeled data. Specifically, the absorbing boundary condition is introduced into the network as a soft regularizer for handling truncated boundaries. To scale up, we consider a sequential training strategy via temporal domain decomposition to improve the scalability of the network and solution accuracy. Moreover, we design a novel surrogate modeling strategy to account for parametric loading, which estimates the wave propagation in semi-infinite domain given the seismic loading at different locations. Various numerical experiments are implemented to evaluate the performance of the proposed PINN model in the context of forward modeling of seismic wave propagation. In particular, we use diverse material distributions to test the versatility of this approach. The results demonstrate excellent solution accuracy under distinctive scenarios.

Physics-constrained neural networks for half-space seismic wave modeling

Forward modeling of seismic waves using physics-informed neural networks (PINNs) has attracted much attention. However, a notable challenge arises when modeling seismic wave propagation in large domains (i.e., a half-space), PINNs may encounter the issue of "soft constraint failure". To address this problem, we propose a novel framework called physics-constrained neural networks (PCNNs) specifically designed for modeling seismic wave propagation in a half-space. The method of images is incorporated to effectively implement the free stress boundary conditions of the Earth's surface, leading to the successful propagation of plane waves and cylindrical waves in a half-space. We analyze the training dynamics of neural networks when solving two-dimensional (2D) wave equations from the neural tangent kernel (NTK) perspective. An adaptive training algorithm is introduced to mitigate the unbalanced gradient flow dynamics of the different components of the loss function of PINNs/PCNNs. Furthermore, to tackle the complex behavior of seismic waves in layered media, a sequential training strategy is considered to enhance network scalability and solution accuracy. The results of numerical experiments demonstrate the accuracy and effectiveness of our approach.

Numerical modeling of seismic wave propagation in loosely deposited partially saturated sands: an application to a mine dump monitoring case

Extensive mine dumps consisting of loosely deposited sands have been created as a result of open-pit lignite mining, with a risk of soil liquefaction under high water saturation and a corresponding initiating event. Soil compaction is one of the feasible methods for reducing the probability of liquefaction. For the monitoring of liquefaction events and the evaluation of compaction work, seismic survey methods with sensitivity to changes in soil saturation and structure may thus complement other methods. Compared to exploration methods for deep systems, the shallow subsurface presents some unique challenges. To this end, an open-source, customizable code based on Biot’s theory was developed in the FEniCS library, which takes into account partial saturation and porosity dependence of stiffness, permeability, and other quantities. Following code verification, a comprehensive investigation of parameter studies is conducted, from which the effects of different factors on wave propagation characteristics were obtained. The numerical model was applied to simulate the expected changes in seismic response following soil compaction. Furthermore, the position of the high saturation area could be detected from the reflection and refraction P waves. The goal of this work is to provide an analysis framework for the assessment of compaction works and monitoring liquefiable soils in mine dumps under conditions of variable saturation due to rising groundwater tables.

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

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

Distinguishing Unique Earthquakes with Overlapping Signals in Oklahoma

AbstractDuring routine operations monitoring Oklahoma earthquakes, we found that certain earthquakes occurred closely both in space and time and had overlapping phases at the recording stations. Through further scrutiny and analysis, we determined that rather than being distinctly different earthquakes, some of the earthquakes exhibited multiphase arrivals and longer than expected coda due to unique ray paths that encounter impedance contrasts such as at the sedimentary rock-basement. Of course, some of these events truly were distinct events, which we term overlapping earthquakes, for which perceived coda duration overlaps and obscures the phase arrivals of the second event due to the source proximity in both time and space. We detail our classification scheme to separate the local earthquakes in Oklahoma as single, overlapping earthquakes, or those associated with multiphase arrivals. We forward model seismic wave propagation in a 2D crustal model and develop a methodology that utilizes waveform correlation to distinguish phases from overlapping earthquakes to those from crustal reverberations. Duration analysis shows a more elongated duration, qualitatively similar to the duration produced by overlapping earthquakes, at the sites where multiphase arrivals are observed.

Love wave fields in a non-local elastic model with reinforced and inhomogeneous media

The primary intent of this paper is to study the propagation of a Love-type wave and its failure in a model consisting of a non-local elastic fiber-reinforced material and a continuously inhomogeneous non-local elastic semi-infinite medium. The inhomogeneity parameters in the lower medium corresponding to the rigidity and density vary quadratically with the z-direction. Expressions for displacement components of the Love-type wave in both the media are derived separately by solving the second-order hyperbolic type differential equation. In the lower medium, we obtain approximate solutions to the displacement component using asymptotic expansions of the Bessel functions. As a result, the dispersion relation obtained by applying appropriate boundary conditions involves inhomogeneity and non-local elastic parameters of both the media. Furthermore, the time-dependent Finite Difference Scheme (FDS) is used to model the propagation of the Love-type wave, and stability conditions are developed to derive the expressions for the phase and group velocities. This scheme is one of the most efficient schemes in modeling seismic wave propagation because of its higher accuracy, applicability to complex structures, and computational benefits. The numerical and graphical discussions are then carried out using MATLAB software to analyze the stability of the Finite Difference Scheme and also to describe the limiting behavior of the Love-type wave due to the non-local elasticity. The developed model provides a more realistic approach to solving the seismological problems concerned with reinforcement due to the inclusion of non-local elastic behavior in a layered media.

A two-step approach combining FK with SE for simulating ground motion due to point dislocation sources

It has been well established that source-path-site effects have had an essential impact on the seismic analysis of major infrastructure over the past few decades, where the accurate and efficient modeling of seismic wave propagation can become a critical issue. In this study, a two-step approach combining the frequency-wavenumber (FK) technique with the spectral element (SE) method is developed based on the concept of domain reduction to simulate three-dimensional (3D) ground motion due to point dislocation sources. In the first step, the seismic responses for the dislocations in a layered half-space crustal model are solved exactly using the FK method and converted into effective input seismic loads around the region of interest. In the second step, the local wavefields of a 3D complex site are finely simulated using the SEM, while absorbing conditions are imposed at five boundaries of the SE model to achieve absorption of the outward scattered energy. The semi-analytical FK method is well-matched to the propagation frequencies because it allows the broadband synthesis associated with engineered systems to be addressed without additional calculations. When the FK is combined with the SEM, it is feasible to efficiently simulate earthquake ground motions covering a broad range of frequencies in the region of interest. Comparisons with the results of existing studies and the traditional SEM show that our approach is accurate and robust for arbitrary geometries and point sources. To further illustrate the efficiency and general applicability of the proposed method, three examples are presented. The results demonstrate that the two-step approach can consider the effects of the seismic source, propagation path, and local site geological conditions on the earthquake ground motions and serve as an effective forward modeling tool for reducing the computational cost.

Seismic waves in medium with poroelastic/elastic interfaces: a two-dimensionalP-SVfinite-difference modelling

SUMMARYWe present a new methodology of the finite-difference (FD) modelling of seismic wave propagation in a strongly heterogeneous medium composed of poroelastic (P) and (strictly) elastic (E) parts. The medium can include P/P, P/E and E/E material interfaces of arbitrary shapes. The poroelastic part can be with (i) zero resistive friction, (ii) non-zero constant resistive friction or (iii) JKD model of the frequency-dependent permeability and resistive friction. Our FD scheme is capable of subcell resolution: a material interface can have an arbitrary position in the spatial grid. The scheme keeps computational efficiency of the scheme for a smoothly and weakly heterogeneous medium (medium without material interfaces). Numerical tests against independent analytical, semi-analytical and spectral-element methods prove the efficiency and accuracy of our FD modelling. In numerical examples, we indicate effect of the P/E interfaces for the poroelastic medium with a constant resistive friction and medium with the JKD model of the frequency-dependent permeability and resistive friction. We address the 2-D P-SV problem. The approach can be readily extended to the 3-D problem.

Accelerating high‐order stencils on GPUs

SummaryFinite‐difference methods based on high‐order stencils are commonly used for modeling of seismic wave propagation, weather forecasting, computational fluid dynamics, convolutional neural networks, and others. Nowadays, the community commonly employs graphics processing units (GPUs) to accelerate such stencil computations. As a result, knowing how to write efficient stencil computations for GPUs is of significant interest. While high‐performance, low‐order stencils on GPUs have been studied extensively in the literature, not all proposed approaches work well for high‐order stencils. Furthermore, coping with boundary conditions used with stencils for seismic modeling makes it challenging to efficiently exploit thread‐level parallelism on GPUs. In this article, we describe several implementations of a 25‐point stencil. We evaluate our stencil code shapes, memory hierarchy usage, data access patterns, and other performance attributes on several modern GPUs and compare them with machine rooflines. On average, our top‐performing kernels achieve six times the performance of a 25‐point stencil code developed in C and mapped to GPUs using OpenACC. Several of our implementations have excellent performance portability across multiple generations of both NVIDIA and AMD GPUs.

Vertical to Horizontal Spectral Ratio (VHSR) Response of Seismic Wave Propagation in a Homogeneous Elastic – Poroelastic Medium Using The Spectral Finite Element Method

<p>Numerical modeling of 2D seismic wave propagation using spectral finite element method to estimate the response of seismic waves passing through the poroelastic medium from a hydrocarbon reservoir has been carried out. A hybrid simple model of the elastic - poroelastic - elastic with a mesoscopic scale element size of about 50cm was created. Seismic waves which was in the form of the ricker function are generated on the first elastic medium, propagated into the poroelastic medium and then transmitted to the second elastic medium. Pororoelastic medium is bearing hydrocarbon fluid in the form of gas, oil or water. Vertical and horizontal component of velocity seismograms are recorded on all mediums. Seismograms which are recorded in the poroelastic and second elastic medium show the existence of slow P compressional waves following fast P compressional waves that do not appear on the seismogram of the first elastic medium. The slow P wave is generated when the fast P wave enters the interface of the elastic - poroelastic boundary, propagated in the poroelastic medium and is transmited to the second elastic medium. The curves of Vertical to horizontal spectrum ratio (VHSR) which are observed from seismograms recorded in the poroelastic and the second elastic medium show that the peak of VHSR values at low frequency correlated with the fluid of poroelastic reservoir. The highest VHSR value at the low frequency which is recorded on the seismogram is above the 2.5 Hz frequency for reservoirs containing gas and oil in the second elastic medium, while for the medium containing water is the highest VHSR value is below the 2.5 Hz frequency.</p>

Applying an advanced temporal and spatial high-order finite-difference stencil to 3D seismic wave modeling

Present finite-difference (FD) algorithms for modeling seismic wave propagation in 3D acoustic media are mainly based on the traditional orthogonal stencil and can achieve spatial high-order but temporal second-order accuracy. Therefore, these approaches may result in visible temporal dispersion and even instability for relatively large time sampling intervals. In this paper, we develop an advanced FD stencil with temporal and spatial arbitrary even-order accuracy to solve the 3D acoustic wave equation. The temporal derivative is approximated by an octahedron stencil with the operator length N together with the traditional second-order FD. Meanwhile, the spatial derivatives are calculated by the orthogonal stencil with the operator length M. The plane-wave theory is introduced to derive the time-space-domain dispersion relation and estimate the corresponding high-order FD coefficients by Taylor expansion. To achieve higher accuracy, we further optimize FD coefficients by applying least-squares (LS) to minimize the relative error of the time-space-domain dispersion relation. Dispersion and stability analyses show that our new schemes have higher modeling accuracy and propagation stability than the conventional method. Numerical examples demonstrate that the proposed FDs can generate greater temporal accuracy on the prerequisite of preserving satisfactory spatial accuracy. Our new temporal high-order FD stencil by incorporating larger time step, shorter operator lengths, LS-based coefficients and variable-length operators can be regarded as a kind of accurate and efficient tool for seismic wave modeling.

МЕХАНІЗМ ВОГНИЩА ЛИТОВСЬКОГО ЗЕМЛЕТРУСУ НА ОСНОВІ ІНВЕРСІЇ ХВИЛЬОВИХ ФОРМ

The aim: Determination of focal mechanism of Lithuanian earthquake of 12.06.15 (t0 = 08:18:26.4; 55.52° N, 21.40° E; hs = 0.9 км.; ML = 2.6) by waveform inversion using direct waves and a limited number of stations. Method: Matrix method is used for modelling of seismic wave propagation in the medium modelled as horizontally layered heterogeneous elastic structure. There were obtained the relations of displacement waves on the free surface that were used for seismic tensor determination using only direct P- and S- waves. Determination of seismic tensor and the focal mechanism on the base of developed method for a point source is described. Thus, based on forward modeling, numerical techniques are developed for the inversion of observed waveforms for the components of moment tensor. Results: In the paper, a method is presented for the focal mechanism determination of Lithuanian earthquake of 12.06.15 (ML = 2.6) by waveform inversion using limited number of stations. The focal mechanism is determined using the data from two stations: PABE, SLIТ and from three stations: PABE, MTSE, SLIТ. These seismic stations are the part of BAVSEN (BalticVirtualSeismicNetwork). Scientific novelty: 1. In the paper, a method is presented for moment tensor inversion for the focal mechanism determination of events with a low seismicity. The East Baltic region (EBR) is the region with low seismicity. 2. The focal mechanism is determined using the data from a limited number of stations. Practical significance: The results of focal mechanism determination can be used to study seismicity for regions with a low seismicity using a limited number of stations.

An adaptive node-distribution method for radial-basis-function finite-difference modeling with optimal shape parameter

The radial-basis-function finite-difference (RBF-FD) method has been proven successful in modeling seismic-wave propagation. Node distribution is typically the first and most critical step in RBF-FD. Regarding the difficulties in seismic modeling, such as node distribution of complex geologic structures, we have designed an adaptive node-distribution method that can generate nodes automatically and flexibly as the computation proceeds with the adaptive grain-radius satisfied dispersion relation and stability condition of seismic-wave propagation. Our method consists of two novel points. The first one is that we adopt an adaptive grain-radius generation method, which can automatically provide a wider scope of grain radius in seismic modeling while satisfying the dispersion relation and stability condition; the second one is that the node-generation algorithm is built by a smoothed model, which significantly improves the modeling stability at a reduced computational cost. Excessive or undesirable shape parameters will create a very ill-conditioned problem. A set of optimal shape parameters for different numbers of neighbor nodes is found quantitatively by minimizing root-mean-square error functions. This optimization method enables us to achieve an improved meshfree modeling process with higher accuracy and practicability and fewer spurious diffractions caused by the transition of different sampling areas. Several numerical results verify the feasibility of our adaptive node-distribution method and the optimal shape parameters.

Imaging the Chicxulub Central Crater Zone from Large-Scale Seismic Acoustic Wave Propagation and Gravity Modeling

Large impact structures are characterized by peak ring and central uplifts with lateral/vertical mass transport during late formation stages. Here we investigate the Chicxulub crater, which has been surveyed by an array of marine, aerial and land-borne geophysical methods. Seismic reflection surveys in its central sector have shown lack of resolution, making it difficult to image the central uplift. We develop an integrated seismic and gravity model for the structural elements, imaging the central uplift and melt and breccia units. The 3D velocity model built from interpolation of seismic data is validated using perfectly matched layer seismic acoustic wave propagation modeling, optimized at grazing incidence using the shift in the frequency domain. Modeling shows that lack of illumination relates to seismic energy that remains trapped in an upper low-velocity zone corresponding to the carbonate sediments, upper melt/breccias and surrounding faulted blocks. After conversion of seismic velocities into a volume of density values, we apply parallel forward gravity modeling to constrain the size and shape of the central uplift, which has a ~ 40 km diameter concave upwards top lying at ~ 3.5–4.5 km depth. The preferred model provides a high-resolution image of crater units and structure. The gravity response of modeled units shows asymmetries in structure and the distribution of breccias, melt and target carbonates. Finally, we apply an adjoint reverse time migration approach for seismic imaging using the density and velocity models built for the acoustic wave propagation and gravity modeling, which allows improved modeling of the crater structure.

Subcell-resolution finite-difference modelling of seismic waves in Biot and JKD poroelastic media

SUMMARY We present a discrete representation of strongly heterogeneous poroelastic medium with the JKD-model of the frequency-dependent permeability and resistive friction, and the corresponding finite-difference (FD) scheme for numerical modelling of seismic wave propagation and earthquake ground motion in structurally complex media. The scheme is capable of subcell resolution, that is, allows for an arbitrary shape and position of an interface in the spatial grid. The medium can have either a zero resistive friction or non-zero constant resistive friction or JKD frequency-dependent resistive friction. The scheme has the same computational efficiency as the scheme for a smoothly and weakly heterogeneous medium (medium without material interfaces) because the number of operations for updating wavefield is the same. Several comparisons with a semi-analytical approach proves the efficiency and reliability of the subcell-resolution FD scheme. An illustrative example demonstrates differences between earthquake ground motion in the Biot's and JKD variants of the model of the surface sedimentary basin. The example indicates that it is desirable to perform an extensive parametric study in order to find out when it is necessary to apply relatively complicated and computationally more demanding JKD model and when much simpler Biot's model is sufficient.

Frequency-domain elastic-wave modeling for polygonal topography using rotated average-derivative difference operators

Modeling of seismic wave propagation in areas with irregular topography is an important topic in the field of seismic exploration. As a popular numerical method for seismic modeling, the finite difference method is nontrivial to consider the irregular free surface. There have been extensive studies on the time-domain finite difference simulations with irregular topography; however, the frequency-domain finite difference simulation considering irregular topography is relatively less studied. The average-derivative approach is an optimal numerical simulation scheme in the frequency domain, which can produce accurate modeling results at a relatively low computational cost. Nevertheless, this approach can only deal with the modeling problems with a flat free surface. To address this issue, we design a new frequency-domain finite difference scheme by introducing the polygonal representation of topography into the average-derivative method. The irregular topography is represented by line segments with various slopes. An extension of the conventional average-derivative difference operator in the local rotated coordinate system is used for formulating the spatial derivatives aligned with the topographic line segments. As a result, new average-derivative difference schemes are obtained for irregular topography. In this way, the average-derivative optimal method is generalized to the model with irregular topography. Numerical examples show the effectiveness of the presented method.

Rational approximation of P-wave kinematics — Part 2: Orhorhombic media

Orthorhombic anisotropy is a modern standard for 3D seismic studies in complex geologic settings. Several seismic data processing methods and wave propagation modeling algorithms in orthorhombic media rely on phase-velocity, group-velocity, and traveltime approximations. The algebraic simplicity of an approximate equation is an important factor in these media because the governing equations are more complicated than transversely isotropic media. To approximate the P-wave kinematics in acoustic orthorhombic media, we have developed a new 3D general functional equation that has a simple rational form. Using the general form, we adopt two versions of rational approximations for the phase velocity, group velocity, and traveltime. The first version uses a simpler functional form and parameter definition within the orthorhombic symmetry planes. The second version is more accurate, using one parameter that is defined out of the symmetry planes. For the phase velocity, we obtain another approximation that is no longer rational but is still algebraically simple, exact for 3D transversely isotropic media, and it is exact within the symmetry planes of orthorhombic media. We find superior accuracy in our approximations compared with previous ones, using numerical studies on multiple moderately anisotropic orthorhombic models. We investigate the effect of the negative anellipticity parameters on the accuracy and find that, in models in which the error of the existing most accurate approximations exceeds 2%, the error of the new approximations remains below 0.2%. The adopted approximations are algebraically simpler and stably more accurate than existing approximations; therefore, they may be considered as attractive alternatives for the existing approximations in many practical applications. We extend the applicability of our approximations by using them to obtain the equations of group direction as a function of phase direction and vice versa, which are useful in wave propagation modeling methods.

Dynamic Modelling of Seismic Wave Propagation due to a Remote Seismic Source: A Case Study

Rockbursts are some of the most dangerous phenomena encountered in underground mines. They can be triggered by an abrupt ground disturbance from a remote seismic source. In this paper, the mechanism of remotely triggered rockbursts around a coal mine drift is analyzed using the synchro-squeezing transform (SST) method and dynamic modeling. Based on seismic monitoring data, the seismic waveform at the source is estimated through empirical scaling law and calibrated in the model. Then, with the SST seismic waveform decomposition codes, the P and S waves are separated and filtered, and then applied in orthogonal oscillating directions at the source simulating their inherently diverse radiation mechanism. A three-dimensional numerical model was constructed with FLAC3D for a case study rockburst of a mine drift at Wudong Coal Mine. The static simulation results show that stress concentration occurs in the roof and rock pillar near the production level, and that the rock mass around the drift is already damaged due to the upper unloading effect and the bending deformation of the roof and pillar. Dynamic modeling confirmed that weak materials cause greater attenuation of particle vibration than hard materials. The rock mass around the drift experienced a significant strain energy release of 2.1 × 108 J and maximum displacements in the working face of 129 mm. The P wave showed a larger contribution to the dynamic disturbance in the horizontal direction than in vertical direction, while the S wave has a predominant proportion in the vertical direction. As expected, particle vibration velocity and displacement on the north side of the drift are greater than those on the south side—under incidents waves propagating from the north side. Based on model results and peak particle velocity (PPV), it was deemed that the computed released energy is large enough to threaten the drift integrity, resulting in drift damage for up to 200 m from the working face. The predicted failure characteristics and potential damage range are consistent with field observations. The proposed approach of synchro-squeezing transform (SST) and dynamic modeling could prove useful in the assessment of damage in rockburst-prone areas—it could further help assess the need for and design of dynamic rock supports.