Geophysical Journal International

SummaryThe estimation of topographic gravity field models has attracted significant interest in recent years due to its growing relevance in Earth sciences. In this study, we present a robust methodology for the computation and comprehensive validation of global, complete spherical Bouguer and isostatic gravity anomalies that are essential for accurately interpreting subsurface mass distributions therefore geological structures. We synthesize these crucial gravitational functionals by leveraging spherical harmonic coefficients from high-resolution global gravity field models and various topographic/topographic-isostatic gravity field models. Our findings underscore the critical role of comprehensive terrain corrections in deriving physically meaningful, complete Bouguer gravity fields. The calculated global anomalies demonstrate strong coherence with established benchmark datasets, such as the World Gravity Map 2012. Residual differences are primarily attributed to variations in input Digital Terrain Models. Comparisons with regional Bouguer datasets reveal systematic biases that are largely explained by differing terrain correction methodologies. After removing this effect, there is a high level of consistency between the calculated global and published regional datasets, highlighting the utility of our global solutions, particularly in regions with sparse terrestrial data. Furthermore, the globally computed isostatic gravity anomalies exhibit significant agreement with both external global and diverse regional datasets, notably without the large systematic biases observed in Bouguer comparisons. This agreement reflects the effectiveness of the combined topographic and isostatic corrections in capturing Earth’s mass balance. This research provides valuable tools for new studies in the geoscience community by offering globally consistent and complete Bouguer and isostatic gravity field anomalies that have been rigorously validated for the ICGEM service.

SummaryDestabilization of volcanic edifices can generate debris avalanches with catastrophic impacts on their environment. We present the first high-resolution muography of Mount Unzen, Japan, conducted to characterize the structure of lava lobes formed on the volcano’s summit and flank during the 1990-1995 eruption. A multi-wire-proportional-chamber-based muon tracking system was operated for 203 days. The obtained high-resolution muographic image shows the internal density structure of Mount Unzen with a spatial resolution of 12 meters. Mean densities were respectively measured as 2,470 kg m−3 and 2,290 kg m−3 for the base rock and a fracture zone, and both were consistent with the results of prior drilling and sampling experiments. The mean density of lava lobes was measured significantly lower value of 1,570 kg m−3, indicating post-eruptive structural weakening. A comparison between the time-series of muographically measured density-lengths and daily precipitation records suggest that rainfall-induced gravitational destabilization did not occur during the observational period. This work demonstrates that long-term (multi-year) muon monitoring of the lava lobes can provide valuable complementary information for volcanic stability assessments.

SummaryThe asthenosphere is a weak layer in the upper mantle that supports the movement of the overriding tectonic plates and facilitates mantle convection. In this study, we compile a global dataset of SS precursors reflected at the base of the asthenosphere, also known as the 220-km discontinuity. The global dataset includes the oceanic SS precursors from Sun & Zhou (2025a) and new measurements with bounce points in continental regions. Similar to the oceanic dataset, the continental SS precursors are observed on about 45% of the SS waves, with bounce points distributed across all tectonic regions — from orogeny belts to stable cratons. We image the depth of the discontinuity at a global scale using finite-frequency tomography. In oceanic regions, the depth of the 220-km discontinuity agree well with the previous study, with discontinuity depth structure characterized by alternating linear bands of shallow and deep anomalies that roughly follow seafloor age contours. In continental regions, the variations are not spatially oscillatory but are instead much broader, with prominent perturbations associated with convergent plate boundaries. The base of the asthenosphere is shallow along the southern border of the Eurasian plate, from the Mediterranean region to Southeast Asia. Shallow discontinuity anomalies are also observed in the continental interiors – in Eurasia, from the northern Tian Shan through Mongolia to eastern Siberia, and in North America east of the Rocky Mountains. These anomalies form a linear structure roughly parallel to the Pacific subduction zones. The average depth of the discontinuity, as well as the velocity contrast across the interface, is globally consistent across both oceans and continents, with an average depth of approximately 251 km and a velocity increase of about 7%. Given that the continental lithosphere has been cooling for much longer than the oceanic lithosphere, the observed consistency in the average depth of the discontinuity implies that secular cooling does not significantly impact the thermal structure at the base of the asthenosphere.

SummaryA very frequent approach for studying lithospheric processes is to deploy temporary seismological networks in dedicated areas and to map the mantle structures with different approaches. One of them is the well-established relative travel time body wave tomography. Different circumstances often lead to a non-uniform deployment of stations both in space and time, and a wish to combine data which have been acquired asynchronously. This is the situation in Patagonia where two distinct seismic experiments provide complementary seismic data over the region covering the Patagonia slab window. Combining these data in one regional relative body wave tomography is however problematic as the two data sets are a priori with respect to two different reference models. In this contribution, we show that the number of finite-frequency relative travel time residuals varies very strongly from station to station for this data set, violating the assumption implicit in relative travel time tomography of a unique reference model due to an even data distribution for all events. We demonstrate the superiority of the inversion using relative sensitivity kernels compared with a traditional approach with absolute kernels and event terms. A resolution test proves how this is crucial for resolving the important issue of the eastern extent of the slab window. In addition, we discuss potential issues related to interference of the direct phases with core phases when measuring finite-frequency travel time residuals by cross-correlation of waveforms in necessarily relatively large time windows. We also briefly outline our preferred strategy for performing crustal correction, keeping in mind that finite-frequency residuals require frequency-dependent crustal corrections.

SummaryP-wave receiver functions (RFs), which utilize P-wave conversions to probe subsurface structures, face significant challenges in sedimentary environments. Specifically, strong reverberations generated by ultra-low-velocity sedimentary layers distort RF waveforms and obscure crustal signals, posing challenges for robust shallow crustal imaging. We develop a novel Bayesian joint inversion framework that simultaneously utilizes three complementary datasets—reverberant receiver functions, dereverberated receiver functions, and surface wave dispersion—to address this challenge. Our approach employs Unscented Kalman Inversion, a derivative-free method that efficiently handles nonlinear joint inversion problems. Synthetic tests demonstrate that our joint inversion recovers sediment thickness and Moho depth with uncertainties of ±0.50 km and ±1.0 km, respectively. Application to real data from the Songliao Basin verifies the approach, successfully reconstructing sediment thickness and Moho depth beneath sedimentary cover. This methodology demonstrates potential for advancing crustal investigations in complex sedimentary settings, such as continental rift basins and oceanic margins, where sedimentary sequences of variable thickness often obscure deeper structures.

SummaryThree-dimensional (3-D) forward modeling of magnetotelluric (MT) data remains a computationally challenging task, particularly when accurate broadband MT responses are simulated for real-world problems that often involve complex multi-scale bathymetry and/or topography. To overcome this challenge, we developed a new efficient numerical approach for 3-D MT forward modeling that combines high-order Nédélec-type finite elements and high-order meshes, allowing us to obtain superior accuracy and account for complex material boundaries and interfaces. Despite gains in accuracy, higher-order FE solvers are often considered impractical owing to higher memory consumption and a more ill-conditioned system. To overcome these limitations, we use an iterative solver accelerated by the Low-Order-Refined (LOR) preconditioner, which uses spectrally equivalent low-order operators, rendering the complexity independent of the polynomial degree. Another key novelty is a matrix-free implementation, where the action of the high-order operator is computed efficiently without explicit matrix assembly. The low-order system is solved using an Auxiliary Space Maxwell (AMS) solver based on a multigrid solver. We demonstrate the efficiency in a series of numerical experiments. Scalability analysis on a 3-D benchmark model demonstrates that the LOR preconditioner significantly outperforms the current state-of-the-art AMS preconditioner in terms of CPU time and memory usage, especially for higher polynomial degrees. Excellent scalability is confirmed by solving a problem with up to 1.5 × 109 degrees of freedom in less than 2 minutes using 16,384 CPU cores, which is, to our knowledge, the largest 3-D MT problem reported to date. We also illustrate that high-order hexahedral meshes allow for accurate discretization of complex geometries, such as topography, with substantially fewer elements than conventional linear meshes. Finally, the capability of the integrated approach is demonstrated on a real 3-D model crossing the ocean trench in the Aleutian subduction zone. The proposed methods pave the way for more efficient and accurate 3-D MT modeling that is crucial for the inversion of complex MT data sets.

SummaryAtmospheric models are based on various types of geophysical data, including lidar and radar. Infrasounds, acoustic waves that can propagate over large distances, have not yet been used in atmospheric models, although they provide valuable information. Besides their sensitivity to atmospheric phenomena such as gravity waves, infrasound also presents the advantage of being omnipresent. Previous studies explored the use of infrasound packet arrival properties for model estimation. However, properties such as arrival times present less information than full waveforms. We aim here to investigate, for the first time, the sensitivity of a full infrasound waveform to model parameters and to use these sensitivities in an inverse problem to recover atmospheric structure. For this purpose, infrasound propagation is modeled by Euler equations (i.e. Navier-Stokes equations in the absence of attenuation effects), and discretization is carried out here using the finite-differences method. Waveform sensitivity to atmospheric parameters is computed through the adjoint method via a novel and optimized double checkpointing-based procedure and validated by comparison with a small perturbation method. As an illustration, these sensitivity kernels are computed for the idealized case of an explosion in Finland, recorded by a CTBT station. These first results demonstrate the high sensitivity of infrasound waveforms to the atmospheric perturbations generated by gravity waves. Moreover, the sensitivity kernels of infrasound waveforms allow us to recover the variations of model parameters by solving an inverse problem. To demonstrate this capability, full waveform non-linear inversions are performed using the Limited Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS): wind and sound speed profiles are inverted for a test case with idealized conditions and a synthetic dataset. These estimates of infrasound sensitivity kernels are closing a knowledge gap that allows the use of infrasound full waveforms to constrain atmospheric models.

SummaryPortland cement remains the most widely used construction material globally, valued for its well-documented properties and performance. However, its production generates substantial CO₂ emissions, mainly due to the decomposition of limestone (CaCO₃) into calcium oxide during clinker formation. In response to these environmental concerns, researchers have been actively exploring ways to lower cement’s carbon footprint and improve its sustainability. One effective strategy involves reducing the clinker content by incorporating supplementary cementitious materials (SCMs). To ensure SCMs enhance performance without compromising safety, it is essential to investigate the properties of blended cements. Natural zeolites have emerged as promising SCMs. Although they do not possess inherent cementitious properties, finely ground zeolites can react with calcium hydroxide in the presence of water, contributing to strength development. This study examines the potential of natural zeolites as SCMs and utilizes the spectral induced polarization (SIP) method to monitor cement hydration and reaction mechanisms. Portland cement mortars containing 25% zeolite were prepared and compared against two reference mixes. Zeolites were analyzed using scanning electron microscopy (SEM) and X-ray diffraction (XRD), while SIP monitoring was conducted continuously over 28 days. Our results reveal that SIP responses are influenced by the specific chemical composition of the mortar. The incorporation of SCMs alters cement chemistry, significantly influencing SIP signals. Over time, we observed an increase in the imaginary conductivity component and a decrease in the real conductivity component. SEM analysis showed the formation of new fibrous mineral habits in zeolite-blended samples, alongside a reduction in pore fluid content. These observations suggest a strong connection between SIP signals and mineralization processes, likely associated with the formation of secondary gels and calcium monosulfoaluminate. The interaction of zeolites with calcium hydroxide promotes the development of calcium aluminate hydrates, which then react with ettringite to form calcium monosulfoaluminate. These results emphasize the importance of studying SIP behavior in cement systems containing SCMs, as assumptions based on ordinary Portland cement may lead to misinterpretations. Our research underscores the potential of SIP as a valuable tool for monitoring cement hydration while offering new insights into the chemical transformations in zeolite-containing mortars. Ultimately, this work contributes to the advancement of more sustainable cement formulations, supporting environmentally responsible construction practices.

SummaryElastic rock physics models are widely used to estimate the saturation of hydrate in isotropic sediments. However, for isotropic media, the influence of heterogeneously distributed hydrate on the P- and S-wave velocities remains unclear, leading to uncertainties in hydrate saturation estimates. To address this issue, in this work we proposed a double-solid-matrix model for predicting the velocities of sediments hosting heterogeneously distributed hydrates. A comparison of simulated velocities of our model and two rock physics schemes designed for homogeneous distributed hydrate (i.e. matrix-supporting and pore-floating models) show that, our model predicts higher S-wave velocity than matrix-supporting and pore-floating models, but yields similar P-wave velocity estimates as matrix-supporting model. We apply our model to two marine hydrate sites in the Cascadia margin: Site 1245 from Ocean Drilling Program Leg 204 and Site U1328 from International Ocean Drilling Program Expedition 311. Two locations yield similar results: velocity estimates from our model are much closer to downhole measurements than matrix-supporting and pore-floating models. Moreover, we estimate in situ hydrate saturation and clay concentration using our model, matrix-supporting model, and pore-floating model independently, and find that (i) hydrate saturations predicted by our model conform better with the saturations from chloride concentration and (ii) clay contents calculated by our model fit the best with results from smear slide analysis. This study demonstrates that our double-solid-matrix model can be an effective tool to understand the effect of heterogeneously distributed hydrates on velocities, as well as obtain accurate hydrate content in marine isotropic sediments.

SummaryThe Hainan volcanic field (HNVF) is one of China’s most active Holocene volcanic areas. Due to a lack of comprehensive geophysical research, questions persist regarding the deep magma system of the HNVF. For example, it is unclear whether the intense seismic activity in its eastern part may be a precursor to renewed volcanic activity. We present new three-dimensional electrical conductivity images, derived from magnetotelluric data, that provide a new understanding of the deep magma system in the HNVF. Our results reveal the presence of multiple sets of low-resistivity structures in both shallow and deep regions. Although once associated with past volcanic activity, a widespread shallow low-resistivity layer on the northwest side of the HNVF is not currently indicative of shallow magma chambers. Instead, a deeper large-volume low-resistivity structure in the western part of the HNVF may represent the current crustal magmatic plumbing system. Our analysis suggests that the intense seismic activity in the east of HNVF lacks corresponding low-resistivity structures, which indicates that there is no direct correlation between seismicity and movement of magma. Recent volcanic eruptions are primarily concentrated near the Changliu-Xiangou fault, which may indicate that the migration of magma has utilized crustal weak zones.