Geophysical Journal International

SummaryLow-frequency tremors provide important constraints on the mechanical and hydrological conditions of subduction plate interfaces, yet their focal depth distributions remain challenging to be resolved compared to those of fast earthquakes. In this study, low-frequency tremors in the Hyuga-nada region are independently relocated using earthquake location techniques with travel-times by converting cross-correlation-derived arrival-time delays into pseudo-arrival times for each single station, without imposing a priori depth constraints. The resulting hypocenter distribution reveals systematic spatial patterns that are robust across multiple observation periods. Hypocenter locations are determined through a multi-stage relocation procedure, in which an initial maximum-likelihood location incorporating station correction terms is followed by a double-difference relocation using a three-dimensional velocity model. Only events meeting stringent quality thresholds on converted arrival counts, azimuthal gap, and initial hypocentral depth uncertainty are retained. Location uncertainties are quantified at each stage using a Bayesian framework and jackknife resampling, which respectively characterize the probabilistic confidence of the solution and its robustness against the influence of individual stations and potentially biased converted arrival times. These steps are iteratively applied to refine source parameters, thereby improving depth resolution and overall location stability. The relocated low-frequency tremors are predominantly concentrated at shallower depths than fast earthquakes and are distributed along the plate interface. Comparison with seismic reflection profiles indicates that low-frequency tremors are generally distributed in the vicinity of strong reflective interfaces along the plate boundary, supporting previous interpretations that associate low-frequency tremors with mechanically weak, fluid-rich layers. The relocated tremor locations further show a clear spatial correspondence with a subducted seamount inferred from magnetic anomaly data, consistent with mechanical and hydrological contrasts associated with seamount subduction, particularly along the downdip flank, which may reflect stronger compression and enhanced fluid drainage. These results demonstrate that low-frequency tremors in the Hyuga-nada region are strongly controlled by the plate interface and structural heterogeneity. The relocation framework presented here provides a unified basis for integrated analyses of fast and slow earthquake processes in subduction zones.

SummaryWave-induced fluid flow (WIFF) across multiple scales causes significant seismic wave dispersion and attenuation over a broad frequency range in partially saturated rocks. Accurately characterizing these effects is crucial for improving hydrocarbon reservoir detection. In this study, we develop a model that integrates mesoscopic Rayleigh bubble oscillation with an extended microscopic squirt-flow mechanism (with varying microcrack aspect ratios), including capillarity effects. Plane-wave analysis indicates that both P1 (fast P) and S waves are strongly influenced by multiscale fluid flow, with the extended squirt-flow mechanism producing a broader range of wave attenuation and dispersion at high-frequency. Numerical analysis shows that the pore structure controls wave-velocity dispersion and attenuation at different scales. The inclusion radius determines the mesoscopic relaxation characteristic frequency, and the microcrack porosity governs squirt flow. Furthermore, the model predictions agree reasonably well with two independent sets of experimental data, validating the applicability of the theory. This study provides valuable insights for understanding seismic wave dispersion and attenuation in partially saturated rocks.

SummaryInterferometric Synthetic Aperture Radar (InSAR) deformation measurements are susceptible to long wavelength errors, yet existing geometric and strain models for estimating three-dimensional (3D) deformation from combined Global Navigation Satellite System (GNSS) and InSAR observations do not explicitly account for these errors. Neglecting long wavelength errors can substantially bias 3D deformation estimates for large earthquakes. We therefore develop a combined adjustment model that introduces additional systematic parameters to represent InSAR long wavelength errors, defines these parameters using GNSS minimum constraints, and optimizes the stochastic model of the heterogeneous observations through least-squares variance component estimation (LS-VCE). We evaluate the model using synthetic tests and observations of the 2011 Tohoku-Oki earthquake. In the synthetic tests, the proposed model constrained the datum deviation to below 1 mm, whereas the corresponding mean datum deviations obtained using the geometric and strain models were 5.40 and 3.80 m, respectively. Relative to the design values, compared with the geometric and strain models, the proposed model reduces the mean absolute error (MAE) for the East (E) and North (N) components by 98.63% and 98.84%, 93.65% and 96.45%, respectively. This indicates that the proposed model can effectively eliminate the interference of long wavelength errors on the 3D deformation parameter estimation. We further applied the proposed model to the 2011 Tohoku-Oki earthquake and compared the resulting deformation with the Jet Propulsion Laboratory (JPL) co-seismic displacement solution. Compared with the geometric and strain models, the proposed model reduces the MAE for the E, N and U components by 89.32% and 89.27%, 83.33% and 85.25%, 37.50% and 82.14%, respectively. Under the adopted fault geometric parameters, the proposed model also recovered candidate strike-slip components that were not resolved by the conventional models. Compared with the geometric and strain models, the proposed model reduces the MAE between forward modeled and JPL co seismic displacements for the E and N components by 98.36% and 98.01%, 95.45% and 94.59%, respectively. Overall, explicitly estimating InSAR long wavelength errors within a combined adjustment framework improves the numerical and physical reliability of GNSS-InSAR 3D deformation estimation.

SummarySeismic velocity profiles are crucial for exploring subsurface features and monitoring seismicity. The depth sensitivity of seismic surface waves depends on frequency and, therefore, dispersion curves, routinely obtained from frequency-wavenumber methods, are used to invert subsurface velocity structures. However, an optimal, scalable array design concept suitable for various target depths and applications remains under discussion. We propose an array design, optimised for passive seismic surveys, based on circular arrays using only the target depth range as input. Exploiting the impact of the array layout on wavenumber resolution and aliasing, we derive new relations for optimized arrays requiring only the target depth range. The concept accounts for different modes and wave types by considering penetration depth in terms of wavelength. For fundamental-mode Rayleigh waves, the maximum inter-station distance should be at least five times the maximum target depth, and the minimum distance smaller than twice the minimum target depth to avoid aliasing.

SummaryGlobal Navigation Satellite System (GNSS) and Gravity Recovery and Climate Experiment (GRACE)/GRACE Follow-On (GFO) observations provide complementary geodetic constraints on terrestrial water storage (TWS) changes, but either data source alone is often insufficient for interpreting water storage response in small humid coastal monsoon regions. We integrated vertical displacements from 49 GNSS stations with GRACE/GFO mascon solutions to estimate monthly TWS changes in Guangdong, China, from 2017 to 2024. An objective-space weighting strategy was used to balance local detail constrained mainly by GNSS with regional mass signals stabilized by GRACE/GFO. The joint solution retains GNSS-supported spatial detail not resolved by GRACE/GFO and reduces boundary artifacts relative to the GNSS-only inversion. It also shows the highest consistency with hydroclimatic water-balance estimates among the geodetic solutions, with a maximum correlation coefficient of 0.85 (95 per cent CI: 0.72–0.93). We further combined the joint solution with hydroclimatic data to interpret water storage response under different moisture conditions. At the seasonal scale, soil moisture storage (SMS) responds nearly synchronously to precipitation, whereas TWS changes peak about one month later. Both responses show weaker amplitudes and shorter lags along the coast, with the opposite pattern inland. For storage anomalies, TWS shows state-dependent memory, mainly 3–7 months under normal and wet conditions but 11–12 months during the prolonged 2021 drought, while SMS remains tied to a much shorter timescale. The joint TWS record also captures three hydrological droughts and two regional flood episodes from 2017 to 2024, including an extreme drought that persisted from December 2020 to January 2022. Together, these results show that joint inversion of GNSS and GRACE/GFO observations provides observation-based constraints on integrated TWS for interpreting storage adjustment, state-dependent storage memory and major hydrological extremes in humid coastal monsoon regions, beyond what can be inferred from meteorological indicators or model-based analyses alone.

SummaryNatural fractured reservoirs or artificially hydraulically fractured networks typically exhibit hierarchical microstructures characterized by porous backgrounds permeated by microcracks and intersected by macroscopic fractures. Prior theoretical frameworks predominantly decouple micro- and macro-scale heterogeneities, thereby overlooking the interplay between pore-scale squirt-flow dissipation and fracture-induced wave attenuation. By integrating microcracked porous wave-motion theory with fracture scattering theory, we present a unified model for characterizing wave propagation in fluid-saturated microcracked porous rocks with intersecting fractures. Five distinct attenuation mechanisms are concurrently incorporated: squirt flow at the microcrack scale, wave-induced fluid flow (WIFF) at both fracture-background and fracture-fracture mesoscales, elastic scattering, and Biot's global flow. The results show that microcrack density primarily affects matrix stiffness and shifts the attenuation of squirt flow to lower frequencies. In contrast, microcrack aspect ratios control the characteristic frequency and its coupling with fracture-related mechanisms. Fluid mobility governs the frequency range of WIFF and Biot flow, with high mobility potentially creating attenuation gaps where scattering-induced velocity drops disappear. The model also accommodates multiple-aspect-ratio microcrack systems and extends to planar fractures as linear-slip interfaces, enabling realistic characterization of anisotropic shale reservoirs. This unified approach provides a comprehensive basis for interpreting seismic dispersion and attenuation in complex fractured formations relevant to hydrocarbon exploration and reservoir monitoring. In parallel fractures and zero microcrack-limiting cases, our model predictions agree well with those of existing models.

SummaryRiver embankments are increasingly afflicted by intense and frequent flood events owing to climate change. The characterization and monitoring of river embankments are hence becoming crucial for climate change adaptation. Electrical resistivity tomography (ERT) is one of the most common geophysical techniques for the characterization of river embankments to identify local anomalies and potentially prevent damages during floods. ERT surveys are usually executed over the embankment crest and interpreted along a longitudinal section. However, 2D interpretation of ERT data can be strongly affected by the embankment morphology along the investigated portion (i.e., the shape of the longitudinal section and the presence of curves). To overcome this issue, we propose 3D inversion of ERT data acquired over embankments, and compare 2D and 3D inversions of both synthetic and field ERT data in order to understand the differences in results between the two approaches. We created a 3D synthetic model with several types of embedded anomalies to realistically represent the embankment geometry. We considered both a straight and a curved longitudinal profile of the embankment. Synthetic data were then generated with two tailored electrode-configuration arrays (i.e., dipole‒dipole and Wenner‒Schlumberger). 2D field data were acquired close to Turin (NW Italy) over an embankment that has often been damaged by floods. The results from 2D and 3D inversions were analyzed and compared. The main outcomes show that, with synthetic and field ERT data from the crest of the embankment, the subsurface was better characterized by means of 3D inversion, especially in presence of 3D eccentric buried features or complex/curved model domains. From the examples of synthetic data, the presence of the embankment slopes and curvatures negatively affected 2D inversion results, which presented more artefacts and higher misfit values than the 3D inversion results. From the field case study, the 3D resistivity volume presented good sensitivity in the transversal direction, thus allowing useful inspection of feature discontinuities or resistivity contrasts. Moreover, for the synthetic data the 3D results were more coherent than the 2D results with the true resistivity and position of the anomalies. For the field data, the 3D results were more coherent than the 2D results with known information from geological and geotechnical data. The more physically consistent 3D solutions may be further exploited for quantitative analysis on the petrophysical properties of the embankment material. In conclusion, 3D ERT inversion is recommended in presence of irregular domains (slopes, curves, strong discontinuities) even with 2D data.

SummaryImplicit full-waveform inversion (IFWI) has shown strong potential for seismic velocity reconstruction due to the implicit regularization of neural networks. However, conventional IFWI based on direct modeling represents multi-scale structures within a unified parameter space, where gradients associated with different wavenumber components may interfere, leading to slow convergence and limited recovery of deep structures. To address this issue, we propose a perturbation-based IFWI framework that decomposes the velocity model into a fixed background and a learnable perturbation. This parameterization reduces multi-scale gradient interference and improves optimization efficiency by restricting updates to the perturbation subspace. In addition, a coordinate normalization strategy is introduced to decouple spatial scale from the network frequency response, enabling stable optimization across different model sizes. Numerical experiments demonstrate that the proposed method achieves more stable convergence, improved structural resolution, and more reliable deep recovery compared to direct modeling IFWI. The approach also shows strong robustness under challenging conditions, including inaccurate initial models, varying spatial scales, and missing low-frequency data.

SummaryWe introduce a new 1D model of shear attenuation in the mantle (QM1D) constructed from measurements of differential attenuation of S and ScS body waves. We build a dataset of over 90,000 measurements of $\delta t^*_{\text{ScS-S}}$ from seismograms with high-quality S and ScS waveforms from deep earthquakes between 20–75° epicentral distance from 1990–2024. To ensure that measurements are robust, we require that values of differential attenuation are similar between the instantaneous frequency matching and waveform matching methods. Using this dataset, we perform a Bayesian inversion to construct 1D profiles of mantle shear quality factor Qμ and its uncertainties. Our dataset has a high sensitivity to the mid-mantle (660–2000 km) corresponding to the turning depths of the S phases, which maps into lower uncertainty in the ensemble distribution at these depths. QM1D suggests that the mantle is more attenuating than PREM on average; but, notably, that there is a strongly attenuating layer around 750 km depth and that there are weakly attenuating layers around 900 km and 2000 km depth. Due to the source-receiver distribution for the phases, our dataset has a sampling bias in the mantle around the edge of the Pacific which may not be representative of the whole mantle. The weakly attenuating zones in the mantle may correspond to a viscosity jump around 900 km depth and a further viscosity jump in the deep mantle around 2000 km depth, whilst the thin strongly attenuating layer may correspond to regions of slab stagnation caused by a potential thin layer of low viscosity beneath the mantle transition zone in the circum-Pacific region.

SummaryWe report the first explicit derivation of the Rikitake dynamo as a direct approximation of mean-field magnetohydrodynamics (MHD), establishing a systematic connection between a classical conceptual reversal model and the governing large-scale field equations. Starting from the mean-field induction and momentum equations, we show that a strongly truncated Galerkin expansion already retains the nonlinear couplings among magnetic induction, diffusion, convective driving, and Lorentz-force feedback required for reversal dynamics. Within this framework, a family of low-dimensional models arise as a limiting approximation of the reduced mean-field system. Placing further constraints on the Galerkin basis functions results in what we call the canonical Rikitake model. We provide an explicit harmonic realization that fulfills these constraints and preserves the minimal nonlinear structure capable of producing irregular oscillations and spontaneous polarity reversals. The reduced parameters acquire clear physical meaning in terms of inductive regeneration, effective magnetic diffusivity, forcing, and magnetic backreaction. Using representative outer-core scalings, we obtain characteristic amplitudes and timescales consistent with geophysical estimates, including observed chron durations. Numerical integrations reproduce several features of geomagnetic variability, including polarity reversals, irregular oscillations, and extended intervals of persistent polarity. The canonical model supports both short and long chrons, including superchrons, whereas the harmonic realization yields only short chrons, indicating that richer basis functions may be needed to recover broader reversal statistics. These results place Rikitake-type models on firm physical grounds and provide a route toward systematically constructed low-dimensional dynamos derived from MHD.

SummaryThis paper has analyzed the spatial and temporal distribution of the local seismic shear-wave splitting parameters at the source region of the Nepal Gorkha MW7.8 earthquake in 2015. There are significant spatial variations of valid results within the aftershock zone. The fast directions change abruptly from west to east, revealing heterogeneity in the physical property or geometry of the upper crustal medium along the strike of the seismogenic fault. Combined with seismic tomography and magnetotelluric sounding, it is inferred that the variation of the upper crust from west to east in this region may represent a response to the tearing of the Indian lithospheric plate. The aftershock area shows large time delays, with their average results decreasing southward, indicating that a large amount of stress accumulated in the upper crust along the seismogenic area during the seismogenic process of the Gorkha earthquake, with the stress gradually weakening from north to south. Furthermore, time delays in the dense aftershock zone remained highly unstable for nearly one year after the Gorkha earthquake, reflecting that the stress accumulated during the seismogenic process of this earthquake failed to fully release with the mainshock and strong aftershocks, and the stress still adjusted violently.

SummaryGeothermal heat flow (GHF) provides a first-order constraint on lithospheric thermal structure, yet continental-scale prediction remains challenging due to sparse, uneven, and strongly nonstationary observations. Most existing machine-learning approaches treat GHF measurements as spatially independent, limiting their ability to capture geological coherence across tectonic provinces. Here, we develop a spatial correlation recurrent neural network (SCRNN) that explicitly incorporates spatial dependencies among multiple variables to predict continental-scale GHF. By employing SCRNN, we produce a high-resolution GHF distribution of the contiguous United States. Higher GHF values are concentrated in tectonically active western regions characterized by crustal thinning, lithospheric extension, and fault systems. Lower GHF values are in the geologically stable central and eastern Craton of North America. SHapley Additive exPlanations analysis indicates that lithosphere-asthenosphere boundary depth, Moho depth, Curie-point depth, and fault systems exert first-order controls on GHF patterns. These results demonstrate that SCRNN captures meaningful controls on continental GHF.

SummarySeismic anisotropy, the dependence of seismic wave speeds on the direction of propagation and/or polarization, places crucial constraints on deformation and convective flow in the lowermost mantle (the D″ layer). Shear waves that diffract along the core-mantle boundary (CMB) are ideally suited for probing this region due to their long horizontal ray paths in the lowermost mantle. However, the inference of D″ anisotropy presents multiple methodological challenges, including the need for accurate upper mantle corrections and the potential for apparent splitting produced by isotropic structures. In this study, we leverage global wavefield simulations to develop an automated measurement procedure that excludes Sdiff waves that are strongly SVdiff-polarized, which can lead to apparent splitting even in the absence of seismic anisotropy. We apply this workflow to a global dataset of ~20 million three-component seismograms. By adhering to strict quality control criteria, we generate a highly regionally consistent dataset of high-confidence measurements. Overall, we obtain ~4, 300 Sdiff splitting measurements that constrain the presence or absence of D″ anisotropy. We project these measurements onto the D″ layer and report Sdiff splitting across 237 global geographic bins (approximately 28 per cent of the D″ layer by surface area) with sufficient data coverage. Of these, 55 bins exhibit clear splitting, many of which are located beneath the northern Pacific. Furthermore, we characterize the directional dependence of splitting using a backazimuth binning algorithm for regions with multi-directional coverage. We make this curated dataset of geographically and directionally binned Sdiff splitting measurements, along with all codes necessary to reproduce it, publicly available.

SummaryThe mechanical properties of Earth’s uppermost layers are seismically relevant because they can modify seismic waves as they pass through them. These properties may significantly modify the intensity of ground shaking during an earthquake. The relatively common presence of fluids in these layers raises the question of whether they can affect ground motion during earthquakes. In general, seismic ground motion results from the complex interaction of P- and S-waves within the subsurface, and it is certainly influenced by the presence of fluids. This work examines the effects of varying the groundwater level on seismic ground motion by studying the propagation of P- and S-waves through a saturated porous layer over an elastic half-space. Our results show that ground motion due to both types of incident waves is affected by the presence of water, with changes in amplitude and resonance frequency. These results highlight the importance of accounting for fluid saturation in both site characterization and site response.

SummaryThe potential for geothermal resources and the tectonic activity of the north-central region of the Gulf of California, Mexico (GoC) offer a good opportunity to investigate the lithospheric structure of this region. We used one year of ambient seismic noise and 11 shallow earthquakes recorded by the Gulf of California Broadband Seismological Network (RESBAN) and stations from the Northwest Mexico Seismic Network (RESNOM) and the National Seismological Service (SSN) to estimate the three-dimensional S-wave velocity structure. We used surface wave travel times at short and long periods, and a two-stage tomography technique to retrieve an S-wave velocity model. We use coarse mesh of 16×40 cells with dimensions of 20×20 km in the horizontal direction and variable up to 60 km in depth. The absolute velocity model shows three areas in depth and significant lateral changes from the north to the center of the GoC associated with crustal and upper-mantle structure. We observed variations in crustal thickness that would indicate a variable Moho depth in the northern GoC and a likely constant depth towards the center of the Gulf. We identified low-velocity anomalies that could confirm the existence of a geothermal system in the north of the Gulf. The velocity models show relative changes associated with a heterogeneous lithosphere from the central region to the north of the GoC, which would indicate a transition from weak to stronger lithosphere. Velocity anomalies explain some of the tectonic differences between the north and the central GoC, as well as geophysical observations, including heat flow measurements and seismicity.

SummaryWe examined source characteristics of the recent earthquake sequence that occurred in February 2025 between the volcanic fields of Santorini, Kolumbo, Anydros, and Amorgos in south-central Aegean. This region lies within the volcanic back-arc domain of the Hellenic subduction zone in the Aegean Sea. During this period, the area experienced intense seismic activity (ten earthquakes 5.0 < Mw < 5.5), with earthquake magnitudes generally increasing towards Amorgos Island. The largest event was the Mw 5.4 earthquake on 10 February 2025, located northeast of Santorini. In this study, we first obtained 4,869 earthquakes that occurred between 2 January 2024 and 23 March 2025. We then derived full moment-tensor (MT) solutions for 30 earthquakes (Mw ≥ 4.5) by inverting seismological datasets. In addition, principal stress orientations were calculated by stress tensor inversion to recognize the most recent crustal stress field in the region. The spatio-temporal distribution of the relocated aftershocks indicates that the seismic activity initiated beneath the Santorini and Kolumbo volcanic centers, and subsequently migrated northeastward towards Anydros and Amorgos. Most earthquakes occurred at shallow crustal depths of approximately 5–15 km. Stress-inversion results and the orientation of T-axes highlight a dominant NW–SE extensional regime that characterizes the combined effects of deep geodynamic processes and shallow crustal dynamics in the back-arc area of the Aegean. The full MT solutions reveal predominantly normal, oblique-normal, and strike-slip faulting mechanisms, accompanied by significant non-double-couple (non-DC) components. These non-DC signatures strongly suggest tensile processes that are likely related to magma intrusion, fluid migration, or hydraulic fracturing. We therefore conclude that the 2025 seismic sequence was largely driven by crack-opening mechanisms and fluid- or magma-induced activity within the crust, consistent with the broader structural framework and volcanic nature of the Santorini–Amorgos region. Overall, this study enhances our knowledge on tectonic characteristics of the 2025 Santorini-Amorgos earthquake activity, which is significant for evaluating and monitoring the potential hazards in the south-central Aegean region.

SummaryJoint inversion of electrical and electromagnetic (E&EM) data can reduce the non-uniqueness of geophysical inversion by exploiting the complementary sensitivities of galvanic and inductive measurements. However, many existing approaches either rely on one-dimensional electromagnetic (1D EM) forward modelling or neglect induced polarization (IP) effects, which can be limiting in the presence of polarizable targets or laterally complex structures. We present a joint inversion framework that combines three-dimensional electromagnetic (3D EM) modelling with two-dimensional direct current resistivity and induced polarization (2D DCIP) modelling. The framework uses a common maximum phase angle (MPA) reparameterization of the Cole-Cole model for both inductive and galvanic data, together with decoupled forward and model meshes that allow each EM system to be modelled on a dedicated 3D forward mesh while all datasets update a common 2D inversion model. The method is evaluated using two synthetic examples and a mineral exploration field dataset. The synthetic examples represent a hydrogeological setting with weak-to-moderate IP effects and a mineral exploration setting with stronger IP effects associated with a finite conductive and polarizable intrusion. For both cases, we compare EM-only inversion, DCIP-only inversion, 1D EM/2D DCIP joint inversion, and 3D EM/2D DCIP joint inversion. The results show that the 3D EM/2D DCIP joint inversion improves the recovery of resistivity and IP structures when 3D EM effects are significant. In the hydrogeological example, the 3D EM kernel improves the imaging of laterally complex clay-rich structures, whereas in the mineral exploration example it reduces artefacts associated with the finite intrusion and improves the separation of distinct polarizable zones. The field application combines airborne EM and ground DCIP data acquired over Fe-Ti-V oxide mineralization in southern Portugal. The joint inversion honours both datasets and provides a single resistivity and IP model compatible with the main geological setting. Overall, the results indicate that 3D EM/2D DCIP joint inversion is useful for interpreting inductive and galvanic datasets across different IP regimes, particularly where finite or laterally complex polarizable structures make 1D EM modelling inadequate.

SummaryGeophysical inversion is inherently non-unique, and efficient uncertainty quantification remains a major challenge for large-scale applications. Conventional Bayesian approaches, such as Markov chain Monte Carlo (McMC), provide rigorous posterior estimates but are often computationally prohibitive, while deterministic inversion methods lack uncertainty characterization. In this study, we develop a Stein variational gradient descent (SVGD) framework for one-dimensional (1-D) geophysical Bayesian inversion. SVGD is a variational inference method that approximates the posterior distribution using an ensemble of interacting particles, which are iteratively updated through gradient-based optimization. The proposed framework is evaluated using airborne electromagnetic (AEM) data in both synthetic and field scenarios. Synthetic tests demonstrate that SVGD effectively recovers subsurface resistivity structures while capturing associated uncertainty. The influence of particle number on posterior resolution and computational cost is analyzed, and comparison with McMC highlights differences in posterior characteristics and computational efficiency. Application to field data from the Mississippi Alluvial Plain shows that the SVGD posterior mean agrees well with deterministic inversion results, while providing the uncertainty estimates. In addition, noise sensitivity tests are conducted to investigate the impact of data errors on inversion behavior. We find that increasing noise level degrades resolution and shifts the posterior toward prior-dominated regimes. Overall, the SVGD framework provides a flexible and scalable approach to probabilistic geophysical inversion, with strong potential for extension to higher-dimensional problems.

SummaryConventional earthquake rupture models in elastic media with scale-independent dissipation predict rupture acceleration up to a limiting speed dictated by seismic wave speeds. However, most earthquakes are inferred to be significantly slower than these theoretical limiting speeds. Previous two-dimensional modeling shows that inelastic deformation of the material surrounding the fault can stabilize rupture propagation at velocities below the limiting speed, owing to scale-dependent dissipation due to the expansion of the off-fault inelastic zone with rupture length. An essential ingredient missing from such models is that large earthquake ruptures become elongated after they saturate the seismogenic width, which prevents a length-dependent growth of the off-fault dissipation zone. Here, we carry out a computational study of elongated ruptures with off-fault plasticity, with a particular focus on steady rupture speeds. We combine a 2.5D (width-averaged) modeling approach and Mohr-Coulomb plasticity to simulate width-limited strike-slip earthquake ruptures propagating in an inelastic medium. We observe steady pulse-like ruptures propagating at speeds below the Rayleigh wave speed. Through a systematic parametric study, we characterize the dependence of the steady rupture speed on the initial fault stress, closeness to plastic failure, and seismogenic zone width. Our simulation results show that off-fault plastic deformation allows the rupture velocity on long faults to saturate at speeds between 6.5 and 99 % of the Rayleigh wave speed, encompassing the common range of speeds of real earthquakes as well as those of tsunami earthquakes. Our results indicate that off-fault plasticity is a viable mechanism for achieving realistic earthquake rupture speeds, even on planar faults with homogeneous stress and strength. The model also provides insights on the factors controlling other earthquake source properties, including pulse width, peak slip velocity, inelastic off-fault deformation and inelastic zone width.

SummaryThe Central Apennines (Italy) host a distributed network of active normal faults capable of generating infrequent but destructive earthquakes with recurrence intervals that often exceed the duration of available historical and paleoseismological records. This observational limitation hampers our understanding of long-term fault behavior, seismic cycle variability, and multi-fault rupture potential. In this study, we apply the MCQsim earthquake-cycle simulator to a three-dimensional fault model comprising 42 major active normal faults in the Central Apennines, generating seven 100 000-year synthetic earthquake catalogs by systematically varying rupture-controlling parameters: fault strength magnitude and heterogeneity, and numerical mesh resolution. We evaluate how these parameters influence seismic productivity, the shape of magnitude–frequency distributions, and the scaling of rupture area and average slip with magnitude. Simulated catalogs are compared with the Italian Parametric Earthquake Catalog (CPTI15) and with empirical magnitude–scaling relationships to assess consistency with observed seismicity. The simulations reproduce the regional magnitude–frequency distribution and yield rupture scaling consistent with empirical relationships, with the strongest sensitivity observed for changes in the difference between static and dynamic friction coefficients and the parameters controlling fault strength heterogeneity. Analysis of interevent times for earthquakes with Mw ≥ 6 shows that recurrence variability is broadly consistent with a quasi-periodic process, while intermittent clustering at multi-centennial timescales is also captured. Although the use of a uniform fault dip and a simplified seismogenic depth boundary, rather than depth-varying fault geometry, may affect the detailed representation of individual ruptures, the overall consistency between synthetic and observed seismicity demonstrates that MCQsim provides a physically grounded framework for investigating long-term fault system behavior and earthquake recurrence variability in the Central Apennines.