Search

Many scientific investigations require that the values of a set of model parameters are estimated using recorded data. In Bayesian inference, information from both observed data and prior knowledge is combined to update model parameters probabilistically by calculating the posterior probability distribution function. Prior information is often described by a prior probability distribution. Situations arise in which we wish to change prior information during the course of a scientific project. However, estimating the solution to any single Bayesian inference problem is often computationally costly, as it typically requires many model samples to be drawn, and the data set that would have been recorded if each sample was true must be simulated. Recalculating the Bayesian inference solution every time prior information changes can therefore be extremely expensive. We develop a mathematical formulation that allows the prior information that is embedded within a solution, to be changed using variational methods, without recalculating the original Bayesian inference. In this method, existing prior information is removed from a previously obtained posterior distribution and is replaced by new prior information. We therefore call the methodology variational prior replacement (VPR). We demonstrate VPR using a 2-D seismic full waveform inversion example, in which VPR provides similar posterior solutions to those obtained by solving independent inference problems using different prior distributions. The former can be completed within minutes on a laptop computer, whereas the latter requires days of computations using high-performance computing resources. We demonstrate the value of the method by comparing the posterior solutions obtained using three different types of prior information: uniform, smoothing and geological prior distributions. [ABSTRACT FROM AUTHOR]

Geophysicists today face the challenge of quickly and reliably interpreting extensive controlled-source electromagnetic (CSEM) data sets to map subsurface conductivity structures within realistic geological environments. An ideal 3-D CSEM inversion algorithm using tetrahedral grids should be capable of distinguishing different resolution requirements between forward modelling and inversion grids, have an optimal parallel strategy that fully exploits the inherent independence of CSEM data sets while also possessing the capability to handle large-scale geo-electrical models, and incorporate conductivity anisotropy which should be a common characteristic in realistic subsurface environments. However, existing tools in the geo-electromagnetic community often fall short of these three demands. Addressing this gap, our study introduces a scalable and parallel anisotropic inversion technique for CSEM data, capitalizing on the potential of unstructured tetrahedral grids. We first apply the tetrahedral longest-edge bisection method to create a refined dense, heterogeneous forward modelling grid from a coarse inversion grid. This refinement, focused on areas around transmitters and receivers, is seamlessly integrated within the coarser inversion grid's topology, enabling precise conductivity mapping and preserving electromagnetic response accuracy during model updates. We further innovate with a source-mesh double-level parallel strategy, utilizing the message passing interface technique for parallel handling of independent CSEM data sets and large-scale geo-electrical models. Externally, we dedicate a processor for inversion model updates employing the Limited-memory Broyden–Fletcher–Goldfarb–Shanno optimization algorithm and divide other processors into groups, each associated with specific transmitting sources and frequencies. Internally, in each group, we employ a domain-decomposition-based scalable and robust iterative solvers using the Auxiliary-Space Maxwell pre-conditioner to parallel quickly calculate the electromagnetic responses from its assigned source-frequency set. Additionally, recognizing the potential for electrical conductivity anisotropy in field data, we incorporate the case of vertical transverse isotropy. We validate the effectiveness of our method through examples, including an isotropic land model with undulating topography, an anisotropic marine model and a real-field data case. Results from both synthetic and field data inversions underscore our method's significant advancements in efficiency and practicality, particularly in addressing large-scale 3-D CSEM data sets inversion challenges in realistic geological environments. [ABSTRACT FROM AUTHOR]

We present a general framework for multiphysics joint inversion of any number of geophysical data sets. Its main feature is the use of the variable splitting approach: an auxiliary multiparameter model space is introduced in which minimization of the coupling and stabilizing functionals is carried out. The use of rediscretization and interpolation to map between this auxiliary space and the model spaces allows the coupled models to have completely different parametrizations. Joint inversion is decoupled into the individual inversion and the coupling-regularization subproblems, each of which can be solved by a different optimization algorithm. For each subproblem, the linking term controlling the distance between the model and the corresponding auxiliary variable takes the form of a quadratic regularization with a reference model. As a result, any existing inversion code supporting such regularization can be integrated without modifications into the developed framework. As a concrete example scheme, we consider an application of the framework to 3-D joint inversion of magnetotelluric, seismic refraction and gravity data. We discuss different coupling functionals, mainly those corresponding to the more universal structural constraints: joint total variation, joint minimum support, cross-gradient, one-way cross-gradient and their combinations for a general multimodel case. The use of coupling based on explicit 'petrophysical' relationship between the properties is also considered. Performance of the developed framework is studied on three synthetic cases: a time-lapse joint inversion of full-tensor gravity gradiometry and seismic data, a joint inversion of magnetotelluric, seismic and gravity data and a joint inversion for electrical resistivity tomography and audio-magnetotellurics. [ABSTRACT FROM AUTHOR]

The Korean infrasound catalogue (KIC) covers 1999–2022 and characterizes a rich variety of source types as well as document the effects of the time-varying atmosphere on event detection and location across the Korean Peninsula. The KIC is produced using data from six South Korean infrasound arrays that are cooperatively operated by Southern Methodist University and Korea Institute of Geoscience and Mineral Resources. Signal detection relies on an Adaptive F-Detector that estimates arrival time and backazimuth, which draws a distinction between detection and parameter estimation. Detections and associated parameters are input into a Bayesian Infrasonic Source Location procedure. The resulting KIC contains 38   455 infrasound events and documents repeated events from several locations. The catalogue includes many anthropogenic sources such as an industrial chemical explosion, explosions at limestone open-pit mines and quarries, North Korean underground nuclear explosions and other atmospheric or underwater events of unknown origin. Most events in the KIC occur during working hours and days, suggesting a dominance of human-related signals. The expansion of infrasound arrays over the years in South Korea and the inclusion of data from the International Monitoring System infrasound stations in Russia and Japan increase the number of infrasound events and improve location accuracy because of the increase in azimuthal station coverage. A review of selected events and associated signals at multiple arrays provides a location quality assessment. We quantify infrasound events that have accompanying seismic arrivals (seismoacoustic events) to support the source type assessment. Ray tracing using the Ground-to-Space (G2S) atmospheric model generally predicts observed arrivals when strong stratospheric winds exist, although the predicted arrival times have significant discrepancies. In some cases, local atmospheric data better captures small-scale variations in the wind velocity of the shallow atmosphere and can improve arrival time predictions that are not well matched by the G2S model. The analysis of selected events also illustrates the importance of topographic effects on tropospheric infrasound propagation at local distances. The KIC is the first infrasound catalogue compiled in this region, and it can serve as a valuable data set in developing more robust infrasound source localization and characterization methods. [ABSTRACT FROM AUTHOR]

Inverse problems are typically tackled using deterministic optimization methods that may become trapped in a local minimum or probabilistic methods that can be computationally demanding. In this study, we explore the potential of the back propagation neural network (BPNN) optimized by the genetic algorithm (GA) for onshore transient electromagnetic (TEM) inversion. The GA is employed to optimize the initial parameters of the BPNN, enhancing its global optimization ability. Once the BPNN optimized by GA (GA-BPNN) is properly trained, it can provide the distribution of subsurface electrical conductivity (σ) in 0.1 s. We train the GA-BPNN using synthetic datasets generated by TEM forward modeling and assess its reliability using both synthetic and field data. Theoretical simulations demonstrate that compared with BPNN, the error of GA-BPNN on the inversion results of six samples is reduced by 23.2%. Furthermore, this method can provide reliable results even in the presence of noise in the TEM response. Finally, applying this inversion method to karst exploration with measured data proves the reliability and robustness of the proposed method. The proposed method can support quasi-real-time imaging of subsurface structures and provides a powerful tool for the interpretation of field TEM data. [ABSTRACT FROM AUTHOR]

Three-dimensional gravity inversion is a process of obtaining the location, shape, and physical property parameters of underground anomaly sources using gravity anomaly data observed on the surface. In recent years, with the rapid development of data-driven methods, the application of deep learning (DL) to 3D gravity inversion has also attracted wide attention and achieved certain results. In this paper, based on the U-Net network, a three-dimensional gravity inversion method using an attention feature fusion mechanism is proposed. Using U-Net as the basic framework, the coarse-grained semantic features and fine-grained semantic features in the encoder and decoder are connected by long hops, and the global and local semantic features are aggregated through the attention feature fusion module, which avoids feature loss in the network training process. Compared with the inversion results of the U-Net network, the proposed method has a higher vertical resolution and effectively alleviates the influence of the skin effect on three-dimensional gravity inversion. Ablation experiments show that the attention feature fusion module is the key to improving the vertical resolution and prediction accuracy of inversion results. Noise experiments show that the inversion network in this study has a strong anti-noise ability and good generalization performance. The experimental results of the inversion network used in the prediction of the SAN Nicolas deposit in Mexico show that the inversion network can clearly predict the basic location and general shape of the sulfur deposit, and the results are in good agreement with the known geological data. [ABSTRACT FROM AUTHOR]

SUMMARY: Inverse problems occur in many fields of geophysics, wherein surface observations are used to infer the internal structure of the Earth. Given the non-linearity and non-uniqueness inherent in these problems, a standard strategy is to incorporate a priori information regarding the unknown model. Sometimes a solution is obtained by imposing that the inverted model remains close to a reference model and with smooth lateral variations (e.g. a correlation length or a minimal wavelength are imposed). This approach forbids the presence of strong gradients or discontinuities in the recovered model. Admittedly, discontinuities, such as interfaces between layers, or shapes of geological provinces or of geological objects such as slabs can be a priori imposed or even suggested by the data themselves. This is however limited to a small set of possible constraints. For example, it would be very challenging and computationally expensive to perform a tomographic inversion where the subducting slabs would have possible top discontinuities with unknown shapes. The problem seems formidable because one cannot even imagine how to sample the prior space: is each specific slab continuous or broken into different portions having their own interfaces? No continuous set of parameters seems to describe all the possible interfaces that we could consider. To circumvent these questions, we propose to train a Generative Adversarial neural Network (GAN) to generate models from a geologically plausible prior distribution obtained from geodynamic simulations. In a Bayesian framework, a Markov chain Monte Carlo algorithm is used to sample the low-dimensional model space depicting the ensemble of potential geological models. This enables the integration of intricate a priori information, parametrized within a low-dimensional model space conducive to efficient sampling. The application of this approach is demonstrated in the context of a downscaling problem, where the objective is to infer small-scale geological structures from a smooth seismic tomographic image. [ABSTRACT FROM AUTHOR]

SUMMARY: We present the first 3-D upper-mantle conductivity models obtained by an inversion of the satellite-derived tidally induced magnetic fields (TIMFs). We primarily use the M $_2$ period, but the potential benefit of the O $_1$ period is also inspected. The inverse-problem solution is found using the recently developed frequency-domain, spherical harmonic finite-element method based on the adjoint approach. We tested two different TIMF data sets derived from the satellite measurements of the Swarm mission and two different regularizations; the solution is either required to be sufficiently smooth or reasonably close to the a priori 3-D conductivity model WINTERC-e Wd-emax. The reconstructed conductivity models are locally compared with the 1-D conductivity profiles from other studies. If we use one of the available TIMF data sets, the smooth reconstructed model gravitates towards Wd-emax and the TIMF-adjusted Wd-emax model is closer to the reference conductivity profiles than the original Wd-emax model. Finally, we use the obtained 3-D conductivity distributions to calculate the corresponding 3-D water distribution in the upper mantle using thermodynamical and compositional models coupled to the electrical-conductivity laboratory measurement of individual mantle constituents. [ABSTRACT FROM AUTHOR]

SUMMARY: The mechanism responsible for the lateral expansion and uplift of the eastern Tibetan Plateau remains a topic of ongoing debate, partly due to discrepancies in the results of seismic velocity and anisotropy. In local earthquake tomography, hypocentral uncertainties can cause significant errors in the tomographic model. However, this issue has received limited attention in previous studies. In this work, we employ the weighted least-squares (WLS) method to solve the tomographic inversion problem. A power exponent coefficient, which is called weighting level, is introduced into the weighting matrix to control the relative contribution of the data with different hypocentral errors to the final tomographic result. Our data set contains high-quality Pg, Pn and Sg arrival times of local earthquakes recorded by the dense Chinese seismic network in eastern Tibet during 2008–2022. We comprehensively analyse the inversion results derived from the WLS inversions with different weighting levels to evaluate the robustness of isotropic velocity anomalies and azimuthal anisotropy. The most robust feature of our results is a striking low-velocity (low-Vp) zone surrounded by high-velocity (high-Vp) anomalies and fault parallel fast-velocity directions (FVDs) of azimuthal anisotropy in the lower crust beneath the western side of the Longmenshan fault zone. Taking into account many previous results of the region, we deem that the low-Vp zone reflects hot and wet upwelling flow from the deep asthenosphere, which ascends to the lower crust along the fault zone. At the NE margin of the Tibetan Plateau, significant low-Vp anomalies exist in the lower crust and the FVDs are consistent with the motion direction of the Tibetan block revealed by GPS (Global Positioning System) observations. We think that lower crustal flow exists beneath NE Tibet, which controls the plateau expansion toward the northeast. A low-Vp anomaly appears at 30 km depth beneath the Sichuan Basin. However, as the weighting level increases, the amplitude of this low-Vp anomaly decreases by more than 6 per cent, suggesting that this low-Vp anomaly has a lower accuracy than the other features. [ABSTRACT FROM AUTHOR]

Full waveform inversion (FWI) has proved to be a reliable tool for high-resolution imaging of lithospheric structures at various depths down to the upper mantle. However, when the size of the model is large, the computational burden is significative and applications are restricted to low frequencies. To tackle this issue, we developed a new 2D time-domain hybrid method to simulate high-frequency teleseismic body waves propagating through a local heterogeneous elastic Earth model: The frequency-wavenumber (FK) integration method is coupled with the staggered grid finite difference method (SGFD). The FK method is used to compute the wavefield due to obliquely incident plane P and SV waves in a 1D multilayered half space that excites a local heterogeneous region. Inside this region, the velocity-stress staggered grid FD method (SGFD) is used to accurately deal with wave propagation in heterogeneous media. Spurious waves that might be generated at the boundaries of the local region are avoided using convolutional perfectly matched layers (CPML). This new hybrid method inherits the low-memory requirements of the FK method and the accuracy, efficiency and easy implementation of the SGFD. The new hybrid method is benchmarked against the analytical FK method for some canonical models and shows good agreement with analytical solutions. Subsequently, our modeling tool is incorporated into a full waveform inversion algorithm adapted for teleseismic configurations to invert the incident P wave and its coda. The inversion is carried out using a gradient approach that is efficiently implemented via the adjoint-state method. The results suggest that our hybrid method and FWI algorithm represent a valuable tool for 2D forward and inverse regional applications using teleseismic data sets. [ABSTRACT FROM AUTHOR]

SUMMARY: A major challenge in seismic tomography consists in quantifying and representing model resolution and uncertainty, particularly at global scales. This information is crucial for interpretations of tomographic images and their technical application in geodynamics. However, due to large computational costs, there have been only few attempts so far to coherently analyse the spatially varying resolving power for a complete set of model parameters. Here, we present a concept for an effective evaluation and global representation of the 3-D resolution information contained in a full set of averaging kernels. In our case, these kernels are constructed using the 'Subtractive Optimally Localized Averages' (SOLA) method, a variant of classic Backus-Gilbert inversion suitable for global tomography. Our assessment strategy incorporates the following steps: (1) a 3-D Gaussian function is fitted to each averaging kernel to measure resolution lengths in different directions and (2) we define a classification scheme for the quality of the averaging kernels based on their focus with respect to the estimated 3-D Gaussian, allowing us to reliably identify whether the inferred resolution lengths are robust. This strategy is not restricted to SOLA inversions, but can, for example, be applied in all cases where point-spread functions are computed in other tomographic frameworks. Together with model uncertainty estimates that are derived from error propagation in the SOLA method, our concept reveals at which locations, resolution lengths and interpretations of model values are actually meaningful. We finally illustrate how the complete information from our analysis can be used to calibrate the SOLA inversion parameters—locally tunable target resolution kernels and trade-off parameters—without the need for visual inspection of the individual resulting averaging kernels. Instead, our global representations provide a tool for designing tomographic models with specific resolution-uncertainty properties that are useful in geodynamic applications, especially for linking seismic inversions to models of mantle flow. [ABSTRACT FROM AUTHOR]

SUMMARY: Minimum-structure, or Occam's style of, inversion introduces a regularization function into the underdetermined geophysical inverse problems to stabilize the inverse problem and mitigate its non-uniqueness. The regularization function is typically designed such that it can incorporate a priori information into the inversion framework, thus constructing models that have more plausible representations of the true Earth's subsurface structure. One type of a priori information is geological orientation information such as strike, dip and tilt angles of the subsurface structure. This type of information can be incorporated into inverse problems through the roughness operators. Designing such roughness operators for inversion frameworks using unstructured tetrahedral meshes is not as straightforward as for inversion frameworks using structured meshes due to the arbitrary and complex geometry of unstructured meshes. Researchers have developed methods which allow us to incorporate geological orientation information into inversion frameworks with unstructured tetrahedral meshes. The majority of these methods consider each cell in a package with its neighbours, hence, the constructed models are not as sharp as desired if the regularization function is measured using an $\ell _1$ -type measure instead of the $\ell _2$ norm. To address this issue, we propose a method that calculates the directional derivatives of physical property differences between two adjacent cells normalized by the distance between the cell centroids. This approach is able to both incorporate geological orientation information into the inversion framework and construct models with sharp boundaries for the scenarios in which the regularization term is quantified by an $\ell _1$ -type measure. This method is an integral-based approach, therefore, the roughness operators are scaled appropriately by the cell volumes, which is an important characteristic for the inversions with unstructured meshes. To assess the performance and the capability of the proposed method, it was applied to 3-D synthetic gravity and magnetotelluric examples. The gravity example was also used to investigate the impact of applying the depth weighting function inside and outside the roughness operators for the scenarios that the model objective function is measured by an $\ell _1$ norm. The examples show that the proposed method is able to construct models with a reasonable representation of the strike and dip directions of the true subsurface model with sharper boundaries if the regularization function is quantified by an $\ell _1$ -type measure. The examples also demonstrate the proposed method behaves numerically well, and has a fast convergence rate. [ABSTRACT FROM AUTHOR]

SUMMARY: Reconciling rock unit boundary geometry is crucial for geological and geophysical studies aiming to achieve a comprehensive 3-D subsurface model. To create a unified 3-D parametrization suitable for both geological modelling and geophysical inversion, an integrated approach utilizing implicit modelling is essential. However, a key challenge lies in encapsulating all pertinent information within the 3-D model, ensuring compatibility with the utilized data sets and existing constraints. In this study, we present a workflow that enables the generation of an integrated 3-D subsurface model primarily using gravity and reflection seismic data sets. Our approach involves a cooperative geophysical inversion workflow, which incorporates the inverted model from the reflection seismic data while leveraging sparse petrophysical information. Despite advances in integrated modelling, the incorporation of implicit modelling approaches in cooperative inversion workflows remains unexplored. In our gravity inversion process, we use a generalized level set method to refine the boundaries of rock units in the prior model. We integrate the inverted model, derived from seismic and other sparse petrophysical data sets, to create a comprehensive 3-D prior model. To enhance the integration of reflection seismic data sets in the level set inversion, we introduce a weighting uncertainty matrix containing constraint terms. This step refines the model's accuracy and ensures greater consistency. Finally, we search for any missing rock units within inverted model through nucleation investigations. The introduced methodology has undergone successful testing in the Boulia region (Southern Mount Isa, Queensland), utilizing two 2-D reflection seismic profiles and regional gravity data sets. This study primarily aims to reconstruct the geometry of major structures within the basement units and the basin at a regional scale. By combining seismic profiles and gravity data sets with constraining information, we are able to create a 3-D model of the area that accurately represents distinct rock units and their boundary geometries. Additionally, relevant legacy data sets and prior modelling results from the region have been incorporated and refined, ensuring that the final model aligns with all available knowledge about the area. [ABSTRACT FROM AUTHOR]

Gravity gradient data can show the structural features of geological bodies in the shallow lithosphere with higher sensitivity and resolution than conventional gravity data. Gravity gradient inversion can be applied to obtain the lithospheric density structures of geological bodies. However, as with gravity data, gravity gradient data have no inherent depth resolution. The methods of gravity gradient depth imaging and gravity gradient inversion are integrated in this study. The depth imaging method is effective for calculations without prior information and iterative computations. As the parameters in the depth weighting function should be chosen from a set of values used in inversion tests of synthetic data, which brings some uncertainties, the depth imaging results of gravity gradient are introduced into the depth weighting function. Several synthetic models are tested to demonstrate the advantages and features of the effective integrated method. Finally, the integrated method is applied to the interpretation of the GOCE satellite gravity gradient tensors over the northeastern margin of the Qinghai–Tibet Plateau. The results reveal that in the crust of the study area, the distribution of density anomalies is more in line with the mechanism of the crustal flow model, in the upper mantle of the study area, the density anomalies are mainly influenced by the high heat flow environment. [ABSTRACT FROM AUTHOR]

In this paper, we use a new workflow to substantiate the characterization of a prominent, deep sediment conductor in the hyperextended Bjørnøya Basin (SW Barents Sea) previously identified in smooth resistivity models from 3-D deterministic inversion of magnetotelluric data. In low-dimensionality environments like layered sedimentary basin, 1-D Bayesian inversion can be advantageous for a thorough exploration of the solution space, but the violation of the 1-D assumption has to be efficiently handled. The primary geological objectives of this work is therefore preceded by a secondary task: the application of a new machine learning approach for handling the 1-D violation assumption for 21 MT field stations in the Barents Sea. We find that a decision tree can adequately learn the relationship between MT dimensionality parameters and the 1-D–3-D residual response for a training set of synthetic models, mimicking typical resistivity structures of the SW Barents Sea. The machine learning model is then used to predict the dimensionality compensation error for MT signal periods ranging of 1–3000 s for 21 receivers located over the Bjørnøya Basin and Veslemøy High. After running 1-D Bayesian inversion, we generated a posterior resistivity distribution for an ensemble of 6000 1-D models fitting the compensated MT data for each 21 field stations. The proportion of 1-D models showing ρ < 1 Ω·m is consistently beyond 80 per cent and systemically reaches a maximum of 100 per cent in the Early Aptian–Albian interval in the Bjørnøya Basin. In hyperextended basins of the SW Barents Sea, the dimensionality compensation workflow has permitted to refine the characterization of the deep basin conductor by leveraging the increased vertical resolution and optimal used of MT data. In comparison, the smooth 3-D deterministic models only poorly constrained depth and lateral extent of the basin anomaly. The highest probability of finding ρ < 1 Ω·m is robustly assigned to the syn-tectonic Early Aptian–Albian marine shales, now buried at 6–8 km depth. Based on a theoretical two phase fluid-rock model, we show that the pore fluid of these marine shales must have a higher salinity than seawater to explain the anomaly ρ < 1 Ω·m. Therefore, the primary pore fluid underwent mixing with a secondary brine during rifting. Using analogue rift systems in palaeomargins, we argue that two possible secondary brine reservoir may contribute to deep saline fluid circulation in the hyperextended basin: (1) Permian salt-derived fluid and, (2) mantle-reacted fluid from serpentinization. [ABSTRACT FROM AUTHOR]

Seismic gradient measurements from distributed acoustic sensors and rotational sensors are becoming increasingly available for field surveys. These measurements provide a wealth of information and are currently being considered for many applications such as earthquake detection and subsurface characterizations. In this work, using a simple 2-D numerical approach, we tackle the implications of such wavefield gradient measurements on full waveform inversion (FWI) techniques using a simple 2-D numerical test. In particular, we study the impact of the wavefield gradient measurement sensitivity to heterogeneities that are much smaller than the minimum wavelength. Indeed, as shown through the homogenization theory, small-scale heterogeneities induce an unexpected coupling of the strain components to the wavefield gradient measurement. We further show that this coupling introduces a potential limitation to the FWI results if it is not taken into account. We demonstrate that a gradient measurement-based FWI can only reach the accuracy of a classical displacement field-based FWI if the coupling coefficients are also inverted. Furthermore, there appears to be no specific gain in using gradient measurements instead of conventional displacement (or velocity, acceleration) measurements to image structures. Nevertheless, the inverted correctors contain fine-scale heterogeneities information that could be exploited to reach an unprecedented resolution, particularly if an array of receivers is used. [ABSTRACT FROM AUTHOR]

Investigating subsurface shear wave velocity (v s) structures using surface wave dispersion data involves minimizing a misfit function that is commonly solved through gradient-based optimization. Sensitivity kernels for model updates are commonly estimated using numerical differentiation, variational methods or implicit functions which however, may involve numerical instability and computational challenges when dealing with complex velocity models and large data sets. In this study, we propose a novel surface wave inversion framework in which error-free gradients are calculated by automatic differentiation (AD) and forward modelling is implemented by convenient computational graphs in the state-of-the-art deep learning framework. The AD-based inversion approach is first validated using two synthetic data sets. Then, the subsurface structures at three distinct locations, namely the Great Plains and the Long Beach in the US and Tong Zhou in China, are also derived using this method with seismic ambient noise data, which show nice consistency with those obtained using traditional methods. With the significantly improved computational efficiency, a great number of initial models can be inverted simultaneously to mitigate the impact of local minima and to estimate the uncertainty in the invert models. We have developed a new surface wave inversion package named ADsurf based on automatic differentiation and computational graphs in the deep learning framework, and its computational efficiency is also compared with the traditional finite-difference-based gradient estimation approach. While a great number of intriguing studies on the geophysical inverse problems have been conducted recently using deep learning for end-to-end mapping, the use of AD provided in the in the deep learning frameworks to assist and expedite the gradient computations are still underexploited in geophysics. Thus, it is expected that various geophysical inverse problems in many different areas beyond the surface wave inversion can also be tackled with this new paradigm in the future. [ABSTRACT FROM AUTHOR]

Seismic full-waveform inversion (FWI) or waveform inversion (WI) has gained extensive attention as a cutting-edge imaging method, which is expected to reveal the high-resolution images of complex geological structures. In this paper, we regard each 1-D signal in the inversion system as a 1-D probability distribution, then use the Jensen–Shannon divergence from information theory to measure the discrepancy between the predicted and observed signals, and finally implement a novel 2-D multiparameter shallow-seismic WI (MSWI). Essentially, the novel approach achieves an implicit weighting along the time-axis for each 1-D adjoint source defined by the classical WI (CWI), thus enhancing the extra illumination for a deeper medium compared with the CWI. By evaluating the inversion results of the two-layer model and fault model, the reconstruction accuracy for S -wave velocity and density of the new method is increased by about 30 and 20 per cent compared with that of the CWI under the same conditions, respectively. The reconstruction performance for P -wave velocity of these two methods is almost equal. In addition, the new 2-D MSWI is also resilient to white Gaussian noise in the data. Numerically, the inversion system has almost the strongest sensitivities to the S -wave velocity and density, performing the poorest sensitivity to the P -wave velocity. Finally, we test the novel method with a detection case for a power tunnel. [ABSTRACT FROM AUTHOR]

Full-waveform inversion (FWI) is an effective method for imaging subsurface properties using sparsely recorded data. It involves solving a wave propagation problem to estimate model parameters that accurately reproduce the data. Recent trends in FWI have seen a renewed interest in extended methodologies, among which source extension methods leveraging reconstructed wavefields to solve penalty or augmented Lagrangian (AL) formulations have emerged as robust algorithms, even for inaccurate initial models. Despite their demonstrated robustness on synthetic data, challenges remain, such as the lack of a clear physical interpretation and reliance on difficult-to-compute least-squares (LS) wavefields. Moreover, the literature lacks a general and through comparison of these source extension methods with each other and with the standard FWI. This paper is divided into three critical parts. In the first, a novel formulation of these methods is explored within a unified Lagrangian framework. This novel perspective permits the introduction of alternative algorithms that use LS multipliers instead of wavefields. These multiplier-oriented variants appear as regularizations of the standard FWI, are suitable to the time domain, offer tangible physical interpretations, and foster enhanced convergence efficiency. The second part of the paper delves into understanding the underlying mechanisms of these techniques. This is achieved by solving the associated non-linear equations using iterative linearization and inverse scattering methods. The paper provides insight into the role and significance of Lagrange multipliers in enhancing the linearization of the equations. It explains how different methods estimate multipliers or make approximations to increase computing efficiency. Additionally, it presents a new physical understanding of the Lagrange multiplier used in the AL method, highlighting how important it is for improving algorithm performance when compared to penalty methods. In the final section, the paper presents numerical examples that compare different methods within a unified iterative algorithm, utilizing benchmark Marmousi and 2004 BP salt models. [ABSTRACT FROM AUTHOR]

Shear wave velocity (Vs) is a fundamental property of elastic media whose estimation from PS converted waves is challenging and requires modelling the boundary where P to S conversion occurs. This paper presents a PS tomography where seismic wave conversion/reflection points correspond to reflectors modelled with the level-set function set to zero [ ϕ (x, z) = 0]. The proposed method aims for stable Vs inversion in a seismic acquisition setting using multicomponent receivers. Synthetic models simulating true Vs , Vp and the location of the geological reflector are used in the study. The inversion starts by locating a flat reflector, ϕ (x, z) = 0, which defines the zone Ω1 between the surface and the reflector, where the initial Vs and Vp fields are also set. To calculate the traveltimes of incident PT (P wave that propagates in Ω1 from source to the reflector), converted PS and reflected PP waves, for both observed and modelled data (forward problem), the methodology proposed by Rawlinson and Sambridge is adopted. This method uses the arrival times of the P waves, Tpt , from the seismic source at each reflector point as secondary sources generating the times Tps and Tpp. These times are calculated as a solution to the eikonal equation by using the Fast Marching method. The PS and PP residual times are minimized by updating Vs , Vp and ϕ (x, z) = 0 through adjoint variables designed from a formulation using Lagrange Multipliers in a variational context. The performance of the algorithm is evaluated for models with synclinal, sinusoidal and monoclinal reflector geometries using numerical tests considering the inversion of: (1) ϕ , given the true values of Vs and Vp ; (2) ϕ and Vs , given the true value of Vp ; (3) ϕ and Vp , given the true value of Vs and (4) the three parameters ϕ, Vs and Vp , simultaneously. Good results are obtained by inverting Vs and ϕ , given the true value of Vp. The simultaneous inversion of the three parameters exhibits promising results, despite the illumination problems caused by the different distribution of the PS, PP and PT time gradients due to the geometry of the reflectors and the acquisition setting (sources–receivers in the same plane). The proposed tomography estimates Vs and reflector positions which could help in statics corrections and improve the lithological characterization of near surface. [ABSTRACT FROM AUTHOR]

We apply the adjoint method to efficiently calculate sensitivity kernels for long-period seismic spectra with respect to structural and source parameters. Our approach is built around the solution of the frequency-domain equations of motion using the direct solution method (DSM). The DSM is currently applied within large-scale mode coupling calculations and is also likely to be useful within finite-element type methods for modelling seismic spectra that are being actively developed. Using mode coupling theory as a framework for solving both the forward and adjoint equations, we present numerical examples that focus on the spectrum close to four eigenfrequencies (the low-frequency mode, 0 S 2, and higher frequency modes, namely 2 S 2, 0 S 7 and 0 S 10 for comparison). For each chosen observable, we plot sensitivity kernels with respect to 3-D perturbations in density and seismic wave speeds. We also use the adjoint method to calculate derivatives of observables with respect to the matrices occurring within mode coupling calculations. This latter approach points towards a generalization of the two-stage splitting function method for structural inversions that does not rely on inaccurate self-coupling or group-coupling approximations. Finally, we verify through direct calculation that our sensitivity kernels correctly predict the linear dependence of the chosen observables on model perturbations. In doing this, we highlight the importance of non-linearity within inversions of long-period spectra. [ABSTRACT FROM AUTHOR]

The increasing need for petrophysical parameter inversion has spurred the advancement of full-waveform inversion (FWI). This study introduces a new inversion method that can directly estimate petrophysical parameters from seismic waveform data. Specifically, we modified the regular elastic wave equation by including Biot poroelastic parameters and designed an inversion workflow based on the framework of elastic FWI. The newly introduced inversion algorithms have demonstrated their effectiveness and accuracy using the Marmousi model and a synthetic model related to the gas hydrate deposits based on the Shenhu area of the South China Sea. Our work demonstrates the feasibility of directly inverting petrophysical parameters and assessing if a deposit contains gas hydrate. [ABSTRACT FROM AUTHOR]

Time-lapse seismic full-waveform inversion (FWI) provides estimates of dynamic changes in the Earth's subsurface by performing multiple seismic surveys at different times. Since FWI problems are highly non-linear and non-unique, it is important to quantify uncertainties in such estimates to allow robust decision making based on the results. Markov chain Monte Carlo (McMC) methods have been used for this purpose, but due to their high computational cost, those studies often require a pre-existing accurate baseline model and estimates of the locations of potential velocity changes, and neglect uncertainty in the baseline velocity model. Such detailed and accurate prior information is not always available in practice. In this study we use an efficient optimization method called stochastic Stein variational gradient descent (sSVGD) to solve time-lapse FWI problems without assuming such prior knowledge, and to estimate uncertainty both in the baseline velocity model and the velocity change over time. We test two Bayesian strategies: separate Bayesian inversions for each seismic survey, and a single joint inversion for baseline and repeat surveys, and compare the methods with standard linearized double difference inversion. The results demonstrate that all three methods can produce accurate velocity change estimates in the case of having fixed (exactly repeatable) acquisition geometries. However, the two Bayesian methods generate significantly more accurate results when acquisition geometries changes between surveys. Furthermore, joint inversion provides the most accurate velocity change and uncertainty estimates in all cases tested. We therefore conclude that Bayesian time-lapse inversion using a joint inversion strategy may be useful to image and monitor subsurface changes, in particular where variations in the results would lead to different consequent decisions. [ABSTRACT FROM AUTHOR]

Density is an important parameter for both geological research and geophysical exploration. However, for model-driven seismic inversion methods, high-fidelity density inversion is challenging due to seismic wave traveltime insensitivity to density and crosstalk that density has with velocity. To circumvent the challenge of density inversion, some inversion methods treat density as a constant value or derive density from velocity through empirical equation. On the other hand, deep learning approaches are completely driven by data and have strong target-oriented characteristics, offering a new way to solve multiparameter coupling problems. Nevertheless, the accuracy of the inversion results of data-driven algorithms is directly related to the amount and diversity of the training data, and thus, they lack the universality of model-driven algorithms. To achieve accurate density inversion, we propose a simultaneous inversion algorithm for velocity and density that combines the advantages of data- and model- driven approaches: A neural network model (U-T), combining the U-net and Transformer architectures, is proposed to construct non-linear mappings between seismic data as inputs and the velocity and density as predictions. Next, the model-driven inversion algorithm uses the U-T prediction as the initial model to obtain the final accurate solution. In the model-driven module, envelope-based sparse constrained deconvolution is used to obtain full-band seismic data, while a variable dominant frequency full waveform inversion algorithm is used to perform multiscale inversion, ultimately leading to accurate inversion results of velocity and density. The performance of the algorithm on the Sigsbee2A and Marmousi models demonstrates its effectiveness. [ABSTRACT FROM AUTHOR]

Full-waveform inversion (FWI) of seismic data provides quantitative constraints on subsurface structures. Despite its widespread success, FWI of data around the critical angle is challenging because of the abrupt change in amplitude and phase at the critical angle and the complex waveforms, especially in the presence of a sharp velocity contrast, such as at the Moho transition zone (MTZ). Furthermore, the interference of refracted lower crustal (Pg) and upper mantle (Pn) arrivals with the critically reflected Moho (PmP) arrivals in crustal and mantle studies makes the application of conventional FWI based on linearized model updates difficult. To address such a complex relationship between the model and data, one should use an inversion method based on a Bayesian formulation. Here, we propose to use a Hamiltonian Monte Carlo (HMC) method for FWI of wide-angle seismic data. HMC is a non-linear inversion technique where model updates follow the Hamiltonian mechanics while using the gradient information present in the probability distribution, making it similar to iterative gradient techniques like FWI. It also involves procedures for generating distant models for sampling the posterior distribution, making it a Bayesian method. We test the performance and applicability of HMC based elastic FWI by inverting the non-linear part of the synthetic seismic data from a three-layer and a complex velocity model, followed by the inversion of wide-angle seismic data recorded by two ocean bottom seismometers over a 70 Ma old oceanic crustal segment in the equatorial Atlantic Ocean. The inversion results from both synthetic and real data suggest that HMC based FWI is an appropriate method for inverting the non-linear part of seismic data for crustal studies. [ABSTRACT FROM AUTHOR]

Material density remains poorly constrained in seismic imaging problems, yet knowledge of density would provide important insight into physical material properties for the interpretation of subsurface structures. We test the sensitivity to subsurface density contrasts of spatial and temporal gradients of seismic ambient noise wavefields, using wave equation inversion (WEI), a form of seismic gradiometry. Synthetic results for 3-D acoustic media suggest that it is possible to estimate relative density structure with WEI by using a full acoustic formulation for wave propagation and gradiometry. We show that imposing a constant density assumption on the medium can be detrimental to subsurface seismic velocity images. By contrast, the full acoustic formulation allows us to estimate density as an additional material parameter, as well as to improve phase velocity estimates. In 3-D elastic media, severe approximations in the governing wave physics are necessary in order to invert for density using only an array of receivers on the Earth's free surface. It is then not straightforward to isolate the comparatively weak density signal from the influence of phase velocity using gradiometric WEI. However, by using receivers both at the surface and in the shallow subsurface we show that it is possible to estimate density using fully elastic volumetric WEI. [ABSTRACT FROM AUTHOR]

We propose a new method to obtain 3-D shear wave velocity structures for regions with undulating topography based on surface wave dispersion data. In our method, we assume that surface waves propagate along a curved free surface, and the dispersion effects of these waves can be modelled as a frequency-dependent ray tracing problem. We use the shortest-path method to address off-great-circle propagation arising from inhomogeneous elastic parameters and topographic variations. We then apply our method to both synthetic and real-data inversions and demonstrate that ignoring topographic effects may significantly distort the inverted images. Finally, we analysed the accuracy of our method and provided a rule-of-thumb principle to quantitatively assess the need to account for topographic effects in surface wave tomography for a region. [ABSTRACT FROM AUTHOR]

More accurate inversion of source fault geometry and slip parameters under the constraint of the Bayesian algorithm has become a research hotspot in the field of geodetic inversion in recent years. In nonlinear inversion, the determination of the weight ratio of the joint inversion of multisource data is more complicated. In this context, this paper proposes a simple and easily generalized weighting method for inversion of source fault parameters by joint geodetic multisource data under the Bayesian framework. This method determines the relative weight ratio of multisource data by root mean square error (RMSE) value and can be extended to other nonlinear search algorithms. To verify the validity of the method in this paper, this paper first sets up four sets of simulated seismic experiment schemes. The inversion results show that the joint inversion weighting method proposed in this paper has a significant decrease in the large residual value compared with the equal weight joint inversion and the single data source joint inversion method. The east–west deformation RMSE is 0.1458 mm, the north–south deformation RMSE is 0.2119 mm and the vertical deformation RMSE is 0.2756 mm. The RMSEs of the three directions are lower than those of other schemes, indicating that the proposed method is suitable for the joint inversion of source parameters under Bayesian algorithm. To further verify the applicability of the proposed method in complex earthquakes, the source parameters of the Maduo earthquake were inverted using the method of this paper. The focal depth of the inversion results in this paper is closer to the focal depth released by the GCMT agency. In terms of strike angle and dip angle, the joint inversion in this paper is also more inclined to the GCMT results. The joint inversion results generally conform to the characteristics of left-lateral strike-slip, which shows the adaptability of this method in complex earthquakes. [ABSTRACT FROM AUTHOR]

In near-surface surveys, shallow-seismic and ground-penetrating radar (GPR) full-waveform inversions (FWIs) have received increasing attention because of their ability to reconstruct high-resolution subsurface models. However, they have different sensitivities to the same targets and thus may yield conflicting geophysical parameter models. To solve this issue, we have developed an indirect joint petrophysical inversion (JPI) integrating shallow-seismic and multi-offset GPR data. These data are used to reconstruct porosity and saturation whereby we use only strong sensitivities between petrophysical and geophysical parameters. To promote its field application, we proposed an input strategy to avoid measuring rock matrix parameters and make indirect JPI more robust. We apply indirect JPI to the field data acquired in Rheinstetten, Germany and find that it reveals the mechanical, electrical and petrophysical properties more reliably than individual inversions. The reconstructed models are assessed by direct-push technology, borehole sample measurements and migrated GPR image. Indirect JPI can fit seismic and GPR observed data simultaneously and provide consistent multiparameter models, which are hard to achieve by FWIs and individual petrophysical inversions. We also find that the method is robust when there are uncertainties in petrophysical a priori information. Overall, the field example proves the great potential of using indirect JPI to solve real-world problems. [ABSTRACT FROM AUTHOR]

Previous studies have shown that it is difficult to determine whether the 2015 Pishan earthquake occurred on a uniform fault or a ramp-flat fault with variable dip angles due to the similar goodness of data fit to coseismic and afterslip models on these two fault models. Here, we first present the InSAR deformation obtained from both ascending and descending orbits, covering the coseismic period and cumulative 5-yr period after the 2015 Pishan earthquake. We then determine the preferred fault geometry by the spatial distributions between the positive Coulomb failure stress change triggered by main shock and the afterslip. Based on the preferred fault model, we finally use a combined model to determine the contributions of elastic and viscoelastic deformation in the post-seismic deformation. We find that the Pishan earthquake prefers to occur on a ramp-flat fault, and the coseismic slip is mainly distributed at a depth of 9–13 km, with a maximum slip of about 1.3 m. The post-seismic deformation is primarily governed by afterslip, as the poroelastic rebound-induced deformation fails to account for the observed post-seismic deformation and the contributions from the viscoelastic relaxation mechanism can be considered negligible in the combined model. Moreover, the modelled stress-driven afterslip and observed kinematic afterslip have good consistency, and the difference between the root mean square error of the two afterslip models is only 4.3 mm. The results from the afterslip model indicate that both of the updip and downdip directions distribute the afterslip, and slip in the updip direction is greater than that of the downdip direction. Meanwhile, the maximum cumulative afterslip after 5 yr is approximately 0.26 m which is equivalent to a released seismic moment of a M w 6.47. [ABSTRACT FROM AUTHOR]

In the context of Bayesian inversion, we consider sequential Monte Carlo (SMC) methods that provide an approximation of the posterior probability density function and the evidence (marginal likelihood). These particle approaches build a sequence of importance sampling steps between gradually tempered distributions evolving from the prior to the posterior PDF. To automate the definition of the tempering schedule, adaptive SMC (ASMC) allows tuning the temperature increments on-the-go. One general challenge in Bayesian inversions is the computational burden associated with expensive, high-fidelity forward solvers. Lower-fidelity surrogate models are interesting in this context as they can emulate the response of expensive forward solvers at a fraction of their cost. We consider surrogate modelling within ASMC and introduce first an approach involving surrogate modelling only, in which either prior samples are used to train the surrogate, or the surrogate model is retrained by updating the training set during the inversion. In our implementation, we rely on polynomial chaos expansions for surrogate modelling, principal component analysis for model parametrization and a ground-penetrating radar cross-hole tomography problem with either an eikonal or finite-difference time-domain solver as high-fidelity solver. We find that the method based on retraining the surrogate during the inversion outperforms the results obtained when only considering prior samples. We then introduce a computationally more expensive multifidelity approach including a transition to the high-fidelity forward solver at the end of the surrogate-based ASMC run leading to even more accurate results. Both methods result in speed-ups that are larger than one order of magnitude compared to standard high-fidelity ASMC inversion. [ABSTRACT FROM AUTHOR]

Muography 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 analysed the altitudinal impact on the data measurements and proposed a simplified method to eliminate it. We also developed the seed algorithm, a novel 3-D reconstruction algorithm, that can reconstruct anomalies located far away from the detectors and avoid their elongation along the observed directions. Benefitting from the seed algorithm, a low-density-contrast ore body 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. [ABSTRACT FROM AUTHOR]

This paper deals with in situ characterization of the small-strain shear-wave velocity VS and damping ratio DS from an advanced interpretation of Multi-channel Analysis of Surface Waves (MASW) surveys. A new approach based on extracting Rayleigh wave data using the CFDBFa method has been discussed in the companion paper. This paper focuses on mapping the experimental Rayleigh wave phase velocity and attenuation into profiles of VS and DS versus depth, which is achieved through a joint inversion procedure. The joint inversion of phase velocity and attenuation data utilizes a newly developed Monte Carlo global search algorithm, which implements a smart sampling procedure. This scheme exploits the scaling properties of the solution of the Rayleigh eigenvalue problem to modify the trial earth models and improve the matching with the experimental data. Thus, a reliable result can be achieved with a limited number of trial ground models. The proposed algorithm is applied to the inversion of synthetic data and of experimental data collected at the Garner Valley Downhole Array site, as described in the companion paper. In general, inverted soil models exhibit well-defined VS profiles, whereas DS profiles are affected by larger uncertainties. Greater uncertainty in the inverted DS profiles is a direct result of higher variability in the experimental attenuation data, the limited wavelength range at which reliable values of attenuation parameters can be retrieved, and the sensitivity of attenuation data to both DS and VS. Nonetheless, the resulting inverted earth models agree well with alternative in situ estimates and geological data. The results stress the feasibility of retrieving both stiffness and attenuation parameters from active-source MASW testing and the effectiveness of extracting in situ damping ratio estimates from surface wave data. [ABSTRACT FROM AUTHOR]

Plate-coupling estimates and previous seismicity indicate that portions of the Makran megathrust of southern Pakistan and Iran are partially coupled and have the potential to produce future magnitude 7+ earthquakes. However, the GPS observations needed to constrain coupling models are sparse and lead to an incomplete understanding of regional earthquake and tsunami hazard. In this study, we assess GPS velocities for plate coupling of the Makran subduction zone with specific attention to model resolution and the accretionary prism rheology. We use finite element model-derived Green's functions to invert for the interseismic slip deficit under both elastic and viscoelastic Earth assumptions. We use the model resolution matrix to characterize plate-coupling scenarios that are consistent with the limited spatial resolution afforded by GPS observations. We then forward model the corresponding tsunami responses at major coastal cities within the western Indian Ocean basin. Our plate-coupling results show potential segmentation of the megathrust with varying coupling from west to east, but do not rule out a scenario where the entire length of the megathrust could rupture in a single earthquake. The full subduction zone rupture scenarios suggest that the Makran may be able to produce earthquakes up to M w 9.2. The corresponding tsunami model from the largest earthquake event (M w 9.2) estimates maximum wave heights reaching 2–5 m at major port cities in the northern Arabian Sea region. Cities on the west coast of India are less affected (1–2 m). Coastlines bounding eastern Africa, and the Strait of Hormuz, are the least affected (<1 m). [ABSTRACT FROM AUTHOR]

This study acquires the coseismic deformation field and the high-frequency dynamic displacement of the M W 7.4 earthquake that occurred in Maduo, China, on 2021 May 22, based on the BeiDou Navigation Satellite System (BDS), and the comparison with the results obtained by the Global Positioning System (GPS) reveals that the two systems are certain differences in their ability to acquire the coseismic deformation field. The maximum difference in the horizontal coseismic deformation is <5 mm, and the maximum difference in the vertical coseismic deformation is 8.7 mm. The dynamic displacement waveforms of the 2021 M W 7.4 Maduo earthquake acquired by BDS and GPS are very similar, which confirms that BDS can acquire ground-shaking images with an accuracy comparable to that of GPS. Based on the empirical relationship equation of the peak ground displacement (PGD) and moment magnitude (M W), this study verifies and calculates both the M W of the 2021 M W 7.4 Maduo earthquake and the error and finds that the M W can be quickly and accurately obtained by using the empirical PGD and M W equations, and this M W value can be used as a supplementary means of calibrating the M W of the large earthquake early warning systems, which can be quickly determined by seismic wave data. Finally, by comparing the slip distributions inverted from the BDS and GPS coseismic deformation fields, this study finds that BDS is equally effective as GPS. [ABSTRACT FROM AUTHOR]

Unexploded ordnance (UXO), often deployed in WWI and WWII, are still present on the seabed posing a risk to seabed construction. UXO surveys identify potential UXO targets which are then investigated and removed where necessary; however, approximately 96 per cent of targets investigated on projects included in this study were non-UXO (false positives). Here we investigate methods of improving statistical confidence in targets picked through the improved classification of magnetic response, utilizing simulation methods and historic data. A simulation of the induced magnetism of UXO was developed to model the most likely dimensions and orientation of a given source item. Modelling results and measured variables were combined in an optimization algorithm, hereon referred to as 'the classifier'. We document the results of integrating the magnetic responses of 3027 targets from 18 offshore wind projects into the classifier. Five of these projects were performed as tests in which the target results were unknown during classification. Finally, the classifier was deployed on a live project, Hollandse Kust Zuid 3 and 4 Offshore Wind Farm, in which it enabled the removal of 67 targets from the from the 1200-target investigation list. The classifier proved to be a valuable additional tool to assist the target rationalization process, reliably and repeatably distinguishing UXO from other items, thereby reducing false positives and repeatably demonstrating no threshold UXO false negatives. Assumptions and limitations of this methodology were highlighted, which must be understood and considered when integrating this approach on future projects. [ABSTRACT FROM AUTHOR]

A large non-double-couple component of a tectonic earthquake indicates that its rupture likely was complex and likely involved multiple faults. Detailed source models of such earthquakes can add to our understanding of earthquake source complexity. The 2007 Martinique earthquake in the Caribbean Sea is one of the largest recent earthquakes with a known large non-double-couple component. It was an intermediate depth intraslab earthquake within the South American plate where it is subducting beneath the Caribbean plate. We applied potency density tensor inversion (PDTI) to teleseismic P waves generated by the 2007 Martinique earthquake to model its source processes and focal mechanism distribution. We identified two focal mechanisms: a strike-slip mechanism with a north–south tension axis (T- axis), and a downdip extension (DDE) mechanism with an east–west T -axis. Rupture by the DDE mechanism was predominant in the northern part of the source region and strike-slip rupture in the southern part. These two focal mechanisms had approximately parallel pressure axes (P- axes) and approximately orthogonal T -axes. The seismic moments released by both types of rupture were almost equal. These results indicate that the 2007 Martinique earthquake had a large non-double-couple component. We identified five subevents with two predominant directions of rupture propagation: two strike-slip subevents propagated to the southeast and three DDE subevents propagated to the east. Although the directions of propagation were consistent for each focal mechanism, each subevent appears to have occurred in isolation. For example, the rupture of one DDE subevent propagated from the edge of the source region back towards the hypocentre. Complex ruptures that include multiple subevents may be influenced by high pore fluid pressure associated with slab dehydration. Our results show that PDTI can produce stable estimates of complex seismic source processes and provide useful information about the sources of complex intermediate depth intraslab earthquakes for which fault geometry assumptions are difficult. [ABSTRACT FROM AUTHOR]

Satellite observations of the geomagnetic field contain signals generated in Earth's interior by electrical currents in the core and by magnetized rocks in the lithosphere. At short wavelengths the lithospheric signal dominates, obscuring the signal from the core. Here we present details of a method to co-estimate separate models for the core and lithospheric fields, which are allowed to overlap in spherical harmonic degree, that makes use of prior information to aid the separation. Using a maximum entropy method we estimate probabilistic models for the time-dependent core field and the static lithospheric field that satisfy constraints provided by satellite observations while being consistent with prior knowledge of the spatial covariance and expected magnitude of each field at its source surface. For the core field, we find that between spherical harmonic degree 13 and 22 power adds coherently to the established structures, and present a synthetic test that illustrates the aspects of the small scale core field that can reliably be retrieved. For the large scale lithospheric field we also find encouraging results, with the strongest signatures below spherical harmonic degree 13 occurring at locations of known prominent lithospheric field anomalies in the northern part of Eastern Europe, Australia and eastern North America. Although the amplitudes of the small scale core field and large scale lithospheric field are likely underestimated we find no evidence that obvious artefacts are introduced. Compared with conventional maps of the core–mantle boundary field our results suggest more localized normal flux concentrations close to the tangent cylinder, and that low latitude flux concentrations occur in pairs with opposite polarities. Future improvements in the recovery of the small scale core field and large scale lithospheric field will depend on whether more detailed prior information can be reliably extracted from core dynamo and lithospheric magnetization simulations. [ABSTRACT FROM AUTHOR]

In geophysical surveys or experiments, recorded data are used to constrain properties of the planetary subsurface, oceans, atmosphere or cryosphere. How the experimental data are collected significantly influences which parameters can be resolved and how much confidence can be placed in the results. Bayesian experimental design methods characterize, quantify and maximize expected information post-experiment—an optimization problem. Typical design parameters that can be optimized are source and/or sensor types and locations, and the choice of modelling or data processing methods to be applied to the data. These may all be optimized subject to various physical and cost constraints. This paper introduces variational design methods, and discusses their benefits and limitations in the context of geophysical applications. Variational methods have recently come to prominence due to their importance in machine-learning applications. They can be used to design experiments that best resolve either all model parameters, or the answer to specific questions about the system to be interrogated. The methods are tested in three schematic geophysical applications: (i) estimating a source location given arrival times of radiating energy at sensor locations, (ii) estimating the contrast in seismic velocity across a stratal interface given measurements of the amplitudes of seismic wavefield reflections from that interface, and (iii) designing a survey to best constrain CO2 saturation in a subsurface storage scenario. Variational methods allow the value of expected information from an experiment to be calculated and optimized simultaneously, which results in substantial savings in computational cost. In the context of designing a survey to best constrain CO2 saturation in a subsurface storage scenario, we show that optimal designs may change substantially depending on the particular questions of interest. We also show that one method, so-called DN design, can be effective at substantially lower computational cost than other methods. Overall, this work demonstrates that optimal design methods could be used more widely in Geophysics, as they are in other scientifically advanced fields. [ABSTRACT FROM AUTHOR]

Bayesian-based data assimilation methods integrate observational data into geophysical forward models to obtain the temporal evolution of an improved state vector, including its uncertainties. We explore the potential of a variant, a particle method, to estimate mechanical parameters of the overriding plate during the interseismic period. Here we assimilate vertical surface displacements into an elementary flexural model to estimate the elastic thickness of the overriding plate, and the locations and magnitudes of line loads acting on the overriding plate to produce flexure. Assimilation of synthetic observations sampled from a different forward model than is used in the particle method, reveal that synthetic seafloor data within 150 km from the trench are required to properly constrain parameters for long wavelength solutions of the upper plate (i.e. wavelength ∼500 km). Assimilation of synthetic observations sampled from the same flexural model used in the particle method shows remarkable convergence towards the true parameters with synthetic on-land data only for short to intermediate wavelength solutions (i.e. wavelengths between ∼100 and 300 km). In real-data assimilation experiments we assign representation errors due to discrepancies between our incorrect or incomplete physical model and the data. When assimilating continental data prior to the 2011 M w Tohoku-Oki earthquake (1997–2000), an unrealistically low effective elastic plate thickness for Tohoku of ∼5–7 km is estimated. Our synthetic experiments suggest that improvements to the physical forward model, such as the inclusion of a slab, a megathrust interface and viscoelasticity of the mantle, including accurate seafloor data, and additional geodetic observations, may refine our estimates of the effective elastic plate thickness. Overall, we demonstrate the potential of using the particle method to constrain geodynamic parameters by providing constraints on parameters and corresponding uncertainty values. Using the particle method, we provide insights into the data network sensitivity and identify parameter trade-offs. [ABSTRACT FROM AUTHOR]

Gravity surveys constitute an important method for investigating the Earth's interior based on density contrasts related to Earth material differentials. Because lithology depends on the environment and the period of formation, there are generally clear boundaries between rocks with different lithologies. Inversions with convex functions for approximating the L 0 norm are used to detect boundaries in reconstructed models. Optimizations can easily be found because of the convex transformations; however, the volume of the reconstructed model depends on the weighting parameter and the density constraint rather than the model sparsity. To determine and adapt the modelling size, a novel non-convex framework for gravity inversion is proposed. The proposed optimization aims to directly reduce the L 0 norm of the density matrix. An improved iterative hard thresholding algorithm is developed to linearly reduce the L 0 penalty during the inner iteration. Accordingly, it is possible to determine the modelling scale during the iteration and achieve an expected scale for the reconstructed model. Both simple and complex model experiments demonstrate that the proposed method efficiently reconstructs models. In addition, granites formed during the Yanshanian and Indosinian periods in the Nanling region, China, are reconstructed according to the modelling size evaluated in agreement with the magnetotelluric profile and density statistics of rock samples. The known ores occur at the contact zones between the sedimentary rocks and the reconstructed Yanshanian granites. The ore-forming bodies, periods, and processes are identified, providing guidance for further deep resource exploration in the study area. [ABSTRACT FROM AUTHOR]

A key initial step in geophysical imaging is to devise an effective means of mapping the sensitivity of an observation to the model parameters, that is to compute its Fréchet derivatives or sensitivity kernel. In the absence of any simplifying assumptions and when faced with a large number of free parameters, the adjoint method can be an effective and efficient approach to calculating Fréchet derivatives and requires just two numerical simulations. In the Glacial Isostatic Adjustment problem, these consist of a forward simulation driven by changes in ice mass and an adjoint simulation driven by fictitious loads that are applied at the observation sites. The theoretical basis for this approach has seen considerable development over the last decade. Here, we present the final elements needed to image 3-D mantle viscosity using a dataset of palaeo sea-level observations. Developments include the calculation of viscosity Fréchet derivatives (i.e. sensitivity kernels) for relative sea-level observations, a modification to the numerical implementation of the forward and adjoint problem that permits application to 3-D viscosity structure, and a recalibration of initial sea level that ensures the forward simulation honours present-day topography. In the process of addressing these items, we build intuition concerning how absolute sea-level and relative sea-level observations sense Earth's viscosity structure and the physical processes involved. We discuss examples for potential observations located in the near field (Andenes, Norway), far field (Seychelles), and edge of the forebulge of the Laurentide ice sheet (Barbados). Examination of these kernels: (1) reveals why 1-D estimates of mantle viscosity from far-field relative sea-level observations can be biased; (2) hints at why an appropriate differential relative sea-level observation can provide a better constraint on local mantle viscosity and (3) demonstrates that sea-level observations have non-negligible 3-D sensitivity to deep mantle viscosity structure, which is counter to the intuition gained from 1-D radial viscosity Fréchet derivatives. Finally, we explore the influence of lateral variations in viscosity on relative sea-level observations in the Amundsen Sea Embayment and at Barbados. These predictions are based on a new global 3-D viscosity inference derived from the shear-wave speeds of GLAD-M25 and an inverse calibration scheme that ensures compatibility with certain fundamental geophysical observations. Use of the 3-D viscosity inference leads to: (1) generally greater complexity within the kernel; (2) an increase in sensitivity and presence of shorter length-scale features within lower viscosity regions; (3) a zeroing out of the sensitivity kernel within high-viscosity regions where elastic deformation dominates and (4) shifting of sensitivity at a given depth towards distal regions of weaker viscosity. The tools and intuition built here provide the necessary framework to explore inversions for 3-D mantle viscosity based on palaeo sea-level data. [ABSTRACT FROM AUTHOR]

Atlantic Canada encompasses geological evidence of the orogenic and rifting episodes that inspired the development of the theory of plate tectonics and the fundamental concept of the Wilson cycle. To provide a regional crustal-scale view that can complement surface mapping studies and sparse seismological investigations, an onshore–offshore 3-D constrained gravity inversion methodology is proposed involving incorporation of topography and an inversion mesh that is laterally variable in terms of its maximum depth extent. A 3-D density anomaly model for the entirety of Atlantic Canada and its environs is generated, with the inverted density distribution structure and extracted isodensity surfaces showing excellent correspondence with independent and co-located controlled source and passive seismic constraints. The full density model and crustal thicknesses from this work are made freely available so that they may be used for further study, for instance as inputs for deformable plate reconstruction modelling. [ABSTRACT FROM AUTHOR]

Analysis of tectonic and earthquake-cycle associated deformation of the crust can provide valuable insights into the underlying deformation processes including fault slip. How those processes are expressed at the surface depends on the lateral and depth variations of rock properties. The effect of such variations is often tested by forward models based on a priori geological or geophysical information. Here, we first develop a novel technique based on an open-source finite-element computational framework to invert geodetic constraints directly for heterogeneous media properties. We focus on the elastic, coseismic problem and seek to constrain variations in shear modulus and Poisson's ratio, proxies for the effects of lithology and/or temperature and porous flow, respectively. The corresponding nonlinear inversion is implemented using adjoint-based optimization that efficiently reduces the cost function that includes the misfit between the calculated and observed displacements and a penalty term. We then extend our theoretical and numerical framework to simultaneously infer both heterogeneous Earth's structure and fault slip from surface deformation. Based on a range of 2-D synthetic cases, we find that both model parameters can be satisfactorily estimated for the megathrust setting-inspired test problems considered. Within limits, this is the case even in the presence of noise and if the fault geometry is not perfectly known. Our method lays the foundation for a future reassessment of the information contained in increasingly data-rich settings, for example, geodetic GNSS constraints for large earthquakes such as the 2011 Tohoku-oki M 9 event, or distributed deformation along plate boundaries as constrained from InSAR. [ABSTRACT FROM AUTHOR]

Recently, the eastern margin of the Tibetan plateau has experienced three moderate earthquakes. Deep crustal structure is a critical factor for understanding the seismotectonic environment and deformations in this region. Although several multidisciplinary geophysical approaches have been used to develop crustal structure models for this area, substantial inconsistencies still persist in the results due to the intricate structural properties of crust. In this study, we introduce a novel approach to construct a 3-D crustal model with assimilated density structure (MADS) at a resolution of 0.2° × 0.2° × 5 km. This model is built by integrating data from various source, including multisource gravity anomalies and diverse available crustal structure reference models. We use a model-assimilated gravity inversion process with Bayesian parameter optimization. First, we combine the terrestrial gravity profiles data set with published global gravity field models to create a data set rich in reliable high-frequency anomaly information. Secondly, we incorporate three reference models derived from different seismological methods as prior constraints for our model. These models encompass seismic tomography, surface wave dispersion and receiver function data. We optimize the hyper-parameters of these constraints using the Bayesian criterion. The results demonstrate that the MADS not only captures significant changes in the crustal density but also discerns subtle variations in the upper and middle crust, thereby providing detailed insights into the morphologies of major faults. For instance, the central section of the Longmenshan fault is revealed as a high-angle deep thrust feature, while the frontal section of the Longmenshan fault appears as a low-angle mid-deep thrust feature, and the Xianshuihe fault exhibits a vertical deep subduction feature. Additionally, our findings indicate a correlation between the locations of moderate-to-large earthquakes in this region and the high density-gradient zones or asperities with high density within MADS. We believe that the insights into density characteristics offered by the new MADS model can shed light on the study of asperities associated with recent moderate earthquakes and enhance our understanding of deformation in this region. [ABSTRACT FROM AUTHOR]

With the rise of low-cost and high-density observation system Distributed Acoustic Sensing (DAS), the effective utilization of high-mode surface wave becomes extremely important due to unique measuring method of DAS. To solve the interference of mode identification of dispersion curve and model dimension division on inversion results, we introduced the fitting degree of the dispersion curve, the model dimension, and the uncertainty estimation of the picked dispersion curve to construct a new objective function, and developed a strategy of adaptive trans-dimensional inversion of multimode dispersion curve based on slime mold algorithm (SMA). The research results show that our objective function can not only satisfy the fitting degree of dispersion curve, but also adaptively select the best model dimension, and does not depend on the mode identification of dispersion curve. Inversion strategy based on SMA algorithm has high flexibility, accuracy, stability, and practicality. Our method develops a new technology for dispersion curve inversion and provides a new idea for DAS system to realize low-cost and high-resolution city underground structure detection. [ABSTRACT FROM AUTHOR]

Bayesian inversions followed by estimations of rare event probabilities are often needed to analyze groundwater hazards. Instead of focusing on the posterior distribution of model parameters, the main interest lies then in the distribution of a specific quantity of interest contingent upon these parameters. To address the associated methodological challenges, we introduce a two‐stage Sequential Monte Carlo approach. In the first stage, it generates particles that approximate the posterior distribution; in the second stage, it employs subset sampling techniques to assess the probability of the rare event of interest. By considering two hydrogeological problems of increasing complexity, we showcase the efficiency and accuracy of the resulting PostRisk‐SMC method for rare event probability estimation related to groundwater hazards. We compare the performance of the PostRisk‐SMC method with a traditional Monte Carlo approach that relies on Markov chain Monte Carlo samples. We showcase that our estimates align with those of the traditional method, but the coefficients of variation are notably lower for the same computational budget when targeting more rare events. Furthermore, we highlight that the PostRisk‐SMC method allows estimating rare event probabilities approaching one in a billion using less than one hundred thousand forward simulations. Even if the presented examples are related to groundwater hazards, the methodology is well‐suited for addressing a wide range of topics in the geosciences and beyond. Plain Language Summary: Risk assessment assumes a crucial role in groundwater management. An illustrative example involves the potential contamination of an aquifer, a process influenced by the hydraulic properties of the subsurface. Our focus lies in hazard analysis, where we only have limited or indirect information about these hydraulic properties. In the context of rare event estimation, the central goal is not the estimation of the model parameters, but assessing the probability of a hazard occurrence contingent upon these parameters. When tackling rare events, it is necessary to rely on specialized methodologies to ensure accurate estimation of occurrence probabilities. We explore Sequential Monte Carlo methods integrated with a technique known as subset sampling for this task. In this approach, we iteratively produce posterior samples and shift them toward the critical region with respect to the hazard under consideration. For two hydrogeological examples, the approach demonstrates high accuracy in estimating the probabilities of rare events, even when their occurrence rate approaches one in a billion. Key Points: We introduce a two‐stage Sequential Monte Carlo method generating a particle approximation of rare events under a posterior distributionBy considering two hydrogeological problems, we showcase the efficiency and accuracy of the method in analyzing groundwater hazardsThe presented Sequential Monte Carlo methodology is applicable to many other problems within the geosciences and beyond [ABSTRACT FROM AUTHOR]

The seismic wave velocity in subsurface formations typically varies both spatially and directionally, resulting in a heterogeneous and anisotropic medium. To successfully image complex structures and potential reservoirs using multichannel seismic data, high-resolution interval velocity models are required. However, most controlled-source seismic recordings lack wide-angle signals, such as diving and refracted waves that propagate through deep formations. Consequently, geophysicists rely on the traveltimes and waveforms of reflections or back-scatterings in finite-frequency and finite-offset seismograms to estimate velocity and account for possible anisotropy. Recently, wave-equation based reflection waveform inversion (RWI) has emerged as an active research topic due to its ability to recover intermediate-wavelength model components that cannot be retrieved in reflection traveltime tomography. Nevertheless, the widely used gradient-type local optimization suffers from multiparameter trade-off and slow convergence of RWI in anisotropic media. Analysis of the Fréchet derivative reveals that the sensitivity of reflections to changes in long-to-intermediate wavelengths of the anisotropic velocities is jointly controlled by the specular reflection on the interface and the radiation pattern of scattered fields concerning the model parameters along the wave paths. Anatomy of the approximate Hessian reveals characteristics of parameter coupling and spatial resolution in the context of multiparameter RWI, inspiring the development of a matrix-free Gauss–Newton algorithm based on the second-order adjoint-state method for large-scale applications. Synthetic and real data examples demonstrate that the proposed approach can appropriately handle overburden heterogeneities and anisotropies, and thus improve imaging of underlying complex structures such as dipping and small-scale faults in the deep parts. [ABSTRACT FROM AUTHOR]

As the main contributor to global sea-level rise, the Greenland Ice Sheet (GrIS) has undergone significant mass change over the last two decades. The satellite mission of GRACE (Gravity Recovery And Climate Experiment) and its follow-on mission (GRACE-FO) provide accurate observations but low-spatial resolution. In contrast, satellite altimetry provides observations at a high-spatial resolution but with large uncertainties, limiting the understanding of glacier-scale mass change. To derive accurate and high-spatial resolution mass change estimates from GRACE/GRACE-FO observations, we present a novel constraint mascon method in which the regularization matrix is constructed with the signal variances from satellite altimetry. Based on the proposed method, we derive a series of high-resolution (25 km × 25 km) monthly mascon solutions from 2002 April to September. The glacier-scale estimates from the input–output method agree better with those from our mascon solutions than those from the global mascons of CSR (Center for Space Research, the University of Texas), JPL (Jet Propulsion Laboratory), and GSFC (Goddard Space Flight Center), with a higher linear regression coefficient of 0.71. Benefitting from the greatly improved spatial resolution, our estimates provide the first accurate monthly glacier-scale mass change estimates from GRACE/GRACE-FO observations over the GrIS, to our knowledge. The results show that 20 of the 260 glaciers contributed to more than 42 per cent of the ice loss in the GrIS from 2002 to 2022. Most strikingly, the mass loss of Jakobshavn Isbrae was the most significant at –18.7 ± 0.05 Gt yr−1, accounting for 7.4 per cent of the total in the GrIS during the study period. Furthermore, we find that the SMB (surface mass balance) and ice-dynamics-related mass changes contribute nearly equally to the observed mass changes, but the corresponding spatiotemporal characteristics differ. SMB contributed the most to the mass change in 2019, while ice dynamics played the most dominant role in 2018. Moreover, the SMB effect is widespread, and ice-dynamics-related mass loss is mainly concentrated in five small areas. [ABSTRACT FROM AUTHOR]