Probabilistic formulations of inverse problems are most often based on Bayes Rule, which is considered a powerful tool for integration of data information and prior information about potential solutions. However, since its introduction it has become apparent that the Bayesian inference paradigm presents a number of difficulties, especially in the phase where the problem is mathematically formulated.

Perhaps the most notable difficulty arises because Bayes Theorem is usually formulated as a relation between probability densities on continuous manifolds. This creates an acute crisis because of a problem described by the French mathematician Joseph Bertrand (1889), and later investigated by Kolmogorov and Borel. According to Kolmogorov's (1933/1956) investigations, conditioning of a probability density is underdetermined: In different parameterizations (reference frames), conditional probability densities express different probability distributions. Surprisingly, this problem is persistently neglected in the scientific literature, not least in applications of Bayesian inversion. We will explore this problem and show that it is a serious threat to the objectivity and quality of Bayesian computations including Bayesian inversion, computation of Bayes Factors, and trans-dimensional inversion.

Another difficulty in Bayesian Inference methods derives from the fact that data uncertainties, and prior information on the unknown parameters, are often unknown or poorly known. Because they are required in the calculations, statisticians have invented

hierarchical methods to compute parameters (known as hyper-parameters) controlling these uncertainties. However, since both the data uncertainties and the prior information on the unknowns are supposed to be known 'a priori', but are calculated 'a posteriori', this creates another crisis, namely a violation of causality. We will take a close look at the consequences of this mixing of 'prior' and 'posterior', and show how it potentially jeopardizes the validity of Bayesian computations.

How to cite: Mosegaard, K.: Inconsistency and violation of causality in Bayesian inversion paradigms, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-17474, https://doi.org/10.5194/egusphere-egu23-17474, 2023.

This work discusses the use of full waveform inversion (FWI) with fully nonlinear estimates of uncertainty, to monitor changes in the Earth’s subsurface due to dynamic processes. Typically, FWI is used to produce high resolution 2D and 3D static subsurface images by exploiting information in full acoustic, seismic or electromagnetic waveforms, and has been applied at global, regional and industrial spatial scales. To avoid the over-interpretation of poorly constrained parts of resulting subsurface images or models, it is necessary to know their uncertainty – the range of possible subsurface models that are consistent with recorded data and other pertinent constraints. Almost all estimates of uncertainty on the results of FWI approximate the model-data relationships by linearisation to make the calculation computationally efficient; unfortunately this throws those uncertainty estimates into question, since their raison d’etre is to account for possible model and data variations which are themselves related nonlinearly.

In a related abstract and associated manuscript we use variational inference to achieve the first Bayesian uncertainty analysis for 3D FWI that is fully nonlinear (i.e., involves no linearisation of model-data relationships: https://arxiv.org/abs/2210.03613 ). Variational inference refers to a class of methods that optimize an approximation to the probability distribution that describes post-inversion parameter uncertainties.

Here we extend those methods to perform nonlinear uncertainty analysis for 4D (time-varying 3D) FWI monitoring of the subsurface. Specifically we apply stochastic Stein variational gradient descent (sSVGD) to seismic data generated synthetically for two 3D seismic surveys acquired over a changing 3D subsurface structure based on the 3D overthrust model (Aminzadeh et al., 1997: SEG/EAGE 3-D Modeling Series No. 1). Iterated linearised inversion of each data set fails to image changes (~1%) in the wave speed of the medium, both when each inversion begins independently from the same (good) reference model, or when the best-fit model from inversion of the first survey’s data was used as reference model for the second inversion. Nonlinear inversion of each data set from the same prior distribution also fails to detect these ~1% changes. However, the changes can be imaged and their uncertainty estimated if variational methods applied to invert data from the second survey are initiated from their final state in the inversion of the first survey data. In addition, the methods then converge far more rapidly, compared to running each inversion independently.

We conclude that the probability distributions describing 3D seismic velocity uncertainty are sufficiently complex that the computations of 3D parameter uncertainty for each survey independently have not converged sufficiently to detect small 4D changes. However, the change in these probability distributions between surveys must be sufficiently small that the final solution found from the first survey could evolve robustly into the second survey solution, such that changes are resolved above the uncertainty using variational methods. Nevertheless, this change must be sufficiently complex that linearised methods can not evolve smoothly from one solution to the next, explaining why linearised methods fail, and highlighting why the estimation of nonlinear uncertainties is so important for imaging and monitoring applications.

How to cite: Curtis, A. and Zhang, X.: On Monitoring Changes in the Earth’s Subsurface using 4D Bayesian Variational Full Waveform Inversion, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-14914, https://doi.org/10.5194/egusphere-egu23-14914, 2023.

Determining if uncertainty quantification is worth it or not is closely related to how that uncertainty is computed and the associated computational cost. For seismic imaging, it is typically done using Markov chain Monte Carlo algorithms (McMC). Solving an inverse problem using McMC means exploring and characterizing the ensemble of all plausible models through more or less point-wise random walk in the data misfit landscape. This is typically done using Bayes’ theorem via the computation of a posterior probability density function. Even though this can sound naively simple, it can come with a significant computational burden given the dimension of the problem to be solved and the expense of the forward solver. This is because as the number of dimensions grow, there are exponentially more possible guesses the algorithm can make, while only a few of these models will be accepted as plausible. More advanced uncertainty quantification methods such as Hamiltonian Monte Carlo (HMC) could be beneficial because they can handle higher dimensions because efficient sampling of the model space through pseudo-mechanical trajectories in the data misfit landscape is expected. In order for an HMC algorithm to efficiently sample the model space of interest and provide meaningful uncertainty estimates, three hyper-parameters need to be tuned for trajectory design: the Leapfrog steps L, the Leapfrog stepsize ε, and the Mass Matrix M. There has been already work showing how one can choose L and ε; however designing the appropriate M is far more challenging. We consider a time-lapse seismic scenario and use a local acoustic solver for fast forward solutions. We then use Singular value decomposition, in the vicinity of the true model, to transform our time-lapse optimal model to a system of normal coordinates and use only a few of the eigenvalues and eigenvectors of the Hessian as oscillators. By doing so, we can efficiently understand the impact of the initial conditions and the choice of M and gain insight on how to design M in the standard system. This gives us an intuitive way to understand the mass matrix, allowing us to determine whether gains from the HMC algorithm are worth the cost of determining the parameters.

How to cite: Malcolm, A., Kotsi, M., Ely, G., and Virieux, J.: Is Hamiltonian Monte Carlo (HMC) really worth it? An alternative exploration of hyper-parameter tuning in a time-lapse seismic scenario, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-8679, https://doi.org/10.5194/egusphere-egu23-8679, 2023.

Due to the nonlinearity of inversion as well as the noise in the data, seismic inversion results certainly have uncertainties. Whether quantifying these uncertainties is useful depends at least in part on the computational cost of computing them.  Bayesian techniques dominate uncertainty quantification for seismic inversion.  The goal of these methods is to estimate the probability distribution of the model parameters given the observed data. The Markov Chain Monte Carlo algorithm is widely employed for approximating the posterior distribution. However, generating the posterior samples by combining the prior and the likelihood is intractable for large problems and challenging for smaller problems. We apply a machine learning method called normalizing flows, which consists of a series of invertible and differentiable transformations, as an alternative to the sampling-based methods. In our work, the normalizing flows method is combined with full waveform inversion(FWI) using a numerically exact local solver to quantify the uncertainty of time-lapse changes. We integrate uncertainty quantification(UQ) and FWI by estimating UQ on the images generated by FWI making it computationally practical. In this way, a reasonable posterior probability distribution is directly predicted and produced by transforming from a normal distribution, measuring the amount and spread of variation in FWI images by sample mean and standard deviation. In our numerical results, the method for calculating the posterior distribution of the model is verified to be practical and advantageous in terms of effectiveness.

How to cite: Sun, C., Malcolm, A., and Kumar, R.: Can Normalizing Flows make Uncertainty Quantification Practical for Time-Lapse Seismic Monitoring, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-8767, https://doi.org/10.5194/egusphere-egu23-8767, 2023.

Antarctic subglacial properties impact geothermal heat, subglacial sedimentation, and glacial isostatic adjustment; critical parameters for predicting the ice sheet's response to warming oceans. However, the tectonic architecture of the Antarctic interior is unresolved, with results dependent on datasets or extrapolation used. Most existing deterministic suggestions are derived from qualitative observations and often presented as robust results; however, they hide possible alternative interpretations.

Using information entropy as a measure of certainty, we present a robust tectonic segmentation model generated from similarity analysis of multiple geophysical and geological datasets. The use of information entropy provides us with an unbiased and transparent metric to communicate the ambiguities from the uncertainties of qualitative classifications. Information theory also allows us to test and optimise the methods and data to evaluate how choices impact the distribution of alternative output maps. We further discuss how this metric can quantify the predictive power of parameters as a function of regions with different tectonic settings.

How to cite: Stål, T., Reading, A. M., Cracknell, M. J., Ebbing, J., Halpin, J. A., Kelly, I. D., MacKie, E. J., Sobh, M., Turner, R. J., and Whittaker, J. M.: Using information entropy to optimise and communicate certainty of continental scale tectonic models, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-4074, https://doi.org/10.5194/egusphere-egu23-4074, 2023.

Geophysical inverse problems (seismic tomography) are often significantly underdetermined meaning that a large range of parameter values can explain the observed data well within data uncertainties. Markov chain Monte Carlo (McMC) algorithms based on Voronoi cell parameterizations have been used for quantifying uncertainty in seismic tomography for a number of years. Since surface waves constrain absolute shear velocities and receiver functions (RFs) image discontinuities beneath receiver locations, joint inversion of both data types based on McMC become a popular method to reveal the structure near Earth's surface with uncertainty estimates.

Joint inversion is usually performed in two steps: first invert for 2-D surface wave phase (or group) velocity maps and then invert 1-D surface wave and RFs jointly to construct a 3-D spatial velocity structure. However, in doing so, the valuable information of lateral spatial variations in velocity maps and dipping discontinuities in RFs may not be preserved and lead to biased 3-D velocity structure estimation. Hence, the lateral neighbors in the final 3-D model typically preserve little of the 2-D lateral spatial correlation information in the phase and group velocity maps.

A one-step 3-D direct inversion based on the reversible jump McMC and 3-D Voronoi tessellation is proposed to improve the above issues by inverting for 3-D spatial structure directly from frequency-dependent traveltime measurements and RFs. We take into account the dipping interfaces according to the Voronoi parameterisation, meaning that back azimuth and incidence angle of individual RFs must be taken into account. We present synthetic tests demonstrating the method. Individual inversion of surface wave measurements and RFs show the limitation of inverting the two data sets separately as expected: surface waves are poor at constraining discontinuities while RFs are poor at constraining absolute velocities. The joint solution gives a better estimate of subsurface properties and associated uncertainties. Compared to two-step inversion which may produce bias propagating between two steps and lose valuable lateral structure variations, the direct 3-D direct inversion not only produces more intuitively reasonable results but also provides more interpretable uncertainties.

How to cite: Ke, K.-Y., Tilmann, F., Ryberg, T., and Mroczek, S.: 3-D joint inversion of surface wave and receiver functions based on the Markov chain Monte Carlo, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-5661, https://doi.org/10.5194/egusphere-egu23-5661, 2023.

Determining the thermochemical structure of the mantle is crucial for understanding its evolution and dynamics. Temperature variations have long been known as important driving forces of mantle convection; however compositional differences can also influence dynamics. Additionally, compositional differences can act as indicators left behind by processes operating in the past. Both aspects have played a role in the ongoing discussions on the Large Low Shear Wave Velocity Provinces (LLSVP), the proposed Bridgmanite Enriched Ancient Mantle Structures (BEAMS) and the fate of subducted oceanic crust.

A prerequisite for determining compositional differences in terms of major oxides with geophysical techniques is a joint determination of several geophysical properties. A single geophysical property (density, velocity) could almost always be explained by temperature or composition variations alone – except in pathological edge cases. The geophysical signature of composition lies in the pointwise relation between properties. This pointwise relation can be distorted by spectral filtering or inversion smoothing and damping.

In this contribution, I parametrize the mantle as a collection of discrete spatial anomalies in terms of seismic velocity and density. Surface wave phase speed and satellite gravity data are used to constrain the anomalies. A transdimensional Monte Carlo Markov Chain method is used to generate ensembles of solutions that try to balance model complexity and data fit. An important aspect of this setup is that the two data sets used are complementary: While satellite gravity data are available (nearly) globally with homogeneous quality, coverage of phase speed data depends on the spatial distribution of seismic stations and large earthquakes. Conversely, the gravity field lacks true depth sensitivity, which surface wave data can provide by combining several frequencies.

I will present synthetic investigations that aim at determining how accuracy and coverage affect the simultaneous recoverability of seismic velocity and density.

How to cite: Szwillus, W.: Sensitivity of surface wave and gravity data to velocity and density structure in the mantle – insights from transdimensional inversion, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-15486, https://doi.org/10.5194/egusphere-egu23-15486, 2023.

Precise determination of earthquake hypocenter (longitude, latitude and depth) and its origin time is of fundamental importance for not only understanding the seismogenic process but also revealing the Earth’s interior structure. Instrumental coverage plays the first-order role in determining earthquake locations. For earthquakes that occurred in the continental interior, it is favorable to have seismic stations with full azimuthal coverage; nonetheless, precise determination of earthquake depth is often challenging due to its tradeoff with earthquake origin time. The situation is even worse for earthquakes that occurred in offshore regions, e.g., Pacific ring of fire, because regional seismic stations are mostly installed on the continent. To deal with those challenges aforementioned, we propose to constrain the earthquake hypocenter by jointly using first arrivals (P and S waves) and depth phase traveltimes. The theoretical travelling times of these phases are precisely and efficiently calculated in 3D velocity model through solving the Eikonal equation. Once the earthquake hypocenter is well constrained, we further improve the accuracy of the origin time. We tested and verified the proposed earthquake location strategy in the Ridgecrest area (southern California), which serves as an end member of continental setting, and central Chile, which serves as another end member of offshore setting. The station coverage is complete in the Ridgecrest area. We have identified and picked first arrivals and sPL phases at local distances. On the contrary, seismic stations are only installed on the continent in central Chile. We have identified and picked first arrivals and sPn phases at regional distances. Determined earthquakes have comparable location accuracy as the regional catalog in the horizontal plane, while the depth uncertainty has been reduced greatly. Our study shows that incorporating depth phases into the earthquake location algorithm together with first arrivals can greatly increase earthquake location accuracy, especially earthquake depth, which will lay the solid foundation for wide-scope topics in earth science studies.

How to cite: Li, T., Chen, J., and Tong, P.: Locating earthquake hypocenter using first arrivals and depth phase in 3D model at local and regional distances, EGU General Assembly 2023, Vienna, Austria, 23–28 Apr 2023, EGU23-2190, https://doi.org/10.5194/egusphere-egu23-2190, 2023.