A review of the influencing factors on teleseismic traveltime tomography

Teleseismic traveltime tomography is an important tool for investigating the crust and mantle structure of the Earth.The imaging quality of teleseismic traveltime tomography is affected by many factors,such as mantle heterogeneities,source uncertainties and random noise.Many previous studies have investigated these factors separately.An integral study of these factors is absent.To provide some guidelines for teleseismic traveltime tomography,we discussed four main influencing factors:the method for measuring relative traveltime differences,the presence of mantle heterogeneities outside the imaging domain,station spacing and uncertainties in teleseismic event hypocenters.Four conclusions can be drawn based on our analysis.(1)Comparing two methods,i.e.,measuring the traveltime difference between two adjacent stations(M1)and subtracting the average traveltime of all stations from the traveltime of one station(M2),reveals that both M1 and M2 can well image the main structures;while M1 is able to achieve a slightly higher resolution than M2;M2 has the advantage of imaging long wavelength structures.In practical teleseismic traveltime tomography,better tomography results can be achieved by a two-step inversion method.(2)Global mantle heterogeneities can cause large traveltime residuals(up to about 0.55 s),which leads to evident imaging artifacts.(3)The tomographic accuracy and resolution of M1 decrease with increasing station spacing when measuring the relative traveltime difference between two adjacent stations.(4)The traveltime anomalies caused by the source uncertainties are generally less than 0.2 s,and the impact of source uncertainties is negligible.

Three-dimensional frequency-domain full waveform inversion based on the nearly-analytic discrete method

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.

A nearly analytic exponential time difference method for solving 2D seismic wave equations

In this paper, we propose a nearly analytic exponential time difference （NETD） method for solving the 2D acoustic and elastic wave equations. In this method, we use the nearly analytic discrete operator to approximate the high-order spatial differential operators and transform the seismic wave equations into semi-discrete ordinary differential equations （ODEs）. Then, the converted ODE system is solved by the exponential time difference （ETD） method. We investigate the properties of NETD in detail, including the stability condition for 1-D and 2-D cases, the theoretical and relative errors, the numerical dispersion relation for the 2-D acoustic case, and the computational efficiency. In order to further validate the method, we apply it to simulating acoustic/elastic wave propagation in mul- tilayer models which have strong contrasts and complex heterogeneous media, e.g., the SEG model and the Mar- mousi model. From our theoretical analyses and numerical results, the NETD can suppress numerical dispersion effectively by using the displacement and gradient to approximate the high-order spatial derivatives. In addition, because NETD is based on the structure of the Lie group method which preserves the quantitative properties of differential equations, it can achieve more accurate results than the classical methods.

A holistic perspective on Earth system science

Earth system can be categorized into three parts, solid Earth system, surface Earth system, and Sun-Earth space system. These three subsystems not only have mutual transmission and coupling relationships in both energy and matter but also involve multiple scales from microscopic to macroscopic. Earth system science is characterized by its globality and unity with a holistic view and a systematic view at multiple scales in both space and time. It focuses not only on the physical, chemical and biological interactions between various geospheres but also on the properties, behaviors, processes, and mechanisms of the entire Earth and its spheres. Although significant progress has been made in the study of internal disciplines of these three subsystems,there is still insufficient understanding of their overall behavior and interactions between individuals, thus facing challenges of different types and levels. The solid Earth system is composed of the crust, mantle, and core. Existing observational techniques struggle to penetrate deep into the mantle, making direct observation and data acquisition difficult;the extreme environments within Earth, such as high temperature, high pressure, and strong magnetic fields, also pose great challenges to observational equipment and scientific experiments. The surface Earth system is an open complex mega-system, in which there are complex interactions and feedback mechanisms among its geospheres(such as atmosphere, hydrosphere, biosphere, pedosphere and lithosphere), leading to difficulties in understanding of its overall behavior and long-term evolution. Biological activities have become increasingly significant in affecting the surface Earth system. The coupling between the internal and external Earth systems becomes more complex. Distinguishing and quantifying the impacts of Earth spherical interactions and biological activities on the surface Earth system is a major challenge. The Sun-Earth space system involves multiple physical processes such as solar activity, Earth's magnetic field, atmosphere, and space weather. Solar activity significantly affects the Earth's space environment, but existing observational and reconstruction methods and prediction models still lack precision and timeliness.Thus it is important to improve the prediction capability of solar activity and reduce the impact of space weather disasters. How to cross different scales and establish coupled models of multiple physical processes is a significant challenge in the study of the Sun-Earth space system. Because the various processes and phenomena within and between these three Earth subsystems often span multiple scales in both space and time and exhibit strong nonlinear characteristics, understanding their behaviors and processes becomes complex and variable, posing great challenges for theoretical modelling and numerical simulation. Therefore,the study of Earth system science requires in-depth interdisciplinary integration to jointly reveal the basic laws and operating mechanisms of Earth system.

An eikonal equation-based earthquake location method by inversion of multiple phase arrivals

The precise determination of earthquake location is the fundamental basis in seismological community,and is crucial for analyzing seismic activity and performing seismic tomography.First arrivals are generally used to practically determine earthquake locations.However,first-arrival traveltimes are not sensitive to focal depths.Moreover,they cannot accurately constrain focal depths.To improve the accuracy,researchers have analyzed the depth phases of earthquake locations.The traveltimes of depth phases are sensitive to focal depths,and the joint inversion of depth phases and direct phases can be implemented to potentially obtain accurate earthquake locations.Generally,researchers can determine earthquake locations in layered models.Because layered models can only represent the first-order feature of subsurface structures,the advantages of joint inversion are not fully explored if layered models are used.To resolve the issue of current joint inversions,we use the traveltimes of three seismic phases to determine earthquake locations in heterogeneous models.The three seismic phases used in this study are the first P-,sPg-and PmP-waves.We calculate the traveltimes of the three seismic phases by solving an eikonal equation with an upwind difference scheme and use the traveltimes to determine earthquake locations.To verify the accuracy of the earthquake location method by the inversion of three seismic phases,we take the 2021 M_(S)6.4 Yangbi,Yunnan earthquake as an example and locate this earthquake using synthetic and real seismic data.Numerical tests demonstrate that the eikonal equation-based earthquake location method,which involves the inversion of multiple phase arrivals,can effectively improve earthquake location accuracy.

熔/流体迁移与2022年泸定Ms 6.8地震

2022年9月5日,鲜水河断裂带东南末端发生泸定Ms 6.8级地震,造成严重的人员伤亡和财产损失.为深入探究该地震构造背景和孕震机理,利用全波形反演方法获得了震源区及周边高精度地震波速度、径向各向异性和泊松比结构.研究表明,泸定地震发生在低波速异常区边缘,在其下方存在显著的低S波速、负径向各向异性(V_(SH)<V_(SV))和高泊松比异常区,此热-流变结构异常区的形成可能与深部地壳和上地幔流体(包括部分熔融)有关,这些熔/流体沿倾斜的Moho面和/或壳内断裂向上迁移,引起中上地壳孕震层的弱化和应力集中,进而导致地震发生,表明深部熔/流体迁移在泸定地震孕育和破裂过程中发挥着重要作用.

Receiver function imaging of dense seismic array and deep dynamic mechanism beneath the eastern South China

The Mesozoic tectonic framework of the eastern South China is mainly controlled by subduction,turning toward,and rollback of the Pacific Plate.Recent studies of receiver function imaging and ambient noise tomography have revealed the“Yshaped”thinnest crustal belt in the eastern South China under the overall extension of the lithosphere.However,the deep dynamic environment and formation mechanisms of the thin crustal belt remain debatable.Here we obtained high-resolution images of the crustal thickness and Poisson’s ratio in the eastern South China Block applying the recently proposed H-κ-c receiver function method,using data recorded by 305 dense portable broadband stations and 219 permanent stations surrounding.Additionally,we discussed the deep dynamic formation mechanism of the“Y-shaped”thinnest crustal belt coupled with two common conversion point stacked images at key locations.Results show that the average crustal thickness of the study area is 33 km(thin crust)and the average Poisson’s ratio is 0.24(low ratio).The overall crustal thinning toward the continental margin is likely because eastern South China was in a back-arc extension environment,which was induced by the rollback of the subducted plate in the Early Cretaceous.The crustal thickness of the“Y-shaped”thinnest crustal belt is<30 km,which is 3-5 km thinner than that outside the zone.The eastern branch is distributed along the trajectory of Nanchang-Ji’an-Ganzhou-Shaoguan-Guangzhou,and the western branch is around the Jianghan-Xiangzhong Basin,both of which intersect in Nanling.The eastern branch of the thin crustal zone indicates the potential location of the Pacific subduction slab breakoff,and the formation mechanism may be related to the interaction of deep-shallow processes,including the upwelling of mantle heat flow through the slab window and transtensional pre-existing faults.We developed a dynamic model that combines subduction-breakoff-rollback processes of the Paleo-Pacific Plate and accompanying deep fluid upwelling to explain the regional extension of the South China lithosphere,the formation mechanism of the thinnest crustal belt,and the distribution of granitic plutons.

THE WASSERSTEIN-FISHER-RAO METRIC FOR WAVEFORM BASED EARTHQUAKE LOCATION

In this paper,we apply the Wasserstein-Fisher-Rao(WFR)metric from the unbalanced optimal transport theory to the earthquake location problem.Compared with the quadratic Wasserstein(W2)metric from the classical optimal transport theory,the advantage of this method is that it retains the important amplitude information as a new constraint,which avoids the problem of the degeneration of the optimization objective function near the real earthquake hypocenter and origin time.As a result,the deviation of the global minimum of the optimization objective function based on the WFR metric from the true solution can be much smaller than the results based on the W2 metric when there exists strong data noise.Thus,we develop an accurate earthquake location method under strong data noise.Many numerical experiments verify our conclusions.

Full waveform inversion based on the ensemble Kalman filter method using uniform sampling without replacement

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 highresolution results and a wider convergence domain. Accordingly, the new method overcomes the disadvantage of conventional FWIs that depend strongly on the initial model.

A unified poroviscoelastic model with mesoscopic and microscopic heterogeneities

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.

Velocity structure and radial anisotropy beneath the northeastern Tibetan Plateau revealed by eikonal equation-based teleseismic P-wave traveltime tomography

The northeastern Tibetan Plateau serves as the frontier for the northeastward expansion of the plateau.In this area,the Tibetan Plateau interacts with the surrounding blocks,such as the Alxa Block,the Ordos Block,the Kunlun-West Qinling belt and the Sichuan Basin.Because of this expansion and interaction,this area suffers from intense deformation.At present,the evolution and deformation mechanisms of the northeastern Tibetan Plateau remain controversial.To provide new insights into these mechanisms,in this study,we conduct tomography of the P-wave velocity and radial anisotropy structures beneath the northeastern Tibetan Plateau.We choose a total of 667 teleseismic earthquakes from August 2006 to October 2020.Waveforms of these earthquakes were recorded by 921 broadband seismic stations in the northeastern Tibetan Plateau and surrounding areas.We first perform cross-correlation on waveforms of each station pair and obtain 770,749 P-wave traveltime differences.Then,we invert the differential traveltime data by applying eikonal equation-based teleseismic tomography.Finally,the P-wave velocity and radial anisotropy structures at depths from 30 to 800 km below the northeastern Tibetan Plateau are obtained.Our tomographic model shows clear low-velocity anomalies and positive radial anisotropy in the lower crust under the northeastern Qilian orogen,the northeastern Songpan-Ganzi belt and the western Qinling fold zone.These features are integrated to demonstrate the existence of lower crustal flow in the study area.Prominent low-velocity anomalies and positive radial anisotropy are found in the uppermost mantle beneath the Qilian orogen,the northeastern Songpan-Ganzi belt and western Qinling fold zone.These characteristics are combined to infer a weak lithosphere and horizontal asthenospheric flow under these tectonic units.Both the Ordos Block and the Sichuan Basin exhibit clear high-velocity anomalies and negative radial anisotropy in the uppermost mantle,thereby reflecting the high mechanical strength of the lithosphere beneath these blocks.High-velocity anomalies are also present in the upper mantle under the northern Chuandian block,potentially implying the northward subduction of the Indian plate.Furthermore,the front of the subducted Indian plate is imaged close to the Xianshuihe fault rather than the Kunlun fault.

Quasi-waveform seismic tomography of crustal structures in the capital circle region of China

Seismic tomography is one of the main tools to explore the interior of the earth.In this study,the quasi-waveform seismic tomographic method is used for the first time to reveal the crustal structures in the capital region of China.3-D highresolution V_P,V_S and the Poisson’s ratio models are generated by inverting 29839 direct P-and 29972 direct S-wave traveltimes selected from 3231 local earthquakes.The results reveal strong crustal heterogeneities.The velocity anomalies at shallow depths are well consistent with surface geologic structures.The relatively low-velocity anomaly layer in the middle crust may be the result of multiple phases of tectonic activity.Earthquakes generally occurred on the boundaries of high-and low-velocity and Poisson’s ratio anomalies.There are obvious low-velocity anomalies below the hypocenters of the Tangshan earthquake and the historical Sanhe-Pinggu earthquake,implying the existence of fluids.The similar velocity structures around the hypocenters of the two earthquakes indicate that the occurrences of the two earthquakes may be related to the same mechanism.The highresolution velocity models provide important observational constraints on the small-scale heterogeneities and dynamic mechanism of the crust in the capital region of China.

Geometry-preserving full-waveform tomography and its application in the Longmen Shan area

Classic L2-norm-based waveform tomography is often plagued by insurmountable cycle skipping problems;as a result,the iterative inversion falls into local minima,yielding erroneous images.According to the optimal transportation theory,we adopt a novel geometry-preserving misfit function based on the quadratic Wasserstein metric(W2-norm),which improves the stability and convexity of the inverse problem.Numerical experiments illustrate that W2-norm-based full-waveform tomography has a larger convergence radius and a faster convergence rate than the L2-norm and can effectively mitigate cycle skipping issues.We apply this method to the Longmen Shan area and obtain a reliable lithospheric velocity model.Our tomographic results indicate that the crystalline crust underlying the Sichuan Basin wedges into the crustal interior of the Tibetan Plateau,and the mid-lower crust of the eastern Tibetan Plateau is characterized by low shear-wave velocities,indicating that ductile crustal flow and strong interactions between terranes jointly dominate the uplift behavior of the Longmen Shan.Furthermore,we find that large earthquakes(e.g.,the Wenchuan and Lushan events)occur not only at the junction between high-and low-velocity regions but also in the transition zone from positive to negative radial anisotropy.These findings improve our understanding of the mechanism responsible for large earthquakes in this region.

Reanimating lithosphere of North China: new insights from full-waveform seismic tomography

The North China Craton(NCC)consists of the western Ordos block and the eastern block of the North China Basin,which are separated by the Taihang Mountains.The eastern NCC experienced multiple stages of tectonic-magmatic activity and internal deformation of the crust in the Mesozoic and Cenozoic[1],while the Ordos block retains the original cratonic characteristics,i.e.,it is stable and lacks seismic activity.Continuous subduction of the western Pacific Plate is responsible for lithosphere reanimating in North China[2],but the detailed morphology of lithospheric thinning patterns remains unresolved.Because rock wave velocity is affected by temperature,and rock anisotropy is mainly dominated by directional alignment of crystals,seismic tomography can be applied to depict the shapes of cooler lithosphere and warmer mantle flow,providing crucial constraints on crust and mantle deformation.

Positions of large earthquakes revealed from velocity heterogeneities and radial anisotropy in the eastern Tibetan Plateau

As the largest and highest plateau on the Earth,the Tibetan Pla-teau is the most important natural laboratory we have for under-standing continent–continent collision and mountain building;the Tibetan Plateau is characterized by conspicuous crustal thick-ening,considerable deformation,complex metamorphism,and sig-nificant uplift[1].Because of the resistance of the rigid Tarim lithosphere in the west,the eastern Tibetan Plateau has been the main region for the occurrence of lateral motions of escaped pla-teau materials caused by the continued penetration of India into Eurasia[2].To compensate for the compression of ongoing conver-gence,the eastward motion of crustal blocks is accommodated by a series of active strike-slip faults(e.g.,Xianshuihe,Kunlun,and AilaoShan-Red River).Based on geologic and geodetic observations,the hypothesis of lower crustal flow has been proposed to reconcile the paradox between limited young upper crustal shortening and anomalously high topography across the eastern plateau margin[3].

A wave propagation model with the Biot and the fractional viscoelastic mechanisms

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.

A Weighted Runge-Kutta Method with Weak Numerical Dispersion for Solving Wave Equations

In this paper, we propose a weighted Runge-Kutta (WRK) method to solvethe 2D acoustic and elastic wave equations. This method successfully suppresses thenumerical dispersion resulted from discretizing the wave equations. In this method,the partial differential wave equation is first transformed into a system of ordinarydifferential equations (ODEs), then a third-order Runge-Kutta method is proposedto solve the ODEs. Like the conventional third-order RK scheme, this new methodincludes three stages. By introducing a weight to estimate the displacement and itsgradients in every stage, we obtain a weighted RK (WRK) method. In this paper, weinvestigate the theoretical properties of the WRK method, including the stability criteria, numerical error, and the numerical dispersion in solving the 1D and 2D scalarwave equations. We also compare it against other methods such as the high-ordercompact or so-called Lax-Wendroff correction (LWC) and the staggered-grid schemes.To validate the efficiency and accuracy of the method, we simulate wave fields in the2D homogeneous transversely isotropic and heterogeneous isotropic media. We conclude that the WRK method can effectively suppress numerical dispersions and sourcenoises caused in using coarse grids and can further improve the original RK methodin terms of the numerical dispersion and stability condition.

A New Earthquake Location Method Based on the Waveform Inversion

In this paper,a new earthquake location method based on the waveform inversion is proposed.As is known to all,the waveform misfit function under the L2 measure is suffering from the cycle skipping problem.This leads to a very small convergence domain of the conventional waveform based earthquake location methods.In present study,by introducing and solving two simple sub-optimization problems,we greatly expand the convergence domain of the waveform based earthquake location method.According to a large number of numerical experiments,the new method expands the range of convergence by several tens of times.This allows us to locate the earthquake accurately even from some relatively bad initial values.

A Weighted Runge-Kutta Discontinuous Galerkin Method for 3D Acoustic and Elastic Wave-Field Modeling

Numerically solving 3D seismic wave equations is a key requirement for forward modeling and inversion.Here,we propose a weighted Runge-Kutta dis-continuous Galerkin(WRKDG)method for 3D acoustic and elastic wave-field mod-eling.For this method,the second-order seismic wave equations in 3D heteroge-neous anisotropic media are transformed into a first-order hyperbolic system,and then we use a discontinuous Galerkin(DG)solver based on numerical-flux formulations for spatial discretization.The time discretization is based on an implicit di-agonal Runge-Kutta(RK)method and an explicit iterative technique,which avoids solving a large-scale system of linear equations.In the iterative process,we introduce a weighting factor.We investigate the numerical stability criteria of the 3D method in detail for linear and quadratic spatial basis functions.We also present a 3D analysis of numerical dispersion for the full discrete approximation of acoustic equation,which demonstrates that the WRKDG method can efficiently suppress numerical dispersion on coarse grids.Numerical results for several different 3D models including homogeneous and heterogeneous media with isotropic and anisotropic cases show that the 3D WRKDG method can effectively suppress numerical dispersion and provide accurate wave-field information on coarse mesh.

A Nearly Analytical Discrete Method for Wave-Field Simulations in 2D Porous Media

The nearly analytic discrete method(NADM)is a perturbation method originally proposed by Yang et al.(2003)[26]for acoustic and elastic waves in elastic media.This method is based on a truncated Taylor series expansion and interpolation approximations and it can suppress effectively numerical dispersions caused by the discretizating the wave equations when too-coarse grids are used.In the present work,we apply the NADM to simulating acoustic and elastic wave propagations in 2D porous media.Our method enables wave propagation to be simulated in 2D porous isotropic and anisotropic media.Numerical experiments show that the error of the NADM for the porous case is less than those of the conventional finite-difference method(FDM)and the so-called Lax-Wendroff correction(LWC)schemes.The three-component seismic wave fields in the 2D porous isotropic medium are simulated and compared with those obtained by using the LWC method and exact solutions.Several characteristics of wave propagating in porous anisotropic media,computed by the NADM,are also reported in this study.Promising numerical results illustrate that the NADM provides a useful tool for large-scale porous problems and it can suppress effectively numerical dispersions.