Geophysical Journal International

SummaryUnderstanding induced seismicity is a critical challenge for seismic risk mitigation, particularly in Europe where geo-resource exploitation increasingly occurs in densely populated areas. The Lacq site in southwestern France represents one of the most significant and long-lasting cases of induced seismicity in Western Europe. Despite decades of observations, the characterization of Lacq seismicity remains limited by insufficient instrumentation and monitoring, leading to poorly constrained earthquake locations and uncertainties regarding whether seismic events occur within or below the reservoir. In this study, we present a high-resolution seismicity catalog covering a five-year period (2021–2026), representing the most comprehensive dataset available for the Lacq area. To achieve this, a temporary seismic network was deployed through a collaboration between GFZ (Potsdam), OMP (Toulouse), and UPPA (Pau). Event detection and location were performed using two independent approaches: a detailed manual picking and relocation workflow, and a fully automated deep-learning-based method. Combined with the use of a pre-existing local 3D velocity model, this dual approach significantly reduces location uncertainties and improves the robustness and completeness of the catalog. The resulting dataset constitutes a valuable resource for future studies.

SummaryHessian–vector products (HVPs) provide matrix-free access to second-order information in full-waveform inversion (FWI), enabling the use of truncated-Newton updates with improved convergence. This may result in enhanced imaging of small-scale scatterers in strongly heterogeneous, high-contrast media. In the time domain, standard adjoint-based HVP derivations are often presented either via perturbation theory or via a Lagrangian (Lagrange-multiplier) framework. Although these routes are formally equivalent, their correctness hinges on a precise treatment of transpose/adjoint placement, boundary/endpoint terms, and the discrete inner product used by the numerical solver. Revisiting these derivations, we show that a commonly used Lagrangian expression for the first HVP contribution can be incorrect when the wave equation operator L is not symmetric, and that rewriting the first term of the HVP, $\mathcal {H}_1$, using L† (or equivalently, swapping the contraction order) requires accounting for the associated boundary terms, where the L† indicates the adjoint of L. To resolve this ambiguity, we provide a complete Lagrangian formulation that explicitly incorporates initial and boundary conditions, together with a Green-identity-based treatment that distinguishes purely algebraic transpose identities from true operator adjoints obtained by integration by parts. We validate the resulting gradient and HVP expressions against reference HVPs computed by automatic differentiation (AD) using an independent JAX computational graph, which we term FDCompGraph. FDCompGraph reproduces the same discrete time stepping, source/receiver injection, and convolutional perfectly matched layer (CPML) treatment as the PDE solver, enabling apples-to-apples validation after matching wavefields under identical discretization and CPML settings. Using the first-order velocity–stress acoustic vertically transversely isotropic (VTI) system with CPML (a representative case where L is non-symmetric), we demonstrate that incorrect handling of $\mathcal {H}_1$ can yield errors up to approximately 400 per cent for HVP in ε and approximately 40 per cent for HVP in δ, whereas the HVPs in Vp and ρ are numerically unaffected. We further assess numerical errors unique to HVP evaluation, which requires perturbed wavefields and their time/space derivatives. In acoustic VTI tests using the damping sponge boundary condition, we observe HVP sensitivity losses of approximately 10 per cent amplitude compared with the results obtained using the CPML boundary condition. These errors may be difficult to diagnose because the commonly used finite-difference check, H δm ≈ g(m + δm) − g(m), is itself a linear approximation, particularly when the gradient g is computed using wavefields affected by the same imperfect boundary absorption. To make accurate HVPs practical at scale, we develop Devito–Hdm, a storage- and I/O-efficient Devito-based engine for co-simulating the required forward, adjoint, tangent, and second-adjoint system, with CPML-consistent checkpointing and partially on-the-fly correlations. As an application, we perform three-parameter VP–ε–δ acoustic VTI FWI using Newton–CG across models with increasing scatterer density and compare it against standard NLCG and L-BFGS. This application shows that incorporating HVP information overall yields a more reliable recovery of small-scale scatterer perturbations at approximately the same computational cost.

SummaryDip is often read from a regularized electrical resistivity tomography (ERT) image without a parameter-specific uncertainty. Image-resolution measures describe how an inversion blurs the resistivity field, but they do not give a noise-dependent limit for the dip of a geological interface. We derive that limit from the 2.5-D forward operator using Fisher information and a Cramér–Rao bound (CRB). The background resistivity is profiled out in log-data space. For a specified structural family, the dip standard error satisfies $\sigma _\theta \ge \sigma /\left\Vert \partial g_0/\partial \theta \right\Vert _{P_{C_0}}$. We also derive a marginal bound for unknown geometry and a first-order bias caused by structural misspecification. Tests with a 30-electrode dipole–dipole array show that the maximum-likelihood estimate attains the CRB when the fitted family is correct. At 3 % noise, the single-parameter floor is 0.14–0.50○ and the geometry-marginalised floor is 0.33–0.72○. Holding a coupled geometry parameter one cell away from its true value instead produces a 3–5○ dip bias. This is 6–25 times the formal floor. An ablation attributes most of the bias to thickness, width, or depth rather than resistivity contrast. The result is conditional on the chosen model family: the CRB measures noise-limited precision within that family, whereas the bias measures the cost of using the wrong geometry. For the cases tested here, improving an independent geometric constraint is therefore more useful than further reducing data noise.

SummaryThe essence of earthquake early warning lies in the rapid and reliable prediction of the ground motion field after an earthquake occurs. While current numerical ground motion prediction methods based on data assimilation avoid the challenge of estimating source parameters, their interpolation processes do not adequately consider the source rupture directivity effects—a key physical factor controlling the spatially heterogeneous pattern of ground motions, resulting in limited prediction accuracy. To address this limitation, we break the traditional framework of numerical ground motion prediction methods and propose a novel real-time method that couples the source and ground motion field. This method incorporates a line-source model, dynamically extracts rupture directivity information, and estimates the line-source strike by comparing observed and predicted ground motion fields in real time, thereby continuously refining subsequent ground motion predictions. Tests using the 2016 Kumamoto earthquake (Mw 7.0) as an example demonstrate that the rupture direction inferred by this method within the early post-event period (within 28s) aligns closely with the results from detailed post-event inversion. After correction for the directivity effects, the accuracy of ground motion prediction is superior to that of traditional numerical prediction methods. The MAE of peak intensity predictions is reduced by 16.1%, 5.4%, and 19.6% for predictions of 5s, 10s, and 20s in advance, respectively.

SummaryGravity curvature is defined as the third-order derivatives of the gravitational potential. It has been proposed as a potential observable for future satellite gradiometry due to its enhanced sensitivity to short-wavelength components of Earth’s gravity field. Previous studies reported that gravity curvature exhibits apparently weaker spatial correlation with Moho geometry than gravity and gravity-gradient observations, raising concerns about its applicability for Moho inversion. In this study, we systematically assess the capability of gravity-curvature data for Moho recovery, with particular emphasis on the contrasting behaviors of spatial- and spectral-domain inversion schemes. The nonlinear spatial inversion is implemented using the iterative regularization, whereas the spectral inversion follows a spherical harmonic Parker-Oldenburg formulation. Synthetic and real-data tests reveal a clear methodological dependence. Spectral inversion demonstrates robust and stable performance, as it primarily exploits the inherent spectral content of the gravitational field and is insensitive to altered spatial patterns. In contrast, spatial inversion is strongly influenced by data-model correlation, with performance improving at higher observation altitudes but degrading under near-surface scenario. At satellite altitudes, with the effect of iterative stabilization, the spatial inversion yields roughly consistent performance to the spectral inversion. These findings highlight that the inversion capability of gravity curvature depends critically on the inversion framework. When spatial distribution and spectral content are properly considered, gravity curvature can provide a promising constraint for density interface inversion by enhancing sensitivity to specific spatial scales.

SummarySeismic imaging of the Earth’s interior has steadily improved in the past 40 years. Full-waveform inversion (FWI) has been the most recent driver of global and local tomography. Although mineral physics experiments link seismic properties to temperature, pressure and composition of the mantle, direct integration of seismic observations and thermodynamic calculations has not yet been achieved. This work illustrates a framework that integrates FWI techniques with mineral physics to directly invert seismic waveforms for mantle temperature and chemical composition. The forward model is based on interpolation of pre-calculated tables of mantle properties, namely density and P- and S-wave speeds. The inverse problem is formulated as an optimization problem. The gradient of the objective function with respect to temperature and composition is computed using the chain rule which allows us to split it into independent seismic and thermodynamic parts. For the seismic part, the gradient is computed with the adjoint method. For the thermodynamic part, gradients of the pre-calculated tables are obtained by analytically deriving the interpolation function, while accounting for the pressure–density coupling, under the assumption of hydrostatic loading. As a proof of concept, we applied this approach on 2-D synthetic models containing thermal and compositional anomalies where P–SV elastic wave propagation is simulated. The optimization problem is solved using the L-BFGS algorithm. We mitigate the trade-off between the effect of temperature and composition by applying a linear transformation to the model parameters before carrying out the minimization and then transforming back once a solution is found. The proposed framework enables the integration of diverse geophysical datasets as well as the incorporation of additional information on the potential origin of certain mantle anomalies based on petrological constraints, which are crucial to tackle the non-uniqueness of the inverse problem.

SummaryInterfacial slip at fluid–solid boundaries can significantly modify hydraulic transport in tight porous media, yet its implications for representative wave-induced fluid-flow mechanisms have not been systematically investigated. In this study, we develop a physically explicit multiscale poroelastic framework to investigate how pore-scale interfacial slip propagates through representative poroelastic relaxation mechanisms. Starting from oscillatory viscous flow in a cylindrical pore with interfacial slip boundary conditions, we derive a slip-modified pore-scale transport model that provides a physically explicit realization of hydraulic transport under slip conditions. Rather than introducing a new wave-induced attenuation mechanism, the derived transport model is consistently incorporated into the Biot, Biot–squirt, and White–squirt formulations to systematically examine how the same pore-scale transport modification influences viscous flow, squirt-flow relaxation, and mesoscopic pressure diffusion. The results show that interfacial slip systematically reduces viscous resistance, enhances dynamic permeability, and primarily manifests as a shift of the characteristic frequencies governing wave-induced fluid flow. Although the same transport modification is introduced into each model, its macroscopic manifestation depends on the governing relaxation mechanism, leading to distinct responses in the Biot, Biot–squirt, and White–squirt frameworks. Comparisons with laboratory measurements on tight carbonate samples demonstrate that the proposed framework reproduces the principal trends of seismic dispersion and attenuation while maintaining physically reasonable model parameters. These results establish a physically explicit connection between pore-scale interfacial transport and multiscale seismic wave propagation, providing a physically consistent basis for incorporating interfacial effects into poroelastic wave-propagation models and for improving the interpretation of seismic responses in tight porous media.

SummaryThe acoustic approximation of seismic data is a widely used approach in marine seismic data processing. Continental shelf margins typically exhibit inclined interfaces. However, existing methods for calculating the P-wave wavefield beneath inclined seabed commonly rely on horizontal assumptions, which introduce errors into seismic data. In this study, an analytical recursive solution for P-wave propagation in inclined seabed formations is derived by establishing boundary condition equations at arbitrarily inclined interfaces using three-dimensional coordinate transformation. By applying three-dimensional coordinate transformation, boundary condition equations are established at arbitrarily inclined interfaces. Based on the acoustic seismic wave differential equations, a set of recursive equations is developed that enables the analytical computation of P-wave propagation in environments with inclined interfaces. This proposed analytical framework provides an effective approach for resolving wavefields in inclined multilayered formations.

SummaryThe Mount Isa Province in northern Australia represents one of the largest exposed Proterozoic crustal regions within the Australian continent and is recognized globally for its world-class mineral deposits. The underlying crustal architecture has been investigated using a series of geophysical profiles informed by geological mapping, gravity, magnetic, magnetotelluric, and seismic data. In this study, we reanalyze legacy passive seismic data recorded across the entire province to generate a high-resolution crustal shear wave velocity model using ambient noise tomography. These results are interpreted alongside reflection seismic profiles acquired along the same survey lines to provide complementary constraints on crustal composition, also with an association of surface geology. Our model reveals distinct velocity contrasts at some major crustal boundaries, including the Gidyea Suture Zone and several major fault zones. Mapped upper crustal intrusive bodies based on seismic reflections appear to cluster around high-velocity anomalies at 5 km depth in our model, and they seem to be linked to the broader higher velocities within the Mount Isa Inlier at 15–20 km depth. These fast, mid-to-lower-crust features probably represent mafic intrusions associated with a trans-lithospheric magmatic system that also involves felsic volcanism and associated metallogenic activity across the broader region. The research results of this area show how passive seismology can be integrated with other geophysical datasets to improve our understanding of fundamental structures for large mineral deposit systems worldwide.

SummaryWe present a comprehensive numerical investigation of the sensitivity of horizontaltovertical spectral ratio (H/V) measurements of ambient seismic noise to the spatial distribution of noise sources and to strong lateral structural contrasts. Our approach relies on the computation of realistic 3D viscoelastic Green’s functions between a dense set of potential surface sources and a few surface receivers, allowing the synthesis of ambient noise wavefields corresponding to arbitrary spatiotemporal source distributions. The study focuses on a simplified 3D representation of the Mygdonian Basin (northern Greece), where significant temporal variations of H/V amplitude, polarisation, and peak frequency have been observed at a site located near a blind fault. Comparisons between the synthetic H/V and the predictions of the elastic Diffuse Field Theory (DFT) show good agreement under isotropic source conditions, confirming the usefulness of DFT as a reference framework in diffuse or quasidiffuse regimes. However, attenuation, lateral heterogeneities, and uneven source distributions lead to substantial deviations from DFT predictions. We find that sourcerelated effects can dominate both amplitude and polarisation patterns, sometimes masking the structural signature of the fault. Nevertheless, the simulations reveal robust spatial trends of H/V polarisation relative to the fault geometry, suggesting that polarisation mapping may help identify lateral heterogeneities when combined with dense, longduration measurements. Our framework provides a physically consistent way to test the limits of the diffuse field assumption and to explore the applicability of ambient noise–based observables in complex geological environments.

SummaryNear-surface geological faults play a crucial role in subsurface systems because they influence fluid flow and transport and can cause geotechnical failure. Topography can sometimes indicate the presence of fault zones, particularly in active tectonic settings where geomorphic features such as triangular facets are preserved. However, particularly for older inactive faults, characteristic surface geomorphic signatures may have been completely removed by erosion, or the fault may be buried beneath younger sediments. Although direct observations from outcrops or borehole data can confirm the presence of a fault very reliably, such data are commonly too sparse and/or distant to precisely constrain fault location in the target area. In this context, near-surface geophysical methods, such as electrical resistivity tomography (ERT) and seismic refraction tomography (SRT), may be useful to locate faults. However, these geophysical data are typically inverted using smoothness-constrained approaches, which can obscure sharp structural features unless the latter are explicitly accounted for in the models. In this study, we present a workflow that enhances the identification of fault zones by jointly inverting collocated ERT and SRT data, considering a parameterization approach based on a layered subsurface model explicitly including a fault. Unlike conventional approaches, we run multiple inversions in which the starting models are initialized with randomly generated fault geometries, guiding the search toward different structural scenarios that yield comparable data misfits. In the absence of additional prior information, these solutions should be interpreted as equally plausible representations of the subsurface. The resulting ensemble of models provides a practical indication of where the structure is well constrained (i.e. features that appear consistently across solutions) and where ambiguity remains (i.e. features that vary among solutions). We showcase the applicability of our workflow with synthetic and field data, where collocated ERT and SRT measurements were collected across expected fault zones. Our results show that the joint inversion of ERT and SRT data consistently narrows the range of plausible fault locations, particularly in the shallow subsurface. Iterative incorporation of prior information (e.g., one vs. two boreholes) leads to significant improvements in model geometry and physical property estimates. However, consistent with the intrinsic ambiguity of geophysical inversion, the results remain non-unique; even models depicting opposite fault types (normal vs. reverse) may fit the data comparably well. This outcome highlights the strength of the workflow in identifying competing structural scenarios that are equally consistent with the observations. Thus, the resulting ensemble provides essential guidance for targeted ground-truthing and follow-up investigations aimed at resolving the remaining ambiguity.

SummaryTo address the challenge of assessing the reliability of three-dimensional (3D) magnetotelluric (MT) inversion results, we have developed a variational inference (VI) inversion framework (VI-MT) based on the Stein Variational Gradient Descent (SVGD) method. Through parallel particle optimization, we efficiently approximate the model posterior distribution, overcoming the computational limitation of traditional Markov Chain Monte Carlo (MCMC) methods. Synthetic tests show that the VI-MT inversion can effectively recover the synthetic model while reducing the tailing effect, and it can approximate the multimodal posterior distribution, providing quantitative uncertainty estimates. Furthermore, the VI-MT inversion is applied to the field MT data collected at the Weishan volcano in northeast China, with the posterior mean model consistent with the deterministic inversion model. A clear low resistivity body of ∼5 Ω·m with small uncertainties is imaged at depths of ∼2-6 km beneath the volcanic crater, suggesting the existence of a shallow magma chamber. Our study shows that the VI-MT inversion based on the SVGD method can efficiently solve the 3D MT Bayesian inversion, providing reliable model uncertainties.

SummaryCharacterizing electronically conductive minerals with spectral induced polarization (SIP) has been a longstanding goal of mineral exploration. A more comprehensive understanding of how mineral texture (e.g. disseminated grain versus veinlet forms) and mineral type influence SIP signals is required, as the polarization mechanisms at the electron conductor-fluid interface and their effect on SIP signatures remain incompletely understood. We performed laboratory measurements on synthetic copper (Cu) and graphite (Gr) samples prepared in disseminated grain and veinlet forms by mixing fixed concentrations of the electron conductor with a non-polarizable sand. Current density dependent SIP (complex impedance) responses were assessed by measuring at multiple current densities, and non-equilibrium was investigated by making measurements as a function of time. SIP responses that remained independent of current density were interpreted as evidence of non-Faradaic processes, with possible minor Faradaic contributions, being dominant. Conversely, current density dependent complex impedance responses were interpreted as an increased contribution from Faradaic polarization processes involving redox reactions at the electrolyte-electron conductor interface. Veinlet Cu samples exhibited strong current density dependent SIP responses and prolonged equilibration, in contrast to the negligible current density dependence of veinlet Gr and disseminated grain Cu and Gr samples. All samples that exhibited negligible current density dependence could be fit to a Pelton-Cole model, whereas samples exhibiting strong current density dependence could only be fit with a more flexible Debye decomposition model. Gr samples displayed a fluid-resistivity dependent peak in the phase response at a specific frequency consistent with mechanistic models, while Cu samples deviated from predicted behaviour. Measurements on rock cores were consistent with observations on synthetic samples, with the SIP response of a Cu bearing rock core exhibiting a clear dependence on current density, whereas the Gr bearing rock core showed negligible current dependence. The observed dependence of the phase and complex resistivity on current density, mineral type, and texture is consistent with differences in the relative importance of Faradaic and non-Faradaic polarization processes. These results suggest that distinguishing current-density-dependent Faradaic responses from current-independent non-Faradaic responses may improve the discrimination of mineralogical and textural characteristics using SIP measurements, although field-scale applications will be limited due to the lower current densities relative to those used in this laboratory study.

SummaryFull Waveform Inversion (FWI) leverages the amplitude and phase of full wavefield seismic data to obtain high-resolution subsurface structures by minimizing an objective function. However, FWI encounters severe cycle-skipping issues when the seismic data lack low-frequency components or the initial velocity model is significantly different from the true velocity model. Furthermore, an incorrect source wavelet can also directly affect the accuracy of inversion results. To overcome these challenges, we introduce the waveform-preserving source-independent wave-equation-based local-scale traveltime inversion (WPS-LTI) method. This approach aims to recover the low-wavenumber velocity structures as an initial velocity model, while reducing the influence of an incorrect source wavelet. In the WPS-LTI method, seismic data are processed via both convolution and deconvolution with reference traces, resulting in modified observed and synthetic seismic data that share the same seismic wavelet. This approach helps avoid waveform distortions typically introduced by conventional source-independent methods that rely solely on convolution. By preserving the waveform characteristics, we can use cross-correlation to compute more accurate local-scale traveltime differences, thus recovering accurate low-wavenumber velocity structures. Tests conducted on the Marmousi model and the Chevron 2014 blind test dataset demonstrate that the WPS-LTI method can significantly alleviate cycle-skipping issues and reduce the dependence on the source wavelet in FWI.

SummaryWith the shale gas exploration and development into deep strata (> 3500 m in depth), weak signal detection becomes more challenging and vital in surface-based microseismic monitoring. Most existing deep learning approaches are trained on synthetic or single-station data, which limits their performance in detecting weak events in field seismic recordings. Here, we propose a deep learning approach based on the computer vision object detection model YOLOv5 via transfer learning and train it on a dataset from field observation to better detect microseismic events in multi-channel seismic data. First, we convert continuous seismic data into waveform images in grayscale with a record length of 20 s, and manually annotate the microseismic events containing the first arrival signals using rectangular boxes. To train the model, the transfer learning strategy is used by using parameters of YOLOs in computer vision to initiate the model. After training and validation, we conducted comprehensive tests to examine the performance of the model. Compared with the single-station short-term average/long-term average (STA/LTA) detector, our method reduces false positives markedly, yet it also delivers higher recall rate than the available multi-station detection scheme. When applied to continuous dataset from a different area, the model achieved good performance with a precision of 0.9832 and a recall of 0.9702, indicating it learns representative features of microseismic signals. All these tests illustrate that the proposed method achieves high detection accuracy, low false positive rates, and robustness. Meanwhile, the method possesses the superiority for real-time monitoring and can be easily updated by adding a few new samples to the training dataset.

SummaryDeep learning methods such as neural networks offer a rapid alternative to geophysical inversion and are increasingly applied to airborne electromagnetic (AEM) data for a range of earth science problems. Such networks rely on both the network architecture as well as suitable training datasets to produce realistic earth resistivity models. Here we develop a Convolutional Neural Network (CNN) to invert AEM data and pay particular attention to the characteristics of the training dataset. As well as covering a realistic range of subsurface property distributions, the training data must emulate the noise and data gap characteristics of the target survey data. We find neural networks trained on noise free data, while performing very well on validation data, perform poorly on real world data, often inducing geologically unrealistic artifacts in the resistivity model. Conversely, training data conditioned with noise and data gap characteristics appropriate for the survey of interest, fit the validation data less well, but produce fewer artifacts when applied to survey data. Choice of loss function metric and data normalisation method is important and may vary between surveys. Despite strong advances in neural network design, developing training data with broad enough characteristics to be generalised across many surveys and data acquisition platforms remains an important challenge.

SummaryAccurate separation of the Earth’s internal and external geomagnetic fields requires robust modeling, especially during the geomagnetic storms. Because the comprehensive CHAOS model is primarily designed to represent the Earth’s magnetic field under geomagnetically quiet periods, residuals of satellite observations relative to the CHAOS model during storms exhibit systematic dependencies on magnetic local time (MLT) and the disturbance levels, highlighting limitations in its representation of external field morphology. To address this limitation, we developed an empirical correction for the CHAOS model that specifically targets the magnetic signatures of the residual ring current at low and mid-latitudes. Using the high-quality scalar magnetic data from Swarm A/B and the Macau Science Satellite-1 (MSS-1) at nighttime hours, between November 2023 and September 2025, we modeled the systematic perturbations remaining after CHAOS background removal. Our separable parameterization explicitly resolves MLT structure while incorporating dependencies on geomagnetic disturbance levels. Event analyses demonstrate that our correction model can successfully capture the pronounced dawn-dusk MLT asymmetries of the external field during storms main and recovery phases. Statistically, the empirical model robustly captures residual ring current signals across all evaluated nighttime sectors, yielding the greatest improvements at dusk. Ultimately, the proposed correction effectively mitigates the activity-dependent biases in the CHAOS model’s external field, yielding a more reliable representation for the accurate interpretation of storm-time magnetic residuals.

SummaryGeophysical joint inversion offers a robust approach for integrating diverse geophysical datasets, leveraging their complementary information to mitigate the nonuniqueness inherent in individual data inversion. This approach thereby enhances model consistency and improves subsurface characterization. In this study, we investigate the effectiveness of three-dimensional (3D) joint inversion of gravity, magnetic and magnetotelluric (MT) data for deep geothermal exploration. We first implemented a reweighting optimization strategy for gravity and magnetic data to increase the vertical resolution of jointly inverted models of potential data. Subsequently, we integrated the optimized gravity and magnetic inversion scheme with MT data using a cross-gradient constraint for comprehensive joint inversion. Synthetic experiments validated the effectiveness of our developed reweighting optimization strategy and joint inversion scheme, demonstrating improved structural resolution and model reliability compared to separate inversions. Finally, we applied the methodology to the Yanggao geothermal field in northern Datong Basin, China, to assess its applicability for deep geothermal exploration. The joint inversion approach, offering an advantage over separate inversions, provides enhanced subsurface characterization by more accurately delineating hydrothermal pathways and potential geothermal reservoirs of the geothermal system within the investigated area.

SummaryLocal seismic topographic amplification of rock slopes is a primary cause of earthquake-induced damage in mountainous engineering. To resolve the scattering of plane SH waves by individual rock slopes, this study introduces an analytical model derived from the Uniform Theory of Diffraction. Conventional wave function expansion methods (WFEM) often struggle with non-closed free boundary conditions across varying elevations and tend to suffer from series divergence under high-frequency incidence. This approach precisely partitions the slope into two semi-infinite ideal wedge domains sharing a common face, reconstructing the generalized geometric wave field through rigorous global ray tracing. Second-order coupled diffraction terms with a slope-diffraction correction are introduced to describe the two-way interaction between the slope crest and toe, while the UTD transition functions ensure continuity of the wavefield across shadow boundaries. Systematic numerical evaluations in both frequency and time domains reveal a frequency-selective interference mechanism with distinct spatial variability. Under horizontal grazing incidence, energy is strongly focused on the leeward face, while the back-slope forms a high-frequency seismic shadow zone due to intense destructive interference between direct and secondary diffracted waves. At certain oblique incidence angles, repeated scattering of broadband waves between the slope crest and toe leads to strong local resonance in the slope region. Findings indicate that the location of maximum dynamic response shifts sensitively with the incident wavefront inclination. The results show that the maximum amplification does not always occur at the slope crest, highlighting the importance of incidence angle in seismic microzonation and site planning for mountainous areas.

SummaryThe statistics of earthquake populations are governed by stress changes and event interactions at various scales which promote non-linear event cascades and seismicity clustering. Earthquake clustering in space and time is observed in both tectonic and volcanically dominated regions; however, notable differences exist for instance during eruptive sequences. Here, we characterize seismicity statistics in California and Hawaii and quantify differences in background rates and clustering. We compare three statistical approaches, i.e, i) Reasenberg declustering, ii) Gamma distribution fit, and iii) Nearest Neighbor Clustering which are applied to seismicity in Hawaii, southern California, and ETAS catalogs with known clustering characteristics. Tests with ETAS-catalogs suggest that background rate variations are more consistently resolvable using the Nearest Neighbor and Gamma distribution algorithm without requiring parameter adjustments to the specific study area. The Nearest Neighbor method tends to overestimate background rates, whereas the Gamma distribution algorithm shows systematic underestimation, suggesting that an ensemble model can improve the performance over any single method. Uncertainties in background rate estimates based on small magnitude events (e.g., M≥2.5) are high, in particular for the Reasenberg method, which overpredicts true rates by up to a factor of 4, and extrapolations to larger magnitudes need to be treated cautiously. Our analysis shows that the combination of the three approaches generates new insights into the processes that govern seismicity clustering. The analysis of southern California seismicity produces consistent results between the different methods, and the distribution of background seismicity rates is indistinguishable from Poissonian rates at magnitudes above 3.5. Conversely, the seismicity in Hawaii exhibits notable non-Poissonian background rate changes. Results from the Nearest Neighbor clustering analysis suggest notable differences between tectonic and volcano seismicity, including background rate changes and more pronounced spatiotemporal clustering around active volcanoes. Hawaiian seismicity is characterized by more localized spatial-temporal clustering, in particular during the 2018 eruption. This difference may be a result of magmatic, hydrothermal and eruptive processes, which make volcanic areas susceptible to stress perturbations and earthquake triggering at small scales.