Geophysical Journal International

SummaryCurrently, gravimetric forward modelling of mass density structures with arbitrary geometries and density distributions typically involves subdividing the mass body into individual geometric elements (such as rectangular prisms), calculating their gravitational contributions that are then summed up to obtain the gravitational attraction of the whole body. To achieve a more accurate approximation of the true geometric shape and density distribution, this rectangular prism model requires fine dividing, which significantly increases computational load and reduces numerical efficiency. To address this issue, we propose the algorithm for gravimetric forward modeling of arbitrary geometric shapes and density distributions in spectral domain that significantly improves numerical efficiency while preserves computational accuracy. The novelty of our proposed algorithm lies in dividing the masses into multiple layers of equal thickness in the vertical direction, providing constant upper and lower bounds. This allows to extended Parker's formulas and apply the Fast Fourier Transform (FFT) to increase numerical efficiency. The algorithm is tested using synthetic models and then used to compute gravitational effects of topography and sediments using real data from Tibet. Results show high accuracy and numerical efficiency than rectangular prism approach.

SummaryWe propose a frequency-domain finite-element (FDFE) method to model the 2D P-SV waves propagating in porous media. This specific finite element method (FEM) is based on the framework of variational principles, which differ from previously widely used FEMs that rely on the weak formulations of the governing equations. By applying the calculus of variations, we establish the equivalence between solving the stress-strain relations, equations of motion and boundary conditions that govern the propagation of P-SV waves, and determining the extremum or stationarity of a properly defined functional. The structured rectangular element is utilized to partition the entire computational region. We validate the FDFE method by conducting a comparison with an analytically-based method for models of a horizontal planar contact of two poroelastic half-spaces, and a poroelastic half-space with a free surface. The excellent agreements between the analytically-based solutions and the FDFE solutions indicate the effectiveness of the FDFE method in modeling the poroelastic waves. Modeling results manifest that both propagative and diffusive natures of the Biot slow P-wave can be effectively modeled. The proposed FDFE method simulates wave fields in the frequency domain, allowing for easily incorporation of frequency-dependent parameters and enabling parallel computational capabilities at each frequency point (sample). We further employ the developed FDFE method to model two simplified poroelastic reservoirs, one with gas-saturated sandstone and the other with oil-saturated sandstone. The results suggest that changing the fluid phase of the sandstone reservoir from gas to oil can substantially impact the recorded solid and relative fluid-solid displacements. The modeling suggests that the proposed FDFE algorithm is a useful tool for studying the propagation of poroelastic waves.

SummaryThe existence of pores, cracks, and cleavage in rocks results in significant non-linear elastic phenomena. One important non-linear elastic characteristic is the deviation of the stress-strain curve from the linear path predicted by Hooke's law. To provide a more accurate description of the non-linear elastic characteristics of rocks and to characterize the propagation of non-linear elastic waves, we introduce the Preisach-Mayergoyz space model. This model effectively captures the non-linear mesoscopic elasticity of rocks, allowing us to observe the stress-strain and modulus-stress relationships under different stress protocols. Additionally, we analyze the discrete memory characteristics of rocks subjected to cyclic loading. Based on the Preisach-Mayergoyz space model, we develop a new non-linear elastic constitutive relationship in the form of an exponential function. The new constitutive relationship is validated through copropagating acousto-elastic testing, and the experimental result is highly consistent with the data predicted by the theoretical non-linear elastic constitutive relationship. By combining the new non-linear elastic constitutive relationship with the strain-displacement formula and the differential equation of motion, we derive the non-linear elastic wave equation. We numerically solve the non-linear elastic wave equation with the finite difference method and observe two important deformations during the propagation of non-linear elastic waves: amplitude attenuation and dispersion. We also observe wavefront discontinuities and uneven energy distribution in the 2-D wavefield snapshot, which are different from those of linear elastic waves. We qualitatively explain these special manifestations of non-linear elastic wave propagation.

SummaryNumerical simulations of seismic wave propagation in heterogeneous 3D media are central to investigating subsurface structures and understanding earthquake processes, yet are computationally expensive for large problems. This is particularly problematic for full waveform inversion, which typically involves numerous runs of the forward process. In machine learning there has been considerable recent work in the area of operator learning, with a new class of models called neural operators allowing for data-driven solutions to partial differential equations. Recent works in seismology have shown that when neural operators are adequately trained, they can significantly shorten the compute time for wave propagation. However, the memory required for the 3D time domain equations may be prohibitive. In this study, we show that these limitations can be overcome by solving the wave equations in the frequency domain, also known as the Helmholtz equations, since the solutions for a set of frequencies can be determined in parallel. The 3D Helmholtz neural operator is 40 times more memory-efficient than an equivalent time-domain version. We employ a Helmholtz neural operator for 2D and 3D elastic wave modeling, achieving two orders of magnitude acceleration compared to a baseline spectral element method. The neural operator accurately generalizes to variable velocity structures and can be evaluated on denser input meshes than used in the training simulations. We also show that when solving for wavefields strictly on the surface, the accuracy can be significantly improved via a graph neural operator layer. In leveraging automatic differentiation, the proposed method can serve as an alternative to the adjoint-state approach for 3D full-waveform inversion, reducing the computation time by a factor of 350.

SummarySignificant compositional differences may exist in the lithospheric mantle and above the core-mantle boundary (CMB) relative to the ambient mantle. The intrinsic density differences may affect the development of thermal boundary layer (TBL) instabilities associated with lithospheric delamination and formation of thermochemical plumes. In this study, we explored the instability of two-layer thermochemical fluid using two different techniques: marginal stability analysis with a propagator-matrix method and finite element modeling. We investigated both the instabilities in lithospheric mantle (i.e., lithospheric instability) and the mantle above the CMB (i.e., plume-forming instability) using a background temperature Tbg(z) with the TBL. For lithospheric instability, we found that two-layer fluid with free-slip boundary conditions mainly undergoes the same three different convective modes (i.e., two oscillatory convection modes and one layered convection regime) as that with no-slip boundary condition reported in Jaupart et al., (2007). However, with free-slip boundary conditions, the transitions between these convection modes occur at larger values of buoyancy number B. Free-slip boundary conditions lead to smaller critical Rayleigh number Rac, but larger convective wavelength and oscillation frequency ωc, compared with those with no-slip boundary conditions. Our numerical modeling results demonstrate that Rac and ωc predicted from the classical marginal stability analyses using Tbg(z) with TBL temperature may have significant errors when the oscillatory period is comparable with or larger than the timescale of lithospheric thermal diffusion that causes Tbg(z) to vary with time significantly. In this case, using a more gently sloped background temperature profile ignoring the TBL temperature, the stability analysis predicts more accurate stability conditions, thus presenting an effective remedy to the stability analysis. For plume-forming instability, because of the reduced viscosity in the hot and compositionally dense bottom layer, the transition to the layered convection occurs at significantly smaller B values, and in the oscillatory convection regime, Rac is larger but ωc is smaller, compared with those for lithospheric instability. Finally, our study provides a successful benchmark of numerical models of thermochemical convection by comparing Rac and ωc from numerical models with those from the marginal stability analysis.

SummaryIn the Northern Andes, partitioning of oblique subduction of the Nazca plate beneath the South American continent induces a northeastward motion of the North Andean Sliver. The strain resulting from this motion is absorbed by crustal faults, which have produced magnitude 7 + earthquakes historically in the Andean Cordillera of Ecuador and southern Colombia. In order to quantify the strain in that area, we derive a high-resolution surface velocity map using InSAR time-series processing. We analyzed 6 to 8 years of Sentinel-1 data and combined different satellite line-of-sight directions to produce a reliable velocity map in the East direction. We use interpolated GNSS data to express the velocity map with respect to Stable South America and remove the long-wavelength pattern due to the post-seismic deformation following the 2016 Mw 7.8 Pedernales earthquake. The InSAR velocity map finds high E-W shortening strain rates along N-S trending structures within the Western Cordillera and the Interandean valley, with little deformation taking place east of them. This result strengthens the previous proposition of a ∼350 km long Quito-Latacunga tectonic block, forming a restraining bend in the overall right-lateral strike-slip fault system accommodating the northeastward escape motion of the North Andean Sliver. However, the high spatial resolution provided by InSAR indicates that previously proposed boundaries for this block need to be revised. In particular, InSAR results highlight high strain rate (>300 nstrain/yr) along undescribed active structures, south and west of the proposed limits for the Quito-Latacunga block, respectively in Peltetec and Ibarra regions. Interestingly, the two areas with the largest strain rates spatially correlate with the proposed areas of large historical earthquakes. Modeling of the InSAR and GNSS velocities in these areas suggests shallow coupling and high slip rates on structures which, previously, were not identified as active. We also demonstrate a slow-down of the shallow aseismic slip on the Quito fault after the Pedernales earthquake, suggesting that stress changes following large megathrust events might trigger transient slip behaviors on crustal faults. The high-resolution strain map provided by this work provides a new basis for future tectonic models in the Ecuadorian and southern Colombian Andes, and will contribute to the seismic hazard assessment in this highly populated area of the Andes.

SummaryThe grid search method is a common approach to estimate the three spatial coordinates of event hypocenters. However, locating events in large search spaces with small grid spacings is computationally prohibitive. This study accelerates the grid searches over large search spaces using a quadratic interpolation technique. We start with the coarse-grid-estimated location, where we have the minimum value of the difference in the traveltimes between S- and P-waves summed over all receivers. Then, we select the neighbouring grid points and build a 3D quadratic function. The unknown coefficients of the 3D quadratic function are computed by solving a system of linear equations. After that, we interpolate the location by solving partial derivatives of the quadratic function. The quadratic interpolation technique performs well on both synthetic and real microseismic data examples, typically leading to similar event locations as those obtained using 10 times smaller grid spacings in all three directions, at a minor additional computational expense, and without the need to generate traveltimes at new spatial positions.

SummaryThe 1934 Mw 8.2 Bihar-Nepal earthquake was one of the devastating earthquakes, which made seismologists realize the importance of proper seismic hazard analysis and design aspects in India. The event occurred way before proper seismic networks were implemented and hence there are no recorded ground motions available for this event. The present study, thus aims to generate possible ground motions for the 1934 Mw 8.2 Bihar-Nepal event. The complex geographical features, ambiguous source information, and lack of ground motion data make the simulation and validation of ground motions very difficult. In this regard, the broadband (BB) ground motions are simulated and validated for the most recent well-documented Himalayan event, i.e., the 2015 Mw 7.9 Nepal earthquake in order to calibrate the model and simulation methodology. For this purpose, the computational model presented by Sreejaya et al. (2023) is extended up to a region of 1000 km × 670 km (longitude 80-89 °E and latitude 23-30 °N) in the Indo-Gangetic Basin to simulate the low-frequency (LF) ground motions using spectral element method (Komatitsch and Tromp 1999). These deterministically simulated LF ground motions are combined with stochastically simulated high-frequency (HF) ground motions based on an improved seismological model following Otarola and Ruiz (2016). The seismic moment and dimensions of the rupture plane presented by Pettanati et al. (2017) are used to generate ten samples for the finite fault source model having different slip distribution along the rupture plane as a random field (Mai and Beroza 2000; 2002). The BB ground motions (0.01–25 Hz) are obtained by merging LF and HF ground motions in the time domain by matching them at a frequency of ∼0.3 Hz. Such BB results are simulated at a grid of stations and at locations where Modified Mercalli Intensity (MMI) intensity values are available. The estimated MMI values and the observed MMI values are compared to emphasize the efficacy of the model. The maximum PGA estimated from the simulated ground motions in horizontal and vertical directions are observed to be 0.48 g and 0.4 g. Further, 5% damped response spectra and spectral amplification are analyzed concerning the sediment depth of the Indo-Gangetic Basin. The results from the study can serve as inputs for dynamic analysis and the design of earthquake-resistant structures across different locations in the Indo-Gangetic Basin.

SummaryFault geometry is a key factor in controling the mechanics of faulting. However, there is currently limited theoretical knowledge regarding the effect of non-planar fault geometry on earthquake mechanics. Here, we address this gap by introducing an expansion of the relation between fault traction and slip, up to second order, relative to the deviation from a planar fault geometry. This expansion enables the separation of the effects of non-planarities from those of planar faults. This expansion is realised in the boundary integral equation, assuming a small fault slope. It provides an interpretation for the effect of complex fault geometry on fault traction, for any fault geometry and any slip distribution. Hence the results are also independent of the friction that applies on the fault. The findings confirm that fault geometry has a strong influence on in-plane faulting (mode II) by altering the normal traction on the fault and making it more resistant to slipping for any fault geometry. On the contrary, for out-of-plane faulting (mode III), fault geometry has a much smaller influence. Additionally, we analyse some singularities that arise for specific fault geometries often used in earthquake simulations and provide guidelines for their elimination. To conclude this study, we discuss the limits of the infinitesimal strain theory when non-planar faults are considered.

SummaryRocks near a fault plane are commonly damaged by multiple earthquake ruptures, forming damage zones. The damage zone is important structures controlling various properties of a fault, yet its fine scale (tens to hundreds of meters) structure is difficult to resolve with surface seismic observations. We propose to use earthquakes that occur at depth within a fault zone as virtual seismometers (VSs) and use surface observations to extract Green's function (GFs) between VS pairs (VSGFs) . This method resembles that of ambient noise tomography and the retrieved VSGFs are related to the structures between event pairs. In this study, we develop the theory about how to extract VSGFs using surface stations deployed across a fault zone. Firstly, we use a half-space model and Fresnel zone analysis to determine the upper and lower limits of the GF frequency band, which is controlled by the station spacing and aperture of a given seismic array. Then, for VS in a fault zone, we demonstrate that the VSGF can be retrieved by linear seismic arrays deployed across the fault, and that the VSGF is equivalent to waves emitted simultaneously from an array of mirror sources of one event and received by the other. Secondly, the half-space result is directly adopted to determine the corresponding frequency band in the damage zone situation. Thirdly, we analyze different combinations of VS pair geometry and conclude that a relatively larger VS distance (much larger than the damage zone width) is more effective to recover damage zone structures for the available frequency bands. In this situation, VSGFs are trapped waves, that is represented by the interference of mirror sources. In such a case, the trapped waves are equivalent to surface waves, which have dispersion features to extract damage zone structures. Finally, we adopt the VSGF method to the Ridgecrest earthquake aftershock monitoring array and use a profile of aftershocks to extract 6 pairs of VSGFs. The spatial variation of VSGFs may reflect the depth-dependent variation of damaged zone. Our analysis shows a promising direction to use VSGFs to extract spatial variations of fault damaged zones.

SummaryJoint inversion, such as the combination of receiver function and surface wave dispersion, can significantly improve subsurface imaging by exploiting their complementary sensitivities. Bayesian methods have been demonstrated to be effective in this field. However, there are practical challenges associated with this approach. Notably, most Bayesian methods, such as the Markov Chain Monte Carlo (MCMC) method, are computationally intensive. Additionally, accurately determining the data noise across different data sets to ensure effective inversion is often a complex task. This study explores the unscented Kalman inversion (UKI) as a potential alternative. Through a data-driven approach to adjust estimated noise levels, we can achieve a balance between actual noise and the weights assigned to different data sets, enhancing the effectiveness of the inversion process. Synthetic tests of joint inversion of receiver function and surface wave dispersions indicate that the UKI can provide robust solutions across a range of data noise levels. Furthermore, we apply the UKI to real data from seismic arrays in Pamir and evaluate the accuracy of the joint inversion through posterior Gaussian distribution. Our results demonstrate that the UKI presents a promising supplement to conventional Bayesian methods in the joint inversion of geophysical data sets with superior computational efficiency.

SummaryThe left-lateral Xianshuihe fault is a seismically active fault system located at the eastern boundary of the Tibetan Plateau. We analyzed the Sentinel InSAR data from 2014 to 2021 to study the temporal and spatial variation of fault creep along the Xianshuihe fault. We applied the InSAR stacking method and the coherence-based SBAS method to derive the Line-Of-Sight (LOS) velocity map and time-series from both the ascending and descending orbits. We studied both the secular component and the time-dependent component of surface deformation from InSAR. We compare the InSAR-derived velocity maps with the GPS-derived velocity field and found that these two independent measurements are consistent. A 200 km long creeping section is identified along the central segment of the Xianshuihe fault. The surface creep rate is measured to be ranging from 0 to 6 mm yr−1. We combined the elastic dislocation model and the InSAR velocity maps to invert for the geodetic fault slip rate and the aseismic slip distribution in the upper crust. The secular fault creep model shows that most of the Xianshuihe fault is creeping at depth. The time-dependent fault creep model indicates that the maximum aseismic slip rate from Bamei to Kangding accelerated from 30 mm yr−1 to 40 mm yr−1 and then decayed to 5 mm yr−1 from 2014 to 2021. The fully creeping segment of the Xianshuihe fault seems to become a partially locked segment in a short time period (a couple of years). We suspect that the acceleration of fault creep from 2017 to 2019 is linked to dynamic triggering by passing seismic waves or fluid migration. Finally, we compare the temporal variation of fault creep with previous studies and discuss the earthquake hazard implications.

SummaryAntarctica has been proposed as a significant source of the meltwater that entered the oceans during Meltwater Pulse 1B (MWP1B) approximately 11,500 years ago. Support for this scenario has been provided by evidence that the deep fjords of coastal Antarctica, which were heavily glaciated at the maximum of glaciation, were deglaciated at this time. Further support for this scenario was provided by the observation that the inter-hemispheric sea level teleconnection associated with significant southern hemisphere deglaciation at this time provided an explanation of the highly non-monotonic relative sea level histories recorded at sites on the coast of Scotland, a region which had also been heavily glaciated at the last glacial maximum. Furthermore, it has been argued that a significant contribution to MWP1B must have also been delivered to the oceans by the abrupt northern hemisphere warming that occurred at the end of the Younger Dryas (YD) cold reversal, which also occurred approximately 11,500 years ago. Our focus in the present paper is to distinguish between these two possible primary sources of MWP1B. The investigation of how local alterations to ice thicknesses are able to explain evidence which has previously been used to argue for an Antarctic dominant MWP1B will lead us to the conclusion that the Laurentide may be primary source of MWP1B.

SummaryThe complete catalog of moment tensor (MT) solutions is essential for a wide range of research in solid earth science. However, the number of reliable MT solutions for small to moderate earthquakes (3.0 ≤ M ≤ 5.5) is limited due to uncertainties arising from data and theoretical errors. In this study, we develop a new procedure to enhance the resolvability of MT solutions and provide more reliable uncertainty estimates for these smaller to moderate earthquakes. This procedure is fully automatic and efficiently accounts for both data and theoretical errors through two sets of hybrid linear-nonlinear Bayesian inversions. In the inversion process, the covariance matrix is estimated using an empirical approach: the data covariance matrix is derived from the pre-event noise and the theoretical covariance matrix is derived from the residuals of the initial solution. We conducted tests using synthetic data generated from the 3D velocity model and interference from background seismic noise. The tests found that using a combination of the non-Toeplitz data covariance matrix and the Toeplitz theoretical covariance matrix improves the solution and its uncertainties. Test results also suggest that including a theoretical covariance matrix when analyzing MT in complex tectonic regions is essential, even if we have the best 1D velocity model. The application to earthquakes in the northern region of the Banda Arc resulted in the first published Regional Moment Tensor (RMT) catalog, containing more than three times the number of trusted solutions compared to the Global Centroid Moment Tensor (GCMT) and the Indonesian Agency for Meteorology Climatology and Geophysics Moment Tensor (BMKG-MT) catalog. The comparison shows that the trusted solutions align well with the focal mechanism of the GCMT and BMKG-MT, as well as with the maximum horizontal stress of the World Stress Map, and tectonic conditions in the study area. The newly obtained focal mechanisms provide several key findings: (i) They confirm that the deformation in the northern and eastern parts of Seram Island is influenced by oblique intraplate convergence rather than by the subduction process; (ii) They validate the newly identified Amahai Fault with a greater number of focal mechanisms; (iii) They reveal an earthquake Mw 4.7 with the same location and source mechanism six years before the 2019 Ambon-Kairatu earthquake (Mw 6.5) which occurred on a previously unidentified fault.

SummaryWithin the last decade, thin ultra-low velocity zone (ULVZ) layering, sitting directly on top of the core-mantle boundary (CMB), has begun to be investigated using the flip-reverse-stack (FRS) method. In this method, pre- and post-cursor arrivals that are symmetrical in time about the ScS arrival, but with opposite polarities, are stacked. This same methodology has also been applied to high velocity layering, with indications that ultra-high velocity zones (UHVZs) may also exist. Thus far, studies using the FRS technique have relied on 1-D synthetic predictions to infer material properties of ULVZs. 1-D ULVZ models predominantly show a SdS precursor that reflects off the top of the ULVZ and an ScscS reverberation within the ULVZ that arrives as a postcursor. 1-D UHVZ models are more complex and have a different number of arrivals depending on a variety of factors including UHVZ thickness, velocity contrast, and lateral extent. 1-D modeling approaches assume that lower mantle heterogeneity is constant and continuous everywhere across the lower mantle. However, lower mantle features display lateral heterogeneity and are either finite in extent or display local thickness variations. We examine the interaction of the ScS wavefield with ULVZs and UHVZs in 2.5-D geometries of finite extent. We show that multiple additional arrivals exist that are not present in 1-D predictions. In particular, multipath ScS arrivals as well as additional postcursor arrivals are generated. Subsequent processing by the FRS method generates complicated FRS traces with multiple peaks. Furthermore, post-cursor arrivals can be generated even when the ScS ray path does not directly strike the heterogeneity from above. Analyzing these predictions for 2.5-D models using 1-D modeling techniques demonstrates that a cautious approach must be adopted in utilization and interpretion of FRS traces to determine if the ScS wavefield is interacting with a ULVZ or UHVZ through a direct strike on the top of the feature. In particular, travel-time delays or advances of the ScS arrival should be documented and symmetrical opposite polarity arrivals should be demonstrated to exist around ScS. The latter can be quantified by calculation of a time domain multiplication trace. Because multiple postcursor arrivals are generated by finite length heterogeneities, interpretation should be confined to single layer models rather than to interpret the additional peaks as internal layering. Furthermore, strong tradeoffs exist between S-wave velocity perturbation and thickness making estimations of ULVZ or UHVZ elastic parameters highly uncertain. We test our analysis methods using data from an event occurring in the Fiji-Tonga region recorded in North America. The ScS bounce points for this event sample the CMB region to the southeast of Hawaii, in a region where ULVZs have been identified in several recent studies. We see additional evidence for a ULVZ in this region centered at 14° N and 153° W with a lateral scale of at least 250 × 360 km. Assuming a constant S-wave velocity decrease of -10 or -20% with respect to the PREM model implies a ULVZ thickness of up to 16 or 9 km respectively.

SummaryUnderstanding the crustal seismic characteristics of tectonically active regions is crucial for seismic hazard assessment. The study conducted in NW Iran utilized surface wave tomography, radial anisotropy, and density information to analyze the complex crustal structure of the region, which is outstanding because of diverse tectonic features, sedimentary basins, and volcanic formations. By selecting a dataset of 1243 events out of over 3,500 earthquakes with M>4, and employing strict data selection criteria (such as SNR, M, Δ), the researchers calculated Rayleigh and Love wave group velocity dispersion curves using Gaussian multiple filters and phase-matched filtering. The tomographic procedure was initiated by excluding data with residuals > 2σ for enhanced stability. Individual inversions were then carried out for local Rayleigh and Love wave dispersion measurements to obtain 1D VSV and VSH models. Radial anisotropy and VS iso were determined through a discrepancy and averaging of the obtained VSH and VSV, respectively. Gravity modeling was also employed alongside surface wave analysis to understand the region's complex geology, revealing insights into upper-middle-lower crust boundaries, subsurface structures, and Moho depths. The study's velocity maps reveal significant findings related to geological units and tectonic features in various regions based on the provided results. Low velocities in the South Caspian Basin (SCB) and Kura Depression (KD) regions are attributed to substantial sedimentary layers, while low velocities, and depth of VS in NW Iran and Eastern Anatolian Accretionary Complex (EAAC) regions suggest the presence of partially molten materials in the upper and middle crust. The Sanandaj-Sirjan Zone (SSZ) region shows a low velocity anomaly in longer periods and greater depths of VS, surrounded by normal to high velocities, indicating a thick middle crust. Analyzing radial anisotropy and VS iso profiles offers insights into upper-middle-lower crust boundaries, subsurface structures, and Moho depths, highlighting middle crust thickening and lower crust thinning beneath the SSZ. The study confirms the gentle subduction of the SCB oceanic-like lower crust beneath NW Iran in the Talesh (TAL) region, with a rigid middle crust. Additionally, cross-sections reveal igneous laccoliths underplate with a VS iso of 3.7 km/s in the volcanic region. The difference observed by subtracting the velocity models at two adjacent depths, combined with parametric test results, indicates that the Sahand volcanic system is clearly identifiable, while the influence of subtle subduction on the Sabalan volcano at depths up to 30 km remains less distinct. The magma chamber beneath Sahand is situated at depths ranging from 18 to 25 km.

SummaryA new automated algorithm for picking the arrival times of the global P-, SH- and SV-wave phases from multi-component seismic waveform data is presented. This picker is based on a sequential approach using autoregressive prediction of the filtered waveform in a sliding time window, the Akaike-Information-Criterion and the Hilbert transform of the original waveform. The quality of the individual picks is computed by combining signal-to-noise ratios and higher order statistics into a single measure. Synthetic tests are used to find values for high and low quality thresholds. The algorithm is applied to a global data set of waveforms from teleseismic events with magnitude 6 or higher that occurred between 1990 and 2019. This resulted in approximately 4 million P-phase arrival times as well as approximately 3 million SH- and SV-phase arrival times each. These automatic picks are compared to approximately 830 000 manual P-picks as well as approximately 70 000 manual S-picks from the ISC-EHB catalogue. An upper bound for the picking errors of the automatic picks is estimated by using high quality picks of neighbouring stations. This upper bound is found to be 0.55s for the P-picks and 4.3s for the S-picks. If only high quality picks are considered, this represents 50 per cent of the P-picks and 25 per cent of the S-picks, then these errors decrease to 0.35s for the P-picks, and 1.5s for the S-picks, respectively. As a by-product of the picking, the dominant periods of the arriving signals are determined as well.

SummaryThe Baihetan Reservoir, the second largest in the world, is located at the intersection of multiple large active fault zones on the eastern boundary of the Sichuan-Yunnan rhombic block. After impoundment on April 6, 2021, many earthquakes occurred around the reservoir area submerged by water. The largest ML 4.7 earthquake in the reservoir area occurred after the water level reached its highest point. But the seismogenic structures and mechanisms of earthquakes in the reservoir area are still unclear. Based on dense array data from the reservoir area, this paper uses the experimental site sub-model (CSES) of USTC-Pickers, transfer learned with “DiTing” dataset of China to obtain a high-precision earthquake catalog that is twice as large as that the manual catalog. This study show that earthquakes in the reservoir region primarily occur on secondary faults of pre-existing ones, characterized by a prominent feature of high dip angles trending northwest to southeast. Combined with the spatiotemperal migration characteristics of earthquakes and the relationship between earthquakes and water levels, we infer that most earthquakes are rapid response type and may be induced by rapid increase in elastic stress. Only the spatiotemporal distribution image of the ML 3.2 earthquakes sequence in the dam site-Toudaogou section conforms to the law of pore pressure diffusion, and belongs to the fast response type, which may be induced by the poroelasiticity coupling mechanism. The ML 3.0 earthquake swarm with deep depths in the Heishui River section belongs to the delayed response type and may be induced by the poroelasiticity coupling mechanism.

SummaryThe recent developments in array-based surface-wave tomography have made it possible to directly measure apparent phase velocities through wavefront tracking. While directionally dependent measurements have been used to infer intrinsic $2\psi $ azimuthal anisotropy (with a 180° periodicity), a few studies have also demonstrated strong but spurious $1\psi $ azimuthal anisotropy (360° periodicity) near major structure boundaries particularly for long period surface waves. In such observations, Rayleigh waves propagating in the direction perpendicular to the boundary from the slow to the fast side persistently show a higher apparent velocity compared to waves propagating in the opposite direction. In this study, we conduct numerical and theoretical investigations to explore the effect of scattering on the apparent Rayleigh-wave phase velocity measurement. Using two-dimensional spectral-element numerical wavefield simulations, we first reproduce the observation that waves propagating in opposite directions show different apparent phase velocities when passing through a major velocity contrast. Based on mode coupling theory and the locked mode approximation, we then investigate the effect of the scattered fundamental-mode Rayleigh wave and body waves interfering with the incident Rayleigh wave separately. We show that scattered fundamental-mode Rayleigh waves, while dominating the scattered wavefield, mostly cause short wavelength apparent phase velocity variations that could only be studied if the station spacing is less than about one tenth of the surface wave wavelength. Scattered body waves, on the other hand, cause longer wavelength velocity variations that correspond to the existing real data observations. Because of the sensitivity of the $1\psi $ apparent anisotropy to velocity contrasts, incorporating such measurements in surface wave tomography could improve the resolution and sharpen the structural boundaries of the inverted model.

SummaryA high-precision and high-resolution vertical velocity for the Chinese mainland is obtained by integrating precise leveling and GNSS data, using a Helmert joint adjustment method. The results show that the surface vertical rates range between -3.0 and 3.9 mm/yr with continuous deformation in most areas, except the obvious subsidence at the rates of -15.0 to -94.2 mm/yr induced by groundwater exploitation in the North China Plain. Particularly, the central and southern Tibet, Tien Shan, Alashan, Ordos, eastern Cathaysia, and Northeast China uplift at the rates of 0.5 – 3.9 mm/yr; the southeastern Tibetan Plateau, Sichuan basin, and Yangtze block are dominated by surface subsidence at the rates of -3.3 to -0.5 mm/yr. Furthermore, the vertical rates vary little between the eastern and western regions of the Chinese mainland despite their pronounced differences in horizontal deformations. The effects of gravity isostasy and non-tectonic factors, including the environmental mass loads, Glacier Isostatic Adjustment (GIA), poroelastic expansion/compression, and mining operations have partially contributed to the vertical deformation of the Chinese mainland. Overall, this velocity reflects the complicated deformation features induced by the multiple geodynamic processes of the Chinese mainland. These geodynamic processes include isostasy, orogenic processes, and geothermal anomalies associated with slab subduction/plate collision.