# Sensitivity of constrained joint inversions to geological and petrophysical input data uncertainties with posterior geological analysis

Jérémie Giraud, Vitaliy Ogarko, Mark Lindsay, Evren Pakyuz-Charrier, Mark

Jessell, Roland Martin

► To cite this version:

Jérémie Giraud, Vitaliy Ogarko, Mark Lindsay, Evren Pakyuz-Charrier, Mark Jessell, et al. Sensitivity of constrained joint inversions to geological and petrophysical input data uncertainties with posterior geological analysis. Geophysical Journal International, 2019, 218 (1), pp.666-688. 10.1093/gji/ggz152. insu-04836399

HAL Id: insu-04836399

https://insu.hal.science/insu-04836399v1

Submitted on 13 Dec 2024

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
<figure><img src="image_2.png" /></figure>

# Sensitivity of constrained joint inversions to geological and petrophysical input data uncertainties with posterior geological analysis

Jérémie Giraud,<fnref n="1" /> Vitaliy Ogarko,<fnref n="2" />,<fnref n="3" /> Mark Lindsay,<fnref n="1" /> Evren Pakyuz-Charrier,<fnref n="1" /> Mark Jessell<fnref n="1" /> and Roland Martin<fnref n="4" />

Accepted 2019 March 24. Received 2019 February 14; in original form 2018 January 4

## SUMMARY

The integration of petrophysical data and probabilistic geological modelling in geophysical joint inversion is a powerful tool to solve exploration challenges. Models obtained from geologically and/or petrophysically constrained inversions are the result of complex interactions between correspondingly diverse data sets. Therefore, it is important to understand how non-geophysical input uncertainty impacts inverted models. In this paper, we propose to study the influence of uncertainty in geological and petrophysical measurements used to derive prior information and constraints onto geophysical inversion. Starting from geological field data from the Mansfield area (Victoria, Australia), we simulate low, medium and high uncertainty levels in geological measurements and petrophysical data, combined into a series of nine realistic case scenarios. This allows us to investigate the impact and propagation of uncertainty from non-geophysical measurements into geophysical inversion. We calculate misfit indicators and reconstruct lithological models *a posteriori* to analyse inversion results. We complement the examination of inverted models with topological analysis of lithological models in order to quantify the geological resemblance between the recovered and reference models. Our work reveals that the influence of uncertainty in geological measurements over the recovered lithological models is significantly stronger than it is for petrophysical data. Our posterior analysis indicates that intermediate petrophysical uncertainty provides optimum results.

**Key words:** Joint inversion; Persistence, memory, correlations, clustering; Statistical methods.

## 1 INTRODUCTION

The premise that modelling a physical system requires the use of all available information about it (Hempel & Oppenheim 1948; Cartwright & McMullin 1984; Nearing et al. 2016) justifies the integration of various disciplines in geophysical inverse modelling. It is motivated by the need to improve a model's geological reliability and to mitigate the inherent limitations of geophysical inversion. With minimum non-geophysical information, geophysical integration efforts focus on inverting several geophysical data sets jointly by enforcing structural similarities between inverted models (Haber & Oldenburg 1997; Gallardo & Meju 2003). Over the past years, the ever growing need to improve subsurface imaging through reduction and analysis of uncertainty increased the interest of the geoscientific community for multidisciplinary integration (Moorkamp 2017; Jessell et al. 2018). Numerous works have explored integration strategies to improve geophysical imaging for mineral, hydrocarbon, near surface and multiscale geophysics, as reviewed by Parsekian et al. (2015), Lelièvre & Farquharson (2016), Linde & Doetsch (2016) and Moorkamp et al. (2016), respectively.

More specifically, recent works highlight the advantages of integrating either petrophysical constraints (Lelièvre et al. 2012a; Carter-McAuslan et al. 2015; Kamm et al. 2015; Zhang & Revil 2015; Heincke et al. 2017; Sun & Li 2017) or geological information or data during inversion (Fullagar & Pears 2007; Lelièvre et al. 2012b, 2015; Lelièvre & Farquharson 2013; Revil et al. 2015; Scholl et al. 2016; Bijani et al. 2017; Lipari et al. 2017; Giraud et al. 2019). Alternatively, Brown et al. (2012), Zhou et al. (2014, 2016), Wiik et al. (2015), Guo et al.

<footnote n="1">¹Centre for Exploration Targeting, School of Earth Sciences, University of Western Australia, 35 Stirling Highway Crawley, Perth 6009, Western Australia, Australia. E-mail: jeremie.giraud@uwa.edu.au</footnote>

<footnote n="2">²International Centre for Radio Astronomy Research (ICRAR), University of Western Australia, Ken and Julie Michael Building, Crawley WA 6009.</footnote>

<footnote n="3">³ARC Centre of Excellence for all Sky Astrophysics in 3 Dimensions (ASTRO 3D).</footnote>

<footnote n="4">⁴CNRS, Géoscience environnement Toulouse, UMR CNRS 5563, Université Paul Sabatier, Observatoire Midi-Pyrénées, 14 Avenue Edouard Belin, 34000 Toulouse, France.</footnote>
(2017) and Giraud et al. (2019) propose the integration of structural information from geophysical or non-geophysical images in geophysical inversion to guide the inverse modelling process.

Concurrently, several works developed methodologies to characterize geodiversity (Lindsay et al. 2013a,b, 2014), the geological counterpart of biodiversity applied to geological modelling. Meanwhile, others focus on the study of geological uncertainty (Wellmann & Regenauer-Lieb 2012; Lark et al. 2013; Park et al. 2013; Kinkeldey et al. 2015; de la Varga & Wellmann 2016; Schneeberger et al. 2017; Schweizer et al. 2017). In particular, Pakyuz-Charrier et al. (2018a,b) produce probabilistic geological models (PGMs) from statistically uncertain inputs. These authors calculate PGMs through a technique called Monte-Carlo Uncertainty Estimator (MCUE). Their methodology relies on the sampling of a probability distribution representing geological measurements and their uncertainty. It is coupled to a geological modelling engine to produce a series of geologically plausible models. From this suite of geological models, the observation probability (i.e. the observed relative frequency) of each lithology is calculated in each model cell across the studied area, thus constituting a probabilistic geological model.

When available, such probabilistic geological models can be utilized to condition petrophysical constraints spatially, which has the potential to improve inversion results (Giraud et al. 2016a,b, 2017). Consequently, models calculated through constrained geophysical inversion encapsulate geological and petrophysical information. In such cases, as it is often observed in inverse modelling, integrated inversion suffers from nonlinearity. There has been little study of the propagation and sensitivity of inversion results to varying uncertainties in both geological and petrophysical input measurements even though it has been identified as a key factor for model evaluation in exploration scenarios (Bosch et al. 2010, 2015). The utilization of inaccurate information to derive petrophysical constraints has been investigated by Sun & Li (2015, 2017) and Carter-McAuslan et al. (2015). Nevertheless, the role of uncertainty in prior information (e.g. how data are scattered around the true value) is only partially addressed as published works mostly consider biased prior information such as systematic errors.

Previous works have investigated the effect of noise in geophysical data inversion (LaBrecque et al. 1996; Fernández-Martínez et al. 2014a,b) and its removal (Yuan et al. 2012; Pilkington & Shamsipour 2014). Understanding the respective influence of petrophysical and geological data uncertainty on the inversion process is elementary to reduce the risk of misinterpretation of results and to sound decision making. Integration approaches considering petrophysical and geological constraints simultaneously are a relatively recent development that present much potential, but that have not been thoroughly studied yet.

In this work, we develop a study that intends to alleviate the scarcity of studies focusing on the propagation of petrophysical and geological uncertainty in integrated inversion. We use a re-designed version of the 3-D inversion platform TOMOFAST3D (Martin et al. 2013, 2018; Ogarko et al. 2019) that integrates both statistical petrophysical constraints and probabilistic geological models. The primary objective of this article is to investigate and understand how uncertainty propagates from the geological and petrophysical input measurements to the recovered lithological model. We perform sensitivity analysis to study how variations in the statistics of petrophysical and geological measurements affect the inverted models and recovered lithologies a posteriori. We interpret inversion results qualitatively through visual inspection and quantitatively by calculating indicators quantifying the discrepancies between the true and inverted models in terms of physical properties (study of misfits) as well as recovered lithological models (topological analysis).

This paper is structured around five sections. In the methodology section (Section 2), we provide a summary of the different steps of the proposed workflow (Section 2.1). This is followed by the formulation of the inverse problem (Section 2.2), where we give a short introduction to our inversion platform. We then provide information about the probabilistic geological modelling procedure, the calculation of local petrophysical constraints and starting models (Section 2.3). Section 3 describes the metrics that assess both inversion quality and recovered lithological models. Section 4 introduces the geological context to the simulated case scenarios and practical information about the synthetic geophysical survey. The results section (Section 5) describes qualitative interpretation and examination of inversion results. We subsequently validate our observations using a series of posterior indicators. The discussion (Section 6) considers the implications of our findings for future work and its practical applications. Section 7 provides the conclusions. Section 8 provides information allowing reproducibility of the work presented in this article.

# 2 INVERSE PROBLEM FORMULATION AND INTEGRATION OF CONSTRAINTS

## 2.1 Integration procedure summary

The main steps of the modelling procedure are summarized in Fig. 1. The inversion algorithm is similar to Giraud et al. (2017), which we extend here to large-scale 3-D problems and complement with the posterior recovery of lithological models and topological analysis. A prerequisite to inversion is to translate petrophysical and geological data into information that can be used in the inverse modelling process. The first step is to calculate the PGM using MCUE and to derive a mixture model (e.g. a weighted sum of probability density functions) representative of the petrophysical measurements. The PGM and the mixture model are then combined to derive local petrophysical constraints and to calculate starting models for inversion (step 2 in Fig. 1). The next step of the workflow is to perform constrained geophysical inversion (step 3). The posterior analysis in step 4 comprises a lithological reconstruction procedure that allows the recovery of a lithological model and subsequent topological analysis. The calculation of a series of indicators provides metrics for the assessment of inversion results in terms of geological plausibility and misfit.
<figure><img src="image_3.png" /><figcaption>Figure 1. Modelling workflow summary.</figcaption></figure>

## 2.2 Cost function and optimization process

### 2.2.1 Problem formulation

The optimized objective function $\Theta$ is derived from the formulation of a density function that relates to the ‘degree of knowledge that we have about the values of the parameters of our system’ (Tarantola & Valette 1982). In our case, it can be expressed as a product distribution of the form

$$ \Theta(d, m) = \Theta_d(d, m) \Theta_m(m) \Theta_c(m), \qquad (1) $$

where subscript $d$ relates to the geophysical measurements $d$ to be inverted, subscript $m$ relates to the inverted model $m$ and subscript $c$ to the petrophysical constraints. In this equation, $\Theta_d(d, m)$ is the density function corresponding to geophysical data misfit, $\Theta_m(m)$ is the prior density function on model $m$ and $\Theta_c(m)$ is the probability density relating to constraints derived from petrophysical information.

Our assumption is that density functions $\Theta_d(d, m)$ and $\Theta_m(m)$ can be expressed as

$$ \left\{ \begin{array}{l} \Theta_d(d, m) = A \exp \left[ -(d - g(m))^T C_d^{-1} (d - g(m)) \right], \quad A \in \mathbb{R}^+ \setminus \{0\} \\ \Theta_m(m) = B \exp \left[ -(m - m_p)^T C_m^{-1} (m - m_p) \right], \quad B \in \mathbb{R}^+ \setminus \{0\}, \end{array} \right. \qquad (2) $$

with

$$ g(m) = \begin{bmatrix} g_G(m^G) \\ g_M(m^M) \end{bmatrix}, \quad m = \begin{bmatrix} m^G \\ m^M \end{bmatrix}^T, \quad d = \begin{bmatrix} d_G \\ d_M \end{bmatrix}, \quad C_d = \begin{bmatrix} C_d^G & 0 \\ 0 & C_d^M \end{bmatrix}, \quad C_m = \begin{bmatrix} C_m^G & 0 \\ 0 & C_m^M \end{bmatrix}, \qquad (3) $$

where superscript $T$ denotes the transpose operator. Superscripts and subscripts $G$ and $M$ refer to the gravity and magnetic problems, respectively. In this case, $\Theta_c(m)$ is expressed using a mixture model where each subpopulation corresponds to a specific rock type. $A$ and $B$ are the normal distributions’ normalization constants, while $C_d$ and $C_m$ are covariance matrices corresponding to data and model weighting, respectively, and $m_p$ is the prior model.

Substituting eqs (1) and (2), we obtain

$$ \Theta(d, m) = AB \exp \left[ -(d - g(m))^T C_d^{-1} (d - g(m)) - (m - m_p)^T C_m^{-1} (m - m_p) \right] \Theta_c(m). \qquad (4) $$
By considering the logarithm of eq. (4), it is straightforward to show that maximizing $Θ(d, m)$ is equivalent to minimizing the following cost function $\theta$:

$$ \theta(d, m) = \underbrace{(d - g(m))^T C_d^{-1} (d - g(m))}_{\text{data term } \phi_d(m)} + \underbrace{(m - m_p)^T C_m^{-1} (m - m_p)}_{\text{model term } \phi_m(m)} + \underbrace{\log (\Theta_c(m)^{-1})}_{\text{constraint } \phi_c(m)}. \quad (5) $$

Generally speaking, the constraint term $\phi_c(m)$ can be utilized to incorporate prior information in the inverse problem and to enforce constraints. It constitutes the link between the two data sets jointly inverted as each element is a function of both the density contrast and magnetic susceptibility. In this method, $\phi_c$ encapsulates local information derived from probabilistic geological modelling combined with the statistics of the physical measurements of the lithologies sampled in the studied area. Calculation of this term is detailed in Section 2.3.2. We assume that $C_d$ and $C_m$ are diagonal covariance matrices for the remainder of this paper. In this work, we weight all the geophysical data points equally. We set $C_d$ as the identity matrix to simplify the estimation of the parameters weighting the different terms. Note that in the general case, values of $C_d$ may vary between the data points accordingly with area, acquisition or measurement-specific information to account for the features of the data set. It remains necessary to estimate $C_m$, the parameters controlling the weight of the petrophysical constraints in inversion, and to determine the relative weight between the gravity and magnetic problems. This procedure is discussed in Section 4.3 as applied to our data set.

### 2.2.2 Inversion algorithm

The form of $\theta(m)$ as per eq. (5) allows us to solve the inverse problem in a least-squares (LSQR) sense using the LSQR algorithm (Paige & Saunders 1982; Chapman & Pratt 1992; Pratt & Chapman 1992; Gerhard Pratt et al. 1998; Martin et al. 2013).

We extend the parallel code TOMOFAST3D for the purpose of this work (Martin et al. 2013, 2018; Ogarko et al. 2019). This extended implementation, ‘Tomofast-x’, follows the object-oriented FORTRAN 2008 standard. The design of Tomofast-x utilizes classes derived to account for the mathematics of the problem. This permits to reduce software complexity, thereby facilitating the addition of new functionalities (see Hammond et al. 2014, and references therein).

## 2.3 Geological conditioning of petrophysical constraints

### 2.3.1 Probabilistic geological modelling

Several studies showed that multiple sources of uncertainty can impact geological modelling, and that quantitative uncertainty estimation of individual geological models is difficult to obtain (Bond et al. 2007, 2015; Alcalde et al. 2017). To mitigate this, Pakyuz-Charrier et al. (2018a,b,c) extend previous works (Jessell et al. 2010; Lindsay et al. 2012, 2013b; Wellmann & Regenauer-Lieb 2012) to generate a collection of geological models reflecting the range of geologically possible models. In this work, we perform probabilistic geological modelling using outputs from MCUE as detailed in Pakyuz-Charrier et al. (2018b) as part of the prerequisite to inversion (step 1 as per Fig. 1).

MCUE propagates input measurement uncertainty to the modelling geological structures via sampling of probability distributions using a Monte-Carlo approach. MCUE perturbs a reference model under the assumption that uncertainty on measurement position can be modelled using a normal distribution and that measured orientation data used to build the model can be modelled using spherical statistics (e.g. foliations through dip and strike converted into a vector in 3-D). This is achieved using the von Mises-Fisher distribution (Davis 2002). In 3-D, the distribution is given by Mardia & Jupp (2008):

$$ vMF(x|y, \kappa) = \frac{\kappa}{4\pi \sinh(\kappa)} e^{\kappa y^T x}, \quad \kappa \in \mathbb{R}^+, \quad (6) $$

where $x$ is the orientation data vector and $y$ is a unit vector representing the mean measurement direction. The term $\kappa$ is called the *concentration parameter*. It characterizes the spread in orientation measurements and is analogous to the inverse of the more familiar standard deviation for normal distributions. In principle, it is possible to assign an individual uncertainty estimate to each measurement depth- or location-dependent uncertainty (e.g. for borehole data or for seismic horizons). In the work presented here, we assume independent surface measurements that present similar levels of uncertainty.

The geological plausibility of model realizations is enforced through application of topological rules and plausibility filters that force the model to honour age relationships and the stratigraphic column. The resulting set of models—typically several hundreds of models—is then combined into a PGM representing the observation probability of all modelled lithologies in each cell. This applies to all lithologies with the exception of the basement, which is assigned in the absence of other units under the hypothesis that it fills space when no other lithology is present.

The results from MCUE are used to condition petrophysical constraints ($\Theta_c(m)$ term in eqs 1–5) and to compute starting models and is detailed in Section 2.3.2 next.
### 2.3.2 Local petrophysical constraints

Petrophysical constraints are applied to inversion through the minimization of the constraint term $\phi_c$ (5) simultaneously to the other terms. We follow concepts introduced and used by Paasche & Tronicke (2007), Doetsch et al. (2010), Lelièvre et al. (2012a) and Sun & Li (2013) to maximize the statistical closeness between inverted properties and measured petrophysical data. We condition the petrophysical constraints geologically to derive local constraints following the procedure introduced in Giraud et al. (2017).

Assuming that the physical properties are normally distributed for each lithology, petrophysical measurements can be represented using a mixture model. Therefore, we can use such statistical description to define a probability distribution to be used to constrain inversion. Let the statistics of petrophysical measurements be represented by a mixture distribution of the form (eqs 7 and 8)

$$p(\mathbf{m}_i) = \sum_{k=1}^{n_1} p_k(\mathbf{m}_i), \qquad (7)$$

$$\text{with} \\ p_k(\mathbf{m}_i) = \omega_k N(\mathbf{m}_i | \mu_k, \sigma_k), \qquad (8)$$

where $n_1$ is the number of lithologies observed in the petrophysical measurements and $k$ denotes the index of the corresponding lithology. The kernels $N(m_i|\mu_k, \sigma_k)$ of the mixture are normal distributions. The normal distribution for the $k$th lithology is characterized by a mean value vector for the considered property, $\mu_k$, and the associated covariance matrix, $\sigma_k$. The relative weights of the corresponding lithologies, $\omega_k = 1..n_1$, are non-negative weights, or mixing coefficients that sum to 1. Note that in our implementation, $\sigma_k$ is a full $2 \times 2$ matrix, thereby accounting for the statistical correlation between the petrophysical properties of rocks.

Introducing the weight $\rho \in \mathbb{R}^+$, we define

$$\phi_c(\mathbf{m}) = \rho \sum_{i=1}^{n_m} (\log p(\mathbf{m}_i) - \log \max p(\mathbf{m}_i))^2 \qquad (9)$$

for optimization in the LSQR problem to be optimized as per eq. (5). The term inside brackets can be interpreted as the summed log-likelihood of the distribution representative of the petrophysical measurements multiplied by a positive scalar weighting the petrophysical constraints term in eq. (5). From eqs (7–9), it becomes clear that such petrophysical constraints constitute soft constraints and favour model changes towards the most likely configuration as indicated by the expression of $p(\mathbf{m})$.

In real-world scenarios, the parameters of $N(m|\mu_k, \sigma_k)$, $\sigma_k$ and $\mu_k$ can be estimated from petrophysical data using an expectation maximization algorithm (e.g. McLachlan & Peel 2000). A normal distribution is commonly expected to describe the physical measurements of rocks. It is clear that the assumption that the statistics of measured petrophysical properties can be represented or approximated using such distribution is not always valid (e.g. bimodal susceptibility distribution of certain rocks), and that the choice of distribution remains case dependent.

The conditioning of petrophysical constraints by geological modelling employs calculation of $\phi_c(\mathbf{m})$ using observation probabilities extracted from the PGM to weight elements of the mixture model defined in eq. (9). We introduce a matrix $\psi$ of dimensions $n_m \times n_1$ containing the observation probabilities of the different rock units. Substituting $\omega_k$ as per eq. (8) with the probability of the respective lithology in each cell of the model provides a local mixture model depending on both geological modelling results and petrophysical measurements. We can rewrite $p_k(\mathbf{m}_i)$ as

$$p_k(\mathbf{m}_i) = \psi_{k,i} N(\mathbf{m}_i | \mu_k, \sigma_k). \qquad (10)$$

This process, which corresponds to the second stage of the workflow (Fig. 1), is illustrated in Fig. 2.

In addition to the calculation of $\phi_c(\mathbf{m})$, the PGM is useful to calculate starting models. The starting model $\mathbf{m}_0$ is set as the average model from the PGM populated accordingly with petrophysical measurements. It is calculated as

$$\mathbf{m}_0 = \psi \mu. \qquad (11)$$

The consequences of the formulation of the petrophysical constraints in eq. (9) and of the starting model in eq. (10) for inversion are the following. Theoretically, inversion will be strongly influenced by the petrophysical constraints in model cells where geology is well determined (e.g. one observation probability dominates the others) and the locally weighted mixture shows well-separated lithologies. In contrast, local petrophysical constraints are similar to global petrophysical constraints in model cells where geology remains undetermined after MCUE. Also note that model updates during inversion are largely driven by geophysical data where the locally weighted mixture model does not exhibit a sufficiently clear distinction between the different lithologies.

## 3 QUALITY INDICATORS

Posterior analysis of inversion results constitutes the fourth and final step of the workflow (Fig. 1). We base our interpretation of inverted models on global and local indicators. All the inversions we present achieve a similar geophysical data misfit (i.e. the root-mean-square error, RMSE), indicating that little insight is gained from comparing models with this particular metric (see Section 5.2.1). We observe that
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Geological conditioning of petrophysical constraints. Modified from Giraud *et al.* (2017).</figcaption></figure>

incorporation of prior information with varying degrees of uncertainty has limited impact on final geophysical data misfit values. Our results are in agreement with Gallardo & Meju (2004, 2007, 2011), Abubakar *et al.* (2012), Gallardo *et al.* (2012), Gao *et al.* (2012), Jardani *et al.* (2013), Molodtsov *et al.* (2013), Moorkamp *et al.* (2013), Juhojuntti & Kamm (2015), Rittgers *et al.* (2016), Sun & Li (2016a, 2017), Demirel & Candansayar (2017) and Giraud *et al.* (2017), who do not observe dramatic data misfit improvements through data integration, and point out that improvements mostly occur in model space.

## 3.1 Global indicators: model misfit and lithological resemblance

Models are assessed using indicators derived from the difference between reference and inverted models and characterize the resemblance between causative bodies and retrieved models. Invoking Laplace's first law of errors (Laplace 1774; Wilson 1923; Stigler 1986), the first indicator we calculate is equivalent to the mean absolute deviation around the reference model, or mean absolute model misfit (MAMM). It is less sensitive to isolated outliers than the more common root-mean-square indicators. Let us express this indicator as follows:

$$ \varphi_m (\mathbf{m}) = \frac{1}{n_m} \sum_{i=1}^{n_m} |m_i^{\text{ref}} - m_i^{\text{inv}}|, \qquad (12) $$

where $m_i^{\text{ref}}$ and $m_i^{\text{inv}}$ correspond to the $i$th cell of the inverted reference and inverted models, respectively.

We also assess inverted models by reconstructing lithologies from inverted physical property models. Lithologies are recovered, for each model-cell assigned with $m_i$, from a membership analysis using eq. (8) similar to Doetsch *et al.* (2010) and Sun & Li (2012), which we restrict locally to lithologies characterized by $\psi_{k, i} > 0$. This allows us to evaluate the lithological resemblance between recovered lithological model $l^{\text{inv}}$ and the reference lithological model $l^{\text{ref}}$ by calculating the proportion of recovered cells that are in agreement. We define what we term the lithological resemblance $r$ as follows:

$$ r (l^{\text{inv}}, l^{\text{ref}}) = \frac{1}{n_m} \sum_{i=1}^{n_m} 1_{l^{\text{ref}}} (l^{\text{inv}}), \qquad (13) $$

where $1_{l^{\text{ref}}}$ symbolizes the indicator function for $l_i^{\text{ref}}$.

## 3.2 Local indicator using the topology of inverted models

The set of cells assigned with recovered lithologies constitute a lithological model characterized by discrete values representing rock types, which can be interpreted geologically. Topology records the geological relationships between adjacent objects within a model, and can be used to represent that model (Perrin & Rainaud 2013). It can also be used as 'tool for quantifying differences between geological models in uncertainty analyses' (Thiele *et al.* 2016b). Topological analysis allows us to interpret and compare results in a way that also accounts for
changes in geological relationships regardless of the scale considered. It is, therefore, well suited to detect local differences between models, and complements indicators reflecting bulk properties (i.e. values averaged from the whole model). This can be useful to identify areas of the recovered models that are geologically implausible and to highlight poorly constrained lithologies. Such analysis of inverted models complements works of Doetsch et al. (2010), Sun & Li (2012), Carter-McAuslan et al. (2015), Martinez & Li (2015), Sun & Li (2015, 2017), Li & Sun (2016) and Melo et al. (2017), who infer lithologies from inverted models but whose study does not extend to the quantitative geological assessment of recovered models. One of the key advantages of topological analysis over qualitative interpretation is that it presents a systematic way to analyse results and to detect features that cannot be measured using mean model indicators and that may be challenging to impossible to detect through visual inspection of 3-D models (Pellerin et al. 2017).

In this case, we perform the topological analysis of recovered lithological models to compare basic geological features of models recovered from inversion results. Following the nomenclature proposed by Thiele et al. (2016b), we use ‘lithological topology’. It is defined in the framework introduced by Burns (1988), who addresses the relationship between adjacent rock volumes (in our case, cells assigned with a given lithology). The different kinds of topological relationships, for instance, ‘contains’, ‘is inside’, ‘meets’, etc., are characterized using the Egenhofer relations (Egenhofer & Herring 1990), generalized to 3-D by Zlatanova (2000).

In the workflow presented here, we limit our investigations to the adjacency of units (i.e. ‘unit A is in contact with unit B’), which is the most common Egenhofer relation used in geology (Thiele et al. 2016b). We represent adjacency following the matrix representation of Godsil & Royle (2001), which can be summarized as follows. For a model consisting of n lithologies, the adjacency matrix can be reduced to a positive strictly lower triangular n_l × n_l matrix M, where element M_{i > j, j} are non-zero if and only if a contact between lithologies i and j is observed.

We attach to each element M_{i > j, j} of the matrix the number of occurrences of the respective contacts to make the adjacency matrix more suitable to the comparison of lithological models. Let the resulting adjacency matrix elements be

$$M_{i > j, j} = n_{ij}, \quad M_{i \le j, j} = 0, \qquad (14)$$

where n_{ij} is the number of contacts between lithology i and j, for node M_{i > j, j}; n_{ij} = 0 indicates that lithology i and j are disjoint.

The synthetic case study uses adjacency matrices to compare recovered lithological models with the reference model. Real-world studies can benefit from comparison of geological scenarios or departures from specific geological models (Perrin & Rainaud 2013; Thiele et al. 2016a; Pellerin et al. 2017). Besides reflecting the global features of a lithological model, the adjacency matrix M is also a local indicator that can change dramatically with small, local alterations in the lithological model, possibly reflecting a change in mathematical properties and geological meaning.

Deviations from the reference lithological model are calculated from the difference between the adjacency matrices of the recovered and reference model and normalized by the highest number of contacts observed in the reference model. Let us define the matrix M^{rel} as follows:

$$M^{\text{rel}} = \frac{1}{\max M^{\text{rel}}} (M^{\text{inv}} - M^{\text{ref}}), \qquad (15)$$

where M^{inv} and M^{ref} represent the adjacency matrices calculated from inversion results and the reference model, respectively.

From eq. (15), M^{rel}_{i,j} > 0 indicates that the recovered lithologies overestimate the contact surface area between lithology i and lithology j. On the contrary, M^{rel}_{i,j} < 0 indicates that the recovered lithologies underestimate the contact surface area between lithology i and lithology j.

# 4 SIMULATED CASE STUDY AND UNCERTAINTY SCENARIOS

## 4.1 Geological context

The unperturbed (or reference) geological model was constructed from contact data and surface orientation (i.e. contact points and foliations, respectively) collected in the Mansfield area (Victoria, Australia) in Geomodeller® (Lajaunie et al. 1997; Calcagno et al. 2008). It presents a Carboniferous sedimentary syncline oriented N170, abutting a faulted contact with a Silurian–Devonian folded basement. Structural complexity was artificially increased to define a problem that is more challenging to the inversion algorithm through the addition of a north–south oriented fault and an ultramafic intrusion to the south. Information used to derive the original model is detailed in the GeoModeller User Manual, Tutorial case study H (Mansfield). Fig. 3(a) shows the unperturbed geological model. The base structural model was made available online by Pakyuz-Charrier (2018; see Section 8 for details about availability of data and models shown here). The utility of this model for testing purposes has been shown by Giraud et al. (2017, 2019) and Pakyuz-Charrier et al. (2018b).

## 4.2 Simulated physical properties

We populate the reference geological model by assigning low density contrast and limited magnetic susceptibility to basin fill (lithologies 3, 5 and 6). Using values from the literature (Hatfield et al. 2002; Hunt et al. 1995; Airo 2005; Clark & Emerson 2006), we assigned higher density contrast and magnetic susceptibilities to ultramafic rocks (lithology 1), diorite (lithology 2) and dolerite (lithology 4). Modelled lithologies and physical properties are provided in Table 1. The reference density and magnetic susceptibility models are shown in Figs 3(b) and (c),
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Reference geological model, reference property model: (b) density contrast (c) magnetic susceptibility, with the corresponding (d) gravity anomaly and (e) magnetic field anomaly. The dashed line square in (d) and (e) shows the outline of the subset of the model used in the first stage of hyperparameter estimation.</figcaption></figure>

<table><caption>**Table 1.** Stratigraphic column showing geological topological relationships and average physical properties. Lithologies are indexed from 1 through 6 by decreasing age. Lithology 1 is the oldest and 6 is the most recent.</caption><thead><tr><th>Lithology index</th><th>Geological relation</th><th>Geological unit</th><th>Density contrast (kg m<sup>-3</sup>)</th><th>Magnetic susceptibility (SI)</th></tr></thead><tbody><tr><td>6</td><td>Sedimentary</td><td>Basin fill 3</td><td>0</td><td>0</td></tr><tr><td>5</td><td>Sedimentary</td><td>Basin fill 2</td><td>110</td><td>0.01</td></tr><tr><td>4</td><td>Intrusive</td><td>Dolerite</td><td>300</td><td>0.025</td></tr><tr><td>3</td><td>Sedimentary</td><td>Basin fill 1</td><td>110</td><td>0.01</td></tr><tr><td>2</td><td>Intrusive</td><td>Diorite</td><td>170</td><td>0.05</td></tr><tr><td>1</td><td>Basement</td><td>Ultramafic rocks</td><td>240</td><td>0.08</td></tr></tbody></table>

respectively. A reference density value of 2670 kg m⁻³ was used.

In real-world scenarios, different lithologies encountered in basins may share similar characteristics in the petrophysical domain, such as low density contrasts and magnetic susceptibilities, and in some cases cannot be differentiated in a density contrast—magnetic susceptibility crossplot. To account for this, we assign the same properties to lithologies 3 and 5, which become undistinguishable.

## 4.3 Geophysical survey

### 4.3.1 Synthetic geophysical data and forward modelling

Magnetic and gravity data were computed to simulate and gravity and magnetic ground surveys. We inverted for the total magnetic field anomaly and the vertical component of the complete Bouguer anomaly. Magnetic data are simulated following the same approach as Bhattacharyya (1964). Gravity data are simulated following Boulanger & Chouteau (2001). We model a magnetic field strength equal to 57 700 nT, which approximates the International Geomagnetic Reference Field for the locality of Mansfield (Victoria, Australia).
<table><caption>**Table 2.** Input uncertainty for geological input data (first block) and physical properties (second block). Cases (a), (b) and (c) correspond to cases scenarios where 95 per cent of the orientation data are contained within 52, 37 and 30 deg aperture cones, respectively. For petrophysics, cases α, β and γ characterize different levels of spread of the measurements around the mean values. Chosen standard deviations for case α (most uncertain case) in terms of magnetic susceptibility are of the same order of magnetic as Sanger & Glen (2003), larger than Törnberg & Sturkell (2005) for density contrast, and conform to Barlow (2004) for both density contrast and magnetic susceptibility. Chosen concentration parameters are in agreement with metrological studies (Bond 2015; Allmendinger et al. 2017; Novakova & Pavlis 2017).</caption><thead><tr><th colspan="3">Geological input data uncertainty</th></tr><tr><th>Uncertainty level</th><th>Orientation data uncertainty<br>(Concentration parameter)</th><th>Contact position uncertainty<br>(Standard deviation, in m)</th></tr></thead><tbody><tr><td>a—most uncertain case</td><td>80</td><td>50</td></tr><tr><td>b—intermediate case</td><td>20</td><td>50</td></tr><tr><td>c—least uncertain case</td><td>25</td><td>25</td></tr><tr><th colspan="3">Physical property uncertainty</th></tr><tr><th>Uncertainty level</th><th>Standard deviation of density contrast (kg m<sup>-3</sup>)</th><th>Standard deviation of magnetic susceptibility (SI)</th></tr><tr><td>α—most uncertain case</td><td>60</td><td>0.02</td></tr><tr><td>β—intermediate case</td><td>45</td><td>0.015</td></tr><tr><td>γ—least uncertain case</td><td>33.75</td><td>0.0113</td></tr></tbody></table>

We add zero-mean normally distributed random noise with standard deviation equal $\sigma_d^G = 0.3$ mGal (e.g. 2.25 per cent of the average amplitude of the data) and to $\sigma_d^M = 10$ nT (e.g. 2.5 per cent of the average amplitude of the data) to gravity and magnetic measurements, respectively, to simulate noise contamination of the data. The synthetic geophysical data we invert are shown in Figs 3(d) and (e) for gravity and magnetic data, respectively. The inverse integrated sensitivities technique of Li & Oldenburg (2000) and Portniaguine & Zhdanov (2002) are used to balance the decreasing sensitivity of potential field data with depth.

The study area is a volume of $16.64 \times 16.64 \times 2.88 \text{ km}^3$ discretized into $128 \times 128 \times 32$ cells (dimensions $130 \times 130 \times 90 \text{ m}^3$), making up a total of 524 288 elements. To avoid dispersion effects, we sampled geophysical data on a $128 \times 128$ grid, for a total of 16 384 measurement points per geophysical data set.

### 4.3.2 Hyperparameter estimation

The values defining the weights of the different terms in the cost function as per eq. (2) through eq. (5) constitute hyperparameters of the inverse problem. We estimate them using the L-curve principle (Hansen & O'Leary 1993; Hansen & Johnston 2001; Santos & Bassrei 2007), seeking geophysical data misfit superior or equal to noise level. We determine $C_m$ and $\rho$ in a two-stage heuristic process that partially avoids the limitations of L-curve analyses in 3-D constrained joint inversion scenarios due to high computation requirements (see Bijani et al. 2017) or the difficulty to visualize hypersurfaces. The first stage involves the estimation of optimum weights for single-physics inversion on a 3-D subset of the full model. It is comprised of $48 \times 48 \times 32$ cells with the corresponding $48 \times 48$ surface data points. The subset we chose is located in the central part of the model. The reference model is made fairly complex in this location by inclusion of an ultramafic intrusion and faulted units. We assume that the corresponding data are representative of the full geophysical data set. The subset comprises approximately 14 per cent of the total volume and data points, thereby dividing the computation time by a factor superior to 10. This allows us to run hundreds of constrained inversions. The outline of the model is shown in top view in Figs 3(d) and (e).

Manual tuning of $C_m$ and $\rho$ for single-physics inversions provides an initial estimate of the optimal parameter values. We refine our search by sampling values spaced at regular intervals in logarithmic scale around these initial estimates. This led us to perform series of 400 and 1032 separate gravity and magnetic inversions, from which we generate the corresponding L-surfaces for the identification of the optimum weights (the corresponding L-surfaces are available in Appendix A1, Fig. A1).

The second stage of hyperparameter estimation pertains to joint inversion. The values obtained for single domain inversion are transferred to the joint inversion of the subset of the full model. The relative weight between the gravity and magnetic problems is obtained through L-curve analysis of 50 different weights. A series of 625 inversions sampling values of $C_m$ and $\rho$ from single domain inversion is then performed to determine optimum values. Finally, fine-tuning of the resulting values using the full model allowed us to obtain the set of optimum hyperparameters that we used to run the inversions shown in this paper.

## 4.4 Case scenario evaluation

In real-world scenarios, statistical geological and petrophysical models derived from geological and petrophysical measurements are dependent upon the quality of data and on the geology of the area. The evaluation of the sensitivity of our inverse modelling procedure to varying uncertainty levels in geological and petrophysical data is performed through the simulation of a range of possible case scenarios. To this end, we simulate three base uncertainty levels for both geological and petrophysical data to simulate extreme and average cases representative of real scenarios. Note that we simulate measurements and uncertainty levels using values reported in the literature (see references provided in the caption of Table 2). We assume that the statistics of physical properties are derived from a representative population in our idealized experiment. Likewise, we assume that measurement uncertainty on orientation data or contact positions are based on outcrop and/or borehole condition
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Simulated case scenario classification matrix representing the tested combinations of uncertainty levels on geological and petrophysical input for inversion as per Table 2. Each row corresponds to a given geological input measurement uncertainty scenario, while each column corresponds to a given petrophysical uncertainty case.</figcaption></figure>

and measuring instrument specifications. Table 2 next gives the tested uncertainty levels on geological and petrophysical measurements, assuming that density contrast and magnetic susceptibility are not correlated (e.g. crosscorrelation equal to zero).

Combining the base uncertainty levels as per Table 2, we obtain a total of nine cases representative of the range of simulated case scenarios (Fig. 4 next). The comparison and analysis of results obtained using these nine cases (Fig. 4) constitute our sensitivity analysis of integrated inversion to geological and petrophysical uncertainty.

Fig. 5 shows the observation probability of the most probable lithology in each model cell, calculated, in each case, from a suite of 300 plausible geological models. It shows the amount of information brought by geological measurements across the different parts of the model. In Fig. 5, one can notice that for case (a) interfaces between geological units are not distinguishable in large portions of the model (darker regions), while case (b) seems less affected and case (c) is the best determined. This clearly illustrates the propagation of input measurement uncertainty to the PGM.

Fig. 6 next shows the resulting mixture model as per eq. (8), where all lithologies are equally weighted. As can be seen in Fig. 6, although the geological model comprises six distinct units, lithologies 3 and 5 present exactly the same density contrast and magnetic susceptibility (Table 1). This simulates cases where geological units cannot be differentiated by the available petrophysical data.

# 5 SENSITIVITY ANALYSIS TO UNCERTAINTY IN PETROPHYSICAL AND GEOLOGICAL INPUT DATA

## 5.1 Inverted models

In this section, we analyse the inverted models qualitatively and examine the differences between reference and inverted models in terms of density contrast for (a, β), (b, β) and (c, β) (Fig. 7) and magnetic susceptibility for (b, α), (b, β) and (b, γ) (Fig. 8). Visual inspection of the recovered lithological models corresponding to the cases shown in Figs 7 and 8 are shown in Figs 9(a) and (b), respectively. The complete set of results for all cases (Fig. 4) in terms of inverted models, differences with the true model and recovered lithological models is given in Appendix A2 (Fig. A2 through A6).

From Figs 7(a) and (b), differences between cases are small. However, visual comparison of inversions (b, α), (b, β) and (b, γ) in Fig. 7(a) indicates that decreasing uncertainty in petrophysical data increases the contrasts in inverted physical properties. This is because the different lithologies are increasingly well differentiated in the crossplot from case α to γ (Fig. 6). To a lesser extent, decreasing petrophysical uncertainty also modifies the geometry of the contacts between structures in the model (small-scale changes in lithological models shown by arrows in Fig. 9a). In particular, changes are more notable in the central part of the model where Fig. 3 respectively show higher geological complexity and uncertainty.

Although the same mixture model is used to define petrophysical constraints for cases (a, β), (b, β) and (c, β), results shown in Fig. 8(a) exhibit only a slight increase in sharpness from (a, β) to (c, β). Visually, it is notable in Fig. 8(a) that decreasing geological input uncertainty results in increased inverted model complexity associated with lower discrepancies with the true model (Fig. 8b). This is particularly notable in areas corresponding to the most uncertain parts of the model as per Fig. 5. Inversion becomes increasingly guided by geology in these regions as large portions of the model undergo a significant decrease in geological uncertainty from case (a) to (c). The result is that geologically conditioned petrophysical constraints guide the inversion to update the model preferentially in the most geologically uncertain parts of the model. The consequence of this is an indirect focusing effect of geophysical inversion on uncertain areas.

The comparison of inversion results obtained for scenarios (b, α), (b, β) and (b, γ) (Fig. 7) and (a, β), (b, β), (c, β) (Fig. 8) shows that larger structures corresponding to ultramafic rocks (lithology 1), diorite (lithology 2) and basin fills (lithologies 3, 5 and 6) are coherent for all cases. It also shows that the differences occur at interfaces and in regions where geological uncertainty varies the most. This translates in pronounced differences between the recovered lithological models in Fig. 9(b), notably in areas marked by an arrow. It also corroborates previous observations, showing that the most notable changes are related to variations in geological uncertainty, while changes related
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Maximum observation probability model for PGMs corresponding to cases (a), (b) and (c) (i.e. decreasing measurement uncertainty) as per Table 2. High values (in white) indicate well-determined areas, while low values materialized by darker shades (e.g. lower values) indicate poorly constrained areas where no lithology dominates.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** simulated global mixture models, with decreasing standard deviations (in eq. 8) from case α to γ as per Table 2. The mean value for each lithology is marked by a cross, of which the colour corresponds to the respective lithology as per Fig. 3(a). Lithologies 5 and 3 completely overlapped because they present the same density contrast and magnetic susceptibility.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Inverted density contrast models (left-hand side rectangle) for cases (b, α), (b, β), (b, γ), and (b) corresponding absolute difference with reference model (right-hand side rectangle).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Inverted magnetic susceptibility models (left-hand side rectangle) and (B) difference with reference model (right-hand side rectangle) for cases (a, β), (b, β), (c, β) from top to bottom.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Recovered lithological model for cases (b, α), (b, β), (b, γ) (left-hand side ransangle) and (b) recovered lithological model for cases (a, β), (b, β), (c, β) (right-hand side rectangle). Note that lithology 5 does not appear as it is merged with lithology 3. The recovered models are overlaid with the boundary of the geological units from the reference model.</figcaption></figure>

to uncertainty in petrophysics are limited to small-scale or low magnitude changes. This observation is also reflected in Fig. 9(a) where the recovered lithological models show largely similar features. The next subsection focuses on quantitative quality indicators and on the topological analysis of recovered lithological models to complete and confirm our interpretations.

## 5.2 Quality indicators

### 5.2.1 Misfits and lithological resemblance

The MAMM $\phi_m$ is derived from the absolute model misfit (eq. 13) and reflects the average differences between inverted and reference models 4. As can be seen in Fig. 10(b), the main control on this indicator is uncertainty in geological input for both density contrast and magnetic susceptibility and it dominates over changes in petrophysical uncertainty. The effect of petrophysics is, therefore, to tune and sharpen the model, while geology strongly influences the structural features recovered by inversion. This confirms observations made in Section 5.1.

Fig. 10 also indicates that the lowest MAMM values for both density contrast and magnetic susceptibility correspond to the intermediate petrophysics case (β) in all three geological uncertainty scenarios. From these observations, we can argue that case (c, β) is the optimum inverted model of this study and that in general, intermediate petrophysical uncertainty similar to case (β) provides optimum results, while inverted model quality is directly related to geological input uncertainty. This interpretation is also valid for lithological resemblance *r* (eq. 12) in Fig. 10(c), for which the main control is geological input uncertainty. A possible explanation to this is that while potential field data inversion is highly affected by non-uniqueness, local entrapment is more likely to occur when the elements of the mixture model present the narrowest standard deviations, and reduces possible changes in the model during inversion. Conversely, normal distributions with broad standard deviations constrain inversion weakly. This interpretation is not contradicted by the inverted data RMSE values (eq. 16 next), which
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Metrics for inversion results analysis: MAMM indicator for density contrast (A) and magnetic susceptibility (B), lithological resemblance r (C), for the different uncertainty scenarios as per Fig. 4. Each colour corresponds to a separate geological uncertainty scenario (a, b, and c; e.g the lines in Fig. 4). Petrophysical uncertainty scenarios (α, b, and γ; e.g. columns in Fig. 4), are differentiated by marker type. The coloured lines connect values attached to the corresponding geological uncertainty scenario for visualization purposes but do indicate a continuum between the points.</figcaption></figure>

<table><caption>**Table 3.** Data RMSE for the tested scenarios.</caption><thead><tr><th rowspan="2"></th><th colspan="9">Uncertainty case scenario</th></tr><tr><th>(a,α)</th><th>(a,β)</th><th>(a,γ)</th><th>(b,α)</th><th>(b,β)</th><th>(b,γ)</th><th>(c,α)</th><th>(c,β)</th><th>(c,γ)</th></tr></thead><tbody><tr><td>Gravity data RMSE</td><td>1.026</td><td>1.028</td><td>1.003</td><td>1.043</td><td>1.045</td><td>1.034</td><td>1.033</td><td>1.039</td><td>1.013</td></tr><tr><td>Magnetic data RMSE</td><td>1.023</td><td>1.01</td><td>1.055</td><td>1.033</td><td>1.014</td><td>1.009</td><td>1.021</td><td>1.032</td><td>1.031</td></tr></tbody></table>

reach similar values for all inverted models (see values in Table 3). We calculate the normalized data RMSE as

$$ \text{RMSE} (\mathbf{m}^P, \mathbf{d}^P) = \left[ \frac{1}{n_d (\sigma_d^P)^2} \phi_d (\mathbf{m}^P, \mathbf{d}_P) \mathbf{d}_P^T \mathbf{d}_P \right]^{1/2} \quad (16) $$

where *P* denotes the problem considered, that is, gravity (G) or magnetic (M) as per eq. (3). The product $\mathbf{d}_P^T \mathbf{d}_P$ provides the sum-of-squares of the data. Values of $\sigma_d^P$ are given in Section 4.3.1, and $\phi_d$ is the corresponding data misfit obtained as explained in Section 2.3.2. Finally, examination of model misfit and lithological resemblance (Fig. 10) clearly illustrates the nonlinearity of the process leading to the lithological model. For instance, MAMM for scenario (c) is significantly lower than for scenario (a) (Figs 10a and b). Meanwhile, *r* increases from 87 per cent to 93 per cent between scenarios (a) and (c) (Fig. 10c). While geological uncertainty seems to be the main influencing factor, visible differences between the simulated scenarios occur primarily in the more geologically complex areas. This makes their geological interpretation challenging and highlights the necessity to study reconstructed lithological models in a geology-related manner.

### 5.2.2 Topology of lithological models

The topology of lithological models is calculated following the technique detailed in Section 3.2. Fig. 11 shows the relative difference between the adjacency matrix of the reference model and that of inverted models, encapsulated in $M^{\text{rel}}$ (eq. 15).

Fig. 11 shows that adjacency matrices consistently exhibit higher relative differences for contacts between lithology 4 and (3+5) (remembering that 3 and 5 are indistinguishable in petrophysical domain, see Fig. 6) than for other lithological relationships. A possible explanation is that lithology 4 (dolerite) has lesser thickness than other lithologies (see Fig. 3) in the more complex and uncertain regions of the model (Fig. 5). Therefore, it may not be recovered by inversion in areas where geophysical inversion cannot resolve it given the constraints applied. In such cases, because lithology (3+5) is intruded by dolerite (lithology 4), the geological meaning and topology of the model are strongly impacted as one may expect the intrusive lithology to be in contact with all the units it cuts. Negative $M^{\text{rel}}$ values for lithologies 4 and (3+5) indicate an underestimate of the number of contacts, suggesting that the intrusion is poorly recovered for cases (a) and (b) but that it is better constrained by cases (c).

The interpretation of Fig. 11 confirms results from other metrics. For geological uncertainty scenarios (a), (b) and (c), the best results are obtained with petrophysics case (β). Comparison of the different adjacency matrices indicates that results obtained with petrophysics case (β) shows better agreement with the reference model. It also demonstrates that geology is the main constraining factor for inversion, and that case (c, β) is optimum in terms of model recovery. Topological analysis as shown in Fig. 11 also corroborates the observations made in Section 5.2.1. in that geological input data uncertainty is the principal factor affecting the recovered models and that petrophysics is tuning the model. Furthermore, this confirms that knowledge of local changes in retrieved models induced by changes in petrophysical uncertainty are important for interpretation, and that the ‘sharpest’ petrophysical distributions do not necessarily lead to optimum results.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Relative difference between the adjacency matrices of the recovered lithological model and of the reference model for the different uncertainty scenarios.</figcaption></figure>

# 6 DISCUSSION

Contrary to intuitive assumptions, a reduction in the spread of physical properties around the mean value measured for each of the lithologies used to derive local constraints does not necessarily improve inverted models. As we have seen in Section 5.2.1 and summary, for each of the PGMs that we tested, petrophysics case $\beta$ leads to better results in terms of model misfit $\phi_m$ and lithological resemblance $r$ than case $\gamma$. In contrast, we observed that a decrease in geological uncertainty is directly related to improvements in the recovered model. This difference resides in that, at least in a controlled environment, better informed geology leads to a starting model that is more representative of the causative model. Importantly, reducing geological uncertainty decreases the effect of non-uniqueness further by reducing the range of models honouring the local constraints. We expect that the positive relationship between geological uncertainty and non-uniqueness will hold true in real-world scenarios based on test performance.

The laws governing geological modelling and geological processes are highly nonlinear. Consequently, the PGM is not simply a smooth version of the reference (or causative, in real-world studies) model as it encapsulates the complexity of all the geologically plausible models from MCUE. This is why topological analysis of the inverted models might reveal the presence of features different or absent from the reference model(s) or from the interpreted geology of the area. Lower geological uncertainty levels reduce the variability of possible models and the associated biases. Consequently, as we have seen, inverted models are strongly impacted by the level of certainty in geological field measurements.

Although the description of physical measurements using normal distributions cannot be applied to all scenarios, our findings regarding the influence of the related uncertainty may be extrapolated, to a certain extent, to the general case. Uncertainty in geological and petrophysical measurements affect inversion largely through their influence on the topography of the constraint term. Consequently, our results can be generalized to other types of distribution describing physical or petrophysical measurements provided that the related uncertainty impacts petrophysical constraints in a similar way. Our utilization of Gaussian mixture models relies on a soft clustering technique that shares a number of characteristics with the fuzzy c-means algorithms (Bordogna & Pasi 2011). We can, therefore, assume that our findings would hold true for inversions utilizing petrophysical constraints relying on soft clustering techniques such as the fuzzy c-means algorithms of Paasche & Tronicke (2007), Carter-McAuslan et al. (2015), Sun & Li (2015, 2016b) and Maag & Li (2018).
The main differences between inversions results shown here are characterized during inverted model analysis since ameliorations brought by constrained joint inversion occur primarily in model space. This confirms observations of Giraud *et al.* (2017), and to a lesser extent, De Stefano *et al.* (2011), Medina *et al.* (2012), Heincke *et al.* (2014, 2017), Linde & Doetsch (2016), Sun & Li (2016a) and Colombo *et al.* (2017), who highlight the interest of joint inversion for improved geological interpretation, and Gallardo & Meju (2004, 2011), Abubakar *et al.* (2012), Gao *et al.* (2012), Jardani *et al.* (2013), Molodtsov *et al.* (2013) and Juhojuntti & Kamm (2015), who stress the fact that improvements brought by joint inversion occur predominantly in model space. Complementing these works, we used a series of metrics to assess the reliability of inverted models for geological interpretation. Our analysis of lithological models complements the work Sun & Li (2016a), who conclude that ‘no additional analysis (of geology differentiation) after inversion is needed’ in a case where geological information is poorly informative or unavailable. This needs to be amended in cases where geological information or modelling is available. As we have demonstrated, geological analysis of recovered lithological models is useful to identify implausible parts of the model or to suggest alternate geological hypotheses. In our case, posterior model examination allows us to highlight poorly constrained areas or lithologies, thus confirming observations of Paasche (2016).

Future developments of the methodology introduced here include the matching of recovered lithological models with the closest geological realization from the suite of models that is used to derive the PGM. A possible extension of the presented methodology that increases the degree of integration between the three disciplines considered here is the incorporation of topological rules from area-specific geological knowledge directly in inversion. This would provide a theoretical framework unifying the inversion strategy we proposed with the works of Lelièvre *et al.* (2015), Miernik *et al.* (2016), Bijani *et al.* (2017) and Fullagar & Pears (2007), where the inverse model is formulated using surface geometry on one hand, and Carter-McAuslan *et al.* (2015), Sun & Li (2015, 2016a, 2017) who use petrophysical information on the other hand. Nevertheless, the applicability of such a modelling scheme is currently limited by the high computational requirements of level sets and wireframe inversion (Lelièvre *et al.* 2015; Lelièvre & Farquharson 2016; Bijani *et al.* 2017), or using hypotheses-limiting in application to specific scenarios (Juhojuntti & Kamm 2015). The same remark applies to the formulation of our joint inversion problem in a multi-objective global optimization scheme, which, while presenting the advantage of being less affected by non-uniqueness, would be limited by high computational costs. Nonetheless, we believe that the implementation of such methodology in 3-D is within reach.

Future works include the application of the method presented here to real-world geophysical data. Ongoing development of Tomofast-x includes the integration of spatial trends, geochemical information in the geological conditioning process and locally varying $C_m$ matrices. Further work also includes tests involving lithology- and location-dependent petrophysical uncertainty to better account for spatial variability of rock properties.

Our topological analysis is a first step towards quantitative posterior geological evaluation of inversion results. It is, however, restricted to the most common and simple geological relationships. Ideally, it should account for as many topological relationships as possible to provide comprehensive topological analysis. In addition, the study presented here could be extended by increasing the complexity of our geological model further through the addition of one or more intrusive bodies not sampled by surface geology or not sampled by petrophysics. In such a case, we expect to observe less sharp contrasts in the physical properties and a topological signature incompatible with the measured geological data in the vicinity of the intrusion.

# 7 CONCLUSIONS

We have studied how petrophysical and geological uncertainty are inputs to and influence geophysical inversion results using a realistic synthetic case study. We have presented the capabilities of our integrated inversion platform and shown how it integrates uncertain petrophysical and geological data in 3-D geophysical inversion.

We explored a relatively new and sparsely documented area of inversion and integrated studies by simulating a series of uncertainty levels in prior information. Our results clearly indicate that in geoscientific integration studies similar to ours in philosophy, the main control on the inverse model’s geological features is geological information. Meanwhile, petrophysics exerts a significantly lower influence that is, as a first-order approximation, restricted to enhancing model resolution. We have also shown that the topological analysis of recovered lithological models is a useful step in geological interpretation. Topological analysis is also crucial for the utilization of inverted models for further modelling because it can reveal poorly constrained geological units, thereby allowing the identification of alternative scenarios and zones to be investigated in more detail.

# 8 DATA, MODELS AND CODE AVAILABILITY

Reference property models, synthetic geophysical data, inversion model and recovered lithological models shown or discussed in this paper are made available by Giraud *et al.* (2018) in an ASCII format usable by Tomofast-x using doi: 10.5281/zenodo.1003105. The source code and the parameter files will be available on Github, project Tomofast-x (Ogarko *et al.* 2019). We used the publicly available structural geological model of the Mansfield area (Victoria, Australia) of Pakyuz-Charrier (2018) as the reference geological structural framework.
ACKNOWLEDGEMENTS

We thank the CALMIP supercomputing mesocentre of Toulouse, France, for their support through Roland Martin's projects #P1138_2017 and #P1138_2018 on the EOS machine (https://www.top500.org/system/178436). We thank P. Kovesi for sharing publicly 'the perceptually uniform colorcore maps' P. Kovesi (Good Colour Maps: How to Design Them, arXiv:1509.03700 [cs.GR] 2015). The authors are thankful to the Australian federal government for granting an International Postgraduate Research Scholarship to JG. They acknowledge the state government of Western Australia for supporting MJ and ML through the Geological Survey of Western Australia, Royalties for Regions and the Exploration Incentive Scheme. MJ is supported by a Western Australian Fellowship. ML is supported by the Mineral Research Institute of Western Australia. The authors finally thank the anonymous reviewers for constructive comments about the manuscript. VO acknowledges the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3-D (ASTRO 3-D) for supporting some of his research efforts.

REFERENCES

Abubakar, A., Gao, G., Habashy, T.M. & Liu, J., 2012. Joint inversion approaches for geophysical electromagnetic and elastic full-waveform data, *Inverse Probl.*, **28**, 055016, doi:10.1088/0266-5611/28/5/055016.

Airo, M., 2005, Regional interpretation of aerogeophysical data: extracting compositional and structural features, *Geol. Surv. Finl.*, **39**, 176-197.

Alcalde, J., Bond, C.E., Johnson, G., Ellis, J.F. & Butler, R.W.H., 2017. Impact of seismic image quality on fault interpretation uncertainty, *GSA Today*, **27**, 4-10.

Allmendinger, R.W., Siron, C.R. & Scott, C.P., 2017. Structural data collection with mobile devices: accuracy, redundancy, and best practices, *J. Struct. Geol.*, **102**, 98-112.

Barlow, M., 2004. Density and susceptibility characterisation of the minerie 100,000 geology sheet: implications for detailed gravity inversion, in *Y2 Project PMD\*CRC*, Commonwealth of Australia, Geoscience Australia.

Bhattacharyya, B.K., 1964. Magnetic anomalies due to prism-shaped bodies with arbitrary polarization, *Geophysics*, **29**, 517-531.

Bijani, R., Lelièvre, P.G., Ponte-Neto, C.F. & Farquharson, C.G., 2017. Physical-property-, lithology- and surface-geometry-based joint inversion using Pareto Multi-Objective Global Optimization, *Geophys. J. Int.*, **209**, 730-748.

Bond, C.E., 2015. Uncertainty in structural interpretation: lessons to be learnt, *J. Struct. Geol.*, **74**, 185-200.

Bond, C.E., Gibbs, A.D., Shipton, Z.K. & Jones, S., 2007. What do you think this is? "Conceptual uncertainty" in geoscience interpretation, *GSA Today*, **17**, 4-10.

Bond, C.E., Johnson, G. & Ellis, J.F., 2015. Structural model creation: the impact of data type and creative space on geological reasoning and interpretation, *Geol. Soc. Lond., Spec. Publ.*, **421**, 83-97.

Bordogna, G. & Pasi, G., 2011. Soft clustering for information retrieval applications, *WIREs Data Min. Knowl. Discovery*, **1**, 138-146.

Bosch, M., Bertorelli, G., Álvarez, G., Moreno, A. & Colmenares, R., 2015. Reservoir uncertainty description via petrophysical inversion of seismic data, *Leading Edge*, **34**, 1018-1026.

Bosch, M., Mukerji, T. & Gonzalez, E.F., 2010. Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: a review, *Geophysics*, **75**, 75A165-75A176.

Boulanger, O. & Chouteau, M., 2001. Constraints in 3D gravity inversion, *Geophys. Prospect.*, **49**, 265-280.

Brown, V., Key, K. & Singh, S., 2012. Seismically regularized controlled-source electromagnetic inversion, *Geophysics*, **77**, E57-E65.

Burns, K., 1988. Lithologic topology and structural vector fields applied to subsurface prediction in geology, in *Proceedings of the GIS/LIS 88*, San Antonio, TX, USA, pp. 26-34.

Calcagno, P., Chilès, J.P., Courrioux, G. & Guillen, A., 2008. Geological modelling from field data and geological knowledge. Part I. Modelling method coupling 3D potential-field interpolation and geological rules, *Phys. Earth planet. Inter.*, **171**, 147-157.

Carter-McAuslan, A., Lelièvre, P.G. & Farquharson, C.G., 2015. A study of fuzzy c-means coupling for joint inversion, using seismic tomography and gravity data test scenarios, *Geophysics*, **80**, W1-W15.

Cartwright, N. & McMullin, E., 1984. *How the laws of physics lie, Am. J. Phys.*, **52**, 474-476.

Chapman, C.H. & Pratt, R.G., 1992. Traveltime tomography in anisotropic media-I. Theory, *Geophys. J. Int.*, **109**, 1-19.

Clark, D.A. & Emerson, D.W., 2006. Notes on rock magnetization characteristics in applied geophysical studies, *Explor. Geophys.*, **22**, 547-555.

Colombo, D., Rovetta, D., Turkoglu, E., Mcneice, G. & Sandoval-curiel, E., 2017. Multiscale hierarchical seismic-CSEM joint inversion for subsalt depth imaging in the Red Sea, in *SEG Technical Program Expanded Abstracts 2017*, pp. 1169-1173, doi:10.1190/segam2017-17558705.1.

Davis, J., 2002. *Statistics and Data Analysis in Geology*, 3rd edn, Wiley.

de la Varga, M. & Wellmann, J.F., 2016. Structural geologic modeling as an inference problem: a Bayesian perspective, *Interpretation*, **4**, SM1-SM16.

Demirel, C. & Candansayar, M.E., 2017. Two-dimensional joint inversions of cross-hole resistivity data and resolution analysis of combined arrays, *Geophys. Prospect.*, **65**, 876-890.

De Stefano, M., Golfré Andreasi, F., Re, S., Virgilio, M. & Snyder, F.F., 2011. Multiple-domain, simultaneous joint inversion of geophysical data with application to subsalt imaging, *Geophysics*, **76**, R69-R80.

Doetsch, J., Linde, N., Coscia, I., Greenhalgh, S.A. & Green, A.G., 2010. Zonation for 3D aquifer characterization based on joint inversions of multimethod crosshole geophysical data, *Geophysics*, **75**, G53-G64.

Egenhofer, M.J. & Herring, J., 1990. Categorizing binary topological relations between regions, lines, and points in geographic databases, 1-28, Tech. rep. 90-12. National Center for Geographic Information and Analysis, University of California, Santa Barbara, Santa Barbara, CA.

Fernández-Martínez, J.L., Pallero, J.L.G., Fernández-Muñiz, Z. & Pedruelo-González, L.M., 2014a. The effect of noise and Tikhonov's regularization in inverse problems. Part II: The nonlinear case, *J. Appl. Geophys.*, **108**, 186-193.

Fernández-Martínez, J.L., Pallero, J.L.G., Fernández-Muñiz, Z. & Pedruelo-González, L.M., 2014b. The effect of noise and Tikhonov's regularization in inverse problems. Part I: The linear case, *J. Appl. Geophys.*, **108**, 176-185.

Fullagar, P. & Pears, G., 2007. Towards geologically realistic inversion, in *Proceedings of Exploration 07: Fifth Decennial International Conference on Mineral Exploration*, Toronto, Canada, pp. 444-460, ed. Milkereit, B.

Gallardo, L.A., Fontes, S.L., Meju, M.A., Buonora, M.P. & de Lugao, P.P., 2012. Robust geophysical integration through structure-coupled joint inversion and multispectral fusion of seismic reflection, magnetotelluric, magnetic, and gravity images: example from Santos Basin, offshore Brazil, *Geophysics*, **77**, B237-B251.

Gallardo, L.A. & Meju, M.A., 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of dc resistivity and seismic data, *Geophys. Res. Lett.*, **30**.

Gallardo, L.A. & Meju, M.A., 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints, *J. geophys. Res.*, **109**.

Gallardo, L.A. & Meju, M.A., 2007. Joint two-dimensional cross-gradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification, *Geophys. J. Int.*, **169**, 1261-1272.

Gallardo, L.A. & Meju, M.A., 2011. Structure-coupled multiphysics imaging in geophysical sciences, *Rev. Geophys.*, **49**, RG1003.

Gao, G., Abubakar, A. & Habashy, T.M., 2012. Joint petrophysical inversion of electromagnetic and full-waveform seismic data, *Geophysics*, **77**, WA3
Giraud, J., Jessell, M., Lindsay, M., Martin, R., Parkyuz-Charrier, E. & Ogarko, V., 2016b. Geophysical joint inversion using statistical petrophysical constraints and prior information, in *ASEG Extended Abstracts 2016, 25th International Conference and Exhibition – Interpreting the Past, Discovering the Future*, Adelaide, pp. 1–6.

Giraud, J., Jessell, M., Lindsay, M., Parkyuz-Charrier, E. & Martin, R., 2016a. Integrated geophysical joint inversion using petrophysical constraints and geological modelling, in *SEG Technical Program Expanded Abstracts 2016*, Tulsa, pp. 1597–1601, doi/10.1190/segam2016-13945549.1.

Giraud, J., Lindsay, M. & Ogarko, V., 2018. Yerrida Basin Geophysical Modeling—input data and inverted models, doi:10.5281/zenodo.1238216.

Giraud, J., Lindsay, M., Ogarko, V., Jessell, M., Martin, R. & Pakyuz-Charrier, E., 2019. Integration of geological uncertainty into geophysical inversion by means of local gradient regularization, *Solid Earth*, **10**, 193–210.

Giraud, J., Pakyuz-Charrier, E., Jessell, M., Lindsay, M., Martin, R. & Ogarko, V., 2017. Uncertainty reduction through geologically conditioned petrophysical constraints in joint inversion, *Geophysics*, **82**, ID19–ID34.

Godsil, C. & Royle, G., 2001. *Algebraic Graph Theory*, Springer.

Guo, Z., Dong, H. & Kristensen, Å., 2017. Image-guided regularization of marine electromagnetic inversion, *Geophysics*, **82**, E221–E232.

Haber, E. & Oldenburg, D., 1997. Joint inversion: a structural approach, *Inverse Probl.*, **13**, 63–77.

Hammond, G.E., Lichtner, P.C. & Mills, R.T., 2014. Evaluating the performance of parallel subsurface simulators: an illustrative example with PFLOTRAN, *Water Resour. Res.*, **50**, 208–228.

Hansen, P.C. & Johnston, P.R., 2001. The L-Curve and its use in the numerical treatment of inverse problems, in *Computational Inverse Problems in Electrocardiography*, pp. 119–142, WIT Press, Johnston, P. WIT.

Hansen, P.C. & O'Leary, D.P., 1993. The use of the L-curve in the regularization of discrete ill-posed problems, *SIAM J. Sci. Comput.*, **14**, 1487–1503.

Hatfield, K.L., Evans, A.J. & Harvey, P.K., 2002. Defining Petrophysical Units of the Palmer Deep Sites from Leg 178. *Proceedings of the Ocean Drilling Program*, 178 Scientific Results, Ocean Drilling Program.

Heincke, B., Jegen, M., Moorkamp, M. & Hobbs, R.W., 2014. Joint-inversion of magnetotelluric, gravity and seismic data to image sub-basalt sediments offshore the Faroe-Islands, in *SEG Technical Program Expanded Abstracts 2014*, Tulsa, pp. 770–775.

Heincke, B., Jegen, M., Moorkamp, M., Hobbs, R.W. & Chen, J., 2017. An adaptive coupling strategy for joint inversions that use petrophysical information as constraints, *J. Appl. Geophys.*, **136**, 279–297.

Hempel, C.G. & Oppenheim, P., 1948. Studies in the logic of explanation, *Phil. Sci.*, **15**, 135–175.

Hunt, C.P., Moskowitz, B.M. & Banerjee, S.K., 1995. Magnetic properties of rocks and minerals, in *Handbook of Physical Constants*, Vol. 3, pp. 189–204, ed. Ahrens, T.J., American Geophysical Union.

Jardani, A., Revil, A. & Dupont, J.P., 2013. Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging, *Adv. Water Res.*, **52**, 62–77.

Jessell, M.W., Ailleres, L. & de Kemp, E.A., 2010. Towards an integrated inversion of geoscientific data: what price of geology? *Tectonophysics*, **490**, 294–306.

Jessell, M.W., Pakyuz-charrier, E., Lindsay, M., Giraud, J. & de Kemp, E., 2018. Assessing and mitigating uncertainty in three-dimensional geologic models in contrasting geologic scenarios, in *Metals, Minerals, and Society*, Vol. 21, pp. 63–74, eds Arribas, A.M. & Mauk, J.L., SEG.

Juhojuntti, N. & Kamm, J., 2015. Joint inversion of seismic refraction and resistivity data using layered models—applications to groundwater investigation, *Geophysics*, **80**, EN43–EN55.

Kamm, J., Lundin, I.A., Bastani, M., Sadeghi, M. & Pedersen, L.B., 2015. Joint inversion of gravity, magnetic, and petrophysical data—a case study from a gabbro intrusion in Boden, Sweden, *Geophysics*, **80**, B131–B152.

Kinkeldey, C., MacEachren, A.M., Riveiro, M. & Schiewe, J., 2015. Evaluating the effect of visually represented geodata uncertainty on

decision-making: systematic review, lessons learned, and recommendations, *Cartography Geogr. Inf. Sci.*, **44**, 1–21.

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

Lajaukie, C., Courrioux, G. & Manuel, L., 1997. Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation, *Math. Geol.*, **29**, 571–584.

Laplace, P-S., 1774. Mémoire sur la probabilité des causes par les événements, *Mem. Acad. R. Sci. Paris*, **6**, 621–656.

Lark, R.M., Mathers, S.J., Thorpe, S., Arkley, S.L.B., Morgan, D.J. & Lawrence, D.J.D., 2013. A statistical assessment of the uncertainty in a 3-D geological framework model, *Proc. Geologists' Assoc.*, **124**, 946–958.

Lelièvre, P., Carter-McAuslan, A., Farquharson, C. & Hurich, C., 2012b. Unified geophysical and geological 3D Earth models, *Leading Edge*, **31**, 322–328.

Lelièvre, P., Farquharson, C. & Bijani, R., 2015. 3D stochastic geophysical inversion for contact surface geometry, in *EGU General Assembly 2015*, Vienna, Austria, Vol. **17**, pp. 3627.

Lelièvre, P., Farquharson, C. & Hurich, C., 2012a. Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration, *Geophysics*, **77**, K1–K15.

Lelièvre, P.G. & Farquharson, C.G., 2013. Gradient and smoothness regularization operators for geophysical inversion on unstructured meshes, *Geophys. J. Int.*, **195**, 330–341.

Lelièvre, P.G. & Farquharson, C.G., 2016. Integrated imaging for mineral exploration, in *Integrated Imaging of the Earth: Theory and Applications*, eds Moorkamp, M., Lelièvre, P.G., Linde, N. & Khan, A., pp. 137–166, Wiley.

Linde, N. & Doetsch, J., 2016. Joint inversion in hydrogeophysics and near-surface geophysics, in *Integrated Imaging of the Earth: Theory and Applications*, pp. 117–135, eds Moorkamp, M., Lelièvre, P.G., Linde, N. & Khan, A., Wiley.

Lindsay, M., Perrouty, S., Jessell, M. & Ailleres, L., 2014. Inversion and geodiversity: searching model space for the answers, *Math. Geosci.*, **46**, 971–1010.

Lindsay, M.D., Aillères, L., Jessell, M.W., de Kemp, E.A. & Betts, P.G., 2012. Locating and quantifying geological uncertainty in three-dimensional models: analysis of the Gippsland Basin, southeastern Australia, *Tectonophysics*, **546–547**, 10–27.

Lindsay, M.D., Jessell, M.W., Ailleres, L., Perrouty, S., de Kemp, E. & Betts, P.G., 2013b. Geodiversity: exploration of 3D geological model space, *Tectonophysics*, **594**, 27–37.

Lindsay, M.D., Perrouty, S., Jessell, M.W. & Aillères, L., 2013a. Making the link between geological and geophysical uncertainty: geodiversity in the Ashanti Greenstone Belt, *Geophys. J. Int.*, **195**, 903–922.

Lipari, V., Urbano, D., Spadavecchia, E., Panizzardi, J. & Bienati, N., 2017. Regularized tomographic inversion with geological constraints, *Geophys. Prospect.*, **65**, 305–315.

Li, Y. & Oldenburg, D.W., 2000. Joint inversion of surface and three-component borehole magnetic data, *Geophysics*, **65**, 540–552.

Li, Y. & Sun, J., 2016. Geology differentiation with uncertainty estimation using inverted magnetization directions, in *SEG Technical Program Expanded Abstracts 2016*, pp. 2159–2164.

Maag, E. & Li, Y., 2018. Discrete-valued gravity inversion using the guided fuzzy c-means clustering technique, *Geophysics*, **83**, G59–G77.

Mardia, K.V. & Jupp, P.E., 2008. *Directional Statistics*, pp. 1–8, Wiley.

Martinez, C. & Li, Y., 2015. Lithologic characterization using airborne gravity gradient and aeromagnetic data for mineral exploration: a case study in the Quadrilátero Ferrífero, Brazil, *Interpretation*, **3**, SL1–SL13.

Martin, R., Monteiller, V., Komatitsch, D., Perrouty, S., Jessell, M., Bonvalot, S. & Lindsay, M.D., 2013. Gravity inversion using wavelet-based compression on parallel hybrid CPU/GPU systems: application to southwest Ghana, *Geophys. J. Int.*, **195**, 1594–1619.

Martin, R., Ogarko, V., Komatitsch, D. & Jessell, M., 2018. Parallel three-dimensional electrical capacitance data imaging using a nonlinear inversion algorithm and Lp norm-based model, *Measurement*, **128**, 428–445.
McLachlan, G. & Peel, D., 2000. *Finite Mixture Models*, Wiley Series in Probability and Statistics, Wiley, pp. 419.

Medina, E., Lovatini, A., Andreati, F.G., Re, S. & Snyder, F., 2012. Simultaneous joint inversion of 3D seismic and magnetotelluric data from the Walker Ridge, *First Break*, **30**, 85–88.

Melo, A.T., Sun, J. & Li, Y., 2017. Geophysical inversions applied to 3D geology characterization of an iron oxide copper-gold deposit in Brazil, *Geophysics*, **82**, K1–K13.

Miernik, K., Bogacz, A., Kozubal, A., Danek, T. & Wojdyla, M., 2016. Pareto joint inversion of 2d magnetotelluric and gravity data – towards practical applications, *Acta Geophys.*, **64**, 1655–1672.

Molodtsov, D.M., Troyan, V.N., Roslov, Y.V. & Zerilli, A., 2013. Joint inversion of seismic traveltimes and magnetotelluric data with a directed structural constraint, *Geophys. Prospect.*, **61**, 1218–1228.

Moorkamp, M., 2017. Integrating electromagnetic data with other geophysical observations for enhanced imaging of the Earth: a tutorial and review, *Surv. Geophys.*, 1–28.

Moorkamp, M., Heincke, B., Jegen, M., Hobbs, R.W. & Roberts, A.W., 2016. Joint inversion in hydrocarbon exploration, in *Integrated Imaging of the Earth: Theory and Applications*, pp. 167–189, eds Moorkamp, M., Lelièvre, P.G., Linde, N. & Khan, A., Wiley.

Moorkamp, M., Roberts, A.W., Jegen, M., Heincke, B. & Hobbs, R.W., 2013. Verification of velocity-resistivity relationships derived from structural joint inversion with borehole data, *Geophys. Res. Lett.*, **40**, 3596–3601.

Nearing, G.S., Tian, Y., Gupta, H.V., Clark, M.P., Harrison, K.W. & Weijs, S.V., 2016. A philosophical basis for hydrological uncertainty, *Hydrol. Sci. J.*, **61**, 1666–1678.

Neumann, J., Watts, M.D. & Hallinan, S., 2016. Geologically constrained 2D and 3D airborne EM inversion through cross-gradient regularization and multi-grid efficiency, in *ASEG-PESA-AIG 2016, 25th International Conference and Exhibition – Interpreting the Past, Discovering the Future*, Adelaid, pp. 1–6.

Novakova, L. & Pavlis, T.L., 2017. Assessment of the precision of smart phones and tablets for measurement of planar orientations: a case study, *J. Struct. Geol.*, **97**, 93–103.

Ogarko, V., Martin, R., Giraud, J. & Jessell, M.W., 2019. Introduction to the Tomofast-x Large-Scale, Three-Dimensional Parallel Integrated and Joint Inversion Platform.

Paasche, H., 2016. Post-Inversion integration of disparate tomographic models by model structure analyses, in *Integrated Imaging of the Earth, Theory and Applications*, pp. 69–91, eds Moorkamp, M., Lelièvre, P.G., Linde, N. & Khan, A., Wiley.

Paasche, H. & Tronicke, J., 2007. Cooperative inversion of 2D geophysical data sets: a zonal approach based on fuzzy c-means cluster analysis, *Geophysics*, **72**, A35–A39.

Paige, C.C. & Saunders, M.A., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares, *ACM Trans. Math. Softw.*, **8**, 43–71.

Pakyuz-Charrier, E., 2018. Mansfield (Victoria, Australia) area original GeoModeller model and relevant MCUE outputs, doi:10.5281/zenodo.848225.

Pakyuz-Charrier, E., Giraud, J., Lindsay, M. & Jessell, M., 2018a. Common uncertainty research explorer uncertainty estimation in geological 3d modelling, in *ASEG Extended Abstracts*, Perth, Vol. 2018, pp. 1–6.

Pakyuz-Charrier, E., Giraud, J., Ogarko, V., Lindsay, M. & Jessell, M., 2018c. Drillhole uncertainty propagation for three-dimensional geological modeling using Monte Carlo, *Tectonophysics*, **747–748**, 16–39.

Pakyuz-Charrier, E., Lindsay, M., Ogarko, V., Giraud, J. & Jessell, M., 2018b. Monte Carlo simulation for uncertainty estimation on structural data in implicit 3-D geological modeling, a guide for disturbance distribution selection and parameterization, *Solid Earth*, **9**, 385–402.

Park, H., Scheidt, C., Fenwick, D., Boucher, A. & Caers, J., 2013. History matching and uncertainty quantification of facies models with multiple geological interpretations, *Comput. Geosci.*, **17**, 609–621.

Parsekian, A.D., Singha, K., Minsley, B.J., Holbrook, W.S. & Slater, L., 2015. Multiscale geophysical imaging of the critical zone, *Rev. Geophys.*, **53**, 1–26.

Pellerin, J., Botella, A., Bonneau, F., Mazuyer, A., Chauvin, B., Lévy, B. & Caumon, G., 2017. RINGMesh: a programming library for developing mesh-based geomodeling applications, *Comput. Geosci.*, **104**, 93–100.

Perrin, M. & Rainaud, J.-F., 2013. *Shared Earth Modeling: Knowledge Driven Solutions for Building and Managing Subsurface 3D Geological Models*, Technip.

Pilkington, M. & Shamsipour, P., 2014. Noise reduction procedures for gravity-gradiometer data, *Geophysics*, **79**, doi:10.1190/geo2014-0084.1.

Portniaguine, O. & Zhdanov, M.S., 2002. 3-D magnetic inversion with data compression and image focusing, *Geophysics*, **67**, 1532–1541.

Pratt, R.G. & Chapman, C.H., 1992. Traveltime tomography in anisotropic media-II. Application, *Geophys. J. Int.*, **109**, 20–37.

Pratt, R.G., Shin, C. & Hicks, G.J., 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion, *Geophys. J. Int.*, **133**, 341–362.

Revil, A., Cuttler, S., Karaoulis, M., Zhou, J., Raynolds, B. & Batzle, M., 2015. The plumbing system of the Pagosa thermal Springs, Colorado : application of geologically constrained geophysical inversion and data fusion, *J. Volc. Geotherm. Res.*, **299**, 1–18.

Rittgers, J.B., Revil, A., Mooney, M.A., Karaoulis, M., Wodajo, L. & Hickey, C.J., 2016. Time-lapse joint inversion of geophysical data with automatic joint constraints and dynamic attributes, *Geophys. J. Int.*, **207**, 1401–1419.

Sanger, E.A. & Glen, J.M.G., 2003. Density and magnetic susceptibility values for rocks in the Talkeetna Mountains and adjacent region, south-central Alaska, Open-File Rep, doi:10.3133/ofr03268.

Santos, E.T.F. & Bassrei, A., 2007. L- and Θ-curve approaches for the selection of regularization parameter in geophysical diffraction tomography, *Comput. Geosci.*, **33**, 618–629.

Schneeberger, R., La Varga, M.D., Egli, D., Berger, A., Kober, F., Wellmann, F. & Herwegh, M., 2017. Methods and uncertainty estimations of 3-D structural modelling in crystalline rocks: a case study, *Solid Earth*, **8**, 987–1002.

Schweizer, D., Blum, P. & Butscher, C., 2017. Uncertainty assessment in 3-D geological models of increasing complexity, *Solid Earth*, **8**, 515–530.

Stigler, S.M., 1986. Laplace's 1774 memoir on inverse probability, *Stat. Sci.*, **1**, 359–363.

Sun, J. & Li, Y., 2012, Joint inversion of multiple geophysical data : a petrophysical approach using guided fuzzy c-means clustering, in *SEG Expanded Abstracts 2012*, Las Vegas, pp. 1–5.

Sun, J. & Li, Y., 2013. Petrophysically constrained geophysical inversion using Parzen window density estimation, in *SEG Technical Program Expanded Abstracts 2013*, pp. 3051–3056, doi:10.1190/segam2013-1163.1.

Sun, J. & Li, Y., 2015. Multidomain petrophysically constrained inversion and geology differentiation using guided fuzzy c-means clustering, *Geophysics*, **80**, ID1–ID18.

Sun, J. & Li, Y., 2016a. Joint inversion of multiple geophysical data using guided fuzzy c-means clustering, *Geophysics*, **81**, ID37–ID57.

Sun, J. & Li, Y., 2016b. Joint clustering inversion of gravity and magnetic data applied to the imaging of a gabbro intrusion, in *SEG Expanded Abstracts 2016*, Dallas, pp. 2175–2179.

Sun, J. & Li, Y., 2017. Joint inversion of multiple geophysical and petrophysical data using generalized fuzzy clustering algorithms, *Geophys. J. Int.*, **208**, 1201–1216.

Tarantola, A. & Valette, B., 1982. Inverse problems = quest for information, *J. Geophys.*, **50**, 159–170.

Thiele, S.T., Jessell, M.W., Lindsay, M., Ogarko, V., Wellmann, J.F. & Pakyuz-Charrier, E., 2016b. The topology of geology 1: topological analysis, *J. Struct. Geol.*, **91**, 27–38.

Thiele, S.T., Jessell, M.W., Lindsay, M., Wellmann, J.F. & Pakyuz-Charrier, E., 2016a. The topology of geology 2: topological uncertainty, *J. Struct. Geol.*, **91**, 74–87.

Törnberg, R. & Sturkell, E.F.F., 2005. Density and magnetic susceptibility of rocks from the Lockne and Tvären marine impact structures, *Meteorit. Planet. Sci.*, **40**, 639–651.

Wellmann, J.F. & Regenauer-Lieb, K., 2012. Uncertainties have a meaning: information entropy as a quality measure for 3-D geological models, *Tectonophysics*, **526–529**, 207–216.
Wiik, T., Nordskag, J.I., Dischler, E.Ø. & Nguyen, A.K., 2015. Inversion of inline and broadside marine controlled-source electromagnetic data with constraints derived from seismic data, *Geophys. Prospect.*, **63**, 1371–1382.

Wilson, E.B., 1923. First and second laws of error, *J. Am. Stat. Assoc.*, **18**, 841–851.

Yuan, S., Wang, S. & Li, G., 2012. Random noise reduction using Bayesian inversion, *J. Geophys. Eng.*, **9**, doi:10.1088/1742-2132/9/1/007.

Zhang, J. & Revil, A., 2015. 2D joint inversion of geophysical data using petrophysical clustering and facies deformation, *Geophysics*, **80**, M69–M88.

Zhou, J., Revil, A. & Jardani, A., 2016. Stochastic structure-constrained image-guided inversion of geophysical data, *Geophysics*, **81**, E89–E101.

Zhou, J., Revil, A., Karaoulis, M., Hale, D., Doetsch, J. & Cuttler, S., 2014. Image-guided inversion of electrical resistivity data, *Geophys. J. Int.*, **197**, 292–309.

Zlatanova, S., 2000. On 3D topological relationships, in *Proceedings 11th Int. Workshop on Database and Expert Systems Applications*, London, UK, pp. 913–919.

# APPENDIX

## A1 L-surfaces for hyperparameter estimation

<figure><img src="image_14.png" /><figcaption>**Figure A1.** Points used for the l-surface analysis for gravity (left) and magnetic inversion (right).</figcaption></figure>
## A2 Complete set of inversion results

### A2.1 Density contrast

<figure><img src="image_15.png" /><figcaption>**Figure A2.** Set of inverted density contrast models arranged as per Fig. 4. The corresponding uncertainty case (Table 2) is given on the figure.</figcaption></figure>

### A2.2 Magnetic susceptibility

<figure><img src="image_16.png" /><figcaption>**Figure A3.** Set of inverted magnetic susceptibility models arranged as per Fig. 4. The corresponding uncertainty case (Table 2) is given on the figure.</figcaption></figure>
## A2.3 Differences between true and inverted density contrast

<figure><img src="image_17.png" /><figcaption>**Figure A4.** Different between true model and inverted density contrast models arranged as per Fig. 4. The corresponding uncertainty case (Table 2) is given on the figure.</figcaption></figure>

## A2.4 Differences between true and inverted magnetic susceptibility

<figure><img src="image_18.png" /><figcaption>**Figure A5.** Different between true model and inverted magnetic susceptibility models arranged as per Fig. 4. The corresponding uncertainty case (Table 2) is given on the figure.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure A6.** Recovered lithological models arranged as per Fig. 4. The corresponding uncertainty case (Table 2) is given on the figure.</figcaption></figure>