Geophysical Journal International

SummaryAround the world the Earth’s crust is blanketed to various extents by sediment. For continental regions, knowledge of the distribution and thickness of sediments is crucial for a wide range of applications including seismic hazard, resource potential, and our ability to constrain the deeper crustal geology. Excellent constraints on the sediment thickness can be obtained from borehole drilling or active seismic surveys. However, these approaches are expensive and impractical in remote continental interiors such as central Australia. Recently, a method for estimating the sediment thickness using passive seismic data, the collection of which is relatively simple and low-cost, was developed and applied to seismic stations in South Australia. This method uses receiver functions, specifically the time delay of the P-to-S converted phase generated at the sediment-basement interface, relative to the direct-P arrival, to generate a first order estimate of the thickness of sediments. In this work we expand the analysis to the vast array of over 1500 seismic stations across Australia, covering an entire continent and numerous sedimentary basins that span the entire range from Precambrian to present-day. We compare with an established yet separate method to estimate the sediment thickness, which utilises the autocorrelation of the radial receiver functions to ascertain the two-way travel-time of shear waves reverberating in a sedimentary layer. Across the Australian continent the new results match the broad pattern of expected sedimentary features based on the various geological provinces. We are able to delineate the boundaries of many sedimentary basins, such as the Eucla and Murray Basins, which are Cenozoic, and the boundary between the Karumba Basin and the mineral rich Mount Isa Province. Contrasts in seismic delay time across these boundaries are upwards of 0.4 s. The delay signal is found to diminish to <0.1 s for older Proterozoic basins, likely due to compaction and metamorphism of the sediments over time. As an application of the method, a comparison with measurements of sediment thickness from local boreholes allows for a straightforward predictive relationship between the delay time and the cover thickness to be defined. This offers future widespread potential, providing a simple and cheap way to characterise the sediment thickness in under-explored areas from passive seismic data.

SummaryEarthquake forecasting poses significant challenges, especially due to the elusive nature of stress states in fault systems. To tackle this problem, we employ features derived from seismic catalogs obtained from acoustic emission (AE) signals recorded during triaxial stick-slip experiments on natural fractures in three Westerly granite samples. We extracted 47 physically explainable features from AE data that described spatio-temporal evolution of stress and damage in the vicinity of the fault surface. These features are then subjected to unsupervised clustering using the K-means method, revealing three distinct stages with a proper agreement with the temporal evolution of stress. The recovered stages correspond to the mechanical behavior of the rock, characterized as initial stable (elastic) deformation, followed by a transitional stage leading to an unstable deformation prior to failure. Notably, AE rate, clustering-localization features, fractal dimension, b-value, interevent time distribution, and correlation integral are identified as significant features for the unsupervised clustering. The systematically evolving stages can provide valuable insights for characterizing preparatory processes preceding earthquake events associated with geothermal activities and waste-water injections. In order to address the upscaling issue, we propose to use the most important features and, in case of normalization challenge, removing non-universal features, such as AE rate. Our findings hold promise for advancing earthquake prediction methodologies based on laboratory experiments and catalog-driven features.

SummarySubduction earthquakes show complex spatial and temporal rupture patterns, exhibiting events of varied sizes, which rupture distinct or overlapping fault segments. Elucidating first-order controlling conditions of rupture segmentation and return periods of large earthquakes is therefore critical for seismic and tsunami hazard estimations. The Chilean subduction zone frequently hosts several Mw > 8 earthquakes, with heterogeneous recurrence rates and locations. Here, we implement 3D quasi-dynamic rate and state frictional models to investigate the role of plate interface geometry on the distribution of interseismic coupling and coseismic ruptures in Central Chile. First, we develop synthetic-parametric models that show how dip and strike variations may increase the probabilities to produce partial seismic barriers, which tend to avoid the production of large earthquake ruptures and modulate rupture lengths. Then, we simulate the subduction seismic cycle processes on Central Chile (25ºS-38°S), imposing depth-dependent frictional properties on a realistic non-planar 3D subduction interface geometry. Similar to results obtained for synthetic-parametric models, after 5000 years of simulation, regions with abrupt dip or strike changes increase the probabilities of stopping coseismic propagation of simulated Mw8.0-9.0 earthquakes. Our simulated earthquake sequences on the Central Chile subduction zone delimit rupture areas that match geometrical interface features and historical earthquakes, results that point to the crucial role of fault interface geometry on seismic cycle segmentation along strike.

SummaryTrain traffic is a powerful source of seismic waves, with many applications for passive seismic imaging. Seismic signals were recorded a few meters away from the railway track. These records display harmonious waveforms below 15 Hz for trains driving at speeds of around 100 km/h. The sensors record an apparent wavelet emitted by the interaction of the axle on a few of the closest sleepers. From this, we build a simple modeling tool using shifted wavelets to simulate a train signal. The relationship involves the varying train speed, the distances between each axle, and the wavelet emitted by each axle. We propose a non-linear deconvolution method to invert this relationship. We use a local minimisation algorithm with gradients derived by the adjoint state method, and use a frequency continuation technique. A linearised picking-based inversion initializes the non-linear inversion. On real data, we apply this automatic workflow to 300 train passages, with an excellent match between the best simulation and the data. We identify the trains decelerating as they enter a train station. We also identify the train type based on inverted wheel spacing with centimetric accuracy. The inverted wavelets are consistent with the assumption that trains emit seismic waves by bending the rail above sleepers, although the theory does not explain why the inverted wavelet is not zero-phase. This automated kinematic inversion algorithm may allow for contact-less railway monitoring, and be used for source characterisation for subsurface monitoring below railway tracks.

SummaryThis study presents a new robust statistical framework, in which to measure relative differences, or deviations from a hypothetical reference value, of Gutenberg-Richter b-value. Moreover, it applies this method to recent seismicity in Italy, to find possible changes of earthquake magnitude distribution in time and space. The method uses bootstrap techniques, which have no prior assumptions about the distribution of data, keeping their basic features. Excluding Central Italy, no significative b-value variation is found, revealing that the frequency-magnitude distribution exponent is substantially stable or that data are not able to reveal hidden variations. Considering the small size of examined magnitude samples, we cannot definitively decide if the higher b-values in Central Italy, consistently founded by all applied tests, have a physical origin or result from a statistical bias. In any case, they indicate short-lived excursions which have a temporary nature and, therefore, cannot be associated solely to spatial variations in tectonic framework. Both the methodological issues and the results of the application to seismicity in Italy show that a correct assessing of b-value changes requests appropriate statistics, that accurately quantify the low accuracy and precision of b-value estimation for small magnitude samples.

SummaryInterferometric synthetic aperture radar (InSAR) decorrelation that creates great challenges to phase unwrapping has been a critical issue for mapping large earthquake deformation. Some studies have proposed a “remove-and-return model” solution to tackle this problem, but it has not been fully validated yet, and therefore has rarely been applied to real earthquake scenarios. In this study, we use the 2023 Mw 7.8 and 7.6 earthquake doublet in Turkey and Syria as a case example to develop an iterative modeling method for InSAR-based coseismic mapping. We first derive surface deformation fields using Sentinel-1 offset tracking and Sentinel-2 optical image correlation, and invert them for an initial coseismic slip model, based on which we simulate InSAR coseismic phase measurements. We then remove the simulated phase from the actual Sentinel-1 phase and conduct unwrapping. The simulated phase is added back to the unwrapped phase to produce the final phase measurements. Comparing to the commonly-used unwrapping method, our proposed approach can significantly improve coherence and reduce phase gradients, enabling accurate InSAR measurements. Combining InSAR, offset tracking and optical image correlation, we implement a joint inversion to obtain an optimal coseismic slip model. Our model shows that slip on the Çardak Fault is concentrated on a ∼100 km segment; to both ends, slip suddenly diminished. On the contrary, rupture on the East Anatolian Fault Zone propagated much longer as its geometry is fairly smooth. The iterative coseismic modeling method is proven efficient and can be easily applied to other continental earthquakes.

SummarySite amplification factors (SAFs) of seismic ground motions are essential in evaluating and estimating seismic hazards. In our previous study, the authors proposed a simple and cost-effective method to estimate a SAF based on a deep neural network (DNN) model and microtremor horizontal-to-vertical spectral ratio (MHVR). Since the previous DNN model was based on the observed SAFs and MHVRs within a limited district in Japan, the applicability of the previous model to non-source regions with different site conditions was limited. This study explored the application of a transfer learning (TL) technique to adapt an existing (pre-trained) DNN model to new regions and a different database. The SAFs obtained through generalized spectral inversion technique (GIT) at the seismic observation stations (K-NET and KiK-net) in Japan were collated as the ground truth for site effects. MHVRs recorded at the stations in several districts of Japan were collected to construct a dataset for the development of the TL model. Subsequently, a TL model was constructed, leveraging the neural network layers and their weights from the pre-trained model while incorporating additional neural network layers to enhance the performance. During the training process, a total data set of 112 sites was divided into training set, validation set, and external test set by 1:1:5. Utilizing a cross-validation strategy, the residuals between pSAFs (pseudo-SAFs) estimated by the TL model and the observed SAFs were analyzed for the external test set containing 80 sites. The results showed that the TL model outperformed the pre-trained DNN model. The cross-validation results demonstrated that almost consistent prediction results were obtained regardless of any combination of 16 sites selected as the training set. Furthermore, by contrasting the influence of varying training set sizes on the performance of the TL model and comparing the TL model to a DNN model with an extended training set, the effectiveness of constructing the model with the limited number of data (16 sites) was ascertained. Finally, the effectiveness and limitations of the TL model were evaluated using MHVRs with peak frequencies falling outside the training set's range.

SummaryAs the high-frequency analog to field-scale earthquakes, acoustic emissions (AEs) provide a valuable complement to study rock deformation mechanisms. During the load-stepping creep experiments with CO2-saturated water injection into a basaltic sample from Carbfix site in Iceland, 8791 AE events are detected by at least one of the seven piezoelectric sensors. Here, we apply a cross-correlation-based source imaging method, called geometric-mean reverse-time migration (GmRTM) to locate those AE events. Besides the attractive picking-free feature shared with other waveform-based methods (e.g., time-reversal imaging), GmRTM is advantageous in generating high-resolution source images with reduced imaging artifacts, especially for experiments with relatively sparse receivers. In general, the imaged AE locations are found to be scattered across the sample, suggesting a complicated fracture network rather than a well-defined major shear fracture plane, in agreement with X-ray computed tomography imaging results after retrieval of samples from the deformation apparatus. Clustering the events in space and time using the nearest-neighbor approach revealed a group of “repeaters”, which are spatially co-located over an elongated period of time and likely indicate crack, or shear band growth. Furthermore, we select 2196 AE events with high signal-noise-ratio (SNR) and conduct moment tensor estimation using the adjoint (back-propagated) strain tensor fields at the locations of AE sources. The resulting AE locations and focal mechanisms support our previously assertion that creep of basalt at the experimental conditions is accommodated dominantly by distributed micro-cracking.

SummaryMuography is an imaging technique that relies on the attenuation of the muon flux traversing geological or anthropogenic structures. Several simulation frameworks help to perform muography studies by combining specialised codes: for muon generation through muon transport to muon detector performance. This methodology is precise but requires significant computational resources and time. We present an end-to-end python-based MUographY Simulation Code, which implements a muography simulation framework capable of rapidly estimating muograms of any geological structure worldwide. This framework considers the generated muon flux as the observation point; the energy loss of muons passing through the geological target; the integrated muon flux detected by the telescope and estimates the 3-dimensional density distribution of the target using Algebraic Reconstruction Techniques. The simulations ignore the relatively small muon flux variance caused by geomagnetic effects, solar modulation, and atmospheric conditions. We validate the code performance by comparing our simulation results with data from other frameworks.

SummaryThe diffusion field theory has been widely used to interpret ambient noise wave fields. Based on this theory, 1D subsurface velocity structure inversion method is developed. However, few studies have referred to the Noise Horizontal to Vertical spectral ratio (NHV) inversion of 3D subsurface velocity structures, and almost no effective 3D NHV inversion tools have been developed. To develop a useful tool for obtaining 3D soil layer velocity structures, we combined the NHV forward calculation formula derived from diffusion field theory with the guided Monte Carlo algorithm and then extended the single-point NHV inversion to multipoint joint inversion through a joint objective function. Subsequently, a new three-dimensional soil layer velocity structure inversion method was proposed. Subsequently, a synthetic 2D case was used to verify the proposed method. Finally, the proposed method was applied to the Xiangtang Array in Tangshan, China, to identify the 3D velocity structures of the site based on noise observations. The results show that the proposed multipoint joint 3D inversion method is effective for identifying 3D underground velocity structures.

SummaryThe physical mechanisms controlling deep-focus earthquakes, or those observed at depths greater than 300 km, remain enigmatic. The leading processes by which deep-focus earthquakes are thought to occur include transformational faulting, thermal runaway, and dehydration embrittlement, but distinguishing observations in support of one or more prevailing mechanisms are needed. In this study, we use a modified back-projection method, data recorded by the Hi-net array in Japan, and a three-dimensional velocity model to produce source images of 19 deep-focus earthquakes within the Izu-Bonin subduction zone. We find that the rupture properties and fault plane orientations of imaged events separate according to reported moment magnitude, indicating the distinct operation of two moment-dependent causal mechanisms of deep-focus earthquakes in this region. We discuss these results in the context of previous observational, laboratory, and numerical studies and emphasize the importance of continued research to validate the dual-mechanism hypothesis both in and outside Izu-Bonin. Such work may not only improve our understanding of the nucleation and propagation of deep-focus earthquakes, but also help clarify slab structure and subduction zone dynamics.

SummaryMuography is a promising and rapidly developing physical prospecting technique based on the attenuation of muon flux. The feasibility and potential of applying muography to mining were presented in a small number of previous case studies. In this work, we aimed to address three challenges that limit the applicability and efficiency of muography in mineral exploration: (1) application to low-density-contrast ore body exploration, (2) analysis of altitudinal impacts on measurements, and (3) precise and efficient reconstruction. We conducted the first case of applying muography to the exploration for low-density-contrast ore bodies. Six muon imaging systems were placed underground to collect surviving muons for roughly half a year. We analyzed the altitudinal impact on the data measurements and proposed a simplified method to eliminate it. We also developed the seed algorithm, a novel three-dimensional reconstruction algorithm, that can reconstruct anomalies located far away from the detectors and avoid their elongation along the observed directions. Benefiting from the seed algorithm, a low-density-contrast orebody and a limonitic siliceous slate structure were reconstructed, demonstrating the sensitivity of this technique in density distinction; discoveries of several mined-out areas are important for accident avoidance; and reconstruction of the stope and scarps served as strong circumstantial evidence of the reliability of the result. The success of this experiment shows the great value of muography in the economic, research and safety aspects of mineral exploration and inspection. Moreover, the overcoming of challenges is helpful for the development of muography, making it an effective and competitive technique in this field.

SummaryThree-dimensional finite element models, which can handle the stress perturbations caused by subsurface mechanical heterogeneities and fault interactions, have been combined with the finite source inversion to estimate the coseismic slip distribution over the fault plane. However, the mesh grid for discretizing the governing equations in the finite element model significantly affects the numerical accuracy. In this study, we performed kinematic finite source inversion with idealized (regular observation point array; M1A–M1D) and regional (GEONET stations in Japan; M2A–M2H) models with different mesh sizes to quantitatively analyse the effect of the mesh grid size around the fault plane on the inverted fault slip distribution. Synthetic observation data vectors obtained from the finest models (M1A and M2A) are compared with those from the coarser models (M1B–M1D and M2B–M2H), which were adopted to construct Green's function matrix. We found that the coarser mesh models derived a smaller surface displacement, leading to a decrease in the norm of Green's function matrix, which in turn influences the fault slip magnitude from the finite source inversion. Finally, we performed the source inversion for the fault slip distribution of the 2011 Mw 9.0 Tohoku-Oki earthquake using the coseismic surface displacements recorded at the GEONET and seafloor stations and finite element modeling. By reducing the mesh size on the fault, we confirmed that the estimated magnitude of fault slip converged to approximately 80 m, which is consistent with the range of fault slip amounts from previous studies based on the Okada model. At least 0.88 million total domain elements and a 6.7 km2 mesh size on the fault plane with an area of 240 × 720 km2 are required for the convergence of the fault slip. Furthermore, we found that the location of the maximum fault slip is less sensitive to the mesh size, implying that source inversion based on a coarse mesh model (i.e. less than 0.5 million elements and > ∼60 km2 mesh size) can quickly provide the rough fault slip distribution.

SummaryA similar waveform search plays a crucial role in seismology for detecting seismic events, such as small earthquakes and low-frequency events. However, the high computational costs associated with waveform cross-correlation calculations represent bottlenecks during the analysis of long, continuous records obtained from numerous stations. In this study, we developed a deep-learning network to obtain 64-bit hash codes containing information on seismic waveforms. Using this network, we performed a similar waveform search for ∼35 million moving windows developed for the 30 min waveforms recorded continuously at 10 MHz sampling rates using 16 acoustic emission transducers during a laboratory hydraulic fracturing experiment. The sampling points of each channel corresponded to those of the 5.8-year records obtained from typical seismic observations at 100 Hz sampling rates. Of the 35 million windows, we searched for windows with small average Hamming distances among the hash codes of 16 channel waveforms against template hash codes of 6057 events that were catalogued using conventional autoprocessing techniques. The calculation of average Hamming distances is 1430–1530 times faster than that of the corresponding network correlation. This hashing-based template matching enabled the detection of 23,462 additional events. We also demonstrated the feasibility of the hashing-based autocorrelation analysis, where similar event pairs were extracted without templates, by calculating the average Hamming distances for all possible pairs of the ∼35 million windows. This calculation required only 15.5 h under 120 thread parallelisation. This deep hashing approach significantly reduced the required memory compared with locality-sensitive hashing approaches based on random permutations, enabling similar waveform searching on a large-scale dataset.

SummaryIrazú and Turrialba are a twin volcanic complex that marks a distinct stop in volcanism along the Central America volcanic arc. We present a new travel-time velocity model of the crust beneath Irazú and Turrialba volcanoes, Costa Rica, and interpret it considering the results of previous ambient noise tomographic inversions. Data were acquired by a temporary seismic network during a period of low activity of the Irazú-Turrialba volcanic complex in 2018-2019. Beneath the Irazú volcano, we observe low P-wave velocities (VP = 5 km s-1) and low velocity ratios (VP/VS = 1.6). In contrast, below the Turrialba volcano, we observe a low S-wave velocities (VS = 3 km s-1) and a high VP/VS (= 1.85) anomaly. We found that locations of low VP and VS anomalies (- 15 per cent) correspond well with shear wave velocity anomalies retrieved from ambient noise tomography. At shallower depths, we observe high VP and VS anomalies (+ 15 per cent) located between the summits of the volcanoes. Sub-vertical velocity anomalies are also observed at greater depths, with high VP and VS anomalies appearing at the lower limits of our models. We propose a complex structure of an intermediate magmatic reservoir, presenting multi-phase fluid states of a liquid-to-gas transition beneath Irazú and a juvenile store of magmatic fluid beneath Turrialba, while shallow fluid transport provides evidence of magmatic-hydrothermal interactions.

SummaryIncreasing magnetisation within mature hydrocarbon reservoirs provides a new technique in identifying oil-water contacts (OWCs) in cored wells with the potential to assess yield thereby reducing the need for further exploration. Authigenic precipitation of magnetic minerals at OWCs may also help locate paleocontacts (PCs), where structural changes to the petroleum system have caused hydrocarbon remigration. This study determines the magnetic characteristics of magnetic enhancements at OWCs and possibly PCs in silliclastic and carbonate reservoirs at the Wytch Farm oil field, Wessex Basin, UK. Increases in saturation magnetisation and susceptibility are observed at the OWC in 11 of the 12 analysed cored reservoirs owing to the increased presence of magnetite and vivianite. Geochemical analysis and shallow reservoirs suggest biogenic and inorganic mineral precipitation is extensive at the OWC depending on iron, sulphur, and phosphorus availability. Similar magnetic characteristics have been observed in magnetic enhancements above the OWC in numerous wells which may represent OWCs before a basin-wide easterly tilt caused hydrocarbon remigration in the Cenozoic. Multiple magnetic enhancements above the OWC in westerly onshore wells, suggest this remigration may have occurred as numerous phases.

SummaryThe drift history of the Lhasa terrane plays an essential role in understanding the tectonic evolution of the Bangong-Nujiang Tethyan Ocean and the Neo-Tethyan Ocean, as well as the evolutionary history of the Tibetan Plateau. Here, a combined rock magnetic, petrographic, and palaeomagnetic study is performed on the Upper Permian–Lower Triassic limestones (∼259–251 Ma) in the western Lhasa terrane. The site-mean direction for the 28 sites is Dg = 32.1°, Ig = 50.3°, kg = 47.9, and α95 = 4.0° in situ and Ds = 342.9°, Is = 32.7°, ks = 43.2, and α95 = 4.2° after tilt-correction, yielding a palaeopole at 68.9° N, 314.4° E with A95 = 4.3°, corresponding to a palaeolatitude of 18.0° ± 4.3° N. The fold tests are not significant because the sampling section shows monoclinic features with minor variations in their bedding attitudes. The palaeopoles for the directions before and after tilt-correction are compared with reliable Late Permian–Palaeogene palaeopoles obtained from the Lhasa terrane. Based on these comparisons, the studied limestones were remagnetized prior to tilting and this remagnetization most likely occurred during the Early Cretaceous. The depositional environment of the limestones may have changed from anoxic to suboxic and oxic during the Early Cretaceous, leading to the oxidation of iron sulphide to authigenic magnetite. Meanwhile, the Late Jurassic–Early Cretaceous convergence between the western Lhasa and Qiangtang terranes may have resulted in tectonic fluid migration and the formation of calcite veins and stylolites in the limestones. This is supported by the presence of small calcite veins and stylolites in some samples, as well as the fact that the framboidal oxides were formerly sulphides (mostly pyrite), implying that the majority of the iron oxides observed in the limestones were authigenic. These processes indicate that chemical remanent magnetization caused by the growth of magnetic minerals related to tectonic fluid migration was most likely the mechanism for the limestone remagnetization.

SummaryIn the immediate vicinity of a source there are strong gradients in the seismic wavefield that are tamed and modified in DAS recording due to combined effects of gauge-length averaging and local stacking on the local strain field. Close to a source broadside propagation effects are significant, and produce a characteristic impact on the local DAS channels. In the presence of topography, of surface or cable, additional effects are introduced that modify the expected signal. All these influences mean that the results of tap tests used to calibrate the channel positions along a DAS cable may give a distorted view of the actual geometry. Such effects can be important for detailed mapping of faulting processes and comparable features.

SummaryAccurate characterization and monitoring strategies are essential for designing and implementing remedial programs for sites polluted with dense non-aqueous phase liquids (DNAPLs). Electrical resistivity tomography (ERT) is a widely used geophysical technique for mapping subsurface features and processes of interest and exhibits desirable characteristics for DNAPL sites due to its ability to gather large volumes of continuous subsurface information in a non-invasive, cost-effective, and time-efficient manner. However, ERT from the surface suffers from poor imaging quality with depth. Enhanced ERT imaging can be obtained via electrodes deployed on the surface and within horizontal boreholes, but so far it has only been investigated for two-dimensional (2D) imaging. This study evaluates the potential of three-dimensional (3D) surface-to-horizontal borehole (S2HB) ERT configurations for imaging 3D DNAPL source zones. Laboratory tank experiments were first conducted with a 3D S2HB ERT configuration, which consisted of a surface grid and a single borehole line of electrodes, being used to monitor DNAPL migration within porous media. Results demonstrate that 3D S2HB ERT with a single borehole provides improved sensitivity at depth, and therefore enhanced imaging compared to conventional 3D surface ERT. Further tank experiments were performed to assess the performance of single borehole S2HB ERT when (i) the distance between surface and borehole is increased, and (ii) additional horizontal boreholes are included. The S2HB ERT with a single borehole significantly outperforms surface ERT at larger depths and performs comparably to S2HB ERT with multiple boreholes. This study suggests that 3D S2HB ERT with a single borehole can provide the enhanced imaging ability needed to map DNAPLs, while also being relatively practical for implementation at field sites.

SummaryAccurate seismic phase detection and onset picking are fundamental to seismological studies. Supervised deep-learning phase pickers have shown promise with excellent performance on land seismic data. Although it may be acceptable to apply them to OBS (Ocean Bottom Seismometer) data that are indispensable for studying ocean regions, they suffer from a significant performance drop. In this study, we develop a generalised transfer-learned OBS phase picker – OBSTransformer, based on automated labelling and transfer learning. First, we compile a comprehensive dataset of catalogued earthquakes recorded by 423 OBSs from 11 temporary deployments worldwide. Through automated processes, we label the P and S phases of these earthquakes by analysing the consistency of at least three arrivals from four widely used machine learning pickers (EQTransformer, PhaseNet, Generalized Phase Detection, and PickNet), as well as the AIC (Akaike Information Criterion) picker. This results in an inclusive OBS dataset containing ∼36,000 earthquake samples. Subsequently, we employ this dataset for transfer learning and utilize a well-trained land machine learning model – EQTransformer as our base model. Moreover, we extract 25,000 OBS noise samples from the same OBS networks using the Kurtosis method, which are then used for model training alongside the labelled earthquake samples. Using three groups of test datasets at sub-global, regional, and local scales, we demonstrate that OBSTransformer outperforms EQTransformer. Particularly, the P and S recall rates at large distances (>200 km) are increased by 68% and 76%, respectively. Our extensive tests and comparisons demonstrate that OBSTransformer is less dependent on the detection/picking thresholds and is more robust to noise levels.