Geophysical Journal International

SummaryElectrical resistivity tomography (ERT) is used to infer the subsurface resistivity structure. ERT requires solving a nonlinear inverse problem that is often approximated as linear to reduce computational time. However, the approximation requires assumptions that cause limitations for the data analysis. Most of the computational time is due to the forward problem that requires solving the Poisson equation. Recently, similar forward problems have been shown to be replaceable with a surrogate model of lower computational cost. We present a geoelectric surrogate based on Fourier Neural Operators (FNO) and demonstrate a successful application in nonlinear inversion. The standard deviation of the surrogate prediction error for unseen samples is <5%. Furthermore, the surrogate reduces computational time by over three orders of magnitude per realization, enabling ERT for previously intractable settings. We apply the surrogate in Markov chain Monte Carlo (MCMC) inversion of simulated data. The results resolve sharp resistivity changes with plausible uncertainties.

SummaryThe Geocenter Motion (GCM) time series captures periodic variations arising from diverse Earth system changes. This study pioneers the use of Successive Variational Mode Decomposition (SVMD) in GCM research, enabling the precise extraction and analysis of these meaningful geophysical signals. SVMD outperformed Singular Spectrum Analysis (SSA) by effectively isolating signals and minimizing interference from components with similar variance contributions. However, a high maximum penalty factor in SVMD may lead to noise-dominated Intrinsic Mode Functions (IMFs). To overcome this limitation, we propose an extraction criterion that utilizes the standard deviation of the correlation coefficient and mean kurtosis as thresholds. Validations with simulations and the real GCM time series demonstrate its superiority over traditional single- and dual-threshold criteria, effectively retaining valuable information while excluding most noise-dominated IMFs. This improved approach is further employed to explore the geophysical driving factors of key periodic variations in the GCM time series, focusing on the annual, semi-annual, 10.5-year, 451-day, ∼160-day, and ∼120-day periods. Multi-source GCM analyses combined with the fingerprint method reveal distinct contributions from the Antarctic and Greenland ice sheets, terrestrial water storage, continental glaciers, and atmosphere-ocean interactions to different periodic signals. This study provides a robust methodology for decomposing GCM and attributing its variations to underlying Earth system changes, advancing our understanding and interpretation of global mass redistribution.

SummaryFull waveform inversion (FWI) is a popular method for subsurface parameter estimation. Despite its effectiveness in building high-resolution velocity models, the quality of the inversion result is significantly dependent on a fairly accurate, smooth initial model, which is often challenging to build. To weaken the influence related to the inaccurate initial model, we propose a deep learning (DL) matching-based FWI framework, namely DLM-FWI, where multiple convolution neural networks (CNNs) are used to construct an adaptive matching filter to better pinpoint the discrepancies between the synthetic and observed data. With the help of the CNN-based matching filter, the synthetic data will be regularized first, leading to intermediate data, and the model update will be conducted by minimizing the misfit between the intermediate and the observed data for improved data-fitting. More importantly, we integrate the whole inversion process into an automatic differentiation (AD) framework, simplifying the implementation of classic FWI. We apply the proposed DLM-FWI method to both synthetic and field datasets to validate its effectiveness. The results demonstrate that compared with classic FWI, DLM-FWI performs better in subsurface model reconstruction when the initial model is far from the global minimum.

SummaryUnderstanding subsurface fluid distribution in volcanic reservoirs is critical for geothermal energy development, critical mineral exploration, and forecasting eruptions. Here, we use travel-time tomography to image the seismic velocity structure beneath Aluto volcano, the first pilot geothermal project in Ethiopia, located in the Main Ethiopian Rift. Using seismic data recorded from January 2012 to January 2014, we invert for the 3D P-wave (Vp), S-wave (Vs), and Vp/Vs ratio. To reduce the non-uniqueness in interpretation, we also compare our results with previously published work on attenuation tomography and magnetotelluric images. Elevated Vp/Vs ratios (at 0 km below sea level (bsl)) around productive geothermal wells suggest high fluid content and/or elevated temperature. Vp/Vs values above 1.8 are observed along the caldera rims and hydrothermal vents, indicating fault and fracture systems as primary fluid conduits. High Vp/Vs below 6 km bsl likely reflects high-temperature areas or the presence of partial melt. In contrast, low Vp/Vs (<1.5), low Vp, and average to high Vs beneath the caldera at around 5 km bsl is interpreted as a crystallised body with over-pressurised gas volume formed during phase separation and transported upward through fractures and fault systems, accumulating at shallower levels. These findings highlight fluid pathways through the caldera rims and faults, with volatile-rich partial melt at greater depth beneath the caldera centre. Travel-time tomography thus offers a valuable constraints on subsurface fluid distribution and is valuable tool in geothermal exploration.

SummaryThis study presents a workflow to monitor spatiotemporal variations of the secondary microseisms using multi-array analysis. We employ ambient-noise cross-correlation beamforming (CC beamforming) across three dense seismic networks with different instrument responses: ANTICS in Albania (nodal-geophone and broadband), Hi-net in Japan (short-period), and SCSN (broadband) in Southern California. Independent of their instrumentation, these networks enable us to track the spatial and temporal evolution of secondary microseism sources in the northern Hemisphere from autumn 2022 to spring 2023. The workflow involves continuous data preprocessing for different instrumented sensors, ambient-noise cross-correlation, beamforming, and beam-power back-projection into a global map. We also propose sliding-window raw-data beamforming (RA beamforming) for the continuous broadband data in this workflow to record the absolute amplitudes of secondary microseisms recorded by ANTICS. Joint CC beamforming analysis across the three different networks improves the resolution of ambient-noise source localization and displays high consistency with the equivalent vertical force at the ocean floor. The results indicate that secondary microseism sources in the northern Hemisphere are predominantly driven by winter storms in the northern Atlantic and northern Pacific. The relative and absolute amplitudes of the beam-power for the northern Atlantic are also extracted from CC beamforming based on geophone sensors and RA beamforming based on broadband instruments from ANTICS, respectively. Both approaches provide robust estimates of microseism strength in the northern Atlantic, with CC beamforming displaying a higher correlation with the modeled ocean floor equivalent forces. This study confirms the feasibility of using cost-effective nodal seismic arrays for detailed monitoring of secondary microseisms and highlights the potential for integrating multi-array seismic data with oceanographic models for an improved understanding of seismic noise generation and propagation.

SummaryMoho topography model of subduction zones is an important component of deep tectonics and an important basis for verifying geodynamic processes. As one of the main factors affecting the accuracy of Moho topography model inverted by gravity method, the selection of inversion parameters suffers from the effect of nonlinear terms, which need to be reduced by constraints. Therefore, we applied the differential evolutionary algorithm to compute the inversion parameters and obtained a refined Moho topographic model of North Manila subduction on this basis. Synthetic tests show that the differential evolution algorithm can effectively mitigate the impact of nonlinear terms. With or without noise, the differential evolution algorithm is effective in finding better inversion parameters compared to the linear regression method. Particularly in Moho density contrast, the average value obtained from multiple runs of differential evolution algorithm still achieved a 54.4% improvement in accuracy. In practical application, the comparison results show that the RMS of the difference between this paper’s model and all seismic control points is 2.37 km with an improvement of at least 35.1%, which proves that this paper’s method is reliable. Furthermore, we examined the impact of various parameters on the method to validate its robustness.

SummaryThe far-field impact of Tibetan Plateau (TP) expansion on cratonic blocks remains enigmatic. We address this for the Ordos Block (OB) by constructing a high-resolution 3-D shear-wave velocity model using ambient noise tomography from an unprecedented dense seismic array (461 stations). Our model reveals: (1) NE-trending high-velocity anomalies at 10–25 km depth correlating with crustal magnetic signatures, providing seismic evidence for late Archean amalgamation of micro-blocks (Jining, Ordos, Xuchang, Xuhuai); (2) TP-induced reactivation manifesting as southwestern OB crustal thickening (50 km) with a high-velocity lower-crustal layer (≥4.0 km/s; 100 km wide, 10 km thick), attributed to TP lower-crustal underthrusting beyond the plateau margin (35.5°–37.5°N), facilitating >200 km strain transfer into the OB interior; (3) Incipient rifting dynamics in the Daihai Rift, where upper-crustal high-Vs (preserved rigidity) overlies mid-lower crust/uppermost mantle low-Vs anomalies (mantle-sourced thermal modification), indicating early-stage rifting driven by combined Pacific plate retreat and TP far-field stresses; (4) Craton-wide segmentation across a fundamental 37.5°N lithospheric boundary demarcating mantle upwelling/crust-mantle interaction (north) from passive TP push-dominated deformation (south). These findings redefine the OB as a strain-partitioned system, where lithospheric heritage controls differential response to plate-boundary forces.

SummaryThe Tanlu fault zone is an NNE-SSW oriented, large and deep strike-slip fault system running through eastern China. To investigate seismotectonics in and around the Tanlu fault zone, we adopt the LOC-FLOW approach to build a high-precision earthquake catalog. Our seismic data were recorded at 120 broadband TanluArray temporary stations with a sampling rate of 40 Hz and 76 broadband permanent stations with a sampling rate of 10 Hz from July 2019 to March 2023. We first conduct a series of experiments around the Luxi uplift and find that a higher sampling rate and a denser array of stations can significantly enhance the earthquake detection ability. Then we utilize both the temporary and permanent stations to conduct earthquake detection and location for the entire study area. As a result, 9648 earthquakes are detected and located in the REAL catalog, 6543 earthquakes in the HypoInverse catalog, and 5619 earthquakes in the HypoDD catalog, representing an increase of 20%, 22%, and 22%, respectively, as compared with the cases when only the TanluArray stations are used. Our location results show that earthquakes are mainly distributed in the Tanlu fault zone and active faults in relevant tectonic units. We collect 322 focal mechanism solutions (M > 2.0) of previous results from 2000 to 2020 to invert the stress field of the whole study region. The results show that the maximum principal stress axis of the whole study area is in the NEE-SWW direction, except that the Huoshan region is in the E-W direction. Along the Tanlu fault zone, the highest seismicity occurs in the Suqian-Weifang segment, and the Suqian seismic gap may be due to the aseismic slip along the fault planes of the Tanlu fault zone.

SummaryIn recent years, physics-informed neural networks (PINNs) have emerged as a powerful tool for seismic traveltime modeling and tomography. However, conventional PINNs do not consider applicable physical priors or quantify the uncertainty of the inverse problem, which is critical for reliable geological interpretation with topographical complexity. Thus, we propose a comprehensive PINN-based framework designed to tackle the critical challenge of inverting for the velocity in models with irregular topography, while also quantifying the inherent uncertainty. Leveraging automatic differentiation, our mesh-free approach directly accommodates complex surfaces without the need for specialized grids. To enhance inversion accuracy and physical consistency, we uniquely incorporate additional physical priors, namely well-log velocity profile and the principle of reciprocity. Furthermore, to address the non-uniqueness of the inverse problem, we integrate Monte Carlo (MC) Dropout to efficiently quantify model uncertainty without architectural modifications. Through 2D and 3D experiments on synthetic and real-world geological models, we demonstrate that our method accurately inverts for velocity structures with highly irregular topography. Results show that the inclusion of physical priors significantly improves model performance, while uncertainty quantification via MC dropout successfully highlights regions of higher uncertainty in the inverted velocity field, aligning with geological complexities in the velocities. This work establishes a robust and practical methodology for accurate and reliable seismic tomography in challenging geological settings.

SummaryShear-wave splitting measurement returns two parameters related to the fabric of the upper mantle: the orientation of the fast polarisation (fast direction), and a measure of the intensity and thickness of the fabric known as the split time. Spatial statistics of compiled splitting measurements indicate that the fast direction is spatially coherent, while the split time is not. We show, through modelling large numbers of noisy measurements, that single-earthquake splitting measurements exhibit a prominent upward bias in split time, the degree of which depends on specifics of the measurement process. Averaging single-event splitting parameters over many measurements does not mitigate this bias; however, stacking of error surfaces from individual measurements does, given sufficient back-azimuthal coverage, while also greatly reducing scatter. Published splitting results use a mix of these two averaging techniques, and this inconsistent bias between studies is likely responsible for the lack of spatial coherence in compiled split-time measurements. We demonstrate this in real data by examining a data set from Alberta, Canada and surrounding areas, for which a recent study published parameter-averaged results. By examining a comparable data set using error-surface stacking, we are able to greatly increase the coherence of the split times while obtaining highly similar fast directions. Our coherent split times are mapped to reveal a zone of strong splitting beneath the active Cordillera, and three zones of moderate to low split time within the cratonic lithosphere.

SummaryUnderstanding the spatial variability of ice content in frozen ground is key to design adequate measures to manage different ecosystems in frame of climate change. To-date investigations in frozen ground require the analysis of borehole data or the collection of multiple geophysical data. Here we propose the use of spectral induced polarization (SIP) as a technique that provides in quasi real-time information about changes in ice content in the subsurface. We demonstrate that exploring the frequency dependence in electrical conductivity and polarization (capacitive) properties in the frequency range between 0.1 and 75 Hz provides direct information about their relative variations in ice content. Our study is based on measurements conducted at nine representative permafrost sites in the European Alps with varying landforms and ice contents, including a pure ice and an unfrozen reference. We use the phase frequency effect (ϕFE) parameter as a parameter describing both the amplitude of the polarization and its frequency dependence to compare the response associated to the different sites. Our results show the lowest ϕFE in sites with low ice content, while increases in this parameter are associated with higher ice content. We evaluate the correlation between SIP parameters and validation ground ice contents for all sites and observe a clear correlation between ϕFE and volumetric ice content. The ϕFE results exhibit distinct landform-specific patterns, with the highest values found in rock glaciers, intermediate values in frozen talus slopes, lower values in bedrock permafrost, and the lowest in unfrozen talus slopes, reducing interpretation ambiguities in electrical resistivity results for assessing ice content.

SummaryTime-domain induced polarization (TDIP) data carry spectral information that can be used for petrophysical interpretation. At the same time, TDIP data can be collected in the field more efficiently than frequency-domain induced polarization (FDIP) data, thanks to the use of square-wave signals. However, TDIP field data are prone to noise, particularly strong near industrial installations and urban areas, above conductive media and in cases where little current is injected. The integral chargeability is a useful parameter to smoothen out the signal but it precludes any spectral interpretation. Debye decomposition (DD) is recognized as one of the best methods for spectral interpretation but the extracted parameters are particularly affected by data noise. More generally, processing TDIP data before further analysis, such as inversion or spectral analysis, is usually necessary for any quantitative interpretation. We propose here an automatic processing algorithm, based on the Kohlrausch-Williams-Watts (KWW) function, which is very close to the Havriliak-Negami model in frequency-domain, that fulfills this need. The processing is completed by an empirical handling of early-time electromagnetic coupling effects to improve the overall performance. The resulting procedure, tested and validated on three datasets that cover a large range of contexts, electrode configurations and acquisition settings, is available as open-source MATLAB scripts. The proposed approach is especially useful for further extracting spectral information from TDIP data through DD. Thanks to the theoretical framework offered by the KWW function, the behavior of the integral chargeability could be investigated in a systematic manner, using both synthetic and field TDIP data. Recommendations could be formulated on how to make use of the spectral information, while keeping the automatic processing transparent and accessible to unexperienced users. This work advances the use of TDIP in the field of environmental geophysics.

SummaryJoint petrophysical inversion is a powerful technique for using multiple geophysical modalities to estimate petrophysical or geotechnical parameters of the subsurface. A precise knowledge of the petrophysical laws for the full model domain is imperative to enable petrophysical coupling. In this work, we investigate the effect of partially invalid petrophysical laws on the inversion of a synthetic data set, using electrical resistivity tomography (ERT) and seismic traveltime data to image a CO2 plume in a Carbon Capture and Storage (CCS) setup. We consider a model consisting of a reservoir and a caprock in which only the reservoir can be described by a petrophysical law. We first apply a conventional (joint) petrophysical inversion (JPI) and show that the use of wrong petrophysical laws leads to systemic artefacts within the parts of the model in which the petrophysical relations are invalid. We then present a new hybrid partially petrophysically coupled joint inversion (P-JPI) approach that combines petrophysical coupling for regions with valid petrophysical laws, and structural coupling, whenever no reliable petrophysical laws are available. The P-JPI approach outperforms tomography based on the individual ERT or seismic data set, as well as joint structural inversion (JSI) based on the cross-gradient functional. The partially petrophysically coupled joint inversion thus enables petrophysical coupling and provides a unique, quantitatively interpretable saturation model for the CO2-plume. We further show that it is possible to detect zones with incorrect petrophysical relations by analysing the difference of the model updates based on the stand-alone data sets. Finally, we combine the detection of zones of incorrect petrophysical laws with the P-JPI to derive an inversion scheme that is independent of prior knowledge of the validity of petrophysical laws. Our novel methods facilitate direct estimation of the petrophysical subsurface parameters from multiple geophysical measurements if petrophysical relations are only available for parts of the model domain and provide means to quantify the spatial extent of regions where the petrophysical relations are valid.

SummaryThe rapid and accurate estimation of strong ground motion is essential for seismic hazard assessment and near-real-time disaster response. Although empirical ground motion models enable fast intensity predictions, they simplify the underlying physics and exhibit large uncertainties. Conversely, physics-based simulations — while capable of more accurately predicting ground shaking — are computationally expensive, making them impractical for large-scale hazard assessments and real-time event response. To overcome these limitations, we introduce a novel two-step machine learning framework that predicts peak ground velocity (PGV) for arbitrary double-couple sources positioned anywhere within a given medium, combining the accuracy of physics-based models with near-instantaneous inference. In the first step, an ensemble of XGBoost predictors, trained on a reciprocal Green’s function database, generates a sparsely sampled PGV map for any input source. In the second step, we refine this map into a continuous spatial prediction. By leveraging Green’s function reciprocity, our approach reduces the required number of simulations in training, lowering both computational cost and storage demands. Our framework provides spatially continuous PGV predictions and inherently accounts for complex 3D geological and topographic effects. It can deliver results within seconds while maintaining accuracy up to the highest frequency captured in the physics-based simulations. This makes PGVnet ideal for applications such as rapid earthquake disaster response, as well as large-scale probabilistic seismic hazard analyses and multi-hazard digital ecosystems. Validated in the geologically complex San Francisco Bay Area, our approach generates PGV maps consistent with physics-based simulations, offering an effective balance between computational speed and accuracy.

SummaryRevil asserted in the comment that none of Macnae’s paper is novel. The paper however introduced a novel method of chargeability prediction as the fraction of pores blocked by metallic particles rather than the prediction using limiting Maxwell effective medium estimates based on volume fractions. This reply presents an analysis of two cases where the novel chargeability prediction proves to be significantly better than earlier methods when applied to laboratory pyrite-clay mixtures and to published petrophysical data from a Co-Cu disseminated mineral deposit. Another item of novelty in the paper was emphasis that the resistivity and conductivity time-constants can differ by orders of magnitude for the high chargeabilities of economic sulphide deposits, a fact not commonly recognised in the literature. Revil asserts in the comment that the term Induced Polarization (IP) should explicitly exclude dielectric effects and analogies as discussed in the paper, an assertion inconsistent with the early literature on IP and the fact that both diffusive and dielectric effects can affect low-frequency data. The reply provides detail on the physical nature of equivalent circuits, and how they can mechanistically model the IP phenomenon and the topology of conductive paths in materials, an issue that I did not emphasize in the paper as I thought it would have been obvious to most readers.

SummaryVolcanic and geothermal areas usually feature highly rugged and variable topography, with both steep and relatively flat areas. Potential field geophysical data being very sensitive to topographic changes, it can be challenging to accurately simulate data in such regions. Traditional modeling approaches for potential field data typically involve discretizing the physical space using rectangular prisms, or alternatively using slightly more sophisticated geometries such as polyhedrons. However, in both cases, these discretizations may lead to an under- or over- estimation of topographic effect on simulated data, unless a very fine discretization around observation points is used, consequently leading to an increase in the computational burden. To address this issue, we introduce new methodological strategies to simulate potential field data in regions with rugged topography. We propose the use of a numerical integration scheme over deformable hexahedral elements that conform to topographic variations, to numerically evaluate the integral equations governing gravity and magnetic anomalies. Density and magnetization properties are defined at grid nodes and interpolated within the elements using polynomial functions. The numerical integration relies on the use of a Gaussian quadrature, combined with a local refinement of the mesh designed to simultaneously increase the precision of the numerical integrals by reducing the effect of singularities and incorporating a more precise description of the topography at the vicinity of measurement points. This refinement is coupled with a discontinuous optimization of the elevation of geometrical nodes lying at the surface of the model by minimizing the distance between the interpolated topographic surface and a point cloud provided by a high-resolution Digital Elevation Model (DEM). The predictions at each observation points are computed with different refined versions of the mesh in order to minimize computational cost. Physical properties within the refined elements are extrapolated from the nodes of the original non-refined mesh in order not to introduce additional unknown physical parameters for the later implementation of this code for the inversion of the potential field data. The approach is first validated by comparing gravity and magnetic predictions with the analytical solution for a simple rectangular prism. We further test the method with realistic simulations using the complex topography of the Krafla geothermal area in northern Iceland. Simulation results are compared against the forward gravity and magnetic responses of a high-resolution benchmark model (discretized with 2 m × 2 m rectangular prisms) computed using the open-source software Tomofast-x (Ogarko et al. 2024). Finally, we use the same high-resolution model to quantify typical errors associated with coarser rectangular prism discretizations (cell sizes of 25 m, 50 m, and 100 m). We demonstrate that these errors can be comparable in magnitude to real geophysical anomalies and are not reproduced when using our proposed numerical approach.

SummaryThe amplitude of the ${{P}_s}$ phase relative to the direct P wave (i.e., the ${{P}_s}/P$ ratio) provides valuable information about the contrasts across the crust-mantle boundary. Understanding how these amplitudes respond to variations in subsurface parameters improve interpretation of lithospheric structures. We examined the sensitivity of eight key parameters, including compressional and shear wave velocities, lower crustal and upper mantle densities, ray parameter, and Moho depth, to receiver function (RF) amplitudes and ${{P}_s}/P$ ratios. Using the synthetic RF (synRF) code of Ammon et al. (1990), we applied a Monte Carlo approach to generate randomized parameter sets, incorporated random noise, and tested both sharp and gradational Moho structures. The results show that lower crustal shear velocity has the strongest influence on all RF amplitudes and ${{P}_s}/P$ ratios, while lower crustal ${{V}_p}$, density, and upper mantle ${{V}_s}$ have moderate effects, and upper mantel ${{V}_p}$ and density show weaker sensitivity. In both sharp and gradational models, lower crustal and upper mantle shear velocities largely control the ${{P}_s}/P$ ratio. Theoretical ${{P}_s}/P$ ratios exhibit higher correlation with observed RFs than synthetic ones. Compared with the Shen & Ritzwoller (2016) model, our analysis yields ∼10% lower uncertainties in lower crustal ${{V}_s}$ and smaller uncertainties in lower crustal density derived from ${{P}_s}/P$ ratios, with consistent results in complex regions such as the northern Rocky Mountains. This study establishes the first quantitative framework linking ${{P}_s}/P$ ratio variability to lower crustal velocity and density while explicitly quantifying parameter sensitivity and uncertainties, clarifying how Moho sharpness and noise affect amplitude stability.

SummaryMacnae (2025) presented a physical interpretation of the Cole Cole Complex Conductivity model in the case of porous materials with sulphides. According to his paper, the Cole Cole parameters determined from such model can be easily interpreted in terms of underlying physics. His model is partly based on the electrochemical polarization model of Wong (1979) to explain the relationship between the chargeability and the volumetric content of sulfide. None of the statements made by Macnae (2025) are however novel. That said, we agree with Macnae (2025) that the Cole Cole complex resistivity relaxation time is quite useless in deciphering the underlying physics of the induced polarization problem.

SummaryEfficient Markov chain Monte Carlo (MCMC) sampling from posterior distributions remains a central challenge in Bayesian geophysical inversion. Recent developments in computational statistics and optimal transport suggest that MCMC efficiency can be improved by reparameterising the sampling problem – specifically, by learning an invertible mapping that recasts the target distribution onto a simpler reference distribution. Here, we introduce a Metropolis–Hastings framework that leverages transport maps parameterised by invertible neural networks. These maps are trained on preliminary MCMC samples from the target distribution and used to propose new samples in a fixed reference space, where proposal design is independent of the target’s structure. The proposed samples are transformed back to the target space via the inverse map, and accepted or rejected according to a modified Metropolis–Hastings criterion. As sampling proceeds, the transport maps are updated, yielding proposals increasingly well adapted to the shape of the target distribution. Across a suite of numerical tests – including a 2-D Rosenbrock distribution, a 3-D earthquake location problem, and Gaussian mixtures up to 16 dimensions – transport-map-driven samplers consistently outperform standard MCMC, reducing integrated autocorrelation times by factors of 2.5 to over 6 (or equivalently, yielding sample sets 2.5–6 times larger for the same number of forward evaluations). This improvement comes at the non-negligible cost of training one or more transport maps, which we quantify systematically. We also provide a quantitative criterion for weighing training cost against sampling speed-up. This shows that transport-map MCMC is advantageous whenever the forward problem is nontrivial, making it a promising approach for Bayesian sampling in geophysics and beyond.

SummaryEarthquake locations and catalogs from routine earthquake monitoring are typically based on manually reviewed arrival-time picks from classical, rule-based automatic pickers. High-performance, deep-learning (DL) pickers can replace this standard approach, rapidly delivering much larger and complete catalogs. A transition to routine monitoring based on DL picks requires that resulting catalogs include all or almost all events identified by current procedures with locations of the same or higher quality. Here we verify these requirements by comparing the performance of DL and manual picking for earthquake relocation and tomographic inversion. We form a reference catalog with a subset of INGV bulletin events and picks from the 2016 Central Italy sequence. This catalog is re-picked using the DL picker PhaseNet trained on the Northern California Earthquake Data Center dataset and on the INSTANCE Italian dataset. We use these three pick sets for high-precision, non-linear earthquake relocation and for 3D tomographic inversion and relocation. Relative to the high-precision relocations using routine picks, those using DL picks show improved organization and clustering, and, in a ground-truth test, smaller hypocenter separation for event pairs with more similar waveforms. The tomographic inversions show statistically better convergence and more organized relocations using the DL picks than with the routine picks. We conclude that DL based monitoring can rapidly produce more consistent picks and higher quality catalogs than standard procedures, while freeing analyst time for improved quality control, assessment, interpretation, and dissemination of information on seismic activity, especially during significant seismic sequences.