Geophysical Journal International

SummaryControlled-source marine seismic experiments are key in advancing our understanding of the Earth’s subsurface structure to study tectonic, magmatic, sedimentary and fluid flow processes. Joint acquisition of Wide-Angle Seismic (WAS) and Multi-Channel Seismic (MCS) streamer data stands as the most robust approach for marine exploration, however effectively mapping subsurface structure remains challenging. The lack of identifiable refractions as first arrivals at short offsets in WAS data creates shallow subsurface illumination gaps up to 6-8 km offsets around Ocean Bottom Seismometers or Hydrophones (OBS/OBH). This inadequate ray coverage, more pronounced in areas with deeper water column and lower seabed velocities, limits the performance of Travel Time Tomography (TTT) techniques. Velocity determination in the sedimentary layer and reflector location are affected, and errors propagate to deeper layers. This study integrates Downward Continuation (DC) to WAS data. Similarly to our former study where DC is applied to MCS data, redatuming WAS data involves solving the acoustic wave equation backward in time. This process virtually repositions the sources to the seafloor, revealing previously masked near-seafloor refractions as first arrivals. This transformation significantly enhances ray coverage in the shallow subsurface, leading to more accurate determinations of both seismic velocity and reflector geometry. By bridging theoretical concepts with a real data application, this study demonstrates the optimization of field seismic data for improved TTT results. This methodology is particularly beneficial for deep water exploration where spatially coincident WAS and MCS are jointly inverted. In such scenarios, DC-processed WAS data provides the refracted phases key for velocity determinations, and that are typically not present in MCS data due to insufficient streamer length relative to the water column depth. Additionally, we contribute to the community by releasing our open-source, High-Performance Computing (HPC) software for WAS data redatuming.

SummaryMexico City is one of the largest cities in North America, facing high seismic hazards and water supply problems. This paper presents an ambient seismic noise tomography of the city's south area in Xochimilco, where large amplifications have already been registered during subduction earthquakes. Eighty-four seismic stations have been installed, and their records processed. The tomography method combines the inversion of the Horizontal-to-Vertical Spectral Ratios (HVSR) and multimodal dispersion curves. The importance of considering a multimodal approach is justified in light of the complex geological setting. The dispersion curves analysis shows that the surface wave energy is divided over the fundamental and the higher modes, particularly between 50 and 300 m/s, and in the whole frequency range analyzed. We observe a spatially continuous decrease of the dominant peak frequency of the HVSRs toward the lake interior but a heterogeneous amplification. By analyzing the velocity profiles associated with the highest amplifications, we discovered that these latter result from the superposition of several resonance peaks. Their coincidence in frequency is due to the overall constant linear gradient velocity in the sedimentary basin crossed by several low-velocity anomalies due to high water content or high-velocity anomalies due to lavas. Although most of the shallow water is trapped in clay sediment, the velocity model also allows for identifying deeper water reserves. All these analyses are of fundamental importance for the correct seismic mitigation in Mexico City but might also be extended to other cities built on top of sedimentary basins.

SummaryIn geodetic slip inversions, resolution decreases rapidly with depth because deformation data are measured at the ground surface. Traditionally, this decrease in resolution has been attributed to the weak deformation signals caused by slips at greater depths. However, this study demonstrated that the primary cause is the stronger smoothing applied to deeper slips compared to shallower ones. This work proposes a method that scales the Green's functions to equalize smoothing effects across depths. The method's effectiveness was validated through both synthetic tests and real earthquake applications. In synthetic tests, it improved recovery of deep slips in both location and amplitude. When applied to the 2008 MW7.9 Wenchuan earthquake, the method produced smaller slips near the ground surface and larger slips at depth. For the land-based deformation inversion of the 2011 MW9.0 Tohoku earthquake, the method resulted in larger shallow slips near the trench and greater sea-floor uplift compared to the conventional inversion, which is valuable for accurate prediction of tsunami wave height. Additionally, this method may also be applicable to other inversions where smoothing is used and observation amplitudes vary with distance.

SummaryLocating earthquakes plays an important role in the study of seismic activity and geological structures. Traditional methods for locating earthquakes mainly rely on waveform matching and travel time fitting. With the development of artificial intelligence technology, machine learning methods have often been applied to locate earthquakes. However, current machine learning approaches may face challenges related to physical constraints. In this study, we build a 3D U-Net network with station distribution constraints to locate earthquakes. To improve the generalizability of the network model, we apply data augmentation techniques, including data shifting, selection, rotation, and fusion. The location results are evaluated using a testing dataset from Oklahoma, showing an average location error of about 5 km. The origin time can be determined based on the earthquake's location and the waveforms recorded by stations through waveform shifting and stacking. This method does not require the complex processing steps of traditional seismic approaches, allowing for rapid earthquake location. Additionally, we apply the network model to data recorded in Southern California through transfer learning for further application. The results show that this new method is stable and generalized, making it applicable to earthquake location problems associated with arbitrary station networks. Furthermore, we discuss the effects of data augmentation, network architecture, and the Gaussian radius of labels on the outcomes. These insights help us better understand machine learning algorithms and improve the application of deep learning in earthquake location.

SummaryEvidence of seismic anisotropy is widespread within the Earth, including from individual crystals, rocks, borehole measurements, active-source seismic data, and global seismic data. The seismic anisotropy of a material determines how wave speeds vary as a function of propagation direction and polarization, and it is characterized by density and the elastic map, which relates strain and stress in the material. Associated with the elastic map is a symmetric 6 × 6 matrix, which therefore has 21 parameters. The 21-dimensional space of elastic maps is vast and poses challenges for both theoretical analysis and typical inverse problems. Most estimation approaches using a given set of directional wavespeed measurements assume a high-symmetry approximation, typically either in the form of isotropy (2 parameters), vertical transverse isotropy (radial anisotropy: 5 parameters), or horizontal transverse isotropy (azimuthal anisotropy: 6 parameters). We offer a general approach to explore the space of elastic maps by starting with a given elastic map T. Using a combined minimization and projection procedure, we calculate the closest Σ-maps to T, where Σ is one of the eight elastic symmetry classes: isotropic, cubic, transverse isotropic, trigonal, tetragonal, orthorhombic, monoclinic, and trivial. We apply this approach to 21-parameter elastic maps derived from laboratory measurements of minerals; the measurements include dependencies on pressure, temperature, and composition. We also examine global elasticity models derived from subduction flow modeling. Our approach offers a different perspective on seismic anisotropy and motivates new interpretations, such as for why elasticity varies as a function of pressure, temperature, and composition. The two primary advances of this study are 1) to provide visualization of elastic maps, including along specific pathways through the space of model parameters, and 2) to offer distinct options for reducing the complexity of a given elastic map by providing a higher-symmetry approximation or a lower-anisotropic version. This could contribute to improved imaging and interpretation of Earth structure and dynamics from seismic anisotropy.

SummaryAutomatic event detection and phase picking are critical for processing the large volumes of data produced by modern seismological instrumentation. Accurate picking is especially challenging in Distributed Acoustic Sensing (DAS) recordings, where data quality can significantly vary along segments of the fiber due to localized environmental noise and coupling issues, reducing signal-to-noise ratios (SNR). Similarly, Ocean Bottom Seismometer (OBS) data quality also suffers from these issues. To improve accuracy under diverse conditions, we developed a novel multiband kurtosis-based picking algorithm, Kurtosis-Value-Picker (KVP), that enhances phase picking for both impulsive and emergent seismic signals. Our approach uses characteristic functions (CFs) calculated with sliding windows across multiple frequency bands. Triggers are identified based on localized kurtosis jumps over a few samples, providing greater sensitivity to emergent signals than traditional finite-difference methods. Each individual CF has its own set of triggers, adding flexibility to phase picking and retaining spectral information. We validate the KVP algorithm using earthquake data recorded with DAS on two land and submarine fiber-optic cables, as well as OBS data. We also compare its performance with a widely-cited, kurtosis-based algorithm, the widely-used FilterPicker algorithm, and the well-known PhaseNet model, using impulsive signals on nearby DAS channels as a ground truth for emergent arrivals. Our results demonstrate that KVP provides accurate picks and is suitable for complex seismic datasets.

SummaryAccording to the classical Mohr-Coulomb-Anderson theory, faults form at an angle from the largest regional compressive stress that is approximately 30° for most rocks. However, real settings are more complex and faults often present orientations inconsistent with the angles predicted by the classical theory applied to the present-day regional stress field. The Ridgecrest region hosts a young fault system that is part of the Eastern California Shear Zone, and the 2019 earthquake sequence unveiled orthogonal ruptures at multiple scales, apparently at odds with the classical brittle failure model. We use the Ridgecrest region as a case study and compare surface ruptures that developed during the 2019 earthquake sequence to the expected orientations derived from classical faulting theory and to observations from rock experiments. We focus on the off-fault secondary fractures that developed coseismically at the northern termination of the mainshock fault. We calculate coseismic stress changes from published slip models superimposed to a background stress field. We find that a combination of tectonic regional stresses oriented with the largest compressive stress at N10E–N14E and very weak intensity of coseismic stresses best captures the orientation of off-fault fractures in the classical Mohr-Coulomb-Anderson framework, with an internal rock friction coefficient μ = 0.6. The secondary fractures also show a scale separation: long fractures are most compatible with shear failure, while short fractures cluster along the direction of the largest horizontal stress. The latter is compatible with either local normal faulting or early tensile failures that would later coalesce to form longer faults, consistent with growth of shear fractures in laboratory experiments. Finally, the different orientations of fractures that developed during and prior to the 2019 events suggest that the tectonic stress has rotated over geological timescales. When accounting for the specificity of the area, orthogonal faulting is thus compatible with brittle fracturing with typical experimental values of rock friction coefficient.

SummarySeismic swarms are known to occur in regions with complex deformation and multiple fault systems. The identification of the affected faults, the evolution of the seismicity and the rupture characteristics are key to understand seismotectonics and seismic hazard of such areas. We here address the Atarfe-Santa Fe earthquake swarm with > 5000 events, including six magnitude 4+ earthquakes, recorded in 2021 in the Granada Basin area (S-Spain). We use continuous data from a dense local recording network and apply deep learning models to pick and associate seismic phase arrivals and construct an automatic event catalogue. A double difference relocation of 3196 earthquakes reveals the seismotectonic fine-structure of the swarm. We identify planar, southwest-dipping structures among the relocated hypocentres, consistent with the NW-SE trending, high-angle normal faults present in this sector of the Granada basin. Earthquakes concentrate between 4-7 km depth along three different lineaments. The distinctive pattern of three equidistant, parallel segments allows an association of the swarm with the Ermita los Tres Juanes, Atarfe and Pinos Puente normal faults. These faults outcrop at the upward extrapolation of the swarm, forming an arrangement of three structures that mimic the geometry of the relocated seismicity at depth. In the course of the swarm, the seismicity jumped from the Ermita los Tres Juanes to the Atarfe fault, in midst of a rapid succession of three magnitude 4 earthquakes within 20 minutes, then migrated laterally along both faults, and later migrated basinwards to the Pinos Puente fault, which produced fewer and smaller events. We estimate apparent source time functions for five earthquakes (Mw 4.1 to 4.4) through the deconvolution of empirical Green's functions from the records, suggesting rupture propagation towards NW, N and E directions. An isochrone back projection of apparent source time functions suggests km-scale ruptures with simple slip distributions, showing lateral and updip components of rupture. Our results shed light on the complexity of this seismic swarm in terms of the fault network involved, the propagation of seismicity across the faults and the variable directions of individual ruptures.

SummarySeismic tomography is routinely used to image the Earth’s interior using seismic data. However, in practice, data limitations lead to discretised inversions or the use of regularisations, which complicates tomographic model interpretations. In contrast, Backus-Gilbert inference methods make it possible to infer properties of the true Earth, providing useful insights into the internal structure of our planet. Two related branches of inference methods have been developed – the Subtractive Optimally Localized Averages (SOLA) method and Deterministic Linear Inference (DLI) approaches – each with their own advantages and limitations. In this contribution, we show how the two branches can be combined to derive a new framework for inference, which we refer to as SOLA-DLI. SOLA-DLI retains the advantages of both branches: it enables us to interpret results through the target kernels, rather than the imperfect resolving kernels, while also using the resolving kernels to inform us on trade-offs between physical parameters. We therefore highlight the importance and benefits of a more careful consideration of the target kernels. This also allows us to build families of models, rather than just constraining properties, using these inference methods. We illustrate the advantages of SOLA-DLI using three case studies, assuming error-free data at present. In the first, we illustrate how properties such as different local averages and gradients can be obtained, including associated bounds on these properties and resolution information. Our second case study shows how resolution analysis and trade-offs between physical parameters can be analysed using SOLA-DLI, even when no data values or errors are available. Using our final case study, we demonstrate that SOLA-DLI can be utilised to obtain bounds on the coefficients of basis function expansions, which leads to discretised models with specific advantages compared to classical least-squares solutions. Future work will focus on including data errors in the same framework. This publication is accompanied by a SOLA-DLI software package that allows the interested reader to reproduce our results and to utilise the method for their own research.

SummaryDetermining the precise pattern of crustal deformation enhances our comprehension of crustal deformation traits and the significant earthquakes. By incorporating 21 additional continuous GNSS stations along with existing ones, we generated an updated GNSS velocity field for the northeastern Tibetan Plateau. Using the back-slip dislocation model, we calculated the average slip rates of three major active faults: the Haiyuan fault, the Liupanshan fault, and the Helanshan fault. Our findings indicated that the regional crustal movement does not conform to the equilibrium principles typically associated with the triple junction-like tectonics. This suggests the existence of a newly active tectonic belt within the Longxi block. Consequently, we proposed a revised block model that incorporates a right-lateral shear zone within the Longxi block to account for the observed crustal deformation in the northeastern Tibetan Plateau. Our study indicates that the right-lateral shear zone significantly contributes to the northeastward expansion of the Tibetan Plateau, accounting for approximately 82 per cent of strain accumulation, while the remaining 18 per cent accumulates along the Liupanshan fault. The revised block model emphasizes the pivotal role of the Haiyuan fault and the right-lateral shear belt as the key tectonic factors shaping the crustal deformation pattern. Our result enables a comprehensive understanding of both the spatial variations observed in the GNSS velocity field and the spatial distribution of significant earthquakes in the region.

SummaryFull waveform inversion (FWI) has enjoyed increased attention the past decade, becoming the state of the art for estimating parameters influencing wave propagation in a medium. However, only a few recent emerging efforts have attempted to tackle the challenge of uncertainty quantification in FWI. In this study, we suggest joining FWI with the Bayesian approach, where we provide a post-processing step with an advantageous starting point defined by the global minimum stemming from a deterministic FWI algorithm. Then, using the local ensemble transform Kalman filter (LETKF), we obtain the uncertainty as a follow-up step to the FWI procedure. Within a probabilistic Bayesian inversion framework, the LETKF uses local seismic data to update sets of variables in the subsurface domain. Seismic data for each shot and receiver in the time-domain is in this way matched with subsurface layers, and assimilated in a sequential manner. The methodology is showcased on a realistic model of the Gullfaks field in the North Sea, where we study effects of various seismic acquisition design set-ups, algorithm and model parameter settings. We investigate how these acquisition designs and parameters influence the uncertainty reduction and bias of the inversion results. We highlight the importance of studying statistical performance metrics to ensure a balance between bias and underestimation of uncertainty.

SummaryModel-based information about the global water cycle, in particular the redistribution of terrestrial water masses, is highly relevant for the understanding of Earth system dynamics. In many geodetic applications, hydrological model results play an important role by augmenting observations with a higher spatio-temporal resolution and gapless coverage. Here we demonstrate the feasibility of the high-resolution, open-source hydrological model OS LISFLOOD to simulate terrestrial water storage (TWS) variations with a spatial sampling of up to about 5 km (0.05○). Validation against data from satellite gravimetry reveals that the choice of the maximum soil depth has a significant impact on long-term trends in TWS, mainly in the deepest soil layer. We find that refining the soil depth definition effectively reduces spurious TWS trends, while preserving accuracy in modeled river discharge. Using the modified model set-up, we show that in many regions TWS from OS LISFLOOD fits better to observations than TWS from the Land Surface Discharge Model (LSDM) routinely operated at the GFZ and used in geodetic applications worldwide. The advantage of the high spatial resolution of the OS LISFLOOD implementation is shown by comparing vertical surface displacements to GNSS observations in a global network of stations. The data set presented here is the first application of OS LISFLOOD to generate quasi-global (regions south of 60○S excluded) daily 0.05○ TWS fields for a 23-year period (2000–2022).

SummaryAnalysis of the 48-year of satellite laser ranging (SLR) data of multiple satellites shows that long-term variations in the Earth's dynamical oblateness represented by the second-degree zonal harmonic J2 is best characterized by the superposition of a quadratic trend, 10.5-year and 18.6-year variations. These variations result from climate-related mass changes, tides, and core flow-induced variations at the core-mantle boundary (CMB). We determined that the global ocean's response to the lunar attraction at the 18.6-year period is near equilibrium, with an amplitude of 0.4735 ± 0.008 cm and an error of ∼11% relative to the modeled amplitude (0.4224 cm), and ∼2±3 degrees of phase lag. The 18.6-year Love Number was found to be 0.01375-i0.00553 with an error of 2% for both the real and imaginary parts and a phase correction π for the imaginary part of the IERS 2010 anelasticity model. The nominal frequency-independent anelasticity Love number, k₂, was determined to be 0.3022 ± 0.0001 for the 18.6-year period, based on a reference frequency of 200 seconds and α = 0.1514 for mantle anelasticity for mantle anelasticity. This study also reveals a significant gravitational signal (3.06×10−11) in J2 obstructs the Earth's mantle anelastic response to the 18.6-year tidal forces, reflecting in phase shift of π the imaginary part of the IERS 2010 Love number. This signal can be characterized by a positive Love number of 0.01106 in the modeling of the variation in J2 coupling with the 18.6-year tide. This signal is possibly produced by the core dynamics, which creates a gravitational signal in J2 with an amplitude of 3.36×10−11 at the decadal time scale and could account for ∼70% of the observed 10.5-year variation.

SummaryTo clarify the 3-D crustal and upper mantle structure of the Bolivian orocline in the Central Andes, we conduct azimuthal anisotropy tomography using newly measured teleseismic fundamental mode Rayleigh-wave phase and amplitude data at periods of 25-110 s. Our tomography shows that the subducting Nazca slab is imaged as a high-velocity zone beneath the study region, except for areas where the Nazca ridge is subducting. Azimuthal anisotropy in the subducting slab generally exhibits trench-parallel fast-velocity directions (FVDs), but it becomes complex in and around the subducting Nazca ridge. Low-velocity anomalies with trench-normal FVDs exist in the mantle wedge beneath active arc volcanoes and backarc regions beneath Altiplano. A significant high-velocity zone with relative weak anisotropy exists in the crust of the overriding plate above the Peruvian flat slab in the study region. In contrast, low-velocity anomalies with trench-parallel FVDs are revealed in the crust beneath Altiplano. Furthermore, a high-velocity zone with depth-varying FVDs appears beneath Eastern Cordillera and its surrounding regions, which may indicate the westward underthrusting cratonic lithosphere. These tomographic features well capture the primary 3-D structure of the middle-lower crust and upper mantle beneath the Bolivian orocline, which results from the subduction of an oceanic lithosphere and the delamination and underthrust of a continental lithosphere, leading to the second-highest plateau on Earth.

SummaryIn this study, we applied the “in-situ Vp/Vs method” to monitor variations of seismic velocity ratio (Vp/Vs) within swarms, providing insights into eruption processes. This method, particularly effective in volcanic regions, estimates Vp/Vs by comparing P- and S-wave arrival times of closely located earthquake pairs, reducing errors from unknown crustal velocity variations and is well-suited for detecting rapid changes associated with volcanic swarms. Our study focused on seismic swarms on the Reykjanes Peninsula, south-west Iceland where, swarms have been frequent since 2017 and led up to eruptions in 2021, 2022, and 2023. We analyzed the entire period (2017–2023) as well as the 2021 swarm separately using data from over 40,000 seismic events recorded by the REYKJANET network. We observed significant decrease in the Vp/Vs ratio before major pre-eruption swarms, compared to the background Vp/Vs value of 1.78. From the 2020 swarm, we observed a lower Vp/Vs of 1.72, but the lowest estimated value was 1.70, associated with the 2021 pre-eruption swarm that preceded Fagradalsfjall's first eruption after 7000 years. Reduced Vp/Vs ratios were also noted before the 2022 and 2023 eruptions, suggesting supercritical fluids in the crust during these stages. We also introduce the concept of “change points” to interpret Vp/Vs variations along the dyke. Change points denote specific locations or times of significant Vp/Vs shifts, potentially indicating subsurface changes such as fluid influx or new fracturing from magma intrusion. Identifying these points allows us to pinpoint key moments when the system undergoes substantial changes, offering insights into eruption timing and location. Focusing on 2021 pre-eruption swarm, interestingly the spatial change point found in a location very close to the eruption site. Temporal analysis identified two main change points: the first corresponding with initial activity in the northern dyke and the second with a shift to the southern segment, ultimately leading to eruption. These points mark stages in magma progression, with each showing an initial rapid Vp/Vs drop that could indicate CO₂-rich fluid infiltration, followed by an increase as magma enters. The in-situ Vp/Vs method's sensitivity to changes in seismic properties makes it a powerful tool for real-time volcanic monitoring. By detecting critical Vp/Vs changes with minimal computational demand, this method has potential for integration with online seismic networks, providing an effective early warning system for volcanic hazards.

SummarySeismic source models that use an elastic relation between pressure decrease, compaction, and stress change have been shown to successfully reproduce induced seismicity in producing natural gas reservoirs undergoing differential compaction. However, this elastic relation is inconsistent with observations of non-linear reservoir compaction in the Groningen field. We utilize critical state mechanics theory to derive a 3D stress-strain framework that is able to house 1D non-linear stress-strain relations typically used for subsidence models, without the need for recalibration of the subsidence model parameters. This is used to adapt the elastic thin sheet stress model that is currently in use as the state-of-the-art for seismicity predictions as part of the hazard and risk assessment of the Groningen gas field. The new thin sheet model has one additional model parameter that modulates the impact of inelastic deformation on fault loading, whilst keeping the intended function of the model calibration from the original elastic thin sheet model intact. The resulting elastic-viscoplastic thin sheet stress model is consistent with previously reported non-linear rate-dependent reservoir compaction in Groningen found from inverting subsidence data and from rock deformation experiments. Our elastic-viscoplastic thin sheet stress model is able to predict ongoing stress increase, and therefore ongoing seismicity, in areas where pressure does not decrease anymore due to shut-in. A pseudo-prospective forecasting exercise indeed shows that the elastic-viscoplastic stress model performs better than the linear elastic stress model. This model addition ensures that the Groningen seismic source model is well suited for predicting seismicity in the post shut-in phase.

SummaryThis paper presents a novel framework for full-waveform seismic source inversion using simulation-based inference (SBI). Traditional probabilistic approaches often rely on simplifying assumptions about data errors, which we show can lead to inaccurate uncertainty quantification. SBI addresses this limitation by learning an empirical probabilistic relationship between the parameters and data, without making assumptions about the data errors. This is achieved through the use of specialised machine learning models, known as neural density estimators, which can then be integrated into the Bayesian inference framework. We apply the SBI framework to point-source moment tensor inversions as well as joint moment tensor and time-location inversions. We construct a range of synthetic examples to explore the quality of the SBI solutions, as well as to compare the SBI results with standard Gaussian likelihood-based Bayesian inversions. We then demonstrate that under real seismic noise, common Gaussian likelihood assumptions for treating full-waveform data yield overconfident posterior distributions that underestimate the moment tensor component uncertainties by up to a factor of 3. We contrast this with SBI, which produces well-calibrated posteriors that generally agree with the true seismic source parameters, and offers an order-of-magnitude reduction in the number of simulations required to perform inference compared to standard Monte Carlo sampling techniques. Finally, we apply our methodology to a pair of moderate magnitude earthquakes in the North Atlantic. We utilise seismic waveforms recorded by the recent UPFLOW ocean bottom seismometer array as well as by regional land stations in the Azores, comparing full moment tensor and source-time location posteriors between SBI and a Gaussian likelihood approach. We find that our adaptation of SBI can be directly applied to real earthquake sources to efficiently produce high quality posterior distributions that significantly improve upon Gaussian likelihood approaches.

SummaryThe thermal structure of the continental crust plays a critical role in understanding its elastic and rheologic properties as well as its dynamic processes. Thermal parameter datasets on continental scales have been used to constrain the crustal thermal structure, including both the direct (e.g., temperature, heat flux, and heat conductivity measured at the surface) and indirect (e.g., seismically derived Mohorovičić discontinuity (Moho) temperature, geomagnetically derived Curie depth) observations. In this study, we present a new continental scale crustal heat generation model with additional information from seismologically-inferred crustal composition. Together with previous direct and indirect thermal parameter datasets in the conterminous United States, we use the new crustal heat generation model to construct a 3-dimensional (3-D) crustal temperature model under a newly developed Bayesian framework. Specifically, we first derive profiles of crustal heat generation based on an empirical geochemical relationship at 1683 locations where seismologically derived crustal composition information is available. Then for each of these locations, the average heat generation values in the upper, middle, and lower crust are combined with other thermal parameters through a Markov Chain Monte-Carlo inversion for a conductive, vertically smooth temperature profile. The results, posterior distributions of temperature profiles, are used to generate a 3-D crustal thermal model with the uncertainties systematically assessed. The new temperature model overall exhibits similar patterns to that from the U.S. Geological Survey National Crustal Model, but also reduces possible biases and the model's dependence on a single thermal parameter.

SummaryThe ∼300km-long rupture of the February 6 2023 Kahramanmaraş earthquake began in the Narlı section of the Karasu trough, a pull-apart basin sandwiched and sheared between the two major strike-slip faults of the region, the East Anatolian Fault (EAF) on the west and the Dead Sea Fault (DSF) on the east. Rupture started where the northern segment of the DSF enters the Narlı Basin with Mw7.0 sub-event and propagated across the basin before making its junction with the EAF. In the seven months preceding the earthquake this basin was the seat of anomalous seismic activity. This activity occurred in bursts interweaved with periods of quiescence. It started near-concomitantly in two clusters located on the opposite edges of the pull-apart basin ∼20 km apart. The organization of this seismicity into families of numerous repeating earthquakes suggests an aseismic process linked to fault healing and rapid reloading in a critically stressed zone. By December 2022, two months before the earthquake, activity had migrated to a cluster located along the path that rupture was to follow during the initial stage of the earthquake. These observations show that the pull-apart basin where rupture started was progressively deforming in a succession of bursts before the earthquake. This regional-scale deformation is closely linked with the transitional nature of geodynamics and kinematics influenced by large-scale fault interactions in the surrounding area. The location of the epicenter near the northern termination of the rupture of the 1822 M7.4 earthquake suggests that the ∼45 km long Narlı sub-rupture which constituted the first stage of this giant earthquake was closing a long-present seismic gap between the DSF and the EAF.

SummaryOcean-bottom seismic acquisitions are gaining widespread popularity across a variety of subsurface applications. However, the high cost of these systems often necessitates receiver geometries with large intervals between ocean-bottom cables or nodes. The upside-down Rayleigh-Marchenko (UD-RM) method has been recently proposed as an effective solution for accurate redatuming and imaging of sparse seabed data. In this paper, we present the first successful application of the UD-RM method to field data. We demonstrate that in the presence of a shallow seabed, an improved data pre-processing workflow is necessary to generate more accurate input wavefields compared to the one produced by the workflow presented in the original paper. To validate the proposed processing workflow, the UD-RM method is initially tested on a synthetic dataset that mimics the Volve field data (referred to as the Volve synthetic dataset); this is followed by its application to a 2D line of the Volve ocean-bottom cable dataset. Subsequently, the field dataset is subsampled by retaining only 25 percent of the total receivers to numerically validate the UD-RM method’s capability to handle sparse receiver arrays. The resulting images reveal that the UD-RM method, when paired with our enhanced data processing workflow, can effectively handle surface-related multiples, internal multiples, and sparse receiver arrays, producing accurate imaging results without the need for costly and labor-intensive multiple removal processes. Finally, we provide theoretical insights and numerical evidence supporting the necessity of source-side deghosting prior to redatuming. While a pre-processing workflow that omits source-side deghosting can offer some practical advantages, we show that this ultimately produces blurrier images compared to those obtained using source-side deghosted input data.