Geophysical Journal International

AbstractDeep learning (DL) approach has gained attention for earthquake (EQ) detection. To alleviate the problem of training data shortage, transfer learning (TL) provides a useful framework to adapt pre-trained models, typically through tuning of model parameters. Nonetheless, the current practice still requires considerable data, which hinders its application where only a small number of data is available. Instead of TL, we propose a novel two-stage of model correction as a solution to this important and ubiquitous problem in EQ detection. In the proposed approach, a pre-trained DL model is directly applied to waveform data in the target domain (first stage), and the cases that are classified as an earthquake signal (i.e., positive cases) are further classified as positives and negatives using a non-DL classification method (second stage). Our classification method for the second stage is based on multiple clustering, which characterizes local waveform patterns in terms of amplitude scale in specific time segments that are inferred in a data-driven manner. This characterization captures complex high-dimensional waveform patterns in a low-dimensional space, which leads to the effective classification of true and false positives. Furthermore, the proposed method is useful when only true positive waveforms are labeled (PU classification). Both synthetic and real data analysis clearly demonstrated effectiveness of unsupervised waveform characterization of the proposed method.

AbstractFull-waveform inversion (FWI) is a powerful imaging technique that produces high-resolution subsurface models. In seismology, FWI workflows are traditionally based on seismometer recordings. The development of fibre-optic sensing presents opportunities for harnessing information from new types of measurements. With dense spatial and temporal sampling, fibre-optic sensing captures the seismic wavefield at metre-scale resolution along the cable. Applying FWI to fibre-optic measurements requires the reformulation of the forward and adjoint problems due to two fundamental differences to seismometer data: i) fibre-optic measurements are sensitive to strain rather than translational motion, and ii) they do not represent the motion at a single spatial point, but instead capture the average deformation over a pre-defined cable segment, known as the gauge length. Within this study, we derive the adjoint sources to perform FWI for data from distributed acoustic sensing (DAS) and integrated fibre-optic sensing (IFOS) that are based on moment tensors. Our formulation incorporates gauge-length effects, direction-dependent sensitivity and complex cable layouts. For the numerical simulations, we use a spectral-element solver that allows us to incorporate surface topography, and coupled viscoacoustic and viscoelastic rheologies. In illustrative examples, we present how our theoretical developments can be used in inversions of synthetic fibre-optic data generated for a realistically curved cable placed on irregular topography. As examples, we invert for source parameters, including moment tensor, location, and origin time for noise-free DAS data, noise-contaminated DAS data, and IFOS data. Further, we present the 3-D imaging results for the three data groups and further analyse the effect of scatterers on the FWI based on DAS data. In all example inversions, we compare how close the found model is to the known ground truth. The codes to produce these results are accessible and ready to be applied to real data inversions.

AbstractHigh-frequency induced polarisation, which measures the complex electrical conductivity in a frequency range up to several hundred kHz, is potentially suitable to detect and quantify ice in the frozen subsurface. In order to estimate ice content from the electrical spectra, a two-component weighted power mean (WPM) model has been suggested and applied to field-scale data. In that model, ice is one of the components, whereas the solid phase, residual liquid water and potentially air form the second component, called “matrix”. Here, we apply the model to laboratory data previously discussed in the literature, with the aim to assess the applicability of the model and to understand the behaviour of the frequency-dependent electrical conductivity. The data were measured on an unconsolidated sediment sample with 20.8% water content from the European Alps, and a consolidated sandstone with 16.6% porosity. Electrical spectra have been measured over a temperature range from approx. - 41 ○C to +20 ○C and a frequency range from 0.01 Hz to 45 kHz. We extend the original WPM model to account for low-frequency polarisation in form of a constant phase angle model. The measured data were fitted with the model by a least-squares inversion algorithm. In order to reduce the ambiguity, we constrained several of the nine underlying parameters by literature values, in particular for the electrical properties of water ice, and the expected ice content according to porosity or water content of the unfrozen sample. Both data sets can be well matched, corroborating the hypothesis that the model is in principle suitable to explain measured data of frozen samples in that frequency range. One important observation is that the mixing parameter, i.e. the power in the WPM model, which is controlled by the geometric arrangement of the two components, depends on temperature. For the unconsolidated sample it even becomes negative at the coldest temperature, which is important because negative shape factors relate to specific geometries. A second observation is that relatively large permittivities of the matrix are required to fit the data, suggesting that processes at the interface between solid/liquid phase and ice, which are not included in the volumetric mixing model, might be relevant and should be considered in future extensions of the model.

SummaryThe dominant forces shaping the unique geometries of flat slabs are still not fully understood. Knowing how the stress field changes with respect to the shape of the slab allows inferences of the dominant forces acting on the slab. In this study we calculated new models of the slab geometry and the intraslab stress field in the Pampean flat slab region of the Chile-Argentina Subduction Zone (latitude ∼25°36°S) where the Nazca Plate subducts together with the aseismic Juan Fernandez Ridge. To build the models, we used a catalog of 1,059 well-located slab earthquakes recorded by the SIEMBRA and ESP temporary seismic arrays and calculated 411 new focal mechanisms that were analyzed together with 407 focal mechanisms from other catalogs. Our results confirmed slab seismicity features such as a reverse dip (i.e. opposite to the subduction direction) of the seismicity band within the flat slab, two bands of descending seismicity, and two regions with an absence of earthquakes. These seismicity patterns express the shape of the slab and its hydration state, with more localized slab dehydration along the inland path of the Juan Fernandez Ridge relative to the surroundings. In one of the regions without earthquakes, the slab is most likely continuous and dry, while in the other one the slab is missing, in agreement with previous works that proposed a hole in the slab visible with other methods. A comparison between the stress field and the local slab dip from both our new model and a previous one (Slab2) indicates that the dominant forces acting on the flat slab are the slab pull and the ridge buoyancy. Finally, the shape of the flat slab is controlled by the geologic migration of the Juan Fernandez Ridge, making the flat slab four times wider than the ridge offshore, and by the competing forces of the slab pull and the ridge buoyancy that creates a notable flexure (bulge) resembling the geometry of the outer rise near the trench.

SummaryIn the analysis of Induced Polarization (IP) data, it is commonly assumed that induction effects (IE) are negligible. However, at higher frequencies, this assumption becomes increasingly invalid, posing challenges for IP measurements. High-frequency induced polarization (HFIP) extends the conventional IP frequency range beyond 100 kHz, allowing estimation of ice content by capturing the characteristic decrease in permittivity of water ice. This study focuses on the interpretation of HFIP data while accounting for IE. We modified an existing one-dimensional simulation code to evaluate HFIP responses over frozen ground with ice, both with and without the influence of IE. Our results demonstrate that IE can distort HFIP measurements in typical permafrost conditions, potentially obscuring the characteristic polarization behavior of water ice. two-dimensional IP inversion codes that account for IE are not routinely available. Even if they were, the simultaneous presence of IE and IP would likely complicate their application. We therefore propose to remove IE from the data, with the aim to use well-established two-dimensional SIP inversion routines. For this purpose, we implemented a one-dimensional inversion routine that includes IE as a frequency-dependent correction factor. The method is based on fitting a layered model to the data. Comparing the responses with and without IE, we calculate a correction factor which is subsequently used to remove IE from the original dataset. The method is conservative in the sense that features not well matched by the one-dimensional inversion are preserved and no information gets lost. We demonstrate the effectiveness of the method with synthetic data, as well as field datasets from an unfrozen site in Germany, and from permafrost peatlands in Scandinavia. We provide evidence from reciprocal measurements that cable coupling effects, not included in the correction procedure, have been effectively minimized by the acquisition system. We further show that high-frequency phase shifts are strongly influenced by IE, and that the correction methodology successfully restores the spectral response. By applying the approach to a measured dataset, we demonstrate that two-dimensional inversion of the corrected data with a well-established code is both feasible and robust. The resulting model deviates markedly from that obtained with uncorrected data, highlighting the critical role of the correction procedure for reliable interpretation.

SummarySubduction zones are tectonic plate boundaries where one tectonic plate is being forced beneath another and known for producing some of the most powerful earthquakes in history. Understanding the regional-scale structural heterogeneity in the subduction zone is crucial for deciphering the genesis of megathrust earthquakes and the factors that control rupture dynamics. The rupture dynamics of the earthquakes in Andaman-Sumatra region are primarily hindered due to the lack of velocity and anisotropy information at the lithospheric level. Compression wave (Pn) waves propagate in the uppermost mantle and can provide velocity and anisotropy constraints for this portion of the mantle along their ray-path. We generated the first comprehensive high-resolution Pn-wave tomography, restricting the turning of diving ray paths less than 50 km deep in the collision zones, to map lithospheric velocity and anisotropy by inverting 65 297 Pn arrivals extracted from 6958 regional events recorded at 384 stations. The results show a strong variation in Pn-wave velocity and fast polarization directions (FPDs) in the entire study region that may be an essential factor for earthquake nucleation. We observed an abrupt change in the Pn-FPDs from Andaman to Sumatra and Sumatra to Java regions. We suggest that well defined plate boundaries between Indian and Capricorn plates, that was earlier reported to be diffused one. Such observation is also supported by the lower spreading rate and higher crustal age of the old Indian oceanic plate in the Andaman region, compared to the faster spreading rate and lower crustal age of the Capricorn oceanic plate in the north Sumatra region. We suggest that the Sunda plate strongly coupled with subducting oceanic slab causes the mega events of 2004 and 2005 supported by trench-normal Pn-FPDs along with rigid-slab nature, and both the rupture propagation terminate near Simeulue Island region due to abrupt material changes beneath it. Our study also reveals the presence of low Pn velocity anomalies in the west of the Andaman Sea, suggesting the existence of a magma reservoir pouring out lava through the steeply torn Indian oceanic slab.

SummaryThe converted wave technique, namely Receiver function (RF), has been routinely employed for estimating one-dimensional velocity models of the Earth’s crust and mantle structures. Physics-driven methods such as inversions are employed to receiver function data to estimate velocity models. Despite their wide utility, the presence of dipping and anisotropic geological structures often complicates this process. To address these complexities, here we introduce a Physics-Guided unsupervised deep learning approach for the inversion of receiver function data. An unsupervised deep learning approach together with implicit neural representations is developed here, for prediction of subsurface model parameters such as the layer thickness, S-wave velocity, anisotropy, trend, plunge, strike and dip without requiring any labelled data. In addition, we incorporate a dynamic layer setting criterion to automatically pick the optimal number of layers required to fit the data, which is found to be particularly effective for identifying pronounced discontinuities like the crust-mantle boundary. For determining the optimal model parameters, neural network output parameters (the model parameters) are used in forward modelling of the receiver functions. Inversion results from both synthetic and real field data from the Indian shield and Hi-CLIMB network suggest that the physics-guided unsupervised deep learning approach is effective in inversion tasks, particularly when dealing with dipping and anisotropic media.

SummaryWe explore the potential of utilizing Distributed Acoustic Sensing (DAS) for Back-projection (BP) to image earthquake rupture processes. Synthetic tests indicate that sensor geometry, azimuthal coverage, and velocity model are key factors controlling the quality of DAS-based BP images. We show that mitigation strategies and data processing modifications effectively stabilize the BP image in less optimal scenarios, such as asymmetric geometry, narrow azimuthal coverage, and poorly constrained velocity structures. We apply our method to the Mw7.6 2022 Michoacán earthquake recorded by a DAS array in Mexico City. We also conduct a BP analysis with teleseismic data for a reference. We identify three subevents from the DAS-based BP image, which exhibit a consistent rupture direction with the teleseismic results despite minor differences caused by uncertainties of BP with DAS data. We analyze the sources of the associated uncertainties and propose a transferrable analysis scheme to understand the feasibility of BP with known source-receiver geometries preliminarily. Our findings demonstrate that integrating DAS recordings into BP can help with earthquake rupture process imaging for a broad magnitude range at regional distances. It can enhance seismic hazard assessment, especially in regions with limited conventional seismic coverage.

SummaryShort-duration instrumental earthquake catalogues, sparsity in seismic stations, and paucity of near-source earthquake data impose strong limitations on data-driven seismic hazard assessment. In contrast, regional-scale multi-cycle earthquake simulations that generate realistic rupture scenarios provide physics-based earthquake rupture forecasts for seismic hazard assessment. In this study, we use the physics-based multi-cycle earthquake simulation engine MCQsim to compute long-term synthetic seismic catalogues that include multi-segment ruptures on a complex-geometry fault system. We expand on previous research by conducting systematic statistical analyses of synthetic earthquake catalogues for the Gulf of Aqaba (GoA, Saudi Arabia) that forms the southern extension of the Dead Sea Fault (DSF). The GoA is characterized by predominantly left-lateral strike-slip faulting and its seismic activity is revealed by past seismic swarms and large-magnitude events, including the 1995 M 7.2 Nuweiba earthquake. To address epistemic uncertainties in the seismic source characterization (SSC) of this region, we explore several modelling realisations by altering the fault system’s geometric configuration and frictional properties. Our simulations reveal that multi-segment ruptures reaching M 7.6 are possible in the GoA. The simulated catalogues reveal source-scaling properties of rupture area and stress drop consistent with global observations and empirical scaling laws. Additionally, the statistical properties of the synthetic catalogues, including the earthquake magnitude-frequency distribution, align with the short-term recorded seismicity in the study region. In summary, our simulation-based study provides insights into the GoA’s seismic behaviour through comprehensive parameter-space exploration and sensitivity analyses that document the possibility for geometrically complex, large multi-segment magnitude earthquakes.

SummaryFull waveform inversion (FWI) is a high-precision subsurface imaging technique that inverts subsurface parameter models by minimizing the discrepancy between observed and synthetic seismic data. However, complicated wave propagation mechanisms, non-convexity of the loss function, and limited seismic acquisition system necessitate the incorporation of sufficient prior and physical constraints to alleviate the ill-posedness and cycle-skipping. Although the implicit FWI (IFWI) can encode implicit spectral bias (i.e., inverting model parameters from low frequencies to high frequencies) to reduce the dependency on an accurate initial model, its limited high-frequency inversion capability results in thousands of iterations for the final results. In this paper, we indicate that the frequency hyperparameters of the sine activation function in IFWI modulate the spectral bias, making a trade-off between inversion accuracy and stability, i.e., lower frequencies yield robust FWI but lower accuracy, while higher frequencies achieve higher accuracy on the premise of an accurate initial model. To improve both the stability and accuracy of IFWI, we propose a novel implicit FWI method with an adaptive Fourier reparameterization strategy (termed FR-IFWI), which explicitly encodes multi-frequency information by reparameterizing the network weights using a learnable coefficient matrix and fixed Fourier frequency bases. The role of learnable matrices in neural networks can evolve from determining frequencies in IFWI to actively selecting frequencies from fixed frequency bases through FR-IFWI, which alleviates the dependence on activation function frequency and obtains more robust and accurate inversion results. Extensive numerical experiments on the modified Marmousi, 2D SEG/EAGE Salt and Overthrust models confirm that FR-IFWI successfully achieves superior inversion efficiency and accuracy compared with conventional FWI and IFWI methods.

SummaryA growing body of literature has contributed to linking the presence of bacteria with SIP signals. Yet, there are still unresolved questions concerning the contribution of cell density and microbial metabolic activity in porous media (soils and sediments) to SIP signals. Moreover, there is continued debate on whether cells themselves polarize or whether a cell-sediment interaction is a prerequisite for the measured responses. This study investigates the SIP response of Shewanella oneidensis MR-1 in isolation, that is, in the absence of a mineral porous medium using two setups (i) cells in aqueous suspension and (ii) alginate bead-packed reactors. Results from experiments conducted with static cell suspensions shed light on the strong control of cell settling that drives erratic, poorly reproducible and difficult to interpret SIP signals. However, incubating cells in bead packed reactors yielded reproducible trends in σ″, with strong (3 – 10 mrad) signals that followed the expected cell growth behaviour. Relating σ″ to measured cell density and metabolic activity (using ATP) highlighted the strongly linked contribution of both activity and cell density and SIP. Here, we report a lower frequency polarization peak between 0.01 and 0.1 Hz in the bead reactors, which we attribute to the polarization of cell colonies in the densely packed reactors. In summary, our findings shed light on the direct contribution of cells and their activity to polarization, in the absence of cell-sediment interactions and provide a novel approach for studying cell polarization in static incubation in a porous environment.

SummaryWe present DASPack, a high-performance, open-source compression tool specifically designed for distributed acoustic sensing (DAS) data. As DAS becomes a key technology for real-time, high-density, and long-range monitoring in fields such as geophysics, infrastructure surveillance, and environmental sensing, the volume of collected data is rapidly increasing. Large-scale DAS deployments already generate hundreds of terabytes and are expected to increase in the coming years, making long-term storage a major challenge. Despite this urgent need, few compression methods have proven to be both practical and scalable in real-world scenarios. DASPack is a fully operational solution that consistently outperforms existing techniques for DAS data. It enables both controlled lossy and lossless compression by allowing users to choose the maximum absolute difference per datum between the original and compressed data. The compression pipeline combines wavelet transforms, linear predictive coding, and entropy coding to optimise efficiency. Our method achieves up to 3 × file size reductions for strain and strain rate data in lossless mode across diverse datasets. In lossy mode, compression improves to 6 × with near-perfect signal fidelity, and up to 10 × is reached with acceptable signal degradation. It delivers fast throughput (100–200 MB s−1 using a single-thread and up to 750 MB s−1 using 8-threads), enabling real-time deployment even under high data rates. We validated its performance on 15 datasets from a variety of acquisition environments, demonstrating its speed, robustness, and broad applicability. DASPack provides a practical foundation for long-term, sustainable DAS data management in large-scale monitoring networks.

SummaryInitial stress exerts a crucial impact on the elastic properties and thus the wave reflection in the layered media. However, the stress effect on wave reflection characteristics in such media remain insufficiently understood. To address this issue, we develop a composite matrix reflectivity method incorporating initial overburden stress (CMRMS) by means of acoustoelasticity theory, enabling accurate modeling of seismic wave propagation in stressed layered media. The proposed method can better simulate multiple reflections, converted waves and transmission loss of seismic waves in layered media, compared to the classic stress-dependent reflection coefficient equation for a single interface. Moreover, our method can degenerate into the existing methods in the cases of no initial stress and single interface, which verifies its correctness. We further extended the CMRMS to elastic and viscoelastic non-welded interfaces using the linear-slip theory and standard linear solid model, respectively. The extended method is used to investigate the impacts of non-welded interface compliance, overburden stress, fluid viscosity and frequency on seismic reflection characteristics within layered model. It is shown that the stress effect magnitude on interface reflection significantly depends on the interface depth, due to cumulative transmission losses from overlying layers. Moreover, increasing either the compliance or the number of overlying non-welded interface significantly reduces the reflection amplitude at deeper interface. Our results show the potential of the proposed composite matrix reflectivity method to consider the joint effects of initial stress, multiple waves and transmission loss in both forward modelling and inverse applications.

AbstractThe low-velocity layer confined by surrounding rocks deep in the subsurface acts as a seismic waveguide. The compressional (P-) and shear (S-) waves propagate in the waveguide are reflected on the top and bottom interface, constructively interfered with to formulate the deep-guided wave. Deep-guided wave has high-frequency contents and notable dispersive features. The dispersion represents the kinematics information and can be used to image the low-velocity structures. The Earth media not only shows elasticity but also attenuates seismic waves. This article presents a theoretical study of deep-guided wave propagation and dispersion analysis in viscoelastic media. We utilize the Thomson-Haskell propagator matrix method to theoretically calculate the phase velocity dispersion and attenuation curves of the deep-guided wave in viscoelastic media. We apply the staggered-grid finite-difference scheme to numerically simulate the elastic wavefield propagation in the shale layer for validation. We have conducted a sensitivity analysis of the dispersion and attenuation curves of the deep-guided wave with respect to different media parameters. The theoretical calculation of dispersion and attenuation curves of deep-guided wave opens the doors for the simultaneous inversion of S-wave velocity and quality factor of the low-velocity layer in the future. Deep-guided waves hold the potential for high-resolution imaging of hydrocarbon reservoirs, geothermal reservoirs, coal seams, saline aquifers, and fault zones.

SummaryA strong earthquake sequence in Storfjorden, south of Svalbard, was initiated by an Mw 6.1 event on 21 February 2008. Earthquake distribution and fault plane solutions indicate that seismic activity is controlled by unmapped NE-SW striking oblique-normal faults, contrasting with the major N-S oriented faults mapped onshore Svalbard. We present a geophysical model derived from an ocean bottom seismometer profile crossing the seismogenic zone to identify structures in the crust and uppermost mantle that potentially control the earthquake source mechanism. Travel-time forward modeling using raytracing, combined with travel-time tomography and gravity-magnetic modeling, reveal distinct crustal domains across the earthquake region. Crystalline crustal P-wave velocities range from 6.1 km/s to 6.7 km/s at the Moho depth in the eastern section. The western profile section exhibits a higher Vp velocity lower crust (6.6–7.0 km/s) with Vp/Vs ratios of 1.75–1.8 and high density (∼3100 kg/m³). Basement depth reaches 8 km in the west, forming a sedimentary basin, and shallows eastward. The Moho remains relatively flat at 29-32 km depth throughout the profile. The N-S oriented Caledonian suture, identified from deep seismic and potential field data, traverses the Storfjorden earthquake zone. The lithological contacts within the suture zone, inferred from the new OBS data, may facilitate seismic failure oblique to the N-S oriented structure, following the regional stress field.

SummaryWe developed a new amplitude correction method for receiver function imaging to analyze velocity contrasts along dipping interfaces. Because receiver function imaging typically assumes a horizontally layered structure, corrections are needed for amplitude and polarity variations of P-to-S converted phases when analyzing dipping interfaces. However, previous studies have not adequately addressed these effects, and improved receiver function analysis is required to better delineate dipping structures, such as subducting plate surfaces and the oceanic Moho. Therefore, we propose formulae that quantify converted S-wave amplitude variations between horizontal and dipping interfaces. This relationship is expressed as a function of the back azimuth, the ray parameter of an incident P wave, and the dip angle and dip direction of a dipping interface, and in this study, the geometry of the dipping interface (dip angle and dip direction) is assumed. We applied these formulae to receiver function imaging using synthetic and observed data and confirmed that the amplitude of seismic discontinuities was successfully reproduced. This method enables the use of numerous receiver functions regardless of the back azimuths of incident P waves, thereby providing more detailed amplitude estimations for dipping interfaces.

SummaryTo better understand the mechanics of injection-induced seismicity, we developed a two-dimensional numerical code to simulate both seismic and aseismic slip on non-planar faults and fault networks driven by fluid diffusion along permeable faults, in an impervious host rock. Our approach integrates a boundary element method to model fault slip governed by rate-and-state friction with a finite-volume method to simulate fluid diffusion along fault networks. We demonstrate the capabilities of the method with two illustrative examples: (1) fluid injection inducing slow slip on a primary rough, rate-strengthening fault, which subsequently triggers microseismicity on nearby secondary, smaller faults, and (2) fluid injection on a single fault in a network of intersecting faults, leading to fluid diffusion and reactivation of slip throughout the network. This work highlights the importance of distinguishing between mechanical and hydrological processes in the analysis of induced seismicity, providing a powerful tool for improving our understanding of fault behavior in response to fluid injection, in particular when a network of geometrically complex faults is involved.

SUMMARYIn mineral exploration, induced polarization and self-potential are two broadly used active and passive geophysical methods, respectively. In the case of ore bodies, both methods are associated with charge distributions associated with a secondary electrical field (induced polarization) and a source current density (self-potential). Both the chargeability and volumetric source current density distributions bring information regarding the shape of ore bodies. Therefore the joint inversion of these datasets is expected to better tomograms of ore bodies. A joint inversion approach is developed to combine both methods. The objective function to minimize includes two independent components plus a cross-gradient joint function. The use of the cross-gradient is justified from the underlying physics of the two geophysical problems at play. The structure of the cost function is tailored to overcome some problems like convergence and parameter determination in the inverse process. Two synthetic tests and a laboratory experiment are used to benchmark the proposed algorithm. We demonstrate that the joint inversion algorithm performs better than the localizations obtained from independent inversion approaches. To refine the interpretation of the shape of ores, we introduce an ore presence index using the chargeability and source current density resulting from the joint inversion algorithm. The K-Medoids clustering algorithm is used to automatically categorize the calculated ore presence index into different clusters. The cluster with larger values successfully identifies the ore bodies associated with strong chargeability and/or volumetric source current density.

SummaryInversion of geomagnetic anomaly data poses an ill-posed problem, and extremal models such as equivalent source layers or point-source distributions can explain observations to the same degree as volumetric magnetisation distributions. However, the spectral characteristics of magnetic anomalies provide fundamental constraints for magnetic source-depth estimation. Specifically, the maximum detectable depth of crustal magnetic sources is dictated by the longest wavelengths present in the field, which correspond to the low-wavenumber bands of the spectrum. This relationship is often analysed through the log power spectrum versus wave number plot, using the slopes of the linear segment for depth estimation. Methods aiming at reconstructing the depth to the bottom of magnetisation from spectral field characteristics are commonly referred to as spectral methods. However, these methods are based on assumptions about the statistical properties of the source distribution and are prone to misinterpretations. Here, we apply sparsity-constrained 3D inversion of magnetic data using an elastic net regularisation to recover the susceptibility distribution and the bottom of magnetisation. We claim that the elastic net (ℓ2ℓ1 norm) regularisation, when properly tuned to balance the solution’s smoothness with sparsity, stabilises the inversion, avoiding extremal magnetisation distributions and generating a geologically plausible source depth distribution that is consistent with the expected source distribution. The ℓ1 norm brings sparsity and high resolution, while the ℓ2 norm brings inversion stability and structural continuity to the final model. From the recovered 3D elastic net sparse inversion model, we extract the depths of all the deepest non-zero susceptibility values and suggest this to be an alternative estimate to the base of magnetisation. Moreover, we suggest that the resulting 3D model has a value in itself and may aid geological interpretation.

SummaryGaining insight into the centennial evolution of the geomagnetic field over the past 2000 years requires the acquisition of reliable palaeomagnetic data from the study of well-dated archaeological materials or rocks. However, despite previous efforts, palaeointensity data from regions south of 30°S are still underrepresented, potentially limiting the accuracy of global geomagnetic field models and their applications. In addition, a comprehensive understanding of the geomagnetic field evolution in South America is particularly relevant, as the recent geomagnetic secular variation has been mainly characterised by the significant growth of the South Atlantic Anomaly over the past three centuries. The evolution of this low-intensity region, currently centered over central South America, is well understood in detail only during the last few centuries, thanks to the availability of direct measurements. For both the geomagnetic and palaeomagnetic communities, understanding its evolution prior to this period remains a challenge. This study presents new palaeointensity estimates from San Juan Province, central western Argentina, based on the analysis of 23 pottery samples dated between the 3rd and 17th centuries CE using radiocarbon and archaeological constraints. We employed the Thellier-Thellier method, incorporating partial thermoremanent magnetisation (pTRM) checks, TRM anisotropy corrections, and cooling rate adjustments, and obtained 11 mean palaeointensity values of good technical quality for central South America. The results are consistent with the limited number of previously reported high-quality palaeointensity data within an area 900 km in radius centered on San Juan, all showing intensity values ranging from approximately 40 to 55 μT. The new data, combined with these previously published high-quality intensities, do not show anomalously low values in intensity in the region between 200 and 1750 CE, suggesting no significant impact of the South Atlantic Anomaly in the region before the past three centuries. Furthermore, the findings suggest the presence of rapid multidecadal variations between 800 and 1100 CE, a behaviour also observed in other regions worldwide, which may point to a global or dipolar origin for these variations. By enhancing the dataset for this latitude range, this work provides new constraints on the geomagnetic field’s past behaviour south of 30°S over South America and contributes to improving future global geomagnetic reconstructions.