Geophysical Journal International

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.

SummaryFluid injection into the subsurface can trigger moderate-magnitude earthquakes days to months after shut-in, complicating hazard assessment. To investigate the governing mechanics, we implemented a fully coupled hydro–mechanical model that couples Darcy flow, poro-viscoelastic deformation and rate-and-state fault friction on a planar fault, allowing two-way feedbacks between pore pressure, volumetric strain and fault slip and the simulation of both aseismic and seismic transients. Compared with decoupled or one-way approaches, the fully coupled formulation generally yields longer post-injection delays, owing to poroelastic stress contributions and a more realistic evolution of volumetric strain. After shut-in, a slow poroelastic redistribution of volumetric compression broadens and migrates along the fault, constructively overlapping regions of elevated shear and reduced effective normal stress. This causes the nucleation of a delayed rupture away from the well, indicating that the point of peak instantaneous pressure does not necessarily coincide with the location of maximum coseismic slip. By scanning permeability and injection rate we construct an empirical injection-rate (IR)–permeability (k) phase diagram that delineates regimes of immediate, delayed and no induced seismicity; this diagram is offered as a conceptual, physics-informed screening tool that requires site-specific calibration. Our results indicate that two-way hydro-mechanical coupling and fault slip evolution should be considered when assessing post-injection seismic hazard and in the design of spatially distributed monitoring.

SummaryTraditional seafloor mapping relies on shipborne soundings which have limited spatial coverage. The Surface Water and Ocean Topography (SWOT) wide-swath altimetry satellite holds the potential for predicting more detailed seafloor topography. In this study, we integrate SWOT gravity data with single-beam shipborne depths to construct seafloor topography models in the Northwestern Pacific using the deep neural network (DNN) method. Compared to shipborne depth checkpoints, the root mean square (RMS) error of the differences between topography model predicted by DNN method and shipborne depths is approximately 97.5 m, improving by 19.5% and 9.9% compared to the gravity-geologic (GGM) method and the Smith and Sandwell (SAS) method respectively. Compared to traditional data, the integration of SWOT gravity data universally enhances prediction accuracy. Furthermore, the DNN method effectively demonstrates superior capability in balancing the characterization of overall structures with the retention of authentic topography features, which we demonstrated in the Mariana region of the NW Pacific Ocean. However, limited by spatial heterogeneity and physical mechanisms, accurate prediction of such complex, fine-scale topography using gravity data remains a significant challenge.

SummarySpatial autocorrelation (SPAC), the azimuthal average of the normalized cross-correlation between equidistant station pairs deployed in a 2-D array, is widely used to image the subsurface structure. However, the rigorous estimate of subsurface structure and its uncertainties as a function of depth using SPAC data is challenging due to the nonlinear relation between the SPAC data and Earth structure as well as the trade-off between depth and velocity. Additionally, data noise is strongly correlated due to data processing (e.g., filtering, stacking from multiple time segments and azimuthal averaging). Most studies do not account for the correlated noise and fix the ratio of compressional-wave velocity (${V}_P$) to shear-wave velocity (${V}_s$) (i.e., ${V}_P$/${V}_s$ ratio) and the number of layers, both of which are typically unknown. To address these challenges, we develop a hierarchical trans-dimensional Bayesian inversion of fundamental mode of SPAC data that properly accounts for the correlated data noise, samples the ${V}_P$/${V}_s$ ratio, and relaxes the number of layers (i.e., model parameterization) to be unknown in the inversion. We further examine the limitation of using only fundamental modes in the inversion. Our synthetic experiments show that the inversion recovers an incorrect model unless we sample the correlated noise and ${V}_P$/${V}_s$ ratio in the inversion. The inversion is then applied to SPAC data acquired at 19 sites across the Tacoma basin in Washington State to characterize the ${V}_s$ and the time-averaged ${V}_s$ over 30-m depth ($V{s}_{30}$). Our results show that the $V{s}_{30}\ $varies from ∼200 m/s to 800 m/s. The $V{s}_{30}$ within the basin is higher in the middle and lower on the east and west sides. We find that these $V{s}_{30}$ values vary with geologic unit. The uncertainties for $V{s}_{30}$ are within 20 m/s in average except for the most eastern site TB28. Additionally, the uncertainties are greater for deeper depths beneath most of the sites as the sensitivity decreases as a function of depth. The ${V}_s$ structure as a function of depth is also complex beneath some sites, possibly because the SPAC curves are affected by higher order Rayleigh modes that are not considered in the inversion. To better constrain the deeper ${V}_s\ $structure, $V{s}_{30}$, and/or other average measures of ${V}_s$ over depth, additional constraints from complementary data, such as ellipticity or geologic data are needed. Moreover, our synthetic experiments show that higher order modes can have significant effect in the inversion results, particularly when there is a low-velocity layer.