Geophysical Journal International

SummaryFault Zone Head Waves (FZHW) are a key diagnostic tool to identify bimaterial interfaces along fault zones. We detect and analyze FZHW recorded in the waveforms from the local MONGAN (MONitoring of the GANos Fault) seismic network along the Ganos section of the North Anatolian Fault Zone, northwestern Türkiye, between October 2017 and July 2019. MONGAN covers the Ganos fault with different inter-station distances ranging from 25 m to ∼4 km. To detect FZHW, an automatic detector is used as a preliminary analysis method followed by manual revision and particle-motion analyses to distinguish between FZHW and direct P waves. FZHWs are predominantly detected at the southern side of the fault. The observed FZHWs have a moveout (∆t) with respect to the direct P arrivals, increasing with distance traveled along the fault and indicating a deep bimaterial interface down to the bottom of the seismogenic crust. The average velocity contrast is estimated to be 5.9% across the fault. Near fault-recordings indicate that the Ganos Fault is offset by ∼250 m with respect to the surface trace obtained from literature. To a lesser extent, FZHW are also observed in the northern stations from the fault, indicating a shallow wedge-shaped low-velocity portion constituted by highly fractured material to either side along the southwestern section of the Ganos Fault between the fast Eocene block to the north and the slow Miocene block to the south. The seismic velocity contrast and geological complexity have important implications for the rupture evolution during future earthquakes on the Ganos fault in that they would progress predominantly westward, away from Istanbul and Tekirdağ. Furthermore, an asymmetric aftershock distribution skewed to the northern block can be expected, with subsequent implications for site-dependent risk there. Our results allow to revise focal mechanism solutions by separating FZHW from direct-P wave for previous Sea of Marmara earthquakes.

SummaryBorehole breakout (BO) has increasingly been utilised to estimate in-situ stress magnitudes given the importance of the stress field in subsurface activities and the limitations of conventional stress measurement techniques. In this study, a new backpropagation neural network model is developed to estimate both maximum and minimum horizontal stress magnitudes from multi-scale BO data. A total of 150 experimental data points from pre-stressed true-triaxial laboratory tests and 44 field data from a mine site in Australia and the literature are collected and employed for model development and validation. Compared to previous studies, the collected dataset is significantly enhanced in both quantity and quality. To address discrepancies in stress magnitudes between experimental and field data, the three principal stresses are normalised by borehole wall strength (BWS). Overall, the model achieves mean absolute percentage errors of below 8% for the maximum horizontal stress and below 20% for the minimum horizontal stress, significantly outperforming the previous model developed for this purpose. Furthermore, these error rates fall within the typical error range (10-20%) of conventional stress measurement techniques, indicating the model's sufficient accuracy for practical applications. Moreover, the effectiveness and generalisability of the model are verified using 166 additional BOs from two mine sites, which are independent of those used in model development. Continuous and detailed stress profiles are established based on these BOs, covering greater depth intervals than the stress measurements from the overcoring method. The results of this study demonstrate that the proposed model can provide reliable and accurate stress estimation, utilising input parameters that can be readily obtained from borehole geophysical logs.

SummaryThe behavior of uniaxial single domain magnetite particles in rock and paleomagnetic experiments was first described in the 1940s by Néel and Stoner and Wohlfarth. Since this time, micromagnetism has allowed us to gain a better understanding of magnetic particles in the single vortex or multi-domain states. By contrast, when describing the behavior of assemblages of single domain particles, simplifying assumptions made in the 1940s are still used today. In particular, most rock and paleomagnetic simulations involve magnetite with a uniaxial anisotropy. These assumptions are not necessary in the modern day, as data on other magnetic minerals has been collected, and modern computers are powerful enough to easily calculate the behavior of multiaxial particles. We present a new software package called the Single Domain Comprehensive Calculator (SDCC). This package can simulate a large number of thermally activated rock and paleomagnetic experiments with distributions of single domain particles. These include acquisition of viscous remanence, thermal demagnetization experiments, hysteresis loops, and paleointensity protocols. The package provides a simple Python scripting interface for users to define custom experiments and run models on a laptop computer. Preliminary simulations run with the SDCC demonstrate that magnetocrystalline anisotropy can have a significant effect on the thermoviscous behavior of single domain particles, despite normally being ignored in models. This highlights a need for further investigation into the behavior of single domain particles.

SummaryIn this study, a modified two-dimensional gravimetric inversion algorithm is presented that is based on the reversible-jump Markov chain Monte Carlo (RJMCMC) method with the Talwani equation. To ensure the validity of the Talwani equation and accurate gravity anomaly calculations, the Shamos-Hoey algorithm is incorporated as an additional acceptance condition to prevent intersections in the model polygon. This improves upon the method proposed by Luo by iteratively refining a polygon model based on gravity anomalies while maintaining physical validity. Additionally, we suggest a revision of the prior density function to better test the proposal models. This method estimates the shape and location of subsurface intrusions, providing valuable insights into subsurface geological structures. This positions the algorithm as a valuable tool for geophysical research.

SummaryNear-field large-amplitude seismograms are essential for the rapid inversion of earthquake source parameters using waveform inversion methods such as the Cut And Paste (CAP) method for disaster assessment. High-rate Global Navigation Satellite System (GNSS) relative positioning (RP) provides precise, rapid, and real-time measurement of near-field large-amplitude displacements. However, RP records motion with respect to a reference station, and the reference station's movements become part of the relative displacement waveforms. Therefore, seismic source parameter estimates may be inaccurate if the reference station's motion is not taken into consideration, and doing so affects some basic assumptions made in the CAP method. To overcome this problem, we develop an expanded differential CAP inversion approach specifically for high-rate GNSS RP (CAP-RP) that accounts for the motion of the reference station. Two methods are proposed to implement CAP-RP: an expanded differential CAP (D-CAP) and an iterative post-processing CAP (P-CAP). We assess the performance of CAP-RP with different datasets, using the July 2019 Mw 6.4 earthquake in California as a case study. Both CAP-RP techniques produce accurate source parameters in synthetic data inversion tests, indicating the feasibility of the strategy. However, P-CAP is more time-efficient than D-CAP, making it the better option. Generally, results from high-rate GNSS RP, broadband seismographs, and their inverted combinations exhibit consistency in observational data inversion testing. Our results also demonstrate that more accurate source parameters can be obtained by combining sensitive far-field broadband seismograph data with large amplitude near-field GNSS RP waveforms.

SummarySpectral induced polarization (SIP) has been suggested as a non-invasive and cost-effective tool to detect and monitor aromatic rich organic matter such as biochar. In our study, we show that SIP can track biochar concentration up to 10% (wt.) in a soil with a clay content of 20%. Assessment of changes in the concentration of biochar was conducted according to double Pelton parameters and the maximal phase determined at 11.7 Hz, a frequency at which a polarization peak is observed in the presence of biochar. All SIP-derived parameters were correlated with the biochar content, with the exception of the relaxation time of the polarization peak occurring at 11.7 Hz, which was related to soil water saturation in previous investigations. Among studied parameters, the phase value that we measured at 11.7 Hz may therefore consist in a simple and reliable methodology to evaluate the biochar content on SIP in our experiment. Several steps are still necessary before a widespread field application notably by considering how modifications in the chemistry of biochar with time can interact with biochar concentration and water saturation to modify polarization processes shaping SIP curves. Beyond the scope of tracking changes in the content of highly aromatic OM – such as biochar - in soils, this study suggests that the degree of aromaticity of OM can play a key role in the SIP response paving the way for wider use of SIP in soil science.

SummaryLong-period seismic events (LPs) are observed within active volcanoes, hydrothermal systems, and hydraulic fracturing. The prevailing model for LP seismic events suggests that they result from pressure disturbances in fluid-filled cracks that generate slow, dispersive waves known as Krauklis waves. These waves oscillate within the crack, causing it to act as a seismic resonator whose far-field radiations are known as LP events. Since these events are generated from fluid-filled cracks, they have been used to analyze fluid transport and fracturing in geological settings. Additionally, they are deemed precursors to volcanic eruptions. However, other mechanisms have been proposed to explain LP seismicity. Thus, a robust interpretation of these events requires understanding all parameters contributing to LP seismicity. To achieve this, for the first time, we have developed a physical model to investigate LP seismicity under controlled-source conditions. The physical model consists of a 30 cm × 15 cm × 0.2 cm crack embedded within a concrete slab with dimensions of 3 m × 3 m × 0.24 m. Using this apparatus, we investigate fundamental factors affecting long-period seismic signals, including crack stiffness, fluid density and viscosity, radiation patterns, and triggering location. Our findings are consistent with the theoretical model for Krauklis waves within a fluid-filled crack.In this study, we examine the interplay between fluid properties and characteristics of waves within and radiated from the crack model. Records from a pressure transducer within the crack model have the same frequency characteristics as the surface sensors, indicating that the surface sensors are recording the crack waves. Because the crack stiffness parameters for all the fluids are relatively high, fluid density variations have a larger effect on the crack wave frequency, with higher density fluids yielding lower resonance frequencies. Similarly, the quality factor (Q) decreases with increasing fluid density. We also find that an increase in fluid viscosity along with the increased fluid density results in a decrease in resonance frequency and Q. Trigger locations at the middle of the crack length and width most effectively resonated the first and second transverse modes. Thus, this physical model can offer new horizons in understanding LP seismicity and bridge the gap between theoretical models and observed LP signals.

SummaryMany geophysical studies require knowledge on the present-day temperature distribution in Earth’s mantle. One example are geodynamic inverse models, which utilize data assimilation techniques to reconstruct mantle flow back in time. The thermal state of the mantle can be estimated from seismic velocity perturbations imaged by tomography with the help of thermodynamic models of mantle mineralogy. Unique interpretations of the tomographically imaged seismic heterogeneity can either be obtained by incorporating additional data sets or requires assumptions on the chemical composition of the mantle. However, even in the case of (assumed) known chemical composition, both the seismic and the mineralogical information are significantly affected by inherent limitations and different sources of uncertainty. Here, we investigate the theoretical ability to estimate the thermal state of the mantle from tomographic models in a synthetic closed-loop experiment. The ‘true’ temperature distribution of the mantle is taken from a 3-D mantle circulation model with Earth-like convective vigour. We aim to recover this reference model after: 1) mineralogical mapping from the ‘true’ temperatures to seismic velocities, 2) application of a tomographic filter to mimic the effect of limited seismic resolution, and 3) mapping of the ‘imaged’ seismic velocities back to temperatures. We test and quantify the interplay of tomographically damped and blurred seismic heterogeneity in combination with different approximations for the mineralogical ‘inverse’ conversion from seismic velocities to temperature. Owing to imperfect knowledge of the parameters governing mineral anelasticity, we additionally investigate the effects of over- or underestimating the corresponding correction to the underlying mineralogical model. Our results highlight that, given the current limitations of seismic tomography and the incomplete knowledge of mantle mineralogy, magnitudes and spatial scales of a temperature field obtained from global seismic models deviate significantly from the true state, even in the idealized case of known bulk chemical composition. The average deviations from the reference model are on the order of 50–100 K in the upper mantle and – depending on the resolving capabilities of the respective tomography – can increase with depth throughout the lower mantle to values of up to 200 K close to the core-mantle boundary. Furthermore, large systematic errors exist in the vicinity of phase transitions due to the associated mineralogical complexities. When used to constrain buoyancy forces in time-dependent geodynamic simulations, errors in the temperature field might grow non-linearly due to the chaotic nature of mantle flow. This could be particularly problematic in combination with advanced implementations of compressibility, in which densities are extracted from thermodynamic mineralogical models with temperature-dependent phase assemblages. Erroneous temperatures in this case might activate ‘wrong’ phase transitions and potentially flip the sign of the associated Clapeyron slopes, thereby considerably altering the model evolution. Additional testing is required to evaluate the behaviour of different compressibility formulations in geodynamic inverse problems. Overall, the strategy to estimate the present-day thermodynamic state of the mantle must be selected carefully to minimize the influence of the collective set of uncertainties.

SummaryGeomagnetic observations at Eskdalemuir observatory in Southern Scotland reveal reduced amplitudes in the vertical component variations compared with the horizontal components for periods of less than an hour. A subsurface high conductivity feature has previously been suggested to account for this anomaly. However, past studies have overlooked the effect of seasonal source changes and impact of solar activity on external geomagnetic field variations. The vertical magnetic transfer function - referred to as the tipper - relates temporal variation in the vertical magnetic field to those in the horizontal magnetic field and is sensitive to lateral electrical conductivity contrasts in the subsurface. Quantifying the seasonal variations in the tipper helps to identify times when external field variations minimally bias tipper estimates, thereby providing a more accurate representation of subsurface conductivity. Ionospheric current systems, particularly during geomagnetic storms, may violate the plane wave assumption underlying tipper estimation at mid-latitudes. This may allude to a more complex source geometry responsible for magnetic field variations. Our study quantifies and proposes a correction for space weather-driven external field contributions to observations for periods shorter than one hour. Using high-quality digital magnetic field data with a one-minute sampling rate from 2001 to 2019, we estimate the tipper at Eskdalemuir, revealing seasonal differences that increase with periods between 1000 s and 10000 s. After finding that tipper estimates during the 2016 time series are least affected by seasonal effects, we used one-second time series and a simple empirical model to quantify the daily variability of the tipper. The model consists of annual and semi-annual terms plus a term proportional to either the F10.7 cm solar flux or geomagnetic Ap index. Neither model fits the data to within the expected error, but the model that uses Ap has better fit. Tipper estimates from temporary site deployments are affected by these seasonal external variations, and we correct those obtained at sites near Eskdalemuir during a recent field experiment using this model.

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.