Geophysical Journal International

SummaryFault systems have geometrically complex structures in nature, such as stepovers, bends, branches, and roughness. Many geological and geophysical studies have shown that the geometrical complexity of fault systems in nature decisively influences the initiation, arrest, and recurrence of seismic and aseismic events. However, a vast majority of models of slip dynamics are conducted on planar faults due to algorithmic limitations. We develop a 3D quasi-dynamic slip dynamics model to overcome this restriction. The calculation of the elastic response due to slip is a matrix-vector multiplication in boundary element method, which can be accelerated by using Hierarchical Matrices. The computational complexity is reduced from O(N2) to O(Nlog N), where N is the number of degrees of freedom used. We validate our code with a static crack analytical solution and the SEAS benchmark/validation exercise from Southern California Earthquake Center. We further employ this method on a realistic fault system with complex geometry that was reactivated during the 2023 Kahramanmaraş–Türkiye doublet earthquakes, generating slip sequences that closely match real observations.

SUMMARYThe spherical harmonic coefficient Level-2 products of the Gravity Recovery and Climate Experiment (GRACE) mission are affected by north-south stripe noise. Toward this end, we have developed a new filter named Variational Mode Decomposition spatial (VMDS) filter that transforms the Equivalent Water Height (EWH) map derived from GRACE level-2 product into a one-dimensional sequence, which is then filtered by using variational mode decomposition. This approach overcomes the limitations of the singular spectrum analysis spatial (SSAS) filter, which well performs in the medium-frequency band but omits the high-frequency NSS noise. We thus put the VMDS filter behind the SSAS filter to utilize the good performance of SSAS in the medium-frequency band and thus propose a combined filter termed SV. A closed-loop simulation demonstrates the better ability of SV to suppress NSS noise and preserve signal at the grid scale compared to the SSAS filter. In the real-world scenario, the SV solution achieves a noise level (46.68 mm of EWH) below that for SSAS and DDK7 solutions (53.52 and 53.68 mm of EWH, respectively) over the ocean at low latitudes. Moreover, the well-documented water level of Lake Victoria and the well-modeled coseismic gravity change of the Mw9.2 2004 Sumatra-Andaman earthquake demonstrate that the SV filter efficiently preserves localized mass evolutions while suppressing north-south stripe noise. Such short-wavelength signals usually miss in highly filtered spherical harmonics (e.g. DDK5 and DDK6) solutions or are significantly inconsistent for various mass concentration solutions.

SUMMARYWe propose a stable and efficient method for high-degree regional lithospheric magnetic field modeling based on spherical Slepian functions, achieving significant improvements in the computational efficiency of solving large linear systems by reducing both complexity and memory requirements. This method leverages the orthogonality of Slepian basis functions on regional domains R to improve the stability of regional modeling. The block-diagonal structure of the normal equation matrix and the sparse representation by Slepian functions are simultaneously exploited to achieve a two-stage matrix compression. Additionally, a Bayesian Information Criterion (BIC)-based strategy is introduced to determine the optimal truncation number for the Slepian basis functions, ensuring a balance between computational efficiency and data fidelity. Our method was validated by directly inverting synthetic regional lithospheric magnetic field data generated from a spherical harmonic model up to degree 1050. The modeling process showed a reduced memory requirement of approximately eight orders of magnitude. The high-degree model successfully reconstructed the desirable magnetic field at different heights, and the residual statistical analysis results showed that the spatial variation of the lithospheric magnetic field was accurately captured. In addition, the proposed method can be readily extended to gravity modeling and other applications that utilize spherical harmonic analysis. The matrix compression techniques adopted provide an ideal framework for parallel computing, showing a wide range of application potential.

SUMMARYSeismic ambient noise tomography has been a powerful tool for seismic imaging, but most existing approaches fail to accurately predict detailed correlation waveforms in the presence of spatially heterogeneous noise distributions. Full-waveform ambient noise inversion allows for high-resolution waveform-based inversion even in substantially heterogeneous noise fields. Unfortunately, the computational cost of this approach is constrained by the noise distribution, rather than the inversion domain: global-scale applications are viable, but smaller-scale applications require a modified approach. We present a general finite-domain full-waveform ambient noise inversion methodology, providing approximation mechanisms for treating out-of-domain propagation from distant noise sources. This makes the problem tractable on much smaller domains. In a numerical example, we demonstrate that this approach enables inversion for both structure within a chosen domain and approximate noise source distribution outside it.

AbstractThe British Isles' lithospheric structure, shaped by a dynamic geological history, remains incompletely understood, particularly regarding anelastic parameters, such as attenuation. In this study, we present a teleseismic attenuation model for the British Isles, using time-domain analysis of teleseismic P-wave data from 28 deep earthquakes. We constructed a 2D differential attenuation map (Δt*) that reveals significant regional variations. Our findings show a weak anticorrelation between Δt* and shear wave velocity at upper mantle depths, suggesting that variations in the lithosphere-asthenosphere system influence this pattern. A high-attenuation zone extends from Scotland across the Irish Sea to southwest England, potentially linked to mantle upwelling associated with the Iceland plume. This model provides new insights into the mantle dynamics beneath the British Isles, offering a crucial reference for future geophysical studies in the region.

SUMMARYInduced seismicity poses a significant challenge to the safe and sustainable development of Enhanced Geothermal Systems (EGS). This study explores the application of machine learning (ML) for forecasting cumulative seismic moment (CSM) of induced seismic events to evaluate reservoir stability in response to fluid injections. Using data from the Cooper Basin (Australia), the St1 Helsinki geothermal project (Finland), and a controlled laboratory injection experiment, we evaluate ML models that integrate catalog and operational features with various frameworks. Results indicate that feature-rich models outperform simpler ones in complex seismic environments like the Cooper Basin and laboratory cases, where seismicity is promoted by earthquake interaction and fault reactivation. However, in scenarios like St1 Helsinki, with minimal event clustering, additional features offer limited predictive benefits. While ML models are promising, several challenges impede reliable forecasting, including data scarcity from operational wells, the extrapolation demands of cumulative output (i.e. CSM), and the difficulty of predicting abrupt CSM increases for large seismic events. Enhancing model robustness requires synthetic data augmentation and improved feature selection capable of capturing diverse reservoir dynamics. These advancements may enable more accurate near real-time forecasts of problematic induced seismic events, informing operational decisions to mitigate seismic risks while maximizing energy extraction, and hence offering a pathway for broader adoption of ML in renewable energy development and management.

SummaryFor the inversion of crustal deformation data, where both the fault geometry and slip distribution need to be estimated, attempts have been made to estimate them simultaneously within the framework of Bayesian inference. In these methods, an analytic expression for the posterior distribution of the parameters cannot be obtained depending on the settings of the probability distribution. The Monte Carlo method is often used to approximate the probability distribution by several samples. However, Monte Carlo sampling becomes computationally expensive as the dimensionality of the parameters increases, making it challenging to apply these methods to inverse analyses targeting large-scale earthquakes with complex rupture patterns, where the parameters to be estimated are high-dimensional. In this study, we developed an efficient algorithm for the simultaneous Bayesian estimation of fault geometry and slip distribution that is robust against high-dimensional parameters, facilitating the estimation of multi-segment fault geometry and complex slip distribution. When comparing the developed method with a conventional method in a numerical experiment where the parameters were high dimensional, it was confirmed that the accuracy and convergence of the estimation results were improved while reducing computational costs. As an application example of the developed method, we estimated fault geometry and slip distribution from the crustal deformation observed during the 2024 Noto Peninsula earthquake, considering models with one-, two-, and three-fault segments. By comparing the results, we demonstrated that the three-segment model is the most plausible, emphasising the importance of improving the estimation methods to accurately estimate multi-segment fault geometry with a complex slip distribution.

SummaryThe characterization of submarine volcanism associated with extension along continental margins is essential for understanding the evolution of rifting basins. In the Okinawa Trough, between the Eurasian continent and the Ryukyu Arc, submarine volcanism remains poorly understood. Here, we conducted detailed magnetic and gravity surveys aboard research vessels off Kume Island in the mid-Okinawa Trough to clarify the distribution of different types of submarine volcanic activities. Equivalent magnetization intensities were estimated from observed magnetic anomalies, and Bouguer gravity anomalies, calculated from the obtained free-air gravity anomalies, were used to estimate crustal thickness. Differences in crustal equivalent magnetization, Bouguer gravity anomaly, and crustal thickness among bathymetric highs allowed us to identify four distinct groups of submarine volcanic edifices reflecting the evolution of rifting. Edifices of Group 1 are affected by basaltic magma intrusions associated with back-arc extension. Those of Group 2 have moderate equivalent magnetizations and Bouguer gravity anomalies compared to other edifices in the study area; the equivalent magnetization intensity in this group is comparable to the average intensity of the oceanic lithosphere. Edifices of Group 3 may have been affected by additional, voluminous magma supplies from active island arc magma sources. Finally, those of Group 4 are characterized by fewer basaltic intrusions compared to other groups and a lower degree of crustal thinning associated with back-arc extension. Moreover, the faults and fractures formed during the development of the Kerama Gap may have affected volcanism in our study area.

AbstractMacroseismic intensity classifies the ground shaking at a locality by comparing the observed effects on humans and buildings with the scenarios characterising each intensity degree according to a macroseismic scale. This practice may involve uncertainty in assessing intensity degree due to several factors. These uncertainties propagate to subsequent elaborations, such as the parameters of pre-instrumental earthquakes determined from macroseismic data. In Italy, more than 60 per cent of the earthquakes in CPTI15, the Italian Parametric Earthquake Catalogue covering the period 1000–2020, rely on intensity data collected in DBMI15, the Italian Macroseismic Database. Their parameters are estimated with the ‘Boxer’ software, which determines the location and magnitude starting from their macroseismic intensity distributions. In this work, we explore the potential impact of possible inaccurate intensity assessments at a single site on macroseismic parameters (i.e. locations and magnitudes) of Italian earthquakes. We select 1108 earthquakes with at least 10 intensity data points from CPTI15 that occurred in the period 1279–2020. For each event, we simulate more than 100 sets of intensity distributions, for a total of 138.327 simulations, by varying the intensities at the sites of ± 1 with a half-degree step starting from the intensity of DBMI15. Each simulated distribution is then parameterised using the same approach adopted by CPTI15 (i.e. Boxer), and the results are compared with the macroseismic epicentre and magnitude of CPTI15. The resulting parameters from all the simulated distributions are coherent with those provided by CPTI15. Locations estimates are within 5 km from CPTI15’s for 55 per cent of the cases, and within 10 km for 83 per cent. The magnitudes of 68 per cent of simulations are within ± 0.2 units of difference from the CPTI15 magnitudes and 87 per cent within ± 0.3 units, similar to the statistical error of instrumental magnitude estimation. Moreover, we treat uncertain intensity values (i.e. 6–7) as equally representative of either the lowest intensity level (i.e. 6–7 as 6) and the highest intensity level (i.e. 6–7 as 7), and we analyse their impact on the parameters. The differences with the CPTI15 magnitude are not significant for either analysis, with more than 97.7 per cent of the earthquakes falling within ± 0.3 units of difference.

AbstractHigh-precision geoid models have traditionally been static, neglecting temporal variations. However, achieving geoid accuracy within 1-2 cm and maintaining dynamic height reference frames necessitates consideration of geoid spatiotemporal variations. GRACE/GRACE-FO and surface mass loading models provide means to estimate geoid changes, but their accuracy and reliability require further validation. This study proposes a method for dynamically maintaining regional height reference frames by integrating GNSS reference stations as core nodes and incorporating time-varying geoid data. This method dynamically corrects station heights by computing normal height variations using GNSS observations and geoid changes. Experiments conducted in Beijing and Shandong derived geoid changes using GRACE/GRACE-FO and surface mass loading, validated against long-term GNSS observations and leveling surveys. Results show a strong correlation (R ≈ 0.9; NSE > 0.4) between geoid changes derived from GRACE and surface mass loading, although amplitude discrepancies of up to 4 mm existed. In 41 experimental cases, accuracy improvement was observed in over 90 per cent of instances following geoid change corrections. In Beijing, 18 out of 26 results achieved accuracy improvements exceeding 20 per cent, five of which surpassed 90 per cent. In Shandong, 11 out of 15 results improved by over 10 per cent, including five exceeding 40 per cent. These findings confirm the feasibility and effectiveness of using GRACE/GRACE-FO and surface mass loading to estimate geoid changes. The proposed method significantly improves the accuracy of dynamic height reference frame maintenance, providing valuable insights for further refinement of geoid models.

SummaryThe Cameroon Volcanic Line (CVL) and other tectonic features in Cameroon remain enigmatic, prompting ongoing debates about their detailed crustal structure, composition, and geodynamic evolution. To shed light on the structural complexities and the underlying crustal processes, we leverage the two-step ambient noise tomography (ANT) method to obtain the 3D shear wave velocity (Vs) and the P-to-S wave speed ratio (Vp/Vs) structure of the crust beneath Cameroon. We start by cross-correlating data recorded at 32 broadband stations from February 2006 to February 2007 to extract Rayleigh wave group dispersion curves on inter-station paths. First, we invert these dispersion curves to obtain group velocity maps across different periods (5-30 s) on a regular grid (0.5o x 0.5o). We then invert the group velocities on each grid node to derive Vs and Vp/Vs as a function of depth. Specifically, we leverage a new evolutionary algorithm called Competitive Particle Swarm Optimization (CPSO) to tightly constrain Vs and Vp/Vs ratios beneath the CVL and surrounding regions. Our inversion results show an anomalously low Vs of ~3.6 km/s in the uppermost crust beneath active volcanic provinces. This low Vs and a high Vp/Vs ratio suggest a mafic composition, possibly due to mafic volcano-plutonic melts driving Cenozoic to modern magmatic activity. Our findings reveal a prominent high-Vs structure at 25-35 km depth, in alignment with the CVL. Characteristic properties, such as the maximum Vs of ~3.9 km/s and Vp/Vs in the range of 1.85-1.88, suggest the presence of cooled mafic materials that have intruded the crust. Our depth cross-sections along the CVL indicate that these mafic intrusions are ubiquitous along the entire CVL. They are spatially separated from the volcano-plutonic structures by a thin middle crust with a Vs of ~3.7–3.8 km/s and a Vp/Vs of 1.70. These properties are indicative of a felsic to intermediate crust, which may be linked to the Neoproterozoic Pan-African Orogeny. We posit that this thinned low Vp/Vs structure may have facilitated the ascent of mafic material, contributing to recent volcanic activity in the region. Conversely, beneath the Oubanguides belt and the Congo craton, these low Vp/Vs structures appear thicker, with mafic intrusions at greater depths. This structural feature suggests a dynamic process involving the pushing and exhumation of lower crustal material by the mafic material. Our images further suggest that an intriguing interaction of crust with deeper structures may be responsible for the intrusions and volcanism observed along the CVL. The findings advance our understanding of the geological and geodynamic complexities associated with the CVL and its origin.

AbstractThe seismic behavior of the near-trench plate interface of the Guerrero seismic gap and other segments of the Mexican subduction zone is likely to play a critical role in the seismic and tsunami hazard of the region. In this context, a detailed study of the near-trench 18 April 2002 Mw 6.7 earthquake that occurred about 55 km off the coast of Guerrero and generated a small tsunami attains particular importance. From an analysis of the teleseismic P waves and S waves, local recordings, and aftershock distribution, we find that the rupture most likely began at a subducted seamount, propagated unilaterally towards NW, parallel to the trench for ~ 54-58 km and a duration of ~ 68-70 s. The moment rate function is highly rugged, with two dominant pulses separated by about 50 s. Although relatively small in magnitude, the earthquake has all the characteristics of a tsunami earthquake: the slip occurs very close to the trench, the rupture speed is slow (~ 1 km/s), the high-frequency radiation is deficient, and, in common with tsunami earthquakes, the moment-scaled radiated energy is low (ER/M0 = 1.45 × 10−6). We confirm that the duration of the event (~ 70 s) is extraordinarily long compared to that expected from scaling relations (~ 12.8 s), consistent with it being the most anomalous of all the events studied in the last 40 years (Duputel et al., 2013). Our results support a conditionally stable upper 15 km of the plate interface in the region reported from recent offshore seismic observations.

SummaryInduced polarization is a geophysical method that can be applied to determine the water content and CEC (Cation Exchange Capacity) of sediments and rocks. We apply here this technique to image the Harmalière landslide in Isère (East of France). This landslide includes a mudflow of glacio-lacustrine silty clays overlying the Toarcian clayrock formation with local occurrences of compact pebbles from a paleochannel of the Drac river. A petrophysical study is used to characterize the properties of the four main lithofacies occurring in the area with a total of 22 samples. We performed complex conductivity measurements in the frequency range 10 mHz-45 kHz at different salinities (NaCl brines and in situ pore waters). We also measured the porosity and CEC of the samples. We calibrate the relationships between surface conductivity, quadrature conductivity, and normalized chargeability with both the porosity and CEC. The relationship between the formation factor and the porosity conforms to Archie's law with a cementation exponent close to 2.0 ± 0.2. In the field, we performed a time-domain induced polarization survey using a 1.26 km-long cable (including roll along of the electrodes) with an electrode spacing of 20 m. The landslide is imaged down to a depth of 220 m. The inversion of the data (788 electrodes, 15.7 km of profiling, 13 012 apparent resistivity and 5539 apparent chargeability data) is done with the least-square technique penalizing the roughness of the tomograms using Occam inversion. The resulting 3D electrical conductivity and normalized chargeability tomograms are analyzed in conjunction with the petrophysical data to image the extension of the lithofacies at the field scale. Furthermore, the water content and CEC of the formations are imaged. We demonstrate that the compacted pebbles of the Drac paleochannel form both a mechanical and hydraulic barrier that is locally breached by the mudflow before entering the Monteynard-Avignonet Lake. This study demonstrates the ability of induced polarization to finely characterize the anatomy of such landslide and image its water content.

AbstractThis study investigates the following three issues in numerical models of the thermal structure of subduction zones, using the Tohoku region in Northeast Japan as an example: (1) a steady state is often assumed in models, (2) quantitative assessment of the uncertainty in the predicted temperatures is lacking, and (3) surface heat flow has been used to constrain many of the models. I found that, at least under the model setting of this study, a steady state may be safely assumed as long as only surface heat flow within 150 km of the trench is used to constrain the model. I used Bayesian inference to predict the thermal structure, with surface heat flow near the trench and the location of the blueschist-out boundary in the oceanic crust as observational constraints. The depth of slab–mantle kinematic decoupling, effective friction coefficient, and rate of radiogenic heat production in the upper island arc crust were constrained simultaneously to be ∼80–100 km, 0.03–0.08, and 1.5–2.16 μW m−3, respectively, although the decoupling depth is sensitive to the assumed location and temperature of the blueschist-out boundary. The uncertainties in slab temperature reach ∼450 K at depths of <100 km and 100 K for greater depths, which are substantial. To reduce these uncertainties, it is necessary to reduce the uncertainty in the input parameters and obtain additional observational constraints.

SummaryWater content and pore fluid pressure increases have been recognized as important drivers of shallow landslides, especially through the role of strong rainfalls promoting gravitational instabilities. Less recognized is the role of vertical hydraulic barriers impeding the flow of ground water at the feet of areas prone to landslides. Induced polarization is a non-intrusive geophysical technique able to image hydraulic properties of the shallow subsurface. Recently developed petrophysical models bridging the gap between hydraulic and electrical properties of soft sediments, soils, and rocks have been developed. Thanks to these relationships, this geophysical method can be used to image the water and clay contents of the formations and their permeability. Therefore, induced polarization can be used to image the occurrence of vertical permeability barriers. We focus our approach on a large landslide that occurred in March 1931 (reactivated in 1971–1972) above Le Châtelard village (Bauges, France). This landslide started inside a kilometer-scale syncline hosting clayey formations and moraines. We performed a 2.2 km profile crossing the syncline and the sliding area including resistivity, induced polarization, and self-potential measurements. In addition, 22 samples were taken from the different formations outcropping at the field site including limestones, sandstones, and clayey formations. The petrophysical investigations are combined with the field data to image the water content and cation exchange capacity as well as their permeability. The dataset shows the existence of a vertical permeability barrier at the bottom of the landslide corresponding to the tight Urgonian limestone formation. We combine the permeability distribution, the resistivity, and self-potential data by forward modeling the groundwater flow and electrokinetic response. We then invert the self-potential measurements to refine the image of the Darcy velocity distribution. The results show a strong upflow of the ground water just above the Manauds canyon where several gravitational instabilities occured in the past.

SummaryHere, we investigate how continental rifts initiate and propagate across cratons by exploring the crustal structure of northwestern tip of the East African Rift System (EARS), hosting the volcanic-rich Edward-George and non-volcanic Albertine-Rhino rifts, and their termination at the Precambrian Aswa Shear Zone. We conducted a derivative analysis of magnetic data, utilized power spectral analyses, and implemented a two-dimensional (2D) forward modeling of gravity data constrained by the seismic results obtained from the region. A magnetic derivative map indicates that the border faults of the Albertine Rift, at regional-scale, trend parallel to the Mesoproterozoic Madi-Igisi fold belt (MIFB) structures, representing the suture zone between two Archean microcratons. Our results show a pronounced thinned crust (∼24–30 km) beneath the southern segments of the rift zone, particularly the Edward-George rift, the Rwenzori Mountains, and the southern Albertine graben, consistent with previous seismic studies. In general, we observe that: 1) the rift system follows the boundary between a broadly thinner crust (21–41 km) to the southeast in Uganda, and thick crust (34–41 km) to the northwest in Congo, and 2) within the rifts, the crustal thickness along the axes exhibits a strong gradient that attenuates northwards beneath the Albertine-Rhino graben. We supplement the geophysical results with field observations of an exhumed Permian ‘Karoo’ rift (Entebbe Graben) in central Uganda, indicating the possible source of inherited thinner crust to the southeast of the Albertine-Rhino Rift. We propose that the northwestern tip of the EARS exploited a cratonic crustal thickness-gradient, assisted by structural inheritance from crustal metamorphic fabrics, and potentially, thermo-mechanical weakening of the deeper crust by partial melts beneath some of the rift segments.

AbstractDeep carbonate reservoirs are characterized by complex stress states and pore structures, and frame compressibility is strongly influenced by differential pressure. We have developed a frequency-dependent acoustoelasticity-squirt flow model based on a threshold pressure to describe the frame compressibility. Based on this model, 3D rock physics templates (RPT) were created and calibrated with laboratory, well log and seismic data from the Longwangmiao Formation in the Gaoshiti-Moxi area, Sichuan Basin, China. The templates are then used to determine the reservoir properties at the well logging and seismic frequencies. The predicted porosity matches well with the actual porosity log curve, and the results are consistent with gas production reports, with high porosity, crack density and threshold pressure (high frame compressibility) indicating high reservoir potential.

SUMMARYArchaeomagnetic study serves as an effective way to investigate the detailed variations of the geomagnetic field during the Holocene. China has accumulated numerous high-quality archaeomagnetic data in recent years, but these data are primarily limited to the period younger than 1500 BCE, which limits our comprehensive understanding of the Holocene geomagnetic field evolution. Here we carried out archaeomagnetic research on 28 unoriented pottery shards and baked clays collected from two Neolithic archaeological sites, with ages spanning from ∼5350 to 5030 BCE, in Anhui province, Eastern China. Rock magnetic analyses show that pseudo-single domain magnetite is the main magnetic carrier among the samples with limited contributions from high coercivity minerals such as hematite. The ideal magnetic behaviour demonstrates the samples are suitable for palaeointensity experiments, as confirmed by the palaeointensity results. A total of 14 samples yielded high-fidelity palaeointensity results, indicating that a fast-changed geomagnetic field with virtual axial dipole moments varying from ∼56 to ∼91 ZAm2 over ∼300 years. This inference is further evidenced by intensities recovered from the double-remanent components of some samples through vector calculation. Our new data provide reliable anchors for the strength variations of the geomagnetic field in Eastern Asia during the poorly constrained time period before 4000 BCE. With the new data in this study, we updated the previous Chinese archaeointensity reference curve and named the new curve as ArchInt_China3, which extends the previous curve forward for ∼500 years. The newly reported archaeointensity data and released reference curve in this study will provide valuable insights into the regional and global geomagnetic field variations during the Holocene, and thus assist in understanding the dynamic processes in the Earth's interior.

SUMMARYInterpretation of palaeomagnetic data requires the detection of magnetofossils in sedimentary rocks and an understanding of their influence on magnetic properties. Subsamples collected from IODP site M0061 lost up to 90 per cent of their initial bulk magnetic susceptibility (MS) during cold-room storage of 4 months, which was attributed to the alteration of single-domain magnetosomal greigite (Fe3S4). To test if the magnetic susceptibility loss affected the anisotropy of MS (AMS) we resampled site M0061 with a Kullenberg piston corer (3 cores), took palaeomagnetic subsamples and undertook time-dependent AMS measurements over 1 year in a controlled cool, humidified environment and exposed to air. Most subsamples possessed an initial normal oblate AMS fabric predicted for laminated sediments (horizontal with respect to the bedding plane) but we also detected a negative trend between the degree of anisotropy (Pj) and MS. In accordance with previous observations, MS decreased over 1 year, which we accredit to oxidation of magnetosomal greigite and conversion into a less magnetic phase (probably FeO(OH)) that does not make a detectable contribution to AMS. These results allowed us to isolate, through application of an AMS tensor subtraction routine, the fabric of the magnetosomal greigite component that had decayed. In subsamples with the largest MS loss over one year, the decayed component had a prolate, inverse AMS fabric (defined as the principal susceptibility axis perpendicular to the bedding place) but relatively low Pj. We conclude that the initial (in-situ) AMS ellipsoid consisted of a mixture of a typical normal, oblate sedimentary fabric and the prolate, inverse magnetosomal fabric. The mostly inverse nature of the separated fabric indicates that the long axis of the magnetosomal greigite (as individual single-domain magnetofossils or chains) must be oriented parallel to the bedding plane, which implies that the magnetosomal greigite was deposited from the water column and contributes to a depositional remanent magnetisation (DRM). Our results indicate that greigite magnetofossils can (i) explain the inverse AMS fabrics that have been reported in similar sedimentary environments and (ii) carry DRM with a median destructive field (MDF) of approximately 20 mT, although this remanence is transient under ambient laboratory conditions and is prone to oxidation.

SummaryEstimating quantitative hydrogeological information from geophysical datasets remains a key challenge in hydrogeophysics, driving the development of innovative inversion methodologies. Among these, coupled hydrogeophysical inversion (CHI) is a promising approach that integrates hydrological and geophysical modeling to improve hydrological property estimations from geophysical observations. However, most CHI applications focus on time-lapse geophysical datasets, while applications using single-time geophysical datasets “historically far more common in hydrological studies” remain scarce. Moreover, CHI also depends on petrophysical relationships, whose accurate calibration is challenging, leading to uncertainties that significantly affect CHI results and must be accounted for. This work proposes Hybrid Bayesian Inversion (HBI) to implement CHI using single time geophysical dataset, which considers coupled hydrological-geophysical modelling constraints and petrophysical relationships, including several calibration uncertainties. HBI is based on solving the hybrid decomposition of subsurface geophysical properties. This decomposition is derived as the sum of the groundwater model, directly predicted by coupled geophysical-hydrological modeling, and the background model, which accounts for the residual geophysical properties not predicted by coupled modeling. The groundwater model is formulated as a stochastic model characterized by probability density functions (PDFs), which enable the derivation of estimations for posterior conditional PDFs of hydrological and geophysical properties. In contrast, the background model is characterized only by estimating its values (e.g. using maximum likelihood estimators) and its PDFs are not determined. The hybrid decomposition in HBI is solved using the Expectation-Maximization (EM) algorithm. This approach divides the iterative solution into sequential steps of Bayesian inversion (E-step) and classic least-squares geophysical inversion (M-step). This formulation allows Bayesian inversion, which is computationally intensive, to focus only on relevant variables linked to hydrological conceptualization (Groundwater model), while classic least-squares geophysical inversion is used to solve the remaining variables (Background model). The HBI methodology was tested using 2D ERT synthetic and experimental data from an unconfined aquifer. In these examples, ERT modelling is integrated with saturated groundwater flow modeling through uncalibrated Archie and CK (electrical conductivity to hydraulic conductivity) petrophysical relationships. Results indicate that even with significant calibration uncertainties in petrophysical relationships, HBI can recover valuable information regarding water table and water conductivity what is not directly derivable from classic least-squares inversion results. Additionally, results derived from experimental data show that HBI can be an effective method to discriminate between low resistivity caused by fine-grain content and the water-saturated zones.