Geophysical Journal International

SummarySeismic scattering waves in random media are usually regarded as noise in conventional seismic imaging, inversion and interpretation. However, the spatial and temporal variation of the scattering energy depends on the stochastic properties of the random media. The extraction of heterogeneity information such as the correlation scale and fluctuation strength from seismic scattering waves remains a challenge. These parameters are inverted from real scattering data by fitting the synthetic envelopes to the observed seismic envelopes. The synthetic envelopes are usually computed using the Monte-Carlo radiative transfer (MCRT) method. However, physical verification of the stochastic parameter inversion based on MCRT theory has not been realized although it is believed to be correct. To this end, we conducted a physical modelling experiment using an ultrasonic acquisition system and recorded the transmitted wavefields through an artificial heterogeneous medium. In this paper, the elastic MCRT method was used to simulate the energy transport, and the correlation length and fluctuation strength of the artificial heterogeneous medium were inverted with a revised objective function, which can better balance the energy level of direct waves and scattering waves in the inversion process. The inversion results of the correlation scale and fluctuation strength match well with true values, suggesting that this method is accurate and reliable. A combination of our physical experiments and the MCRT theory gives strong proof that this inversion method is correct. Therefore, it can be used with confidence to estimate the properties of the heterogeneities from the ‘undesired’ scattering waves, both in the oil/gas exploration and earth structure investigation.

SummaryThe attenuation operator t* represents the total path attenuation and characterizes the amplitude decay of a propagating seismic wave. Calculating t* is typically required in seismic attenuation tomography. Traditional methods for calculating t* require determining the ray path explicitly. However, ray tracing can be computationally intensive when processing large datasets, and conventional ray tracing techniques may fail even in mildly heterogeneous media. In this study, we propose a modified fast sweeping method (MFSM) to solve the governing equation for t* without explicitly calculating the ray path. The approach consists of two main steps. First, the traveltime field is calculated by numerically solving the eikonal equation using the fast sweeping method. Second, t* is computed by solving its governing equation with the MFSM, based on the discretization of the gradient of t* using an upwinding scheme derived from the traveltime gradient. The MFSM is rigorously validated through comparisons with analytical solutions and by examining t* errors under grid refinement in both simple and complex models. Key performance metrics, including convergence, number of iterations, and computation time, are evaluated. Two versions of the MFSM are developed for both Cartesian and spherical coordinate systems. We demonstrate the practical applicability of the developed MFSM in calculating t* in North Island, and discuss the method’s efficiency in estimating earthquake response spectra.

SummaryVarious methods for determining the magnetic field at the core-mantle boundary (CMB) from the observed geomagnetic core field have been explored over recent decades. These include the harmonic downward continuation of surface data and the stabilised iterative upward continuation. The instability of the inverted poloidal magnetic field at the CMB for a radial conductivity structure is complemented by the non-uniqueness of determining the toroidal magnetic field at the CMB for a laterally inhomogeneous conductivity model. We reformulate this unstable and non-unique inverse problem as an iterative upward continuation approach, in which the magnetic field at the CMB is successively updated. The uniqueness of the inverse solution is ensured by the initial choice of the toroidal magnetic field at the CMB, while the stability is achieved by stopping the iterations once the desired tolerance is reached between the spectral index of the updated solution and that obtained from numerical geodynamo simulations. We consider two significantly different radial electrical conductivity models of the lower mantle, each with conductance near 108 S: conductivity model A, based on external electromagnetic sounding, which includes a significant conductivity increase in a 10 km thick layer above the CMB, and conductivity model B, characterized by a gradual conductivity increase determined from the Voigt-Reuss-Hill average of the bridgmanite-ferropericlase aggregate, with an additional conductivity increase in the 300 km thick D” layer associated with post-perovskite. Models A and B bracket the lower and upper bounds of conductivity structures derived from thermal and compositional constraints below 1600 km depth. We find that the differences between the magnetic field components at the CMB inverted for models A and B are approximately 1-2 per cent of the total field. To explore lateral variations, we construct a synthetic model of the Pacific and African superplumes by simplifying their geometric shapes, estimating the temperature increase within the plumes and allowing mantle mineral activation energies to vary only with temperature. Our results show that, in the regions of the superplumes, the poloidal and toroidal magnetic fields at the CMB change by approximately 12,000 nT and 2,500 nT, respectively. The changes in the horizontal poloidal field at the CMB are comparable in magnitude to those resulting from substituting model A with model B. However, the changes in the radial field inverted for the three-dimensional plume conductivity model are significantly larger than those arising from replacing model A with model B.

SummarySeismic surface wave tomography, particularly when leveraging dense array data, has become a widely used method for investigating shallow subsurface velocity structures. The shallow structures are usually characterized by rapid seismic velocity changes (i.e. seismic interfaces) due to variations in rock properties, sedimentary environments, or tectonic features. However, the commonly used grid-based parameterization of the velocity field in surface wave tomography often struggles to accurately constrain such interface geometries. In addition, traditional surface wave inversion methods typically rely on 1D inversion at individual stations using dispersion curves, followed by interpolation to construct 2D or 3D models. This approach can sometimes introduce spurious features and reduce model reliability. To address these limitations, we propose a geological and level-set parameterization approach for surface wave tomography, allowing for the explicit consideration of interface structures in inversion. This method is then combined with the Ensemble Kalman Inversion to optimize subsurface structures. Synthetic tests demonstrate that integrating 3D interface parameterization in tomography significantly enhances the reliability of the velocity model and the recovery of interface geometries. Applying this approach to the Woxi gold mine region in China yielded inversion results that closely align with existing borehole data. This study highlights the advantages of level-set parameterization for 3D interface imaging in seismic tomography, underscoring its potential in subsurface mineral exploration.

SummaryIn order to better understand the regional tectonics of western part of Africa (WA) and adjacent islands, joint inversion (Jinv) of body wave and surface wave measurements is conducted to construct new sets of crustal models. Teleseismic P-wave receiver function, receiver function horizontal-to-vertical ratio and Rayleigh wave ellipticity are jointly inverted based on a fast simulated-annealing scheme. All three types of observables are derived from single-station recordings and are primarily sensitive to structures beneath the station. The integration of these datasets through Jinv allows for complementary constraints, thereby improving the resolution of crustal velocity structures and the characterization of velocity variations with depth. We present improved and some new crustal structure parameters including bulk crustal VP/VS ratio, crustal thickness (H) estimates, and shear-wave velocity (VS) models beneath 25 broad-band seismic stations across inland, coastal, and island terrains. Using an improved approach involving the correction of misorientation error effect from seismic waveform data, the data quality is well-enhanced leading to improved resolutions of structures across the different terrains. Results from H-k and crustal models showed a general northward thinning from Congo Craton (> ∼48 km) towards the Lower Benue Trough (∼15 km), and from coastal terrain along Gulf of Guinea (< ∼44 km) towards Mauritanian Belt (> ∼16 km). Compared to other terrains, the islands show very thin depth to the Moho, but higher than the global estimates. In the Mauritanian-Senegal Basin, sharp differential in crustal thickness and Jinv results at neigbouring G.SOK and G.MBO are observed, where slower Vs revealed a LVZ anomaly at G.SOK in contrast with faster Vs at G.MBO—which could be due to local subsidence from sediment loading, or uplift from tectonic activities. In the upper-middle crust, the Jinv imaged structures with faster VS characteristic of felsic to intermediate bulk crustal composition beneath inland terrain (West Africa Craton, Congo Craton, Hoggar), attributed to highly depleted and stable nature of the cratonic lithosphere, contributing to faster VS compared to other terrains. Low velocity structures underlying the island stations are attributed to partial melts and high temperature materials, indicative of volcanic and Basaltic composition. Similarly, the low velocity structures deciphered beneath coastal stations G.SOK and AF.EDA could be related to the structures in their adjacent areas of Tenerife and the Cameroon Volcanic Line, respectively. The nbroad range of VP/VS (∼1.58–1.85) ratio along the coastal terrains demonstrates its complexity; from the low VP/VS which may be attributed to indurated or low porosity sedimentary materials, and high VP/VS —typical of cracks, fluids inundated sedimentary or volcanic materials. Island terrain are associated with higher bulk VP/VS indicative of volcanics and Mafic-Basaltic materials, with the low velocity zones (LVZs) suggestive of the presence of magmatic materials. These broad crustal configuration highlights the complexity and provides new insight for developing more accurate regional model for western Africa and its adjacent islands, and global reference models in future studies.

SummaryThe Xiangshan volcanic basin locates in southeast China hosts the world’s third-largest volcanogenic uranium deposit. However, the structure of the volcanic system remains poorly resolved, limiting insights into the uranium mineralization. To address this, we conducted a joint inversion of gravity and magnetic data collected in the basin. Our inversion results reveal a southeast-dipping porphyroclastic lava conduit beneath the peak of Mount Xiangshan, characterized by low density and high magnetic susceptibility. A southwest-dipping volcanic conduit has also been identified beneath the rhyodacite crater in the Shutang area of the western basin. It connects to the porphyroclastic lava conduit in the deep. Both of these volcanic conduits are controlled by an EW-trending, low-density basement fault zone. This spatial relationship indicates that the volcanic eruptions in the western basin share a common subvolcanic plumbing system, which collectively acted as principal pathways for ore‑forming hydrothermal fluids and uranium enrichment. These results underscore the role of volcanic-intrusive architecture in controlling the mineralization processes in the Xiangshan volcanic basin.

SummaryThe South American continent (SAC), a region of pronounced geodynamic and hydrological activity, exhibits crustal deformation and gravity field anomalies driven by the interplay of tectonic forces and surface/subsurface mass redistribution. While previous studies have mainly focused on gravity changes driven by terrestrial water storage (TWS), mass variations of the solid Earth remain inadequately addressed. In this study, we resolve deep-seated mass transport Gravity Recovery and Climate Experiment (GRACE) satellite gravimetry, hydrological model outputs, GPS-derived vertical crustal motions, and glacial isostatic adjustment (GIA) correction. Our results reveal an internal mass variation of 0.21 ± 0.45 cm yr -1 in equivalent water height (EWH), independent of surface hydrological contributions. Interpreting this signal as predominantly driven by crust-mantle boundary (Moho) displacements, we estimate an average Moho depth uplift rate of 0.37 ± 0.80 cm yr -1 across SAC, based on the crust–mantle density contrast. The Moho interface depth variations exhibit significant spatial heterogeneity. Through uncertainty analysis, four distinct regions (A, B, C, and D) are identified: Region A exhibits Moho uplift and Region B exhibits subsidence, with part contributions from the isostatic adjustment. Key uncertainties in these estimates stem from sedimentation effects and the accuracy of current observations or models. Subsidence in Region C and uplift in Region D are related to the co-seismic and post-seismic effects of the 2010 Chile earthquake. These findings underscore the significance of solid Earth mass flux in active continental regions and unravel the mechanisms governing crust-Moho mass redistribution.

SummaryRecent advancements in high-resolution Digital Elevation Models (DEMs) derived from LIDAR and satellite radar technologies have added a new dimension to the determination of height systems and geoid models. However, their benefits are limited by simplified assumptions inherited from past practices. In mountainous areas, taking into consideration of topography as the Bouguer plate or employing inaccurate terrain corrections can constitute to a problematic approach. Even though the gravity reduction procedures mentioned above have been enhanced in geoid determination studies, the Helmert orthometric heights based on them are still used in some countries such as Türkiye and Taiwan. It is inevitable that this contradiction will negatively affect geoid modeling studies that are intended to be verified or combined with GNSS/leveling data. Another issue arises by ignoring density variations of topographic masses. Through a comparative analysis, this study reviews combined and individual impacts of terrain roughness and density variations on geoid models in the Konya Closed Basin (KCB) and the Auvergne regions, with a focus on their distinctive topographical characteristics. Using 1″ DEMs of the SRTM mission and 30″ UNB_TopoDensT lateral density models, we reveal that terrain corrections in gravity reductions significantly affect geoid heights, with deviations of up to 11.9 cm in KCB and 4.2 cm in Auvergne. Incorporating lateral density models has resulted in geoid height discrepancies of up to 26.8 cm in KCB and 6.7 cm in Auvergne. A validation strategy implemented through GNSS/leveling paths showed that terrain corrections markedly improved geoid model accuracy, particularly in relation to elevation. However, the contribution of the UNB_TopoDensT model to geoid accuracy is questionable in terms of accuracy. Notably, applying density values below 2.4 g·cm⁻³ in high-altitude regions can lead to disruptive effects on geoid determination. This result is underscoring of the need on a realistic modeling of topographical densities in high elevated and rugged terrains. A further conclusion that emerged from these analyses is that gravimetric geoid models should be verified by rigorous orthometric heights, which are observed to fit them better at the 1-2 cm level, instead of the Helmert orthometric heights.

AbstractThis study introduces a novel method for performing 3D inversion of magnetotelluric (MT) data. The proposed method, referred to as the radiation boundary scheme, employs a two-step simulation strategy for the computation of both forward and adjoint responses. One of the key advantages of the scheme is its ability to handle arbitrarily shaped inversion domains, thereby optimizing the number of unknown model parameters by discarding model parameters that are not constrained by the data. Consequently, it significantly improves accuracy and computational speed as compared to traditional inversion algorithms. The effectiveness of the developed algorithm is demonstrated through a comprehensive analysis of 3D inversion using synthetic and continental-scale (SAMTEX) MT data. The method’s efficiency facilitates a detailed analysis of large-scale MT data inversion. Through numerical experiments, it is observed that using a coarse mesh for inversion, the resolution is compromised in the shallower part, resulting in inferior imaging and, consequently, affecting the estimation of resistivity value in the deeper subsurface. The detailed numerical experiments indicate that performing a fine-scale inversion on a small portion of the survey data utilizing a coarsely inverted model may run into a local minimum. Hence, caution should be exercised in employing such an approach. Instead, the investigations suggest simultaneously executing a fine-scale inversion on the entire data set. The forward/adjoint problem can be solved with a two-order higher tolerance as compared to the conventional finite-difference-based inversion algorithm. Therefore, the proposed algorithm holds significant value for the MT inversion of large data sets.

SummaryMineral exploration is frequently centred around delineating discrete geological units with, typically, sharp boundaries that could represent economic targets. In the case of complex resistivity (CR) inversions, the choice of regularisation and model parameterisation significantly impacts the inversion’s ability to delineate targets. Initially, however, a prudent researcher may not wish to bias their inversion towards sharp distinct units without prior justification. Here we explore how a suite of regularisation approaches to the CR inverse problem allows to encompass different classes of prior beliefs. We present these as a progression as more information becomes available regarding the likelihood of distinct geological units. The most weakly informed approach with respect to the delineation of geological units we consider is the classic ℓ2-type regularisation, tend to produce smeared-out fuzzy images. However this is typically not what is expected for distinct geological units, and we compare this with schemes that increasingly resolve sharp boundaries. We test a range of ℓ1-type regularisations, which have been frequently touted in the geophysics and optimisation literature as being well-suited for such tasks. We experiment with using a so-called overcomplete parameterisation of the CR field, which aims to separate smooth background and sharp foreground features. These ℓ1 schemes are shown to produce generally sharper images than ℓ2. In the most informed case, where strong assumptions can be made about the local geology, we represent the CR field as a foreground ellipse in a homogenous background. This approach significantly reduces the size of the parameter space, and tends to have a simple geometric interpretation. While the anomaly parameterisation has some unique challenges, we show it clearly resolves distinct units compared to both the ℓ2 and ℓ1 regularisations. Applications first to synthetic data and then to field data from Century Zinc Deposit in northern Australia, demonstrate the progression from weakly informed to strongly informed regularisation and parameterisation and the sharpness of the recovered geological units.

SummaryFollowing reanalysis of data from 8 seismic networks that operated in the region surrounding the Pampean flat slab during the past several decades, we generated 3D images of Vp, Vs, and Vp/Vs from a combination of arrival times of P and S waves from local earthquakes, and Rayleigh wave dispersion curves from both ambient noise and existing shear wave models. Among the robust features in these images is a low velocity, root-like structure that extends beneath the high Andes to a deflection in the flat slab, which suggests the presence of an overthickened Andean crust rather than a hypothesized continental lithospheric root. Most of the larger scale features observed in both the subducted Nazca plate and the overriding continental lithosphere are related to the intense seismic activity in and around the Juan Fernandez Ridge Seismic Zone (JFRSZ). Vp/Vs ratios beneath, within, and above the JFRSZ are generally lower (∼1.65–1.68) than those in the surrounding Nazca and continental lithosphere (∼1.74–1.80). While the higher continental lithosphere ratios are due to reduced Vs and likely a result of hydration, the lower JFRSZ related ratios are due to reduced Vp and can be explained by increased silica and CO2 originating from beneath the slab, perhaps in concert with supercritical fluid located within the fracture and fault networks associated with the JFR. These and related features such as a region of high Vp and Vs observed at the leading edge of the JFRSZ are consistent with a basal displacement model previously proposed for the Laramide flat-slab event, in which the eroded base of the continental lithosphere accumulates as a keel at the front end of the flat slab while compressional horizontal stresses cause it to buckle. An initial concave up bend in the slab facilitates the infiltration of silica and CO2-rich melts from beneath the slab in a manner analogous to petit spot volcanism, while a second, concave down bend, releases CO2 and supercritical fluid into the overlying continental lithosphere.

SummaryTomographic inversion of traveltime picks from both P-wave and S-wave wide-angle seismic data acquired along and across the Louisville Ridge Seamount Chain (LRSC) provides key insights into its magmatic construction and subsequent subduction-related deformation. Our P-wave velocity-depth models reveal that each seamount along the LRSC comprises an intrusive mafic-ultramafic core that rises within the crust to within 1–2 km of the seabed summit (P-wave velocity, Vp = 5.5–6.5 km s−1; S-wave velocity, Vs < 3.6 km s−1), with each underlain by a crustal root ∼4–5 km thick. Notably, Canopus seamount comprises two adjacent eruptive centres, and our modelling shows that the more northern is currently being internally deformed as it rides up (ascends) the Tonga-Kermadec Trench (TKT)-related plate bending outer rise. Lateral variation in Vs within models along and across the LRSC also primarily reflects subduction-related deformation, with low-velocity regions corresponding to large-scale faulting constrained within the crust. Comparison of pre- and post-LRSC-TKT collision forearc crustal structure indicates that bulk Vp properties recover within ∼50 kyr, whereas Vs structure retains it fault-related fabric for at least ∼740 kyr. Vp/Vs ratios (1.75–1.85) confirm a magmatic origin for all LRSC seamounts, with evidence of localized water-filled cracks due to seawater infiltration along faults, particularly beneath the TKT-ward side of the Osbourn seamount. Estimated water content within the upper crust ranges from 12–15 per cent by weight, decreasing to < 10 per cent in the mid-lower crust, with no evidence of > 12 per cent water content within the Pacific crust being subducted. In comparison with post-collision subduction further north, where the observed upper mantle velocity suggests up to 30 per cent water content, our models suggest that, although deformed and faulted as part of subduction, the LRSC appears more resistant to this deformation than the background Pacific crust adjacent. Our findings provide new constraints on the mechanical and compositional evolution of the LRSC, both prior to and during its collision with the overriding Indo-Australian plate.

SummaryThis research had an initial goal to quantitatively fit and then separate an induced polarization (IP) contribution to extensive ground electromagnetic (EM) data from the Girrilambone area, NSW. A secondary goal identified during the study was to explain why inversion of data from two different EM systems covering the same area each consistently predicted different IP time-constants and chargeabilities. The mineral exploration area was originally surveyed by a 6.25 Hz central loop SIROTEM survey measuring dB/dt. The area was later resurveyed with 1 Hz base-frequency Slingram survey using a Landtem B field sensor. The targets were economic sulphides at depth, with expected signatures being slowly decaying EM responses of small amplitude. Most of the data was affected by inductive IP effects of negative sign, with potential late-delay time EM responses of positive sign obscured. The Girrilambone area surveyed includes the Tritton Mine, discovered in 1995 as a result of the 6.25 Hz SIROTEM survey. To enable the subtraction of IP effects from the EM data, our primary goal, we used the EM data to predict Cole-Cole IP parameters that are consistent with documented values associated with extensive in-situ regolith clay resulting from weathering. The data sets were inverted using a polarisable thin-sheet model that estimated regolith conductivity-thickness or conductance S, chargeability m, IP frequency dependence c and conductivity IP time constant τσ. The thin sheet model was generally able to fit the observed responses, with the fitted IP contribution subtracted from the observed data to produce an ‘IP corrected’ data set of EM data more suitable for the detection of slow decays indicative of sulphide targets. The 6.25 Hz dB/dt data was however modelled with quite different parameters to the1 Hz B field data. The 6.25 Hz IP conductivity time constant was smaller by a factor of 10 while the chargeability was smaller by a factor of more than 2. This initial goal of the research was achieved in that subtraction of the fitted IP contributions in either case improved the capability to identify deeper conductive targets. We are confident that the systematic differences in fitted IP conductivity time constant and chargeability are not due to data or system description error, or to inversion constraints. We conclude that TEM systems will not accurately estimate intrinsic IP conductivity time-constants as rigorously defined from wideband laboratory physical property measurements but rather estimate an IP time-constant whose characteristic frequency (inverse of IP time constant) lies within the bandwidth of the TEM system used. Further, the chargeability estimate will reflect only that fraction of polarizable material whose response is within the bandwidth of the system.

SummaryThe magnitudes of earthquakes are generally described by an empirical relation called the Gutenberg-Richter law. This relation corresponds to a well-known statistical distribution, i.e. the exponential distribution. In this work, we verify the validity of the Gutenberg-Richter law using a 44-year-long worldwide seismic catalog of strong (Mw ≥ 6.5) events, by testing the exponentiality and the independence of the magnitudes. Moreover, we suggest a new way to visualize the distribution of the magnitudes, which complements the classical magnitude frequency distribution plot.

SummaryThe Cole-Cole (CC) equation was introduced as an empirical formula that in many cases matched laboratory measurements of dispersive frequency-dependent dielectric, conductive or resistive physical properties. The CC formula has four parameters, a high and low frequency limit, one of which is often normalised into a polarizability η or a chargeability m, a time constant τ and a frequency dependence c. The chargeability can be directly related to the volume fraction of polarizable material in a uniform conductive background. In an electrically resistive host, the chargeability is given by the fraction of pores blocked by sulphides or other electronically conductive polarizers in the Pelton conceptual model. Fundamentally, an electrochemical model predicts that uniformly sized, non-reactive polarizers would produce exponential or Debye decays during the return to equilibrium after excitation, while uniformly-sized, chemically-reactive materials such as sulphides should exhibit Warburg decays. Numerically, the frequency dependence c can be matched and predicted from the standard deviation of a log-normal size distribution of polarizable material. Using a simple circuit analogy, excitation by a Voltage or a Current source can produce very different decay time-constants. This analogy and mathematical analysis predict that the time constants from fitting a CC conductivity model (${{\tau }_\sigma }$) can be very different from a CC resistivity model (τρ). A published laboratory sulphide measurement example with ${{\tau }_\rho }$ = 10 days, m = 0.96 and c = 0.2 corresponds to ${{\tau }_\sigma }$ = 88 ms, or 7 orders of magnitude shorter in time. Further, the physical circuit analogy confirms the electrochemical model that while the different time constants are mathematically independent of m and c, ${{\tau }_\sigma }$ is an intrinsic time-constant fundamentally related to the grain size of polarizable material and as such a far better parameter than τρ to use for mineral discrimination studies. Published tabulations of fitted parameters need to specify c as well as the chargeability and CC time-constant to allow for transformation between conductivity and resistivity time-constants. The maximum phase time constant ${{\tau }_\phi } = \sqrt {{{\tau }_\sigma }{{\tau }_\rho }} $ is a fitting rather than experimentally measurable parameter, that is still dependent on m and c, but less dependent than ${{\tau }_\rho }$ as previously discussed in the literature.

SummaryThe northeastern Nam Co region, historically impacted by significant earthquakes such as the 1951 Ms 8.0 Beng Co and 1952 Ms 7.5 Gulu events, exhibits intricate crustal deformation patterns shaped by strike-slip and extensional fault systems. This study integrates multi-source geodetic data to analyze three-dimensional deformation patterns and fault interaction. In this region, we established 16 new GNSS campaign-model observation stations. By aligning these with existing GNSS data within a unified reference frame, we obtained a high-resolution horizontal velocity field. Additionally, ascending and descending InSAR deformation velocity fields were derived utilizing Sentinel-1A data and SBAS-InSAR technology. By fusing the GNSS and InSAR velocity fields, we extracted and analyzed the three-dimensional deformation velocity field. Utilizing an enhanced back-slip dislocation model that accounts for fault dip angles, we inverted the slip behaviors of three major faults and investigated their tectonic transition patterns. The deformation field reveals distinct kinematic behaviors among these faults. Specifically, the Beng Co fault demonstrates dextral strike-slip motion, increasing from 4.8 ± 0.1 mm/yr in west to 5.4 ± 0.1 mm/yr east, accompanied by thrusting at a rate of 3.5 ± 0.1 mm/yr. Notably, the locking depths deepen eastward from 12.6 ± 0.6 km to 17.4 ± 0.8 km. In contrast, the Dong Co and northern Yadong-Gulu faults exhibit sinistral strike-slip rates of 1.0 ± 0.1 mm/yr and 4.2 ± 0.1 mm/yr, respectively, paired with maximum extensional rates of 2.0 ± 0.3 mm/yr and 5.2 ± 0.5 mm/yr. Collectively, these faults form a stable strike-slip to extension coupled system, modulating regional crustal deformation through kinematic interactions. This study quantifies a tectonic transition model, elucidating how strike-slip faulting evolves into extensional structures in central Tibetan Plateau. Our findings contribute to a deeper understanding of strain partitioning and intracontinental deformation mechanisms within the central Tibetan Plateau.

SummarySeismic activity induced by underground engineering projects often involves complex causal mechanisms, and represents significant hazards, including ground subsidence, disruption of surface and underground water systems, ecological damage, structural damage to buildings, and even casualties. Consequently, induced seismicity has become an important topic in the risk assessment and protective measures for underground engineering projects. During the construction of the Hongtu Tunnel on the Dafenghua Expressway in Guangdong, China, a series of earthquakes occurred nearby, with the biggest of magnitude ML = 3.7, alongside significant water inflows at multiple locations. This study analyzed seismic network data from 2017 to 2022 around the tunnel area to investigate the potential relationship between the seismic swarm and tunnel construction and uncover the underlying mechanisms. After velocity model corrections and double-difference relocation, the earthquakes were primarily distributed at depths of 1∼4 km. Three concealed, steeply dipping NE-trending faults, each 3∼7 km in length, were identified based on the earthquake distribution. The swarm began about one month after the onset of water inflows in the tunnel and grew significantly after the peak daily inflow, culminating in the ML 3.7 mainshock. A strong spatiotemporal correlation was observed between the seismic swarm and the water inflows. During the first year of the swarm, the seismicity displayed migration characteristics consistent with pore pressure diffusion, with an initial diffusion depth of approximately 2 km and a diffusion rate of 0.0039∼0.0446 m²/s, and best fit by the classical parabolic diffusion model (α = 0.5). After 2021, the earthquakes occurred more consistently, mainly exhibiting stress-triggering characteristics. Over time, the seismicity gradually extended to greater depths, with focal mechanisms changing from normal faulting to strike-slip faulting. The local stress field shifted from extensional to shear, which reflected the sustained influence of pore pressure diffusion on fault activation. Fluid diffusion not only initially activated the faults but also promoted repeated fault slip during the seismic swarm, indicating that prolonged water inflow significantly altered fault activity patterns and the regional stress field. This study is the first to reveal the phenomenon of long-distance induced seismicity caused by tunnel water inflow and the role of pore pressure diffusion in triggering such events, which offers new insights into the safety of underground construction and the study of fluid-related geological processes.

SummaryThe phenomenon of transient strain from seismic waves triggering earthquakes is a robust observation. However, our understanding as to why seismic waves can trigger earthquakes remains incomplete. In this study, we use particle simulations to investigate the response of sheared granular matter to dynamic strain perturbations in order to better understand the dynamic triggering of earthquakes. In our simulation, an unstable slip is triggered when the dynamic strain above a threshold (critical strain) is applied to the system. We show that the critical strain is of the same order (10–6 to 10–9) as those in some experimental and observational studies. This enhanced response is observed at resonance wavelengths. Resonant vibration decreases the shear modulus of the granular system, and accordingly the shear strength is reduced, leading to unstable slip. This modulus softening is due to the increase in slipping contacts between particles. The relevance of simulation results to natural earthquake faults is discussed as to whether seismic waves can satisfy the resonance condition.

SummaryAs seismic migration is increasingly applied to more and more complex media, more sophisticated imaging techniques are required to generate accurate images of the subsurface. Currently, the best results for imaging are achieved by Least-Squares Migration (LSM) methods, such as Least-Squares Reverse Time Migration (LS-RTM) and Full-Wavefield Migration (FWM). These methods iteratively update the image to minimize the misfit between the forward modelled wavefield and the recorded data at the surface. However, a key challenge for these techniques is the speed of convergence. To accelerate the speed of convergence, pre-conditioning is commonly applied. The most common preconditioner is the reciprocal of the Hessian operator. However, this operator is computationally expensive to calculate, making it difficult to apply directly. In this paper, we present a novel, alternative, preconditioner for FWM. This preconditioner is based on applying Galerkin projections to a linear system, which projects the system onto a set of known basis vectors. To find an appropriate set of basis vectors for this approach we apply Proper Orthogonal Decomposition (POD) to a set of partial solutions of the linear system. The resulting method gives an approximation to the pseudo-inverse based on these basis vectors. To test this technique, which we name Model-Order Reduced FWM (MOR-FWM), we apply it to the synthetic Marmousi model as well as to field data from the Vøring basin in Norway. For these examples, we show that MOR-FWM yields an improved data-misfit compared to the standard FWM approach. In addition, we show that the result for the field data case can be improved by normalizing the partial solutions before applying POD.

SummarySeismic interferometry, applied to continuous seismic records, yields correlation wavefields that can be exploited for information about Earth’s subsurface. The coda of the correlation wavefield has been described as multiply scattered waves that are highly sensitive to crustal heterogeneity and its changes. Therefore, the coda of consecutive correlation wavefields allows to monitor velocity variations over time to detect weak changes in the medium at depth. Ocean microseisms, generated by ocean-land interactions, are the dominant continuous source of seismic energy at frequencies below 0.5 Hz. It is well-understood that these oceanic sources are not homogeneously distributed over Earth and change over the seasons, which commonly results in asymmetric correlation wavefields from seismic data. The impact of these seasonal changes on the coda of the correlation wavefield is typically considered negligible. In contrast, we demonstrate that oceanic noise sources and their changes directly impact the composition of the coda. We compute correlation wavefields between several master stations throughout Europe and the Gräfenberg array in Germany. We beamform these correlation wavefields, in the microseism frequency band, to detect coherent waves arriving at the Gräfenberg array. We perform this analysis for a two-year period, which enables us to compare variations in source direction over the seasons. We find seismic waves arriving from dominant sources to the North-Northwest of Gräfenberg in boreal winter (with slownesses corresponding to surface waves) and towards the South in summer (with slownesses corresponding to body waves) throughout the entire correlation wavefield, including its late coda. Beamforming the original recordings before cross-correlation confirms that the seasonally dominant source regions are directly detected also in the correlation wavefield coda. We derive that seismic waves propagating from isolated microseism source regions will be present in correlation wavefields even if the master station, or ”virtual source”, used for correlation recorded no physical signal at all. The findings we present raise concerns about velocity monitoring approaches relying on the coda being comprised exclusively of scattered waves. Our results also suggest that higher-order correlations do not achieve an effectively more homogeneous source distribution, and instead may even enhance such bias.