Efficient Markov chain Monte Carlo (MCMC) sampling from posterior distributions remains a central challenge in Bayesian geophysical inversion. Recent developments in computational statistics and optimal transport suggest that MCMC efficiency can be improved by reparametrizing the sampling problem-specifically, by learning an invertible mapping that recasts the target distribution onto a simpler reference distribution. Here, we introduce a Metropolis-Hastings framework that leverages transport maps parametrized by invertible neural networks. These maps are trained on preliminary MCMC samples from the target distribution and used to propose new samples in a fixed reference space, where proposal design is independent of the target's structure. The proposed samples are transformed back to the target space via the inverse map, and accepted or rejected according to a modified Metropolis-Hastings criterion. As sampling proceeds, the transport maps are updated, yielding proposals increasingly well adapted to the shape of the target distribution. Across a suite of numerical tests-including a 2-D Rosenbrock distribution, a 3-D earthquake location problem and Gaussian mixtures up to 16 dimensions-transport-map-driven samplers consistently outperform standard MCMC, reducing integrated autocorrelation times by factors of 2.5 to over 6 (or equivalently, yielding sample sets 2.5-6 times larger for the same number of forward evaluations). This improvement comes at the non-negligible cost of training one or more transport maps, which we quantify systematically. We also provide a quantitative criterion for weighing training cost against sampling speed-up. This shows that transport-map MCMC is advantageous whenever the forward problem is non-trivial, making it a promising approach for Bayesian sampling in geophysics and beyond.

Inference problems within the geosciences vary considerably in terms of size and scope, ranging from the detection of changepoints in 1D time/depth models, to the construction of complex 3D or 4D models of the Earth. Solving an inverse problem typically requires fusing various classes of data, each associated with its own forward model. The choice of an appropriate inference method is itself not obvious. An investment of much time effort and is required in software development and education. Many researchers have developed bespoke inversion and parameter estimation algorithms tailored to their specific needs. Associated software is then typically bespoke to the particular application, often requiring significant investment by new researchers to master with minimal documentation. This is entirely understandable as generalisation and ongoing support of inference codes requires significant time and effort that is frequently beyond the primary objectives of the research. As a result the all important experimentation often required to choose an appropriate inversion method for a new data set or domain, is often not practical. Furthermore design choices made in existing software implementations often dictate those by subsequent researchers and influence the scientific direction taken. The Common Framework for Inference, CoFI, is an open source project which aims to capture inherent commonalities present in all types of inverse problems, independent of the specific methods employed to solve them. CoFl is codifies the definition of an inference problem and then provides an interface to reliable and sophisticated third-party packages, such as SciPy and PyTorch, to tackle inverse problems across a broad range. The modular and object-oriented design of CoFI, supplemented by a comprehensive suite of tutorials and practical examples, ensures its accessibility to users of all skill levels, from experts to novices. This not only has the potential to streamline research and promote best practice but also to support education and STEM training. This poster gives an overview of CoFl through domain relevant examples, from optimisation to probabilistic sampling. With a focus on CoFI’s modular approach we hope to foster collaboration centred around interaction by expanding the set of inference algorithms and domain-relevant examples.

Geophysical inversion often involves choices and assumptions. Which combination of forward model, noise model, regularisation, and solution algorithm is best suited to a given dataset is rarely obvious in advance. Practitioners typically design inversion schemes, and software, around a single class of method based on their choices. In these circumstances, moving from a parameter estimation approach to an ensemble or Bayesian one requires substantial recoding effort. This discourages cross-fertilization of methods between sub-fields that share a common mathematical structure.With these issues in mind we present CoFI, the Common Framework for Inference, an open source Python package that separates inverse problem definition from algorithm selection. CoFI provides three objects (classes): BaseProblem, which captures the forward model, data, noise model, and any regularization or prior; InversionOptions, which specifies the algorithm; and Inversion, which executes the inversion. Because problem and algorithm are held separately, switching between solvers requires only a change to InversionOptions, with BaseProblem unchanged. CoFI connects user-written forward models to existing solver tools and software, including SciPy, emcee and BayesBay, without re-implementing inference methods.We illustrate CoFI through three worked examples. For the educator, a polynomial regression and sea-level change-point problem introduces three solution styles, from direct matrix solvers through trans-dimensional Bayesian sampling, within a single teaching sequence. For those who carry out inversion studies, an ambient-noise surface-wave tomography of Australia demonstrates progression from regularised linear inversion through nonlinear iterative methods to Bayesian uncertainty quantification. For those who develop inference methods, a Himmelblau function example shows how to register a new inference algorithm and integrate it into the standard CoFI interface.CoFI is accompanied by a growing library of worked examples and tutorials. We invite contributions from researchers, developers, and teachers, and hope that the framework will become a shared resource for the geoscientific community.

Trans-dimensional Bayesian sampling has been applied to subsurface imaging and other inference problems across the Earth Sciences. A particular style of Markov chain Monte Carlo (McMC) method, known as reversible-jump has been used almost universally in such studies. This algorithm allows sampling across variably dimensioned model parameterizations. However, for practical reasons, it is limited to cases where the number of free parameters differ in a regular sequence between alternate models, usually by addition or subtraction of a single variable. Furthermore, jumps between model dimensions rely on bespoke mathematical transformations, which are bespoke to each class of application. As a result, implementations are dependent on the choice of model parameterization employed. A framework for Trans-conceptual Bayesian sampling, which is a generalization of trans-dimensional sampling, is presented. Trans-C Bayesian sampling allows exploration across a finite, but arbitrary, set of conceptual models, that is ones where the number of variables, the type of model basis function, nature of the forward problem, and assumptions on the measurement noise statistics, may all vary independently. The new framework avoids parameter transformations and thereby lends itself to development of automatic McMC algorithms, that is where the details of the sampler do not require knowledge of the parameterization. Algorithms implementing Bayesian conceptual model sampling are presented and illustrated with examples drawn from geophysics, using real and synthetic data. Comparison with reversible-jump illustrates that trans-C sampling produces statistically identical results for situations where the former is applicable, but also allows sampling in situations where trans-D would be impractical to implement.

Over the past twenty years, Trans-dimensional Bayesian Inference has become a popular approach for Bayesian sampling. It has been applied widely in the geosciences when the best class of model representation, e.g. of the subsurface, is not obvious in advance, or the number of free variables undecided. Making arbitrary choices in these areas may result in sub-optimal inferences from data. In trans-D, one typically defines a finite number of model states, with differing numbers of unknowns, over which Bayesian Inference is to be performed using the data.A key attraction of Trans-D Bayesian Inference is that it is designed to let the data decide which state, as well as which configurations of parameters within each state, are preferred by the data, in a probabilistic manner. Trans-D algorithms may hence be viewed as a combination of fixed dimensional within-state sampling and simultaneous between-state sampling where Markov chains visit each state in proportion to their support from the data.In theory, each state may be completely independent, involving different classes of model parameterization, with different numbers of unknowns, data noise levels, and even different assumptions about the data-model relationship. Practical considerations, such as convergence of the finite length Markov chains between states, usually mean that each state must be closely related to each other, e.g. differing by a single layer in a 1-D seismic Earth model. In addition, since the form of the necessary Metropolis-Hastings balance condition depends on the mathematical relationship between the unknowns in each state, then implementations are often bespoke to each class of model parameterization and data type. To our knowledge there exists no automatic trans-D sampler where one can define arbitrary independent states, together with a prior and Likelihood, and simply pass to a generalised sampling algorithm, as is common with many fixed dimensional MCMC algorithms and software packages. A second limitation in trans-D sampling is that since implementations are bespoke within a class of model parameterizations, within-state sampling is typically performed with simplistic and often dated algorithms, e.g. Metropolis-Hastings or Gibbs samplers, thereby limiting convergence rates. Over the past 30 years fixed dimensional sampling has advanced considerably with numerous efficient algorithms available and many conveniently translated into user friendly software packages, almost all of which have not been used within a trans-D framework due to a lack of a way to conveniently deploy them in a trans-D setting.In this presentation we will address all of these issues by describing the theory under-pinning an ‘Independent State’ (IS) Trans-D sampler, together with some illustrative examples. In this algorithm class, sampling may be performed across states that are completely independent, containing arbitrary numbers of unknowns and parameter classes. In addition, the IS-sampler can conveniently take advantage of any fixed dimensional sampler without the need to derive and re-code bespoke Markov chain balance conditions, or specify mechanisms for transitions between model parameters within different states. In this sense it represents a general purpose automatic trans-D sampler.

Mineral exploration is frequently centred around delineating discrete geological units with, typically, sharp boundaries that could represent economic targets. In the case of complex resistivity (CR) inversions, the choice of regularization and model parametrization significantly impacts the inversion's ability to delineate targets. Initially, however, a prudent researcher may not wish to bias their inversion towards sharp distinct units without prior justification. Here, we explore how a suite of regularization approaches to the CR inverse problem allows to encompass different classes of prior beliefs. We present these as a progression as more information becomes available regarding the likelihood of distinct geological units. The most weakly informed approach with respect to the delineation of geological units we consider is the classic $\ell _2$-type regularization, tend to produce smeared-out fuzzy images. However, this is typically not what is expected for distinct geological units, and we compare this with schemes that increasingly resolve sharp boundaries. We test a range of $\ell _1$-type regularizations, which have been frequently touted in the geophysics and optimization literature as being well suited for such tasks. We experiment with using a so-called overcomplete parametrization of the CR field, which aims to separate smooth background and sharp foreground features. These $\ell _1$ schemes are shown to produce generally sharper images than $\ell _2$. In the most informed case, where strong assumptions can be made about the local geology, we represent the CR field as a foreground ellipse in a homogenous background. This approach significantly reduces the size of the parameter space, and tends to have a simple geometric interpretation. While the anomaly parametrization has some unique challenges, we show it clearly resolves distinct units compared to both the $\ell _2$ and $\ell _1$ regularizations. Applications first to synthetic data and then to field data from Century Zinc Deposit in northern Australia, demonstrate the progression from weakly informed to strongly informed regularization and parametrization and the sharpness of the recovered geological units.

The lack of versatile tools for Bayesian inference presents a significant challenge to researchers in geophysics, who often resort to developing bespoke codes to address specific classes of inverse problems. In this study, we present BayesBay, a Python package for generalized transdimensional and hierarchical Markov chain Monte Carlo sampling. Leveraging object-oriented programming principles, BayesBay facilitates the definition of Bayesian sampling problems across a range of applications. This includes joint inversions of multiple data sets with different forward functions and unknown noise properties, as well as complex parameterizations involving multiple parameters with unknown dimensionality and/or spatially varying priors. We illustrate BayesBay from both a technical and a practical perspective. The first two applications are common in geophysics: a 2D tomographic problem and a joint inversion for the 1D subsurface structure. The third involves partition modeling and requires a sophisticated parameterization with two nested levels of transdimensionality. In all cases, BayesBay recovers known solutions, highlighting its potential to address a broad range of inverse problems.

The Earth's interior structure must be inferred from geophysical observations collected at the surface. Compared to just a few decades ago, the amount of geophysical data available today is voluminous and growing exponentially. Dense seismic networks like USArray, AlpArray, and AusArray now enable joint inversions of various geophysical data types to maximise subsurface resolution at scales ranging from local to continental. However, the practical application of joint inversions faces several challenges: Various geophysical techniques typically probe different scales and depths, complicating the choice of an appropriate discretisation for the Earth's interior. Different geophysical observables may respond to physical properties that are not directly related (e.g., density and electrical conductivity), making the construction of self-consistent parameterisations a non-trivial task. Without a comprehensive understanding of noise characteristics, standard methods require assigning weights to different data sets, yet robust choices remain elusive. Capable of overcoming these recognised challenges and allowing estimates of model uncertainty, probabilistic inversions have grown in popularity in the geosciences over the last few decades, and have been successfully applied to specific modelling problems. Here, we present BayesBridge, a user-friendly Python package for generalised transdimensional and hierarchical Bayesian inference. Computationally optimised through Cython, our software offers multi-processing capabilities and runs smoothly on both standard computers and computer clusters. As opposed to existing software libraries, BayesBridge provides high-level functionalities to define complex parameterisations, with prior probabilities (defined by uniform, Gaussian, or custom density functions) that may or may not be dependent on depth and/or geographic coordinates. By default, BayesBridge employs reversible-jump Markov chain Monte Carlo for sampling the posterior probability, with the option of parallel tempering, but its low-level features enable effortless implementations of arbitrary sampling criteria. Utilising object-oriented programming principles, BayesBridge ensures that each component of the inversion -- such as the discretisation, the physical properties to be inferred, and the data noise -- is a self-contained unit. This design facilitates the seamless integration of various forward solvers and data sets, promoting the use of multiple data types in geophysical inversions.

Inference problems within the geosciences vary significantly in size and scope, ranging from the detection of data trends through simple linear regressions, to the construction of complex 3D models representing the Earth’s interior structure. Successfully solving an inverse problem typically requires combining various types of data sets, each associated with its own forward solver. In the absence of established software, many researchers and practitioners resort to developing bespoke inversion and parameter estimation algorithms tailored to their specific needs. However, this practice does not promote reproducibility and necessitates a substantial amount of work that is frequently beyond the primary objectives of the research. Our aim with CoFI (pronounced: coffee), the Common Framework for Inference, is to capture inherent commonalities present in all types of inverse problems, independent of the specific methods employed to solve them. CoFI is an open-source Python package that provides a link to reliable and sophisticated third-party packages, such as SciPy and PyTorch, to tackle inverse problems of a broad range. The modular and object-oriented design of CoFI, supplemented by our comprehensive suite of tutorials and practical examples, ensures its accessibility to users of all skill levels, from experts to novices. This not only has the potential to streamline research but also to support education and STEM training. This poster presentation aims to give an overview of CoFI’s main features and usage through practical examples. Moreover, we hope to foster collaboration and invite contributions on inference algorithms and domain-relevant examples.

<p>Many Earth systems cannot be observed directly, or in isolation. Instead, we must infer their properties and characteristics from their signature in one or more datasets, using a variety of techniques (including those based on optimization, statistical methods, or machine learning). Development of these techniques is an area of focus for many geoscience researchers, and methodological advances can be instrumental in enhancing our understanding of the Earth.         </p> <p>In our experience, progress is substantially hindered by the absence of infrastructure facilitating communication between sub-disciplines. Researchers tend to focus on one area of the earth sciences — such as seismology, hydrology or oceanography — with only slow percolation of ideas and innovations from one area to another. Indeed, silos often exist even within these subfields. Testing new ideas on new problems is challenging as it requires the acquisition of domain knowledge, an often difficult and time-consuming endeavour with uncertain returns. Key questions that arise include: What is a relevant field data set, and how has it been processed? Which simulation package is most appropriate to predict the data? What would a 'good' model look like and what should it be able to resolve? What is the current best practice?</p> <p>To address this, we introduce the ESPRESSO project — a collection of Earth Science Problems for the Evaluation of Strategies, Solvers and Optimisers. It aims to provide  access to a suite of ‘test problems’, spanning a wide range of inference and inversion scenarios. Each test problem defines appropriate dataset(s) and simulation routines, accessible within a standardised Python interface. This will allow researchers to rapidly test new techniques across a spectrum of problems, share domain-specific inference problems and ultimately identify areas where there may be potential for fruitful collaboration and development. ESPRESSO is envisaged as an open, community-sourced project, and we invite contributions from across the geosciences.</p>

Regularized least-squares tomography offers a straightforward and efficient imaging method and has seen extensive application across various fields. However, it has a few drawbacks, such as (i) the regularization imposed during the inversion tends to give a smooth solution, which will fail to reconstruct a multi-scale model well or detect sharp discontinuities, (ii) it requires finding optimum control parameters, and (iii) it does not produce a sparse solution. This paper introduces ‘overcomplete tomography’, a novel imaging framework that allows high-resolution recovery with relatively few data points. We express our image in terms of an overcomplete basis, allowing the representation of a wide range of features and characteristics. Following the insight of ‘compressive sensing’, we regularize our inversion by imposing a penalty on the L1 norm of the recovered model, obtaining an image that is sparse relative to the overcomplete basis. We demonstrate our method with a synthetic and a real X-ray tomography example. Our experiments indicate that we can reconstruct a multi-scale model from only a few observations. The approach may also assist interpretation, allowing images to be decomposed into (for example) ‘global’ and ‘local’ structures. The framework presented here can find application across a wide range of fields, including engineering, medical and geophysical tomography.

SUMMARY Concerns raised by Okazaki & Ueda (2022) on the paper by Sambridge et al. (2022) are addressed. Two issues are discussed and some new numerical results presented. The first concerns whether the properties of the Wasserstein time-series misfit introduced in our earlier paper will translate to model space non-uniqueness in a seismic waveform inversion setting. It is argued that this is unlikely, given the special conditions, which must exist between all observed/predicted seismic waveform pairs for non-uniqueness to result. The second issue discussed is the efficacy of using the Sliced Wasserstein algorithm of Bonneel et al. (2015) as an alternate to the marginal Wasserstein algorithm, as proposed by Okazaki & Ueda (2022). It is argued that for optimization-based waveform fitting, the Sliced Wasserstein algorithm is a viable alternate provided care is taken to ensure that conditions arise which do invalidate analytical derivative expressions of the resulting Wasserstein misfit. In practice, this would likely mean recasting the 2D Optimal Transport problem posed in our earlier paper onto unstructured grids.

We present a theoretical framework that links Fermat’s principle of least time to optimal transport theory via a cost function that enforces local transport. The proposed cost function captures the physical constraints inherent in wave propagation; when paired with specific mass distributions, it yields shortest paths in the considered media through the optimal transport plans. In the discrete setting, our formulation results in physically significant optimal couplings, whose off-diagonal entries identify shortest paths in both directed and undirected graphs. For undirected graphs with positive edge weights, commonly used to parameterize seismic media, our method provides solutions to the Eikonal equation consistent with those from the Dijkstra algorithm. For directed negative-weight graphs, corresponding to transportation cost matrices with negative entries, our approach aligns with the Bellman–Ford algorithm but offers considerable computational advantages. We also highlight potential research directions. These include the use of sparse cost matrices to reduce the number of unknowns and constraints in the considered transportation problem, and solving specific classes of optimal transport problems through the Dijkstra algorithm to enhance computational efficiency.

In this chapter, we survey some recent developments in the field of geophysical inversion. We aim to provide an accessible general introduction to the breadth of current research, rather than focussing in depth on particular topics. In particular, we hope to give the reader an appreciation for the similarities and connections between different approaches, and their relative strengths and weaknesses.

When working with hydrological data, the ability to quantify the similarity of different datasets is useful. The choice of how to make this quantification has a direct influence on the results, with different measures of similarity emphasising particular sources of error (for example, errors in amplitude as opposed to displacements in time and/or space). The Wasserstein distance considers the similarity of mass distributions through a transport lens. In a hydrological context, it measures the “effort” required to rearrange one distribution of water into the other. While being more broadly applicable, particular interest is paid to hydrographs in this work. The Wasserstein distance is adapted for working with hydrographs in two different ways and tested in a calibration and “averaging” of a hydrograph context. This alternative definition of fit is shown to be successful in accounting for timing errors due to imprecise rainfall measurements. The averaging of an ensemble of hydrographs is shown to be suitable when differences among the members are in peak shape and timing but not in total peak volume, where the traditional mean works well.

We propose a new approach to measuring the agreement between two oscillatory time series, such as seismic waveforms, and demonstrate that it can be employed effectively in inverse problems. Our approach is based on Optimal Transport theory and the Wasserstein distance, with a novel transformation of the time series to ensure that necessary normalisation and positivity conditions are met. Our measure is differentiable, and can readily be employed within an optimization framework. We demonstrate performance with a variety of synthetic examples, including seismic source inversion, and observe substantially better convergence properties than achieved with conventional $L_2$ misfits. We also briefly discuss the relationship between Optimal Transport and Bayesian inference.

SUMMARY Monte Carlo methods are widespread in geophysics and have proved to be powerful in non-linear inverse problems. However, they are associated with significant practical challenges, including long calculation times, large output ensembles of Earth models, and difficulties in the appraisal of the results. This paper addresses some of these challenges using generative models, a family of tools that have recently attracted much attention in the machine learning literature. Generative models can, in principle, learn a probability distribution from a set of given samples and also provide a means for rapid generation of new samples which follow that approximated distribution. These two features make them well suited for application to the outputs of Monte Carlo algorithms. In particular, training a generative model on the posterior distribution of a Bayesian inference problem provides two main possibilities. First, the number of parameters in the generative model is much smaller than the number of values stored in the ensemble, leading to large compression rates. Secondly, once trained, the generative model can be used to draw any number of samples, thereby eliminating the dependence on an often large and unwieldy ensemble. These advantages pave new pathways for the use of Monte Carlo ensembles, including improved storage and communication of the results, enhanced calculation of numerical integrals, and the potential for convergence assessment of the Monte Carlo procedure. Here, these concepts are initially demonstrated using a simple synthetic example that scales into higher dimensions. They are then applied to a large ensemble of shear wave velocity models of the core–mantle boundary, recently produced in a Monte Carlo study. These examples demonstrate the effectiveness of using generative models to approximate posterior ensembles, and indicate directions to address various challenges in Monte Carlo inversion.

Summary The package pyprop8 enables calculation of the response of a 1-D layered halfspace to a seismic source, and also derivatives (‘sensitivity kernels’) of the wavefield with respect to source parameters. Seismograms, seismic spectra, and measures of static displacement (e.g. GPS, InSAR and field observations) may all be simulated. The method is based on a Thompson-Haskell propagator matrix algorithm, described in O’Toole & Woodhouse (2011) and O’Toole et al. (2012). The package is entirely written in Python, dependent only on the mainstream libraries numpy (Harris et al., 2020) and scipy (Virtanen et al., 2020). As such, it is lightweight and easy to deploy across a variety of platforms, making it particularly suited to use for teaching and outreach purposes.

The relationship between Delta S-36 and Delta S-33 in Archean sedimentary pyrites has been used to evaluate early geologic processes, including photochemical reactions in the anoxic atmosphere, biological activity and thermochemical alteration during sediment deposition. We have applied statistical methods to quadruple S isotope analyses of Archean sedimentary pyrites, using data compiled from the literature. Most of the best-fit lines, on plots of Delta S-36 against Delta S-33, have Archean reference array-like Delta S-36/Delta S-33 slopes that vary between -1.5 and - 0.9. Rigorous statistical tests were conducted to calculate the probability of the best-fit lines passing through the origin. Seventeen of 23 Delta S-36-Delta S-33 regression lines, which pass our reliability filter of R-2 >= 75% and Delta S-33 range >= 2 parts per thousand, have positive intercepts on the Delta S-36 axis, and 13 of these have a probability of < 5% of a zero intercept on the Delta S-36 axis. The observed Delta S-36/Delta S-33 slopes and the non-negative intercepts, which requires at least two mass-independent fractionation source reactions to operate simultaneously, can be produced by UV radiation in the atmosphere at low SO2 partial pressures by combining collision-induced intersystem crossing in the SO2 photoexcitation band (240-340 nm), with the self-shielding effect in the SO2 photolysis band (190-220 nm). The two SO2 photochemical processes must occur simultaneously in a single atmospheric reservoir in order that the fraction contributed by the end-member process remains constant across the full range of Delta S-33 values. We call this process simultaneous fractionation. We applied a two-end-member model to calculate the fraction of S contributed by the SO2 photoexcitation end-member (f) needed to produce the observed Delta S-36/Delta S-33 gradients and variable intercepts on the Delta S-36 axis in the Archean sedimentary pyrites, when the other end-member is SO2 photolysis with the self-shielding. The simplest explanation for variations in f, and therefore variations in Delta S-36/ Delta S-36 gradients, is that it is controlled by changes in the partial pressure of SO2 in the atmosphere.