Search

In this paper, we derive a new poroelastic wave equation in triple-porosity media and develop a weighted Runge-Kutta (RK) discontinuous Galerkin method (DGM) for solving it. Based on Biot’s theory and Lagrangian formulas, we obtain 3D Biot’s equations in a heterogeneous anisotropic triple-porosity medium. We also summarize poroelastic wave equations in single- and double-porosity media. The traditional two-phase theory is a special case of the porous theory. Compared with single- or dual-porosity wave equations, our triple-porosity wave equation generates more accurate wavefield information. The isotropic and anisotropic cases are considered. Subsequently, we formulate the new poroelastic equation into a first-order hyperbolic conservation system, which is suitable to be solved by DGM. An optimized local Lax-Friedrichs flux and an implicit-weighted RK time discretization scheme are used for this computation. We use two types of mesh elements — quadrilateral and unstructured triangular elements. We find that there are two kinds of slow P waves, P1 and P2 waves, in double-porosity media, whereas three kinds of slow P waves, P1, P2, and P3 waves, exist in three-porosity media. We also study the analytical and numerical solutions of propagation velocities for different waves in isotropic media without dissipation using the Jacobian matrix of DGM and provide a comparison of field variables about three types of wave equations. Finally, we conduct a series of examples to quantitatively investigate the propagation properties of seismic waves in isotropic and anisotropic multiporosity media computed by DGM. The slow P wave in multiporosity media with dissipation decays rapidly, which also will lead to phase distortion. Numerical results verify the correctness and applicability of our proposed new equation and indicate that the weighted RK DGM is a stable and accurate algorithm to simulate wave propagation in poroelastic media.

SUMMARYThe Qinling belt is a transitional zone lying among three units: the North China block (NCB), the South China block (SCB) and the northeastern Tibetan Plateau (NETP). Owing to the interaction of these units, complex deformation has occurred in the Qinling belt. Although many studies have been conducted to understand the deformation mechanism in the Qinling belt, some key issues are still under debate, such as whether middle-lower crustal flow exists beneath the western Qinling belt (WQB). High-resolution images of subsurface structures are essential to shed light on the deformation mechanism. In this paper, high-resolution images of the velocity structure and azimuthal anisotropy beneath the Qinling belt are obtained by using an eikonal equation-based traveltime tomography method. Our seismic tomography inverts 38 719 high-quality P-wave first arrivals from 1697 regional earthquakes recorded by 387 broad-band seismic stations. In the WQB, our tomography results show low-velocity anomalies but relatively weak anisotropy in the middle-lower crust. These features suggest that middle-lower crustal flow may not exist in this area. In the central Qinling belt (CQB), we find low-velocity anomalies in the middle-lower crust; however, the fast velocity directions no longer trend E–W but vary from NNE–SSW to N–S. These characteristics can be ascribed to the convergence and collision between the NCB and the SCB. In addition, we find strong low-velocity anomalies in the uppermost mantle beneath the CQB, which may indicate delamination of the lower crust. In the southern Qinling belt, we observe significant high-velocity anomalies in the upper crust beneath the Hannan–Micang and Shennong–Huangling domes. These high-velocity anomalies indicate a mechanically strong upper crust, which is responsible for the arc-shaped deformation process of the Dabashan fold. Based on the P-wave velocity and azimuthal anisotropic structures revealed by the inversion of high-quality seismic data, the deformation of the Qinling belt is affected mainly by the convergence between the NCB and the SCB rather than by the middle-lower crustal flow from the Tibetan Plateau.

We use seismic waveform adjoint tomography to constrain the shear-wave attenuation models of the crust and upper mantle in Sichuan–Yunnan region of China. On the base of 3D velocity model with high accuracy by the previous work, we use the adjoint approach to efficiently construct the anelastic structures. Spectral element method with graphic processing unit acceleration is implemented in our work. We use an envelope-based misfit function and develop a mini-batch gradient descent algorithm for model update. We have applied the adjoint tomography algorithms to 41 seismic events, including 1911 high-quality three-component displacement seismic records from 2009 and 2017 in Sichuan–Yunnan region. Synthetic tests show that the attenuation model is well resolved. The generated shear-wave attenuation model reveals detailed structural characteristics of the upper mantle in Sichuan–Yunnan region. Some notable features are observed, such as an obvious strong-attenuation zone in northern Yunnan, which provides evidence for the existence of high-attenuation middle and lower crustal channels.

Tight porous media have extremely low porosities and permeabilities, which can induce low-fluid mobility and a significant non-Darcy flow phenomenon. We build a unified mathematical model considering the non-Darcy seepage and the viscoelasticity mechanism, which includes two new parameters related to the non-Darcy flow mechanism. We call this the non-Darcy viscoelastic (NDV) model. Numerical experiments reveal that the two parameters have significant effects on the dispersion and attenuation of the P and S waves. We find that the dispersion and attenuation of the waves depend on the initial amplitude. A comparison is conducted between experimental data for tight porous media in the low (seismic) frequency band and the numerical results of the NDV model. The comparison demonstrates that the NDV model can predict the velocity dispersion in tight oil and gas reservoirs more accurately.

It is of great importance to predict reservoir petrophysical properties, such as porosity and saturation based on elastic properties, including density and wave propagation velocity, through quantitative seismic interpretation. A joint two-step inversion (JoiTSI) method is proposed for predicting porosity and saturation with a forward petrophysical model. As the forward model, the Biot/squirt model with multiphase fluids (gas, oil, and water) describes the nonlinear relationships between petrophysical and elastic properties. The first step of the two-step inversion is to invert porosity and saturation via a semianalytical approximation approach. The second step is a joint inversion at high porosities where porosity and saturation are simultaneously optimized using a hybrid genetic algorithm. In these steps, maximum a posteriori estimations are performed based on normal prior distributions, where the posterior distributions may differ between two steps because of the potential variability in reservoir properties. The numerical tests indicate that this JoiTSI method is more efficient and more accurate than the traditional inversion methods. In addition, Poisson’s ratio, P-wave impedance, and S-wave impedance are found to give better inversion results than those based on more commonly used attributes (density, P-wave impedance, and S-wave impedance). In the oilfield seismic example, the inversion accuracy of porosity and saturation is verified by actual well logs, and the major oil and water layers are correctly identified.

SUMMARY The Tengchong–Baoshan region (TBR) is located at the margin of the southeastern Tibetan Plateau (SETP) and is subject to the complex interaction among multiple tectonic blocks. Therefore, the TBR is approximate for investigating the mechanism responsible for the interaction among the SETP and its surrounding blocks. Based on eikonal equation-based seismic tomography, we obtain high-resolution P-wave velocity and azimuthal anisotropy at depths of 0–40 km beneath the TBR. We discover a large magma reservoir (with horizontal scale of ∼50 km and vertical scale of ∼15 km) in the middle-lower crust beneath the Tengchong volcanic zone, which is directly fed by hot materials from the uppermost mantle. We also detect a shallow magma chamber in the uppermost crust above the large reservoir. These magma reservoirs at different depths constitute a multiscale magma system. Furthermore, we observe a low-velocity belt beneath the Weixi–Qiaohou–Weishan and Nantinghe faults and postulate that it corresponds to lower crustal flow and hot mantle upwelling. Additionally, we infer that two M > 7.0 earthquakes that have occurred in Baoshan block nucleated in the transition zone between low- and high-velocity anomalies. Based on the azimuthal anisotropy tomography results, we consider that incoherent deformation in the crust causes the uneven accumulation of stress, which facilitated the occurrence of these two events. The 2021 Ms 6.4 Yangbi earthquake also occurred in the transition zone between low- and high-velocity anomalies. Based on the P-wave velocity tomography results, we suggest that hot material decreases the fault strength and leads to the instability of the fault.

The discontinuous Galerkin (DG) method is a numerical algorithm that is widely used in various fields. It has high accuracy and low numerical dispersion with advantages of easy handling boundary conditions and good parallel performance. In this study, we develop an efficient parallel weighted Runge–Kutta discontinuous Galerkin (WRKDG) method on unstructured meshes for solving 3D seismic wave equations. The DG method we use is based on the first-order formulation of a hyperbolic system with an explicit weighted Runge–Kutta time discretization. We describe the numerical scheme and parallel implementation in detail, and analyze the stability condition and numerical dispersion and dissipation. We carry out a convergence test on unstructured meshes and investigate the parallel efficiency of the implementation of the WRKDG method. The speedup curve shows that the method has good parallel performance. Finally, we present several numerical simulation examples, including acoustic and elastic wave propagations in isotropic and anisotropic media. Numerical results further verify the effectiveness of the WRKDG method in solving 3D wave propagation problems.

Although teleseismic waveform tomography can provide high-resolution images of the deep mantle, it is still unrealistic to numerically simulate the whole domain of seismic wave propagation due to the huge amount of computation. In this article, we develop a new three-dimensional hybrid method to address this issue, which couples the modified frequency–wavenumber (FK) method with the 3D time–space optimized symplectic (TSOS) method. First, the FK method, which is used to calculate the semianalytical incident wavefields in the layered reference model, is modified to compute the wavefields efficiently with a significantly low-memory requirement. Second, 3D TSOS method is developed to model the seismic wave propagating in the local 3D heterogeneous domain. The low memory requirement of the modified FK method and the high accuracy of the TSOS method make it feasible to obtain highly accurate synthetic seismograms efficiently. A crust–upper mantle model for P-, SV-, and SH-wave incidences is calculated to benchmark the accuracy and efficiency of the 3D optimized FK-TSOS method. Numerical experiments for 3D models with heterogeneities, undulated discontinuous interfaces, and realistic model in eastern Tibet, illustrate the capability of hybrid method to accurately capture the scattered waves caused by heterogeneities in 3D medium. The 3D optimized FK-TSOS method developed shows low-memory requirement, high accuracy, and high efficiency, which makes it be a promising forward method to further apply to high-resolution mantle structure images beneath seismic array.

We have developed a Runge-Kutta discontinuous Galerkin (RKDG) method for solving wave equations in isotropic and anisotropic poroelastic media at low frequencies. First, the 2D Biot’s two-phase equations are transformed into a first-order system with dissipation. Then, the system is discretized by using the discontinuous Galerkin method with a third-order Runge-Kutta time discretization. The numerical stability conditions for solving porous equations are also investigated. We test several examples to validate our method in isotropic and anisotropic poroelastic media. Comparisons of seismic responses with the finite-difference method on fine grids show the correctness of this method. Moreover, the numerical results indicate that the RKDG method can provide clear fast P-, slow P-, and S-waves for anisotropic poroelastic media on coarse meshes. Also, a two-layer porous model, a poroelastic-elastic model with horizontal interface, and an isotropic-anisotropic poroelastic model with sinusoidal interface demonstrate that our method can deal with complex wave propagation. Therefore, the simulation results show that RKDG is an accurate and stable method for solving Biot’s equations.

The nearly analytic discrete (NAD) method is a kind of finite difference method with advantages of high accuracy and stability. Previous studies have investigated the NAD method for simulating wave propagation in the time-domain. This study applies the NAD method to solving three-dimensional (3D) acoustic wave equations in the frequency-domain. This forward modeling approach is then used as the “engine” for implementing 3D frequency-domain full waveform inversion (FWI). In the numerical modeling experiments, synthetic examples are first given to show the superiority of the NAD method in forward modeling compared with traditional finite difference methods. Synthetic 3D frequency-domain FWI experiments are then carried out to examine the effectiveness of the proposed methods. The inversion results show that the NAD method is more suitable than traditional methods, in terms of computational cost and stability, for 3D frequency-domain FWI, and represents an effective approach for inversion of subsurface model structures.

Reverse time migration (RTM) is widely used in the industry because of its ability to handle complex geologic models including steeply dipping interfaces. The quality of images produced by RTM is significantly influenced by the performance of the numerical methods used to simulate the wavefields. Recently, a weighted Runge-Kutta discontinuous Galerkin (WRKDG) method has been developed to solve the wave equation, which is stable, explicit, and efficient in parallelization and suppressing numerical dispersion. By incorporating two different weights for the time discretization, we have obtained a more stable method with a larger time sampling. We apply this numerical method to RTM to handle complex topography and improve imaging quality. By comparing it to the high-order Lax-Wendroff correction method, we determine that WRKDG is efficient in RTM. From the results of the Sigsbee2B data, we can find that our method is efficient in suppressing artifacts and can produce images of good quality when coarse meshes are used. The RTM results of the Canadian Foothills model also demonstrate its ability in handling complex geometry and rugged topography.

Energy loss in porous media containing fluids is typically caused by a variety of dynamic mechanisms. In the Biot theory, energy loss only includes the frictional dissipation between the solid phase and the fluid phase, resulting in underestimation of the dispersion and attenuation of the waves in the low frequency range. To develop a dynamic model that can predict the high dispersion and strong attenuation of waves at the seismic band, we introduce viscoelasticity into the Biot model and use fractional derivatives to describe the viscoelastic mechanism, and finally propose a new wave propagation model. Unlike the Biot model, the proposed model includes the intrinsic dissipation of the solid frame. We investigate the effects of the fractional order parameters on the dispersion and attenuation of the P- and S-waves using several numerical experiments. Furthermore, we use several groups of experimental data from different fluid-saturated rocks to testify the validity of the new model. The results demonstrate that the new model provides more accurate predictions of high dispersion and strong attenuation of different waves in the low frequency range.

Northeastern Tibet is still in the primary stage of tectonic deformation and is the key area for studying the lateral expansion of the Tibetan plateau. In particular, the existence of lower crustal flow, southward subduction of the Asian lithosphere, and northward subduction of the Indian lithosphere beneath northeastern Tibet remains controversial. To provide insights into these issues, a high-resolution 3D radially anisotropic model of the lithospheric structure of northeastern Tibet is developed based on adjoint tomography. The Tibetan plateau is characterized as a low S-wave velocity lithosphere, in contrast with the relatively high S-wave velocities of the stable Asian blocks. Our tomographic result indicates that the low-velocity zone (LVZ) within the deep crust extends northeastward from Songpan–Ganzi to Qilian, which is interpreted as a channel flow within the crust. The upper mantle of Alxa and Qinling–Qilian are dominated by a rather homogeneous LVZ, which is inconsistent with the hypothesis that the Asian lithospheric mantle is being subducted southward beneath northeastern Tibet. Furthermore, high-velocity regions are observed in the southern Songpan–Ganzi region at depths ranging from 100 to 200 km, indicating that the northward-subducting Indian plate has probably reached the Xianshuihe fault.

We have developed a new discontinuous Galerkin (DG) method to solve the 2D seismic wave equations in isotropic and anisotropic media. This method uses a modified numerical flux that is based on a linear combination of the Godunov and the centered fluxes. A weighting factor is introduced in this modified numerical flux that is expected to be optimized to some extent. Through the investigations on the considerations of numerical stability, numerical dispersion, and dissipation errors, we develop a possible choice of optimal weighting factor. Several numerical experiments confirm the effectiveness of the proposed method. We evaluate a convergence test based on cosine wave propagation without the source term, which shows that the numerical errors in the modified flux-based DG method and the Godunov-flux-based method are quite similar. However, the improved computational efficiency of the modified flux over the Godunov flux can be demonstrated only at a small sampling rate. Then, we apply the proposed method to simulate the wavefields in acoustic, elastic, and anisotropic media. The numerical results show that the modified DG method produces small numerical dispersion and obtains results in good agreement with the reference solutions. Numerical wavefield simulations of the Marmousi model show that the proposed method also is suitable for the heterogeneous case.

The dispersive and dissipative properties of numerical methods are important for numerical modeling. We have evaluated a numerical dispersion-dissipation analysis for two discontinuous Galerkin methods (DGMs) — the flux-based DGM (FDGM) and the interior penalty DGM (IP DGM) for scalar wave equation. The semidiscrete analysis based on the plane-wave analysis is conducted for quadrilateral and triangular elements. Two kinds of triangular elements are taken into account. The fully discrete analysis for each method is conducted by incorporating a classic third-order total variation diminishing (TVD) Runge-Kutta (RK) time discretization. Our results indicate that FDGM produces smaller numerical dispersion than IP DGM, but it introduces more numerical dissipation. Notably, the two methods have different local convergence orders for numerical dispersion and dissipation. The anisotropy properties for different mesh types can also be identified. Several numerical experiments are carried out that verify some theoretical findings. The experiments exhibit that the numerical error introduced by FDGM is less than that introduced by IP DGM whereas the storage and calculation time of FDGM are greater than that of IP DGM. Overall, our work indicates that when both methods adopt third TVD RK time discretization, the computational efficiency of FDGM is slightly larger than IP DGM.

Input files for the experiments (examples of modes I, II, III, IV) performed in the paper "How do pre-existing weak zones and rheological layering of the continental lithosphere influence the development and evolution of intra-continental subduction?" using the codeASPECT.

The North China craton (NCC) was dominated by tectonic extension from late Cretaceous to Cenozoic, yet seismic studies on the relationship between crust extension and lithospheric mantle deformation are scarce. Here we present a three dimensional radially anisotropic model of NCC derived from adjoint traveltime tomography to address this issue. We find a prominent low S-wave velocity anomaly at lithospheric mantle depths beneath the Taihang Mountains, which extends eastward with a gradually decreasing amplitude. The horizontally elongated low-velocity anomaly is also featured by a distinctive positive radial anisotropy (VSH > VSV). Combining geodetic and other seismic measurements, we speculate the presence of a horizontal mantle flow beneath central and eastern NCC, which led to the extension of the overlying crust. We suggest that the rollback of Western Pacific slab likely played a pivotal role in generating the horizontal mantle flow at lithospheric depth beneath the central and eastern NCC. Published version . Tis work was supported by the Joint Earthquake Research Program of the National Natural Science Foundation and the China Earthquake Administration (No. U1839206) and the National R&D Program on Monitoring, Early Warning and Prevention of Major Natural Disaster, China (Grant No. 2017YFC1500301).

The wave-induced fluid flow (WIFF) is considered to be the main cause of dispersion and attenuation of seismic waves in fluid-saturated porous media. Among numerous theories, the mesoscopic and microscopic heterogeneities are considered to be the primary mechanisms causing the WIFF. Furthermore, in most rocks, the mesoscopic and microscopic heterogeneities exist simultaneously and can cause obvious transitions of the fast P-wave velocity, which means it is necessary to consider the influence of the two mechanisms on the dispersion and attenuation simultaneously. Numerous results have shown that the dispersions and attenuations caused by these two mechanisms can be approximated in terms of the Zener model. To combine the two mechanisms into a unified model, we introduce a new generalized Zener model into the Biot poroelasticity theory to obtain a new poroviscoelastic model. Comparisons between the numerical results and two groups of experimental data further confirm the validity of our new model.

The perfectly matched layer (PML) is an efficient artificial boundary condition that has been routinely implemented in seismic wave modeling. However, the effective combination of PML with symplectic numerical schemes for solving seismic wave equations has rarely been studied. In a companion paper, we have developed a complex-frequency-shifted convolutional PML (CPML) with a nonconstant compression grid parameter for solving the time-domain second-order seismic wave equation. Subsequently, we combine this CPML with two classes of symplectic methods to formulate symplectic partitioned Runge-Kutta (SPRK) + CPML and nearly analytic SPRK (NSPRK) + CPML, both of which are properly synchronized. To further investigate their validity, the two algorithms are then applied to acoustic and elastic wave simulations in typical geologic models, including a heterogeneous acoustic model, several isotropic and orthotropic elastic models, and an isotropic elastic model with a free-surface boundary. Relevant numerical results demonstrate the effectiveness of our CPML and combination algorithms. Specifically, the numerical accuracy and stability of the CPML that we develop are greatly improved compared with the classic split-field PML. Moreover, the final model with the free-surface boundary condition indicates that the nonconstant grid-compression parameter can eliminate the unstable modes at the free surface in the PML domain. The (N)SPRK + CPML that we propose is prospective for future application in other complex models and wave-equation-based migration and inversion.

To improve the computational accuracy and efficiency of long-time wavefield simulations, we have developed a so-called symplectic interior penalty discontinuous Galerkin (IPDG) method for 2D acoustic equation. For the symplectic IPDG method, the scalar wave equation is first transformed into a Hamiltonian system. Then, the high-order IPDG formulations are introduced for spatial discretization because of their high accuracy and ease of dealing with computational domains with complex boundaries. The time integration is performed using an explicit third-order symplectic partitioned Runge-Kutta scheme so that it preserves the Hamiltonian structure of the wave equation in long-term simulations. Consequently, the symplectic IPDG method combines the advantages of discontinuous Galerkin method and the symplectic time integration. We investigate the properties of the method in detail for high-order spatial basis functions, including the stability criteria, numerical dispersion and dissipation relationships, and numerical errors. The analyses indicate that the symplectic SIPG method is nondissipative and retains low numerical dispersion. We also find that different symplectic IPDG methods have different convergence behaviors. It is indicated that using coarse meshes with a high-order method produces smaller errors and retains high accuracy. We have applied our method to simulate the scalar wavefields for different models, including layered models, a rough topography model, and the Marmousi model. The numerical results show that the symplectic IPDG method can suppress numerical dispersion effectively and provide accurate information on the wavefields. We also conduct a long-term experiment that verifies the capability of symplectic IPDG method for long-time simulations.

Full waveform inversion (FWI) has been increasingly more and more important in seismology to better understand the interior structure of the Earth. FWI, by taking advantage of both the traveltime and amplitude in the data, provides high-resolution model parameters of the earth which can produce images with high resolution. However, this inversion method conventionally suffers from non-uniqueness due to many local minima of the objective function and large computing costs. In this study, we propose a new FWI method in a semi-random framework by integrating the ensemble Kalman filter and uniform sampling without replacement. Numerical results demonstrate that the new method can achieve high-resolution results and a wider convergence domain. Accordingly, the new method overcomes the disadvantage of conventional FWIs that depend strongly on the initial model.