# 1D geological imaging of the subsurface from geophysical data with Bayesian Evidential Learning

Hadrien Michel, Frédéric Nguyen, Thomas Kremer, Ann Elen, Thomas Hermans

► To cite this version:

Hadrien Michel, Frédéric Nguyen, Thomas Kremer, Ann Elen, Thomas Hermans. 1D geological imaging of the subsurface from geophysical data with Bayesian Evidential Learning. Computers & Geosciences, 2020, 138, pp.104456. 10.1016/j.cageo.2020.104456. hal-04342450

HAL Id: hal-04342450

https://hal.science/hal-04342450v1

Submitted on 5 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 1D geological imaging of the subsurface
from geophysical data with Bayesian
Evidential Learning.

**Authors**<fnref n="1" />:

Hadrien MICHEL<fnref n="1" />,²,³, Frédéric NGUYEN<fnref n="1" />, Thomas KREMER<fnref n="1" />, Ann ELEN⁴, Thomas HERMANS²

(1) University of Liège, Urban and Environmental Engineering Department, Faculty of Applied Sciences, Liège, Belgium, (2) Ghent University, Department of Geology, Ghent, Belgium, (3) F.R.S.-FNRS (Fonds de la Recherche Scientifique), Brussels, Belgium, (4) KU Leuven, Department of Earth and Environmental Sciences, Leuven, Belgium

**Corresponding author:**

Hadrien MICHEL
University of Liège
Allée de la découverte, 9
B-4000, Liège (Belgium)
E-mail: hadrien.michel@uliege.be

**Declaration of interest:** none

**Link to the code:** https://github.com/hadrienmichel/BEL1D

**Highlights:**

- A new approach for 1D imaging from geophysical data.

- A set of MATLAB graphical user interfaces is provided.

- Estimation of the uncertainty on the model parameters.

- Consistent results for SNMR.

- BEL1D can be generalized to any type of 1D geophysical data.

<footnote n="1">¹ Authorship statement: Hadrien Michel developed the BEL1D codes and produced the results presented in this paper. Thomas Kremer has significantly contributed and improved all the SNMR related elements developed in this paper. The Mont Rigi SNMR dataset was acquired and pre-processed by Ann Elen for her Master Thesis. Frédéric Nguyen and Thomas Hermans both co-supervised the whole work and are at the basis of methodological development.</footnote>
**Abstract:**

Imaging the subsurface of the Earth is of prime concern in geosciences. In this scope, geophysics offers a wide range of methods that are able to produce models of the subsurface, classically through inversion processes. Deterministic inversions lack the ability to produce satisfactory quantifications of uncertainty, whereas stochastic inversions are often computationally demanding. In this paper, a new method to interpret geophysical data is proposed in order to produce 1D imaging of the subsurface along with the uncertainty on the associated parameters. This new approach called Bayesian Evidential Learning 1D imaging (BEL1D) relies on the constitution of statistics-based relationships between simulated data and associated model parameters. The method is applied to surface nuclear magnetic resonance for both a numerical example and field data. The obtained results are compared to the solutions provided by other approaches for the interpretation of these datasets, to demonstrate the robustness of BEL1D. Although this contribution demonstrates the framework for surface nuclear magnetic resonance geophysical data, it is not restricted to this type of data but can be applied to any 1D inverse problem.

**Keywords:** 1D imaging, uncertainty, Bayesian Evidential Learning, machine learning, SNMR, MATLAB.
1. Introduction

Imaging the subsurface of the Earth through geophysical methods allows reducing costly invasive investigations or is sometimes the only way to investigate the Earth. Behind the geophysical images are models conceptualized in terms of physical properties (density, seismic velocities or resistivity) which may be converted to properties of interests (hydraulic conductivity, porosity or temperature) using petrophysical relationships. Decision-making, whether planning additional investigations for a scientific project, exploiting a geo-resource, or constructing a building, is often a driver for geophysical methods (Lawyer et al., 2001). Uncertainty assessment on those geophysical models is therefore a key concern to geophysicists, as it provides a much more consistent information in decision-making processes (Scheidt et al., 2018).

Computation methods to process geophysical data into models (other than simple mapping of the data) are either deterministic, often using the regularization approach, or stochastic, often using the Bayesian framework. Deterministic methods are limited in terms of uncertainty quantification and often rely on linear data error propagation based on a linearized approximation if needed (Aster et al., 2013; Kemna et al., 2007). On the other hand, stochastic methods may compute sets of models based on numerous inversions of perturbed datasets, or on the exploration of a prior model space using Markov chain Monte Carlo (McMC) methods (Aster et al., 2013; Sambridge, 2002).

Recently, due to the rise of computer performances, stochastic methods have been applied to geophysics and are able to quantify more properly uncertainty. Most algorithms rely on Markov chains Monte Carlo (McMC) methods in order to explore the prior model space (Sambridge, 2002). Those methods offer the advantage to provide a much more complete appraisal of the ensemble of possible solutions. Nevertheless, stochastic methods are often
linked to high computation costs, since they require a large number of forward model runs to sample the posterior distribution in the prior model space, explaining why they are still often limited to scientific or numerical studies (e.g., Irving and Singha, 2010; Trainor-Guitton and Hoversten, 2011). Approximation of the posterior can be obtained (under some assumptions) through ensemble-based inversion (e.g., Bobe et al., 2019). Even more recently, machine learning also emerged as a potential contender to replace the inversion process, typically yielding to deterministic-like results due to the nature of neural networks and such (e.g., Laloy et al., 2019, 2017; Yang and Ma, 2019).

In this paper, we present Bayesian Evidential Learning (BEL) (Scheidt et al., 2018) as an innovative method for the unidimensional geophysical imaging of the subsurface. Bayesian Evidential Learning (BEL) is a general framework to handle data in geosciences, from the optimization of the data acquisition to the data interpretation, prior falsification, global sensitivity analysis, and the prediction of relevant model responses (Scheidt et al., 2018).

We will refer to the approach presented here as Bayesian Evidential Learning 1D imaging (BEL1D). It uses statistics-based relationships between forecast variables and data learned from realizations sampled from a prior distribution. As opposed to inversion methods, this approach does not require the time-consuming inversion of the dataset, but rather numerous (and parallelized) runs of the faster forward model. So far, BEL has been used to predict the dynamical response of subsurface models from well monitoring data (Satija and Caers, 2015), push-pull tests (Hermans et al., 2019) or time-lapse geophysical experiments (Hermans et al., 2018), i.e. model predictions mostly characterized by smooth variations in time and/or in space. BEL was also applied to show the non-uniqueness of the gravity inverse problem (Phelps et al., 2018).

In this contribution, we extend BEL to 1D geophysical imaging. Such 1D imaging arises in many geophysical methods where forward models assume a succession of homogeneous
horizontal layers such as surface nuclear magnetic resonance (SNMR e.g., Behroozmand et al., 2015), surface seismic waves (e.g., Socco et al., 2010), electromagnetic surveys (e.g., Hanssens et al., 2019; Li et al., 2018), or vertical electrical sounding (e.g., Jha et al., 2008).

The main contributions of this paper compared to previous works are:

1. The application of BEL1D for static geophysical imaging. Here, we predict directly the posterior distribution of spatially distributed geophysical model parameters.

2. The estimation of the data-prediction relationship using kernel density estimation instead of Gaussian regression, because the former can be applied when the relationship deviates from linearity and Gaussian assumptions.

3. The validation of the approach by comparing it with state-of-the-art McMC

4. The link made between distance-based global sensitivity analysis and the canonical correlation analysis, yielding similar results.

5. The analysis of the number of samples in the prior.

This paper first introduces our implementation of BEL: Bayesian Evidential Learning 1D imaging (BEL1D). The method is then validated within a numerical example. To demonstrate the broad applicability of our approach, we also provide a validation within a non-geophysical related inverse problem (supplementary material: "Testing BEL1D on a synthetic mind experiment: oscillations of a pendulum"). Then, we present an example of BEL1D applied to field data, providing evidence that the method is already mature. Finally, we discuss on elements governing the accuracy of BEL1D, such as the number of models sampled in the prior and the choice of the bandwidth for kernel density estimation.
71

## 2. BEL for 1D imaging: methodology

72

73

74

75

76

77

78

79

80

In this work, the focus is on computing geophysical models and assessing the uncertainty of the inferred models (Phelps et al., 2018). This interpretation part of BEL is also known as the prediction-focused approach (PFA), introduced by Scheidt et al. (2015). Contrary to deterministic approaches, this method does not rely on the stabilization of the ill-posed inverse problem, such as regularization for example, that imposes non-realistic constraint to the solution. It rather relies on the constitution of statistical relationships between predictors (the set of parameters of interest to the end-user) and data. These relationships are originating from models and (numerical) simulations of the geophysical problem that reflect the available prior knowledge.

81

82

83

84

BEL has already been applied previously to geophysical data – time-lapse electrical resistivity – by Hermans et al. (2016). Hermans et al. (2016) demonstrated the applicability of BEL to predict variations of subsurface physical properties with time-lapse electrical resistivity data. More broadly, they also demonstrated that BEL was a possible tool for geophysicists.

85

86

87

88

89

90

91

92

93

94

However, Hermans et al. (2016) presented a scheme that required numerical simulations of groundwater flow and transport and petrophysical relationships specific to hydrogeophysical monitoring (i.e. time-lapse data), limiting the extend of the geophysical prior model. Solving the time-lapse problem in geophysics is easier because we can invert for changes in the model (e.g., Kemna et al., 2002; Nguyen et al., 2016). In Hermans et al., (2016), a petrophysical relationship and the background distribution of resistivity were taken as known, which simplified the inference of the posterior distribution. Static data are also characterized by a higher noise level and thus more challenging to process (LaBrecque et al., 1996; Lesparre et al., 2017). This is also the first application where BEL is used to directly estimate spatially distributed model parameters. In our case, the forecast variables are static geophysical models,
directly linked to the geophysical experiment, hence no petrophysical relationship is required.

Such situations generally bear more prior uncertainty and wider variability. Moreover, our aim

is to propose a fast generation of subsurface models along with uncertainties. The adapted

BEL1D method consists of six steps (Figure 1):

1. Definition of the prior uncertainty based on prior field knowledge, generation of

prior models (yellow box in Figure 1) and associated data (blue box in Figure 1)

by forward modeling

2. Reduction of the dimensionality of the data using principal component analysis

(PCA)

3. Constitution of statistical relationships between the model parameters and the

reduced data (green cloud in Figure 1), using canonical correlation analysis

(CCA)

4. Generation of posterior distributions for the model parameters in reduced space

by constraining the bivariate distributions to field data (red box in Figure 1)

using kernel density estimators (contrary to previously used Gaussian process

regression (Hermans et al., 2016))

5. Sampling of the constituted (non-Gaussian) distributions

6. Back-transformation of the samples into the original space, delivering a set of

1D models of the subsurface that are constrained to the knowledge of the

geophysical data (purple box in Figure 1).

Each step is explained in detail in the following sections.
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of the BEL1D process. See text for abbreviation.</figcaption></figure>

Although BEL1D can be applied to any 1D inverse problem, for the clarity of the paper, we briefly introduce SNMR, the specific method chosen to illustrate the BEL1D method. SNMR is a geophysical technique dedicated to the detection and characterization of aquifers, whose main advantage and particularity is that the measured signal (an exponential decay) is directly linked to the amount of water present in the subsurface. For more details about the nature and origin of the SNMR signal, we refer to Behroozmand et al. (2015). The main goal of SNMR experiments is to retrieve the distribution of the water content and the relaxation time with depth. The latter depends on the pore space geometry and can be related to hydrodynamic parameters such as the porosity or the permeability. The acquisition of SNMR data on the field consists in the injection of an alternative current in a transmission loop (typically circular or in an eight-shape) directly followed by the measurement of the induced current in the reception loop (often the same as the transmission loop). The injection is tuned at the resonance frequency
of water (called the Larmor frequency), perturbing the quantum state of protons contained in
water molecules. Once the injection is stopped, the protons are relaxed back to their original
state, inducing a current with a decreasing amplitude with time in the receiver loop. This
response is called the free-induction decay (FID) and constitutes the collected data. The
injection is repeated with different pulse moments, sounding different zones of the subsurface.

2.1. Step 1: prior realizations and associated data

In BEL1D, the models are described in a classic manner with a finite number of layers
$N_L$, the last one being a half-space. For each layer, the $N_p$ physical properties related to the
geophysical experiment and layer thickness are assigned. This means that a model can be
described with a low number of parameters: $(N_L \times (N_p + 1)) - 1 = q$. Defining the prior
model space $f_H(\mathbf{m})$, i.e. the space describing the prior knowledge of the subsurface physical
properties before data acquisition, consists in assigning a distribution to the values of each of
those $q$ parameters. The shapes of those distributions are unconstrained and must best represent
the prior knowledge of the survey site, either originating from previous experiments or from
general geological and geophysical considerations. Once the prior model space is defined, we
generate random models within the prior boundaries, hereafter called prior realization:
$$ \mathbf{m} (q \times 1) \in f_H(\mathbf{m}). $$

For example, the prior model space for SNMR experiments is defined using both the
water content (in %) and relaxation time (in ms) distributions along depth. This leads to a finite
number of layers, of unknown thickness, and of unknown water content and relaxation time,
thus to the description of a model with only five parameters for a 2-layer model.

For all prior realizations, we compute their associated response using the same
acquisition parameters as the experimental conditions. Doing so, we produce the synthetic data
vector: $d(l \times 1)$. To perform this operation, we use the (non-linear) forward operators ($K$) of
the geophysical problem of interest such that

$$d = K(m) \qquad (1)$$

For more details about the forward operator associated with the SNMR method, we refer
to Hertrich et al. (2007). Often, the dimensionality of the data ($l$) is relatively large. For SNMR,
the data space is composed of tens of time-decaying signals (one for each pulse moment) with
sampling rates on the order of 100 Hz measured for about 1 seconds, thus resulting in several
thousands of dimensions.

## 2.2. Step 2: reduction of the dimensionality of the data using principal component analysis

At this stage, we have generated the “model space” which contains the prior realizations,
and the “data space” which contains the synthetic geophysical data sets associated with each
model. To explore the statistical relationships that exist between these two spaces, one must
first ensure that their dimensionality is sufficiently low (typically about 10 dimensions per
space) to allow for reasonable computational costs. In the case of SNMR, as in most other cases
if the prior models space remains simple, the model space will generally have a reasonable
number of dimensions, thus no dimension reduction is required ($m^f = m(q \times 1)$). As stated
previously, the data space is generally large and cannot be handled properly under reasonable
computational costs. In this case, the Principle Component Analysis (PCA) technique can be
used to transform the data space and reduce its dimensionality while preserving its original
variability.

The PCA method seeks for linear combinations of variables that maximize the
variability in the first dimensions using eigenvalue decomposition (Krzanowski, 2000). This
operation results in PCA scores, the values of the transformed data. This way, the first
dimensions are the most informative about the dataset, and the others only represents a small
amount of the variability. Those later can be discarded, to gain memory and ease following computations. In the SNMR context, the same process is applied to the original data space and performs very well since the dimensions are often highly correlated. The obtained reduced data space has a dimensionality $k \ll l$ and is noted $d^f$ ($k \times 1$). Note that in other geophysical problems, this PCA step may not be necessary if the dimensions of both the model and data space are low or may also be necessary for both the model and data spaces if their dimensions are too high. Noise analysis using Monte Carlo simulations (Hermans et al., 2016) provides covariance between the different PCA components representing the data ($C^f$). This translates in uncertainties on the PCA scores of the test data. The detailed procedure for noise propagation is presented in subsection 2.7.

2.3. Step 3: constitution of statistical relationships between the models parameters and the reduced data, using canonical correlation analysis

BEL1D relies on the constitution of statistics-based relationships between earth models and simulated data. This step is the core of BEL1D. Several methods are suitable to derive such a relationship, but we choose canonical correlation analysis (CCA) for its simplicity. In essence, CCA transforms the (possibly PCA-reduced) model space ($m^f$) and the PCA data space ($d^f$), so that the canonical correlation (an approach formula of the correlation) between the resulting CCA spaces ($m^c$ – Equation 2 – and $d^c$ – Equation 3) is maximized (Krzanowski, 2000). Each of the dimensions of the CCA spaces are orthogonal, meaning that there is no redundancy between the dimensions. We obtain:

$$d^c = d^f A^T \qquad (2)$$

$$m^c = m^f B^T \qquad (3)$$

where $A$ is a matrix of dimensions $m \times k$ and $B$ is a matrix of dimensions $m \times q$ ($m$ is the minimum between $k$ and $q$). Graphically, we can observe the statistical relationship between $d^c$ and $m^c$ as both spaces have the same dimensions ($1 \times m$). As long as $k$ is larger than $q$, it
200 is possible to back-transform the models from the CCA model space to the original model space
201 using the inverse of **B** (Hermans et al., 2016).

202 2.4. Step 4: generation of posterior distributions for the model parameters in CCA
203 space using kernel density estimators

204 In canonically correlated space, the correlation between dimensions is maximized.
205 Hence, it is possible to produce a meaningful statistical description of the bi-variate distribution.

206 This enables to determine a posterior distribution in the CCA model space constrained by the
207 knowledge of field data. This posterior distribution represents the probability density function
208 for the model parameters in the CCA model space. To approximate the posterior distribution
209 corresponding to the observed data $f_H(\mathbf{m}^c|\mathbf{d}_{\text{obs}}^c)$, we use a kernel density estimator (KDE) with
210 a Gaussian kernel (Wand and Jones, 1993) similar to Hermans et al. (2019) where they faced
211 non-linearity issues when validating BEL for hydrogeology. Gaussian regression (Hermans et
212 al., 2016; Satija and Caers, 2015) was tested; but in many cases, CCA did not yield perfectly
213 linearly correlated relationships with Gaussian prior distributions, two necessary conditions to
214 apply Gaussian regression. We thus implement KDE for its robustness to estimate posterior
215 distribution in a large variety of situations. KDE computes, for a given location in space, the
216 sum of the contribution of each point from the CCA space. We use a multi-Gaussian kernel
217 centered on the points with bandwidths chosen accordingly to the point density (see subsection
218 5.2). The resulting distributions in the CCA model space are not constrained to any given
219 distribution of known shape. The process is illustrated in the center of Figure 1.

220 2.5. Step 5: sampling of the constituted distributions

221 The posterior distribution is not constrained to a known distribution. Therefore,
222 sampling is done through the inverse transform sampling method (Devroye, 1986). This
223 sampling procedure benefits from the properties of the cumulative distribution function that
links a uniformly distributed variable to any distribution. The sampler thus transforms a value
sampled in a uniform distribution to the corresponding value in the original distribution via the
inverse of the cumulative distribution function. We can now easily generate a set of samples
from the posterior distributions in the reduced space:

$$
m_{post}^{c} \in f_{H}(m^{c} | d_{obs}^{c}) \qquad (4)
$$

2.6. Step 6: back-transformation of the samples into original space

The set of samples in CCA model space can be back-transformed into the original model
space, using the inverse of the **B** matrix (q × q). We thus have the posterior in the original
physical space:

$$
m_{post} = m_{post}^{f} = B^{-1} m_{post}^{c} \quad (5)
$$

Here, in the SNMR context, since the model space was not reduced using principle
component analysis, the models in PCA space are the same as the ones in original space.
Because Gaussian regression or KDE with Gaussian kernel has no limit on the value of sampled
parameters, a few samples might occasionally be located outside the prior model space, in
particular when the model is at the extremity of the prior. Those models are removed from the
sampled space as they do not correspond to the definition of the problem.

238 2.7. Noise propagation in BEL1D

In order to investigate the effect of noise on BEL1D, one can analyze the variation of
the PCA scores to the estimated field noise level (if the latter can be estimated). If the scores of
the first *k* components are not varying in a significant way, then, one can assume that noise has
no impact on the imaging process. On the other hand, if the noise impacts the scores
significantly, BEL1D considers this impact by propagating uncertainty (represented by *C*<sup>*f*</sup> in
Equation 6) on the reduced observed data in the canonically correlated space, leading
potentially to more complex computation of the posterior distributions (Hermans et al., 2016).
This uncertainty is propagated in the CCA dimension according to linear error propagation from equation (2):

$$C^c = AC^f A^T \quad (6)$$

To estimate $C^f$, we explore the PCA scores and their modifications when noise is present in the dataset. To do so, we take the data associated to each of the prior realizations as a basis and compute the PCA scores when one model is perturbed with a noise level similar to the data. The operation is repeated for a significant sub-sample (in the latter example, we used 50 randomly selected models) of the prior samples allowing to derive the PCA score covariance matrix $C^f$.

Contrary to Gaussian regression where the addition of covariance to the observed reduced data only slightly increases the complexity of the formula for the sampling parameters (Hermans et al., 2016), the addition of uncertainty on the reduced observed data highly hinders performances of the Kernel density estimation. Assigning uncertainty to the observed data (X in Figure 2) requires the computation of the kernel density estimation of the reduced posterior distributions for multiple sets of randomly sampled $d_{obs}^c + \epsilon$, then processing to sampling and back transformation from those sets of distributions to constitute models of the posterior distributions. This can be avoided by enlarging the bandwidth of the kernel density estimator (Bowman and Azzalini, 1997) in the X direction (in the CCA data space) according to $C^c$. This accounts for the uncertainty on the exact position of the data in the reduced space and therefore accounts for noise, without adding complexity to the computations, as is shown in Figure 2.
<figure><img src="image_3.png" /><figcaption>Figure 2: Example of the effect of noise on the computation of the kernel density estimator for the probability density functions. On the left, a random dataset that is highly correlated and the associated true value of X. On the right, the estimated distributions for the values of Y constrained to the known value of X and its associated error of 1. The computation time t is given for the different cases.</figcaption></figure>

## 3. Numerical example

In this section we will present results obtained by applying BEL1D to numerical examples of SNMR. The model that we will use corresponds to a classical context with two layers: a first unsaturated layer and a saturated half-space (Table 1). A more complex case is presented in section 4, and a non-geophysical related inverse problem is presented in supplementary material (“Testing BEL1D on a synthetic mind experiment: oscillations of a pendulum”).

<table><caption>Table 1: Characteristics of the synthetic model</caption><thead><tr><th></th><th>Thickness (m)</th><th>Water content (%)</th><th>Relaxation time (ms)</th></tr></thead><tbody><tr><td>Layer 1</td><td>5</td><td>5</td><td>100</td></tr><tr><td>Half-space</td><td>/</td><td>25</td><td>200</td></tr></tbody></table>

### 3.1. Computation of the model space and data space

The data are simulated for SNMR using the MRSMatlab toolbox (Müller-Petke et al., 2016). For the simulated experiment, the circular transmitter is the same as the receiver (the
loops have a 50 m diameter with 2 turns). To simplify the computations, the earth is set as
resistive, neglecting the effect of electrical conductivity in signal attenuation and phase shift.
Therefore, the signal is purely real. However, the process presented in this paper is not limited
to real-valued data. SNMR is known to be impacted significantly by electromagnetic (EM)
noise (e.g., Behroozmand et al., 2016). Different types of EM noise can affect the SNMR data.
However, when most of these components have been filtered out, the SNMR signal will remain
contaminated by a Gaussian noise distribution, whose magnitude cannot be reduced further. To
represent this, we chose to add a Gaussian noise level to the dataset, corresponding to a mean
signal-to-noise ratio of 10 (35 nV), which is a relatively high value in usual SNMR surveys.

<figure><img src="image_4.png" /><figcaption>Figure 3: Simulated data for the SNMR experiment. The black curves represent the noisy dataset and the red curves the noise-free one.</figcaption></figure>

As can be observed in Figure 3, the initial amplitude of the signal is around 500 nV, and
signal-to-noise ratio of 10 produces quite noisy signals, but the decaying behavior remains
observable.

The prior model space reflects the prior knowledge of the study site. In our case, we will
assume that it is known that the subsurface is composed of one layer resting over an infinite
half-space. Then, the rest of the properties will be defined as uniformly distributed variables varying in the intervals described in Table 2. The prior intervals are not necessarily centered on the actual value of the parameter. This is to better simulate real situations, where the user will only have rough information about the subsurface constitution. 5000 prior realizations are generated using a Latin-hypercube sampler (McKay et al., 1979). The forward model that was used to simulate the test data is reused here to generate the synthetic data set.

<table><caption>Table 2: Prior model space for the SNMR experiment</caption><thead><tr><th rowspan="2"></th><th colspan="2">Thickness [m]</th><th colspan="2">Water content [%]</th><th colspan="2">Relaxation time [ms]</th></tr><tr><th>Minimum</th><th>Maximum</th><th>Minimum</th><th>Maximum</th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>Layer 1</td><td>2.5</td><td>7.5</td><td>3.5</td><td>10</td><td>5</td><td>350</td></tr><tr><td>Half-space</td><td>/</td><td>/</td><td>10</td><td>30</td><td>5</td><td>350</td></tr></tbody></table>

## 3.2. Global sensitivity analysis and dimension reduction

We performed a distance-based global sensitivity analysis (Park et al., 2016) whose results are presented in Figure 4 (Left). The parameters corresponding to the first layer ($W_1$
309

and $T_{2,1}^*$) are not sensitive to the SNMR signal in this specific experimental configuration and
should therefore lead to a poorly reduced uncertainty on those parameters.

310

<figure><img src="image_5.png" /></figure>

311

312

Figure 4: Left: multivariate sensitivity analysis performed on the simulated data (the red bars correspond to parameters that have an effect on the data and the blue ones to parameters that have a negligible effect). Right: combinations of parameters representing the canonically correlated model space dimensions (the values are scaled for visual purposes)

313

314

315

<figure><img src="image_6.png" /></figure>

316

317

Figure 5: SNMR - Data space versus model space in the canonically correlated space. The blue dots represent the prior models and data spaces and the orange line represents the exact value of the test data in the reduced space.

318

319
From the initial 11.904 dimensions in the dataspace, PCA enabled a reduction to only 5 dimensions, while managing to keep more than 90% of the variability. The canonical correlation analysis is shown in Figure 5. It is observed that (1) the correlation between the models and the data is larger for the first 3 CCA dimensions, (2) the correlation is not linear from the 2nd dimension and (3) from the 4th dimension, the points are scattered and not correlated. The fact that the last dimensions are scattered (note that each dimension is a linear combination of the model parameters) informs us on the inability of the data to efficiently predict all the parameters simultaneously. This is confirmed by the linear combinations that constitute the CCA model space dimensions (Figure 4, Right): we observe that the three last reduced dimensions are mainly linear combinations of the insensitive values (Figure 4, Left). Note that non-linear statistical techniques (e.g., Gorban et al., 2008; Lawrence, 2012) could potentially lead to better characterization of the data-model relationship.
<figure><img src="image_7.png" /><figcaption>Figure 6: Results of the BEL1D process on the synthetic noise-free data from the 2 layers model. The obtained distributions are the posterior distributions estimated from the noisy dataset in red. The sampled prior model space is represented by the blue distributions. The test values are represented with dashed red lines in the histogram plots. Results of the McMC algorithm DREAM are presented in yellow.</figcaption></figure>

333

334

335

336

337

338

339

340

341

342

BEL1D is applied to the simulated noisy dataset. Figure 6 shows the prior and posterior distributions of parameters obtained through BEL1D. We observe that the uncertainty on the half-space parameters is dramatically reduced through BEL1D. On the other hand, the insensitive parameters $T_{2,1}^*$ and $W_1$ are poorly reduced. This insensitivity to some of the parameters is also clearly identified in the pendulum example (Supplementary Material). As
can be observed, the test value of the parameter is always contained in the posterior distribution.

The multivariate analysis of the parameters shows that most of the parameters are independent.

However, some correlations are observed for parameters that are related to the same layer ($T_{2,2}^*$ and $W_2$ for example). Figure 7 (Left) presents another representation of the posterior results, where the data RMS of posterior samples are color-coded, and where each model is a set of linked parameters. The presented results show that a large part of the posterior models has a very low RMS error. The trend in the water content set through the prior model space is respected. In contrast, the trend for the relaxation time is not resolved by BEL1D. The analysis of the RMS and the comparison with DREAM (Figure 7, Right) shows that is related to the insensitivity of the data and not to the method itself.

<figure><img src="image_8.png" /><figcaption>Figure 7: Left: results of the BEL1D process on the synthetic noisy data from the 2 layers model. Each displayed model is presented with its associated RMS error value. In order to enhance visualisation, each colour is associated to the same number of models. The dashed white lines represent the test values and the dashed black lines represent the extent of the prior model space. The similar graph is also presented for DREAM results (Right).</figcaption></figure>

Finally, we applied the DREAM McMC algorithm (Vrugt, 2016) to the same dataset.
We applied a Gaussian likelihood with 15 generations and 10000 chains and a burn-in period of 5000 chains. However, to achieve convergence towards reasonable uncertainties, changes to the definition of the dataset are required. We had to use gate integration, similarly to what is
applied in MRSMatlab (Müller-Petke et al., 2016) in order to attribute higher loads to the first part of the decay curves, thus lowering the impact of noise on the estimation. For the presented results, we use 5 gates for each pulse moment. However, we must note that the results of DREAM are highly sensitive to the number of gates: with 1 or 2 gates, nearly no reduction of uncertainty is observed, whereas with a too large number of gates (to the limit with no gate integration) false posterior distributions are observed (the test values are outside of the predicted posterior). The results are presented alongside the BEL1D results in Figure 6. They show that BEL1D tends to have a slightly lower performance than DREAM. This is related to the approximation of the problem in a lower dimension space where the actual data misfit is never computed. Applying a cutoff on the RMS values observed in BEL1D would result in the same distribution (Figure 7).

<figure><img src="image_9.png" /><figcaption>Figure 8: Impact of noise on the PCA scores (the first 5, kept after dimensionality reduction), covariance matrix for the PCA scores impacted by noise (left) and the corresponding covariance in CCA space (right). This first matrix is estimated through the mean of the covariance matrices observed between the noise-free scores and the noisy ones.</figcaption></figure>

To estimate the posterior distributions, we needed to account for the impact of noise in the process. The matrix $C^f$ presented in Figure 8 shows very low values for the covariance. However, once transformed into CCA space ($C^c$, Equation 6), we observe that the effect of
382

noise, even if limited, is still observable. From a default bandwidth of 0.01 in the X direction for the kernel density estimation, we computed bandwidths ranging from 0.017 to 0.059.

383

384

Those increased bandwidths explain the presence of some outliers in the dataspace (Figure 7, Left**Érreur** ! **Source du renvoi introuvable**.). When observing the CCA space (Figure 5), it is noticeable that the second dimension (mainly representing the decay time of the second layer) shows models that are outside the main distribution but close to the position of the reduced test data. When using the bandwidth that takes into account the impact of noise, those models are englobed into the posterior. Nonetheless, those models correspond to large RMS error and are easily identified (Figure 7).

385

386

387

388

389

390

391

An example with more than two layers is presented in supplementary material ("Application of BEL1D to a four-layer model").

392

393

## 4. Case study: Mont Rigi (Belgium)

394

The Mont Rigi is located in the Belgian Fagnes region, in the Eastern part of the country.

395

This site presents an ideal case study for SNMR as it is remote and far from any electromagnetic noise sources, due to its natural reserve classification.

396

397

Geologically, the site is characterized by a metric peat layer on top of a Cambrian bedrock (La Venne formation) known as an aquiclude with a very low water content (Gilson et al., 2017). In contrast, peat is known to present very high water contents with observed total porosities around 90% (Wastiaux, 2008). According to previous GPR exploration (Wastiaux and Schumacker, 2003), the peat layer at the experimental site should have a thickness between 2.5 and 4.5 meters. The SNMR response should therefore allow to easily distinguish the two layers with a properly designed experiment. The designed on-site experiment consisted of a classical coincident loops transmitter/receiver couple with a diameter of 20 meters.

398

399

400

401

402

403

404
The raw data have been preprocessed to lower the impact of noise with despiking, harmonic modelling, reference noise cancellation and despiking again (Müller-Petke et al., 2016). The signals have been further truncated to 0.2 seconds to decrease the impact of noise, mostly present at large timeframes. The resulting noise amplitude is quite low for SNMR (18 nV). Then, the signals have been inverted using the QT Inversion (Müller-Petke and Yaramanci, 2010), to constitute a benchmark to compare to the results of BEL1D. We used a smooth-mono model description (smooth models with a single relaxation time for a given depth) and a regularization parameter of 6000 for both water content and relaxation time using the L-curve criteria (Aster et al., 2013). Then, BEL1D was applied with the uniformly distributed prior model space described in Table 3. We sampled 5000 prior realizations and produced 5000 models for the posterior.

<table><caption>Table 3: Prior model space for the Mont Rigi SNMR experiment</caption><thead><tr><th rowspan="2"></th><th colspan="2">Thickness [m]</th><th colspan="2">Water content [%]</th><th colspan="2">Relaxation time [ms]</th></tr><tr><th>Minimum</th><th>Maximum</th><th>Minimum</th><th>Maximum</th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>Layer 1</td><td>0</td><td>7.5</td><td>30</td><td>80</td><td>0</td><td>200</td></tr><tr><td>Layer 2</td><td>/</td><td>/</td><td>0</td><td>15</td><td>100</td><td>400</td></tr></tbody></table>

The CCA space analysis (Figure 9) shows that we could expect a significant reduction of uncertainty for some parameters, as the first three dimensions show high correlations. The first dimension, dominated by the water content of the half-space shows an especially narrow relationship between the reduced data and model space. On the other hand, the last dimension, mostly represented by the water content of the first layer, shows no specific correlation. We should therefore expect a much better reduction for the parameters of the half-space than for the first layer.
<figure><img src="image_10.png" /><figcaption>Figure 9: Mont Rigi – Data versus model space in the CCA space. The blue dots represent each prior realizations and the orange line the position of the field data.</figcaption></figure>

The noise level (after denoising) impact on BEL1D was assessed and resulted in significant changes in bandwidths for the kernel density estimation, between 0.04 and 0.12 instead of the default 0.01 value.
<figure><img src="image_11.png" /><figcaption>Figure 10: SNMR at Mont Rigi - Prior (blue) and posterior (orange) distributions of the parameters.</figcaption></figure>

The analysis of the posterior parameters model space (Figure 10) shows that the reduction of uncertainty from the prior model space to the posterior is significant and produces consistent results, when compared to the geological context and the QT inversion solution (dashed gray lines in Figure 11). The thickness of the peat layer tends to be smaller than the maximum estimate (4.5 m) from the GPR interpretation (Wastiaux and Schumacker, 2003). The use of a narrower prior taking into account the GPR information would prevent such a behavior. This is illustrated by the decrease of the maximum value of the thickness for the first
layer when the water content of the aforementioned layer increases. This behavior arises from
the non-uniqueness of the SNMR response stated above.

In the posterior models produced by BEL1D (Figure 11), a large majority shows a root
mean square error lower than the noise level in the dataset, implying that the set of probable
models provides an efficient estimation of the site parameters. The RMS errors of the posterior
samples show that they all explain the data to a similar level; the uncertainty is thus intrinsic to
the data set and the non-unicity of the solution. The non-uniqueness of the solution is easily
observable, with two distinct behaviors: one with a median water content and a thick first layer
and another with a very high water content but a thinner first layer. Those two cases show
similar RMS errors, hence are explaining the data at the same level. The probability of each
posterior model is shown in the RMS error colorbar for which each color represents the same
number of models. Therefore, it is observed that the low RMS error models also correspond
with the most probable ones.

<figure><img src="image_12.png" /><figcaption>Figure 11: SNMR at Mont Rigi - Posterior models distributions with their associated misfit and results of the QT inversion (dashed gray line). The black line in the RMS error graph represents the level of noise in the original dataset.</figcaption></figure>
456 5. Discussion

457 5.1. Impact of the number of sampled models in the prior model space

458
459
460
461
462
463
464
465
466
467
468
469
The optimum number of models in the prior model space sample is a tradeoff between
precision of the posterior models distributions and ease of computation: the larger the number
of samples, the larger the memory requirement and the harder the computations. We ran BEL1D
on the previous synthetic SNMR dataset with 10, 50, 100, 500, 1000, 5000, 10 000, 50 000 and
100 000 models in the prior samples. For each case, we ran BEL1D 10 times and aggregate the
results. It is observed that the RMS error distribution tends to stabilize above 1000 models,
suggesting that this value is optimal in our case (Figure 12), especially for the minimum and 5th
percentile values. This observation is confirmed by the evolution of the parameter distribution
characteristics (Figure 12). In those graphs, the mean values of the parameters are normalized
with the test values, and the standard deviations are normalized with the standard deviation in
the prior model space. We see that the parameter distributions are stabilized above 1000 models,
even if the values are not minimal.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Figure 12: Impact of the number of prior samples on the posterior distribution. The top graph presents RMS evolution, while the two bottom graphs show (from left to right) the normalized mean value and the normalized standard deviation.</figcaption></figure>

## 5.2. Impact of the chosen kernel density bandwidth on the posterior distributions

To assess the impact of the kernel density bandwidth on the posterior distribution, we ran 10 times 10 different values ranging from $10^{-5}$ and 1. In each case, 5000 models were sampled in the prior (Latin-hypercube sampler) and 5000 models were predicted, to ensure that no side effect is linked to other parameters. A similar analysis to the one performed in the previous subsection is then performed (Figure 13).

We observe that the bandwidth has a strong impact on the resulting distributions. Too low values result in erratic behaviors of the posterior distributions. On the other end, too large bandwidths result in the smoothing of the posterior distributions, as is observed in the evolution
of the standard deviation where this value increases sharply for the sensitive parameters.
483
Therefore, the optimal bandwidth is a value located in-between the two extremes, where the
484
reduction of uncertainty is significant for the sensitive parameters and the other parameters are
485
stable. This occurs around $10^{-2}$ in our case.
486

In our implementation of the kernel density estimator, the effect of a too small
487
bandwidth is avoided by verifying that enough points (at least 1% of the sampled prior) are
488
present in the direct surrounding of the reduced data ($d_{obs}^c \pm 3 \times bandwidth$). To perform the
489
test presented in this bandwidth analysis, this safeguarding measure was disabled. In the case
500
the default bandwidth is too small, then, the bandwidth is multiplied by 2 and the verification
501
is made again until the criterion is satisfied and only then, the kernel density estimation is
502
performed. We ran BEL1D with a bandwidth set to $10^{-5}$ and the safeguard algorithm enabled.
503
The obtained result is here similar to the one obtained with a suitable bandwidth as is shown by
504
the diamonds in the graph presenting the evolution of the normalized standard deviation in
505
Figure 13.
506
<figure><img src="image_16.png" /><figcaption>Figure 13: Impact of the bandwidth in the kernel density estimation on the posterior distribution. The top graph presents RMS evolution, while the two bottom graphs show (from left to right) the normalized mean value and the normalized standard deviation. This latter also presents results obtained with the safeguarding measure enabled (diamonds) for the most extreme case: a bandwidth of 1E-5.</figcaption></figure>

## 6. Conclusion

We propose an adaptation of the Bayesian Evidential Learning framework as a new tool for the unidimensional interpretation of geophysical data called BEL1D. This approach offers an alternative to stochastic inversions, as it requires fewer and simpler computations, thus lowering the CPU cost. This approach is also presented as an alternative to deterministic inversions as its behavior is superior both in terms of ease of implementation and completeness of the results (handling properly non-uniqueness, quantifying the uncertainty, etc.). We therefore developed this method theoretically, as a new general algorithm.
The framework does not require the use of regularization parameters or any other form of bias, but rather requires the proper definition of the prior model space. Contrary to other Bayesian methods, the use of a large prior model space (hence, introducing less bias in the posterior information) does not involve heavier computations to converge as the algorithm is always building in one step the posterior model space that corresponds to the observed data (see, 2.4 Step 4: generation of posterior distributions for the model parameters in CCA space using kernel density estimators). This method can therefore be considered unbiased if the prior model space is sufficiently large. BEL1D is thus a computationally efficient alternative to stochastic inversions as it only requires *n* forward models corresponding to the prior samples. Moreover, in BEL1D noise is considered via an imbedded theoretically based method, contrary to classically used high regularization parameters that lack of theoretical basis to be properly assessed and are therefore difficult to estimate.

When compared to McMC methods (such as DREAM), one of the key advantages of BEL1D is the tractability of all operations. When McMC classically applies random changes to the models and select them randomly, BEL1D constitutes a correlation between models parameters and the data. This means that one could seek for the origin of all the posterior models in the prior model space. This latter enables, for example, to efficiently sample multi-modal distribution or to understand the origin of odd models in the posterior distribution (see the numerical example), impossible to track back in McMC algorithms. The other key advantage is the simpler tuning to convergence of BEL1D compared to DREAM (McMC). In DREAM, it was required to change the dataspace to enable convergence towards reasonable (or even correct) posterior distributions. In BEL1D, as long as the prior is correctly assumed, the only parameter that the user needs to take care of is the number of models to sample in this space. As shown in the discussion, this parameter can be selected by observing the stabilization of the posterior distributions. In DREAM, the application without changes in the dataset
representation lead to either nearly no reduction of uncertainty (using the RMS as a proxy of the probability) or to a reduction towards level of uncertainty close to deterministic values but without the true model inside the posterior (using all data points as proxy of the probability). Therefore, applying this method to an unknown dataset might lead to erroneous interpretations, whereas BEL1D, even if producing slightly larger uncertainties should always predict consistent posterior spaces.

We demonstrated that BEL1D is efficient for the interpretation of SNMR data. We presented results for noisy datasets originating from a numerical model and field. In both cases, we validate the results from BEL1D, either against the test values from the model and McMC methods – DREAM (Vrugt, 2016) – for the numerical example or other geophysical experiments and classical state-of-the art inversion technics for the field case. On the other end, the careful analysis of the RMS error of the proposed models indicates that the range of uncertainty delivered by BEL1D is intrinsic to the non-unicity of the solution for geophysical inverse problems: many models explain the data to the same level.

In term of performance, the computation cost of BEL1D is directly proportional to the efficiency of the forward model and the number of samples in the prior model space. It must be noted that the prior samples being independent, the simulation of synthetic data can be fully parallelized. The performance of BEL1D can thus be easily estimated from the available computing facilities. The cost is therefore much smaller than other stochastic methods requiring ten to hundreds thousands of runs to converge towards a solution. Moreover, computation timing could be dramatically reduced by pre-field forward modelling operations for large geophysical campaigns where the prior would be similar. Prior realizations and the associated datasets would then be reused (PCA and CCA operations also), letting only the kernel density estimation, sampling and back-transformation operations to be applied after field acquisition,
559 leading to extremely rapid imaging. This is an extreme advantage compared to other stochastic
560 methods that rely on the data misfit (such as McMC) to sample the posterior.

561 This application, along with the pendulum case presented in supplementary material,
562 demonstrate the versatility of BEL1D and its applicability to multiple cases with very few
563 adaptations. Provided that a forward model already exists, the adaptation of the codes to any
564 other 1D geophysical method would be straightforward. For example, the forward operator
565 from HV-Inv (García-Jerez et al., 2016) could be used for the interpretation of H/V data. System
566 calls may also help further extend the applicability of BEL1D. This may be the case for the use
567 of all the 1D forward operators implemented in pyGIMLi (Rücker et al., 2017) – vertical
568 electrical sounding, frequency- and time-domain electromagnetic and magnetotelluric.

569 In its present version, BEL1D is dedicated to the 1D inversion of geophysical data in
570 layered Earth model where the contrast between layers is supposed to be sharp. In the future,
571 we plan to extend the method to smoothly varying systems accounting for a large number of
572 thin, correlated layers.

573 7. Acknowledgment

574 The authors would like to thank the contributors to the development of MRSmatlab, the
575 cornerstone of SNMR data interpretation in BEL1D. Hadrien Michel is a fellow funded by the
576 F.R.S.-FNRS.

577 8. Computer codes availability

578 The different MATLAB codes that have been used to produce the results presented in
579 this paper are available at github.com/hadrienmichel/BEL1D. They are released with this paper
580 under the BSD-2-Clause license. The codes consist in a series of graphical user interfaces for
581 the interpretation of SNMR data (BEL1DSNMR), among others. All the functions are
582 standalone (no specific toolbox required) in MATLAB but can run faster when the Parallel

583 Computing toolbox is available.

584 The main codes have been developed by H. Michel (ULiège, Belgium):

585 Hadrien MICHEL

586 Urban and Environmental Engineering Department

587 Allée de la découverte, 9

588 B-4000 Liège

589 Belgium

590 The forward models that is used is an Open Source MATLAB toolboxes for SNMR

591 (MRSMatlab: Müller-Petke et al., 2016).

592 In order to run the codes, a MATLAB license is required. Minimum requirements for

593 the computations are the ones for your MATLAB license, at the exception of RAM, for which

594 a minimum of 8 Go is recommended and more could be required for more complex

595 computations. If using the parallel capabilities of MATLAB (parallel computing toolbox), a

596 large number of cores is preferred and, globally, a 4-core (or more) processor is advised.

597 9. References

598 Aster, R., Borchers, B., Thurber, C., 2013. Parameters estimation and inverse problems, 2nd

599 ed. Academic Press.

600 Behroozmand, A.A., Auken, E., Fiandaca, G., Rejkjaer, S., 2016. Increasing the resolution and

601 the signal-to-noise ratio of magnetic resonance sounding data using a central loop
602 configuration. Geophysical Journal International 205, 243-256.

603 https://doi.org/10.1093/gji/ggw004

604 Behroozmand, A.A., Keating, K., Auken, E., 2015. A Review of the Principles and Applications

605 of the NMR Technique for Near-Surface Characterization. Surveys in Geophysics 36,

606 27-85. https://doi.org/10.1007/s10712-014-9304-0

607 Bobe, C., Van De Vijver, E., Keller, J., Hanssens, D., Van Meirvenne, M., De Smedt, P., 2019.

608 Probabilistic One-Dimensional Inversion of Frequency-Domain Electromagnetic Data

609 Using a Kalman Ensemble Generator. arXiv:1911.03275 [physics].
610 Bowman, A.W., Azzalini, A., 1997. Applied Smoothing Techniques for Data Analysis: The
611 Kernel Approach with S-Plus Illustrations, Oxford Statistical Science Series. Oxford
612 University Press, Oxford, New York.

613 Devroye, L., 1986. Non-uniform random variate generation. Springer, New York.

614 García-Jerez, A., Piña-Flores, J., Sánchez-Sesma, F.J., Luzón, F., Perton, M., 2016. A computer
615 code for forward calculation and inversion of the H/V spectral ratio under the diffuse
616 field assumption. Computers & Geosciences 97, 67–78.
617 https://doi.org/10.1016/j.cageo.2016.06.016

618 Gilson, M., Briers, P., Ruthy, I., Dassargues, A., 2017. Elsenborn - Langert - Dreiherrenwald.

619 Gorban, A.N., Kégl, B., Wunsch, D.C., Zinovyev, A.Y. (Eds.), 2008. Principal Manifolds for
620 Data Visualization and Dimension Reduction, Lecture Notes in Computational Science
621 and Enginee. Springer Berlin Heidelberg, Berlin, Heidelberg.
622 https://doi.org/10.1007/978-3-540-73750-6

623 Hanssens, D., Delefortrie, S., Bobe, C., Hermans, T., De Smedt, P., 2019. Improving the
624 reliability of soil EC-mapping: Robust apparent electrical conductivity (rECa)
625 estimation in ground-based frequency domain electromagnetics. Geoderma 337, 1155–
626 1163. https://doi.org/10.1016/j.geoderma.2018.11.030

627 Hermans, T., Lesparre, N., De Schepper, G., Robert, T., 2019. Bayesian evidential learning : a
628 field validation using push-pull tests. Hydrogeology journal in press.

629 Hermans, T., Nguyen, F., Klepikova, M., Dassargues, A., Caers, J., 2018. Uncertainty
630 Quantification of Medium-Term Heat Storage From Short-Term Geophysical
631 Experiments Using Bayesian Evidential Learning. Water Resources Research 54, 2931–
632 2948. https://doi.org/10.1002/2017WR022135

633 Hermans, T., Oware, E., Caers, J., 2016. Direct prediction of spatially and temporally varying
634 physical properties from time-lapse electrical resistance data. Water Resources
635 Research 52, 7262–7283. https://doi.org/10.1002/2016WR019126

636 Hertrich, M., Braun, M., Gunther, T., Green, A.G., Yaramanci, U., 2007. Surface Nuclear
637 Magnetic Resonance Tomography. IEEE Transactions on Geoscience and Remote
638 Sensing 45, 3752–3759. https://doi.org/10.1109/TGRS.2007.903829

639 Irving, J., Singha, K., 2010. Stochastic inversion of tracer test and electrical geophysical data
640 to estimate hydraulic conductivities. Water Resources Research 46.
641 https://doi.org/10.1029/2009WR008340

642 Jha, M.K., Kumar, S., Chowdhury, A., 2008. Vertical electrical sounding survey and resistivity
643 inversion using genetic algorithm optimization technique. Journal of Hydrology 359,
644 71–87. https://doi.org/10.1016/j.jhydrol.2008.06.018

645 Kemna, A., Kulessa, B., Vereecken, H., 2002. Imaging and characterisation of subsurface solute
646 transport using electrical resistivity tomography (ERT) and equivalent transport models.
647 Journal of Hydrology 267, 125–146.

648 Kemna, A., Nguyen, F., Gossen, S., 2007. On linear model uncertainty computation in electrical
649 imaging. Presented at the SIAM Conference on Mathematical and Computational Issues
650 in Geosciences, Santa Fe.

651 Krzanowski, W.J., 2000. Principles of multivariate analysis: A user's perspective, Revised
652 Edition. ed, Oxford Statistical Series. Oxford University Press, New York.

653 LaBrecque, D.J., Miletto, M., Daily, W., Ramirez, A., Owen, E., 1996. The effects of noise on
654 Occam's inversion of resistivity tomography data. Geophysics 61, 538–548.

655 Laloy, E., Hérault, R., Lee, J., Jacques, D., Linde, N., 2017. Inversion using a new low-
656 dimensional representation of complex binary geological media based on a deep neural
657 network. Advances in Water Resources 110, 387–405.
658 https://doi.org/10.1016/j.advwatres.2017.09.029
659 Laloy, E., Linde, N., Ruffino, C., Hérault, R., Gasso, G., Jacques, D., 2019. Gradient-based
deterministic inversion of geophysical data with generative adversarial networks: Is it
feasible? Computers & Geosciences 133, 104333.
660
661
662 https://doi.org/10.1016/j.cageo.2019.104333

663 Lawrence, N.D., 2012. A Unifying Probabilistic Perspective for Spectral Dimensionality
664 Reduction: Insights and New Models 30.

665 Lawyer, L.C., Bates, C., Rice, R., 2001. Geophysics in the affairs of mankind: a personalized
666 history of exploration geophysics - 2nd edition, Society of Exploration Geophysics. ed.
667 Tulsa, Oklahoma.

668 Lesparre, N., Nguyen, F., Kemna, A., Robert, T., Hermans, T., Daoudi, M., Flores-Orozco, A.,
669 2017. A new approach for time-lapse data weighting in electrical resistivity
tomography. Geophysics 82, E325–E333. https://doi.org/10.1190/geo2017-0024.1

671 Li, J., Lu, X., Farquharson, C.G., Hu, X., 2018. A finite-element time-domain forward solver
672 for electromagnetic methods with complex-shaped loop sources. GEOPHYSICS 83,
673 E117–E132. https://doi.org/10.1190/geo2017-0216.1

674 McKay, M.D., Beckman, R.J., Conover, W.J., 1979. A Comparison of Three Methods for
675 Selecting Values of Input Variables in the Analysis of Output from a Computer Code.
676 Technometrics 21, 239. https://doi.org/10.2307/1268522

677 Müller-Petke, M., Braun, M., Hertrich, M., Costabel, S., Walbrecker, J., 2016. MRSmatlab —
678 A software tool for processing, modeling, and inversion of magnetic resonance
sounding data. GEOPHYSICS 81, WB9–WB21. https://doi.org/10.1190/geo2015-
679
680 0461.1

681 Müller-Petke, M., Yaramanci, U., 2010. QT inversion — Comprehensive use of the complete
682 surface NMR data set. GEOPHYSICS 75, WA199–WA209.
683 https://doi.org/10.1190/1.3471523

684 Nguyen, F., Kemna, A., Robert, T., Hermans, T., 2016. Data-driven selection of the minimum-
685 gradient support parameter in time-lapse focused electric imaging. GEOPHYSICS 81,
A1–A5. https://doi.org/10.1190/geo2015-0226.1

687 Park, J., Yang, G., Satija, A., Scheidt, C., Caers, J., 2016. DGSA: A Matlab toolbox for
688 distance-based generalized sensitivity analysis of geoscientific computer experiments.
689 Computers & Geosciences 97, 15–29. https://doi.org/10.1016/j.cageo.2016.08.021

690 Phelps, G., Scheidt, C., Caers, J., 2018. Exploring viable geologic interpretations of gravity
691 models using distance-based global sensitivity analysis and kernel methods.
692 GEOPHYSICS 83, G79–G92. https://doi.org/10.1190/geo2017-0742.1

693 Rücker, C., Günther, T., Wagner, F.M., 2017. pyGIMLi: An open-source library for modelling
694 and inversion in geophysics. Computers & Geosciences 109, 106–123.
695 https://doi.org/10.1016/j.cageo.2017.07.011

696 Sambridge, M., 2002. Monte Carlo methods in geophysical inverse problems. Reviews of
697 Geophysics 40. https://doi.org/10.1029/2000RG000089

698 Satija, A., Caers, J., 2015. Direct forecasting of subsurface flow response from non-linear
dynamic data by linear least-squares in canonical functional principal component space.
Advances in Water Resources 77, 69–81.
700
701 https://doi.org/10.1016/j.advwatres.2015.01.002

702 Scheidt, C., Li, L., Caers, J., 2018. Quantifying Uncertainty in Subsurface Systems, Wiley-
703 Blackwell. ed.

704 Scheidt, C., Renard, P., Caers, J., 2015. Prediction-focused subsurface modeling: investigating
705 the need for accuracy in flow-based inverse modeling. Mathematical Geosciences 47,
173–191.
706
707 Socco, L.V., Foti, S., Boiero, D., 2010. Surface-wave analysis for building near-surface velocity

708 models — Established approaches and new perspectives. GEOPHYSICS 75, 75A83-

709 75A102. https://doi.org/10.1190/1.3479491

710 Trainor-Guitton, W., Hoversten, G.M., 2011. Stochastic inversion for electromagnetic

711 geophysics: Practical challenges and improving convergence efficiency.

712 GEOPHYSICS 76, F373–F386. https://doi.org/10.1190/geo2010-0223.1

713 Vrugt, J.A., 2016. Markov chain Monte Carlo simulation using the DREAM software package:

714 Theory, concepts, and MATLAB implementation. Environmental Modelling &

715 Software 75, 273–316. https://doi.org/10.1016/j.envsoft.2015.08.013

716 Wand, M.P., Jones, M.C., 1993. Comparison of Smoothing Parameterizations in Bivariate

717 Kernel Density Estimation. Journal of the American Statistical Association 88, 520.

718 https://doi.org/10.2307/2290332

719 Wastiaux, C., 2008. Les tourbières sont-elles des éponges régularisant l'écoulement ? Bulletin

720 de la Société géographique de Liège 50.

721 Wastiaux, C., Schumacker, R., 2003. Topographie de surface et de subsurface des zones

722 tourbeuses des réserves naturelles domaniales des Hautes-Fagnes (No. C60/1-2-3).

723 Université de Liège, Service de phytosociologie, flore et végétation des Hautes-Fagnes.

724 Yang, F., Ma, J., 2019. Deep-learning inversion: A next-generation seismic velocity model

725 building method. GEOPHYSICS 84, R583–R599. https://doi.org/10.1190/geo2018-

726 0249.1

727