Geophysical Journal International

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.

SummaryElectromagnetic induction logging-while-drilling tools measure real-time electromagnetic data to characterize the electrical properties of formations, thereby supporting geosteering decisions. The multiscale structures formed by heterogeneous fractures within complex reservoirs lead to electromagnetic field discontinuities, which in turn limit the accuracy of reservoir evaluation. To elucidate the mechanisms by which multiscale fracture conductivity and spatial distribution influence the attenuation of the electromagnetic field, a discontinuous Galerkin method based on vector basis functions is developed, using the open-source finite element library deal.II, for modelling complex fracture clusters. To reduce the number of elements associated with the explicit volumetric discretization of fractures and truncate the computational domain, the impedance transition boundary condition that treats the fractures as element interfaces and the perfectly matched layer approach are employed. A goal-oriented adaptive mesh refinement strategy based on element-wise residuals and interelement field jumps is adopted to flexibly refine the mesh near the transmitter–receiver system and fractures, since the complex distribution of field discontinuities induced by fracture clusters renders empirical mesh refinement inadequate. The numerical results validate the accuracy and performance of the algorithm. For a layered fracture model, the numerical accuracy is assessed by comparison with the semi-analytical solutions from the open-source software empymod. In interlaced fracture models with different in-plane dimensions, the flexibility and adaptability in modelling fractures are demonstrated by the mesh distributions after adaptive refinement. For fracture cluster models, an advantage in computational efficiency over the finite element method for complex fracture networks is demonstrated as fracture complexity increases. These results demonstrate the applicability to electromagnetic induction logging-while-drilling modelling of fractures in conductive-to-moderately resistive sedimentary formations, with improved computational efficiency and modelling flexibility.

SummaryGlobal Navigation Satellite System (GNSS) technology and Gravity Recovery and Climate Experiment (GRACE) satellite gravimetry provide essential observational constraints for monitoring terrestrial water storage (TWS) changes. In this study, we develop a spatiotemporally constrained joint inversion model that inverts GNSS and GRACE observations to estimate reliable TWS changes over Brazil from January 2008 to July 2016. We evaluate the performance of the spatiotemporally constrained joint inversion model using both closed-loop simulation and independent GNSS surface displacement observations. The simulation results indicate that the joint inversion results using spatiotemporal constraints outperform GNSS-only solutions, GRACE-only solutions, yielding improved accuracy and reliability than joint inversion results with spatial constraints. The corresponding standard deviations are 45.94 mm, 57.44 mm, 52.53 mm, and 48.88 mm, respectively. Analysis from the measured data shows that vertical displacement time series simulated from TWS changes derived from spatiotemporally constrained joint inversion indicate higher consistency with GNSS-observed surface displacement time series at nine GNSS stations compared with the GNSS-only, GRACE-only (from CSR-M), and the joint inversion results with spatial constraints. The corresponding standard deviations and correlation coefficients are 1.93 mm & 0.92, 2.54 mm & 0.87, 4.45 mm & 0.74, and 1.94 mm & 0.91, respectively. Meanwhile, the joint inversion results from measured GNSS and GRACE data using spatiotemporal constraints show lower uncertainty and higher stability than those of GNSS-only and joint inversion results with spatial constraints. The proposed inversion model provides an alternative means to fully exploit the potentials of GNSS and GRACE technologies for integrated monitoring of TWS changes.

SummaryReanalysis of Laterolog-7 (LL7) data in thin-bed sequences suffering from severe resistivity suppression caused by shoulder-bed averaging and mud filtrate invasion remains important. We developed a forward modelling and inversion workflow based on the Finite Difference Method (FDM) in cylindrical coordinates (r, z), reducing the problem to an efficient two-dimensional solver. The LL7 bucking current focusing condition is satisfied at every forward step, ensuring correct reproduction of the dynamic instrument response under high resistivity contrasts. The bucking current ratio log (n1/n2) serves as a boundary-detection indicator to initialize the layered-earth model and reduce non-uniqueness. Formation resistivity profiles are recovered through an Inexact Gauss-Newton (IGN) scheme that avoids explicit Jacobian formation, maintaining computational efficiency while guaranteeing a valid descent direction. Synthetic validation demonstrates recovery of the resistivity structure of a single-bed model (0.6 m bed of 30 Ωm with 5 Ωm invaded zone). The inversion successfully recovers the true resistivity of 30 Ωm and background of 10 Ωm. A multi-bed test shows the inversion resolves bed boundaries at decimeter scale, recovering layers as thin as 20 cm from a smeared raw log. Application of the inversion to a real LL7 dataset from the Cuu Long Basin, Vietnam, confirms field applicability. The inversion reconstructs sharper resistive peaks throughout the thin-bed interval, providing an improved representation of true formation resistivity.

SummaryPetrophysical characterization in near-surface environments remains challenging due to the high heterogeneity, sharp contrasts in fluid content, and complex hydrogeological structures. This work introduces a novel petrophysical joint inversion (PJI) framework that integrates SH-wave seismic refraction tomography (SRT) with electrical resistivity tomography (ERT) and, optionally, time-domain induced polarization (IP) data, replacing the conventional approach based on P-wave velocity and Wyllie’s equation with the Castagna relation for S-wave velocity. By reducing sensitivity to pore-fluid effects, the framework preserves robust saturation estimates while significantly improving porosity reconstruction in both its spatial distribution and absolute values. A synthetic example and two field case studies with contrasting hydrogeological layering demonstrate the capability of our method to provide physically consistent, fully quantitative, and high-resolution characterization of complex near-surface scenarios, extending the scope of PJI beyond classical applications. Additionally, in settings where the IP contribution is limited, the strong sensitivity of S-waves to site lithology allows for a comprehensive characterization without including IP data, which are required in traditional approaches to fully characterize complex scenarios, thereby reducing survey costs and overall complexity.

SumamryThe introduction of full waveform inversion has been the starting point of four decades of investigation in seismic imaging, which aim to reduce the ill-posedness of this nonlinear and quantitative imaging method. In this review, we propose a comprehensive overview of these investigations. Beyond simply cataloging proposed approaches, this work aims to classify them, identify similarities and differences between methods, establish connections with common concepts and practices in full waveform inversion, and emphasize the few core ideas underlying this long-term research effort. We introduce a general mathematical framework that encompasses all proposed methods, which we use throughout this study to support their presentation, analysis, and comparison. The concluding discussion should promote a better understanding of the current state of research on this topic. In particular, we examine which methods have been applied to field data and which have not, highlighting the main difficulties and the requirements such techniques must fulfill to transition to such applications.

SumamryElectrical resistivity models derived from magnetotelluric (MT) measurements across the Khangai Dome, an intracontinental plateau in central Mongolia, have identified several intriguing features, including a locally thinned lithosphere and extended fluid-rich domains in the lower crust. However, due to the limited spatial coverage of these studies, whether those features extend beyond central Mongolia could not be assessed. This is an important question with major scientific implications for understanding the geodynamic evolution of the region. The surrounding areas also host significant mineral deposits and geothermal resources, whose origin and interpretation depend on the characterization of deep structures to be fully understood. In this study, we report on MT data from 378 new locations to the west, east, and north of the Khangai Dome, acquired between 2020 and 2023. A new 3-D electrical resistivity model of the region was derived using a combination of the new and previously acquired datasets. This new model covers an area of approximately 900 × 1250 km2, substantially expanding the spatial coverage compared to previous models. The key methodological novelty of this study is the implementation of a recently developed, scalable, and open-source 3-D inverse solver based on the integral equation approach and non-local parametrization. The 3-D model fits the observed data remarkably well, and the recovered features show good agreement with those from previously published models obtained with different solvers and subsets of the data. The model reveals an upper crustal dichotomy, with very high resistivity in the northern region and low resistivity in the southern region. The boundary between these regions follows the Main Mongolian Lineament in the central and eastern regions and the Ikh-Mongol arc system in the western region. In the lower crust, the resistivity is generally lower, and several long, laterally extended, very low-resistivity anomalies are observed. These anomalies extend westward from central Mongolia but are bounded to the east by a high-resistivity anomaly whose location coincides with the Mogod fault zone. The aforementioned features highlight the control of major tectonic boundaries on the electrical resistivity distribution. An upper-mantle low-resistivity anomaly below the Khangai Dome, previously interpreted as a locally thinned lithosphere with an upwelling asthenosphere, is imaged not only beneath central Mongolia but also, to some extent, in the western region. The model also provides new information on the structure of major features across central-western Mongolia, such as the extensive Bulnay fault. Taken together, these new constraints on lithospheric properties and architecture advance our understanding of the mechanisms that shaped the region and its subsequent evolution.

SummaryReciprocal measurements are standard for data quality control in electrical resistivity surveys, but their potential as an interpretation tool is rarely exploited. Under the assumption of linear subsurface behaviour, normal and reciprocal readings should yield the same results, but when a systematic offset is present it can be indicative for additional processes. Here we present a time-domain induced polarization (TDIP) dataset from a volcanic hydrothermal system (Reykjanes, Iceland) where a strong, systematic offset between normal and reciprocal measurements is observed. The discrepancy is present as a positive shift of the reciprocal decay curves (>100 mV/V). While both normal and reciprocal inversions resolve a strong IP anomaly in the southern part of the profile, a substantial offset exists between the two ( >20 mS/m). The offset is spatially confined and persistent over 100 consecutive days, the anomaly itself coincides with a strong IP response attributed to the mutual occurrence of clay minerals and disseminated iron sulphides and oxides. Besides random noise, we evaluate four possible mechanisms for the offset: (1) polarization of the electrodes, (2) inherent sensitivity differences between normal and reciprocal configurations, (3) stray currents from the nearby powerplant and (4) nonlinear IP effects where the response depends on the current density. Given the geological context and the localized nature of anomaly, a non-linear effect is plausible but the ambiguity between systematic noise and a genuine subsurface signal cannot be fully resolved with field data alone. This study demonstrates that the systematic analysis of normal-reciprocal misfit in TDIP can serve as a widely accessible tool for identifying signals that, whether geologic, or anthropogenic in origin, are overlooked by standard data quality workflows.

AbstractTomographic reconstruction involves fundamental trade-offs between resolution, data coverage, and parameterisation. Conventional approaches commonly rely on regularised pixel-wise or cell-based parameterisations. These regularisations are valuable and can themselves be geology-informed when their assumptions match the target setting. However, they often encode structural information only indirectly through generic penalty terms and regularisation weights. As a result, they can be less suited to questions in which the objective is to test a specific structural hypothesis. We develop a Bayesian framework for tomographic reconstruction and uncertainty quantification that formulates such hypotheses explicitly and assesses them through posterior uncertainty. As an illustrative application of the framework, we consider ambient noise tomography. The phase velocity reconstruction is reformulated to address a targeted geological hypothesis by introducing a parameterised structural prior. The prior is constructed from a Whittle–Matérn latent Gaussian field, combined with spectral dimensional reduction via a truncated Karhunen–Loève expansion and a differentiable nonlinear pushforward mapping that promotes approximately piecewise-constant velocity structures. The prior parameters, including the correlation length, smoothness, pushforward slope, and velocity bounds, provide interpretable hyperparameters related to expected lateral scale, structural regularity, transition sharpness, and plausible phase-velocity contrasts. In this formulation, the prior is used as a parameterised structural assumption informed by geological context, rather than as hard geological knowledge. The resulting parameterisation substantially reduces dimensionality, making both maximum a posteriori estimation and Hamiltonian Monte Carlo sampling computationally feasible for the synthetic and observed datasets considered here. Synthetic experiments demonstrate accurate recovery of prescribed velocity contrasts and show that posterior uncertainty patterns reflect acquisition geometry and data coverage. Application to a ∼1000-sensor dense array dataset produces structures consistent with the imposed geology-informed prior assumptions, while providing spatially resolved uncertainty estimates that identify regions of robust inference and remaining ambiguity. The framework offers a proof-of-concept formulation for goal-oriented, uncertainty-aware seismic tomography when prior information supports an explicit structural hypothesis that can be encoded through the proposed Whittle–Matérn latent-field parameterisation.

SummarySubterranean rock masses are in true triaxial stress (TTS) fields, which can preferentially close internal microcracks, thereby inducing elastic anisotropy. However, the seismic response under TTS is poorly understood. Here, an analytical PP-wave reflection coefficient for triaxial stress-induced anisotropic media modelled by microcrack closure mechanism is proposed for seismic inversion. Considering a micromechanical model in which microcracks are represented by stress-dependent compliances, the stress magnitude and orientation are firstly incorporated to account for the elastic anisotropy resulting from triaxial stress-induced microcrack closure. A good agreement is obtained between the model predictions and the existing laboratory measurements. Based on weak anisotropy assumption, we then deduce the effective stiffness tensor of the triaxial stress-induced anisotropic media. Three sets of stress-related anisotropy indicators (SRAIs) are introduced to quantify the microcrack closure effect and anisotropy magnitude induced by triaxial stress. Furthermore, a linearized PP-wave reflection coefficient equation for triaxially stressed isotropic media is derived using the scattering theory. Numerical results validate the feasibility and accuracy of the proposed formula, and reveal the influence of TTS on the PP-wave amplitude variation with angle and azimuth (AVAz). Finally, due to the highly ill-conditioned AVAz inverse problem in triaxial stress-induced anisotropic media, a model and data driven inversion approach is proposed through building on a convolutional neural network. We stepwise estimate the isotropic elastic parameters and SRAIs based on azimuthal seismic amplitude difference inversion strategy. Tests on both synthetic and real seismic data indicate that the nonlinear AVAz inversion framework outperforms conventional approaches in terms of stability and accuracy. Our study allows the construction of elastic properties in triaxial stress-induced anisotropic media, and may provide new insights into determining in-situ stress from seismic data.

SummaryBy leveraging a significantly expanded GNSS dataset (997 vectors, 335 new), we construct a high-resolution crustal strain rate model for the Southeastern Tibetan Plateau (SETP) that approximately doubles the resolution of previous studies. Our model delineates a prominent J-shaped high-strain zone that follows the Xianshuihe-Xiaojiang fault system (XXFS) and extends ∼600 km southwestward. This feature, combined with unsupervised Euler pole clustering (EPC) machine learning results, provides further evidence for a mid-to-lower crustal linkage between the XXFS and faults southwest of the Red River fault (RRF), suggesting the active XXFS is propagating across this ancient boundary (RRF) at depth. Within the SETP interior, our model reveals widespread extensional deformation characterized by a systematic clockwise rotation of its principal axes from north to south, alongside discrete rotational domains partitioned by major faults. Shear strain patterns within the SETP indicate that major strike-slip faults, predominantly sinistral, drive regional strain localization. Moreover, the close association between modeled high-strain zones and historical M ≥ 6 earthquakes suggests that our high-resolution model can serve as a refined, physically-based perspective for regional seismic hazard assessment.

SummaryWhile Bayesian inference is widely applied to magnetotelluric, seismic, and magnetic datasets, its application to gravity anomalies remains relatively unexplored. To address this, we introduce InDIA (Inversion of Density Interface and its Application), a Python-based Bayesian inference framework utilizing a structural Gaussian process. InDIA is developed to estimate discretized subsurface layer depths and densities while explicitly quantifying uncertainty. Unlike previous techniques restricted to continuous depth profiles and localized issues, InDIA integrates diverse prior information to effectively resolve both local and regional gravity anomalies. This approach successfully mitigates the multi-parameter challenges and convergence issues common in local gradient-based optimization techniques, such as Adam. The superiority of this algorithm has been validated through various synthetic models (involving multi-prism configurations, two-layer models with lateral density variations, heterogenous subsurface model with lateral and vertical density variation along with faulted dipping models) featuring both constant and variable density distributions—whether lateral, vertical (prior-constrained), or both—incorporating Gaussian noise, Random-walk noise, Systematic noise and Salt and Pepper noise to replicate real conditions. Furthermore, the framework’s field applicability is demonstrated using two real gravity datasets. First, it tackles multi-layered subsurface profiling in the West Korea Basin. Second, it simultaneously resolves depth and prior-constrained vertical density variations in the Godavari Basin, India. In both field applications, the inverted parameters are geologically viable and closely align with previously established models, confirming InDIA as a reliable tool for gravity data inversion.

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.