Keywords: . . . . . . . . . 4.1 Convergence statistics for probabilistic Inverse theory, Probability distributions, Tomography, Interfaces. xiii

field data inversion. . . . . .

Résumé

De nombreux processus physiques, chimiques et biologiques du sous-sol ont lieu proche de la surface, au niveau des interfaces où les propriétés en jeu sont contrastées. Ainsi, la caractérisation de la profondeur et de la géométrie de telles interfaces ainsi que leurs incertitudes revêtent un intérêt majeur dans de nombreux domaines des sciences de la Terre. Les méthodes géophysiques sont intrinsèquement sensibles aux discontinuités dans les propriétés physiques souterraines, mais l'interprétation des observations recueillies et la reconstruction d'une image du sous-sol sont difficiles. La principale limitation provient du fait que les problèmes d'inversion sont mal posés. L'inversion déterministe résout cette problématique par régularisation, le plus souvent en utilisant des contraintes de lissage qui gomment toutes les interfaces naturellement présentes. L'inversion probabiliste modélise plutôt la solution inconnue comme une variable aléatoire en la décrivant par sa fonction de densité de probabilité à posteriori. Dans ce cadre, la distribution de probabilité à priori a une forte influence sur les réalisations des modèles a posteriori et sur la complexité du problème. De plus, sa formulation lorsque peu de connaissance à priori est disponible n'est pas anodine. Nous proposons ici une formulation (et une solution) probabiliste au problème d'inversion d'un ou plusieurs jeux de données géophysiques pour inférer les interfaces en présence de sous-domaines hétérogènes, lorsque les connaissances préalables sont rares. Dans le but de proposer une formulation alternative pour une distribution de probabilités à priori "non informative" en cas de paramétrage par discrétisation spatiale, nous développons un algorithme d'échantillonnage qui suppose une distribution à priori uniforme sur les mesures de variabilité spatiale du modèle, au lieu du choix classique des distributions log-uniformes et non corrélées pour les paramètres du modèle. Nous démontrons la capacité de cette technique d'inversion à échantillonner des réalisations et des statistiques de modèle postérieur satisfaisantes sur des jeux de données synthétiques et basés sur le champ. La méthode est théoriquement valable, mais sa mise en oeuvre numérique est limitée aux cas dans lesquels la valeur des paramètres à priori est limitée à des gammes relativement étroites. Ainsi, dans notre implémentation qui sépare les mises à jour des paramètres géométriques (interfaces) et physiques (valeurs de propriétés physiques), basé selon un échantillonnage de Gibbs, nous contraignons les propriétés géométriques et physiques de manière à favoriser des transitions spatiales graduelles en utilisant les méthodes empiriques de Bayes. Nous étendons enfin l'algorithme empirical-Bayes-within-Gibbs développé pour inverser conjointement plusieurs jeux de données géophysiques afin de réduire l'ambiguïté inhérente à l'interprétation de mesures individuelles. Le couplage entre les modèles capables d'expliquer les différents jeux de données géophysiques est réalisé en considérant une interface commune. Des test synthétiques et basés sur le terrain démontrent que la méthode échantillonne l'interface cible de manière plus précise que les résultats obtenus par inversion déterministe et que les résultats sont encore améliorés par l'inversion conjointe de deux jeux de données. Background and Motivation

Imaging the Earth

The imaging and exploration of the Earth's shallow subsurface is of great importance for human society. The pressure on natural resources is increasing (in terms of water supply, fossil fuels, hydrocarbons, geothermal reservoirs, etc.) and it is forecast to continue growing as the earth's population increases (United Nations , 2017). To aid in the study of these resources and the impact of human activity on them, geophysics can provide useful information through remotely sensed data that are sensitive to different subsurface physical properties. This enables the study of shallow structures for engineering purposes, the exploration of economically useful geological deposits, or the management of near surface biological, chemical and physical processes that affect local environments. The popularity of geophysical measurements lies particularly in the large spatial coverage and their ability to see through, below, and into solid materials using non-invasive (or minimally invasive) techniques. This leads to important applications in a variety of fields where investigations may otherwise require intrusive or destructive methods.

Furthermore, the inherent sensitivity of geophysical observations to discontinuities in the physical properties of the Earth makes geophysics a very useful tool for mapping subsurface boundaries. These boundaries are of increasing interest in different fields of the Earth sciences. This can be seen, for instance, in the characterization of the interface between saltwater and freshwater in coastal areas, where its position and dynamics are needed to study and predict mixing processes [START_REF] Bear | Seawater Intrusion in Coastal Aquifers-Concepts, Methods and Practices[END_REF][START_REF] Werner | Seawater intrusion processes, investigation and management: Recent advances and future challenges[END_REF]. Other examples include the discontinuity between frozen and unfrozen ground for permafrost monitoring purposes [START_REF] Kneisel | Advances in Geophysical Methods for Permafrost Investigations[END_REF] and the boundary between regolith and intact bedrock, which is of crucial importance for landslide studies [START_REF] Lanni | Simulated effect of soil depth and bedrock topography on near surface hydrologic response and slope stability[END_REF] and for the characterization of the Earth's critical zone (the interface between geosphere, atmosphere, biosphere and hydrosphere that provides the ecosystem for most of the life on Earth; Brantley et al., 2011).

Inferring the causes of observed measurements

Once the geophysical data have been collected, the information contained in these needs to be interpreted, which implies inferring the unknown causes for the observed indirect measurements. Inversion attempts to address this problem by estimating parameters describing the investigated physical system that explain the observations [START_REF] Mosegaard | Probabilistic Approach to Inverse Problems[END_REF]. The use of the term inversion stresses the need to reverse the natural direction of causality of predictive models (i.e., forward models), that are generally well-posed and therefore characterized by a unique and stable solution (small changes in the causes lead to small changes in the consequences; [START_REF] Hadamard | Sur les problèmes aux dérivées partielles et leur signification physique[END_REF].

Given a physical system, described by a set of M model parameters, m = {m 1 , m 2 , • • • , m M }, the forward or direct problem consists of using a physical theory to predict the N theoretical data values,

d si m = {d si m 1 , d si m 2 , • • • , d si m
N }, that we would observe as outcome of the geophysical measurements. The (generally) nonlinear forward operator, g(• ), consists of a set of equations that causally associate the observations d to a model m: g(m) = d.

(

The forward model, together with the choice of parameterization of the subsurface, allows for a mathematical description of the investigated system that can then be used for testing predictions [START_REF] Gupta | Reconciling theory with observations: Towards a diagnostic approach to model evaluation[END_REF]. Such quantitative models are generally a simplification of reality, partly because of lack of knowledge about how to express the processes characterizing the subsurface, but also because of conscious strategies aiming to decrease computational efforts. These simplifications are based on assumptions regarding the importance of different subsurface mechanisms on the observations. In order to account for these simplifications and the inherent uncertainties in the observed measurements, the forward problem can be reformulated as:

g(m) = d + e, (1.2)
where e is a vector of dimension N, containing different sources of uncertainties [START_REF] Arendt | Quantification of Model Uncertainty: Calibration, Model Discrepancy, and Identifiability[END_REF]. In addition to measurement noise, this term should also take into account model discrepancies, which results from inadequately describing the true physical system [START_REF] Brynjarsdóttir | Learning about physical parameters: The importance of model discrepancy[END_REF], and modeling error, which typically refers to the result of working with a model of lower complexity and accuracy than the sometimes hypothetical high-fidelity natural simulator that supposedly produced the data [START_REF] Bayarri | A framework for validation of computer models[END_REF]. The importance of modeling, and especially the effect of a chosen model prediction strategy on the estimation of parameters, has been demonstrated in several studies [START_REF] Refsgaard | A framework for dealing with uncertainty due to model structure error[END_REF][START_REF] Bond | What do you think this is? "Conceptual uncertainty"in geoscience interpretation[END_REF].

In contrast to the forward problem, inversion moves in the opposite direction, that is, making inferences about the physical systems, mathematically represented through model parameters m (causes), from the observations d (consequences). This leads to problems that in most cases are ill-posed: their solution is neither unique, nor stable [START_REF] Carrera | Estimation of Aquifer Parameters Under Transient and Steady State Conditions: 2. Uniqueness, stability and solution algorithms[END_REF]. Thus, if any set of model parameters can be found that is able to explain the observations, then an infinite number of parameters-sets would exist and arbitrarily small errors in the measurement data may lead to indefinitely large errors in the solutions [START_REF] Kabanikhin | Definitions and examples of inverse and ill-posed problems[END_REF]. In other words, if the information gathered through the data is somehow "imperfect", then it is subjected to multiple plausible interpretations whose validity strongly depends on the initial assumptions [START_REF] Glendinning | Stability, Instability and Chaos[END_REF].

In geophysics, the ill-posed nature of the inverse problem can be solved by either reformulating it with a nearby well-posed one (deterministic approach), or explicitly incorporating additional knowledge or beliefs about the subsurface to supplement the noisy observations (probabilistic approach). In the first case, the solution to the inverse problem is given by the optimized set of model parameters that best explains the observations under certain constraints. The reformulation of the problem is carried out through regularization, which, in geophysical inversion generally penalizes the complexity of the solution [START_REF] Constable | Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data[END_REF], thus improving the conditioning of the problem [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF][START_REF] Tikhonov | Solution of ill-posed problems[END_REF]. Regularization can be interpreted as an implicit form of incorporating prior information in the formulation of the inversion problem. The choice of which feature to penalize (e.g., the distance from a reference model, the variability within model parameters, etc.), is in fact an implicit assumption of some knowledge or belief about the investigated subsurface before considering the observed measurements [START_REF] Calvetti | Inverse problems: From regularization to Bayesian inference[END_REF]. The deterministic formulation of the inverse problem involves assuming a locally linear relationship between data and model parameters, so that it is a suitable approach for linear, quasi-linear or "linearizable" problems. Moreover, the only source of uncertainty considered is the noise contained in the experimental measurements, while the errors embedded in the computed model predictions are assumed negligible. The noise in the observations can then be propagated to the estimates of model parameters, allowing variations around the optimal model [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF].

The probabilistic approach models the unknown solution as a random variable that is described through its probability distribution, obtained by combining all the known "states of information" (Tarantola & Valette, 1982). Thus, the probabilistic approach allows the description of multiple sources of uncertainties (i.e., noise in the data but also the uncertainties in the model predictions), which influence the degree of variability for the possible solutions to the problem [START_REF] Mosegaard | Resolution analysis of general inverse problems through inverse Monte Carlo sampling[END_REF]. The inversion problem can then be recast in the form of combining the knowledge gathered by the measurements with the prior information (or assumptions) on the investigated subsurface through Bayesian Inference. The term Bayesian refers to a branch of probability theory that allows the modeling of uncertainties in the outcomes of interest by combining different sources of knowledge and observational evidence (Bayes, 1763); while inference, or model evaluation, is the process of updating probabilities of outcomes based upon the relationships between the model and available evidence [START_REF] Jensen | An Introduction to Bayesian Networks[END_REF].

The choice of inversion routine implemented and especially the explicit or implicit assumptions made about the prior knowledge on the inferred parameters, strongly influence the inversion results, especially when the collected observations do not constrain the model parameters adequately [START_REF] Carrera | Estimation of Aquifer Parameters Under Transient and Steady State Conditions: 2. Uniqueness, stability and solution algorithms[END_REF]. It is thus often useful to consider multiple inversions with different underlying assumptions, for both verifying which features are well defined by the observations and to obtain an ensemble of rival models that can be used for decision support [START_REF] Ferre | Revising the Relationship Between Data, Models and Decision-Making[END_REF].

The probabilistic framework and Bayesian inference

The first implementations of a probabilistic framework for geophysical applications date to the late 1960s [START_REF] Westwater | Statistical information content used in indirect sensing[END_REF]Backus, 1970a,b,c), but Tarantola & Valette (1982) were the first to formalize this approach in general terms for geophysical inversion problems. Afterwards, there have been many contributions to both the theoretic and practical development within this framework (Duijndam, 1988a,b;[START_REF] Mosegaard | Monte Carlo sampling of solutions to inverse problems[END_REF][START_REF] Sambridge | Monte Carlo methods in geophysical inverse problems[END_REF][START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF][START_REF] Calvetti | Introduction to Bayesian scientific computing-Ten lectures on subjective computing[END_REF][START_REF] Hansen | Probabilistic integration of geoinformation[END_REF][START_REF] Vrugt | Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation[END_REF].

In probabilistic inversion, once the choice of model prediction strategy and parametrization is settled, it is possible to use a density function to provide a baseline description of the unknown model parameters. The prior distribution, ρ(m) describing all the a priori information available about the subsurface, is updated in the light of the observations, resulting in a posterior distribution, ρ(m|d), that summarizes all information included in the formulation of the inversion problem [START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF]. Bayes theorem can then be used to obtain such posterior distribution inasmuch it describes how to update the probabilities of hypotheses (m) given the evidence (d) (Bayes, 1763):

ρ(m|d) = ρ(m)ρ(d|m) ρ(d) , (1.3)
where ρ(d) is the evidence, or marginal probability of the data given a certain conceptual model, and ρ(d|m) is the conditional probability of the observations given the model parameters. The likelihood function describes the relative probability (i.e., degree of confidence in the occurrence of an event) for the proposed m to explain d and it is proportional to this conditional probability: L(m) ∝ ρ(d|m). Therefore, if the model parametrization is fixed and we measure the fit between observations and model predictions with the likelihood function, it is possible to extract all the relevant statistical information from the unnormalized distribution:

ρ(m|d) ∝ ρ(m)L(m). (1.4) This distribution is described analytically only for a few specific cases (e.g., linear problems based on the least-squares formalism). When this is not possible, the inversion problem solution is obtained by sampling model realizations distributed accordingly to ρ(m|d). In geophysics this is commonly addressed with the use of Markov chain Monte Carlo methods [START_REF] Chib | Understanding the metropolis-hastings algorithm[END_REF]Gilks et al., 1996;[START_REF] Gamerman | Markov Chain Monte Carlo, Stochastic Simulation for Bayesian Inference[END_REF] In summary, the main steps needed to formulate and solve the inversion problem with Bayesian inference (Fig. 1.1) are:

• the definition of the prior distribution,

• the formulation of the likelihood function,

• the sampling strategy used to explore the posterior distributions. Defining the prior distribution is of great importance in Bayesian inference, especially in the presence of few or noisy data or in high parameter-dimensions [START_REF] Scales | Prior information and uncertainty in inverse problems[END_REF]. When considering geophysical problems, the models describing the investigated subsurface include both geometrical (e.g., position and size of geological bodies, interfaces, layering structure, etc.) and physical parameters (e.g., electrical conductivity, mass density, porosity, etc.). The prior distribution describes all the knowledge available about these parameters, independently of the particular observations considered as data within the formulation of the inversion problem. It can be constructed based on existing measurements and our beliefs (conceptualizations, structural assumption) about the subsurface. Notably, the a priori knowledge concerning subsurface spatial variability (i.e., model parameter correlation) has a large influence on the posterior model realizations. When such variability is known (or assumed) it can be formulated in terms of either specific parametrization choices, or by considering priors that account for higher-order spatial statistics.

BAYES THEOREM

Geostatistics uses random functions to describe the spatial variability within earth models and can be implemented to define different degrees of spatial correlation within the prior distribution, with the hope of generating realizations that are geologically realistic [START_REF] Hansen | Probabilistic integration of geoinformation[END_REF]. The underlying random function, used as a model, can be based either on two-point statistics, as for multi-Gaussian fields (Figs. 1.2c-g) where the models are fully characterized by an expected value and a covariance function or matrix [START_REF] Dietrich | Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix[END_REF][START_REF] Laloy | Probabilistic inference of multi-Gaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction[END_REF], or on multiple-point statistics (Figs. 1.2h,i). In this latter case, the models have no parametric description and are instead inferred from a training image [START_REF] Guardiano | Multivariate Geostatistics: Beyond Bivariate Moments[END_REF][START_REF] Strebelle | Conditional simulation on complex geological structures using multiple point statistics[END_REF][START_REF] Mariethoz | Multiple-Point Geostatistics[END_REF][START_REF] Cressie | Image analysis with partially ordered Markov models[END_REF].

The use of higher-order statistics as priors are useful when the investigated subsurface is well-characterized. However, when the a priori knowledge is scarce or not reliable, it is desirable to assume as little as possible about the subsurface. When the model parameterization is done through spatial discretization, such assumptions are often implemented with the use of uncorrelated log-uniform priors over a range of reasonable values of the model parameters [START_REF] Scales | To Bayes or not to Bayes?[END_REF][START_REF] Sambridge | Monte Carlo methods in geophysical inverse problems[END_REF]. Nevertheless, the spatially uncorrelated uniform distribution provides maximum entropy, which in this case translates for a high level of "disorder/disorganization" of the parameter-values within a model (see Fig. 1.2b). Finally, this results in posterior model realizations that can be acceptable if interpreted in terms of their average, but if taken singularly, they are too variable to be geologically realistic. This is the case because the posterior realizations are consistent with the prior model assumed: within the sampling routine, whenever two sampled models have the ability to explain the data within a similar discrepancy, the prior contribution favors the model with spatial variability closer to the one maximizing the entropy. Moreover, [START_REF] Hansen | Inverse problems with non-trivial priors: efficient solution through sequential Gibbs sampling[END_REF] demonstrated how the use of such prior leads to hard inverse problems, practically unsolvable except when the number of model parameters is small. ors. The sequential re-simulation sampler provides an eff cient and non-biased approach to sample complex prior information, that is well designed for use with the generalized Metropolis algorithm.

Synthetic case study: Application of non-linear inverse problems with complex prior information

Figure 1 is a channel-based training image from which we generate an unconditional realization, Figure 2, using the single normal equation simulation algorithm, SNESIM [START_REF] Strebelle | Conditional simulation on complex geological structures using multiple point statistics[END_REF]). This will be our reference velocity model for a synthetic cross borehole inversion problem. Travel times are computed, traveling from the 20 sources located to the left in Figure 2, to 40 receivers located to the right. 3% Gaussian noise is added to the synthetic travel-time delay data, Figure 3, and used as observed data. We now consider solving this inverse problem, using the generalized Metropolis algorithm to generate samples of the a posteriori probability distribution, given the observed data, the assumed noise model, and an assumed prior model. We consider a number of prior models based on both 2-point and multiple-point based random models. All prior models are assumed to have the correct mean and variance, as obtained from the training image in Figure 1. The f rst 6 prior models are based on 2-point random models. The pure nugget model assumes no spatial correlation, and thus all model parameters are a-priori considered uncorrelated. The 'Gau( 1 2017) proposed an alternative description of priors within geophysical Bayesian inversion when little (or no) a priori knowledge is available. Working on parameterizations based on a spatially regular discretization, the basic idea is to change the focus by assuming a uniform prior on global summary statistics describing the model spatial variability, instead of on the uncorrelated model parameters. This choice implies that all the possible prior model realizations are not equally likely, but that is equally likely to sample a smooth model realization (i.e., with little spatial variability) as a completely "scattered" one (i.e., with large spatial variability). The outcome of this change in prior description are model realizations that better resemble subsurface property distributions. Moreover, this approach can be applied to any summary statistic chosen to describe the model parameters and it can handle different structure prior pdfs.

In a previous work, [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF] proposed an alternative solution to the formulation of prior pdfs in case of scarce or unreliable a priori knowledge. They introduced an explicit prior constraint on the model variability within an empirical Bayes (EB) inversion framework [START_REF] Casella | An introduction to empirical Bayes data analysis[END_REF]. This method is computationally more efficient because it doesn't explore the full space of possible model variability, focusing the search in the neighborhood of values assumed to be more reasonable for the investigated subsurface. The likelihood function describes how well the model prediction, g(m), explains the observed data, d, given a statistical description of the measurement noise and modeling inadequacies [START_REF] Cordua | Monte-Carlo full-waveform inversion of cross-hole GPR data using multiple-point geostatistical a priori information[END_REF]: the closer the model prediction is to the experimental data, the larger the value of its likelihood. The form taken by the likelihood function depends on what is hypothesized/known about the residual errors (here referred to as the combination of measurement noise and modeling inadequacies). The form given to this function, and especially the approach used to define it, has been the center of an animated debate especially within the hydrology community [START_REF] Mantovan | Hydrological forecasting uncertainty assessment: Incoherence of GLUE methodology[END_REF][START_REF] Beven | So just why would a modeler choose to be incoherent[END_REF][START_REF] Mcmillan | Rainfall-runoff model calibration using formal likelihood measures within a Markov chain Monte Carlo sampling scheme[END_REF].

In the formal approach, the form of the likelihood function is derived through an assumed statistical model for the residual errors. For example, the standard least squares approach for parameter estimation assumes mutually independent (uncorrelated) residuals identically distributed according to a Gaussian pdf with zero mean and a constant variance [START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF]. In this case the likelihood function assumes its "classical" form: .5) where σ i represents the standard deviation of the i-th residual. For computational reasons (namely, to avoid integer overflow) it is favorable to evaluate the natural logarithm of the likelihood, which is expressed as:

L(m) = N i =1 1 2πσ 2 i exp - 1 2 (g i (m) -d i ) 2 σ 2 i , ( 1
l (m) = - N 2 l og (2π) - 1 2 l og ( N i =1 σ 2 i ) - 1 2 φ, (1.6) where φ = N i =1 g i (m) -d i σ i 2
represents the misfit between observations and model predictions. Another simple alternative often adopted in geophysics is to assume that the data follow a two-sided exponential pdf, which results in the Laplacian likelihood function [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF]:

L(m) = N i =1 1 2σ 2 i exp -2 |g i (m) -d i | σ . (1.7)
While the Gaussian model is based on using the L 2 norm as a measure of the residual length, the Laplacian model is based on the L 1 norm, and it is a more "robust" method, since it better tolerates few outliers [START_REF] Claerbout | Robust modeling with erratic data[END_REF]).

An advantage of the formal approach is that the error assumptions are stated explicitly and their validity can be verified a posteriori [START_REF] Stedinger | Appraisal of the generalized likelihood uncertainty estimation (GLUE) method[END_REF]. Nevertheless, the problem (and main argument against it in the debate) is that it relies too strongly on assumptions on the residuals errors that do not hold in many applications, leading to bias in estimated parameter values and uncertainties [START_REF] Thyer | Critical evaluation of parameter consistency and predictive uncertainty in hydrological modeling: A case study using Bayesian total error analysis[END_REF]. On the other hand, informal likelihood functions have been proposed as a pragmatic approach to uncertainty estimation in the presence of complex residual error structures. An example is the generalized likelihood uncertainty estimation methodology (GLUE), where the likelihood function is specified a priori without linking it to an underlying error model [START_REF] Beven | The future of distributed models: model calibration and uncertainty prediction[END_REF]. The informal approach is attractive in situations where traditional error assumptions are violated because the modeler has flexibility when specifying the form of the likelihood function. However, since it makes no explicit reference to the underlying error model, its assumptions are implicit and cannot be verified a posteriori.

Finally, an alternative to integrating information through the likelihood function (either within the formal or not formal approach) is the approximate Bayesian computation (ABC). In this framework, model realizations sampled from the prior pdf are kept as populations of the posterior pdf whenever the chosen summary statistics of the model predictions are within predefined distances from those of the observations [START_REF] Marjoram | Markov chain Monte Carlo without likelihoods[END_REF][START_REF] Turner | A tutorial on approximate Bayesian computation[END_REF].

Extracting information from the posterior distribution is generally addressed through Monte Carlo simulations and, in particular Markov chain Monte Carlo (MCMC; [START_REF] Brooks | Handbook of Markov chain Monte Carlo[END_REF]. This choice has been naturally developed in geophysics to face two related issues: the presence of multiple secondary maxima in the posterior distribution (representing possible solutions), and the non-linearity and multi-dimensional character of the inverse problem [START_REF] Mosegaard | Monte Carlo sampling of solutions to inverse problems[END_REF]. In fact, in the case of highly nonlinear and high-dimensional problems the posterior pdf has a complex multimodal shape and the definition of "central estimators" (e.g., mean or median)

and "estimators of dispersion" (i.e., variance and covariance matrix) have little meaning [START_REF] Mosegaard | Probabilistic Approach to Inverse Problems[END_REF]. Therefore it is necessary to implement global sampling methods to gather information on the complete shape of the posterior distribution. In this framework, the simplest approach is an exhaustive search, where all the models within a chosen subspace are visited, but in case of high-dimensional and non-linear problems, also limiting the number of misfit calculations is necessary. Moreover, a large-dimensional model space tends to be 'empty", therefore, the probability of sampling the target distribution rapidly tends to zero with the increase in dimensionality. MCMC methods are a common choice in geophysics because they directly simulate the posterior pdf, filtering proposed values from the prior distribution of model parameters in order to obtain a sample of values, referred as a chain, distributed according to the posterior and from these calculate Bayesian estimates [START_REF] Koren | Monte Carlo estimation and resolution analysis of seismic background velocities[END_REF][START_REF] Gouveia | Bayesian seismic waveform inversion: Parameter estimation and uncertainty analysis[END_REF].

Historically, the starting point of MCMC lies in the work of [START_REF] Metropolis | The Monte Carlo method[END_REF] and others in the 50s, together with advancements of [START_REF] Hastings | Monte Carlo sampling method using Markov chains and their applications[END_REF] and [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF]. The idea behind MCMC is to sample a target distribution through the construction of a computationally realizable random Markov process (i.e., memoryless process). A special property of the chain is that the probability of sampling a certain future state of the process (i.e., a new sample within the chain) depends only on the present state and not on the sequence of events that preceded it, which is the socalled Markov property ("the future is independent of the past given the present"; [START_REF] Markov | Theory of Algorithms[END_REF]. This means that the sampling algorithm is completely defined by the transition probabilities between states of the chain. There exists many MCMC sampling algorithms, two of which are often implemented in geophysics are Metropolis-Hastings (Metropolis et al., 1953;[START_REF] Hastings | Monte Carlo sampling method using Markov chains and their applications[END_REF] and Gibbs [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF] samplers.

Given a target pdf that is hard to sample from, P (m), and a transition kernel, T (m j |m i ), describing the probability to move in the model space from the state i to j, Metropolis-Hastings generates a sequence of iterates {m 1 , m 2 , . . . , m M } which approximates a unique stationary pdf:

lim i t er →∞ π(m) = P (m). (1.8)
In the case of Bayesian inference, the target distribution is the posterior:

P (m) = ρ(m|d) ∝ L(m)ρ(m).
The existence and uniqueness of the sampled pdf is guaranteed because Metropolis-Hastings satisfies [START_REF] Chib | Understanding the metropolis-hastings algorithm[END_REF]:

1. Detailed balance (or reversibility), which is a sufficient condition for a random walk to asymptotically reach a stationary pdf and it requires each transition to be reversible:

T (m pr op |m cur r )L(m cur r )ρ(m cur r ) = T (m cur r |m pr op )L(m pr op )ρ(m pr op ), (1.9) where the subscript curr refers to the current position of the chain and prop to the proposed one. 2. Ergodicity of the Markov process, which requires that every state must be aperiodic (the system does not return to the same state at fixed intervals), positive recurrent (the expected number of steps for returning to the same state is finite) and irreducible (each state is accessible in a finite number of steps).

One challenge is then to define a proper transition kernel, sure that the sequence of samples drawn will converge to the target distribution. When the proposal probability q, which describes how the chain moves in the model space, satisfies the reversibility condition (eq.1.9), then the transitional kernel coincides with the proposal distribution: T (m pr op |m cur r ) ≡ q(m cur r -→ m pr op ). In practice, this is often not the case, therefore, the Metropolis-Hastings introduces the acceptance probability 0 < α(m cur r , m pr op ) < 1 that the move happens in order to achieve detailed balance:

α M H = mi n 1, L(m pr op )ρ(m pr op )q(m pr op → m cur r )
L(m cur r )ρ(m cur r )q(m cur r → m pr op ) .

(1.10)

The proposal pdf (q) generates model perturbations at each proposal step. In many applications, such a distribution is chosen to be symmetric, which means that there is the same probability to move from one state to another and in the other direction: q(m pr op → m cur r ) = q(m cur r → m pr op ). In this case the Metropolis-Hastings acceptance ratio simplifies to the Metropolis ratio: .11) Gibbs sampling is a special case of the Metropolis-Hastings algorithm, where each iteration consists of M sub-steps, one for each model parameter: the mth sub-step perturbs only the m-th parameter, and it has its own transitional kernel.

α M = mi n 1, L(m pr op )ρ(m pr op ) L(m cur r )ρ(m cur r ) . ( 1
MCMC simulations are often started at random points in the model space that are drawn from the prior pdf. These points are generally far from high density regions of the posterior distribution, therefore in the early stages of the simulations the sampled models are unlikely to produce representative samples from the target pdf. This initial part of the chain is referred as the "burn-in" period and the remaining part is the stationary part of the chain, or the part where the chain has "converged in distribution " [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF]. Operationally, effective convergence of Markov chain simulations is reached once the inferred quantities of interest do not depend on the starting point of the chain [START_REF] Brooks | General Methods for Monitoring Convergence of Iterative Simulations[END_REF]. [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF] proposed an analysis based on the comparison of inferences from different chains. Their analysis results in a potential scale reduction factor which gives a measure of how much the posterior distribution has been explored by the Markov chains: when it is close to 1, it means that each of the x chains is close to the target distribution.

Designing reliable and computationally efficient sampling techniques is fundamental, especially for high dimensional problems (such as inverse problems of spatially distributed parameters), where running the forward solver can be time-consuming.

The trade-off between exploration and exploitation is of central importance to evaluate different algorithms and to design new sampling methods. Exploration refers to the ability of the chain to gather more information, while exploitation refers somehow to the opposite: the ability to choose where to sample using the information from previous samples [START_REF] Vrugt | Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation[END_REF]. Therefore, a more explorative MCMC sampler would use a wide proposal pdf, which might lead to a high rejection rate for the candidate points and consequentially to a slow convergence to the target distribution. On the other hand, a more exploitative MCMC sampler would use a narrow proposal pdf, which would generally result in accepting almost all the proposal candidate points, but with a short distance moved, so that the chain can converge to the target pdf only after a large number of updates [START_REF] Vrugt | Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation[END_REF]. Moreover, the choice of the sampling algorithm should also facilitate the description of the subsurface accounting for the target phenomenon/physical property of the investigation. In this field, the work of [START_REF] Iglesias | Well-posed Bayesian geometric inverse problems arising in subsurface flow[END_REF] investigates inverse problems for determining subsurface permeability from hydraulic head measurements in case of sharp discontinuities due to geometrical subsurface features (as for instance channel structures, layers, faults, etc..). In such a problem, both the geometry and the physical characteristics of the explored area are to be inferred.

The authors thus introduce a Metropolis within Gibbs algorithm, which separates the effect of parameters describing the geometry (i.e., location and shape of regions where discontinuities in subsurface permeability arise due to the presence of different geologic facies) from those describing spatial variability of the target physical property. This results in a more precise description of the subsurface in case of sharp boundaries.

A general issue with geophysical observations, independent of the assumptions made for their interpretation, is that their information content is limited. Even if geophysical surveys can result in vast datasets, only a finite number of model parameters can be independently inferred from them [START_REF] Backus | Uniqueness in the inversion of inaccurate gross earth data[END_REF]. There are multiple reasons for this limited resolution of subsurface properties: from restrictions on the possible acquisition geometries, to the noise-contamination of the observations and the undergoing physics related to the measurement produced (as for instance diffusion or wave propagation). Beyond awareness of these limitations, complementary datasets about the same investigated subsurface can be incorporated and interpreted together, thus reducing the inherent ambiguity affecting geophysical observations (Moorkamp et al., 2016).

Joint inversion provides a formal approach to integrate multiple datasets with the aim of better constraining the model results [START_REF] Vozoff | Joint Inversion of Geophysical Data[END_REF][START_REF] Moorkamp | A framework for 3-D joint inversion of MT, gravity and seismic refraction data[END_REF][START_REF] Hellman | Structurally coupled inversion of ERT and refraction seismic data combined with cluster-based model integration[END_REF]. Combined interpretation of data from several geophysical methods can be addressed in three main ways [START_REF] Linde | Joint Inversion in Hydrogeophysics and Near Surface Geophysics[END_REF]:

i) Joint interpretation of the different datasets refers to when observations are separately inverted for, resulting in different models that are afterwards interpreted together. This approach is useful especially when the inversion results of one geophysical data type are used as prior information for the inversion of others [START_REF] Saunders | Constrained resistivity inversion using seismic data[END_REF][START_REF] Doetsch | Zonation for 3D aquifer characterization based on joint inversions of multi-method crosshole geophysical data[END_REF]. Nevertheless, it does not reduce the inherent ambiguity on the parameters estimation and it is generally not considered as joint inversion since the data misfit of the different observations with the model predictions are never considered simultaneously.

ii) Joint inversion of separate data sets for a common parameter. This is often accomplished through petrophysical relationships, used to transform the primary property fields into geophysical properties that the observations are sensitive to (Hetrich & Yaramanci, 2002;Jardaniet al., 2010). The main attractive feature of joint inversion by petrophysical coupling is that it allows formulating the inverse problem in terms of the target properties of primary interest (e.g., porosity, permeability, lithology). Nevertheless, this approach is valuable only when global search inversion methods are used since local inverse formulations in fact are very sensitive to errors in the petrophysical model and it is likely that model artifacts will be introduced to compensate for these errors. Moreover, it introduces additional sources of uncertainty in the inversion problem through the petrophysical relation themselves.

iii) Coupled inversion refers to when separate datasets constrain each other, generally through structural coupling which seeks multiple distributed models that share common interfaces or have a similar model gradients [START_REF] Gallardo | Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints[END_REF].

The assumption in this case is that the proposed models, besides explaining the different collected observations, must also have similar structure as they are describing the same subsurface area. This strategy is also the most robust coupling for the observations [START_REF] Linde | Joint Inversion in Hydrogeophysics and Near Surface Geophysics[END_REF].

Bayesian inference grants a natural framework to incorporate information obtained from different sources [START_REF] Reid | Bayesian joint inversions for the exploration of earth resources[END_REF], nevertheless, for coupling the models explaining the different observations it is necessary to include a constraint between them. If we consider two sets of observations d A and d B respectively sensitive to the physical properties parameterized with m A and m B , the posterior distribution on m A and m B can be written as:

ρ(m A , m B |d A , d B ) ∝ ρ(m A )ρ(m B |m A )L(m A )L(m B ).
(1.12)

Here ρ(m B |m A ) contains the information about the constraints between models. The Metropolis within Gibbs algorithm previously introduced allows spontaneously such constraint is in a structural form. This sampling method accounts for the subsurface geometry, parameterizing and inferring it within the inversion routine, that updates alternately the geometric and physical parameters within the Gibbs framework. The geometric parameter can then be used as structural constraints between the different models when jointly inverting multiple datasets. Figure 1.3 schematically represents a flow chart of joint inversion of two (or more) geophysical datasets that can be applied to any type of MCMC samplers (besides the Metropolis proposed in the original algorithm presented by [START_REF] Iglesias | Well-posed Bayesian geometric inverse problems arising in subsurface flow[END_REF]. The choice of the sampler (e.g., Metropolis-Hastings, Metropolis, etc...) would simply modify the acceptance ratio α, that, for the sake of simplicity, is herein noted only as proportional to the likelihood ratio.

1.2

Geophysical methods

Among the various available geophysical methods, the inversion algorithms proposed in this thesis were tested on observations from crosshole ground penetrating radar (GPR), electrical resistivity tomography (ERT) and seismic refraction considering both field datasets and synthetic test cases.

Ground Penetrating Radar

Ground penetrating radar uses electromagnetic (EM) fields to probe the subsurface.

Compared with other technologies, GPR is a relatively new geophysical exploration tool: one of its earliest successful applications was the measure of ice thickness on polar ice sheets in 1960s [START_REF] Waite | Gross errors in height indication from pulsed radar altimeters operating over thick ice or snow[END_REF]. Since then, there have been rapid developments in hardware, measurement and analysis techniques, and the method has been extensively used in many applications, such as archeology, civil engineering, forensics, geology and utility detection [START_REF] Daniels | Ground penetrating radar 2nd edition (Radar, Sonar &amp[END_REF]. For ground penetrating radar data, the tomographic inversion is based on electromagnetic waves transmitted through the medium. The traveltime and amplitude of propagating EM wave traveling through the Earth's subsurface are controlled mainly by the electric properties of the medium [START_REF] Born | Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light[END_REF]. In terms of EM wave transmission and reflection, the most important property is , which describes polarization effects resulting from bounded charges in the medium.

When considering GPR tomography, a typical setup for a crosshole experiment consists of a transmitting dipole-type antenna located in a borehole and a corresponding receiving antenna in a neighboring one [START_REF] Davis | Ground-penetrating radar for high-resolution mapping of soil and rock stratigraphy[END_REF]. Considering different depths in the two boreholes, it is possible to gather multi-offset traveltimes, representing the times it takes for the first arrivals of high-frequency electromagnetic pulses. The measurements are then taken for all possible combinations of transmitter and receiver positions to ensure that the medium between the boreholes is sampled by a large number of rays to enable a sufficient angular coverage. The highest resolution and data cumulative sensitivity is found in the central part of the tomogram (Day- [START_REF] Day-Lewis | Applying petrophysical models to radar travel time and electrical resistivity tomograms: Resolution-dependent limitations[END_REF].

The simplest ray-based modeling possible of first-arrival times is based on the assumption that velocity variations along the ray-path are relatively small, so that the EM wave energy propagates dominantly along straight rays [START_REF] Holliger | Ray-based amplitude tomography for crosshole georadar data: a numerical assessment[END_REF]. Within this assumption, reconstructing tomograms requires solving a linear system of equations [START_REF] Peterson | Applications of algebraic reconstruction techniques to crosshole seismic data[END_REF]. If the region between a transmitter (T x ) and a receiver (R x ) is discretized into pixels (Fig. 1.4), and the slowness p (inverse of the velocity v) of the pixel j is assumed to be constant over the area of a pixel, the relation between propagation velocity and total traveltime for a particular ray-path is given by [START_REF] Hanafy | Ground-penetrating radar tomography for soilmoisture heterogeneity[END_REF]:

t i = R i da v(x, y) = R i p(x, y)da, i = 1, 2, . . . k (1.13)
where t i is the total traveltime along k different ray-paths, v(x,y) is the EM wave propagation velocity along each ray-path, da is the differential length along the ray, and the integration is along the particular ray-path R i . The discrete form of Eq. 1.13 is easily implemented when the ray-paths are assumed to be straight lines through the finite difference (FD) approximation: (1.14) where ∆a i j is the length of the ray i that crosses the pixel j, and m is the total number of pixels intersected by ray i. The assumption of straight ray-paths is seldom valid in real problems and in such case the eikonal equation that considers non-linear effects could be considered. Moreover, ray methods are based on high-frequency approximations and only account for a small fraction of the information contained in the radar traces and are restricted to resolving relatively large-scale features. To address this issue, Ernst et al. (2007a,b) introduce a full-waveform inversion scheme that is based on a finite-difference time-domain solution of Maxwell's equations [START_REF] Klotzsche | 3D characterization of high-permeability zones in a gravel aquifer using 2D crosshole GPR full-waveform inversion and waveguide detection[END_REF].

t i = m j =1 ∆a i j p j , i = 1, 2, . . . m,

Electrical Resistivity Tomography

Electrical resistivity tomography collects information about subsurface electrical resistivity by injecting direct (DC) electric currents into the ground and measuring electric voltages at different locations. The data acquisition is performed with quadripole geometry, where the current flowing between two electrodes is causing the voltage measured between a second pair of electrodes through a high-impedance voltmeter. This geometry is necessary to avoid an impact of contact resistances (depending on ground moisture, contact area, etc.). The high-impedance voltmeter draws virtually no current so that the voltage drop between the electrodes is negligible. Thus, the contact resistance still limits the current flow, but does not affect the resistivity calculations.

A geometric factor is then needed to convert the readings obtained into apparent resistivity. The observations are collected through measurements of multiple electrode configurations (many tenths to hundreds) with different spacing and central points, yielding a 2D dataset along a profile line which is usually plotted as pseudo section (see Fig. The basics of ERT has its roots in the resistivity sounding methods, that were first developed by Conrad Schlumberger in 1912 [START_REF] Kunetz | Principles of Direct Current Resistivity Prospecting[END_REF]. Nevertheless, it is only with the rapid development of multiple-electrode equipment in the 1990s [START_REF] Dahlin | On the automation of 2D resistivity surveying for engineering and environmental applications[END_REF] that ERT studies multiply in different applications, as in engineering problems [START_REF] Kunetz | Principles of Direct Current Resistivity Prospecting[END_REF], soil science (Samouëlian et al., 2005) or monitoring hydrogeological processes [START_REF] Robinson | Advancing process-based watershed hydrological research using near-surface geophysics: a vision for, and review of, electrical and magnetic geophysical methods[END_REF]. ERT is largely used for the imaging of the shallow subsurface (tens to hundreds meters depth), for which it provides comparatively high resolution images.

In the theoretical case of a homogeneous earth, the expression for the potential at a surface location, caused by a point source located at the surface, is: .15) where I is the intensity of the current injected at the point source, l is the distance between the source and the location where the potential is considered, and r is the resistivity of the homogeneous earth. In reality the earth is not homogeneous, therefore the quantity r is replaced with the apparent resistivity (r a ). This is defined as the resistivity of a homogeneous, isotropic ground that would give the same voltage-current relationship as the one measured . Usually what it is measured in ERT surveys is the electrical impedance V p I , which is then transformed into apparent resistivity by means of the geometric factor k. In eq.1.15, this factor is simply k

V p = r I 2πl , ( 1
= 1 2πl
, but more generally is a function of both the electrode layout and surface geometry. Finally, the apparent resistivity pseudo-sections are the observations from which it is possible to infer a subsurface resistivity model by inversion.

Electrical resistivity (the inverse of conductivity) quantifies how strongly a material opposes the flow of electric current. In most rocks and soils, electrical current is carried by movements of ions in the pore water, with the actual mineral matrix practically being an isolator. Thus, important factors affecting subsurface resistivity are water saturation, salinity, temperature, porosity and the connectivity of the water phase [START_REF] Lesmes | Relationships between the Electrical and Hydrogeological Properties of Rocks and Soils[END_REF]. Moreover, the clay content strongly influences the electrical resistivity through its contribution to surface conductivity occurring at the mineralwater interface related to the electrical double layer [START_REF] Revil | The Self-Potential Method: Theory and Applications in Environmental Geosciences[END_REF]. Among the many petrophysical relationships available, a quantitative formulation linking the different contributions to the global conductivity of the subsurface is [START_REF] Linde | Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data[END_REF]: .16) where σ eff is the effective conductivity of the matrix-water system, σ s is the surface conductivity and σ w is the pore water conductivity. F = φ -m is the electrical formation factor which depends on the porosity, φ, of the matrix and on the cementation exponent m, S w is the water saturation and n is the saturation exponent, or Archie's second exponent, which is inversely related to the connectivity and tortuosity of the water phase.

σ eff = 1 F S n w σ w + (F -1)σ s , ( 1
An inherent problem related to resistivity measurements in general (sounding as well as tomography) is a phenomenon described as equivalence or suppression [START_REF] Koefoed | Geosounding Principles I, Resistivity Sounding Measurements[END_REF]. These terms refer to the fact that strongly different resistivity distributions in the subsurface may lead to apparent resistivity curves (in case of sounding) or fields (in case of ERT) which differ so slightly they cannot be distinguished within the accuracy of the measurements. Moreover, ERT observations have very limited sensitivity to the depth and thickness of deep subsurface layers, yet to simulate the model response of such measurements it is necessary to consider models that extend outwards from the survey area and to greater depths in order to avoid boundary effects on the numerically computed electrical field [START_REF] Oldenburg | Estimating depth of investigation in dc resistivity and IP surveys[END_REF].

Refraction seismics

Refraction seismics is based on the analysis of first arrival traveltimes of critically refracted seismic waves. A seismic wave is energy transmitted by vibration of rock particles, and the aim in exploration seismology is the reconstruction of subsurface wave velocity fields. The beginning of seismic methods can be dated to the experiment of Mallet (1847) , who was the first to use an artificial source for generating seismic waves.

Refraction seismics has undertaken significant developments during the First World War when the research focused on the location of heavy artillery from the study of the waves generated by the recoil of firing guns. This work was further developed by Ludger Mintrop, who obtained the first patent for a portable seismograph in 1919 [START_REF] Keppner | Ludger Mintropp[END_REF]. Afterwards, the advance in technology provided during Second World War and especially the beginning of the computing revolution in the late 70s early 80s brought major development for this methodology. Even though seismic methods are mainly used in hydrocarbon exploration [START_REF] Berg | Seismic stratigraphy II: An integrated approach to hydrocarbon exploration[END_REF], they have a large number of applications, as for depth-to-bedrock detection [START_REF] Saas | Bedrock detection and talus thickness assessment in the European Alps using geophysical methods[END_REF], characterization of fracture system in the subsurface [START_REF] Liu | Seismic Fracture Characterization. Concepts and Practical Applications[END_REF] and they are also widely used in detecting the water table or saturated aquifer thickness [START_REF] Nur | The effect of saturation on velocity in low porosity rocks[END_REF].

A seismic wave in solid material can travel in two ways: either through the oscillation of molecules backward and forward in the direction of the energy transport (i.e., pressure wave -P wave-) or perpendicular to this direction (i.e., shear wave -S wave-). Seismic refraction focuses principally on P waves and it is based on the principle that, when seismic waves impinge on a boundary across which there is a contrast in velocity, their direction of propagation changes. The amount of this change in the direction of the wave propagation depends on the contrast in seismic velocity across the boundaries accordingly to Snell Law [START_REF] Blackstock | Fundamental of physical acoustics[END_REF]:

sin Θ i sin Θ r f = V 1 V 2 , V 2 > V 1 , (1.17)
where Θ i is the angle of incidence, Θ r f is the angle of refraction and V 1 and V 2 are respectively the velocity of the upper and lower layer. When sin

Θ i = V 1 V 2
, then the refracted wave travels parallel to the interface and some of the acoustic energy will return to the surface as head waves leaving the interface at the same angle of incidence: the critical angle Θ i c (see Fig.1.6). This critical refraction is the base of refraction seismics which collects the first arrival times of such waves in order to infer a 1D/2D/3D P-wave subsurface velocity model.

𝛩 ic 𝛩 ic V 2 > V 1 V 1 geophones source

Critical refracted wave

Direct wave

Head wave Figure 1.6: Seismic critical refraction ray-path scheme. In figure Θ i c is the critical angle and V 2 > V 1 are the P-waves propagation velocities of the two layers.

When the composition of the explored subsurface is complex, as in the case of fractured rocks or partially saturated soils, seismic properties can be quantitatively interpreted through an equivalent medium representation. Here, the multiphase (minerals, water, air) properties and their connectivities are replaced locally (at the scale of a representative elementary volume) by an upscaled homogeneous medium with the same macroscopic properties [START_REF] Liu | Seismic Fracture Characterization. Concepts and Practical Applications[END_REF]. Nevertheless, within quite broad limits, the velocity of a mixture of different materials can be obtained by averaging the transit times (reciprocals of velocities) through the pure constituents, weighted accordingly to the relative amounts present. This principle can be used also when one of the constituents is a liquid [START_REF] Wyllie | Elastic wave velocities in heterogeneous and porous media[END_REF][START_REF] Wyllie | AN experimental investigation of factors affecting elastic wave velocities in porous media[END_REF]):

1 .18) where V eff , V f and V m are the P-wave velocities of the saturated rock, the pore fluids and the rock matrix (mineral matrix), respectively. When considering dry rocks, the pore space is saturated with air rather than water and equation 4.5 can not be applied. Nevertheless, dry materials generally have low P-wave velocities. Moreover, poorly consolidated water-saturated materials have velocities slightly higher than that of water, such that the water table is generally a prominent seismic interface.

V eff = φ V f + 1 -φ V m , ( 1
First-arrival refraction methods use only a small portion of the information contained in the seismic traces and strongly depend upon there being a general increase in velocity with depth. In case of velocity inversion (i.e., the lower medium has a lower P-wave velocity than the overlaying one), the refracted wave will bend towards the normal. This gives rise to the so-called "hidden layer" phenomenon [START_REF] Banerjee | Hidden layer problem in seismic refraction work[END_REF]. There are different scenarios where refraction seismics is unable to detect subsurface layers.

The presence of such layers can sometimes be recognized from secondary arrivals, but this is only occasionally possible, especially because refracted waves are strongly attenuated in thin layers. Also, if the velocity decreases at an interface, critical refraction cannot occur and no refracted energy returns to the surface.

Modeling tools

Within probabilistic inversion, the efficiency of the forward code is crucial and can be sometimes favored over precision. When considering high-dimensional inversion problems, a significant exploration of the posterior distribution can require up to a few millions of iterations. Thus, an efficient forward code is necessary in order to run the MCMC chain(s) in a reasonable time.

GPR

Concerning the GPR forward model, we considered electromagnetic energy to propagate along ray paths. Therefore the forward model was evaluated through the finite difference approximation described in equation 1.14, resulting in an efficient forward operator defined on regular grids [START_REF] Podvin | Finite difference computation of travel times in very contrasted velocity models: a massively parallel approach and its associated tools[END_REF].

ERT

For the computation of the ERT forward response the modeling is demanding because it requires numerical computation of the electric field. For this purpose, both the finite difference (FD -Mufti, 1976-) and finite element (FE -Coggon, 1971-) methods have been developed. A main limitation of the first method is the restriction to orthogonal grids, which limits the ability of FD to reproduce non-orthogonal geometries, whereas FE approaches are not subject to such drawbacks.

FE methods are based on Hamilton's principle of minimization of an energy function (in this specific case, EM energy). Subdividing the modeling domain into disjunct elements, allows to solve a weak formulation of this minimization problem within each element. This finally results in a system of linear equations that is solved by splitting the computation in two: a background potential (primary field) and a secondary one, caused by the conductivity deviations from the background homogeneous model (i.e., singularity removal; [START_REF] Lowry | Singularity removal: a refinement of resistivity modeling techniques[END_REF].

For computing the ERT model predictions, we used the forward operator from the Boundless Electrical Resistivity Tomography (BERT) python library [START_REF] Rücker | Three-dimensional modelling and inversion of dc resistivity data incorporating topography-I[END_REF][START_REF] Günther | Three-dimensional modelling and inversion of dc resistivity data incorporating topography-II. Inversion[END_REF]. BERT forward operator allows efficient numerical computation of 22 the electrical potential with finite-element methods in 1,2 and 3 dimensions. One of its main innovations is the incorporation of tetrahedral/triangular unstructured meshes for model discretization (Si, 2015;[START_REF] Shewchuk | Thetraedral mesh generation by delaunay refinement[END_REF]. This allows a flexible description of arbitrary model geometries and efficient local refinement, therefore, saving computing resources. In fact, refinement of regular grids leads to an excessive increase in the number of nodes, which consequentially expands the numerical effort. Unstructured meshes on the other hand facilitate the local refinement of the grid within specific regions without requiring a significant increase in the number of nodes. Thus, it is possible to choose a very fine mesh in zones of varying potential gradients (e.g., close to the electrodes or at strong conductivity contrasts) and to increase the cell size toward the boundaries of the domain, where the electrical field is smoother. Finally, within the BERT library a numerically calculated geometric factor k is used to consider complex topography.

Refraction Seismics

The first arrival traveltimes for seismic refraction were simulated with the physics refraction class of pyGIMLI (Geophysical Inversion and Modelling Library in Python; [START_REF] Rücker | pyGIMLi: An open-source library for modelling and inversion in geophysics[END_REF]. Specifically, the forward operator is based on the calculation of the shortest traveltime path from source to receivers through a network that represents the earth to approximate seismic ray-paths [START_REF] Moser | Shortest path calculation of seismic rays[END_REF][START_REF] Dijkstra | A note on two problems in connection with graphs[END_REF]. This approximation is based on Fermat's principle [START_REF] Schuster | An Introduction to the Theory of Optics[END_REF], which states that seismic rays follow the path that gives the shortest traveltimes between points. This principle is then implemented through the analogy between seismic ray-path and the shortest path in a network. Within this method, the investigated subsurface is in fact represented by a network consisting of points (nodes), each of which is connected with a restricted number of other points within its neighborhood. A connection within two nodes is then weighted as the traveltime of a seismic wave between them and, by virtue of the reciprocity principle, it is invariant to the directions. The main advantage of the shortest path method is its simplicity, capacity for simultaneous calculation of the first arrival times and the associated ray paths of each mode without missing any "target" (receiver) in a complex geological structure [START_REF] Zhou | Shortest path' ray tracing for most general 2D/3D anisotropic media[END_REF]. On the other hand, by forcing the seismic ray-paths to follow the connection of a network, errors may rise in the geometry and therefore traveltime along the rays due to space and angle discretization [START_REF] Sethian | Fast Marching Methods[END_REF].

pyGIMLI is an open source library, firstly written in C++ (GIMLI) and recently translated in Python, that provides tools for modeling and inversion of various geophysical but also hydrogeological methods (pyBERT, the Python implementation of BERT, is built within pyGIMLI). The modeling component of this software supplies discretization tools and the numerical basics for both FD and FE solvers in 1D, 2D and 3D on structured and unstructured meshes. Moreover, allowing the computation of different geophysical model predictions, pyGIMLI gives a facilitating framework for the model coupling necessary for the joint inversion of multiple geophysical observations.

Objectives

The aim of this thesis was to address the choice of priors within Bayesian inversion of geophysical datasets when little is known about the investigated subsurface variability, and to infer unknown interfaces, together with their uncertainties, in the presence of significant subsurface heterogeneity. Specifically, we firstly focused on the effects of using log-uniform uncorrelated prior distributions on the posterior model realizations when parameterizing the subsurface through spatial discretization (Chapter 2). The twofold objective was to raise awareness of the consequences carried by this choice of prior pdf (commonly implemented in the case of scarce a priori knowledge about the investigated subsurface) and to propose an alternative description. Particularly, describing the prior in terms of model spatial variability allows a characterization that is easily understood in geoscience and especially for researchers used to classical regularized deterministic inversion methods. In this sense, re-defining a uniform prior pdf in terms of spatial model variability implies the same a priori probability of sampling smooth models as highly variable ones, finally leaving the observations to carry the information about subsurface structure. More generally, this approach can help to understand which are the features of the subsurface that are effectively constrained by the observations.

The second topic concerns the way we describe and sample the inversion model parameters, especially when the focus of the investigation are on geometrical features of the subsurface (specifically the detection of interfaces; Chapter 3). The central point in this case, besides inferring the depth and topography of interfaces from geophysical observations was the quantification of uncertainty. The probabilistic framework is a powerful tool for uncertainty quantification and its practical implementation in case of geometric parameters was introduced by [START_REF] Iglesias | Well-posed Bayesian geometric inverse problems arising in subsurface flow[END_REF]. With the aim of sampling target interfaces in the presence of subsurface heterogeneities, we adapted their Metropolis-within-Gibbs algorithm to an inversion framework that takes into account also weak prior information.

The third topic we addressed is the inherent ambiguity of geophysical inversion results in describing all the model parameters due to the limited information content based on a single set of observations (Chapter 4). It is well known that joint inversion of multiple datasets can help reducing such ambiguity, therefore, we extended the methodology proposed to probabilistically infer subsurface interfaces to jointly invert multiple geophysical datasets. This was a straightforward development of the algorithm since model coupling between the different physical properties can be achieved by considering a common interface.

Outline

The next three chapters present work that has been published, submitted or is in preparation for submission in peer-reviewed journals. Namely:

• In Chapter 2 we present a novel approach for MCMC inversion of geophysical data.

Here, the prior distributions are re-defined in terms of model spatial variability (S(m)) instead of marginal model parameters. This re-definition, requires an algorithm that accounts for the asymmetry of the proposal distribution for S(m) within the Metropolis-Hastings acceptance ratio (eq.1.10). In the absence of an analytical expression for such proposals, one of the main constriction of this work has been the numerical estimation of pdf ratios from small samples of their population with empirical kernel density estimation or method of moments. The approach was evaluated with synthetic and field-based ground penetrating radar observations, showing in both cases that the posterior model realizations and statistics obtained are significantly more satisfactory than those obtained when assuming uncorrelated model parameters or using explicit penalties on model structure within an empirical Bayes framework. The same concepts and methodology are reformulated for a wider readership in Appendix A.

• In Chapter 3 we describe the implementation of an empirical-Bayes-within-Gibbs algorithm which explicitly parameterizes and infers both interface geometry and spatial heterogeneity of physical properties. We assumed that reliable prior knowledge was absent and constrained both the interface and the physical properties to favor smooth spatial transitions in an empirical Bayes framework.

Here, we couldn't apply the methodology presented in Chapter 2, because of the considerable wider ranges for the model parameters values in this case, which makes the task of numerically approximating the proposal pdfs for model variability computationally too demanding. We evaluated the methodology on synthetic and field surface-based ERT datasets, with the aim of inferring regolith-bedrock interfaces. The proposed algorithm improves the accuracy in detecting the interface location when compared to maximum gradient methods applied to deterministic inversion results. Moreover, the introduction of the interface as a parameter to infer for gives the possibility to quantify its uncertainty and leads to a dramatic improvement in the estimations of bedrock properties when compared with smoothness-constrained deterministic inversions and MCMC-based inversions without an explicit interface.

• The field case dataset study in Chapter 3 was first presented by St. Clair et al. (2015), where ERT and seismic refraction surveys were carried out in order to investigate the influence of surface topography on bedrock geometry and properties.

In Chapter 4 we extended the empirical-Bayes-within-Gibbs method to probabilistically and jointly invert both datasets. The model coupling was achieved by considering a common interface geometry, allowing a partial resolution of the inherent ambiguity in the interpretation of one single set of geophysical observations. Through a synthetic test case it is shown that the joint inversion is more accurate in sampling the target model compared to inversion results obtained from single datasets. Moreover, the introduced methodology is applicable to other method combinations and it could easily be extended to three or more datasets.

Finally, Chapter 5 concludes with a summary of this thesis work, some remarks on the limitations of the different methodologies presented and possible outlooks.

Chapter 2

On structure-based priors in Bayesian geophysical inversion

Giulia de Pasquale and Niklas Linde.

Published 1 in Geophysical Journal International and herein slightly adapted for this thesis.

1 de Pasquale, G. and N. Linde (2017). On structure-based priors in Bayesian geophysical inversion.

Geophysical Journal International, 208(3) 1342-1358, doi:10.1093/gji/ggw458

Abstract

Bayesian methods are extensively used to analyze geophysical datasets. A critical and somewhat overlooked component of high-dimensional Bayesian inversion is the definition of the prior probability density function that describes the joint probability of model parameters before considering available datasets. If insufficient prior information is available about model parameter correlations, then it is tempting to assume that model parameters are uncorrelated. When working with a spatially gridded model representation, this leads to posterior realizations with far too much variability to be deemed realistic from a geological perspective. In this study, we introduce a new approach for structure-based prior sampling with Markov chain Monte Carlo that is suitable when only limited prior information is available. We evaluate our method using model structure measures related to standard roughness and damping metrics for l 1 -and l 2 -norms.

We show that our structure-based prior approach is able to adequately sample the chosen prior distribution of model structure. The usefulness and applicability of the methodology is demonstrated on synthetic and field-based ground penetrating radar data. We find that our method provides posterior model realizations and statistics that are significantly more satisfactory than those based on underlying assumptions of uncorrelated model parameters or on explicit penalties on model structure within an empirical Bayes framework.

Introduction

In Bayesian inference, the model parameters are treated probabilistically. Bayes theorem is used to combine a prior probability density function (pdf) and a likelihood function that quantifies the agreement between proposed model parameter values and the observed data [START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF]. The corresponding posterior pdf is only defined analytically for specific situations, such as, for linear inverse problems based on least-squares formalism. For non-linear problems, it is necessary to rely on sample-based methods, such as Markov chain Monte Carlo (MCMC) methods [START_REF] Gamerman | Markov Chain Monte Carlo, Stochastic Simulation for Bayesian Inference[END_REF]Gilks et al., 1996;[START_REF] Sambridge | Monte Carlo methods in geophysical inverse problems[END_REF][START_REF] Mosegaard | Monte Carlo sampling of solutions to inverse problems[END_REF]. Assumed or known spatial correlations are summarized either through parametrization choices (the number of unknowns and the type of basis function used; e.g., [START_REF] Lochbühler | Probabilistic electrical resistivity tomography of a CO 2 sequestration analog[END_REF] or by considering prior pdfs defined in terms of, for example, multi-Gaussian fields, multiple-point statistics or partially ordered Markov models [START_REF] Linde | Geological realism in hydrogeological and geophysical inverse modeling: a review[END_REF][START_REF] Hansen | Probabilistic integration of geoinformation[END_REF]. Such geostatistical priors are useful when dealing with well-characterized sites, or when the geological setting is well understood [START_REF] Huysmans | Application of multiple-point geostatistics on modeling groundwater flow and transport in cross-bedded aquifer (Belgium)[END_REF][START_REF] Caers | Multiple-point geostatistics: a quantitative vehicle for integrating geologic analogs into multiple reservoir models[END_REF]). Nevertheless, focusing on a specific class of spatial random field models might provide overly optimistic estimates of posterior model uncertainty and the morphology of individual posterior model realizations will be strongly affected by the prior pdf [START_REF] Hansen | Probabilistic integration of geoinformation[END_REF].

Another option is to simply state that if there is insufficient information about model parameter correlations, then it might be wiser to simply assume that model parameters are uncorrelated [START_REF] Linde | Distributed Soil Moisture from Crosshole Ground-Penetrating Radar Travel Times using Stochastic Inversion[END_REF]. Unfortunately, when working with a spatially gridded model of moderately high dimensions (some 100's of unknowns), this leads to posterior model realizations with far too much variability (the entropy is too high) to be deemed realistic from a geological perspective [START_REF] Linde | Distributed Soil Moisture from Crosshole Ground-Penetrating Radar Travel Times using Stochastic Inversion[END_REF][START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF]. In fact, actual smoothly varying subsurface fields will never (in practice) be part of the inferred posterior pdf. This is disturbing as the choice of ignoring parameter correlations is often motivated in terms of making the least assumptions possible.

To address this issue, [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF] used an explicit prior constraint on model structure using classical roughness measures that are common in deterministic inversions [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF][START_REF] Degroot-Hedlin | Occam's inversion to generate smooth, twodimensional models from magnetotelluric data[END_REF]. The appropriate weight to this prior constraint was estimated within an empirical Bayes (EB) framework [START_REF] Casella | An introduction to empirical Bayes data analysis[END_REF][START_REF] Malinverno | Expanded uncertainty quantification in inverse problems: Hierarchical Bayes ad empirical Bayes[END_REF]. [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF] applied the methodology to electrical resistance tomography (ERT) and radio magnetotelluric (RMT) data. Visually, the posterior realizations offered a clear improvement compared with an uncorrelated prior, but this came at the cost of model realizations that were overly smooth and marginal probability estimates that were slightly biased. [START_REF] Lochbühler | Summary statistics from training images as prior information in probabilistic inversion[END_REF] tackled the problem differently, by using prior pdfs of summary statistics (e.g., total model roughness) that were estimated from a large set of training images that describe the expected sub-surface structure. This approach leads to less bias provided that reliable training images are available.

In this study, we introduce a formal Bayesian method that defines the prior pdf in terms of global summary statistics describing the model structure. We consider uniform priors of model structure and highlight the main differences and advantages compared with the recent work by [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF] and [START_REF] Lochbühler | Summary statistics from training images as prior information in probabilistic inversion[END_REF]. We explain why structure-based spatial priors offer an useful complement to uncorrelated parameters or more traditional geostatistical models when addressing subsurface investigations with limited amounts of prior information. The examples considered are related to near-surface applications, but the approach is equally applicable to crustal or mantle-scale problems.

In order to sample a prior pdf of a given global measure of model structure, we rely on the Metropolis-Hastings (M-H) algorithm for computing the acceptance ratio within an MCMC framework [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling method using Markov chains and their applications[END_REF]. Given a target pdf that is hard to sample from and a model proposal pdf, the M-H algorithm generates a sequence of iterates whose distribution approximates the target pdf once the chain converges. This happens because the M-H acceptance ratio satisfies the two conditions that ensure that any Markov process have a unique stationary pdf: detailed balance and ergodicity.

Instead of considering the posterior pdf of the model parameters as the target pdf, we formulate the prior and its sampling in terms of global measures of model structure.

The numerical challenge of our approach lies in estimating the proposal pdfs in terms of the chosen model structure metric. Indeed, proposal pdfs that sample individual model parameters symmetrically will result in asymmetric proposal pdfs in terms of model structure and it is absolutely essential to account for this asymmetry. The absence of an analytical expression for the proposal pdf in terms of model structure requires that we sample it numerically and estimate it empirically. The numerical challenge is non-trivial and requires an accurate approximation of the ratio between the probability (in terms of model structure) of moving from the proposed state to the current one and of the corresponding reverse step. This ratio depends strongly on the current model structure (and the individual values of model parameters) so it must be evaluated at each MCMC step. After demonstrating the ability of our algorithm to adequately sample the desired prior pdf for various global model structure metrics, we apply the Bayesian inversion algorithm to a simple synthetic test case and to a field-based crosshole groundpenetrating radar (GPR) experiment.

In this paper, we first introduce structure-based prior sampling (SBP in Section 2.3) and present the results when considering MCMC inversion of a synthetic and a field-based crosshole GPR experiment (Section 2.4). In Section 2.5, we discuss the advantages and drawbacks of this methodology with respect to uncorrelated prior (UP) MCMC inversion and with the methodologies by [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF] and [START_REF] Lochbühler | Summary statistics from training images as prior information in probabilistic inversion[END_REF]. Finally, in Section 2.6 we provide our conclusions.

Methodology

Bayesian inversion

Within this study, we define the model parameters on a squared l × l Cartesian grids. The subsurface model, m, is then described by a vector of l × l = M variables, while the N observations are contained in the data vector d. These are related to the model through a forward operator g(-):

d = g(m) + , (2.1)
where is an N-dimensional vector containing measurement errors.

Bayes Theorem (Tarantola & Valette, 1982) states that the posterior pdf of the model parameters conditional to the observations is given by:

ρ(m|d) = ρ(m)ρ(d|m) ρ(d) , (2.2)
where ρ(m) is the prior pdf, ρ(d) is the evidence, which is a constant when the parametrization is fixed (as in our case), and L(m|d) ≡ ρ(d|m) is the conditional probability of the observations given the model, also called the likelihood function. This function quantifies the agreement between data predictions based on proposed model realizations and actual observations (the larger the value of the likelihood, the closer are the predictions and observations). Under the assumption of normally distributed and uncorrelated measurements errors, the log-likelihood is [START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF]:

l (m|d) = - N 2 l og (2π) - 1 2 l og ( N i =1 σ 2 i ) - 1 2 φ, (2.3) where φ = N i =1 g i (m) -d i σ i 2
represents the data misfit and σ i is the standard deviation of the i-th measurement. For non-linear inverse problems it is not possible to find an analytical formulation for the posterior pdf, ρ(m|d), and it is common to rely on MCMC methods [START_REF] Sambridge | Monte Carlo methods in geophysical inverse problems[END_REF].

Metropolis-Hastings algorithm

There exists many MCMC samplers, but we rely herein on the Metropolis-Hastings algorithm [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling method using Markov chains and their applications[END_REF]. In the standard formulations, the transition probability is proportional to the acceptance ratio α: .4) In this formulation, q is the proposal pdf, which generates model perturbations at each proposal step, while subscript prop refers to the proposed model state and curr to the current state of the Markov chain.

α = mi n 1, L(m pr op |d)ρ(m pr op )q(m pr op → m cur r ) L(m cur r |d)ρ(m cur r )q(m cur r → m pr op ) . ( 2
A basic characteristic of MCMC samplers is the so-called Markov property, which states that the probability of visiting the proposed state m pr op only depends on the current state m cur r and not on previously visited states. This implies that the sampling algorithm is completely defined by the transition matrix:

T (m pr op |m cur r ) ≡ P r (m cur r → m pr op ), (2.5) where P r (m cur r → m pr op ) is the probability of moving from the current state to the proposed state, and the transition matrix is such that T (m pr op |m cur r ) ≥ 0 and

m i T (m pr op |m cur r ) = 1.
The elements of this transition matrix are the transition probabilities and they define the sampling method used.

The Metropolis-Hastings algorithm is a powerful and widely used MCMC method to sample from unknown pdfs. Given a target pdf P (m) that is hard to sample from and a proposal pdf q(m cur r -→ m pr op ), the M-H algorithm generates a sequence of iterates {m 1 , m 2 , . . . , m T } which approximates a unique stationary pdf (π(m) = P (m)) after the burn-in period. The burn-in corresponds to the initial phase of the chain, in which the sampled states are still influenced by the arbitrary starting point and is, therefore, not considered part of the sampled posterior pdf. The Metropolis-Hastings algorithm is able to sample a unique stationary pdf, π(m), because it fulfills the following two conditions [START_REF] Robert | Introducing Monte Carlo Methods with R[END_REF]):

1. Detailed balance (or reversibility) is a sufficient condition for a random walk to asymptotically reach a stationary pdf. This principle requires each transition (m cur r -→ m pr op ) to be reversible and it can be stated as follows: given a transition matrix T (m pr op |m cur r ) (eq. 2.5), a stationary distribution P (m) satisfies T (m pr op |m cur r )P (m cur r ) = T (m cur r |m pr op )P (m pr op ).

(2.6) 2. Ergodicity of the Markov process requires that every state must be aperiodic (e.g., the system does not return to the same state at fixed intervals), positive recurrent (e.g., the expected number of steps for returning to the same state is finite) and irreducible (e.g., each state is accessible in a finite number of steps). This condition guarantees the uniqueness of the stationary pdf π(m).

A central aspect of MCMC theory is to define transition kernels, such that the sequence of samples drawn will converge to the target pdf π(m). Let us consider below the case of prior sampling, in which π(m) = ρ(m). If for some m cur r and m pr op we have that ρ(m cur r )q(m cur r -→ m pr op ) > ρ(m pr op )q(m pr op -→ m cur r ), (2.7) then this implies that the chain moves too often from m cur r to m pr op and too rarely in the other direction. To counteract this tendency, the Metropolis-Hastings algorithm reduces the number of moves from m cur r to m pr op , to achieve detailed balance. This is done by introducing a probability (0 < α(m cur r , m pr op ) < 1) that the proposed move is executed. The resulting transition from m cur r to m pr op is defined as:

T (m pr op |m cur r ) = q(m cur r -→ m pr op )α(m cur r , m pr op ).

(2.8)

The probability α(m cur r , m pr op ) is calculated to ensure that T (m pr op |m cur r ) satisfies the detailed balance criterion (eq. 2.6), as:

ρ(m cur r )T (m pr op |m cur r ) = ρ(m pr op )T (m cur r |m pr op ), ρ(m cur r )q(m cur r -→ m pr op )α(m cur r , m pr op ) = ρ(m pr op )q(m pr op -→ m cur r )α(m pr op , m cur r ).
(2.9)

The inequality in eq. ( 2.7) indicates that the move from m pr op to m cur r is not made often enough. Setting α(m pr op , m cur r ) = 1, one obtains from eq. ( 2.9):

ρ(m cur r )q(m cur r -→ m pr op )α(m cur r , m pr op ) = ρ(m pr op )q(m pr op -→ m cur r ), α(m cur r , m pr op ) = mi n 1, ρ(m pr op )q(m pr op -→ m cur r )
ρ(m cur r )q(m cur r -→ m pr op ) , (2.10) which is the acceptance ratio given in eq. ( 2.4) when no data are considered.

Structure-based priors

Most MCMC applications in geophysics rely on a symmetric proposal pdf: q(m pr op -→ m cur r ) = q(m cur r -→ m pr op ), which reduces the Metropolis-Hastings algorithm to the simpler Metropolis algorithm, in which the proposal pdf is not considered in the computation of the acceptance ratio. The acceptance ratio is then given by: .11) which reduces to likelihood ratios when considering uniform and uncorrelated prior pdfs:

α = mi n 1, L(m pr op |d)ρ(m pr op ) L(m cur r |d)ρ(m cur r ) , ( 2
ρ(m) = U [a, b] .
Our interest lies in sampling a given prior pdf of a global measure of model structure.

To do so, we need to account for the asymmetry of the proposal step q(S(m cur r ) -→ S(m pr op )), which arises when using a symmetric model proposal pdf on the individual model cells. The M-H ratio to sample ρ(S) is given by:

α = mi n 1, ρ(S(m pr op ))q(S(m pr op ) -→ S(m cur r )) ρ(S(m cur r ))q(S(m cur r ) -→ S(m pr op )) , (2.12)
where, in the numerator, ρ(S(m pr op )) is the prior probability of the chosen model structure metric, and q(S(m pr op ) -→ S(m cur r )) is the probability of proposing a given model structure when using an underlying symmetric proposal pdf q(m pr op -→ m cur r ) for the individual model parameters. The proposal pdf in terms of model structure is often asymmetric and depends strongly on m cur r in the denominator of the ratio, and on m pr op in the nominator. The absence of an analytical expression for the proposal pdfs in eq. ( 2.12), requires us to sample it numerically and estimate it empirically. This estimation must be repeated at each proposal step as the ratio q(S(m pr op ) -→ S(m cur r ))

q(S(m cur r ) -→ S(m pr op )) depends on m cur r and m pr op .

Our proposed structure-based implementation of the Metropolis-Hastings acceptance ratio (eq. 2.12), can be applied to any prior pdf of model structure (e.g., a Gaussian distribution with a known mean and standard deviation). Herein, we consider sampling from a uniform structural prior, but it is trivial to extend the method to other choices. A uniform structural prior is the least restrictive and it is the numerically most challenging case as it gives equal weights to all possible levels of model structure.

Empirical estimation of probability density functions

We consider two alternative approaches to sample the proposal probability ratio of model structure in eq. ( 2.12) from a finite sample of variables. In the non-parametric case, no assumptions are made about the underlying pdf, whereas in the parametric case samples are assumed to be drawn from a known distribution type (e.g., normal, gamma, beta...).

Kernel density estimation (KDE) is a frequently used non-parametric approach to estimate the pdf of a random variable [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. If the vector X = [X 1 , X 2 , . . . , X n ] is an independent and identically distributed sample drawn from a distribution with an unknown density f (x), its KDE is given by:

f (x) ≈ f (x) = 1 nh n i =1 K ( x -X i h ), (2.13)
where K is the kernel (or weighting function), which is usually chosen from symmetric functions of unitary area ( K (t )d t = 1), x is the point where the function needs to be evaluated and h > 0 is the bandwidth (or smoothing parameter). Both the kernel shape and the bandwidth dimension are crucial parameters and there is a vast literature on the topic [START_REF] Pedretti | An automatic locally-adaptive method to estimate heavily-tailed breakthrough curves from particle distributions[END_REF][START_REF] Shimazaki | Kernel Bandwidth optimization in spike rate estimation[END_REF][START_REF] Engel | An iterative bandwidth selector for kernel estimation of densities and their derivatives[END_REF][START_REF] Kuruwita | Density estimation using asymmetric kernels Bayes bandwidth with censored data[END_REF]. After conducting trials with different configurations for the case of a known underlying pdf, we found that the most suitable setting is to use a triangular kernel and the adaptive bandwidth (AB) method with sensitive parameter α = 0.5 [START_REF] Pedretti | An automatic locally-adaptive method to estimate heavily-tailed breakthrough curves from particle distributions[END_REF].

In the parametric approach, the underlying pdf is approximated by finding the bestfitting parameter values that define a specific distribution of a given pdf family. One way to estimate these parameter values is to use specific relations between the distribution's parameters and the sample's moments (method of moments; e.g., Frase 1958). In our case, preliminary investigations suggest that q(S(m cur r ) -→ S(m pr op )), as well as q(S(m pr op ) -→ S(m cur r )), are often asymmetric with both positive and negative skewness (e.g., the third standardized moment that quantifies asymmetry of the pdf ).

For the parametric approach, we fit the unknown pdf with the gamma distribution, which allows different degrees of symmetry (it also reduces to the normal distribution in the limit when the shape parameter goes to infinity):

f (x; k, θ) = 1 θ k Γ(k) x k-1 e -x θ for x > 0. (2.14)
In this equation, k > 0 represents the shape parameter and its value determines the skewness of the distribution, while θ > 0 is the scale parameter, whose value determines the statistical dispersion of the pdf. The gamma distribution is defined to have positive skewness, but negatively skewed pdfs can be estimated by reflecting the sampled data with respect to their median value. In implementing the method of moments we consider also an additional shifting parameter µ 0 . The support of the gamma pdf is [0, ∞], while the proposal pdfs we need to fit are defined on a smaller support ([a, b], with a > 0 and b < ∞). While the finite upper limit doesn't affect significantly the pdf estimation, the non-zero lower limit has a larger impact and requires the introduction of µ 0 in the parametric approximation of the pdfs. This additional shift parameter requires the pdf to be computed as f (x -µ 0 ; k, θ). The solution is obtained by solving three equations in three variables using direct relationships between the skewness and the shape parameter (eq. 2.15), the variance and the shape and scale parameters (eq. 2.16) and finally we introduce the shifting parameter in the mean definition (eq. 2.17): (2.16)

Skewness = 2 k , (2.15) V ar i ance = E (x -µ) 2 = kθ 2 ,
Mean = µ 0 + µ = µ 0 + kθ.
(2.17)

By computing the first three moments from the sample, we can obtain a unique estimate of the underlying gamma function.

Algorithm to compute the structure-based model proposal ratio

To ensure that the Metropolis-Hastings algorithm samples from a specified pdf of model structure, we need to estimate the proposal probability ratio (eq. 2.12). We use four steps to estimate the structure-based proposal probability, herein described for the forward step q(S(m cur r ) -→ S(m pr op )):

1. starting from m cur r = {m 1 , ...., m M }, draw the proposed model m pr op from the model proposal distribution (q(m cur r -→ m pr op ) in eq. 2.11) and compute the measure of structure S(m pr op ); 2. again from m cur r = {m 1 , ...., m M }, draw P new realizations from q(m cur r -→ m pr op ); 3. for each of the P new models, compute the corresponding model structure S p = S(m p ) ∀p ∈ {1, 2, . . . , P } ; 4. the obtained {S 1 , • • • , S P } are considered samples from the proposal pdf and we apply one of the methods described in section 2.3 to estimate the probability of draw S(m pr op ) from the sampled distribution.

The algorithm is equivalent for the backward step q(S(m pr op ) -→ S(m cur r )), but we do not need step (i) in this case, since we are interested in the probability of drawing the same amount of structure as S(m cur r ). For computational reasons, the sample size is kept relatively small (P = 1000). An example of generated samples and the corresponding estimated model structure-based proposal pdf is given in Fig. 2.1.

Measures of model structure

To illustrate our method, we consider four measures of model structure that are widely used in deterministic inversion [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF]. We use l 1 -and l 2 -norms related to model roughness and damping. To introduce these measures, we describe below each grid element (e.g., model parameter) by two indices. If we consider the model grid as a matrix, then the first index relates to the row and the second one to the column. With this notation the implemented measures of model structure are:

• Roughness (l 1 -norm): sum of the absolute differences between neighboring parameters

S R1 = |i -l |≤1 | j -k|≤1 |m i , j -m l ,k | (i , j ) = (l , k).
(2.18)

• Roughness (l 2 -norm): sum of the squared differences between neighboring parameters

S R2 = |i -l |≤1 | j -k|≤1 (m i , j -m l ,k ) 2 (i , j ) = (l , k). (2.19)
• Damping (l 1 -norm): sum of the absolute differences between each model parameter and a prior reference value (m r e f )

S D1 = i j |m i , j -m r e f |.
(2.20)

• Damping (l 2 -norm): sum of the squared differences between each model parameter and a prior reference value (m r e f )

S D1 = i j (m i , j -m r e f ) 2 .
(2.21)

Empirical Bayes

To evaluate the performance of our proposed structure-based prior MCMC method, we compare the results with the approach of Rosas-Carbajal et al. (2014), in which empirical Bayes is used to constrain model structure in favor of smooth spatial transitions. The empirical Bayes method [START_REF] Casella | An introduction to empirical Bayes data analysis[END_REF] uses hyper-parameters to describe a family of prior pdfs (e.g., the standard deviation of a normal distribution is treated as an unknown). The hyper-parameters are estimated together with the regular model parameters during the inversion. Consequently, they require both a prior and a proposal pdf.

Following [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF], the logarithmically constrained function on model structure is given by

l og (c m,2 )(m) = -Q log(2πλ 2 ) - S R2 (m) 2λ 2 , (2.22)
in case of a model structure defined in terms of roughness using an l 2 -norm, while in case of an l 1 -norm the expression is:

l og (c m,1 )(m) = -Q log(2λ) - S R1 (m) λ .
(2.23)

The corresponding equations for S D2 and S D1 are obtained by replacing S R2 and S R1 in eqs (2.22) and (2.23), respectively. Here, Q indicates the total rank of the model structure operator and λ is a hyper-parameter to invert for (together with other model parameters). The smaller the λ, the more weight is given to the prior constraint on model structure.

Estimating MCMC convergence

MCMC simulations are generally initiated at a random point in the prior parameter space, which is most likely far from regions of high posterior density. Therefore, the early stages of the Markov chains are generally not considered to be part of the target pdf. This initial part of the chain is referred as the burn-in period and indicates the part of the chain that is still dependent on the starting point. Estimating the burn-in time is necessary to decide on what part of the sampled chains that should be discarded.

We estimate the burn-in as the iteration when the estimated log-likelihood stabilizes at a value that correspond to the expected measurement uncertainty [START_REF] Cowles | Markov chain Monte Carlo convergence diagnostics: A comparative review[END_REF].

We rely on the potential scale reduction factor [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF] to estimate how many samples that are needed after burn-in to adequately sample the posterior pdf. In this study, we periodically compute this R-statistic using the last 50% of the samples within each chain. Following common practice, convergence to a stationary pdf is declared when R is less than 1.2 for all model parameters.

Results

To evaluate the influence of our proposed model structure-based prior, we apply our methodology to a synthetic test case and to a field example. In both cases, the dataset consists of first-arrival travel times from crosshole ground penetrating radar experiments. In this kind of experiments, information about radar wave speed (or its reciprocal: slowness) are obtained by recording the arrival times of radar waves transmitted in one borehole and received in an adjacent one. The forward problem consists of solving the eikonal equation using a finite-difference scheme [START_REF] Podvin | Finite difference computation of travel times in very contrasted velocity models: a massively parallel approach and its associated tools[END_REF].

Synthetic example

Our synthetic example considers a simple toy model on a 2D square domain with dimension of 5 m (Fig.

Prior structure-based sampling

We first evaluate the ability of our methodology to sample a uniform prior pdf for the four measures of structure defined in section 2.3: S R1 , S R2 , S D1 and S D2 . Our proposed structure-based prior Metropolis-Hastings algorithm is theoretically valid, but it is at the mercy of the finite size (P=1000) evaluations of the structure-based proposal distribution. To allow for ad hoc corrections to the numerical limitation of the empirical density estimation, we introduce an exponent correction ν to the computed ratio between the structure proposal pdfs (eq. 2.12). The acceptance ratio for a uniform structure-based prior sampling then becomes:

α = mi n 1, q(S(m pr op ) -→ S(m cur r ) q(S(m cur r ) -→ S(m pr op ) ν . (2.24)
For the prior sampling, the uncorrelated prior acceptance ratio is 1, while for the empirical Bayes case, we only consider the ratio of the constraint function. Finally, for our structure-based prior approach, the acceptance ratio is given by the estimated proposal probability ratio.

The sampled model structure for the different acceptance ratios are depicted for the parametric (Figs 2.3a-d) and non-parametric (Figs 2.3e-h) empirical structure proposal pdf estimation. When using the standard uncorrelated prior approach, we find (as expected) that the sampled model structure pdf is very narrow and do not sample low or high model structure values. The empirical Bayes method only samples models with very small model structure values. The results of our proposed structure-based prior method with ν = 1 in eq. ( 2.24) show a very significant widening of the sampled model structures. Adding a small empirically-estimated exponential correction (ν ∈ [1.02, 1.07]) leads to an even wider distribution (note that the y-axes in Fig. 2.3 are logarithmic).

Clearly, the structure-based prior approach does not perfectly match the ideal uniformstructure prior. This is perhaps not surprising given the numerical nature of our method to estimate the proposal ratios.

Inversion results

The synthetic dataset is generated by calculating the forward response of the synthetic model depicted in Fig. 2.2 and contaminating the simulated data with 1% Gaussian noise (a typical error level in crosshole GPR studies). For each model structure metric, we consider 5 Markov chains running in parallel with the structure-based prior Metropolis-Hastings acceptance ratio:

α SB P = mi n 1, exp(l (m pr op |d) -l (m cur r |d)) × q(S(m pr op ) -→ S(m cur r ) q(S(m cur r ) -→ S(m pr op ) ν , (2.25)
with l (m|d) being the log-likelihood function defined in eq. ( 2.3).

The results are compared with the uncorrelated prior approach: For our model structure-based prior approach, we first compare the inversion results with and without the exponent correction (e.g. ν = 1 or ν as in Fig. 2.3). The sampled posterior pdfs of model structure (Fig. 2.4) almost overlap, which suggests that the influence of the exponent correction is minor on the inferred posterior pdf. Therefore, in the following we kept ν = 1 in eq. ( 2.25).

α U P = mi n 1,
For the empirical Bayes approach, we follow Rosas-Carbajal et al. (2014) in defining the prior on the model structure weight (λ): Jeffreys pdf with a range between one quarter and four times the optimal value found when eqs (2.22) and (2.23) are maximized for the true model.

The results in terms of the sampled posterior model structure are shown in Fig. 2.5 for the three different acceptance ratios (eqs 2.25-2.27). For all measures of model structure, the uncorrelated prior overestimates the actual model variability severely. The empirical Bayes approach behaves differently depending on the norm used to quantify the model variability: the method is able to sample the correct model structure when using the l 1 -norm, while it underestimates the structure when using an l 2 -norm. For our structure-based prior approach, the actual model structure is sampled in all cases. Random posterior model realizations are shown in Fig. 2.7 in case of an S R2 measure of model structure. As for the prior realizations, the uncorrelated prior results (Fig. 2.7a) show too much model variability, while the empirical Bayes results (Fig. 2.7b) are overly smooth. In case of our structure-based prior approach (Figs 2.7c and d), the ability to sample the correct structure is reflected by model realizations that appear closer to the actual model (c.f., Fig. 2.2). Fig. 2.8 depicts the mean and standard deviation of the posterior realizations for the different acceptance ratios. Our structure-based prior approach shows mean values that are closer to the true model (c.f., Fig. 2.2) and provides standard deviations with magnitudes that are in-between the other two methods, with a slightly lower variability in case of non-parametric pdf estimation. Convergence analyses for the different acceptance ratios are summarized in Table 2.1, where we list the number of iteration after which R< 1.2 for all model cells, the relative data misfit root mean squared error (r mse D ), and the average acceptance rate after burn-in. Moreover, for estimating the quality of the posterior model realizations, we also computed the root mean squared error (r mse M ) between the actual model, m ac t (Fig. 2.2), and the posterior ones, m (Fig. 2.7):

r mse M = m ac t -m M 2 .
(2.28)

The values of r mse M (see Table 2.1) demonstrate that our structure-based prior approach samples models that are significantly closer to the underlying one compared with those sampled by the other methods (r mse M is 50% higher when using empirical Bayes and 100% higher when using the uncorrelated prior approach). Furthermore, the R-convergence is significantly reduced from uncorrelated to structure-based prior approach (uncorrelated prior R-convergence time is 1.3 to 2 times longer) and from uncorrelated prior to empirical Bayes (2.2 to 3.4 times). We also find that the r mse D values are smaller in case of the uncorrelated prior approach (5% smaller than empirical Bayes results and 2% smaller than when using our structure-based prior approach). The mean r mse D values for our structured-based prior approach (between 1.04 and 1.05) are very close to the r mse D calculated from the forward response of the true model (1.04). 3.9 × 10 5 25 1.05 1.9

Field example

We now apply our proposed methodology to a crosshole GPR experiment performed at a field site close to the Thur River in northern Switzerland [START_REF] Coscia | 3D crosshole ERT for aquifer characterization and monitoring of infiltrating river water[END_REF][START_REF] Klotzsche | 3D characterization of high-permeability zones in a gravel aquifer using 2D crosshole GPR full-waveform inversion and waveguide detection[END_REF]. This site has been investigated thoroughly, for example, to delineate lithologic sub-units within the gravel aquifer by joint inversion of seismic, ERT and crosshole GPR data [START_REF] Doetsch | Zonation for 3D aquifer characterization based on joint inversions of multi-method crosshole geophysical data[END_REF]. For the GPR experiments, the data were acquired at a 0.4 ns sampling rate using a RAMAC 250 MHz system, which at this site had a center frequency of 100 MHz with significant energy between 50-170 MHz. The first-arrival travel times ranged between 50 and 116 ns, with estimated picking errors of approximately 1%. The vertical spacing of sources and receivers was 0.5 m and 0.1 m, respectively. Moreover, in order to ensure fully symmetric radar coverage, source and receiver antennas were interchanged and the experiment repeated for each borehole pair.

We focus only on one borehole pair and N=878 first-arrival GPR travel times (e.g. 2D crosshole GPR inversion). The two boreholes considered are C2 and B2 (following the naming convention of [START_REF] Coscia | 3D crosshole ERT for aquifer characterization and monitoring of infiltrating river water[END_REF]. Fig. 2.9 shows the GPR inversion model obtained for this area by [START_REF] Doetsch | Zonation for 3D aquifer characterization based on joint inversions of multi-method crosshole geophysical data[END_REF] using a classical deterministic inversion.

Inversion results

For In agreement with the synthetic case, we find that the sampled posterior pdfs of model structure show a strong overlap, which indicates that the influence of the exponent cor-rection on the inferred posterior pdf is minor. The results based on the non-parametric method (c.f., Fig. 2.10d) are analogous. Similarly to the synthetic example (c.f., Fig. 2.5b), we find that the structure-based prior approach provide a posterior distribution with a wider class of model structure than the other methods. Random prior model realizations are shown in Fig. 2.11. As expected, we find a large model variability for the uncorrelated prior approach (Fig. 2.11a) and a high correlation between model parameters for the empirical Bayes approach (Fig. 2.11b). With our structure-based prior method (Figs 2.11c and d) it is possible to sample different degrees of model parameter correlation. In Fig. 2.12 we show random posterior model realizations. The uncorrelated prior results are too variable to be deemed geologically realistic (Fig. 2.12a), while the empirical Bayes realizations tend to provide very smooth transitions and almost no small-scale variability (Fig. 2.12b). Our structure-based approach provides an in-between behavior (Figs 2.12c and d), in which geologically realistic structure is imaged with more details. Fig. 2.13 depicts the means and standard deviations of the posterior realizations for the different acceptance ratios. Similar to the synthetic case, we find that our structure-based prior approach provides standard deviations that are in-between the other two methods, with a slightly lower variability in case of non-parametric pdf estimation.

In Table 2.2 we summarize the convergence analysis of the different MCMC inversion approaches. As in the synthetic case, the R-convergence for the uncorrelated prior is significantly larger than for the structure-based prior (1.3 times longer) and the empirical Bayes approach (1.5 times longer). Finally, the relative r mse D values are smaller in case of the uncorrelated prior approach (30% smaller than for the empirical Bayes results and 10% smaller than for the structure-based prior approaches).

Discussion

We have introduced a new approach for structure-based prior sampling in MCMC inversion. The measures of model structure considered are related to standard roughness and damping measures based on l 1 and l 2 -norms. We have shown that the introduced approach is able to adequately sample the prior distribution of model structure. We focused on uniform structural priors to give equal weight to all possible levels of model structure, but our approach can easily be applied to other prior pdfs of model structure (it is enough to assign the prior probabilities to ρ(S(m pr op )) and ρ(S(m cur r )) in eq. 2.12). To quantify the efficiency of the structure-based prior method, we compared it with the uncorrelated prior and empirical Bayes approach (e.g., Rosas-Carbajal et al.

2014).

We presented examples for synthetic and field-based crosshole GPR data, but the methodology can be applied to any kind of geophysical data set or spatial scale.

Figs 2.7 and 2.12 show that the introduced structure-based prior approach returns posterior model realizations that are more similar to the underlying model (in the case of synthetic dataset) or to the expected subsurface properties (in the case of field data) than results based on an uncorrelated prior or empirical Bayes. This happens as our structure-based prior approach samples a larger range of model structure. In both empirical Bayes and uncorrelated prior methods, the range of the sampled model structure is limited: in the first case towards higher values (e.g., more heterogeneous models) and in the second case towards lower ones (e.g., more homogeneous models).

When the actual structure of the subsurface is outside of the sampled ranges, then the two approaches are practically unable to sample the true model. To illustrate this, Fig. 2.14 shows that when using an uncorrelated prior on a 10 × 10 discretized grid, the a priori probability to sample a structure value that is half of the sampled mean is approximately 10 -100 . The structure-based prior approach scans a much larger range of model structure (see Figs 2.3,2.10a and b). Conceptually closer to our structure-based prior approach is the technique introduced by [START_REF] Lochbühler | Summary statistics from training images as prior information in probabilistic inversion[END_REF]. They formulate an indirect prior on model parameters by using statistics on model structure values derived from training images. This approach uses an indicator function to ensure that the Markov chains sample model structures within the range found in the training images realizations. Nevertheless, this method does not ensure detailed balance in sampling the model structure, which leads to inversion results that are slightly biased towards higher model structure values.

Our structure-based prior approach is built on solid theory, but its performance is affected by the numerical implementation. The main drawback is related to the empirical pdf estimation schemes used to compute the proposal structure distribution ratio. We have found that using a small exponent correction (see eq.

Conclusions

The prior pdf has a strong influence on moderately high-dimensional (100s of parameters) Bayesian inversion results. Despite its importance, there is little guidance in the geophysical literature about how to choose meaningful priors when prior knowledge is weak. We argue that prior models that are described in terms of global model structure are easy to define and understand, especially for researchers with a background in classical regularized deterministic inversion. We propose a theoretically solid approach to enable prior sampling in terms of a pre-defined model structure metric. The method is demonstrated for l 1 -and l 2 -norm measures related to model roughness and damping, but the framework can be applied to any metric of model structure. For both synthetic and field-based crosshole GPR data, we find that our structure-based prior approach provides posterior model realizations and statistics that better resemble the expected subsurface properties than those obtained by assuming that individual model parameters are uncorrelated or from empirical Bayes implementations that favors smoothly varying models. The proposed framework is applicable to 3D inversion, other geophysical data types and scales.

Chapter 3

Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data Giulia de Pasquale, Niklas Linde, Joseph Doetsch and W. Steven Holbrook.

Accepted1 in Geophysical Journal International and herein slightly adapted for this thesis.

ing heat exchange with the atmosphere and indirectly as its melting is responsible for greenhouse gas emissions [START_REF] Dafflon | Geophysical estimation of shallow permafrost distribution and properties in an ice-wedge polygon-dominated Arctic tundra region[END_REF]. Geophysical detection and monitoring of this interface is most commonly achieved using seismic, electromagnetic and electrical methods [START_REF] Minsley | Airborne electromagnetic imaging of discontinuous permafrost[END_REF]. Yet another example is landslide studies, in which near-surface geophysics plays a central role when investigating slip and failure surfaces [START_REF] Jongmans | Geophysical investigation of landslides: a review[END_REF] as they often correspond to the contact between regolith and bedrock, with its contrast in physical properties. Imaging with focus on the depth to the bedrock and its topography can also provide useful constraints in the study of landslide triggering mechanisms [START_REF] Lanni | Simulated effect of soil depth and bedrock topography on near surface hydrologic response and slope stability[END_REF][START_REF] Reneau | Size and location of colluvial landslides in a steep forested landscape[END_REF][START_REF] Sanchez | Relationships between tectonics, slope instability and climate change: Cosmic ray exposure dating of active faults, landslides and glacial surfaces in the SW Alps[END_REF]. Moreover, bedrock mapping is important when characterizing the critical zone (CZ): the near-surface and external layers of the Earth that include most of the life-sustaining resources on our planet (Brantley et al., 2011). Here, knowledge about the depth to the regolith-bedrock interface (i.e., the lower boundary of the critical zone) is necessary to model the balance between soil formation and erosion, which is important for understanding landscape evolution, sediment transport and budgets (St. Clair et al., 2015;[START_REF] Heimsath | The soil production function and landscape equilibrium[END_REF][START_REF] Rempe | A bottom-up control on fresh-bedrock topography under landscapes[END_REF][START_REF] Parsekian | Multiscale geophysical imaging of the critical zone[END_REF]. Non-invasive geophysical investigation techniques have been increasingly applied to map regolith-bedrock interfaces and, more specifically, electrical resistivity tomography (ERT) has been commonly used because of its comparatively high resolution over the relevant depth range and the strong sensitivity of electrical conductivity to variations of relevant hydrogeological and geological subsurface properties [START_REF] Chambers | Bedrock detection beneath river terrace deposits using three-dimensional electrical resistivity tomography[END_REF][START_REF] Saas | Bedrock detection and talus thickness assessment in the European Alps using geophysical methods[END_REF].

Even though geophysical datasets are sensitive to such subsurface discontinuities, inferred interface locations and geometries are usually obscured by common inversion routines that seek smoothly-varying subsurface models. Geophysical inversion methods seeks to provide quantitative information about physical properties from indirect geophysical observations. It is generally an ill-posed problem, which implies that a solution based on data alone is neither stable nor unique [START_REF] Kabanikhin | Definitions and examples of inverse and ill-posed problems[END_REF]. Hence, when working with finely discretized property fields it is necessary to include a model regularization term (deterministic approach) or defining a prior probability density function (pdf) on the model parameters (Bayesian inference). In deterministic inversions, the regularization term often imposes smooth spatial variations of the Earth's subsurface by penalizing a measure of model variability [START_REF] Constable | Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data[END_REF][START_REF] Degroot-Hedlin | Occam's inversion to generate smooth, twodimensional models from magnetotelluric data[END_REF]. Because of the nature of the constraints and the non-constraining nature of the data, sharp interfaces are generally smeared out within a thick transition zone. In practice, this implies that it is up to the interpreter to identify the most likely interface geometry [START_REF] Christensen | Interpretation attributes derived from airborne electromagnetic inversion models using the continuous wavelet transform[END_REF].

To decrease such subjectivity, [START_REF] Hsu | Bedrock detection using 2D electrical resistivity imaging along the Peikang River, central Taiwan[END_REF] propose an automatic approach to bedrock edge detection, which analyzes 2-D smoothness-constrained deterministic inversion models obtained from an ERT dataset using image analysis. This is achieved by a gradient method (i.e., Laplacian edge detection), which searches for zero values in the second derivative of the resistivity image in both the horizontal and vertical direc-tion. The main issue with such a methodology is the possible interference from other sources of subsurface heterogeneity that are unrelated to the large-scale interfaces of interest. To reduce the complexity of such a methodology, especially when considering 3D subsurface models, [START_REF] Chambers | Bedrock detection beneath river terrace deposits using three-dimensional electrical resistivity tomography[END_REF] only consider the maximum resistivity gradient in the vertical direction to locate the bedrock beneath a river terrace. Another approach that uses a guided fuzzy clustering algorithm has been proposed by [START_REF] Ward | Distributionbased fuzzy clustering of electrical resistivity tomography images for interface detection[END_REF]. Also in this case, the edge detection of bedrock interfaces is based on analysis of ERT inversion results.

Instead of imposing smoothness by model regularization, it might be more reasonable to postulate that physical properties of the subsurface are characterized by a layered structure (homogeneous layers of varying topography with distinctly different properties).

Clearly, targeting the inversion towards the automatic detection of interfaces instead of producing smooth minimum-structure models and then, automatically or manually, identifying the expected boundaries from these smoothly varying images has its advantages. Within such a framework, Auken & Vest Christiansen (2004) propose a twodimensional deterministic inversion scheme with lateral smoothness constraints and sharp vertical boundaries. Similarly, [START_REF] Juhojuntti | Joint inversion of seismic refraction and resistivity data using layered models -Applications to groundwater investigation[END_REF] propose a method for joint inversion of seismic refraction and resistivity data using sharp-boundary models with few layers. Hence, this method produces models with laterally-varying properties within the layers, but without vertical variations within each layer. Another scenario occurs when interface locations are well known, but the fields of physical layer properties are of interest. In this case, it is possible to pre-determine unit geometries and infer physical properties and variations within each unit [START_REF] Doetsch | Constraining 3-D electrical resistance tomography with GPR reflection data for improved aquifer characterization[END_REF][START_REF] Uhlemann | Four-dimensional imaging of moisture dynamics during landslide reactivation[END_REF][START_REF] Vignoli | Examples of improved inversion of different airborne electromagnetic datasets via sharp regularization[END_REF].

In this paper, we propose a probabilistic formulation and solution to the geophysical inverse problem of inferring an interface separating two heterogeneous sub-domains. We rely on Bayes' theorem to combine the prior distribution about model parameters with a likelihood function that describes the probability of observing the collected data given proposed model parameter values and a data noise model [START_REF] Sambridge | Monte Carlo methods in geophysical inverse problems[END_REF]. The resulting posterior distribution is generally not known analytically. Instead, it must be sampled numerically using, for instance, Markov chain Monte Carlo (MCMC) methods [START_REF] Sambridge | Monte Carlo methods in geophysical inverse problems[END_REF]. When performing MCMC inversion for highdimensional distributed models, the choice of the prior pdf has a very strong influence on inversion results [START_REF] Hansen | Inverse problems with non-trivial priors: efficient solution through sequential Gibbs sampling[END_REF][START_REF] Hansen | Probabilistic integration of geoinformation[END_REF][START_REF] De Pasquale | Changing the Prior Model Description in Bayesian Inversion of Hydrogeophysics Dataset[END_REF]. To favor posterior model realizations with smoothly varying properties, a possible approach is to include structure constraints [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF][START_REF] Besag | Bayesian computation and stochastic systems[END_REF][START_REF] Chen | Stochastic inversion of magnetotelluric data using sharp boundary parametrization and application to geothermal site[END_REF] within an empirical Bayes (EB) inversion framework [START_REF] Casella | An introduction to empirical Bayes data analysis[END_REF].

Here, we seek to probabilistically infer subsurface interface geometry and heterogeneities within the sub-domains that are delimited by the interface. To do so, we follow [START_REF] Iglesias | Well-posed Bayesian geometric inverse problems arising in subsurface flow[END_REF] and implement an empirical-Bayes-within-Gibbs methodology that separates the interface and physical property updates within the MCMC scheme. More specifically, we address the particular problem of using ERT datasets to infer the regolith-bedrock interface, at a site in the Calhoun Critical Zone Observatory in South Carolina, USA. Both the interface and the physical properties of the sub-domains are constrained to favor smooth spatial variations and to honor pre-defined property bounds. We demonstrate our methodology using synthetic and actual surface-based ERT datasets. We stress that the methodology is general and that it could be used to study other types of interfaces (see examples above) and that other types of geophysical methods (e.g., seismic refraction, magnetotellurics) could be used.

We first introduce the theoretical background of the proposed empirical-Bayes-within-Gibbs inversion algorithm (Section 3.3), followed by the results obtained when applying it to two synthetic and a real dataset (Section 3.4). We then discuss our approach and findings (Section 3.5) before we conclude (Section 3.6).

Methodology

To address the challenge of probabilistic inference of subsurface interfaces and subdomain heterogeneities, we present an inversion algorithm which combines modelstructure constraints within an empirical Bayes framework with Metropolis-within-Gibbs MCMC sampling.

Structure-constrained empirical Bayes inversion

By combining a prior probability density function (pdf) of the model parameters (ρ(m)), summarizing our a priori information about the subsurface, and the likelihood function, which expresses the probability of the proposed model vector (m) given the available data (d): ρ(d|m) ∝ L(m) [START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF], Bayes Theorem expresses the posterior pdf of the model parameters given the data:

ρ(m|d) = ρ(m)L(m) ρ(d) , ( 3.1)
where ρ(d) only acts as a normalizing constant when the model parametrization is fixed. Analytical solutions to equation (3.1) are not available for non-linear problems and it is thus necessary to numerically sample from the posterior distribution. To do so, MCMC methods are often used to construct Markov chains in the model space, whose steady state distribution corresponds to the posterior distribution [START_REF] Gamerman | Markov Chain Monte Carlo, Stochastic Simulation for Bayesian Inference[END_REF], with the sampling algorithm being completely defined by the transition probabilities. Metropolis-Hastings (M-H) is a commonly used MCMC sampler [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling method using Markov chains and their applications[END_REF], which defines the transition probability from different states of the chains as being proportional to the acceptance ratio:

α M H = mi n 1, L(m pr op )ρ(m pr op )q(m pr op → m cur r ) L(m cur r )ρ(m cur r )q(m cur r → m pr op ) , (3.2) where q is the proposal pdf describing the probability to generate a given model perturbation at each proposal step, and the superscripts prop and curr refer, respectively, to the proposed and the current state of the Markov chain. Acceptance ratios used in MCMC codes are typically reformulated using a logarithmic formulation to avoid floating-point under-and overflow. Most MCMC applications in geophysics rely on a symmetric proposal pdf, q(m pr op → m cur r ) = q(m cur r → m pr op ), which leads to the simpler Metropolis acceptance ratio:

α M = mi n 1, L(m pr op )ρ(m pr op ) L(m cur r )ρ(m cur r ) , (3.3)
for which there is no evaluation of the proposal pdfs at each step of the chain.

When the prior knowledge about the subsurface is limited, it is common to consider uncorrelated parameters with uniform probability over a given parameter range. This choice of priors leads to posterior model realizations that are too spatially variable to be geologically realistic when applied to high-dimensional distributed models [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF][START_REF] Hansen | Inverse problems with non-trivial priors: efficient solution through sequential Gibbs sampling[END_REF][START_REF] Hansen | Probabilistic integration of geoinformation[END_REF][START_REF] De Pasquale | On structure-based priors in Bayesian geophysical inversion[END_REF]). An alternative is to also constrain global model structure by penalizing model variability. The appropriate weight given to these constraints can be obtained through an empirical Bayes approach.

Empirical Bayes (EB) inference [START_REF] Casella | An introduction to empirical Bayes data analysis[END_REF] uses hyper parameters describing the prior distribution (e.g., its standard deviation, integral scales). The hyper parameters are then estimated together with the regular model parameters during the inversion process. When using MCMC, the hyper parameters require both prior distributions and proposal pdfs. For the pdf describing model structure used herein, the hyper-parameter λ defines the mean deviation of the exponential model:

c(S(m), λ) = (2λ) -Q exp - ||S(m)|| 1 λ , (3.4)
where Q indicates the total rank of the model structure operator and the measure of model structure S(m) is computed with an l 1 -norm (a corresponding formulation is used to define interface roughness) as further explained in section 3.3. Similar to the model regularization weights used in deterministic inversion, a high λ strongly penalizes model structure. Considering symmetric model proposal pdfs for the model parameters and the hyper-parameter, the acceptance ratio in case of structure-constrained empirical Bayes MCMC inversion is:

α E B = mi n 1,
L(m pr op )ρ(m pr op )ρ(λ pr op )c(S(m pr op ), λ pr op ) L(m cur r )ρ(m cur r )ρ(λ cur r )c(S(m cur r ), λ cur r ) , (3.5) where ρ(λ) is the prior pdf of the hyper parameter.

In this study, we considered ρ(λ) and ρ(m) to have log-uniform and uncorrelated uniform pdfs, respectively, such that the acceptance ratio simplifies to the ratio involving likelihood and constraint functions only.

Empirical-Bayes-within-Gibbs

We conceptualize the subsurface as being composed of two overlapping domains, m C Z and m b , that are separated by an interface, I. We first draw two uncorrelated physical property fields that are parametrized with regular grids that cover the whole investigated area (Fig. 3.1a). In the examples that follow, we used 50 × 10 grid with cell sizes of 8 × 8 meters. Herein, the physical properties considered are the logarithm of electrical resistivity r [Ωm], that is discretized throughout the sub-domains representing the critical zone (CZ, m C Z ) and the bedrock (m b ). The forward operator is defined on an unstructured mesh that covers the investigated area (Fig. 3.1b) and extends on the sides and below to account for boundary conditions. In the examples, the unstructured mesh used for forward simulations is discretized with 1578 cells. We linearly interpolate the sampled electrical resistivity values into the unstructured mesh (Fig. 3.1c). In parallel, we draw an interface from the corresponding prior pdf, ρ(I).

In 2D, the interface defines a line delineated by a set of connected nodes within the mesh (Fig. 3.1d). The actual resistivity field (m) used for forward simulations is built by mapping each of the interpolated physical property values into the appropriate sub-domain defined by the interface (Fig. 3.1e):

m = G(m C Z , m b , I), (3.6)
where G is the mapping operator.

Iglesias et al. (2014) introduced a Metropolis-within-Gibbs method in which the model proposals are symmetric and the interface and physical properties are updated alternately within the Gibbs framework [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF]. In our implementation, we additionally constrain the interface and the resistivity fields to favor smoothness, resulting in a slightly modified empirical-Bayes-within-Gibbs algorithm in which we alternately propose an update to the interface or to the resistivity field of the domains (Fig. 3.2 shows a flowchart of our algorithm).

The interface proposal updates are implemented as follows. Initially, one vertex on the current interface is chosen at random and removed from the interface definition. Then, one of its two adjacent vertices in the current interface is picked at random. When this vertex (anchor) is located on the left, one of its adjacent vertices in the mesh is drawn at random under the constraint that it is located to the right and that it is not the previously removed vertex. From this point, the shortest path is sought to a vertex belonging to the current interface under the constraint that all horizontal increments in the path are to the right. The vertices of this path are now part of the proposed interface geometry. When the anchor is located to the right, the model proposal procedure is adapted in a straightforward manner. In this way, we ensure symmetry in our model proposals and avoid a possible situation with the bedrock overlying itself. The probability of accepting the proposed interface is given by: The electrical resistivities are described using priors defined on different ranges:

α I = mi n 1, L(m cur r C Z ,
ρ(m C Z ) = ρ(m b
), as bedrock is generally more resistive than the overlying CZ. To update model properties, we choose a certain fraction of model parameters to update, randomly disposed within the domain, and sample the new values from a Gaussian pdf centered on the previous value. Moreover, with the aim of augmenting the posterior model space exploration, we alternate the choice of step size for the proposal distributions by switching between applying a randomly drawn standard deviation value from a Gaussian pdf and applying it uniformly over the whole domain or scaling it inversely to the square root of the cumulative sensitivity of the investigated subsurface. Model proposals that use the same standard deviation throughout the model domain help to converge to the target misfit, while the proposals based on the cumulative sensitivities allows for improved exploration in regions of low sensitivity. The cumulative sensitivities are here obtained from the model result of a deterministic inversion, but the outcomes are expected to be very similar if using the sensitivities obtained from a uniform model. Indeed, such a scaling of the standard deviations does not impact the symmetry and stationarity of the proposal distribution throughout the MCMC procedure, hence, it is a theoretically-valid model proposal. The model proposals are implemented for both domains at the same time and they can involve parameters that are not mapped into the resistivity field used for the forward computation (e.g., a parameter describing CZ resistivity at the bottom of the model domain). The corresponding acceptance ratio for the physical properties update is: First, we initialize the model (as described in Fig. 3.1) and evaluate its likelihood, physical properties structure and interface roughness. We then sequentially evaluate proposed interface and physical property updates in an iterative scheme until the chosen length of the chain is reached. Here, u ∼ U (0, 1) represents a randomly chosen number with uniform probability between 0 and 1, while the acceptance ratios α I and α P are described in eqs. 3.7 and 3.8, respectively.

α P = mi n 1, L(m

Measures of model structure and interface roughness

As mentioned above, we consider two different measures of model structure: one quantifying the spatial variability of the subsurface resistivity field (i.e., S(m) in eq. 3.8) and another quantifying the roughness of the interface (i.e., R(I) in eq. 3.7). In the first case, we consider a measure of model structure that is widely used in deterministic inversion [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF], namely l 1 -norms of model roughness, and we apply it to the inversion parameter grid (i.e., the regular grid parametrization represented in Fig. 3.1a). We follow de Pasquale & Linde (2017) by describing each grid element (i.e., model parameter) with two indices: m i , j . Considering the grid analogously to a matrix, the measure of electrical property variability can be expressed as:

S(m) = |i -l |≤1 | j -k|≤1 |m i , j -m l ,k | |i -l | + | j -k| = 1. (3.9)
To quantify the interface roughness, we use a corresponding l 1 metric of the differences between z-coordinates (e.g., depth) of the interface vertex, according to the following equation:

R(I) = n-1 l =1 |z l +1 -z l -1 |, (3.10)
where l = 0 . . . n refers to the nodes intersected by the interface I, with index l growing with the profile length (i.e., from left to right in Fig. 3.1d). Initial tests (not shown) suggested superior inversion results when considering a central difference measure compared with a forward difference measure.

ERT forward modeling

In the following test cases, we consider 2D electrical resistivity tomography datasets. Electrical resistivity tomography is widely used to determine the spatial distribution of electrical resistivity. In practice, apparent resistivity data are obtained using a large number of resistance measurements made from spatially-distributed four-point electrode configurations [START_REF] Dahlin | 2D resistivity surveying for environmental and engineering applications[END_REF]. In order to simulate the ERT forward response needed to compute the likelihood of each model realization within the Markov chain, we used the Boundless Electrical Resistivity Tomography (BERT) library, which is an open source software package widely used for deterministic inversion [START_REF] Rücker | Three-dimensional modelling and inversion of dc resistivity data incorporating topography-I[END_REF]Günther et al., 2006). Using BERT, we can compute the 2.5D forward response by relying on a finite element scheme implemented on unstructured meshes (Si, 2015). Using unstructured meshes allow for a more efficient local refinement (e.g., close to electrodes) and are ideally suited to account for surface topography and internal boundaries.

Results

To evaluate the performance of our proposed algorithm, we consider a 2D surface-based ERT dataset acquired above a regolith-bedrock interface. More specifically, we refer to the data from the Calhoun Critical Zone Observatory (South Carolina, USA) presented by St. Clair et al. (2015), in which ERT and seismic refraction surveys were used to investigate the influence of surface topography on bedrock geometry and properties.

The ERT survey was acquired using an Advanced Geosciences Inc's (AGI) SuperSting instrument with 56 electrodes spaced 5 m apart with a dipole-dipole electrode configuration. The profile length of 400 m was obtained by using one roll-along in which 50 % of the electrodes stayed in place. Consequently, St. Clair et al. (2015) used 84 different electrode positions. To save computational time and since our interest is in the deeper CZ structure, we only consider every second electrode position and we removed the configurations with a maximal electrode separation of 30 m. Moreover, the raw data were filtered by St. Clair et al. (2015) to eliminate negative apparent resistivity values. This resulted in 42 electrodes and 645 data being used. A noise description based on a relative Gaussian error of 3.7 % was used in order to reproduce the same data misfit as the inversion results by St. Clair et al. (2015).

Synthetic test cases

To test the methodology in a controlled environment, we first consider two synthetic examples with the same electrode configuration and surface topography as for the real dataset. Both examples present the same interface geometry, but the first example features homogeneous sub-domains and the second has heterogeneities within the layers.

To demonstrate that the empirical-Bayes-within-Gibbs algorithm can sample the target interface, we test the method on the model shown in Fig. 3.3(a). In this case, the interface divides the domain in two homogeneous domains: an upper layer with resistivity of 1800 Ωm and a lower layer with resistivity of 15000 Ωm. Fig. 3.3(b) represents the interpolation (through the nearest neighbor method) of the input resistivity model on the mesh used for the forward computations within the MCMC chains. Note that this mesh is different from the mesh used to generate the synthetic observations. The best interface we can sample is represented by a black line in Figs. 3.3(c) and (e).

We contaminate the simulated observations with 3.7 % uncorrelated Gaussian noise. The mean model sampled by the empirical-Bayes-within-Gibbs chain in Fig. 3.3(c) provides a very close approximation of the interface; moreover, the posterior realization corresponding to the maximum likelihood (Fig. 3.3d) exactly reproduces the interface of the interpolated target model. We contrast these results with the deterministic inversion result, in which the interface is obtained from the maximum vertical gradient method (i.e., the blue line in Fig. 3.3e ) using the approach described in [START_REF] Chambers | Bedrock detection beneath river terrace deposits using three-dimensional electrical resistivity tomography[END_REF].

From these results, it is evident that the empirical-Bayes-within-Gibbs approach is more successful in sampling the target interface. To quantify this, we compute the average distance between the target interface (I t ar g et ) and the sampled (in case of empirical-Bayes-within-Gibbs) or inferred one (in case of deterministic inversion result):

D = x max
x mi n |I t ar g et -I sampl ed /i n f er r ed | x maxx mi n , (3.11) where x mi n and x max are the horizontal limits of the model domain. For the deterministic inversion this distance is 2.9 m, while when evaluating the measure on a sample of posterior realizations, we find that the empirical-Bayes-within-Gibbs algorithm yields an average distance of 0.4 ± 0.1 m only. After verifying that the algorithm samples the target interface for this simple example, we introduce significant heterogeneities within the sub-domains. In Fig. 3.4(a), we show the corresponding resistivity model used to generate the new synthetic dataset together with the electrode positions. The underlying mesh is constructed to exactly accommodate the defined interface, while the two heterogeneous sub-domains are obtained by stationary multivariate Gaussian process generations through circulant embedding of the covariance matrix (Dietrich & Newsam, 1997). Fig. 3.4(b) represents the interpolation (through nearest neighbor method) of the input resistivity model on the mesh used for the forward computations within the MCMC chains. This image gives an idea of the best possible representation of the interface that can be sampled within the empirical-Bayes-within-Gibbs inversion. As before, the simulated data were subsequently contaminated with 3.7 % uncorrelated Gaussian noise (Fig. 3.4c).

Within our inversion routine, we used intervals for the electrical resistivity fields according to typical ranges of CZ and bedrock in Fig. These values were chosen to ensure an acceptance rate close to 25 % and to accommodate both smaller and larger model updates. We alternate between using a uniform value as standard deviation for all the model parameter updates and reweighing them inversely to the square root of the cumulative sensitivities (Fig. 3.4d). Based on extensive testing, we find that alternating model updates in this way enables improved exploration in regions of low sensitivity while still ensuring that the residuals of the model predictions have the same chi-square mean as the assumed noise term (i.e., that we reach the target misfit). The prior on the interface is defined by the way the interface can be constructed within the unstructured mesh used for computing the model response: nodes defining an interface have to be connected between each other through edges of the forward mesh and the horizontal increment must be either always positive or always negative within one realization (i.e., as in Fig. The "posterior" distributions are constructed based on the sampled realizations within the stationary part of the chains (i.e., after burn-in when the chains start to sample proportionally to the target distribution). To evaluate the burn-in length for each chain, we rely on the Geweke diagnostic [START_REF] Geweke | Evaluating the Accuracy of Sampling-Based Approaches to the Calculation of Posterior Moments[END_REF]. This method proceeds by testing if the mean of the log-likelihood of the first part (20 % in our implementation) of the supposedly stationary section of the MCMC chain can be assumed to be the same as the later part (last 50 % in our implementation). By incrementally increasing the initial part of the chain that is discarded from the analysis, it is possible to estimate the burn-in period. Here, this diagnostic results in a burn-in of 3 × 10 5 , 2.5 × 10 5 and 4 × 10 5 iterations for the first, second and third chain, respectively. Figs. 3.6(a) and (b) show the sampled pdfs of model structure for CZ and bedrock resistivity fields, respectively, while Fig. 3.6(c) represents the interface roughness. The marginal distribution of the hyperparameter used as structure-constraint weights are depicted in Figs. 3.6(d)-(f), together with the optimal values for the true model. The mean acceptance rate for the interface updates (Fig. 3.6g) is lower (average value of 8 %) than for the resistivity field updates (Fig. 3.6h; average value of 25 %). Finally, Fig. 3.6(i) represents the log-likelihood evolutions for the three chains. The potential scale reduction factor (?) is commonly used to assess if a sufficient number of posterior samples have been obtained in order to adequately describe the posterior distribution. For each model parameter, it compares the between-chain and the within-chain variance. In practice, it is generally agreed that a value of 1.2 is sufficient to declare convergence. Considering the last 5 × 10 5 iterations, the potential scale reduction factor is found to be below 1.2 for only 16 % of the CZ parameters and 45 % of the bedrock parameters. Therefore, these results indicate that we did not sufficiently sample the posterior pdf.

Nevertheless, we have seen that the different chains produce similar posterior model realizations (Fig. 3.5). The mean values (Figs. 3.7a, d and g) and normalized standard deviations of electrical resistivity (i.e., divided by the mean resistivity values; Figs. 3.7b, e and h) are comparable for all the chains. This strong similarity is also seen in the probability maps of the interface locations (Figs. 3.7c, f and i), and in the vertical resistivity profiles (Figs. 3.7l,m and o). For all chains we see that the empirical-Bayes-within-Gibbs samples the underlying interface well, except for the left side of the domain (first 50 m of the profile length) and around 370 m along the profile, where the inferred bedrock interface is too shallow. Moreover, all the three chains show the highest normalized standard deviation in the vicinity of the inferred interface. In this region, the interface is varying throughout the MCMC chain and the same cell is alternately associated with the CZ (more conductive) or bedrock (more resistive).

In order to evaluate the quality of the sampled posterior model realizations, we computed the model discrepancy (Fig. 3.8b): where m M is the vector of the log 10 resistivity values of the true model (Fig. 3.6b) and m is the mean of the posterior model realizations (Fig. 3.8a) for the three Empirical-Bayes-within-Gibbs chains. For comparison purposes, we also present the inversion model result obtained by a traditional l 1 -norm smoothness-constrained deterministic inversion that is fitted to the same error level (Fig. 3.8c) and the corresponding model discrepancy (Fig. 3.8d). We find that the MCMC inversion strongly overestimates the resistivity of the upper 40 m at the beginning of the profile, while the deterministic inversion consistently underestimates the resistivity of the bedrock, this is particularly visible in the vertical resistivity profiles represented in Figs. 3.8(e)-(h). The mean absolute discrepancy in Fig. 3.8(b) is lower (0.19) than for the deterministic inversion (0.32 in Fig. 3.8d). Fig. 3.8(c) shows also the interface obtained from the deterministic inversion using the maximum vertical gradient method. Also in this case, it is seen that the estimated interface is the worst on the left side of the parameter domain.

w r = m M -m, ( 3
The distances to the true interface (i.e., eq. 3.11) are comparable for the two types of inversions: for the deterministic inversion the average distance is 3.9 m, while the distance for the different empirical-Bayes-within-Gibbs chains are 3.6 ±0.1 m, 4.5 ±0.2 m and 4.1 ±0.1 m. The error in the MCMC inversion is almost entirely due to the poor performance in the first 50-100 m of the profile. Note that the problematic zones are on

Field data

As already discussed above, we filter the ERT dataset from the Calhoun Critical Zone Observatory (CZO) in South Carolina (USA). Indeed, when considering all 84 electrodes, the chains were unable to converge to the target data misfit. We attribute this to significant small-scale heterogeneity in the near-surface below our model discretization size (see Fig. 3.1a). We could have refined the inversion grid, but we chose instead to decrease the data sensitivity to small-scale near-surface variability by not considering neighboring electrodes and therefore, removing the shallowest pseudo depth level from the data. This choice is motivated by our focus on the bedrock interface and that a model refinement would have led to unnecessarily high computational times.

n agreement with the synthetic test case, we sampled the electrical resistivity, model structure and the interface roughness weights from log-uniform distributions with the structural constraints expressed in Eq. (3.4). Here, the prior range of the hyper parameters are defined broadly enough to avoid boundary effects (see [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF]. To establish the prior range on resistivities, we first determined a global resistivity range (i.e., over the whole subsurface), which was obtained from the deterministic inversion results when considering the full dataset (i.e., considering the full 84-sensor dataset). Our target interface is the one between regolith (weathered/fractured bedrock) and unweathered bedrock, which is clearly seen at 20-40 m depth in seismic and resistivity data presented by St. Clair et al. (2015). In those results, unweathered bedrock occurs where seismic velocity increases to more than 4 km/s and resistivity increases to more than 4000 Ω m. We had to run different tests to establish an upper limit for the CZ resistivity range and a lower limit for the bedrock. [2000,30000] Ωm, were chosen to ensure the convergence to the target data misfit. The physical parameter and interface updates, as well as the priors and proposal pdfs for the hyper-parameters, are implemented in the same way as for the After MCMC inversion, we find that the inferred model structure for the CZ (Fig. 3.10a), the bedrock (Fig. 3.10b) and the interface (Fig. 3.10c) are overall similar for the three chains, as are the hyper-parameters (Figs. 3.10d-f). In agreement with the synthetic results, the acceptance rate for the interface updates (Fig. 3.10g) is lower (average value of 6 %) than for the resistivity field updates (Fig. 3.10h; average value of 17 %). According to the Geweke diagnostic, the burn-in is 5 × 10 5 , 4.5 × 10 5 and 3.5 × 10 5 iterations for the first, second and third chain, respectively. As for the synthetic test case, we could declare formal convergence for only a percentage of the model parameters: 14 % for the CZ and 23 % for the bedrock. This implies that even if we obtain a set of models that are able to explain the observations well, the chains do not fully explore the posterior distribution.

Inversion results in terms of mean, normalized standard deviation of resistivity and bedrock probability map of the "posterior" model realizations are shown in Fig. 3.11. The results show a very similar behavior among the different chains, especially in the upper 40 meters (where the resistivity data are able to constrain the model parameters). This is also seen in the vertical resistivity profiles in Figs. 3.11(l)-(o), where the consistency in the inferred interface locations is manifested by the jumps in resistivity. For comparison purposes, we also run the empirical Bayes framework with three chains for which we did not consider the subdivision of the subsurface. That is, we infer one log-resistivity field using a prior pdf that spans the range of the previously defined sub-domains. Similarly to previous results, we start to sample the posterior probability distribution after 10 5 iterations and we are again unable to fully sample the posterior distribution. The corresponding approximate posterior results in terms of mean, normalized standard deviation and vertical resistivity profiles are represented in Fig. 3.12. Finally, we compared these results to the ones obtained by smoothness-constrained deterministic inversion for the same data misfit (Fig. 3.13c). For both inversions we inferred an interface using the vertical gradient method (i.e., the black lines in Fig. 3.13c and d).

Discussion

We have presented the first inversion approach that uses ERT data to probabilistically infer interface properties in the presence of heterogeneous sub-domains. The results are promising, but there are also certain limitations. Initially, we attempted to address the inverse problem using the more general formulation presented by de Pasquale & Linde ( 2017), but we were unable (despite significant testing) to obtain stable results. We attribute this to the wider prior ranges in the present work compared to the successful case-studies considered by de Pasquale & Linde (2017). Both Rosas-Carbajal et al. ( 2014) and de Pasquale & Linde (2017) have highlighted that the type of structure-constrained inversion routine implemented herein tends to favor model realizations with too little complexity (i.e., too smooth models), especially where the data are weakly sensitive.

The empirical-Bayes-within-Gibbs method successfully sampled the modeled interface in the synthetic example with homogeneous properties (Fig. 3.3). When we consider the synthetic test case with heterogeneous sub-domains (Fig. 3.4), however, the inversion results highlight that the target interface is not part of the sampled interface locations along the first 50 m and around 370 m of the profile length (see Figs. 3.7a,d,g and 3.8a).

The true model in Fig. 3.4(a) presents a strong resistor (40 m depth and 50 m along the profile line) which is situated below a conductive area. At 370 m, the true model has a high resistivity bedrock at the bottom, which is overlain by a conductive deep-CZ anomaly and a higher resistivity feature at the surface. Due to the inherent equivalence problem of DC resistivity data [START_REF] Koefoed | Geosounding Principles I, Resistivity Sounding Measurements[END_REF], the resulting forward response can, thus, also be explained by the posterior realizations we obtained (i.e., a large area with intermediate resistivities for both cases). Such results are favored in the empirical-Bayes inversion because they are less complex (in terms of the variability within the two subdomains) than the underlying true model. These inconsistencies persist even when using homogeneous starting models. Moreover, Fig. 3.8(c) shows that the interface inferred using the maximum gradient method on the deterministic inversion model is shallower than the target one in the beginning of the profile. This inefficiency of the different inversion routines, stresses the inability of the ERT data to resolve this area. Nevertheless, the model results are satisfactory in the central part of the profile and remain an advancement compared with smoothness-constrained deterministic inversion results (see Fig. 3.8). Perhaps the most dramatic improvement is found in the bedrock domain, whose resistivity values are always severely under-estimated by the deterministic inversion.

Even if the posterior model realizations are similar across chains, we are unable to declare convergence of the MCMC chains for all the model parameters (?); this implies that each individual chain has not sufficiently sampled the posterior pdf. Consequently, the spread of the parameters is likely larger than those inferred in each chain. By combining the sampled posterior realizations of the three chains, we partially reduce this issue. Inversion for sharp boundaries between subsurface physical properties is hard as it makes the likelihood surface highly irregular (Tavassoli et al., 2005), especially for high-dimensional problems and large data sets with high signal-to-noise ratios. More advanced MCMC algorithms (e.g., parallel tempering) could improve the situation [START_REF] Laloy | Probabilistic inference of multi-Gaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction[END_REF]. We have found that scaling of the model proposal updates using model sensitivities improve exploration of the posterior distribution. Adaptive Metropolis (Haario et al., 2001) and its variants (e.g., [START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing[END_REF], in which an appropriate proposal distribution is determined for each model parameter might further improve exploration. We also suspect that some of the problems encountered in this study are further enhanced by the fact that ERT data have no (or very limited) inherent sensitivity to the depth and thickness of subsurface layers [START_REF] Parker | The inverse problem of resistivity sounding[END_REF][START_REF] Oldenburg | Estimating depth of investigation in dc resistivity and IP surveys[END_REF].

Another possibility would be to address this problem with a transdimensional inversion algorithm [START_REF] Sambridge | Transdimensional inverse problems, model comparison and the evidence[END_REF][START_REF] Belhadj | New parameterizations for Bayesian seismic tomography[END_REF], in which the complexity of the interface geometry and the physical property fields are determined by treating the model dimension of sub-domains and the number of interface locations as unknowns.

The results obtained for the Calhoun field data are summarized in Fig. 3.13, where we show the models obtained for the different inversion routines considered. Regardless of the inversion approach, we find that the upper part of the domain (first 10-20 meters of depth) is well defined. However, the inferred bedrock interface and bedrock properties are highly dependent on the inversion approach. As expected, we find that the smoothness-constrained deterministic inversion result displays the least features (Fig. ] Ω m. The interfaces we inferred through the maximum vertical gradient method is represented with a black line in Fig. 3.13(c). When considering the subsurface as a single domain, the probabilistic inversion provides large fluctuations in resistivity values (Fig. 3.13b) including geologically-unrealistic conductors at depth. The inferred interface (black line in Fig. 3.13b) shows unrealistic depth variations especially at the far left end of the domain. In contrast, the proposed empirical-Bayes-within-Gibbs algorithm ensures that no conductors are found at depth (Fig. 3.13a), while it still allows for a significant variability within the CZ and bedrock sub-domains. Moreover, the interface is part of the model parameterization and is, hence, clearly described in terms of its depth, geometry and uncertainties. We note that it is only this inversion that clearly identify a deepening of the inferred bedrock interface at the two topographic lows. Note further that the small-scale irregularity of the inferred interface in Fig. 3.13(a) is a consequence of the mesh discretization. It could be removed in a post-processing step by smoothing the interface over a length scale corresponding to twice the average node separation or by using a finer mesh. The vertical resistivity profiles represented in Figs. 3.13(d)-(g) underline the different characteristics of the inversion results: the deterministic inversion results are generally smoother and they present a smaller resistivity range, while the probabilistic inversion with single subsurface domain show larger fluctuations. Finally, the empirical-Bayes-within-Gibbs vertical profiles clearly display the resistivity jump at the interface location and the widest resistivity range.

Our results are now compared with those of St. Clair et al. (2015) regarding the depth and shape of the regolith-bedrock transition. Those authors located the transition between CZ and bedrock at the site between 10 and 40 meters depth, with a shape that mirrors topography (deep under ridges and shallow under valleys). This pattern, which has been attributed to the effect of compressive tectonic stress on fracture opening (St. Clair et al., 2015;[START_REF] Moon | A model of three-dimensional topographic stresses with implications for bedrock fractures, surface processes, and landscape evolution[END_REF] is most apparent in seismic velocity (Fig. 4 of St. Clair et al. 2015). Our results indicate similar depths, but rather interface patterns that follows the topography (Figs. 3.11 a, d, and g). For instance, all the inversion results suggest the presence of a low-resistivity zone in the valley at 300 m along the profile line (m ≈ 2000Ω m). Extensive testing with many alternative prior resistivity ranges and starting models illustrate that such a zone is needed if the data are to be fitted with a realistic-looking model. A possible explanation for this subsurface resistivity pattern is the presence of deep fracture zones or faults underlying the stream valleys. To better understand this difference in behavior, we suggest that joint inversion of the seismic and ERT data would be of great interest.

Conclusions and Outlook

Determining subsurface interface geometries together with reliable uncertainty quantification is important in various Earth science settings. This objective is generally not achievable when using deterministic smoothness-constrained inversions or when inferring for interface locations while ignoring subsurface heterogeneities within layers. In this work, we introduce an empirical-Bayes-within-Gibbs MCMC inversion algorithm that explicitly parameterize and infer both interface geometry and spatial heterogeneity of physical properties. Our synthetic and field-based results consider 2D surface-based ERT data aiming at inferring the regolith-bedrock interface, but the methodology is general and, since it is independent of the forward solver, it can be extended to other observation types (e.g., magnetotellurics or seismic refraction), 3D datasets and parametrization choices. For a synthetic heterogeneous test case, we find that the interface location is well-resolved in the central part of the profile, but less so on the sides where the data are less constraining. The introduction of the interface in the inversion leads to a dramatic improvement in the estimations of the bedrock properties when compared with smoothness-constrained deterministic inversions and MCMC-based inversions without an explicit interface. Nevertheless, longer MCMC chains or adaptations of more advanced sampling methods (e.g., parallel tempering) is needed to fully explore the posterior distribution. For the field example at the Calhoun Critical Zone Observatory, our new method suggested a significant deepening of the bedrock geometry coinciding with topographic lows. The time needed for each iteration is mainly determined by the forward operator (in the case of 2D ERT, this time scale linearly with the number of current electrodes) and the number of nodes in the unstructured forward mesh, which affects the interface update performance (also in this case, the time scales linearly with the number of nodes). A logical extension of this work is to perform joint inversion of ERT and seismic refraction data. This would help to better resolve the regolith-bedrock interface and it would overcome some of the equivalence issues that are inherent with ERT datasets. The model coupling could be achieved by considering a common interface geometry and possibly also parameter correlation of the physical property fields (e.g., to enforce that highly resistive bedrock is likely to correspond to regions of high seismic velocity).

Abstract

The depth and geometry of the bedrock/regolith interface influence many near-surface processes and they can be imaged using non-invasive geophysical surveys that are sensitive to discontinuities in physical properties. For instance, refraction seismic applications rely on the strong correlation between seismic velocity and porosity, while DC resistivity data are sensitive to electrical conductivity, which in turn depends on porosity and clay content. Combining these two data types is beneficial since their spatial resolution patterns and sensitivities to relevant target properties differ. Here, we propose a probabilistic joint inversion framework that uses DC resistivity and first-arrival seismic travel times to infer a common interface geometry separating two heterogeneous sub-domains. A synthetic test case demonstrates that the method is more accurate in sampling the target interface than inversion results obtained using single datasets.

We then applied the joint inversion method to field data from the Calhoun Critical Zone Observatory in South Carolina, where it resolves apparent contradictions in the individual inversion results. It produces a bedrock topography that mirrors surface topography and offers clear indications of a more fractured bedrock below the surface topography lows.

Introduction

The depth to the bedrock and its topography influence many surface and near-surface processes [START_REF] Rempe | A bottom-up control on fresh-bedrock topography under landscapes[END_REF], such as, water runoff to channels [START_REF] Onda | The role of subsurface water flow paths on hillslope hydrogeological processes, landslides and landform development in steep mountains of Japan[END_REF] and water chemistry [START_REF] Anderson | Weathering profiles, mass-balance analysis, and rates of solute loss: Linkages between weathering and erosion in a small, steep catchment[END_REF]. Furthermore, the bedrock/regolith interface is a key boundary when studying erosion rates and certain atmospheric processes related to plant water uptake and transpiration [START_REF] Jones | Water-holding characteristics of weathered granitic rock in chaparral and forest ecosystems[END_REF]. Its characteristics are crucial for modeling soil thickness, which depends on the balance between production at depth (i.e., where weathering of bedrock takes place) and erosion at the surface [START_REF] Heimsath | The soil production function and landscape equilibrium[END_REF]; an interplay which ultimately determines the critical zone (CZ) thickness (Brantley et al., 2011). Moreover, this boundary is a key controlling factor for landscape evolution and landslide investigations [START_REF] Lanni | Simulated effect of soil depth and bedrock topography on near surface hydrologic response and slope stability[END_REF][START_REF] Reneau | Size and location of colluvial landslides in a steep forested landscape[END_REF]. For these reasons, its imaging is beneficial for a wide range of application areas. Hence, geophysics has become a common tool in geomorphological research [START_REF] Schrott | Geophysical applications in geomorphology[END_REF], because it allows for non-invasive and relatively fast surveys that are sensitive to discontinuities in physical properties [START_REF] Saas | Bedrock detection and talus thickness assessment in the European Alps using geophysical methods[END_REF]. Moreover, geophysical measurements provide complementary information about subsurface heterogeneities that is complementary to those offered by remote sensing, geomorphological or geological mapping (Olona etal., 2010;[START_REF] Tye | Using integrated near-surface geophysical surveys to aid mapping and interpretation of geology in an alluvial landscape within a 3D soil-geology framework[END_REF][START_REF] Hirsch | A Comparison of Electrical Resistivity, Ground Penetrating Radar and Seismic Refraction Results at a River Terrace Site[END_REF].

One geophysical technique that is traditionally used to infer interface geometries is the seismic refraction method, which uses the first arrival travel times of critically refracted body waves to reconstruct P-wave velocity models of the subsurface (Zhuo, 2014). Given the strong correlation between seismic velocity variations and contrasts in density and porosity, the seismic refraction method has been used to detect subsurface interfaces; for instance, the thickness of weathered layers or fractured zones [START_REF] Lee | Seismic refraction surveys for predicting the intensity and depth of weathering and fracturing in granitic masses[END_REF].

Another geophysical method that is widely used in this respect is electrical resistivity tomography (ERT; [START_REF] Chambers | Bedrock detection beneath river terrace deposits using three-dimensional electrical resistivity tomography[END_REF][START_REF] Hsu | Bedrock detection using 2D electrical resistivity imaging along the Peikang River, central Taiwan[END_REF]. ERT provides comparatively high-resolution images of the shallow subsurface (down to depths of tens to a few hundreds of meters) with changes in electrical resistivity being strongly related to changes in water and clay content [START_REF] Dahlin | 2D resistivity surveying for environmental and engineering applications[END_REF]. Their difference in spatial resolution and sensitivity to target properties (e.g., water content) has motivated the combined use of seismic refraction and ERT in numerous studies [START_REF] Berge | Joint inversion of geophysical data for site characterization and restoration monitoring[END_REF][START_REF] Linder | Zonal cooperative inversion of crosshole P-wave, S-wave and georadar traveltime data sets[END_REF][START_REF] Hellman | Structurally coupled inversion of ERT and refraction seismic data combined with cluster-based model integration[END_REF][START_REF] Juhojuntti | Joint inversion of seismic refraction and resistivity data using layered models -Applications to groundwater investigation[END_REF]. It is well established that joint inversions aiming at inferring the two physical fields (electrical resistivity and P-wave velocity) can improve resolution and decrease interpretation ambiguity compared to individual inversions [START_REF] Gallardo | Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints[END_REF].

In this paper, we focus on the seismic refraction and DC resistivity datasets that were used by St. Clair et al. (2015) to infer the geometry of the interface between the critical zone (CZ) and the intact bedrock at the Calhoun Critical Zone Observatory in South Carolina. In their work, the underlying motivation was to better understand the influence of surface topography on bedrock weathering. To do so, St. Clair et al. (2015) inverted the two datasets separately before interpreting them qualitatively. Here, we jointly invert the two datasets using structural constrains that impose a common interface geometry of the two physical property fields (Kozlovskaza et al., 2007;[START_REF] Linde | Joint Inversion in Hydrogeophysics and Near Surface Geophysics[END_REF]. To achieve this, we adapt a probabilistic inversion framework for single-method data aiming at inferring an interface separating two heterogeneous sub-domains [START_REF] De Pasquale | Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data[END_REF]. The changes made to this work include the consideration of two data types with model coupling offered by a common interface. We investigate the added value offered by the joint inversion results compared to those obtained by individual inversions. Based on these results, we seek to complement the findings and interpretations by St. Clair et al. (2015).

Methodology

Probabilistic joint inversion with interface constraints [START_REF] De Pasquale | Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data[END_REF] presented a probabilistic formulation and solution to the inverse problem of using one geophysical data type to infer the location and geometric shape of an interface separating two heterogeneous sub-domains. In this empirical-Bayes-within-Gibbs algorithm, the interface geometry and physical property field updates are performed alternately within a Markov chain Monte Carlo (MCMC) scheme [START_REF] Sambridge | Monte Carlo methods in geophysical inverse problems[END_REF]. Furthermore, the interface geometry and the physical properties of the sub-domains are constrained to favor smooth spatial transitions and pre-defined property bounds. The algorithm by de Pasquale et al. ( 2019) is here adapted to joint inversion in which the interface is considered common for two physical property fields. The subsurface is parameterized by two overlapping domains that are divided by an interface, which in 2D is a line intersecting connected nodes in the common unstructured mesh that is used for forward computations. The joint probabilistic inversion as adapted to joint inversion is summarized below in four steps for the problem of delineating the CZ, while a more detailed treatment of algorithmic details is provided by [START_REF] De Pasquale | Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data[END_REF]. (1) Chain initialization: the initial model realizations of seismic P-wave velocity (m v 0 ) and electrical resistivity (m r 0 ) are obtained from the initial property fields of both the CZ (m r CZ 0 , m vCZ 0 ) and bedrock (m r B 0 , m vB 0 ) domains over the whole model domain (Figs. 4.1a and b) and the initial interface (I 0 , Fig. 4.1c). We initialize the physical properties on a regular grid by randomly drawing values from uncorrelated log-uniform distributions. These fields are linearly interpolated on the unstructured mesh used for the forward computations. The interface is sampled as a set of connected nodes within the unstructured mesh and is used to map the field that is to be used for forward simulations (i.e., CZ above and bedrock below the sampled interface as shown in Figs. 4.1d and e). Once m v 0 and m r 0 are built, we evaluate their likelihoods [START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF]: L(m v 0 , I 0 |d RS ) and L(m r 0 , I 0 |d ERT ), where d RS and d ERT are the first-arrival travel times and apparent resistivities, respectively. Following de [START_REF] De Pasquale | Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data[END_REF], we also evaluate corresponding constraint functions that are used to quantify model parameter variability, c(S v,r 0 ), and interface roughness, c(R 0 ). Their values are used in the inversion to penalize model structure with weights that are inferred using the empirical Bayes approach as described by [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF]. For both constrain types, we rely on a total variation measure.

(2) Interface update: after updating the interface, following a principle of minimum change, both the physical parameter fields (seismic velocity and electrical resistivity) are re-mapped accordingly to the proposed interface. The update is then accepted or rejected with probability:

α I = mi n 1, L(m v cur r , I pr op |d RS )L(m r cur r , I pr op |d ERT )c(R pr op ) L(m v cur r , I cur r |d RS )L(m r cur r , I cur r |d ERT )c(R cur r ) , (4.1)
where the subscripts prop and curr refer, respectively, to the proposed and current state of the Markov chain. (3) Physical property updates: both physical property fields are perturbed for each sub-domain and mapped according to the interface obtained in step (2). The resulting resistivity and velocity field proposals are accepted or rejected independently from each other, with probabilities:

α ERT = mi n 1, L(m r pr op , I cur r |d ERT )c(S r pr op ) L(m r cur r , I cur r |d ERT )c(S r cur r ) , (4.2)
α RS = mi n 1, L(m v pr op , I cur r |d RS )c(S v pr op ) L(m v cur r , I cur r |d RS )c(S v cur r ) . ( 4.3)
The procedure is then repeated starting from step 2. (4) Chain finalization: after "burn-in" (i.e., the number of iterations needed to start sampling proportionally to the posterior distribution; [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF] and when the stationary part of the chain has reached the desired length, the iteration loop over steps 2 and 3 terminates and the MCMC realizations represent our estimate of the posterior distribution of the investigated subsurface.

A schematic flowchart of the algorithm is represented in Fig. 4.2.

Figure 4.2:

Schematic flow chart of the proposed probabilistic joint inversion algorithm with a common interface. First, we initialize the electrical resistivity (m r 0 ) and seismic P-wave velocity (m v 0 ) models (as described in Fig. 4.1), evaluate their likelihoods (L) and the corresponding structural constraint functions. We then sequentially evaluate proposed interface and physical property updates in an iterative scheme until the chosen length of the chain is reached. Here, u ∼ U (0, 1) represents a randomly drawn number with uniform probability between 0 and 1, while the acceptance probabilities α I , α E RT and α RS are described in eqs. 4.1, 4.2 and 4.3, respectively. For more details, please refer to de [START_REF] De Pasquale | Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data[END_REF].

DC resistivity is a geophysical method that is used to image the subsurface soil and rock resistivity by observations made with electrodes at the Earth's surface or in boreholes [START_REF] Binley | DC Resistivity and Induced Polarization Methods[END_REF]. Voltage measurements are made between pairs of potential electrodes in response to current flowing between two current electrodes. The geometry of the current injection and potential electrode pairs are varied; typical set-ups involve many tenths (or hundreds) of electrodes and a few hundreds or thousands of data. These data are then used in an inversion process to infer the spatial distribution of electrical resistivity in the subsurface; a process that is often referred to as ERT.

We use the Boundless Electrical Resistivity Tomography (BERT) library [START_REF] Rücker | Three-dimensional modelling and inversion of dc resistivity data incorporating topography-I[END_REF][START_REF] Günther | Three-dimensional modelling and inversion of dc resistivity data incorporating topography-II. Inversion[END_REF] to simulate the DC resistivity forward response needed to compute the likelihood of each model proposal within the Markov chain. Using BERT, we can compute the 2.5D forward response by relying on a finite element scheme implemented on unstructured meshes (Si, 2015). Unstructured meshes allow for efficient local refinement and is suitable to account for surface topography and internal boundaries.

Seismic refraction tomography often relies on first-arrival travel times corresponding to direct and critically-refracted seismic waves to reconstruct seismic P-wave (i.e., compressional wave) velocity models. When seismic waves impinges on velocity boundaries, they undertake a change in their direction of propagation. At a critical angle that depends on the velocity contrast, head waves are created that moves along the interface at the speed of the faster lower-lying layer velocity and emits refracted waves. The time at which these waves arrive at receiver positions (i.e., first-arrival travel times) are the main observations used in seismic refraction surveys. Besides the density of the rock matrix (related to mineralogical composition and texture), seismic velocities are affected by porosity and pore-fluids, as well as confining stress, pressure and temperature [START_REF] Schmitt | Geophysical Properties of the Near Surface Earth: Seismic Properties[END_REF].

IIn order to simulate the first-arrival travel times, we used the physics refraction class of pyGIMLI (Geophysical Inversion and Modelling Library in Python; [START_REF] Rücker | pyGIMLi: An open-source library for modelling and inversion in geophysics[END_REF]. Specifically, the forward operator used is based on the calculation of the fastest paths from sources to receivers along the elements in the unstructured mesh [START_REF] Moser | Shortest path calculation of seismic rays[END_REF][START_REF] Dijkstra | A note on two problems in connection with graphs[END_REF]. This algorithm by Dijkstra is very popular, but also imperfect, since the traveled distance is always overestimated given that rays can only move along the discretized elements. It is used herein because of its speed, but we will see later that the related modeling errors have adverse effects on the joint inversion results.

Petrophysical relationships

In most rocks and soils, electrical current is due to electrolytic conduction of ions. As a consequence, the main factors affecting subsurface resistivity are water content, salinity, porosity and the connectivity of the water phase [START_REF] Lesmes | Relationships between the Electrical and Hydrogeological Properties of Rocks and Soils[END_REF]. Another important factor affecting electrical resistivity is clay type and content through its contribution to enhanced conduction within the electrical double layer at the mineral-water interface. Among many petrophysical relationships, we consider below the relationship proposed by [START_REF] Linde | Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data[END_REF], which is based on volume-averaging in the high-salinity limit (the pore contribution to overall conduction is much larger than the surface contribution), to highlight how electrical conductivity (the inverse of electrical resistivity) is affected by properties and state variables of interest:

σ e f f = 1 F S n w σ w + (F -1)σ s . (4.4)
Here, σ e f f (S/m) is the effective conductivity of the matrix-water system at the scale of a Representative Elementary Volume (REV), σ s the surface conductivity, σ w the pore water conductivity, F = φ -m the electrical formation factor which depends on the porosity, φ, and on the cementation exponent m, S w is the water saturation and n is the saturation exponent, which is related to the tortuosity of the water phase.

To quantitatively interpret seismic properties in fractured rock it is common to rely on an equivalent medium representation. Here, the multiphase (minerals, water, air) properties and their connectivities are replaced locally (at the scale of a REV) by an upscaled homogeneous medium with the same macroscopic properties [START_REF] Liu | Seismic Fracture Characterization. Concepts and Practical Applications[END_REF]. In this field of rock physics [START_REF] Mavko | Rock physics handbook: Tools for seismic analysis in porous media[END_REF], there are many formulations with different underlying assumptions. They range from macroscopic models of heterogeneous rocks that do not account for microscopic heterogeneities (i.e., pore shape, connectivity and fluid distribution) [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low frequency range and II. higher frequency range[END_REF][START_REF] Gassmann | Über die Elastizität poröser medien. Vierteljahresschrift der naturforschenden[END_REF], to inclusion-based models which upscale microscopic features to describe the macroscopic elastic behavior of the rock [START_REF] Xu | Modelling the effects of fluid communication on velocities in anisotropic porous rocks[END_REF][START_REF] Hudson | Wave speeds and attenuation of elastic waves in material containing cracks[END_REF][START_REF] Chapman | Frequency dependent anisotropy due to meso-scale fractures in the presence of equant porosity[END_REF]. A rough indication of the impact of porosity in saturated rocks is given by Wyllie's mixing law [START_REF] Wyllie | Elastic wave velocities in heterogeneous and porous media[END_REF][START_REF] Wyllie | AN experimental investigation of factors affecting elastic wave velocities in porous media[END_REF]):

1

V e f f = φ V f + 1 -φ V m , ( 4.5)
where V e f f , V f and V m are the P-wave velocities of the saturated rock, the pore fluids and the rock matrix (mineral matrix), respectively. All common rock physics models and observations agree that seismic velocities in fractured rocks is lower than in intact rocks, with the contrast depending on fracture porosity, connectivity and filling material.

For instance, the P-wave velocity in water is five times greater than in air, therefore, water saturated porous/fractured rocks presents a higher elastic wave velocity than unsaturated rocks [START_REF] Barton | Rock Quality, Seismic Velocity, Attenuation and Anisotropy[END_REF][START_REF] Nur | The effect of saturation on velocity in low porosity rocks[END_REF].

Results

Our field-based case-study considers the datasets collected at the Calhoun Critical Zone Observatory in South Carolina and interpreted by St. Clair et al. (2015). Here, DC resistivity and seismic refraction surveys were used to investigate the influence of surface topography on bedrock geometry. The DC resistivity survey was acquired using 56 electrodes spaced 5 m apart with a dipole-dipole electrode configuration. The 400 m long profile was obtained by using one roll-along in which 50 % of the electrodes stayed in place, which resulted in a total of 84 different electrode positions. Since our interest is in the deeper CZ structure, we could save computational time by only considering every second electrode position and removing configurations with a maximal electrode separation of 30 m. This resulted in 42 electrodes and 645 measurements being used.

For the seismic refraction survey, 96 geophones were used with a spacing of 2.5 m and the shots spaced 10 m apart. This configuration of 237.5 m in length was repeated once after moving the set-up and complemented with off-line shots to link the two acquisitions, resulting in 192 geophone positions. We manually picked the first arrival travel times on each trace after processing the data with a zero-phase Butterworth filter (15 Hz 8dB/octave low cut and 150 Hz 24 dB/octave high cut) to improve the signalto-noise ratio. Particularly, the signal was weak at far offsets and it was particularly noisy at the end of the profile. Travel time curves were plotted for quality control purposes, whereby erroneous picks were ascertained and deleted from the data through comparison with neighboring travel-time curves. This process resulted in 3896 travel times that were kept for inversion. For the ERT dataset, we used a noise description based on an uncorrelated relative Gaussian error of 3.7 %, while for the refraction seismic observations we considered an absolute error of 0.003 seconds. These choices were made to obtain similar data misfits as for the inversion results presented by St. Clair et al. (2015), thereby, facilitating comparisons.

Synthetic test case

To demonstrate the probabilistic joint inversion with interface constraints, we first consider a synthetic test case in which DC resistivity and seismic refraction surveys are simulated for a test model. The resulting data are noise-contaminated and subsequently used to jointly invert for the interface between regolith and intact bedrock. In the synthetic test case, we consider identical survey layouts and the same surface topography as for the Calhoun CZO. Moreover, the physical properties are assumed to have a loguniform prior with the same prior ranges as later used for the Calhoun field data 4.3(d) represents its interpolation on the mesh used for the forward computation. Fig. 4.3(e) shows the ERT inversion model result obtained by a l 1 -norm mimicking smoothness-constrained deterministic inversion, while Fig. 4.3(f ) shows the corresponding deterministic inversion model result for the refraction seismic observations. In Figs. 4.3(e) and (f) we indicate the interface obtained by using the maximum vertical gradient method on the deterministic inversion results (following the procedure described by [START_REF] Chambers | Bedrock detection beneath river terrace deposits using three-dimensional electrical resistivity tomography[END_REF]. Moreover, the deterministic inversion model results are depicted with a transparency that is inversely proportional to the data coverage. This highlights that the DC resistivity data have virtually no coverage after 400 m along the profile, simply because no electrodes were located in this area. Also, the coverage of the seismic refraction measurements deteriorates significantly after 400 m. Therefore, for the analysis and interpretation of the inversion results, we focus only on the results obtained in the first 400 m. Finally, in Figs. 4.3(g-h) we show a measure of discrepancy between the target model and the inversion results:

w r = m r t r ue -m r , (4.6)
w v = m v t r ue -m v ,
where m r t r ue and m v t r ue are the vectors of the resistivity and velocity values of the true model (Fig. 4.3c for the resistivity and Fig. 4.3d for the velocity). The mean absolute value of the model discrepancy is 0.28 for the resistivity and 0.12 for the velocity inversion model results. In both deterministic inversions, the data were fitted to the same error levels as those used for the subsequent probabilistic inversions. In fact, the assumed errors in the inversions were increased by 1.1 times with respect to the standard deviation used to contaminate the DC resistivity observations and by 1.7 times for the seismic refraction observations. This inflation of the noise level was needed to partially account for the modeling error related to the change between the unstructured mesh used for creating the synthetic observations and the one used for the forward computations within the MCMC chains. The factors were determined by comparing the residuals between the forward simulations for the two meshes. To assess the added value of the proposed probabilistic joint inversion algorithm, we also run the empirical-Bayes-within-Gibbs algorithm by de [START_REF] De Pasquale | Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data[END_REF], that is, considering the DC resistivity and refraction seismic datasets individually. For each case (i.e., two individual and one joint inversion), we run three independent MCMC chains for 5 × 10 5 iterations with each iteration taking ≈ 0.4 seconds when inverting the seismic dataset, ≈ 1.7 seconds when inverting the DC resistivity dataset and ≈ 2 seconds for the joint inversion. We start all nine MCMC chains with random uncorrelated initial models and use Geweke analysis [START_REF] Geweke | Evaluating the Accuracy of Sampling-Based Approaches to the Calculation of Posterior Moments[END_REF] to asses the burn-in period of the chain (the convergence statistics for the different inversion schemes are summarized in Table 4.1). Considering the last 2x10 5 iterations, the potential scale reduction factor ( R; [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF]) is below 1.2 for the majority of the model parameters of the velocity fields for both the bedrock (100 % of the parameters for individual and 98 % for joint inversion ) and the CZ sub-domains (69 % for seismic refraction and 61 % for joint inversion). The resistivity field parameters have R < 1.2 for the majority of the bedrock parameters (70 % for individual and 68 % for joint inversion) but not for the CZ sub-domains (20 % for individual and 24 % for joint inversion). Since this threshold value is typically used to assess if the MCMC chains have sampled the posterior sufficiently, we find that we did not sufficiently sample the posterior pdf, in particular, the CZ resistivity fields. Table 4.1: Convergence statistics for the probabilistic inversion of synthetic datasets. ERT stands for probabilistic inversion of DC resistivity data with the empirical-Bayes-within-Gibbs approach, RS refers to empirical-Bayes-within-Gibbs inversion of the seismic refraction data and Joint refers to the probabilistic joint inversion of the two datasets with interface constraints. The sub-scripts (1,2,3) indicate the three different chains used for each inversion routine. AR stands for the averaged acceptance rate of the chains after burn-in and the sub-scripts I, v and r refer respectively to the interface, velocity and resistivity fields. Finally the burn-in columns enumerate the number of iterations needed for each chain to converge to the target posterior distribution [START_REF] Geweke | Evaluating the Accuracy of Sampling-Based Approaches to the Calculation of Posterior Moments[END_REF]. In Fig. 4.4 we display the last posterior model realization for each of the nine MCMC chains, while in Fig. 4.5 we plot the corresponding vertical resistivity and velocity profiles at 50 m, 150 m, 250 m and 350 m along the profile length, together with target model values. For each inversion case, the three individual MCMC chains show similar posterior model realizations and vertical profiles. Here, the consistency in the inferred interface locations is manifested by similar depths and magnitudes of the jumps in resistivity and/or velocity. Fig. 4.6 shows the mean, normalized standard deviation, probability maps of the inferred interface location and model discrepancies (eq. 4.6).

Chain AR I [%] AR v [%] AR r [%] Burn-in v [n°of iterations] Burn-in r [n°of
These results indicate that the probabilistic inversion of a single dataset is able to sample interfaces that are partially coherent with the target one. In the case of individual inversion of the DC resistivity dataset (Fig. 4.6a), the chains are found to sample the underlying interface well in the central area of the investigated domain (i.e., between 100 m and 300 m) but not on the sides (i.e., where there is less coverage). The probabilistic inversion of the seismic refraction dataset (Fig. 4.6i) samples the proper bedrock topography, but it places the interface slightly deeper than the true one. Finally the joint inversion model results (Figs. 4.6d and g) combines the resolution ability of the two data types and obtain an interface that is more consistent in terms of shape and depth. This is also seen by studying the discrepancies between the target models and the mean of the posterior realizations. In case of individual inversion of ERT data,Fig. 4.6(n) shows a thick area of over-estimated resistivity (hence the negative relative error) on the left and right side of the domain, reflecting the too shallow interface that was sampled by the MCMC chains. The mean absolute value of model discrepancy is in this case 0.24. For the individual refraction seismic inversion, Fig. 4.6(p) shows a layer of underestimated velocity (hence the positive relative error) along the whole profile, mimicking the target interface shape and expressing a deeper interface than the target one. This discrepancy leads to a mean absolute model discrepancy of 0.074. Finally, for both resistivity and velocity, the probabilistic joint inversion mitigates these issues. The improvement is particularly evident in the resistivity model discrepancy (Fig. 4.6o), where the area of over-estimation of the inferred field is significantly reduced. An improvement is also seen at 80 m and 300 m along the profile length (i.e., the valleys) for which the velocity model discrepancy (Fig. 4.6q) reduces significantly. The resulting mean absolute model discrepancies are 0.17 for resistivity and 0.066 for velocity.

Calhoun CZO

We now turn our attention to the results obtained when considering the field data from the Calhoun CZO. In Fig. 4.8, we show the (a) DC resistivity and (b) seismic refraction data together with the resulting (c-d) l 1 -norm mimicking smoothness-constrained deterministic inversion results and the corresponding interface obtained using the maximum vertical gradient method. The deterministic inversion model results are plotted with transparency that is inversely proportional to the data coverage, which shows the inability of both datasets to resolve subsurface properties after 400 m along the profile. For this reason, the model results shown hereafter are cut at 400 m.

For each individual or joint inversion approach, we run three independent MCMC chains for 5 × 10 5 iterations. The initial models are again randomly sampled from log-uniform uncorrelated property fields. The convergence statistics of the different probabilistic inversion schemes are summarized in Table 4.2. Using the last 10 5 MCMC iterations, we find that the potential scale reduction factor is only below 1.2 for most model parameters describing the velocity field of the bedrock sub-domain (99 % of the parameters for both individual and joint inversions). However, the CZ velocity fields have only 39 % (individual inversion) and 40 % (joint inversion) of the parameters with R < 1.2. For the resistivity fields, only 11-12 % of the CZ and 36-37 % of the bedrock parameters have R < 1.2 regardless of if a joint or individual inversion is carried out. This implies that only an incomplete sampling of the posterior distributions were achieved for most of the model parameters considered. This suggests that model uncertainty is underestimated when using the results of one MCMC chain only. However, this problem is partly mitigated by merging the results obtained by the three MCMC chains. Table 4.2: Convergence statistics for the probabilistic inversion of the Calhoun CZO data. ERT stands for probabilistic inversion of DC resistivity data with the empirical-Bayes-within-Gibbs approach, RS refers to empirical-Bayes-within-Gibbs inversion of the seismic refraction observations and Joint refers to the probabilistic joint inversion of the two datasets with interface constraints. The sub-scripts (1,2,3) indicate the three different chains used for each inversion routine. AR stands for the averaged acceptance rate of the chains after the burn-in and the sub-scripts I, v and r refer respectively to the interface, velocity and resistivity fields. Finally the burn-in columns enumerate the number of iterations needed for each chain to converge to the target posterior distribution [START_REF] Geweke | Evaluating the Accuracy of Sampling-Based Approaches to the Calculation of Posterior Moments[END_REF]. Figs. 4.9b,d,f with Figs. 4.9a,c,e and Figs. 4.9h,l,n with Figs. 4.9g,i,m), we see that the bedrock topography is mainly constrained by the seismic refraction observations. The DC resistivity data has the strongest sensitivity to heterogeneity within the bedrock layer. This is seen in the well-defined low resistivity zone of the bedrock around 300 m while the velocity fields of the bedrock are more homogeneous. Finally, around 80 m along the profile, all the posterior velocity models obtained by joint inversion show a slightly deeper bedrock interface than the ones inferred from the inversion of the seismic observations alone. In Fig. 4.11 we show the inversion results in terms of mean, normalized standard deviation and probability maps of the interface locations for the three MCMC chains for each inversion case. In Fig. 4.12, we plot the vertical resistivity and vertical profiles obtained from l 1 -norm mimicking smoothness-constrained deterministic inversion, as well as the individual and joint inversions with interface constraints.

Chain AR I [%] AR v [%] AR r [%] Burn-in v [n°of iterations] Burn-in r [n°of

Discussion

As an extension of the empirical-Bayes-within-Gibbs algorithm by de Pasquale et al.

(2019), we have introduced a probabilistic joint inversion method to infer an interface separating two heterogeneous sub-domains. Our synthetic example demonstrates that the joint inversion provides more accurate interface estimates than those based on DC resistivity or seismic refraction data alone. Furthermore, a general improvement is obtained with respect to deterministic inversion results in terms of the inferred interface geometry and physical property fields, particularly with respect to bedrock properties. This is seen by studying the model discrepancies with respect to the true model for the deterministic (Fig. 4.3g) and probabilistic (Fig. 4.6n) inversions of the DC resistivity data from which it is evident that the deterministic inversion severely underestimates bedrock resistivity. For the seismic data, the deterministic inversion (Fig. 4.3h) overestimates the velocity in the central area of the bedrock sub-domain (between 100 and 300 meters along the profile) while the errors for the probabilistic case (Fig. 4.6p) are the largest in the vicinity of the inferred interface. The improved ability to locate the interface was quantified by a measure of the distance between the target and the inferred interface (eq. 4.7). For the joint inversion results, this distance is at least 30 % smaller than those obtained by individual deterministic or probabilistic inversions, thereby, highlighting the added value of jointly inverting the two types of geophysical data.

When inverting the synthetic datasets it was necessary to inflate the assumed data errors by 10 % for the ERT data and 70 % for the seismic data. This was needed to account for modeling errors arising from the discrepancies between the forward simulations when using the mesh used to generate the observations (Figs. 4.3a and b) and the one used for forward simulations within the MCMC runs (Figs. 4.3c and d). This error inflation was calculated by computing the weighted root mean squared errors between the datasets and the model predictions obtained from the models depicted in Figs. 4.3(c) and (d).

Indeed, different meshes have different possible interface geometries that the MCMC chains can sample, as they are bound to follow the pre-defined node connections. The true model has a rather smooth interface and the mesh conforms to this surface, while the target interface defined by the mesh used in the MCMC forward simulations is more irregular. The travel time path calculations used are restricted to ray-paths that follow the elements connecting the nodes in the network. This introduces errors in the ray geometry and it will overestimate the traveltimes [START_REF] Moser | Shortest path calculation of seismic rays[END_REF]. The algorithm used is in fact said to be inconsistent with the problem it seeks to approximate since it uses a network solution for a continuous problem [START_REF] Sethian | Fast Marching Methods[END_REF]). An alternative that could improve the accuracy of the forward computation (and reduce the modeling error) is, for instance, to use the fast marching approach which directly approximates the solution of the underlying partial differential equation [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Fomel | A variational formulation of the fast marching eikonal solver[END_REF]. Fig. 4.13 shows posterior model realizations whose model predictions fit the synthetic datasets to the added data errors (i.e., without inflation to account for model errors). The effects are small for the individual probabilistic inversions of DC resistivity (Fig. 4.13a) and seismic (Fig. 4.13b) datasets. However, strong artifacts appear when we jointly invert the two datasets (Figs. 4.13c and d). Indeed, a much more complex interface is needed to accommodate the model errors in both the seismic and DC resistivity data. We expect that a more detailed quantification and inclusion of model errors [START_REF] Hansen | Accounting for imperfect forward modeling in geophysical inverse problems -Exemplified for crosshole tomography[END_REF] would further improve our results. We leave this topic for future research. This problem does not appear as acute for the Calhoun CZO inversions, probably because of the rather conservative choice of the data error models. The results obtained by jointly inverting the Calhoun datasets indicate sharp interfaces at locations that agree well with the more gradual transitions in St. Clair et al. (2015).

In agreement with St. Clair et al. (2015), the transition between the CZ and bedrock has a shape that mirrors surface topography, with bedrock rising almost to the surface below the valleys. Our results clearly indicate two regions of lower resistivities within the bedrock domain that are located below the topographic depressions (i.e., around 80 m and 350 m along the profile; Fig. 4.11d). Extensive testing indicates that it is impossible to fit the ERT data without introducing such lower-resistivity regions. We also note that corresponding features are also visible in Fig. 4a in St. Clair et al. (2015). These lower-resistivity zones are likely to be caused by a higher fracture intensity. In Fig. 4.14, we display the resistivity and velocity values corresponding to the bedrock at the inferred interface between bedrock and regolith. The resistivity values drop below the topographic depressions, while the seismic velocities slightly increase in magnitude.

We attribute the increase in seismic velocity to model errors (as explained above) that is likely to lead to over-estimated velocities that compensate for ray-paths that are too long. Since there is no seismic evidence of mechanical weakening in the bedrock below the topographic lows, we expect that the decrease in resistivity is caused by a small increase in fracture porosity and that these fractures may by clay-filled. Induced polarization data or a borehole would shed more light about the presence of clay-filled fractures and lead to more conclusive findings.

Conclusions

The geometry and depth to the regolith-bedrock interface impact many surface and near-surface processes. de Pasquale et al. (2019) presented a probabilistic inversion framework to infer the posterior distribution of an interface separating two heterogeneous sub-domains. We have extended this framework to joint inversion of DC resistivity and seismic refraction datasets in presence of a common interface with unknown geometry. The synthetic example demonstrates that the interface inferred by joint inversion is significantly more accurate than those obtained by individual deterministic or probabilistic inversions of either DC resistivity or seismic refraction data. When applied to field-data from the Calhoun CZO in South Carolina, we find in agreement with St. Clair et al. (2015) that the bedrock mirrors the surface topography. Furthermore, the sampled posterior resistivity fields suggests that bedrock below the valleys in surface topography has an enhanced fracture porosity. Since there is no evidence of a coincident mechanically-weak zone in seismic P-wave velocity, we suggest that the fractures are clay-filled. In the travel time forward solver used, the ray-paths are restricted to the elements linking the nodes leading to ray-lengths that are always overestimated. This bias creates artifacts in the joint inversion results in terms of an overly variable bedrock interface. This problem can be reduced by assuming too high data errors and we suggest that a more accurate forward solver or a more advanced treatment of model errors would lead to further improvements. Our methodology is applicable to other method combinations, application areas and it could easily be extended to three or more datasets.

Chapter 5

Conclusions

Inferring the position and geometry of unknown subsurface interfaces is a challenging, yet crucial task to investigate numerous phenomena of interest in the Earth sciences. Geophysics data contains useful informations for characterizing the Earth near-surface but they are subject to multiple interpretations. When addressing the inversion problem in a probabilistic framework, the shape and formulation of the prior distribution have strong influence on the posterior model realizations, especially where the observation do not constrain the model parameters. When the subsurface is parameterized via spatial discretization, with a moderately-high (order of hundreds) number of model parameters, the choice of priors reflecting limited knowledge about the subsurface is not trivial. In geophysics it is often suggested to use uncorrelated log-uniform priors, which increases the complexity of the inversion problem and leads to posterior realizations that are not geologically realistic. Moreover, even though standard geophysical investigation methods provide useful information for detecting subsurface boundaries, the inversion routines available do not target at the same time interface locations and geometries and the heterogeneous physical property fields as part of their solution. Also, available methods for detecting interfaces from inversion results (e.g., maximum gradient methods) do not provide uncertainty quantification. In this sense, formulating the inversion problem of detecting unknown subsurface interfaces in the presence of heterogeneous sub-domains in probabilistic terms is a reasonable choice, since the solution is given in the form of a probability distribution and, therefore, is inherently described by its uncertainty. To address these problematics, in this thesis we proposed a probabilistic formulation and solution to the inverse problem of using multiple geophysical datasets to infer interfaces in the presence of heterogeneous sub-domains when little a priori knowledge is available.

In Chapter 2 we proposed a theoretically solid approach to enable prior sampling in terms of a pre-defined global model structure metric within MCMC inversion. The measures of model structure considered are related to standard roughness and damping based on l 1 and l 2 -norms. These are normally used within classical regularized deterministic inversion and are, therefore, easy to define and understand for most researchers in geophysics. We focused on structural uniform priors, in order to give equal weight to all possible levels of model structure and, therefore, to propose an alternative prior formulation to the uncorrelated log-uniform pdf in case of scarce prior knowledge. Nevertheless, the method is general and can be applied to other density functions or model summary statistics. We demonstrated the methodology on synthetic and field-based GPR datasets, showing that the structure-based prior approach is able to adequately sample the chosen prior distribution of model structure and to provide satisfactory posterior model realizations and statistics. The structure-based prior approach is built on solid theory, re-defining the Metropolis-Hastings acceptance ratio in terms of the chosen measure of model structure. Nevertheless, its performance is affected by the numerical implementation: the need to account for the asymmetry of the model structure proposal distributions requires to numerical estimation of the proposal pdf ratios, which is the main challenge and limitation of this algorithm. The approach we proposed is to empirically estimate the two proposal pdfs from a limited number of samples at each step of the MCMC chain. Here, the trade-off between efficiency and precision lies in the number of samples we produce in order to approximate the pdfs. Moreover, when considering inversion problems with moderately high dimension and wide ranges for the model parameter values, exhaustive sampling of the proposal pdfs on the model structure is too computationally demanding. Thus, our approximation scheme cannot be applied, as in the case of synthetic and field studies presented in Chapters 3 and 4. A more efficient and precise approximation of pdf ratios could significantly improve the performance of the structure-based prior inversion. Keeping the same approximation approach, one first improvement (in terms of computational time) could be achieved by parallelizing the generation of samples from the model structure proposal pdfs. Nevertheless, in case of wide ranges for the parameters values, the tails of the distributions are extremely hard to sample.

In Chapter 3 we presented the empirical-Bayes-within-Gibbs algorithm as a probabilistic inversion method to infer interface properties in the presence of heterogeneous sub-domains when little prior knowledge is available. Because of the limitations described above, we could not implement the structure-based prior approach and we resolve to constrain the global measure of model structure by penalizing model variability, with appropriate weights obtained via an empirical Bayes approach. Moreover, for computational reasons, we avoided the re-meshing of the unstructured grid used for the ERT model predictions, constraining the degrees of freedom for the inferred interface to the set of nodes of one subsurface discretization choice. We applied the proposed inversion algorithm on synthetic and field surface-based ERT datasets, with the aim of inferring regolith-bedrock interfaces. For both synthetic and field-based cases we were not able to sufficiently explore the posterior pdf for all the model parameters. Also, for the synthetic test case with heterogeneous layers we do not correctly infer the target interface at all locations. In both cases, the limitations related to the high complexity of the problem considered. Inversion for sharp boundaries between subsurface physical properties makes the likelihood surface highly irregular, especially for high-dimensional problems and large data sets with high signal-to-noise ratios, making it difficult for the MCMC chains to explore the posterior space exhaustively. Moreover, the local interface misdetection for the synthetic case with heterogeneous sub-domains is a consequence of the inherent equivalence problem of DC resistivity data and of the type of structureconstrained inversion routine applied, which tends to favor model realizations with too little complexity, especially where the data are weakly sensitive. Implementing alternative posterior sampling strategies (as for example parallel tempering or transdimensional inversion) could increase the MCMC exploration of the posterior model space, while to improve the interface detection a possible solution is to jointly invert additional observations in order to better constrain the model parameters.

In Chapter 4 we extended the empirical-Bayes-within-Gibbs framework to jointly invert ERT and seismic refraction datasets in presence of a common interface with unknown geometry. The joint inversion of the two sets of observations is proved to be more accurate in inferring the physical property fields, the location and geometry of the interface. Nevertheless, in case of the synthetic test case, we needed to inflate the modeled noise on the observations to account for modeling errors arising from the change in the mesh used to generate the observations. As pointed out in Chapter 3, the interfaces we are able to sample depend on the mesh, as they are bound to follow pre-defined node connections to avoid re-meshing at each proposed updated of the interface. For the synthetic test case presented, the target interface defined within the mesh used for the forward simulation is more irregular than the modeled one. This affects especially the algorithm used for computing the seismic refraction first arrival traveltimes, which is restricted to ray-paths that follow the elements connecting the nodes and therefore introduces errors in ray geometry that overestimate the traveltimes. We expect that improved results could be obtained by implementing a more accurate forward solver to simulate refraction seismic observations (e.g., fast marching approach) or a more detailed quantification of the model errors. Moreover, to increase the possible interface geometries the algorithm can sample, further work could be addressed to allow re-meshing at each interface proposal or to define the interface externally (e.g., polynomial formulation) and then interpolate the sub-domains it defines into the mesh used for model response computation.

Ultimately, prior assumptions have a strong influence on inversion model results. These assumptions can be included in the formulation of the inversion problem either implicitly or explicitly. When the inversion problem is reformulated in deterministic terms, the solution sought is the one able to explained the observations with less complexity. Through regularization in fact, it is implicitly assumed the model result to be smooth or similar to a chosen reference model. On the other hand, within the probabilistic formulation, prior assumptions are explicitly formulated through the shape and formulation of the prior pdf (e.g., in terms of model parameters or variability within a model). Also, including the geometry and location of subsurface interfaces as inversion parameter is an explicit prior assumption. Within this thesis the comparison of the results from various inversion routines show the strong influence that prior assumptions have on the resulting model(s). Therefore, once the data are collected, an investigation strategy that can be useful is to compare inversion model results obtained from routine characterized by different implicit or explicit prior assumptions. This can help to identify which are the model features that are constrained by the data and help proposing rival models for decision support.

A.1 Introduction

Geophysical imaging is largely employed to determine parameters and monitor processes for hydrological studies. The popularity of geophysics measurements lies in the large investigation coverage and minimal invasiveness. Once the data are acquired, we need to invert the collected observations to resolve an image of the subsurface, which is then used to estimate quantitative representations of the distribution of Earth properties. The attributes of these models depend on the inversion method used, the quality and amount of acquired data, the prior knowledge, and the degree of understanding of the investigated system.

Hydrogeophysical inversion is usually ill-posed. That is, the information is somehow "imperfect", so it is subject to multiple plausible interpretations. As a result, our interpreted Earth property distributions are neither unique nor stable [START_REF] Kabanikhin | Definitions and examples of inverse and ill-posed problems[END_REF].

To solve the inverse problem, we then need to include all of the available information about the subsurface.

The workflow of most hydrologic investigations begins with data collection. But, there is growing recognition that the nature of the problem, and the questions that are to be answered, should influence data collection and model development efforts [START_REF] Ferre | Revising the Relationship Between Data, Models and Decision-Making[END_REF][START_REF] Kikuchi | Towards Increased Used of Data Worth Analyses in Groundwater Studies[END_REF][START_REF] White | Forecast First: An Argument for Groundwater Modeling in Reverse[END_REF]. Similarly, the assumptions made in defining the inverse problem, and therefore the form taken by the solution, also impact the value of a prior knowledge and the way that it should be considered. The simplest inverse approach seeks a deterministic solution. These approaches assume negligible uncertainties in the observations and in the theoretically calculated model predictions, so that it is possible to find one optimized solution (the "best fit"to the data). Generally, this is achieved discarding all the other models through a regularization term, which penalizes the complexity of the solution and represents the a priori information. One common regularization approach assumes that the Earth's subsurface structure changes smoothly in space, so the regularization term penalizes models with large spatial variability [START_REF] Constable | Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data[END_REF].

With increased interest in quantifying prediction uncertainties, there has been growing use of probabilistic inverse methods (Woodbury & Ulrych, 2000). These approaches aim to describe the uncertainties on the solution together with the solution itself. For example, in Bayesian inference [START_REF] Stone | Bayes' Rule: A tutorial Introduction to Bayes Analysis[END_REF] uncertainty is express in terms of probabilities. Via Bayes' theorem, the probability of a hypothesis (e.g., a model of the Earth's subsurface property distribution) is updated based on the ability of the model to simulate observations as more information becomes available. The inversion problem is in this case postulated in terms of combining "state of information"on the model parameters, which are treated as random variables so that their "state of information"can be described through a probability density distribution (Tarantola & Valette, 1982).

All inverse problems begin with the choice of the complexity of the system [START_REF] Guthe | Defensible Model Complexity: A Call for Data-Based and Goal-Oriented Model Choice[END_REF]. The complexity includes the choice of the physical laws used to represent the system (forward model) and the values of the parameters needed to characterize the system with the chosen model (system parameterization). In Bayesian inversion, once this framework is defined, it is possible to formulate the a priori information on the unknown model parameters by a density distribution (prior distribution). This a priori information includes both existing measurements and our beliefs (conceptualization, structural assumptions) about the subsurface. Bayes' theorem is then used to combine the prior distribution with a likelihood function that describes the probability of observing additional collected data given the proposed model parameter values (our hypothesis). The resulting posterior distribution on the model parameters summarizes all the information available before the next phase of data collection [START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF]. This posterior density is only defined analytically for few specific cases. Generally, it is necessary to rely on sample-based methods such as Markov chain Monte Carlo (MCMC) [START_REF] Chib | Understanding the metropolis-hastings algorithm[END_REF][START_REF] Mosegaard | Monte Carlo sampling of solutions to inverse problems[END_REF], which result in an ensemble of models approximating the posterior distribution of the system properties.

In both the deterministic and the probabilistic approaches, the absence of reliable a priori information tends to lead to a bias in the solution. For example, if the subsurface is assumed to be smooth and this is imposed through the regularization constraint, deterministic inversion can result in an over-smooth image of the subsurface. On the other hand, in the probabilistic framework it is often seen as logical, or at least as more objective, to assume "uninformative"prior distributions. This is achieved by assuming that parameters are uncorrelated with uniform probability over the range of reasonable values. This choice of "uninformative"prior on the model parameters can also introduce a bias, maximizing the level of disorganization of the parameters and resulting in posterior model realizations that are too spatially variable to be geologically realistic. This may be acceptable if the models are only interpreted in terms of their average behavior; but it is misleading to pose the models as an ensemble of rival models [START_REF] Ferre | Revising the Relationship Between Data, Models and Decision-Making[END_REF] that can be used for decision support.

Selection of geologically unrealistic models occurs because the assumption of an uninformative prior prevents algorithms from sampling a large enough range of possible model structures. Therefore, it might be useful to constrain model proposals not only in terms of the parameter values, but also in terms of a global measure of model structure (e.g., roughness/smoothness or distance from a reference model) [START_REF] Pirot | Training Images as Advocacy Models to Describe Structural Variability of Spatio-Temporal Parameters[END_REF]. If available information is limited, it may be more effective to assume a uniform prior on the model structure rather than a uniform prior on the model parameters. This choice implies that all possible prior model realizations are not equally likely, but that is equally likely to sample a model realization with little spatial parameters variability (i.e., low structure) as one with a higher one. de Pasquale & Linde (2017), demonstrate that changing the focus of the model prior description from the parameters to the structure results in a posterior ensemble of models that better resemble realistic subsurface property distributions. Importantly, these more realistic models can be regarded as reasonable competing images of the domain (i.e., rival models).

A.2 A Synthetic Example

The synthetic test case presented in de Pasquale & [START_REF] De Pasquale | On structure-based priors in Bayesian geophysical inversion[END_REF], is useful to illustrate the effects that this change in the description of the prior model ensemble has on the inversion results. This example considers a crosshole ground penetrating radar (GPR) experiment, where the information about radar wave speeds is obtained by recording the first-arrival times of GPR waves transmitted in one borehole and received in another.

It is assumed that the information on the subsurface is unreliable or insufficient, so the prior is chosen to be uniform and therefore uninformative with regard to the parameter values or the global model structure. The subsurface is parametrized with a regular 10×10 Cartesian grid, where each cell (e.g., model parameter) is assigned a value of radar wave propagation speed.

The focus of standard Bayesian inversion is on the model parameters, which in this case are described with an uncorrelated log uniform distribution [START_REF] Jeffreys | Theory of probability[END_REF] within the range of 50 to 100 m/µs. The alternative to this classic Bayesian formulation is to sample a uniform distribution from the global model structure. In this example, the structure is defined by a Euclidean measure of global model roughness (zero in the case of uniform models). But, this could be substituted with any measure of parameter spatial variability within the model. In a broader sense, the posterior ensemble of models is in this case characterized by plausible structure behaviors. This property is crucial in many application, as for risk assessment, where many predictions of interest are in fact controlled by subsurface structure (e.g., faults or preferential flow paths).

In addition to improved ability to identify structure, shifting the focus of the prior description also reduces the computational effort required to sample the posterior distribution. This results because the ensemble of models that need to be sampled to represent the posterior ensemble in a statistically reliable manner is smaller when the models are described in terms of their structure than in the case of the standard parameter-focus Bayesian formulation. In other words, many random parameter distributions would have to be sampled to 'happen' to recover the structured subsurface if the prior assumption on the model structure is a random spatial distribution. and I thank them for that. I thank Sarah for all the dancing, Dafne for the warm welcomes, Samuel and Erika for their open-mindedness, Anna for her calm and deep listening, Flo, Tendai and all the people that gave me time to play music and/or climb together. I thank Gloria for her open-arms every time our path crosses, Laura and Brig for all the long skype talks here and there and Olivier, Jerome, Elliot, Diana and Justine for being an extension of home in the far Bretagne.

A crucial point of my personal growth in this Lausanne experience is Collectif R with its "little sister" UNIL Sans Frontiers. Graziella, Aline, Valentina, Bernard, Annabel, Stephan, Mathiew, Nicole, Judith, Pauline, Anne, Emilie, Rahel, Gabriel, Hamish, Axelle, Yann, Lucile, Gäel, Maila, Giulia... they are just some of the amazing people that have made and make these two realities exist. Together, besides being my personal "French course", they showed me that it is possible to fight for what I believe it is right and proved me that "we" is more than the sum of "us" and that solidarity is not just a word but is a way to live and to built a better society.

Family has always been the biggest support in my life and I'll never be thankful enough for all it gave me. I thank my father Nicola for supporting me in my education and his wife Vlada for supporting both of us in moments of need. I thank my aunt Chiara and uncle Antonio for always had been there and help filling my time with laugh and love. I thank Alexis, my brother in arms, that with a random conversation on his balcony brought me here in the first place and for being an example in many aspect of life I aim to grow into. Finally, I am deeply grateful to Dondi, Gió and Romana, three women that have been my life pillars and my first society. In the role of mother, sister and grandmother, they have guided me from the beginning till they could, giving me the solid roots from which all of this come from and I owe them all.

This project has been possible thank to the Swiss National Foundation, which financed four years of scientific research.

Lausanne, February 2019 G

Mots clefs :

 clefs Théorie de l'inversion, Distribution de probabilité, Tomographie, Interface.

Figure 1 . 1 :

 11 Figure 1.1: Schematic of Bayesian inference.

Figure 1 :

 1 Figure 1: Training image used to generate the reference model. Black channel structures have a velocity of 0.09 m/ns. The background velocity (white) has a velocity of 0.13 m/ns.

Figure 1 . 2 :

 12 Figure1is a channel-based training image from which we generate an unconditional realization, Figure2, using the single normal equation simulation algorithm, SNESIM[START_REF] Strebelle | Conditional simulation on complex geological structures using multiple point statistics[END_REF]). This will be our reference velocity model for a synthetic cross borehole inversion problem. Travel times are computed, traveling from the 20 sources located to the left in Figure2, to 40 receivers located to the right. 3% Gaussian noise is added to the synthetic travel-time delay data, Figure3, and used as observed data. We now consider solving this inverse problem, using the generalized Metropolis algorithm to generate samples of the a posteriori probability distribution, given the observed data, the assumed noise model, and an assumed prior model.We consider a number of prior models based on both 2-point and multiple-point based random models. All prior models are assumed to have the correct mean and variance, as obtained from the training image in Figure1. The f rst 6 prior models are based on 2-point random models. The pure nugget model assumes no spatial correlation, and thus all model parameters are a-priori considered uncorrelated. The 'Gau(1)', 'Gau(3)' and 'Gau(8)' prior models are based on a Gaussian covariance model with an isotropic range of 1m, 3m and 8m respectively. The SGSIM prior is 6

Figure 1 . 4 :

 14 Figure 1.4: Medium-pixel geometry and ray-path from transmitter (T x ) to receiver (R x ). For the ray i the distance traveled in a pixel j is denoted by a i j . Figure from: Hanafy & al Hagrey (2001)

Figure 1 . 5 :

 15 Figure 1.5: Example of plotting system for ERT pseudo section. V represents the high-impedance voltmeter measuring the voltage generate by injecting the current I into the subsurface.

Figure 2 . 1 :

 21 Figure 2.1: Examples of samples and estimates of S R2 model structure proposal pdfs. (a), (c) and (e) Forward step, q S(m cur r ) -→ S(m pr op ) . (b), (d) and (f) Backward step, q S(m pr op ) -→ S(m cur r ) . The models considered are 2D models discretized with an equally spaced 10 × 10 Cartesian grid. The histograms are based on P=1000 proposed model structure evaluations. The pdf estimators are the non-parametric kernel density estimation (KDE) and the parametric method of moments (Moments). The current model structure value (purple cross), the proposed model structure value (black cross) and the skewness are indicated.

  2.2). At x = 0 m, we have 19 sources equally spaced in depth, while 19 receivers (also equally spaced in depth) are disposed on x = 5 m. The model has a uniform background characterized by a radar wave speed of v = 75 m µs and a three-layer structure in the middle with a central area of higher wave speed (v = 90 m µs ) and slower upper and lower layers (v = 60 m µs ).

Figure 2 . 2 :

 22 Figure 2.2: A simple synthetic test model. Red crosses represents the 19 equally spaced sources while the blue crosses represents the 19 equally spaced receivers for the crosshole GPR experiment. Colors represent the radar wave speed (v).

Figure 2 . 3 :

 23 Figure 2.3: Prior density of model structure in case of UP (uncorrelated prior), EB (empirical Bayes) and SBP (structure-based prior) sampling. Structure-based prior results are obtained by approximating the proposal pdf of model structure with (a)-(d) parametric, and (e)-(h) nonparametric approaches. Yellow histograms represents the structure-based prior results after applying the exponent correction on the structural proposal ratio (ν in eq. 2.24) and the red histograms are the resulting prior pdf sampled when no correction is implemented (e.g., ν = 1). The orange area indicates the overlap of these two sampled distributions. (a) and (e) Results in case of roughness measure computed with l 1 -norm. (b) and (f ) Results in case of roughness measure computed with l 2 -norm. (c) and (g) Results in case of damping measure computed around m r e f = 75 m µs with l 1 -norm. (d) and (h) Results in case of damping measure computed around m r e f = 75 m µs with l 2 -norm.

  exp(l (m pr op |d)l (m cur r |d)) ; (2.26) and with the empirical Bayes structural-constrains by Rosas-Carbajal et al. (2014): α E B = mi n 1, exp l (m pr op |d)l (m cur r |d) + l og (c m pr op )l og (c m cur r ) , (2.27) where c m indicates either c m,2 (eq. 2.22) or c m,1 (eq. 2.23).

Figure 2 . 4 :

 24 Figure 2.4: Posterior density of model structure for structure-based prior acceptance ratios with parametric approximation of structure proposal pdf and S R2 measure of model structure with and without exponent correction (e.g., ν = 1.04 and ν = 1 in eq. 2.25, respectively). The actual model structure value for the known underlying model (Fig. 2.2) is indicated by a black cross.

Figure 2 . 5 :

 25 Figure 2.5: Posterior density of model structure for the different acceptance ratios after MCMC inversion of the synthetic dataset. (a) and (b) Results for roughness measures respectively, with l 1 -and l 2 -norms. (c) and (d) Results for damping measures computed around m r e f = 75 m µs with l 1 -and l 2 -norms respectively. The actual model structure value for the known underlying model (Fig. 2.2) is indicated by a black cross.

Figure 2 . 6 :

 26 Figure 2.6: Random prior model realizations based on the different sampling methods for the synthetic data example. (a) Uncorrelated prior random realizations. (b) Empirical Bayes random realizations with S R2 constraints. (c) and (d) Structure-based prior random realizations for an S R2 measure of model structure when the ratio in eq. (2.24) is approximated with a (c) parametric or (d) non-parametric approaches.

Figure 2 . 7 :

 27 Figure 2.7: Random posterior model realizations for the different acceptance ratios after MCMC inversion of the synthetic dataset. (a) Uncorrelated prior random realizations. (b) Empirical Bayes random realizations with S R2 constraints. (c) and (d) Structure-based prior random realizations for an S R2 measure of model structure when the ratio in eq. (2.25) is approximated with a (c) parametric or (d) non-parametric approaches.

Figure 2 . 8 :

 28 Figure 2.8: Means and standard deviations of the model realizations of MCMC inversion of the synthetic dataset. (a) and (b) Uncorrelated prior results. (c) and (d) Empirical Bayes results with S R2 constraints. (e)-(h) Structure-based prior results when the model structure proposal pdfs are estimated with (e) and (f) parametric, (g) and (h) non-parametric approach for an S R2 measure of model structure.

Figure 2 . 9 :

 29 Figure 2.9: Sliced plane of 3D GPR deterministic inversion of the dataset acquired close to the Thur River in northern Switzerland (modified after Doetsch et al. 2010 ).

  full MCMC inversion of the field data, we consider a domain of 6 m× 5.25 m discretized in 24 × 21 squared model cells with a 0.25 m edge length. The forward code utilized a finer mesh with 0.05 m edge lengths. We use 5 MCMC chains running in parallel and for each model proposal step, we randomly update 25% of the model parameters by adding perturbations from the current model with independent Gaussian realizations having standard deviation of 10 0.0035 m µs . To avoid model proposals outside of their prior range, we use reflection boundary handling. The marginal prior distributions are Jeffreys priors in a range from 60 m µs to 90 m µs and the model structure consideredfor the structure-based prior and empirical Bayes approaches is the roughness measure computed with l 2 -norm (S R2 ). In case of empirical Bayes constraints, the choice of the λ prior pdf is a Jeffreys prior with a range between 10 -3 and 10 -1 .The sampled prior model structure values for the different methods are depicted in Figs2.10(a) and (b) for the parametric and non-parametric empirical structure proposal pdf estimation, respectively. The prior sampling with our structure-based prior method with ν = 1 shows a very significant widening of the sampled model structures compared to the uncorrelated prior and empirical Bayes methods. Nevertheless, implementing an empirically-estimated exponential correction (ν = 1.1) leads to an even more uniformlooking distribution. Fig.2.10(c) depicts the posterior pdf of model structure for the case of parametric structural proposal approximation with and without exponent correction.

Figure 2 . 10 :

 210 Figure 2.10: Prior and posterior density of model structure for the different acceptance ratios after MCMC inversion of the real data. (a) and (b) Prior distributions of S R2 , in case of UP (uncorrelated prior), EB (empirical Bayes) and SBP (structure-based prior) sampling when the ratio in eq. (2.24) is approximated with a (a) parametric or (b) non-parametric approaches . (c) Posterior model structure pdf in case of parametric approximation of proposal ratio with and without exponent correction. (d) Posterior model structure pdfs for the different acceptance ratios.

Figure 2 . 11 :

 211 Figure 2.11: Random prior model realizations for the different sampling methods for the field data example. (a) Uncorrelated prior random realizations. (b) Empirical Bayes random realizations with S R2 constraints. (c) and (d) Structure-based prior random realizations for an S R2 measure of model structure when the proposal ratio in eq. (2.24) is approximated with a (c) parametric or (d) non-parametric approaches.

Figure 2 . 12 :

 212 Figure 2.12: Random posterior model realizations for the different sampling methods for MCMC inversion of the field data. (a) Uncorrelated prior random realizations. (b) Empirical Bayes random realizations with S R2 constraints. (c) and (d) Structure-based prior random realizations for an S R2 measure of model structure when the proposal ratio in eq. (2.25) is approximated with a (c) parametric or (d) non-parametric approaches.

Figure 2 . 13 :

 213 Figure 2.13: Means and standard deviations of the model realizations of MCMC inversion of the field data. (a) and (b) Uncorrelated prior results. (c) and (d) Empirical Bayes results with S R2 constraints. (e)-(h) Structure-based prior results when the model structure proposal pdfs are estimated with (e) and (f) parametric, (g) and (h) non-parametric approach for an S R2 measure of model structure.

Figure 2 . 14 :

 214 Figure 2.14: Sampled prior density of S R2 model structure when assuming an uncorrelated prior. (a) Linear y-axis, (b) logarithmic y-axis. Results are for 10 × 10 discretized model grid with parameters randomly chosen from uncorrelated Jeffreys priors with range within 50 and 100 m µs .The black cross corresponds to the probability density to sample a model structure which is half of the sampled mean values.

  2.24) helps to better sample the target prior structure distribution (Figs 2.3, 2.10a and b). Luckily, the posterior realizations are only weakly sensitive to this ad hoc correction (Figs 2.4 and 2.10c).

Figure 3 . 1 :

 31 Figure 3.1: Proposed model setting for inversion of geophysical datasets in the presence of an unknown interface, which divides the investigated subsurface in two sub-domains. (a) Physical properties, representative of the CZ (m C Z ) and bedrock (m b ) in the present example, defined on regular grids over the whole domain of investigation. (b) Unstructured mesh used for forward computation and (c) the corresponding interpolation of the regular grids (a) onto this mesh. (d) The interface I, defined by connected nodes of the forward computation mesh, is used to assign (e) the electrical properties used in the forward calculations for a given model (m = G(m C Z , m b , I)).

  quantifies the interface roughness and the likelihood ratio is computed using the model obtained after updating the interface (m cur r C Z , m cur r b , I pr op ) and the model accepted in the previous step of the chain (m cur r C Z , m cur r b , I cur r ).

,

  structure (S(m)) is computed for both domains (i.e., not only for the sub-domains found below or above the interface). We use different constraint functions for CZ and the bedrock, and sample two different hyper-parameters (λ C Z and λ b ). The likelihood ratio is computed between the model obtained by updating the two resistivity fields I cur r ) and the current model of the chain (m cur r C Z , m cur r b , I cur r ).

Figure 3 . 2 :

 32 Figure 3.2:Flow chart of the proposed empirical-Bayes-within-Gibbs algorithm. First, we initialize the model (as described in Fig.3.1) and evaluate its likelihood, physical properties structure and interface roughness. We then sequentially evaluate proposed interface and physical property updates in an iterative scheme until the chosen length of the chain is reached. Here, u ∼ U (0, 1) represents a randomly chosen number with uniform probability between 0 and 1, while the acceptance ratios α I and α P are described in eqs.3.7 and 3.8, respectively.

Figure 3 . 3 :

 33 Figure 3.3: (a) Resistivity model and electrode positions (black dots) used to generate the synthetic dataset with homogeneous sub-domains and (b) interpolation of (a) on the unstructured mesh used for ERT forward modeling. (c) Mean of posterior realizations from the empirical-Bayes-within-Gibbs chain and (d) model realization corresponding to the maximum likelihood. (e) Deterministic inversion result. In (c) and (e), the black line indicates the interface we aim to infer, while the blue line in (e) is the interface obtained using the maximum vertical gradient in the deterministic inversion result.

  3.4(b): m C Z ∈ [100, 3000] Ωm and m b ∈ [5000, 30000] Ωm. As described in Section 3.3, the resistivity fields are updated by randomly choosing 10 % of the corresponding model parameters (i.e., varying 10 % of the resistivity values in the regular grids in Fig. 3.1a). At each m C Z and m b model proposal step, the model parameters are updated using a Gaussian distribution centered on the current model with a standard deviation between 0.002 and 0.02: σ ∼ U (0.002, 0.02).

Figure 3 . 4 :

 34 Figure 3.4: (a) Resistivity model and electrode positions (black dots) used to generate the synthetic dataset with heterogeneous sub-domains and (b) interpolation of (a) on the unstructured mesh used for ERT forward modeling. (c) Noise-contaminated synthetic dataset obtained from (a) (645 apparent resistivity data points). (d) Distributed weights applied to the model parameter update.

  3.1b). For the prior definition of model structure or interface roughness weights used to constrain the model complexity (i.e., the hyperparameters λ C Z , λ b and λ I introduced in Sections 3.3 and 3.3), we follow[START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic data: methodology, model constraints and joint inversion with electrical resistivity data[END_REF] and use log-uniform priors with ranges between one quarter and six times the values found when eq. 3.4 is maximized for the true model. In the case of interface roughness, the maximization is done considering the value obtained by the interface of the model depicted in Fig.3.4(b). For CZ and bedrock, we use the same value that was obtained by considering the structure of the full model (i.e., computed from the model represented in Fig.3.4c). For the hyper-parameters, we consider Gaussian proposal distributions with a constant standard deviation of σ λ = 0.1.To sample the posterior space, we run three independent chains for 10 6 iterations with each iteration taking ≈ 0.6 seconds. The starting model of the three chains are presented in Figs.3.5(a), (d) and (g), where the logarithmic resistivity fields are sampled with uniform uncorrelated prior probability. The corresponding model predictions have weighted root mean squared errors (WRMSE) of 28, 34 and 35, which is much higher than the target value of 1 on which the posterior distribution is centered. That is, the misfits of the initial models are some 30 times larger than the posterior realizations. The last posterior model realizations(Figs. 3.5b, e, h), as well as the maximum likelihood models (Figs.3.5c, f, i) present resistivity fields and bedrock topography that are largely consistent with each other and with the underlying true model (Fig.3.4b).

Figure 3 . 5 :

 35 Figure 3.5: (a), (d) and (g) The three initial models used for the synthetic test case with heterogeneous sub-domains. (b), (e) and (h) The last corresponding realizations of each MCMC chain. (c), (f) and (i) Corresponding maximum likelihood realizations.

Figure 3 . 6 :

 36 Figure 3.6: Behavior of the three independent chains (each colors refers to a specific chain) for the synthetic example with heterogeneous sub-domains. Sampled model structure for (a) CZ, (b) bedrock and (c) regolith-bedrock interface, with the black cross referring to the values for the actual interface. Sampled hyper parameters penalizing (d) CZ, (e) bedrock and (f) interface roughness, with black crosses referring to the optimal hyper parameter values for the true model. Acceptance rate of (g) the interface and (h) the physical property update. (a)-(f ) share the same y-label, so does (g) and (h). (i) Log likelihood evolution during the course of the inversion with target log-likelihood (dashed black line)

Figure 3 . 7 :

 37 Figure 3.7: Posterior model results for empirical-Bayes-within-Gibbs inversion of the synthetic dataset with heterogeneous sub-domains. (a), (d) and (g) represent the mean of the posterior model realizations for each chain with the actual interface superimposed (black line). (b), (e) and (h) show the standard deviation of the electrical resistivity divided by the corresponding mean values. (c), (f) and (i) show the probability of each cell of the mesh being part of the bedrock. (l), (m), (n) and (o) are vertical resistivity profiles at 50, 150, 250 and 350 meters along the profile; we show in blue, red and green the mean model results (solid lines) of each chain with the inferred posterior ranges (shadow areas), while the black dashed lines represent the synthetic model.

Figure 3 . 8 :

 38 Figure 3.8: (a) Mean of the posterior model realizations for the three empirical-Bayes-within-Gibbs chains and (b) its model discrepancy compared to the true model. (c) Deterministic inversion result and its (d) model discrepancy. The black line in (a) and (c) it indicates the interface we aim to infer (i.e., Fig. 3.4b), the gray line in (a) is the mean of the inferred posterior interfaces, while in (c) represents the interface obtained with the maximum vertical gradient method based on the deterministic inversion result. (e), (f ), (g) and (h) show vertical resistivity profiles respectively at 50, 150, 250 and 350 meters along the profile. The red solid lines represent the mean model of empirical-Bayes-within-Gibbs inversion and the red shadow areas indicate the inferred posterior ranges, the blue dotted lines are obtained from the deterministic inversion, while the black dashed lines represent the underlying synthetic model.

  The final ranges of m C Z ∈ [80, 10000] Ωm and m b ∈

Figure 3 . 9 :

 39 Figure 3.9: (a), (d) and (g) The three initial models used for the field-data inversion using empirical-Bayes-within-Gibbs sampling. (b), (e) and (h) The last corresponding realizations of each MCMC chain. (c), (f) and (i) Corresponding maximum likelihood realizations.

Figure 3 .

 3 Figure 3.10: Characteristics of the three independent chains (each colors refers to a specific chain) for the field example. Sampled model structure for (a) CZ, (b) bedrock and (c) regolithbedrock interface. Sampled hyper parameters penalizing (d) CZ, (e) bedrock and (f) interface roughness, with black crosses referring to the optimal hyper parameter values for the true model. Acceptance rate of (g) the interface and (h) the physical property update. (a)-(f ) share the same y-label, so does (g) and (h). (i) Log likelihood evolution during the course of the inversion with target log-likelihood (dashed black line)

Figure 3 . 11 :

 311 Figure 3.11: Posterior model results for empirical-Bayes-within-Gibbs inversion of the field dataset. (a), (d) and (g) represent the mean of the posterior model realizations for each chain. (b), (e) and (h) show the standard deviation of the electrical resistivities divided by the corresponding mean values. (c), (f ) and (i) show the probability for each cell of the forward mesh to be part of the bedrock. (l), (m), (n) and (o) are vertical resistivity profiles at 50, 150, 250 and 350 meters along the profile. Blue, red and green indicate the mean model of each chain (solid lines) with the inferred posterior ranges (shadow areas).

Figure 3 . 12 :

 312 Figure 3.12: Posterior model results for one-domain empirical Bayes inversion of the field dataset. (a), (c) and (e) represent the mean of the posterior model realizations for each chain, while (b), (d) and (f) show the standard deviation of the electrical resistivities divided by the corresponding mean values.(g), (h), (i) and (l), vertical resistivity profiles respectively at 50, 150, 250 and 350 meters along the profile. Blue, red and green indicate the mean model of each chain (solid lines) with the inferred posterior ranges (shadow areas).

Figure 3 . 13 :

 313 Figure 3.13: Comparison of field data inversion results. (a) Mean of the posterior model realizations for the three empirical-Bayes-within-Gibbs chains, (b) mean of posterior model realizations for the three one-domain empirical Bayes chains and (c) deterministic inversion result. In (a) the black line is the mean of the inferred posterior interfaces, while in (b) and (c) it represents the interfaces obtained by the maximum vertical gradient method. (d), (e), (f) and (g) show vertical resistivity profiles at 50, 150, 250 and 350 meters along the profile. Red solid lines are used to represent the mean of the empirical-Bayes-within-Gibbs inversion result and the red shadow areas the inferred posterior ranges, the black dashed lines with the gray shadow areas represent the one-domain empirical Bayes inversion results and the blue dotted lines are obtained from the deterministic inversion.

  3.13c) and a smaller range of resistivity values: m ∈ [150, 7000] Ω m compared with the probabilistic inversion results (Figs. 3.13a and b) with m ∈ [100, 25000

Figure 4 . 1 :

 41 Figure 4.1: Proposed model setting for joint inversion of DC resistivity and refraction seismic datasets in the presence of an unknown interface (I), which divides the investigated subsurface into two sub-domains. (a) Electrical resistivity fields of the CZ and bedrock. (b) Seismic P-wave velocity fields of the CZ and bedrock. (c) The interface, defined by connected nodes in the forward computation mesh, is used to assign (d) the electrical resistivity (m r ) and (e) P-wave velocity (m v ) used for forward calculations.

  Fig. 4.3(a) represents the electrode positions and the resistivity model used to generate the synthetic ERT datasets, while Fig.4.3(c) shows the interpolation (through the nearest neighbor method) of the input resistivity model on the mesh used for the forward computations. Fig.4.3(b) displays the geophone positions and the P-wave velocity model, while Fig.4.3(d) represents its interpolation on the mesh used for the forward computation. Fig.4.3(e) shows the ERT inversion model result obtained by a l 1 -norm mimicking smoothness-constrained deterministic inversion, while Fig.4.3(f ) shows the corresponding deterministic inversion model result for the refraction seismic observations. In Figs.4.3(e) and (f) we indicate the interface obtained by using the maximum vertical gradient method on the deterministic inversion results (following the procedure described by[START_REF] Chambers | Bedrock detection beneath river terrace deposits using three-dimensional electrical resistivity tomography[END_REF]. Moreover, the deterministic inversion model results are depicted with a transparency that is inversely proportional to the data coverage. This highlights that the DC resistivity data have virtually no coverage after 400 m along the profile, simply because no electrodes were located in this area. Also, the coverage of the seismic refraction measurements deteriorates significantly after 400 m. Therefore, for the analysis and interpretation of the inversion results, we focus only on the results obtained in the first 400 m. Finally, in Figs.4.3(g-h) we show a measure of discrepancy between the target model and the inversion results:

Figure 4 . 3 :

 43 Figure 4.3: (a) Resistivity model and electrode positions (black dots) used to generate the synthetic DC resistivity dataset and (c) interpolation of (a) on the unstructured mesh used for ERT forward modeling. (b) Velocity model and geophone positions (small black dots) used to generate the synthetic refraction seismic dataset and (d) interpolation of (b) on the unstructured mesh used for the refraction seismic forward modeling. Deterministic inversion results based on the (e) DC resistivity and (f) seismic refraction observations. In (e) and (f) the transparency is inversely proportional to the data coverage, the black line represents the interface obtained by the maximum vertical gradient method and the red line is the interface we aim to infer. Model discrepancies for deterministic inversion of the (g) DC resistivitiy and (h) seismic refraction observations.

Fig. 4 .

 4 Fig.4.7 shows vertical profiles of the velocity and resistivity fields obtained for the different probabilistic and deterministic inversion approaches plotted together with the true model parameter values. Both the individual and joint probabilistic inversions present jumps in the physical property values at locations that are mostly in agreement with the target interface location. In case of deterministic inversion results, the physical property values smoothly increase with depth. This results in significant underestimation of the bedrock properties.To quantify the ability of the different inversion approaches to infer the target interface (I t ar g et ), we compute the average distance between I t ar g et and the sampled (in case of probabilistic approaches) or inferred one (in case of deterministic inversion results):

Figure 4 . 4 :

 44 Figure 4.4: Last posterior model realizations from probabilistic inversion of synthetic datasets. (a), (c) and (e) represent resistivity model realizations from each MCMC chain for individual ERT inversion, while (b), (d) and (f) show such realizations in case of joint inversion. (g), (i) and (m) are P-velocity model realizations from each MCMC chain for individual refraction seismic inversion, while (h), (l) and (n) show such realizations in case of joint inversion. The inferred interface of each realization is represented with a black line.

Figure 4 . 5 :

 45 Figure 4.5: Vertical resistivity and P-wave velocity profiles obtained by probabilistic inversion of the synthetic datasets. Blue, red and green indicate the mean model of each chain (solid lines) with the inferred posterior ranges (shadow areas) and the black dashed lines represent the target model. (a)-(d) Inferred resistivity values from probabilistic inversion of DC resistivity data. (e)-(h) Inferred resistivity and (i)-(n) velocity values from probabilistic joint inversion. (o)-(r) Inferred velocity values from probabilistic inversion of individual seismic refraction data.

Figure 4 . 6 :

 46 Figure 4.6: Summary of posterior realizations for the three MCMC chains used for each synthetic inversion case. Mean of the posterior model realizations for (a) individual DC resistivity and (i) seismic refraction inversions; (d) and (g) for joint inversion. The mean inferred interfaces are represented with a black line, while the interface we aim to infer is shown in blue. Standard deviation of the resistivities, divided by the corresponding mean values of (b) individual DC resistivity inversion and (e) joint inversion. Standard deviation of the velocities, divided by the corresponding mean values for probabilistic inversion of (h) individual refraction seismic data and (l) joint inversion. Probability for each cell of the forward mesh to be part of the bedrock in case of probabilistic inversion of individual (c) DC resistivity (m) and refraction seismic data and in case of (f) joint inversion. Model discrepancies for the resistivity fields inferred from probabilistic inversion of (n) individual ERT data and (o) joint inversion. Relative model errors for the velocity fields inferred from probabilistic inversion of (p) individual seismic refraction data and (q) joint inversion.

Figure 4 . 7 :

 47 Figure 4.7: Vertical resistivity and P-wave velocity profiles from inversion of the synthetic datasets. Red solid lines represent the mean of the probabilistic joint inversion results and the red shadow areas the inferred posterior ranges, the green solid lines with the green shadow areas represent individual dataset probabilistic inversion results, the blue dotted lines are obtained from deterministic inversion and the black dashed lines are the target profiles. (a)-(d) Inferred resistivity values and (e)-(h) inferred velocity values.

Figure 4 . 8 :

 48 Figure 4.8: (a) ERT observations and (b) refraction seismic first-arrival travel times from the Calhoun CZO field site. Deterministic inversion results based on (c) DC resistivity and (d) seismic refraction observations. In (c) and (d) the models transparency is inversely proportional to the data coverage and the black line represents the interface obtained by the maximum vertical gradient method.

Fig. 4 .

 4 Fig.4.9 presents the last posterior realization obtained for each chain and inversion case, while in Fig.4.10 we show the corresponding vertical resistivity and velocity profiles at 50 m, 150 m, 250 m and 350 m along the profile. The posterior model realizations obtained by inverting the DC resistivity dataset alone were already presented in de[START_REF] De Pasquale | Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data[END_REF], who stressed the similarity between the model results in terms of both posterior realizations and vertical profiles. This is similarly seen in the posterior velocity realizations obtained by inversion of the refraction seismic data alone(Figs. 4.9g,i,m) and in the vertical profiles, where the inferred interface location is manifested by velocity jumps that appear at very similar positions. Particularly, all the three chains suggest an interface that mirrors the topography, with bedrock rising almost to the surface at locations where the surface topography creates valleys at around 80 m and 350 m. The posterior models obtained by joint inversion are represented in terms of both resistivity (Figs.4.9b,d,f) and velocity (Figs. 4.9h,i,n) fields. As for the single-dataset inversion results, the upper part of the domain is well defined for both property fields (i.e., upper 40 m). Moreover, comparing the joint inversion results with the individual inversion results(Figs. 4.9b,d,f with Figs. 4.9a,c,e and Figs. 4.9h,l,n with Figs. 4.9g,i,m), we see that the bedrock topography is mainly constrained by the seismic refraction observations. The DC resistivity data has the strongest sensitivity to heterogeneity within the bedrock layer. This is seen in the well-defined low resistivity zone of the bedrock around 300 m while the velocity fields of the bedrock are more homogeneous. Finally, around 80 m along the profile, all the posterior velocity models obtained by joint inversion show a slightly deeper bedrock interface than the ones inferred from the inversion of the seismic observations alone.

Figure 4 . 9 :

 49 Figure 4.9: Last posterior model realizations from probabilistic inversion of the Calhoun datasets. (a), (c) and (e) represent resistivity model realizations from each of the MCMC chains for individual DC resistivity inversion, while (b), (d) and (f) show such realizations in case of joint inversion. (g), (i) and (m) are P-velocity model realizations from each of the MCMC chains for individual refraction seismic inversion, while (h), (l) and (n) show such realizations in case of joint inversion. The inferred interface of each realization is represented with a black line.

Figure 4 . 10 :

 410 Figure 4.10: Vertical resistivity and P-wave velocity profiles from probabilistic inversion of the Calhoun datasets. Blue, red and green indicate the mean model of each chain (solid lines) with the inferred posterior ranges (shadow areas). (a)-(d) Inferred resistivity values from probabilistic inversion of DC resistivity data alone. (e)-(h) Inferred resistivity and (i)-(n) velocity values from probabilistic joint inversion. (o)-(r) Inferred velocity values from probabilistic inversion of the refraction seismic dataset alone.

Figure 4 . 11 :

 411 Figure 4.11: Summary of posterior realizations for the three MCMC chains used for the inversion of the Calhoun datasets. Mean of the posterior model realizations when inverting for individual (a) DC resistivity and (i) refraction seismic dataset; (d) and (g) when jointly inverting the two datasets. The mean inferred interfaces are represented with a black line. Standard deviation of the resistivities, divided by the corresponding mean values for probabilistic inversion of (b) individual DC resistivity data and (e) for joint inversion. Standard deviation of the velocities, divided by the corresponding mean values for probabilistic inversion of the (h) individual refraction seismic dataset and (l) for joint inversion of DC resistivity data and refraction seismic dataset. Probability for each cell of the forward mesh to be part of the bedrock in case of probabilistic inversion of individual (c) DC resistivity, (m) refraction seismic data (m) and (f) joint inversion.

Figure 4 . 12 :

 412 Figure 4.12: Vertical resistivity and P-wave velocity profiles obtained by inversion of the Calhoun datasets. Red solid lines are used to represent the mean of probabilistic joint inversion results and the red shadow areas the inferred posterior ranges, the green solid lines with the green shadow areas represent individual dataset probabilistic inversion results and the blue dotted lines are obtained from deterministic inversion. (a)-(d) Inferred resistivity values and (e)-(h) Inferred velocity values.

Figure 4 . 13 :

 413 Figure 4.13: Posterior model realizations from probabilistic inversion of synthetic datasets when the model errors are not taken into account. (a) Resistivity model realization inferred from inversion of individual DC resistivity data and (b) velocity model realizations inferred from inversion of individual refraction seismic data. (c) Resistivity and (d) velocity model realizations inferred from joint inversion. The inferred interface of each realization is represented with a black line and the target interface with a blue line.

Figure 4 . 14 :

 414 Figure 4.14: Inferred mean values and ranges of (a) resistivity and (b) velocity at the interface location as obtained from probabilistic joint inversion of the Calhoun datasets. In both plots, the black line represents surface topography.

  Fig.A.1 shows the results of these two alternative descriptions of the prior. In the case of parameter focused prior description (marginal uniform uncorrelated priors), the sampled model structure distribution is very narrow and unable to represent low or high model structure values when considering a finite number of samples. This behavior is then reflected in the prior model realizations (Fig.A.1b), which are similar in terms of spatial parameter variability within the model. When the focus of the prior description is instead shifted on the global measure of structure of the models, it is possible to sample a much larger range of structure values, resulting in more diverse spatial parameter variability within the prior ensemble of models (Fig.A.1c).

Figure A. 1 :

 1 Figure A.1: (a) Sampled prior distribution on S R2 (model structure: Euclidean measure of global model roughness) when considering marginal uncorrelated priors on the model parameters and on the model structure. Random prior model realizations in case of marginal uncorrelated uniform priors on the model parameters (b) and on the model structure (c).

Fig.A. 2 (

 2 Fig.A.2(a) shows the toy-model assumed to represent a structured subsurface. This model is used to simulate the times for the radar waves to travel from transmitter to receiver antennas situated at the opposite edges of the domain. The information carried by these simulated measurements are combined with the prior distribution via Bayes' theorem, to obtain a sample of the posterior distribution. This resulting posterior ensemble of models is strongly influenced by the choice of the priori model description. In Fig.A.2(b) the posterior results are represented in terms of sampled model structure, showing that the parameter-focused prior description severely overestimates the actual

Figure A. 2 :

 2 Figure A.2: (a)Synthetic toy model on 2D square domain with lateral dimension of 5m. At x=0m there are 19 GPR transmitter antennas equally spaced in depth, while at x=5m are disposed 19 receiver antennas, also equally spaced in depth. (b)Sampled posterior distribution on S R2 (model structure: Euclidean measure of global model roughness) when considering marginal uncorrelated priors on the model parameters and on the model structure. Random posterior model realizations in case of marginal uncorrelated uniform priors on the model parameters(c) and on the model structure(d).

  Nevertheless, both Rosas-Carbajal et al. (2014) and de Pasquale & Linde (2017) highlighted that this type of structure-constrained inversion routine tends to favor model realizations with too little complexity, especially where the observations are weakly sensitive.

  Flow chart for Markov chain Monte Carlo joint inversion, following the principle of the Metropolis within Gibbs presented by[START_REF] Iglesias | Well-posed Bayesian geometric inverse problems arising in subsurface flow[END_REF]. d A and d B are the observations, respectively sensitive to the subsurface physical properties parametrized with m A and m B , and to the subsurface geometry parametrized with Geom.

				Initialization
				7 3 4			7 3 ?
						9:;7 3
				2 33 4 = 2(7 3 4 , 9:;7 3 |= 4 )	2 33 ? = 2(7 3 ? , 9:;7 3 |= ? )
				Geometry update
				9:;7 3		9:;7 @
			2 3@ 4 = 2(7 3 4 , 9:;7 @ |= 4 )	2 3@ ? = 2(7 3 ? , 9:;7 @ |= ? )
				Accept the update? 2 33 4 2 33 ? D ∝ 2 3@ 4 2 3@ ?	YES	Update geometry: 9:;7 3 = 9:;7 @
						NO
				Keep 9:;7 3
				Physical properties updates
			7 3 4	7 @ 4			7 3 ?	7 @ ?
			2 @3 4 = 2(7 @ 4 , 9:;7 3 |= 4 )			2 @3 ? = 2(7 @ ? , 9:;7 3 |= ? )
	Update 7 4 : 7 3 F = 7 @ 4	YES	Accept the update? D ∝ 2 @3 4 2 33 4			Accept the update? D ∝ 2 @3 ? 2 33 ?	YES	Update 7 ? : 7 3 ? = 7 @ ?
				NO			NO
				Keep 7 3 4			Keep 7 3 G
			END	YES	Desired length of the chain reached?	NO
	Figure 1.3:

Table 2 . 1 :

 21 Convergence diagnosis for the synthetic data inversion. UP stands for uncorrelated prior implementation, EB refers to the empirical Bayes method and SBP to our structure-based prior sampling (where par. stands for parametric approximation of model structure proposal pdf and non-par. stands for the non-parametric approximation). Structure indicates model structure metric (Section 2.3). The R-convergence shows at which iteration the potential scale reduction factor gets smaller than 1.2. AR refers to the average acceptance rate of the chains after the burn-in, r mse D is the data misfit root mean squared error and r mse M is the root mean squared model error defined in eq. (2.28).

	Method	Structure R-convergence AR [%] mean r mse D mean r mse M Std r mse M
	UP		5.8 × 10 5	39	1.03	3.2
	EB	S R1	2.6 × 10 5	26	1.08	2.3
		S R2	1.7 × 10 5	21	1.09	1.3
		S D1	2.1 × 10 5	32	1.07	2.6
		S D2	1.7 × 10 5	24	1.08	1.7
	SBP: par.	S R1	3.1 × 10 5	31	1.04	2.9
		S R2	2.9 × 10 5	36	1.04	1.9
		S D1	4.6 × 10 5	24	1.05	2.1
		S D2	4.4 × 10 5	23	1.05	1.9
	SBP: non-par.	S R1	3.8 × 10 5	30	1.04	1.6
		S R2	3.4 × 10 5	35	1.04	1.5
		S D1	4.2 × 10 5	26	1.05	2.1
		S D2

Table 2 . 2 :

 22 Convergence diagnosis for the field data inversion. UP stands for uncorrelated prior implementation, EB refers to the empirical Bayes method and SBP to our structure-based prior sampling (where par. stands for parametric approximation of model structure proposal pdf and non-par. stands for the non-parametric approximation). Structure indicates the model structure metric (Section 2.3). The R-convergence shows at which iteration the potential scale reduction factor gets smaller than 1.2. AR refers to the averaged acceptance rate of the chains after the burn-in and r mse D is the data misfit root mean squared error.

	Method	Structure R-convergence AR [%] mean r mse D
	UP		9.6 × 10 5	42	0.70
	EB	S R2	6.4 × 10 5	27	0.92
	SBP: par.	S R2	7.1 × 10 5	38	0.79
	SBP: non-par.	S R2	7.7 × 10 5	36	0.76

de Pasquale, G., N. Linde, J.Doetsch and W.S: Holbrook (2019). Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data. Geophysical Journal International, doi: 10.1093/gji/ggz055

de Pasquale, G., N. Linde and A. Greenwood. Joint probabilistic inversion of ERT and seismic refraction data: a case-study for bedrock/regolith interface delineation.

de Pasquale, G. (2017). Changing the Prior Model Description in Bayesian Inversion of Hydrogeophysics Dataset. Groundwater, 55(5) 651-655, https://doi.org/10.1111/gwat.12547

Acknowledgements

There is some kind of chaotic flow in what brought me here and now: a moving frame with a lot of people animating it. I owe them all, and even though there is only my name at the beginning of this manuscript, lots of this work is mine as much as theirs, therefore, it feels right to add their names at the end of it.

Abstract

Geophysical methods provide remotely sensed data that are sensitive to subsurface properties and interfaces. Knowledge about discontinuities is important throughout the Earth sciences: for example, the saltwater/freshwater interface in coastal areas drive mixing processes; the temporal development of the discontinuity between frozen and unfrozen ground is indicative of permafrost development; and the regolith-bedrock interface often plays a predominant role in both landslide and critical-zone investigations. Accurate detection of subsurface boundaries and their geometry is challenging when using common inversion routines that rely on smoothness constraints that smear out any naturally-occurring interfaces. Moreover, uncertainty quantification of interface geometry based on such inversions is very difficult. In this paper, we present a probabilistic formulation and solution to the geophysical inverse problem of inferring interfaces in the presence of significant subsurface heterogeneity. We implement an empirical-Bayes-within-Gibbs formulation that separates the interface and physical property updates within a Markov chain Monte Carlo scheme. Both the interface and the physical properties of the two sub-domains are constrained to favor smooth spatial transitions and pre-defined property bounds. Our methodology is demonstrated on synthetic and actual surface-based electrical resistivity tomography datasets, with the aim of inferring regolith-bedrock interfaces. Even if we are unable to achieve formal convergence of the Markov chains for all model parameters, we demonstrate that the proposed algorithm offers distinct advantages compared to manual or algorithm-based interface detection using deterministic geophysical tomograms. Moreover, we obtain more reliable estimates of bedrock resistivity and its spatial variations.

Introduction

Near-surface geophysical datasets are inherently sensitive to physical properties of the subsurface and their spatial discontinuities. The underlying interests and reasons for mapping subsurface boundaries in the Earth sciences is discipline-dependent. For instance, knowledge about saltwater/freshwater interfaces and their evolution in coastal areas is crucial for informed groundwater management [START_REF] Bear | Seawater Intrusion in Coastal Aquifers-Concepts, Methods and Practices[END_REF][START_REF] Kim | Tidal effects on variations of fresh-saltwater interface and groundwater flow in a multilayered coastal aquifer on volcanic island (Jeju Usland, Korea)[END_REF][START_REF] Werner | Seawater intrusion processes, investigation and management: Recent advances and future challenges[END_REF]. Both aquifer overexploitation and catastrophic events (e.g., hurricanes, tsunamis) cause seawater intrusion that leads to contamination of freshwater reserves. The mixing zone between lighter freshwater and denser seawater is very difficult to study using traditional borehole monitoring data [START_REF] Carrera | Computational and conceptual issues in the calibration of seawater intrusion models[END_REF]. Geophysical methods that are sensitive to electrical conductivity are suitable to detect and study this interface and its dynamics [START_REF] Nguyen | Characterization of seawater intrusion using 2D electrical imaging[END_REF][START_REF] Falgas | Monitoring freshwater-seawater interface dynamics with audiomagnetotelluric data[END_REF]). Another interface of importance is the one between frozen and unfrozen ground [START_REF] Kneisel | Advances in Geophysical Methods for Permafrost Investigations[END_REF]. The evolution of permafrost impacts climate both directly by influenc-

Chapter 4

Joint probabilistic inversion of ERT and seismic refraction data: a case-study for bedrock/regolith interface delineation Giulia de Pasquale, Niklas Linde and Andrew Greenwood.

Submitted 1 in Journal of Applied Geophysics

Appendix A

Changing the Prior Model Description in Bayesian Inversion of Hydrogeological Dataset

Giulia de Pasquale Published 1 in Groundwaterl and herein slightly adapted for this thesis.

A.3 Conclusions

Modelers seek tools that allow them to avoid subjectivity in model construction. But, as with the influence of bias on decision making [START_REF] Ferre | Revising the Relationship Between Data, Models and Decision-Making[END_REF], there are often hidden biases that can even affect these objective tools. One example is the impact of the assumption of the nature of the prior information in Bayesian inversion. While it may seem most unbiased to assume an "uninformative" prior distribution on the model parameters, this actually impacts the distribution of model structures that are included in the posterior model ensemble. Changing the focus of the a priori model description from the parameter values (traditional approach) to a measure of global model structure alleviates this problem. By making this relatively simple change, more diverse model structures are selected, which may greatly improve our ability to identify potentially important system behaviors with the model ensemble. As an added benefit, this change also reduces the computational effort needed to generate the model ensemble.

Niklas Linde has been a great mentor, and I'm not just talking about his incredible ability of exploring and imagining alternative solutions to recurrent problems, or the patient in correcting my spelling mistakes and dealing with all the "Giulia's factor" I carefully hide within the codes. He also gave me the chance to explore and learn and question: basically to grow both as a researcher and as a person. I thank him for that. There is a certain number of researchers that collaborated with me and I will never thank them enough, because finding an answer together is always more efficient than by oneself. I thank David Ginsbourger for the mathematical inputs and discussions, Joseph Doetsch and Florian Wagner for giving me the hints I needed to get started with pybert and pygimli, and especially to keep answering the most trivial questions all the way through. I thank Daniel Fernàndez-Garcia for sharing codes and literature on the kernel density estimation, W. Steven Holbrook for sharing field-data and insights on them and Andrew Greenwood for helping with their processing. I also thank the committee members Thomas Mejer Hansen and Grégoire Mariéthoz for the interesting discussion during my defense and the president of the jury, Christian Kull for his availability.

Géopolis is a concentration of inspiring human beings and I had the chance to spend four years wandering around it, collecting a great deal of conversations and stories. I thank Jürg and Alejandro for the time shared discussing small and big issues in the office and Carlotta for all the rainy adventures during field courses. I thank all my research group for the scientific feedbacks and the coffee breaks and all the "donkeys" for filling the time in Géopolis and elsewhere with laugh and friendship. I thank Sté for have been my personal translator to proper French, Gogol for all the concerts and festivals he brought or follow me into, Nath for being a source of inspiration, Nico for helping me to find my motivation at some point, Benjy for the running-hugs and Mike for the bear-like ones. I thank Vjeran for the great high fives, Christophe for introducing me to the art of parking security, Femi for showing me that you can always use a bit more of your energy and Seb for confirming it. I thank Ignes for the good mood, Leanne for her crazy ideas, Mor1 for the good vibes towards everything and everyone. I thank Martino for the bouldering sessions, Josh for his contagious enthusiasm, Zoneibe for the smiles, Maga and Nadja for listening to everything I needed to say, even and especially when I didn't want to say it.

Outside the scientific encounters, often developed in friendship, I had the chance of crossing paths with Michael, Aki, Bernard and Basile, who have been home for me in avenue de la Harpe • de Pasquale, G., N. Linde, J. Doetsch and S. Holbrook (2018). Probabilistic inversion of geophysical datasets for subsurface interface detection. European Geosciences Union General Assembly 2018, Oral presentation. • de Pasquale, G., N. Linde, J. Doetsch and S. Holbrook (2018). Probabilistic joint inversion of ERT and seismic refraction data to detect subsurface interfaces. PhD conference "Water-Earth Systems", Zurich.