Geophysical Journal International

SummaryCharacterizing ore deposits or mining dumps in terms of mineral content and grain size remains a challenge. Since the 1950s the Induced Polarization (IP) method has been successfully applied in ore prospecting. However, reliably interpreting field survey data requires comprehensive laboratory studies to establish a link between the IP parameters from empirical or phenomenological models and the type and quantity of ore minerals. In this study, we use numerical electrical networks to replicate the complex electrical resistivity spectra observed in experiments on sand-pyrite-water mixtures. A network consists of a 3D assembly of resistors, representing the saturated pore space, and leaky capacitors simulating the electrical behaviour of ore minerals. A sophisticated fitting procedure enables the precise determination of resistor and capacitor parameters, ultimately leading to strong agreement between measured and synthetic IP spectra. The results obtained from the 3D network align well with the classical Pelton model, which is based on a simple equivalent circuit. Our findings indicate that the network's chargeability depends on the fraction of capacitors in the system (i.e. the number of capacitors divided by the number of capacitors and resistors), and that the Pelton time constant of the measured spectra is closely related to the resistor and capacitor parameters. We argue that a 3D approach offers a more realistic framework, paving the way for future studies on the effects of ore grain size distribution, and the spatial arrangement of ore grains.

SummaryThe most widely used method to derive global geomagnetic field models for historical and longer timescales has long been regularized non-linear least squares inversion. It is based on spherical harmonics for the spatial part and cubic B-splines for the temporal dynamics. Recently, different versions of Bayesian inversion have been applied for this purpose. Early literature on the traditionally used formalism states the inverse problem in a Bayesian setting and discusses uncertainty estimation via the posterior covariance, but this view was lost in subsequent studies in the geomagnetic community. Here we aim to provide both geomagnetic field modellers and users of such models with a comparative view of the methods to enable them to better evaluate strengths and weaknesses of different models. We first describe the connection between regularized least squares and Bayesian inversion in general form in a linear, one-dimensional setting. A fully Bayesian perspective allows interpreting the regularization term as a form of prior and offers new ways of comparing models from both approaches. We then discuss the particular case of geomagnetic field modelling. We find that in comparison to Bayesian modelling approaches the prior corresponding to the widely used regularization does not imply reasonable field properties and does not lead to meaningful uncertainty estimates.

AbstractThis study presents a refined interpretation of the Jeokjung-Chogye Basin (JCB), a confirmed meteorite impact structure in South Korea, by integrating high-resolution gravity data, microtremor measurements, and borehole information. A total of 1 700 gravity stations including 1 000 newly acquired in 2023, were used alongside horizontal-to-vertical spectral ratio (HVSR) analysis and well-log constraints to characterize subsurface structures. To isolate the impact-induced deformation from overlying sedimentary effects, gravity stripping was applied to remove the signal from post-impact ejecta deposits. The residual gravity field was analysed using dip-curvature mapping and Euler deconvolution, which revealed concentric ring structures with displaced centres. These asymmetries, corroborated by 3D forward gravity modelling using IGMAS+, suggest a northeast-to-southwest impact trajectory with an oblique incidence angle of approximately 45°, contrasting with earlier estimates of ~55° from east to west. The final 3D density model achieves a strong correlation with observed anomalies (R ≈ 0.95) and successfully resolves variations in the autochthonous and basement layers.

SummaryThe mantle transition zone (MTZ) plays an important role in the global material circulation, slab dynamics, and seismogenesis of deep earthquakes in subduction zones. Here we construct fine MTZ structures of both P and SH waves beneath the northeast Asia continental margin using improved grid-search waveform modelings, based on high-quality triplicated waveforms of three deep Kuril earthquakes recorded by the China Digital Seismograph Network (CDSN). We find a high-velocity anomaly (HVA, average 3.3 per cent δVp and 2.3 per cent δVs) with a thickness of 130-138 km in the lower MTZ. The HVA hosts a top interface with positive velocity contrasts (δVp: 4.3 per cent, δVs:3.2 per cent), while the 660-km discontinuity (660) shows reduced velocity contrasts (δVp: 3.6 per cent, δVs: 5.1 per cent) and negligible depressions of less than 10 km. The HVA we detected likely implies the thickened stagnant Pacific slab that may alter localized heat exchanges between the MTZ and lower mantle. The increased Vp/Vs ratio (∼1.85) indicates a water-rich state (∼0.42 wt per cent) inside the stagnant slab, evidencing the deep water transportation by the slab subduction. We infer that the interior localized dehydration of hydrous minerals within the stagnant slab may trigger the large outboard 1990 Mw 7.2 Sakhalin (Kuye) earthquake. Our results can provide more insight into slab dynamics and seismogenesis of deep earthquakes in northeast Asia.

SummaryWe used broadband ocean bottom seismometer (BBOBS) data from the RHUM-RUM experiment to derive the compliance function and estimate the shear velocity (Vs) structure of the subsurface at several sites beneath the Indian Ocean. The primary objective is to map the geological features of poorly explored marine regions, utilizing the compliance function, a measure of seafloor deformation in response to infragravity pressure signals at low frequencies (0.003 to 0.04 Hz). Compliance is the transfer function between vertical displacement and pressure, which is most sensitive to subsurface shear velocities. Our analytical process involves several data processing steps, including the removal of glitches, filtering out seismic events, minimizing tilt effects, calibrating pressure gauges, searching over the frequency and coherence domains to determine the optimal data window, and performing depth-velocity inversion using Monte Carlo method, specifically the Metropolis-Hastings algorithm. We present the ’ComPy’ software, which automates these processing steps for seafloor compliance analysis. The data, recorded over 13 months in 2012-2013 over a large region stretching from La Reunion Island to the Central Indian Ridge (CIR) and the South-West Indian Ridge (SWIR) (water depths of 3 to 5 km), confirm the stability of the compliance function over time. Depth-velocity inversions of the derived compliance measurements, using the Metropolis-Hastings algorithm, illuminate the Vs structure of the oceanic crust down to 8 km. Low Vs anomalies in the crust at the SWIR are consistent with significant serpentinization of a crustal component of tectonically exhumed mantle-derived peridotites.

SummaryIn this study, we present a new algorithm and accompanying software for 3D gravity inversion of density structures. The algorithm combines the strengths of the Bayesian approach-which incorporates prior model information through a variogram model-with the advantages of the Tikhonov regularization framework to address the challenge of depth resolution. We also provide a detailed derivation of the procedure for calculating and fitting the 3D experimental variogram, which serves as a fundamental input to the algorithm. The software implementing the proposed algorithm was developed using the widely adopted computational programming language Matlab. To evaluate its effectiveness, we conducted four representative experiments, ranging from simple to complex scenarios. The synthetic results demonstrate that incorporating model covariance constraints yields a more localized and better-focused density distribution compared to results obtained without such constraints. Additionally, we tested the algorithm's robustness by introducing noise into the observation data. The results show that the proposed method is resistant to noise and maintains strong performance. Finally, we applied the algorithm and software to real field data and compared the results with those from previous studies. The comparison confirms that our method is capable of producing reliable, high-resolution 3D density models, with the added advantage of integrating prior information.

SummaryReliable earthquake forecasts depend on modellers making suitable choices regarding data selection and model implementation. We explore some of these choices in order to construct short-term forecasts for independent background events in Italy, using a Log-Gaussian Cox Process to describe the spatial intensity of seismicity as a function of physical spatial covariates and a random field. We explore correlations between a range of physical covariates available for Italy, and investigate how they might contribute to observed spatial seismicity patterns. We find that historic smoothed seismicity, strain rates and elevation are most useful in describing instrumentally-observed seismicity. We use point process intensity models constructed with combinations of these covariates with past seismicity to generate Bayesian earthquake forecasts, constructing simulated catalogue forecasts of future earthquake catalogues in Italy, discussing modelling choices we make along the way. Finally, we test the performance of these forecasts in a pseudo-prospective manner for the 2010-2021 period and using historical data to assess how well they might perform for extreme events. While the inclusion of historic events improves forecasts, models without historic seismicity also perform well, indicating a degree of stationarity in the process. We demonstrate that combinations of covariates can add predictive power over modern catalogue data alone, and provide a framework for future modelling.

SummaryThe rapid development of the global economy brings frequent artificial seismic events related to industrial activities. The identification of faults revealed by natural earthquakes, along with the temporal and spatial distribution of b-values derived from these seismic events, is influenced by their occurrence. Correcting identification of seismic events has significant influence on seismic risk assessment. In this paper, we present a tri-branch convolutional neural network (CNN) model, known as EQTypeNet, to automatically classify three classes of events (natural earthquakes, explosions, and collapses) in China, using waveform, spectrogram, and amplitude ratio features as inputs. The EQTypeNet has developed in two steps: 1. Binary classification: Mitigating data imbalance by supplementing non-natural events outside the test set area. In the classification of explosions and earthquakes, we obtained an individual station macroF1 of 0.99 and a classification accuracy of 98.7 per cent. In the classification of collapses and earthquakes, we obtained an individual station macroF1 of 0.99 and a classification accuracy of 99.5 per cent. 2. Ternary classification of natural earthquakes, explosions, and collapses: Enhancement of the model's classification in small areas by transfer learning. In the transfer learning of Shanxi and Northeast China, the F1-score of individual stations reached 0.98 and 0.97, with corresponding accuracies of 98.6 per cent and 96.8 per cent. The EQTypeNet applies to different magnitude ranges for seismic event discrimination and shows a satisfactory application performance in a large area. It will be expected to achieve higher performance in the ternary classification by further supplementing the training with suitable physical features.

SummaryIn February 2025, a strong seismic crisis occurred 35 km northeast of Santorini, an active volcano located in South Cyclades, Greece, a region of distributed extensional active faulting. The GNSS data shows an inflation of the volcano since August 2024. We model it with a magma source of 7.7 106 m3 located 3.1 km under the north-central floor of the caldera, near the inflation centre of 2011–2012. After 24 January 2025, the seismic activity, until then localised within the caldera, shifted offshore Santorini and increased with eight Mw ≥ 5 events and ground motion exceeding 5 mm at Syros at 110 km from the epicentres. The GNSS data is consistent with a model of dislocation involving a south-east dipping normal fault located between the Kolumbo submarine volcano and the Anydros islet, 18 km long, 12 km wide, with a tip at 7.5 km depth and ~3 m of slip. The ~17.5 1018 Nm corresponding seismic moment, much greater than the ~1.5 1018 Nm of the recorded earthquakes, reveals the probable occurrence of a slow-slip earthquake of equivalent magnitude 6.8. This event might have been triggered by a small dyke injected from the inferred source beneath Santorini. However, the subsidence recorded at Santorini and Anydros is incompatible with the hypothesis of a large dyke injected beneath Kolumbo-Anydros.

SummarySeasonal variations in Nepalese seismicity have been reported with varying degrees of confidence. We re-investigate these claims by analysing 20 years of Nepalese seismicity before the 2015 Gorkha earthquake, as detected by the Nepalese national network, and focusing on earthquakes located along the eastern and central sections of the Nepalese Main Himalayan Thrust. Using several declustering techniques, we find no statistically robust evidence of seasonal seismicity in the studied record, regardless of magnitude threshold above completeness. This suggests that previously reported seasonality may be restricted to the western section of the Nepalese orogeny, may be an artefact, or may indicate that nucleation times of earthquakes are longer than the year. We also investigate potential annual variations in the Gutenberg-Richter b-value, given its recent observed modulation by transient stressing. Additionally, we use large-scale mass redistribution estimated from the monthly gravity field retrieved from the Gravity Recovery And Climate Experiment and Follow-On (GRACE/-FO) missions, to resolve stress variations at depth induced by transient surface loads. We find that the mean annual b-value peaks when seasonal Coulomb stress rates reach their minimum value at the height of the summer rainy season. When considering the combined effect of tectonic and seasonal loading, this corresponds to a recurring period of stress reversal, when Coulomb stress momentarily decreases. This suggests that periodic clamping of the Main Himalayan Thrust reduces the likelihood of earthquakes growing to larger magnitudes in accordance with hierarchical rupture models. The susceptibility of b-value to stress variations of roughly 0.1 points.kPa−1 is consistent with recent estimates of b-value sensitivity to transient loading, although it remains high when compared to the stress-dependency associated with both static differential stress, and with long-term evolution during the seismic cycle. This discrepancy points to the large impact of stress transients on the dynamics of seismic rupture.

SummaryThe goal of this paper is to present considerations, steps and tools to perform statistical analysis of InSAR time series by relying on multichannel singular spectrum analysis (M-SSA). We apply these tools to Sentinel-1 InSAR time series processed for Pacaya Volcano in Guatemala in two steps. First, we produce, in a data-adaptive way, estimates of data points to obtain evenly sampled time series. The resulting time series are then decomposed using M-SSA into long-periodic nonlinear trends and oscillatory modes providing a sparse representation of the signals present in the data. Combining M-SSA that includes varimax rotation with power spectrum analysis augments the physical interpretability of the InSAR dataset presented herein. Monte Carlo SSA hypothesis testing further helps estimate the statistical significance of the M-SSA modes with respect to a red-noise null hypothesis. The dominant frequencies of the main oscillatory modes retained correlate with frequency peaks of the seasonal variability of the regional hydrological system, as determined from correlograms of rainfall time series. The spatial patterns of the significant modes correlate with three types of geological structures present at Pacaya volcano: the volcanic edifice, the 2010 and 2014 lava flows, and a collapse scarp dividing the volcanic edifice into an eastern and western part. These findings suggest that, when including the complementary tools presented herein, M-SSA is able to provide a reliable statistical picture of InSAR datasets and that the main M-SSA modes are geophysically meaningful.

SummaryMachine-learning-based phase pickers have been successfully leveraged to build high-resolution earthquake catalogs using seismic data on land. However, their performance when applied to ocean bottom seismic (OBS) data remains to be evaluated. In this study, we first adopt three machine-learning-based phase pickers—EQTransformer, Pickblue, and OBSTansformer—to build three earthquake catalogs for the 350-km-long Blanco oceanic transform fault (BTF) based on a year-long OBS deployment. We then systematically compare these catalogs with an existing catalog which utilized a traditional workflow. Results indicate that the Pickblue-based catalog documents more events and/or provides better-constrained locations than the other catalogs. The different performances of the three phase pickers suggest that detailed assessment of catalogs built using automatic workflows is necessary to prevent misinterpretations, especially when applied to regions without training samples. The Pickblue-based catalog reveals seismicity gaps in three extensional segments of BTF, which likely represent aseismic slip zones affected by seawater infiltration. Furthermore, most earthquakes are shallower than the 600 °C isotherm predicted by a half-space conductive cooling model, except for the Blanco Ridge segment which has hosted 80 per cent of the Mw > 6.0 earthquakes along BTF since 1976. These Blanco Ridge deep earthquake clusters can be explained by hydrothermal cooling or the serpentinization of mantle peridotite due to seawater infiltration along conduits created by the deeper ruptures of large earthquakes. Our analyses also demonstrate the importance of careful examination of automatically produced earthquake catalogs since mislocated events can lead to very different interpretations of fault slip modes from seismicity distribution.

SummarySeismic wave multiples have historically been treated as noise when performing velocity model building or imaging and have, therefore, typically been attempted removed from the data during preprocessing. Recently, there has been interest in developing methods to utilise their information to better illuminate areas poorly covered by primary seismic wave energy. These often either predict and iteratively apply information from reflected wave multiples or incorporate multiple information by retaining surface reflected waves in a full/reflection waveform inversion framework. Such approaches have provided insight into the behaviour of reflected wave multiples, both alone and when mixed with multiples of other wave types. Much less is known about the behaviour of diving wave multiples alone. In velocity model building using full waveform inversion, diving waves are a significant source of the low-frequency information needed to update low wavenumbers. Here, we derive equations for the ray paths, traveltime, relative geometric spreading, and acoustic reflection coefficient of diving wave multiples when the vertical velocity increases linearly with depth in both an acoustic isotropic model and an acoustic factorized vertical transverse isotropic anisotropic model with constant anisotropy. These are used to investigate the effect of increasing the number of surface reflections. We also study their interference by deriving equations for the amplitude spectrum and traveltime of a wavetrain consisting of a sum of diving wave multiples, where the number of surface reflections increases linearly. Our results show that the properties of the diving wave multiples converge to those of the direct wave as the number of surface reflections approaches infinity. For the relative geometrical spreading in the acoustic factorized vertical transverse isotropic anisotropic model, this leads to a new equation describing the direct wave when the vertical velocity increases linearly with depth. This equation has an anisotropic dependence solely on the η parameter, indicating that η has the largest anisotropic influence on the relative geometrical spreading of the diving wave and its multiples. Both the amplitude spectrum and traveltime of the wavetrains show signs of dispersion. The dispersion primarily arises from the traveltime differences between the diving wave and the diving wave multiple that undergoes one surface reflection. In the amplitude spectrum, this causes notches to appear at different offsets that shift downward in frequency as the offset increases. Adding diving wave multiples that undergo more reflections to the sum significantly changes the amplitudes at short offsets. At mid to long offsets, this primarily affects the amplitude around existing minimum and maximum points, with new small local peaks and troughs additionally being induced at long offsets. The notches in the amplitude spectrum cause complex phase behaviour at different offsets, resulting in difficulties computing the traveltime beyond their appearance. Before this behaviour emerges, we observe low-frequency wavetrains to arrive earlier than those with higher frequencies. We attribute this to the lower-frequency wavetrains sampling the higher seismic velocities deeper in the subsurface. Overall, our results show a significant difference between properties in the acoustic isotropic model and the acoustic factorized vertical transverse isotropic anisotropic model.

SummaryEfficient simulation of controlled-source electromagnetic (CSEM) methods in complex geological settings is critical for accurate subsurface imaging and resource evaluation. This study presents an advanced adaptive finite element algorithm leveraging non-conforming arbitrary hexahedral meshes and combining the flexibility of adaptive refinement with the accuracy of high-order basis functions to solve challenges in electromagnetic field modelling. We present a high-order finite element algorithm for electromagnetic geophysics that supports adaptive refinement on arbitrary hexahedral meshes, a capability not commonly available in existing implementations. The proposed algorithm resolves sign conflict issues inherent in non-conforming meshes by employing global vertex indices and adjusting constraint matrices, ensuring the tangential continuity of the electromagnetic field across the computational domain. Adaptive refinement is guided by goal-oriented error estimation, enabling precise solving for the electromagnetic response while maintaining computational efficiency. Numerical examples validate the algorithm’s accuracy and performance. In a steel-cased well model, the computational accuracy of our algorithm is verified by comparison with the 3-D cylindrically symmetric algorithm of the open-source software SimPEG. Furthermore, it is demonstrated that arbitrary hexahedral meshes can achieve accurate modelling of the steel casing with a reduced mesh count. The application to the Vallès Basin model demonstrates substantial computational savings using arbitrary hexahedral meshes compared to tetrahedral meshes. Finally, the algorithm successfully models a large-scale oil-gas reservoir, illustrating its robustness in addressing complex geological scenarios and its potential for monitoring dynamic processes such as water intrusion during resource extraction. These results establish the proposed method as a reliable and efficient forward modelling tool for controlled-source electromagnetic applications, providing significant improvement in flexibility, accuracy, and computational efficiency.

SummaryThe base of the Earth's lower mantle is characterized by large seismic velocity anomalies, known as large low-velocity provinces (LLVPs) (Garnero et al., 2016). There are several hypotheses related to the origin of LLVPs, such as remnants of Earth's early differentiation (Labrosse et al., 2007; Lee et al., 2010) and buried relics of proto-Earth's mantle after the Moon-forming giant impact (Yuan et al., 2023). However, the geodynamical implications, such as the role of LLVPs as driving mechanisms of plumes or subducted slabs, are not well resolved because some observations of the polarization of seismic velocity at LLVPs use the azimuthal anisotropy of shear wave splitting. Here, we combine new observations of antipodal PKPab seismic waves with the adjoint method to perform an inversion of the radially anisotropic Vp structure at the base of the lower mantle. We have carefully examined antipodal stations with sufficient signal to noise ratios for both the vertical and horizontal components over the past 30 years and selected 23 source‒receiver pairs with epicentral distances greater than 178.0 degrees and Mw values less than 7.0. We calculate synthetic seismograms with an accuracy of 6.9 s and perform an inversion of the radially anisotropic Vp structure at the base of the lower mantle by the adjoint method. The results of our inversion show that vertically polarized Vp is dominant within the LLVPs of the Pacific and African regions. These features are characterized by relatively small spots of high vertically polarized Vp anomalies, which may be interpreted as the locations of ascending mantle plumes inside LLVPs in the Pacific region.

AbstractSeismic inversion translates seismic data into subsurface elastic property models, enabling geophysicists to better understand underground rocks and fluids. Due to the inherently ill-posed nature of this inverse problem, accurately capturing the uncertainty associated with the solution is essential for reliable interpretations. Traditional Bayesian inversion methods, such as Markov Chain Monte Carlo (MCMC) and Laplace approximations, have been employed for this purpose but face significant limitations in terms of scalability and computational efficiency for large-scale problems. Combined with deep learning, Variational Inference (VI) has emerged as a promising alternative, striking a balance between computational efficiency and flexibility (i.e. the ability to approximate complex posterior distributions). However, selecting an appropriate proposal distribution remains a key challenge, as it directly influences the quality of the estimated posterior distribution. In this study, we extend IntraSeismic, an implicit neural representation (INR)-based framework for seismic inversion applications, to Bayesian inversion using VI with different parameterizations of the proposal distribution. We introduce two methods: B-IntraSeismic (BIS), which uses a mean-field Gaussian proposal, and B-IntraSeismic with Conditional Normalizing Flows (BIS-Flow), which utilizes a mean-field unparameterized proposal distribution to better capture deviations from Gaussianity in the posterior distribution. These methods are evaluated on a synthetic dataset (Marmousi) and two field data (Volve and Sleipner). Our results indicate that both BIS and BIS-Flow can accurately capture structural details and produce high-resolution mean models and standard deviation maps. BIS-Flow is also shown to be able to model complex posterior distributions, offering a more comprehensive characterization of uncertainty while maintaining computational feasibility.

AbstractThe phenomenon of elastic wave conversions, where acoustic, pressure (P-) waves are converted to elastic, shear (S-) waves, and vice-versa, is commonly disregarded in seismic imaging. This can lead to lower-quality images in regions with strong contrasts in elastic parameters. While a number of methods exist that do take wave conversions into account, they either deal with P- and S-waves separately, or are prohibitively computationally expensive, as is the case for elastic Full-Waveform Inversion. In this paper an alternative approach to taking converted waves into account is presented by extending Full Wavefield Migration (FWM) to account for wave conversions. FWM is a full-wavefield inversion method based on explicit, convolutional, one-way propagation and reflection operators in the space-frequency domain. By applying these operators recursively, multi-scattering data can be modelled. Using these operators, the FWM algorithm aims to reconstruct the reflection properties of the subsurface (i.e. the ‘image’). In this paper, the FWM method is extended by accounting for wave conversions due to angle-dependent reflections and transmissions using an extended version of Shuey’s approximation. The resulting algorithm is tested on two synthetic models to give a proof of concept. The results of these tests show that the proposed extension can model wave conversions accurately and yields better inversion results than applying conventional, acoustic FWM.

SummaryGlobal seismology mainly uses seismic waves propagating in the sagittal plane along the great circle path (GCP). However, heterogeneities in the mantle laterally deviate the path of seismic signals, which arrive out-of-plane (OOP) at arrays of sensors at teleseismic distances. Detection and back-projection of these signals have, in the past, provided independent evidence for the location of distant subducted slabs in the deep mantle, complementing global tomographic imaging. To infer physical properties of these subducted slabs, 3D waveform modeling of OOP waves for a finite-thickness slab is needed but still missing. In this study, we conduct a series of synthetic tests using a spectral element solver. We test the detectability of OOP signals and, by progressively adding complexities, we evaluate to which extent these signals can be used to infer physical properties of the modeled slab. We carry out three-component array analysis and investigate focal mechanism dependency. Our results show that the transverse component might be the best candidate for such studies, also for P-to-P OOP signals. Vertical and radial component recordings are usually dominated by P-SV energy arriving from the earthquakes along the GCP, which masks possible OOP signals. Contrary, the transverse component filters out any P-SV energy arriving directly from the source and, owing to its intrinsic directionality, allow for higher resolution measurement of P-to-P OOP signals. This is especially the case prior to the arrival of the S-wavefield. We pick a series of OOP arrivals which are back-projected using a multi-phase trial-and-error approach, that is considered successful only when different OOP seismic phases converge to the modeled (true) structure. We retrieve the location of the slab, its bottom and top edges, and its thickness in the lower mantle. These inferences are tested against varying topography, orientation and size of the modeled slab. The insights gained with modeling are confirmed with real data examples, supporting higher resolution mapping of 3D mantle structure based on OOP seismology.

SummaryWe investigate the lithospheric structure of the Southwest Iberian margin along an active seismic profile southwest of São Vicente Cape, ranging from the southern Tagus Abyssal Plain to the westernmost part of the Gulf of Cadiz. This profile, approximately 320 km long, intersects almost perpendicularly three major thrust faults: the Tagus Abyssal Plain, Marquês de Pombal and Horseshoe faults. The crustal structure, derived from spatially coincident wide-angle seismic (WAS) and multichannel seismic (MCS) data, was validated and constrained using gravimetric data. Joint travel-time inversion of refracted phases identified in WAS and reflected seismic phases from both WAS and MCS records were used to build a detailed two-dimensional P-wave velocity (Vp) structure. The resulting model reveals a Vp distribution with abrupt lateral velocity and structural variations, characterized by a rugged basement top and sharp changes in crustal thickness. Three main lithospheric domains consisting of continental, oceanic, and exhumed mantle affinity were identified from south to north. The travel-time inversion of the deepest reflected seismic phases reveals four major southeast-dipping reflectors, likely corresponding to major regional thrust faults with significant seismic and tsunamigenic potential. Integrating the modelled and interpreted seismic results with the locations of recent well-constrained earthquakes suggests that the Marquês de Pombal and Tagus Abyssal Plain extend deeper than previously thought, with fairly high seismic activity in the deep levels. This has significant implications for their seismogenic potential and should be taken into account for accurate assessment of seismic hazards in the region.

SummaryFirst-arrival traveltime tomography (FATT) is a widely used method for characterizing near-surface velocity structures in geotechnical engineering and resource exploration. We introduce an improved version of FATT, named first-arrival adjoint tomography (FAAT), which involves the joint inversion of first-arrival absolute traveltimes and differential traveltimes. Unlike absolute traveltimes, differential traveltimes, derived from common sources or receivers, offer heightened sensitivity to fine-scale structures near the receivers or sources, respectively. This dual sensitivity makes FAAT particularly effective in imaging highly heterogeneous media. However, the proximity of rays associated with differential traveltimes can lead to the instability of the inversion. To overcome this challenge, we simultaneously incorporate absolute traveltimes and differential traveltimes to update the velocity model. This approach improves the stability of the inversion process, leading to improved resolution of inverted results. We specifically employ the fast-sweeping method to solve the factored eikonal equation, providing robust solutions in models with complex geological structures. Furthermore, we address the inverse problem by computing the gradient of data misfit using the efficient adjoint-state method. Through numerical testing, we validate the effectiveness of FAAT in comparison to that using only absolute or differential traveltimes. Finally, we apply the proposed FAAT method to near-surface site characterization at Bishan-AMK Park in Singapore. Compared with FATT and validated against borehole data, FAAT demonstrates its ability to establish more reliable velocity models, revealing finer details and substantially improving geological interpretation.