Petrophysically and geologically guided multi-physics inversion using a dynamic Gaussian mixture model

by   Thibaut Astic, et al.

In a previous paper, we introduced a framework for carrying out petrophysically and geologically guided geophysical inversions. In that framework, petrophysical and geological information is modelled with a Gaussian Mixture Model (GMM). In the inversion, the GMM serves as a prior for the geophysical model. The formulation was confined to problems in which a single physical property model was sought, with a single geophysical dataset. In this paper, we extend that framework to jointly invert multiple geophysical datasets that depend on multiple physical properties. The petrophysical and geological information is used to couple geophysical surveys that, otherwise, rely on independent physics. This requires advancements in two areas. First, an extension from a univariate to a multivariate analysis of the petrophysical data, and their inclusion within the inverse problem, is necessary. Second, we address the practical issues of simultaneously inverting data from multiple surveys and finding a solution that acceptably reproduces each one, along with the petrophysical and geological information. To illustrate the efficacy of our approach and the advantages of carrying out multi-physics inversions, we invert synthetic gravity and magnetic data associated with a kimberlite deposit. The kimberlite pipe contains two distinct facies embedded in a host rock. Inverting the datasets individually leads to a binary geological model: background or kimberlite. A multi-physics inversion, with petrophysical information, differentiates between the two main kimberlite facies of the pipe. Through this example, we also highlight the capabilities of our framework to work with interpretive geologic assumptions when minimal quantitative information is available. In those cases, the dynamic updates of the Gaussian Mixture Model allow us to perform multi-physics inversions by learning a petrophysical model.




Joint geophysical, petrophysical and geologic inversion using a dynamic Gaussian mixture model

We present a framework for petrophysically and geologically guided inver...

A probabilistic approach to emission-line galaxy classification

We invoke a Gaussian mixture model (GMM) to jointly analyse two traditio...

Dynamic Weights in Gaussian Mixture Models: A Bayesian Approach

In this paper we consider a Gaussian mixture model where the mixture wei...

Skewed Distributions or Transformations? Modelling Skewness for a Cluster Analysis

Because of its mathematical tractability, the Gaussian mixture model hol...

Directivity Modes of Earthquake Populations with Unsupervised Learning

We present a novel approach for resolving modes of rupture directivity i...

Enhance Feature Discrimination for Unsupervised Hashing

We introduce a novel approach to improve unsupervised hashing. Specifica...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Mineral deposits, or other geologic features, are characterized by different physical properties, and hence multiple geophysical surveys can be used to delineate them. For example, kimberlites have signatures that depend upon density, magnetic susceptibility, and electrical conductivity. They are often discovered through data collected in airborne surveys and appear as a circular low gravity anomaly, high magnetic response, and sometimes negative electromagnetic transient response (Macnae, 1995; Keating and Sailhac, 2004; Bournas et al., 2018). Although a joint interpretation of several individually inverted datasets can significantly improve our understanding of the subsurface (Oldenburg et al., 1997; Devriese et al., 2017; Kang et al., 2017; Giuseppe et al., 2014; Paasche et al., 2006; Paasche, 2016; Martinez and Li, 2015; Melo et al., 2017), multiple cases studies have shown that multi-physics inversions can reveal information that was not accessible through individual geophysical dataset inversions (Doetsch et al., 2010; Jegen et al., 2009; Kamm et al., 2015; Lelièvre and Farquharson, 2016). An extensive compilation of methods and their applications can be found in Moorkamp et al. (2016). Multi-physics inversions require a coupling term that mathematically describes a relationship between the different physical property models responsible for the geophysical data. Coupling methods generally use one or a combination of structural or physical property relationships.

The first frameworks for joint inversion focused on linking geophysical models through their structural similarities. Haber and Oldenburg (1997) defined the structure of a model in terms of the absolute value of its spatial curvature and compared different models to see if variations occurred at the same locations. This idea was further developed by Gallardo and Meju (2003) with the introduction of the concept of cross-gradient between geophysical models. This approach has become commonly used, and both Gallardo and Meju (2011) and Meju and Gallardo (2016) provide in-depth reviews of the method and its application. However, this strategy assumes that contrasts in physical properties are co-located, and that is not valid for many problems. This drawback can be overcome by using other coupling methods.

The second coupling approach uses physical property relationships to link geophysical models. Some of the earliest works used empirical constitutive formulae as their physical properties constraint (Afnimar et al., 2002; Hoversten et al., 2006; Chen et al., 2007). De Stefano et al. (2011) combined this approach with the above mentioned cross-gradient method for sub-salt imaging. Moorkamp et al. (2011) compared the constitutive relationship and the cross-gradient approaches on a 3-D synthetic example combining magnetotelluric, gravity, and seismic data. They concluded that, overall, a cross-gradient approach was preferable compared to using constitutive equations because deviations from the constitutive relations resulted in artefacts in the inverted models; in those situations, the cross-gradient method gave consistent satisfactory results. They also pointed out that the cross-gradient method relies on fewer assumptions about the models than the constitutive equations. Some stochastic frameworks have also been proposed that leverage geostatistical tools to define relationships between physical properties. Chen and Hoversten (2012) used a similar coupling approach as in Bosch (2004) by building a rock-physics model from borehole data to jointly invert seismic and controlled-source electromagnetic data. Shamsipour et al. (2012) used the geostatistical techniques of cokriging and conditional simulation to jointly invert gravity and magnetic data assuming that the auto- and cross-covariances of the density and magnetic susceptibility follow a linear model of coregionalization. On the deterministic side, of which the framework we present belongs, recent frameworks use clustering techniques such as the fuzzy C-means (FCM) algorithm, which was first used in Paasche and Tronicke (2007) and further expanded in Lelièvre et al. (2012). This approach adds a clustering term to the objective function, which allows more flexible relationships between physical properties. Beyond the addition of the FCM term to the objective function, Sun and Li (2015) introduce an iterative update to the cluster centre throughout an inversion for a single physical property, a technique they call guided FCM; this introduced the concept of uncertainties for physical properties into the inversion. In Sun and Li (2016), they generalize this work to consider multiple physical properties, and further in Sun and Li (2017), they added tools to their approach to consider various types of correlations between physical properties (linear, quadratic, etc.). Giraud et al. (2017) represented the petrophysical information as a fixed GMM and focused on reducing uncertainties in stochastic geological modelling by linking potential field inversions and geological models through petrophysical information. To this end, they added, to the Tikhonov objective function, a sum of least-squares differences between the GMM function, evaluated at the current model, and reference values representing the likelihood of their prior knowledge. In Giraud et al. (2019), they modified their formulation of this coupling term to work with a least-squares difference between the log-likelihood of the GMM and their reference values. In both formulations, they required extensive and fixed quantitative petrophysical and geological information. In practice, that information is not often available and might be only qualitative.

In Astic and Oldenburg (2019)

, we presented a petrophysically and geologically guided inversion (PGI) framework that generalized concepts presented in previously published researches. We showed how geological and petrophysical knowledge represented as a univariate GMM could be incorporated in a voxel-based geophysical inversion through a single smallness term. Each contrasting geological unit was represented by a univariate Gaussian distribution, which summarized its physical property signature. Geological information was included in the GMM through its proportions in a manner similar to that of

Giraud et al. (2017). The log-likelihood of the GMM was then used to regularize the geophysical inversion; this is analogous to the approach taken by Grana and Della Rossa (2010); Grana et al. (2017). Incorporating both petrophysical and geological information into a single smallness term in the regularization has several advantages. First, it does not require adding a term in the classic Tikhonov formulation. Second, this approach brings the petrophysical data to the same level as the geophysical data, which allows us to define a misfit, with a target value, between the geophysical model and the petrophysical and geological data. The iteration steps were decomposed into a suite of cyclic optimization problems across the geophysical, petrophysical, and geological data. The petrophysical step formalized the idea of learning the physical property mean values described in Sun and Li (2015)

and generalized it to enable the variances and proportions of the GMM to be learned as well during the inversion. These updates to the GMM parameters allowed us to work with partial petrophysical information. The geological step built at each iteration a ”quasi-geologic model”

(Li et al., 2019), based on the current geophysical and petrophysical models. We applied the PGI approach on synthetic and field data, but the analysis was restricted to single datasets and a single physical property.

This study extends the PGI framework to perform multi-physics inversions, involving several physical properties. We show how geological and petrophysical information represented as a multivariate GMM can be used to couple multiple voxel-based geophysical inversions through a single smallness term in the regularization. The updates to the means, covariances and proportions of the GMM are extended to inversions with multiple physical properties, which expands the work of Sun and Li (2016). Tools for handling various types of relationships between physical properties are designed; this further develops ideas presented in Sun and Li (2017). These capabilities, and the gains that they generate, are demonstrated on inversions of synthetic gravity and magnetic data.

In this paper, we start by reviewing the key concepts of the PGI approach and generalize them to the case of multiple physical properties and governing equations for performing multi-physics inversions. We then delineate our strategy to address the numerical challenges of finding a solution to the inverse problem that fits each dataset and, at the same time, adequately fits the petrophysical data. Finally, we demonstrate the advantages of performing joint inversions, with various levels of prior knowledge, by using a synthetic model of the DO- Tli Kwi Cho kimberlite pipe, Northwest Territories, Canada (Jansen and Witherly, 2004).

Figure 1: A graphical representation of the PGI framework, modified from Astic and Oldenburg (2019). Each diamond box is an optimization process that takes data (shown in rectangular boxes) as well as information provided by the other processes.

2 Background and motivation for a multi-physics PGI framework

In this section, we present the conventions we use for notation, provide background on the univariate PGI framework, and motivate its extension to multiple physical properties.

2.1 Notation conventions

In terms of notation for parameters, we use the following conventions:

  • Lowercase italic symbols are used for scalar values, such as the trade-off parameter .

  • Bold lowercase symbols designate vectors, such as the geophysical model


  • Bold uppercase symbols designate matrices, such as the weight matrix .

For a multi-physics inversion, the geophysical model is likely to contain multiple physical properties. Several surveys might be associated with the same physical property (e.g. gravity and gravity gradiometry) or one survey might depend on several physical properties (e.g. electromagnetic surveys depend on both electrical conductivity and magnetic permeability). We thus adopt the following notations for the geophysical model indices:


A row of the matrix represents all the physical properties that live in the same location. A column represents a single physical property at all the cells of the mesh. For clarity, we are consistent throughout this study with the index notation. The index always refers to the cell number, from to . The vector then denotes all of the physical properties at the th cell:


Likewise, we denote the vector model for a single physical property on the whole mesh with the index with a superscript:


Lastly, the geophysical model at iteration of an inversion is denoted with parentheses .

2.2 The Tikhonov inverse problem and its PGI augmentation

The geophysical inverse problem can be posed as an optimization process where the goal is to find a geophysical model that minimizes an objective function . Using the same formulation of the inverse problem as in Oldenburg and Li (2005), the geophysical optimization problem takes the form:


In equation 5, the vector is the geophysical model, which represents physical properties on a mesh. The term is the geophysical data misfit. The regularization is composed of a smallness term , that penalizes model-values different from a reference model, and smoothness terms , which penalize variations between adjacent cells; those terms are weighted by scaling parameters . The trade-off parameter is a positive scalar that adjusts the relative weighting between the regularization and the data misfit. A value of is sought so that the data misfit is below an acceptable target misfit (Parker, 1977). The scaling parameters and weight the relative importance of the smallness and smoothness terms. In the Tikhonov inversion, first introduced in Tikhonov and Arsenin (1977), each term of the objective function takes a least-squares form. In particular, the smallness term, which is essential in the PGI framework, can be written as:


where is a reference model and the matrix represent local weights.

The PGI approach can be considered as an augmentation of the well-established Tikhonov inversion. This is detailed in Astic and Oldenburg (2019), where we start from a probabilistic formulation of the least-squares inverse problem (Tarantola, 2005), to then include petrophysical and geological information in the form of a GMM. The resulting term is analogous to the smallness, but the reference model and the weights are updated at each iteration. The minimization of the geophysical objective function is labelled Process in the PGI framework (Fig. 1).

2.3 Multivariate GMM: modelling multiple physical properties

To extend the approach presented in Astic and Oldenburg (2019), we represent the petrophysical signature of each geological unit () as a multivariate Gaussian probability distribution, denoted by . The Gaussian function representing the probability distribution of the physical properties of interest for each unit is defined by its mean (vector of size ), and its covariance matrix (matrix of size ), plus its proportion .

The multivariate Gaussian Mixture Model (GMM) simply sums the Gaussian probability distribution representing each known unit, weighted by their proportion:


The variable holds the GMM global variables . With some modifications, the GMM can also represent nonlinear relationships (such as polynomial as presented in Onizawa et al. (2002)). We present those modifications in Appendix A, along with an example of an inversion with various nonlinear relationships between two physical properties.

The geological classification (or membership) is denoted . It is defined as the most probable geological unit, given a set of values for the physical properties:


The categorical variable is key for building a ”quasi-geologic model”

(Li et al., 2019) from the physical properties model obtained by inversion. This corresponds to Process in the PGI framework (Fig. 1).

2.4 Motivation for simultaneously inverting multiple physical properties

Figure 2: Example of a two-dimensional Gaussian mixture model with three rock units. The background is coloured according to the geological classification evaluated by equation 8. A thicker contour line indicates a higher iso-probability density level. On the left and bottom panels, we provide the 1D projections of the total and individual probability distributions for each physical property, and the cumulative histograms of the fictitious samples of each rock unit.

In Fig. 2, we present an example of a GMM with three distinct rock units characterized by two physical properties (two-dimensional distributions). The background is coloured according to the geological identification that would be made at each location using equation 8. The bottom and left panels represent the marginal GMM probability distribution for each physical property individually. We notice that, while all three units are distinct in the two-dimensional space, they can overlap significantly when only considering one physical property at the time. For physical property , rock units and are distinct while rock unit is indistinguishable from rock unit . For physical property , rock units and are now distinct while unit is indistinguishable from either rock units or . This highlights that it is only by jointly inverting for several physical properties that we might be able to uniquely identify three rock units. Units that might not be distinguishable in one survey may be in another one, and by simultaneously working with both physical properties in an inversion, we are able to explore the (or multi)-dimensional physical property space in the centre panel of Fig. 2. In the following section, we show how to use this probability distribution that links the various physical properties as a priori information to regularize the multi-physics inverse problem.

3 Extension of the PGI framework to multi-physics inversions

3.1 Definition of the GMM prior

Astic and Oldenburg (2019) developed a GMM smallness prior to include petrophysical and geological information in inversions involving only one type of geophysical survey with only one physical property to recover. In equation 9, we propose a generalized version of the GMM smallness prior that is designed to couple multiple physical properties and incorporate geological information. Note that now the means are vectors that are the size of the number of different physical properties, and the scalar variance becomes a full positive-definite matrix of the same size. The parameters are both spatially (index ) and lithologically (index ) dependent.



  • is the number of distinct rock units.

  • is the number of active cells in the mesh.

  • represents the physical property values at location .

  • is the a priori probability of observing rock unit at location . It can be either constant over the whole area, then denoted by , or locally determined by a priori geological knowledge.

  • contains the means of the physical properties of rock unit .

  • is the covariance matrix of the physical properties of rock unit .

  • is a weighting term at location , used for example to include depth or sensitivity weighting. We define it from weights ; this is a scalar value for each cell and physical property . is defined as a diagonal matrix made of the combination of all the weights at the specific location : , with the same notation convention as for the model . This allows the weighting to be different for each physical property. We can thus weight each physical property according to the survey on which it depends.

  • holds the GMM global variables .

This GMM probability distribution (equation 9), representing the current geological and petrophysical knowledge, is used to define the smallness term in the regularization. In the next subsection, we use this multivariate GMM prior to develop a modified objective function for the inverse problem.

3.2 The multi-physics PGI geophysical objective function

In Astic and Oldenburg (2019), we demonstrated how to use the negative log-likelihood of a univariate GMM as the smallness term in the Tikhonov inverse problem. The resulting smallness term could be approximated by a least-squares misfit between the current model and a reference model , which was updated at each iteration, as was the smallness matrix . Those dynamic reference model and smallness matrix updates were determined based on the current geophysical model and the geological and petrophysical prior information. The goal of the least-squares approximation was to enable the use of the PGI framework with compiled codes working with the Tikhonov formulation.

Generalizing the result obtained in Astic and Oldenburg (2019) to multiple physical properties, we use the negative log-likelihood of the GMM defined in equation 9 to obtain a single smallness term that couples all of the model parameters. That smallness term can be approximated by the following least-squares misfit:


where is the upper triangular matrix from the Cholesky decomposition of the precision matrix .

Note that our implementation can handle either the log-likelihood of the GMM or its least-squares approximation. When the petrophysical signature of the rock units are not known, it is possible to learn the parameters of the GMM (Astic and Oldenburg, 2019) (Fig. 1, Process ). The extension of that learning process to multiple physical properties can be found in Appendix B.

3.3 Petrophysical target misfit

The PGI smallness expresses a misfit between the petrophysical and geological information and the geophysical model. To measure the goodness of fit and define a stopping criterion for the petrophysical misfit, Astic and Oldenburg (2019) defined a measure and its target value . This measure is similar to the PGI smallness term but without the weights in equation 13. The same approach can be taken here to define the value for the multivariate case:


where and are the same as in equations 11 and 12.

Looking at the term (equation 14) from a probabilistic point of view (Tarantola, 2005; Astic and Oldenburg, 2019), each variable follows a multivariate Gaussian of dimension , with mean , and covariance matrix . Thus, the variable follows a multivariate Gaussian variable with mean and an identity covariance matrix. Thus, each term of the sum in

follows a chi-squared distribution and we can apply Pearson’s chi-squared test

(Pearson, 1900). The target misfit value is defined as the expectation of :


with being the number of active cells in the mesh and being the number of physical properties. This generalizes the result obtained in Astic and Oldenburg (2019) (). This is a similar approach to the definition of a target misfit for the geophysical data as given in Parker (1977).

Our algorithm stops when all of the target misfits, geophysical and petrophysical, are achieved. In the next section, we present our strategy for handling multiple geophysical data misfits as well as an additional petrophysical misfit, each with its target value we seek to reach.

4 Numerical considerations for reaching multiple target misfits

We have multiple geophysical data misfits that we wish to fit. The inclusion of petrophysical and geological data with PGI adds another data misfit term that also needs to reach its target misfit (section 3.3). In this section, we provide our strategies for choosing and dynamically adjusting the various parameters of the objective function to find a solution that fits all the data. An algorithm that summarizes the whole framework is provided in Appendix C.

4.1 Objective function with multiple geophysical data misfits

The objective function we seek to minimize for the multi-physics inversion process (Fig. 1, Process ) takes the form:


The data misfit term now contains multiple geophysical data misfits, each defined as a weighted least-squares norm. is the data misfit of the th survey, where the forward operator generates the predicted data for that survey from . The notation denotes the subset of model parameters associated with the th survey. Note that multiple surveys can be associated with the same physical property, for example, gravity and gravity gradiometry both depend on density contrasts. Similarly, one survey can be sensitive to several physical properties; for example, electromagnetic surveys are sensitive to electrical conductivity and magnetic susceptibility. The data measured by the th survey is symbolized by and the uncertainty on those measurements by the matrix . Each data misfit is weighted by a scaling parameter . Those scaling parameters are important for balancing the various geophysical data misfits and finding a solution that fits all of them. Our approach for updating these parameters is developed later in this section.

The regularization is still composed of the smallness and smoothness terms. The smallness term is our coupling term, which is defined in equation 10. The smoothness terms, one for each direction and physical property, are represented by . In the smoothness terms (equation 19), the smoothness operators (usually first or second-order difference) are represented by the matrix , and weights (sensitivity or depth) are represented by the matrix .

Equation 17 is an intricate objective function that is the sum of many quadratic regularization terms, each of which is multiplied by an adjustable constant. Finding values for these constants and carrying out a nonlinear inversion to produce a model that acceptably fits the data, and is a good candidate for representing information from the true geologic model, is numerically challenging.

In the following sections, we present our approach for estimating and updating these parameters throughout the inversion. The smoothness scaling parameters

are the only values that we keep constant. The scaling parameter weights the importance of the petrophysical misfit term, and the trade-off parameter influences the importance of the regularization (which contains the petrophysical misfit) relative to the geophysical data misfits. The geophysical data misfit scaling parameters are used to adjust the relative importance of each geophysical dataset. Values of , , and are sought so that each geophysical data misfit is below or equal to its target misfit , along with a value of the petrophysical data misfit that is less than or equal to its target misfit .

4.2 The regularization scaling parameters and

Two types of scaling parameters act on the regularization terms; they are the trade-off parameter and the parameters. In our implementation, we keep the parameters acting on the smoothness terms constant while we update and to reach a suitable solution to the PGI problem. Next, we outline our strategies for each of these parameters.

4.2.1 Fixed parameters: the smoothness scaling parameters

In the smoothness terms, the scaling parameters control the relative importance of spatial derivative terms in the regularization. Each set of assigned values will yield different outcomes. This is often a way in which model space can be explored (e.g. preferential smoothness in some directions, see Williams (2008); Lelièvre et al. (2009)). They are generally specified a priori, and we keep them fixed in our objective function throughout the inversion process. In addition to the common practical considerations for choosing the smoothness parameters (Oldenburg and Li, 2005; Williams, 2006), we use the to weight each physical property, by dividing it by the square of its expected maximum amplitude (available through the GMM means if provided). This helps equalize the contribution of the smoothness terms to the objective function value when parameters have widely different scales (like density, log- electrical conductivity and magnetic susceptibility contrasts). The scaling of the physical properties in our extended smallness term (equation 10) is taken care of by the covariance matrices of the GMM.

4.2.2 Adjusted parameters: and

In the case of a single geophysical data misfit with a petrophysical misfit, Astic and Oldenburg (2019) developed a strategy for cooling and warming to find a solution to the inverse problem that reaches the target values of both misfits. This approach is still appropriate in the multi-physics inversion framework and is what we use in this study (step 6 in algorithm 1). For the multi-physics case, we alter it in the following way: when all geophysical data misfits are equal or below their target value, the strategy for warming the scaling parameter on the coupling term is:


4.3 Balancing the geophysical data misfits with the scaling parameters

Our goal is to develop a strategy for scaling multiple geophysical data misfits so that each geophysical dataset is adequately fit. We propose an strategy where the scaling parameters in equation 18

are successively updated. Our approach has a heuristic foundation and does not incur the significant computational cost often associated with optimization-based approaches, and generalize to any number of geophysical data misfits. Before presenting its details, we first outline some strategies that others have taken in addressing this problem.

4.3.1 Review of previous strategies for balancing various geophysical data misfits and coupling terms

Several approaches for weighting multiple geophysical data misfits have been proposed in the recent literature. Some frameworks do not follow any prescribed strategy for updating the scaling parameters of the geophysical data misfits. This is the case for the approaches proposed by Sun and Li (2016, 2017) and Sosa et al. (2013); both keep those scaling parameters constant. Sun and Li use values of unity, while Sosa et al. normalize each geophysical data misfit by its number of data. In our experience, keeping the weights constant has led us to overfit some surveys while underfitting others. Other frameworks have adopted the approach of running their joint inversions for multiple combinations of parameters. For three geophysical data misfits, Moorkamp et al. (2011) adopted a manual check-and-guess approach to adjust the parameters. For two geophysical data misfits, Giraud et al. (2019) ran a subset of their inversions hundreds of times for various combinations of scaling parameters before choosing values based on the L-curve principle (C. and P., 1993; Hansen, 2000; Santos and Bassrei, 2007). They then manually ”fine-tuned” those values using the full joint inversion problem. To avoid the issue of having to choose multiple appropriate scaling parameters, Bijani et al. (2017)

developed a compromise between deterministic and stochastic optimizations for joint inversions. They adopted a ”Pareto Multi-Objective Global Optimization” strategy with genetic algorithms that generate populations of candidate models that ”simultaneously minimize multiple objectives in a Pareto-optimal sense,” rather than working with a fully aggregated objective function. This approach was still computationally expensive and limited to small 2-D studies in the paper.

To limit the number of multiple runs of the same inversion, Lelièvre et al. (2012) devised a rigorously defined, but computationally expensive, strategy for dynamically balancing two geophysical data misfits. Their approach relies on adjusting the trade-off parameter until the two data misfits are Pareto-optimal. Next, they adjust the relative weights of the two surveys to fit both geophysical surveys. They then reinforce the importance of the coupling term before going into another round of adjustments of the trade-off and surveys weights parameters. The approach developed in Astic and Oldenburg (2019) for a single geophysical data misfit, but with a petrophysical data misfit, is related to the work of Lelièvre et al. (2012). Both focus first on fitting the geophysical data misfit terms and then adjusting the coupling term. On the contrary, the strategy presented in Gallardo and Meju (2004) favoured the cross-gradients coupling over the geophysical misfits.

Here we define a practical, computationally inexpensive, heuristic strategy for balancing the geophysical data misfits as well as the coupling term. We design this strategy to work for any number of surveys, and thereby generalize the work of Lelièvre et al. (2012). For the full algorithm, the reader can refer to Appendix C.

4.3.2 Strategy for updating

We now define our strategy for weighting the multiple geophysical data misfits to reach all the target values. We use the scaling parameters defined in equation 18. We dynamically update each geophysical misfit scaling parameter based on its current misfit and target values, compared to the other surveys.

We start with a set of initial scaling parameters that sums to unity. To ensure the progress of all data misfits, while limiting the possibilities of overfitting any given term, we update the scaling parameters during the inversion. Our approach is philosophically similar to what is proposed in Astic and Oldenburg (2019) for and in order to balance the geophysical data misfit and the petrophysical misfit at each iteration in the inversion, and generalized ideas proposed in Lelièvre et al. (2012) to more than two geophysical data misfits. If one geophysical data misfit reaches its target value before the others, we use the ratio of its current value with its target to warm the scaling parameters of the other geophysical data misfit terms. We then normalize the sum of the scaling parameters to be equal to unity again. This is to keep the importance of the total term relatively similar before and after adjusting the scaling parameters . If several surveys are below their respective targets, we simply use the median of the ratios to warm the scaling parameters of the still unfit surveys. Thus, at any iteration of the geophysical inverse problem, if an ensemble of surveys has reached their respective targets, we warm the scaling parameters of the remaining surveys that are not yet fit as (step in algorithm 1):

then we normalize the sum:

An example of convergence curves for the data misfits and evolution of the dynamic scaling parameters is proposed in Fig. 8 for a multi-physics PGI with full petrophysical information.

Our strategy has proven to be insensitive to the initialization of the scaling parameters for linear problems. To demonstrate this point, we show in Fig. 3 the evolution of the scaling parameters for three multi-physics PGI with full petrophysical information. The synthetic example presented in section 6.4 is run with various initializations . The outcomes of all these three PGI were similar to the result we show in Fig. 7. The scaling parameters associated with the magnetic and gravity misfits, respectively, all finish at approximately the same value even though the initializations are very different. The final scaling values are about for the gravity data misfit and around for the magnetic data misfit. This is an appealing property as it reduces the need to fine-tune the initialization of the scaling parameters , which can be costly for large-scale inversions.

Figure 3: Evolution curves of the scaling parameters with the proposed strategy for three multi-physics PGI with full petrophysical information with different initializations for . The color of each line corresponds to the geophysical misfit: blue for gravity and red with markers for magnetic. The style of the lines corresponds to one of the three inversions ( in each inversion).

5 Numerical implementation

We implemented our framework as part of the open-source software

SimPEG (Cockett et al., 2015; Heagy et al., 2017). As such, we are able to share both the software environment and the scripts to reproduce the examples shown in this paper (Astic, 2020). In this section, we highlight some key points of our implementation to encourage the use of this work and future collaborations. A more detailed tutorial is provided in Appendix D, which lays out a pseudocode sketch of the implementation.

SimPEG is designed to be a modular, extensible framework for simulations and inversions of geophysical data. In particular, two features enabled us to focus our implementation efforts on the PGI framework, while using tools provided by the open-source community (such as the forward operators etc.):

  1. the composability of objective functions in the data misfit and the regularization terms,

  2. the use of directives, which orchestrate updates to components of the inversion at each iteration.

The first point enables the implementation of joint inversions. In the code, each misfit term is a Python object that has properties, such as the weights used to construct , and methods, including functions to evaluate the misfit given a model as well as derivatives for use in the optimization routines. To construct a joint inversion with gravity and magnetic data, we first define each misfit term independently and then sum them. We use operator-overloading in Python so that when we express the addition of two objective functions in code, the evaluation of this creates a combo-objective function. This is an object that has the same evaluation and derivative methods as the individual data misfits, and thus readily inter-operates with the rest of the simulation and inversion machinery in SimPEG.

To the second point, directives are functions that are evaluated at the beginning or end of each iteration in the optimization. They are the mechanism we use to make updates to components of the inversion, including the data misfit scaling parameters (equation 21) or smallness weights and reference model (equations 12 and 13), as well as for evaluating the target misfits and stopping criteria for the inversion.

6 Example: The DO- kimberlite pipe

In this section, we illustrate the joint PGI approach on synthetic gravity and magnetic data based on the DO- kimberlite pipe (Jansen and Witherly, 2004), which is composed of two different kimberlite facies. We compare standard Tikhonov inversions of the individual geophysical datasets and independent PGIs of the gravity and magnetic data with the multi-physics PGI approach. Both the Tikhonov inversions and single-physics PGIs produce models that enable only a binary distinction: kimberlite or host rock. Only the multi-physics PGI allows us to identify the two kimberlite facies as distinct from the background host rock.

Scripts and Jupyter notebooks to reproduce the examples presented in this study are available through GitHub at (Astic, 2020).

6.1 Setup

Figure 4: Setup: DO- synthetic example: (a) Plan map, East-West and North-South cross-sections through the synthetic geological model. The dots represent the data locations for the gravity and magnetic survey; (b) Cross-plot and histograms of the physical properties of the synthetic model; (c) Synthetic ground gravity data; (d) Synthetic total amplitude magnetic data.

The DO- Kimberlite pipe (Northwest Territories, Canada) is part of a complex known as the Tli Kwi Cho (TKC) kimberlite cluster (Jansen and Witherly, 2004) (Fig. 4). The pipe has two distinctive kimberlite units that are embedded in a background consisting of a granitic basement covered by a thin layer of till (Fig. 4a). The first pipe unit is a pyroclastic and volcanoclastic kimberlite (called PK/VK), which has a weak magnetic susceptibility and a very high density contrast. The second unit is a hypabyssal kimberlite (called HK), which has a strong magnetic susceptibility and a weak density contrast. The Tikhonov inversions of the field gravity and magnetic datasets have been documented in Devriese et al. (2017).

For this example, we use simulated surface gravity and airborne magnetic data modelled from a synthetic model of the DO-

pipe. The forward and inversion mesh is a tensor mesh with 375 442 active cells; each has a pair of density-magnetic susceptibility values. The smallest cells are cubes with a

m edge length. All chosen values for the surveys and geological units are consistent with observations documented in Devriese et al. (2017). For the PK/VK unit, we assume a magnetic susceptibility of and a density contrast with the background of . For the HK unit, the magnetic susceptibility is set to and the density contrast to (Fig. 4b). We forward modelled the data over a grid of 961 receivers, at the surface for the gravity survey and at the height of m for the airborne magnetic survey (Fig. 4

c and d). We added Gaussian noise to the gravity and magnetic data with standard deviations of

mGal and nT, respectively. These standard deviations are input into the data weighting matrices .

For each inversion, we added bound constraints so that the sought density contrast values are null or negative, and the magnetic susceptibility contrast values are null or positive. We used the sensitivity of each survey to define the weights. Each physical property is weighted by the sensitivity of its associated survey. Sensitivity-based weighting is a common practice for potential fields inversions (Li and Oldenburg, 1996, 1998; Portniaguine and Zhdanov, 2002; Mehanee et al., 2005). The initial model is the background half-space for all inversions.

6.2 Tikhonov inversions

Figure 5: DO- gravity and magnetic Tikhonov inversion results. (a) Plan map, East-West and North-South cross-sections through the recovered density contrast model; (b) Plan map, East-West and North-South cross-sections through the recovered magnetic susceptibility contrast model; (c) Cross-plot of the density and magnetic susceptibility models. The points are coloured using the density and the magnetic susceptibility contrast values (white for the background (BCKGRD), blue for density contrast only, red for magnetic susceptibility contrast only, and purple for co-located significant density and magnetic susceptibility contrasts). The side and bottom panels show the marginal distribution of each physical property, with the best fitting univariate Gaussian (proba. stands for probability, and hist. stands for histogram). Those two univariate Gaussian distributions are used to compute the multivariate Gaussian showed in the background of the cross-plot; (d) Plan map, East-West and North-South cross-sections coloured based on the combination of density and magnetic susceptibility contrasts recovered by Tikhonov inversions (same colourmap as used in (c)).

We first run the individual inversions of the gravity and magnetic data using the well-established Tikhonov approach described in Section 2.2. The results, shown in Fig. 5, are relatively smooth. The gravity inversion (Fig. 5a) provides an approximate outline of the pipe. The magnetic inversion (Fig. 5b) shows a body centred on HK, but it is too diffuse to delineate a shape. Fig. 5(c) shows the cross-plot of the recovered density and magnetic susceptibility models. Each point is coloured based on both its density and magnetic susceptibility values: white when both contrasts are low, with a blue-scale for a significant density contrast only, with a red-scale for only a significant magnetic susceptibility, and with a purple-scale when both contrasts are significant. We observe the expected continuous Gaussian-like distribution of the model parameters (in the region allowed by the bound constraints). Petrophysical signatures are not reproduced, with notably the strongest density contrasts being co-located with the highest magnetic susceptibility values (mesh cells coloured in purple in Fig. 5c and d); this is in contradiction with the setup where PK/VK, the unit with a low density, is distinct from HK, the unit with high magnetic susceptibility. In both gravity and magnetic inversions, the two anomalous units are indistinguishable from each other. To highlight this, we show an overlap of the two inversions in Fig. 5d). We coloured each point relative to its density and magnetic susceptibility values, as in Fig. 5c). This juxtaposition highlights that combining both models does not show structures that seem closer to the ground truth. Post-inversions classification would give highly variable results, depending on the thresholds chosen to delineate units.

We next move to a PGI approach and include petrophysical information. We start by inverting each geophysical dataset individually, and we assess what gains are made before moving to a multi-physics PGI.

6.3 Single-physics PGIs

Figure 6: Results of the individual PGI. (a) Plan map, East-West and North-South cross-sections through the density model recovered using the petrophysical signature of PK/VK; (b) Plan map, East-West and North-South cross-sections through the magnetic susceptibility model obtained using the HK petrophysical signature; (c) Cross-plot of the inverted models. The 2D distribution in the background has been determined by combining the two 1D distributions used for density and magnetic susceptibility PGIs, respectively. With only one anomalous unit in each case, there is still four possible combinations; (d) Plan map, East-West and North-South cross-sections through the quasi-geologic model built from the density and magnetic susceptibility models, see cross-plot in (c).

We apply the PGI framework developed in Astic and Oldenburg (2019) to invert each geophysical dataset individually. The results are shown in Fig. 6. For the prior petrophysical distribution, we use the true value for the means of each unit. For the petrophysical noise levels, we assign standard deviations of of the highest known mean value for each physical property except for the background for which we assign (see the 1D left and bottom panels for the distribution of each physical property, respectively, in Fig. 6c). For the proportions, we also used the true values. We acknowledge that proportion values could be difficult to estimate in practice. However, in our experiments, the values of the global proportions have not had a significant impact on the inversion result. The use of locally varying proportions can, however, guide the reproduction of particular features (Giraud et al., 2017; Astic and Oldenburg, 2019). All the GMM parameters are held fixed in those single-physics inversions.

In carrying out the inversions, we found that both magnetic and gravity data can be explained individually by assuming a single unit, either PK/VK or HK. Each dataset, gravity or magnetic, can be fit by either reproducing the signature of PK/VK or HK, or any value in between. The difference in physical property contrast is compensated for by a difference in the volume of the recovered body. The two kimberlite facies are thus indistinguishable when we consider one geophysical survey at the time. Adding a third cluster in either inversion to represent the second kimberlite facies does not help, as it only gives the algorithm more ”choices” that are not supported by the data. For conciseness, we choose to show here the gravity result recovered using only the petrophysical signature of the PK/VK unit (Fig. 6a), which is the most responsible for the gravity anomaly; for the magnetic inversion, we show the model obtained by only using the petrophysical signature of HK, which is the unit that is the most responsible for the magnetic response (Fig. 6b). The additional models (gravity inversion with HK’s density signature and magnetic inversion with PK/VK’s magnetic signature) are shown in Appendix E; these demonstrate the discrepancy in the recovered volumes of each unit between the magnetic and gravity inversions. For example, explaining the gravity anomaly with only a body with the same density as HK leads to a very large body, bigger than the volume of that same unit recovered through the magnetic inversion. The same reasoning applies to the PK/VK body.

The gravity PGI using the PK/VK unit petrophysical signature (Fig. 6a) gives useful information about the depth and delineation of the pipe that was not available from the Tikhonov inversion. The magnetic PGI using the HK petrophysical signature (Fig. 6b) places a body around the HK unit location but misses its elongated shape. From the petrophysical perspective, both the gravity signature of PK/VK and the magnetic signature of HK are individually well reproduced. However, the combination of the density and magnetic susceptibility contrasts recovered by the individual PGIs (cross-plot in Fig. 6c) is very far from the desired multidimensional petrophysical distributions. Even by assuming just two units for each inversion (background and kimberlite) as we did, there are still four different combinations of density and magnetic susceptibility values. In this specific case, looking at Figs 6c) and d), there is: 1) a cluster representing the background with both weak density and magnetic susceptibility contrasts (coloured in white); 2) a cluster with a large density contrast and a low susceptibility (coloured in blue); this is close to the petrophysical signature of the PK/VK unit; 3) a cluster with high magnetic susceptibility and very small density contrasts (coloured in orange); This would be the HK unit; 4) a cluster that has both high magnetic susceptibility and high-density contrasts, that we identify as ’undefined’ in the figures. This last cluster does not correspond to any unit signature and occupies a large volume. This hinders our ability to resolve two clear kimberlite facies from the inversions. Therefore, this motivates us to move to a multi-physics inversion approach to take advantage of the density-magnetic susceptibility relationships in the inversion and finally delineate two distinct kimberlite facies.

6.4 Multi-physics PGI with petrophysical information

Figure 7: Results of the multi-physics PGI with petrophysical information. (a) Plan map, East-West and North-South cross-sections through the recovered density contrast model; (b) Plan map, East-West and North-South cross-sections through the magnetic susceptibility contrast model; (c) Cross-plot of the inverted models. The colour of the points has been determined by the clustering obtained from this framework joint inversion process. In the background and side panels, we show the prior joint petrophysical distribution with true means we used for this PGI; (d) Plan map, East-West and North-South cross-sections through the resulting quasi-geologic model.

We now apply our multi-physics PGI framework to jointly invert the gravity, magnetic data with petrophysical information (means, covariances, proportions for each unit, background, PK/VK, HK, of the GMM). The parameters are again used to include the appropriate sensitivity weighting for each method and physical property. We use the same uncertainties that we used for the individual PGIs. The off-diagonal elements of the covariance matrices are set to null, which just means we assume no correlations between the density and magnetic susceptibility variations within a single rock unit. The GMM parameters are still held fixed. The results are shown in Fig. 7.

The improvement is significant. The joint inversion succeeds in recovering two distinct kimberlite facies that reproduce the provided petrophysical signatures. The quasi-geologic model (Fig. 7d) is geologically consistent and does not introduce erroneous structures. The surface outline of the pipe is well recovered. The vertical extension is similar to that of the true model. We also now have indications of the elongated shape and tilt of the HK unit. The magnetic susceptibility of HK is slightly overestimated but still within the acceptable margins defined by the petrophysical noise levels we assigned (Fig. 7c).

To illustrate the behaviour of our heuristic strategy for the update of the objective function scaling parameters, we provide in Fig. 8 the convergence curves of the three data misfits (gravity, magnetic, and petrophysical), and the evolution of the various dynamic scaling parameters, for the multi-physics PGI with full petrophysical information. All target values are reached after iterations, and the PGI stops.

Figure 8: Convergence curves for the three misfits, and evolution curves for the dynamic scaling parameters during the multi-physics PGI with petrophysical information shown in Figure 7. (a) Gravity and magnetic geophysical data misfits and their targets (same number of data); (b) Convergence curves of the petrophysical misfit, defined in equation 14, and its target value, defined in equation 16; (c) Evolution of the scaling parameters; (d) Evolution of the scaling parameter; (e) Evolution of the trade-off parameter .

This result was obtained by providing the petrophysical means of the rock units in the GMM. In the next inversion, we devise our approach for using the multi-physics PGI framework when quantitative information is not available.

6.5 Multi-physics PGI with limited information

We have illustrated the gains made by the multi-physics PGI framework when extensive and quantitative a priori information is provided. We now investigate how to perform multi-physics inversions when a priori information to design the coupling term is not available. Astic and Oldenburg (2019)

emphasized the benefits of learning a GMM during the inversion process to compensate for uncertain, or unknown, petrophysical information. At each iteration, the GMM parameters are determined by running a Maximum A Posteriori Expectation-Maximization (MAP-EM) clustering algorithm

(Dempster et al., 1977). The MAP-EM algorithm estimates compromise values for the GMM parameters based on the prior GMM parameters, weighted by confidence parameters in this prior knowledge, and the current geophysical model. We generalize the learning process of the GMM parameters to a multidimensional case in Appendix B.

In the next example, we demonstrate how learning the means of the kimberlite units iteratively through the inversion allows us to still perform multi-physics inversions without providing physical property mean values. This is done by acting on the confidence parameters in the means. Confidence values of zero indicate that the mean is fully learned from the inversions, while infinite confidences fix the mean to its prior value. In all the inversions with limited information, we fix the means of the background for both physical properties to their true values (zero); this is a usual assumption in Tikhonov inversions of potential fields that the background has a zero contrast. We keep the covariances of the GMM fixed and similar to what we used previously. The covariance matrices define our petrophysical noise levels and how spread each petrophysical signature can be.

6.5.1 Employing qualitative information with PGI

Once the Tikhonov inversions have been run (see Fig. 5), it already appears likely that the gravity and magnetic anomalies are mostly generated by two distinct bodies, as the centres of the two recovered anomalous bodies (density and magnetic susceptibility) are at different locations. Lacking quantitative petrophysical information, the multi-physics PGI framework allows us to formulate the following ”interpreter’s assumption”: one kimberlite unit is responsible for the gravity anomaly, while a second one is responsible for the magnetic response. Employing this assumption is made possible in our framework by defining the confidences in the means of each unit as vectors. This allows us to act on each physical property mean value of each unit. For example, for the kimberlite unit that is assumed to be responsible for the magnetic response, we set the confidence in its magnetic susceptibility to zero; the MAP-EM algorithm decides its value at each iteration based solely on the current geophysical model. On the contrary, its mean density contrast is kept fixed at zero by setting the confidence in this mean value to infinity. The same procedure is applied for the kimberlite unit that is assumed to be responsible for the gravity response. The initialization of the density contrast and magnetic susceptibility mean values, for the kimberlite units responsible of the gravity and magnetic response respectively, has little impact on the inversion result, so long as the initial guess is reasonable. It is also common practice, in general, to run clustering algorithms multiple times from various initializations before choosing a specific outcome (Dempster et al., 1977; Murphy, 2012). In the specific result shown in Fig. 9, the density contrast and magnetic susceptibility mean values for the respective kimberlite rock units were initialized at and SI. Similar results were obtained with other initializations ( and SI etc., but not with and SI for all units). Because of the weak dependency of the result with regard to the initialization, we choose not to show the initial value in Fig. 10, which presents the evolution of the means throughout the multi-physics PGI with qualitative information.

Figure 9: Results of the multi-physics PGI without providing the means of the physical properties for the kimberlite facies, and assuming a low-density unit and a magnetized unit; (a) Plan map, East-West and North-South cross-sections through the recovered density contrast model; (b) Plan map, East-West and North-South cross-sections through the magnetic susceptibility contrast model; (c) Cross-plot of the inverted models. The colour of the points has been determined by the clustering obtained from this framework joint inversion process. In the background and side panels, we show the learned petrophysical GMM distribution; (d) Plan map, East-West and North-South cross-sections through the resulting quasi-geologic model.

The result of the multi-physics PGI, with no petrophysical information but with the assumption of distinct low-density and magnetized units, is shown in Fig. 9. Three distinct clusters (background, low-density unit, magnetized unit) are well recovered. The final learned means are respectively for the low-density unit and SI for the magnetized unit.

This multi-physics inversion has several advantages over any of the single-physics inversions. First, by bringing in a qualitative, geologic assumption, we are able to delineate two units by avoiding the overlap of low density and high susceptibility anomalies. Second, we get a sense of the dip of the HK unit. None of those two achievements was reached by the Tikhonov inversions or the single-physics PGIs, even with petrophysical information.

The means of the GMM are learned iteratively following the constraints defined by our assumptions: the background has a fixed contrast of zero in both physical properties, one rock unit is responsible for the gravity response with a null magnetic contrast, and one rock unit is responsible for the magnetic response with a null density contrast. The evolution of the estimations of the means throughout the inversion is shown in Fig. 10. All target values are reached after iterations, and the PGI stops.

Figure 10: Evolution of the learned means of the GMM throughout the multi-physics PGI with qualitative information for the three assumed rock units (background, low-density kimberlite and magnetic kimberlite) shown in Figure 9. The background mean values, the density of the magnetic rock unit, and the magnetic susceptibility of the low-density rock unit are kept fixed. Initialization has a low impact on the learned mean values, and thus the values at iteration 0 are not shown in the plot. (a) Evolution of the density contrast mean values ; (b) Evolution of the magnetic susceptibility mean values.

Finally, we note that distinguishing two bodies, each mostly responsible for a particular geophysical response, is not automatically required from the geophysical datasets themselves. Indeed, we present in Appendix F a multi-physics PGI with no petrophysical information and only two clusters, where the confidences are all zeros for the kimberlite unit; the background mean is kept fixed at zero contrast for both physical properties. A single anomalous body is able to fit both gravity and magnetic datasets, with a learned mean and an acceptable spread according to the set covariance matrices. This further highlights the gains made possible by the opportunities to simply incorporate a qualitative, geologic assumption within the inversion framework.

6.6 Do- example summary

From the petrophysical perspective, the density-magnetic susceptibility cross-plot for the standard Tikhonov inversions is very different from the expected distribution (Fig. 7c). The smoothness of the recovered models and physical property distributions does not allow us to delineate and distinguish between the two kimberlite facies clearly. Using PGI, individual datasets can both be reproduced using a single kimberlite facies. The individual gravity PGI gives us more information about the depth and delineation of the main PK/VK body. The individual magnetic PGI yields a reasonable estimate for the depth of the HK unit but misses its elongated shape. The two individually recovered geologic models are, however, incompatible when they are combined because of the significant overlap of the PK/VK and HK unit. The multi-physics PGI without petrophysical information produces a geological model that distinguishes between the two kimberlite facies. It also begins to give us information about the elongated shape and dip of the HK unit; this result was not achieved by any of the single-physics inversions, not even by the ones that included petrophysical information. However, the accuracy of the boundaries of the bodies is affected by the lack of petrophysical information. The result is improved by providing petrophysical information to the multi-physics PGI, which yields our best recovered model.

7 Discussion

We have expanded the PGI framework developed in Astic and Oldenburg (2019) to carry out joint inversions, and we have proposed a strategy to balance any number of geophysical data misfits along with a coupling term. In our experiments, this strategy appeared to be critical to fitting data from multiple surveys as well as petrophysical data. Finally, we have used a synthetic example to demonstrate the capabilities of the multi-physics PGI framework.

With regards to the iterative learning of GMM means when limited information is available, considering the confidences as vectors is an important contribution of our framework and it advances the approach of Sun and Li (2016). In their approach, the updates to the means are controlled per unit only, without differentiating the physical properties that are well-known from the undocumented ones. They can either learn the means of a unit or keep it fixed, whereas the framework I present is capable of learning specific components of the GMM means for each unit, as demonstrated in the DO- example.

We demonstrated examples of multi-physics inversions with potential fields, which are linear problems. In previous works (Astic and Oldenburg, 2019)

, the PGI approach was applied to nonlinear electromagnetics problems (magnetotelluric, direct-current resistivity, and a field frequency-domain electromagnetics dataset), but it considered only individual surveys depending on a single physical property. We plan to implement this approach for performing multi-physics inversions with electromagnetic methods. This will also be an opportunity to test the robustness of our reweighting strategies for multi-physics inversions with nonlinear geophysical problems. Areas such as the DO-

kimberlite pipe (Devriese et al., 2017; Fournier et al., 2017; Kang et al., 2017), with many different types of geophysical surveys available, are prime candidates for the application of the PGI approach to refine the image of the subsurface structures by integrating more datasets and physical properties in a single inversion.

As we apply the PGI framework to more complex problems, the handling of various types of relationships between physical properties is required. Linear relationships are straight-forward to implement with Gaussian distributions through the covariance matrix, which can define tilted, elongated probability distributions. In Appendix A, we discuss how to account for nonlinear relationships. Such nonlinear relationships are found, for example, between density and seismic velocities (Onizawa et al., 2002). Our framework is flexible enough that different relationships can be included for each rock unit. While modest in size, our example in Fig. 12

is, to our best knowledge, the first one in the literature with such diverse relationships in a single inversion. For the moment, our framework assumes that those nonlinear relationships are given. An interesting avenue of research would be to develop the mathematics for the learning of those nonlinear relationships, along with the other GMM parameters (such as defined in Appendix


Our PGI framework is composed of three regularized optimization problems (Fig. 1). Astic and Oldenburg (2019) laid the mathematical foundation of the framework, with an emphasis on the inclusion and the learning of the petrophysical signature (Process in Fig. 1). In the current study, we focused on the coupling of several geophysical surveys with various physical properties, thus extending the Process in Fig. 1. The third and last process, the geological identification, is still one where there is much room for advancement. Astic and Oldenburg (2019) showed, with a direct-current resistivity example, the efficacy of the proportions when they are locally set to zero or unity. Proportions values of zero or unity are ”constraining,” in the sense that they forbid local occurrences of certain units, rather than just favouring it. While intermediate values of the global proportions appear to have minimal impact on the inversions in our experiments, further studies are required to address the importance and the effects of intermediate (strictly between zero and unity) local proportions. The approach taken by Giraud et al. (2017) combines local proportions from stochastic geological modelling with various warm-started initial models. The integration, extrapolation, and learning of geological information, is part of our current active research. More types of geological information, such as dips, contacts or strikes, also need to be formalized within our framework approach. Combining the PGI smallness term with approaches including prior structural information in the smoothness terms (Fournier and Oldenburg, 2019; Lelièvre and Oldenburg, 2009; Brown et al., 2012; Yan et al., 2017; Giraud et al., 2019) is also to be investigated.

This framework has been implemented as part of the open-source SimPEG (Simulation and Parameter Estimation in Geophysics) project (Cockett et al., 2015). This has two major advantages. First, this enables the reproducibility of the approach by making freely available online, in GitHub repositories, the software environment and Python scripts to recreate the example (, Astic (2020)). Second, by providing a common environment, it allows the implementation of the framework to be readily used with any type of geophysical surveys (such as EM surveys) or discretization (such as OcTree meshes) that are supported in the source code, and facilitate the collaboration with others researchers (Oldenburg et al., 2019).

8 Conclusion

We have expanded the PGI framework to use petrophysical and geological information, represented as a Gaussian Mixture Model, as a coupling term to perform multi-physics joint inversions. We described our strategies for handling multiple geophysical target misfits as well as a petrophysical target misfit. We presented our efforts to make the implementation modular, extensible, and shareable. Finally, we demonstrated, through the DO- kimberlite pipe synthetic example, the gains that can be made by including various types of information into a single inversion. Only a joint approach for inverting the potential field datasets allows us to delineate two kimberlite facies and to reproduce their petrophysical signature.

9 Acknowledgments

We sincerely thank the open-source software SimPEG community whose work has considerably facilitated the research presented here. Special thanks go to Dominique Fournier, for the implementation of the potential fields operators, Seogi Kang, for contributions to the mapping module (Kang et al., 2015) and Rowan Cockett for the pioneering development of SimPEG. We thank Condor Consulting Inc., Peregrine Diamonds Ltd. and Kennecott for making the DO- geologic model available for our research.


  • P. Afnimar, K. Koketsu, and K. Nakagawa (2002) Joint inversion of refraction and gravity data for the three-dimensional topography of a sediment–basement interface. Geophysical Journal International 151 (1), pp. 243–254. External Links: Document Cited by: §1.
  • T. Astic and D. W. Oldenburg (2019) A framework for petrophysically and geologically guided geophysical inversion using a dynamic Gaussian mixture model prior. Geophysical Journal International 219 (3), pp. 1989–2012. External Links: Document,, ISSN 0956-540X, Link Cited by: Appendix B, Appendix B, Appendix B, Appendix B, Appendix B, Figure 1, §1, §2.2, §2.3, §3.1, §3.2, §3.2, §3.2, §3.3, §3.3, §3.3, §4.2.2, §4.3.1, §4.3.2, §6.3, §6.5, §7, §7, §7, Algorithm 1.
  • T. Astic (2020) Collection of scripts for forward modeling and joint inversion of potential fields data. Zenodo. Note: 2020-01-31 External Links: Document, Link Cited by: §5, §6, §7.
  • R. Bijani, P. G. Lelièvre, C. F. Ponte-Neto, and C. G. Farquharson (2017) Physical-property-, lithology- and surface-geometry-based joint inversion using pareto multi-objective global optimization. Geophysical Journal International 209 (2), pp. 730–748. External Links: Document,, ISSN 0956-540X, Link Cited by: §4.3.1.
  • M. Bosch (2004) The optimization approach to lithological tomography: Combining seismic data and petrophysics for porosity prediction. Geophysics 69 (5), pp. 1272–1282. External Links: Document,, Link Cited by: §1.
  • N. Bournas, A. Prikhodko, K. Kwan, J. Legault, V. Polianichko, and S. Treshchev (2018) A new approach for kimberlite exploration using helicopter-borne tdem data. In SEG Technical Program Expanded Abstracts 2018, pp. 1853–1857. External Links: Document,, Link Cited by: §1.
  • V. Brown, K. Key, and S. Singh (2012) Seismically regularized controlled-source electromagnetic inversion. Geophysics 77 (1), pp. E57–E65. External Links: Document,, Link Cited by: §7.
  • H. C. and O. P. (1993) The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems. SIAM Journal on Scientific Computing 14 (6), pp. 1487–17 (English). Note: Copyright - Copyright] © 1993 Society for Industrial and Applied Mathematics; Last updated - 2012-06-29 External Links: ISBN 10648275, Link Cited by: §4.3.1.
  • J. Chen, G. M. Hoversten, D. Vasco, Y. Rubin, and Z. Hou (2007) A Bayesian model for gas saturation estimation using marine seismic AVA and CSEM data. Geophysics 72 (2), pp. WA85–WA95. External Links: Document,, Link Cited by: §1.
  • J. Chen and G. M. Hoversten (2012) Joint inversion of marine seismic AVA and CSEM data using statistical rock-physics models and Markov random fields. Geophysics 77 (1), pp. R65–R80. External Links: Document,, Link Cited by: §1.
  • R. Cockett, S. Kang, L. J. Heagy, A. Pidlisecky, and D. W. Oldenburg (2015) SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications. Computers and Geosciences 85, pp. 142–154. External Links: Document, ISSN 00983004, Link Cited by: §5, §7.
  • M. De Stefano, F. Golfré Andreasi, S. Re, M. Virgilio, and F. F. Snyder (2011) Multiple-domain, simultaneous joint inversion of geophysical data with application to subsalt imaging. Geophysics 76 (3), pp. R69–R80. External Links: Document,, Link Cited by: §1.
  • A. P. Dempster, N. M. Laird, and D. B. Rubin (1977) Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, series B 39 (1), pp. 1–38. Note: With discussion External Links: ISSN 0035-9246, Link Cited by: Appendix B, §6.5.1, §6.5.
  • S. G. R. Devriese, K. Davis, and D. W. Oldenburg (2017) Inversion of airborne geophysics over the DO-27/DO-18 kimberlites — Part 1: Potential fields. Interpretation 5 (3), pp. T299–T311. External Links: Document,, Link Cited by: §1, §6.1, §6.1, §7.
  • J. Doetsch, N. Linde, I. Coscia, S. A. Greenhalgh, and A. G. Green (2010) Zonation for 3D aquifer characterization based on joint inversions of multimethod crosshole geophysical data. Geophysics 75 (6), pp. G53–G64. External Links: Document,, Link Cited by: §1.
  • D. Fournier, S. Kang, M. S. McMillan, and D. W. Oldenburg (2017) Inversion of airborne geophysics over the DO-27/DO-18 kimberlites — Part 2: Electromagnetics. Interpretation 5 (3), pp. T313–T325. External Links: Document,, Link Cited by: §7.
  • D. Fournier and D. W. Oldenburg (2019) Inversion using spatially variable mixed p norms. Geophysical Journal International 218 (1), pp. 268–282. External Links: Document,, ISSN 0956-540X, Link Cited by: §7.
  • L. A. Gallardo and M. A. Meju (2003) Characterization of heterogeneous near-surface materials by joint 2D inversion of dc resistivity and seismic data. Geophysical Research Letters 30 (13). External Links: Document Cited by: §1.
  • L. A. Gallardo and M. A. Meju (2004) Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints. Journal of Geophysical Research: Solid Earth 109 (B3). External Links: Document,, Link Cited by: §4.3.1.
  • L. A. Gallardo and M. A. Meju (2011) Structure-coupled multiphysics imaging in geophysical sciences. Reviews of Geophysics 49 (1). Note: doi: 10.1029/2010RG000330 External Links: Document, ISSN 8755-1209, Link Cited by: §1.
  • J. Giraud, M. Lindsay, V. Ogarko, M. Jessell, R. Martin, and E. Pakyuz-Charrier (2019) Integration of geoscientific uncertainty into geophysical inversion by means of local gradient regularization. Solid Earth 10 (1), pp. 193–210. External Links: Document, Link Cited by: §7.
  • J. Giraud, E. Pakyuz-Charrier, M. Jessell, M. Lindsay, R. Martin, and V. Ogarko (2017) Uncertainty reduction in joint inversion using geologically conditioned petrophysical constraints. Geophysics 82 (6), pp. 1–16. External Links: Document, ISSN 0016-8033, Link Cited by: §1, §1, §6.3, §7.
  • J. Giraud, V. Ogarko, M. Lindsay, E. Pakyuz-Charrier, M. Jessell, and R. Martin (2019) Sensitivity of constrained joint inversions to geological and petrophysical input data uncertainties with posterior geological analysis. Geophysical Journal International 218 (1), pp. 666–688. External Links: Document,, ISSN 0956-540X, Link Cited by: §1, §4.3.1.
  • M. G. D. Giuseppe, A. Troiano, C. Troise, and G. D. Natale (2014) k-Means clustering as tool for multivariate geophysical data analysis. An application to shallow fault zone imaging. Journal of Applied Geophysics 101, pp. 108 – 115. External Links: Document, ISSN 0926-9851, Link Cited by: §1.
  • D. Grana and E. Della Rossa (2010) Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion. Geophysics 75 (3), pp. O21. External Links: Document, ISBN 00168033, ISSN 00168033 Cited by: §1.
  • D. Grana, T. Fjeldstad, and H. Omre (2017) Bayesian Gaussian Mixture Linear Inversion for Geophysical Inverse Problems. Mathematical Geosciences 49 (4), pp. 493–515. External Links: Document, ISSN 1874-8953, Link Cited by: §1.
  • E. Haber and D. Oldenburg (1997) Joint inversion: a structural approach. Inverse Problems 13 (1), pp. 63. External Links: Document, Link Cited by: §1.
  • P. C. Hansen (2000) The L-Curve and its Use in the Numerical Treatment of Inverse Problems. In in Computational Inverse Problems in Electrocardiology, ed. P. Johnston, Advances in Computational Bioengineering, pp. 119–142. Cited by: §4.3.1.
  • L. J. Heagy, R. Cockett, S. Kang, G. K. Rosenkjaer, and D. W. Oldenburg (2017) A framework for simulation and inversion in electromagnetics. Computers & Geosciences 107, pp. 1–19. Cited by: §5.
  • G. M. Hoversten, F. Cassassuce, E. Gasperikova, G. A. Newman, J. Chen, Y. Rubin, Z. Hou, and D. Vasco (2006) Direct reservoir parameter estimation using joint inversion of marine seismic AVA and CSEM data. Geophysics 71 (3), pp. C1–C13. External Links: Document,, ISSN 0016-8033, Link Cited by: §1.
  • J. Jansen and K. Witherly (2004) The Tli Kwi Cho kimberlite complex, Northwest Territories, Canada: A geophysical case study. In 2004 SEG Annual Meeting, Cited by: §1, §6.1, §6.
  • M. D. Jegen, R. W. Hobbs, P. Tarits, and A. Chave (2009) Joint inversion of marine magnetotelluric and gravity data incorporating seismic constraints: Preliminary results of sub-basalt imaging off the Faroe Shelf. Earth and Planetary Science Letters 282 (1-4), pp. 47–55. External Links: Document Cited by: §1.
  • J. Kamm, I. A. Lundin, M. Bastani, M. Sadeghi, and L. B. Pedersen (2015) Joint inversion of gravity, magnetic, and petrophysical data — A case study from a gabbro intrusion in Boden, Sweden. Geophysics 80 (5), pp. B131–B152. External Links: Document,, Link Cited by: §1.
  • S. Kang, R. Cockett, L. J. Heagy, and D. W. Oldenburg (2015) Moving between dimensions in electromagnetic inversions. In SEG Technical Program Expanded Abstracts 2015, pp. 5000–5004. External Links: Document,, Link Cited by: Appendix D, §9.
  • S. Kang, D. Fournier, and D. W. Oldenburg (2017) Inversion of airborne geophysics over the DO-27/DO-18 kimberlites — Part 3: Induced polarization. Interpretation 5 (3), pp. T327–T340. External Links: Document,, Link Cited by: §1, §7.
  • P. Keating and P. Sailhac (2004) Use of the analytic signal to identify magnetic anomalies due to kimberlite pipes. Geophysics 69 (1), pp. 180–190. External Links: Document,, Link Cited by: §1.
  • P. G. Lelièvre, C. G. Farquharson, and C. A. Hurich (2012) Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration. Geophysics 77 (1), pp. K1–K15. External Links: Document,, Link Cited by: §1, §4.3.1, §4.3.1, §4.3.2.
  • P. G. Lelièvre and C. G. Farquharson (2016) Integrated Imaging for Mineral Exploration. In Integrated Imaging of the Earth, pp. 137–166. External Links: Document,, ISBN 9781118929063, Link Cited by: §1.
  • P. G. Lelièvre, D. W. Oldenburg, and N. C. Williams (2009) Integrating geological and geophysical data through advanced constrained inversions. Exploration Geophysics 40 (4), pp. 334–341. External Links: Document, ISBN 0812-3985, ISSN 08123985 Cited by: §4.2.1.
  • P. G. Lelièvre and D. W. Oldenburg (2009) A comprehensive study of including structural orientation information in geophysical inversions. Geophysical Journal International 178 (2), pp. 623–637. External Links: Document,, ISSN 0956-540X, Link Cited by: §7.
  • Y. Li, A. Melo, C. Martinez, and J. Sun (2019) Geology differentiation: A new frontier in quantitative geophysical interpretation in mineral exploration. The Leading Edge 38 (1), pp. 60–66. External Links: Document,, Link Cited by: §1, §2.3.
  • Y. Li and D. W. Oldenburg (1996) 3-D inversion of magnetic data. Geophysics 61 (2), pp. 394–408 (English). External Links: Document, ISBN 0016-8033, ISSN 1070485X Cited by: Appendix A, §6.1.
  • Y. Li and D. W. Oldenburg (1998) 3-D inversion of gravity data. Geophysics 63 (1), pp. 109–119. External Links: Document,, Link Cited by: §6.1.
  • J. Macnae (1995) Applications of geophysics for the detection and exploration of kimberlites and lamproites. Journal of Geochemical Exploration 53 (1), pp. 213 – 243. Note: Diamond Exploration: Into the 21st Century External Links: Document, ISSN 0375-6742, Link Cited by: §1.
  • C. Martinez and Y. Li (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 (2), pp. SL1–SL13. External Links: Document,, Link Cited by: §1.
  • S. Mehanee, N. Golubev, and M. S. Zhdanov (2005) Weighted regularized inversion of magnetotelluric data. In SEG Technical Program Expanded Abstracts 1998, pp. 481–484. External Links: Document,, Link Cited by: §6.1.
  • M. A. Meju and L. A. Gallardo (2016) Structural Coupling Approaches in Integrated Geophysical Imaging. In Integrated Imaging of the Earth, pp. 49–67. External Links: Document,, ISBN 9781118929063, Link Cited by: §1.
  • A. T. Melo, J. Sun, and Y. Li (2017) Geophysical inversions applied to 3D geology characterization of an iron oxide copper-gold deposit in Brazil. Geophysics 82 (5), pp. K1–K13. External Links: Document,, Link Cited by: §1.
  • M. Moorkamp, B. Heincke, M. Jegen, A. W. Roberts, and R. W. Hobbs (2011) A framework for 3-D joint inversion of MT, gravity and seismic refraction data. Geophysical Journal International 184 (1), pp. 477–493. External Links: Document,, ISSN 0956-540X, Link Cited by: §1, §4.3.1.
  • M. Moorkamp, P. G. Lelièvre, N. Linde, and A. Khan (2016) Integrated Imaging of the Earth: Theory and Applications. American Geophysical Union (AGU). External Links: Document,, ISBN 9781118929063, Link Cited by: §1.
  • K. P. Murphy (2012) Machine learning: a probabilistic perspective. The MIT Press. External Links: ISBN 0262018020, 9780262018029 Cited by: Appendix B, §6.5.1.
  • D. W. Oldenburg, L. J. Heagy, S. Kang, and R. Cockett (2019) 3D electromagnetic modelling and inversion: a case for open source. Exploration Geophysics 0 (0), pp. 1–13. External Links: Document,, Link Cited by: §7.
  • D. W. Oldenburg, Y. Li, and R. G. Ellis (1997) Inversion of geophysical data over a copper gold porphyry deposit: A case history for Mt. Milligan. Geophysics 62 (5), pp. 1419–1431. External Links: Document,, Link Cited by: §1.
  • D. W. Oldenburg and Y. Li (2005) Inversion for Applied Geophysics: A Tutorial. In Near-Surface Geophysics, D. K. Butler (Ed.), pp. 89–150. External Links: Document,, Link Cited by: §2.2, §4.2.1.
  • S. Onizawa, H. Mikada, H. Watanabe, and S. Sakashita (2002) A method for simultaneous velocity and density inversion and its application to exploration of subsurface structure beneath Izu-Oshima volcano, Japan. Earth, Planets and Space 54 (8), pp. 803–817. External Links: Document, ISSN 1880-5981, Link Cited by: §2.3, §7.
  • H. Paasche, J. Tronicke, K. Holliger, A. G. Green, and H. Maurer (2006) Integration of diverse physical-property models: Subsurface zonation and petrophysical parameter estimation based on fuzzy c -means cluster analyses. Geophysics 71 (3), pp. H33–H44. External Links: Document,, Link Cited by: §1.
  • H. Paasche and J. Tronicke (2007)

    Cooperative inversion of 2D geophysical data sets: A zonal approach based on fuzzy c-means cluster analysis

    Geophysics 72 (3), pp. A35–A39. External Links: Document,, Link Cited by: §1.
  • H. Paasche (2016) Post-Inversion Integration of Disparate Tomographic Models by Model Structure Analyses. In Integrated Imaging of the Earth, pp. 69–91. External Links: Document,, ISBN 9781118929063, Link Cited by: §1.
  • R. L. Parker (1977) Understanding Inverse Theory. Annual Review of Earth and Planetary Sciences 5 (1), pp. 35–64. External Links: Document,, Link Cited by: §2.2, §3.3.
  • K. Pearson (1900) On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 50 (302), pp. 157–175. External Links: Document,, Link Cited by: §3.3.
  • F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay (2011) Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research 12, pp. 2825–2830. Cited by: Appendix D.
  • O. Portniaguine and M. S. Zhdanov (2002) 3‐D magnetic inversion with data compression and image focusing. Geophysics 67 (5), pp. 1532–1541. External Links: Document,, Link Cited by: §6.1.
  • E.T.F. Santos and A. Bassrei (2007) L- and -curve approaches for the selection of regularization parameter in geophysical diffraction tomography. Computers & Geosciences 33 (5), pp. 618 – 629. External Links: Document, ISSN 0098-3004, Link Cited by: §4.3.1.
  • P. Shamsipour, D. Marcotte, and M. Chouteau (2012) 3D stochastic joint inversion of gravity and magnetic data. Journal of Applied Geophysics 79, pp. 27 – 37. External Links: Document, ISSN 0926-9851, Link Cited by: §1.
  • A. Sosa, A. A. Velasco, L. Velazquez, M. Argaez, and R. Romero (2013) Constrained optimization framework for joint inversion of geophysical data sets. Geophysical Journal International 195 (3), pp. 1745–1762. External Links: Document,, ISSN 0956-540X, Link Cited by: §4.3.1.
  • J. Sun and Y. Li (2015) Multidomain petrophysically constrained inversion and geology differentiation using guided fuzzy c-means clustering. Geophysics 80 (4), pp. ID1. External Links: Document, /gsw/content_public/journal/geophysics/80/4/10.1190_geo2014-0049.1/2/geo2014-0049.pdf, Link Cited by: §1, §1.
  • J. Sun and Y. Li (2016) Joint inversion of multiple geophysical data using guided fuzzy c-means clustering. Geophysics 81 (3), pp. ID37–ID57. External Links: Document, Link Cited by: §1, §1, §4.3.1, §7.
  • J. Sun and Y. Li (2017) Joint inversion of multiple geophysical and petrophysical data using generalized fuzzy clustering algorithms. Geophysical Journal International 208 (2), pp. 1201–1216. External Links: Document, ISSN 1365246X Cited by: §1, §1, §4.3.1.
  • A. Tarantola (2005) Inverse Problem Theory and Methods for Model Parameter Estimation. Vol. 89, siam. External Links: Document,, Link Cited by: §2.2, §3.3.
  • A. N. Tikhonov and V. Y. Arsenin (1977) Solutions of ill-posed problems. V. H. Winston & Sons, Washington, D.C.: John Wiley & Sons, New York. Note: Translated from the Russian, Preface by translation editor Fritz John, Scripta Series in Mathematics External Links: MathReview (M. Z. Nashed) Cited by: §2.2.
  • N. C. Williams (2006) Applying UBC-GIF potential field inversions in greenfields or brownfields exploration. In Proceedings of the Australian Earth Sciences Convention. ASEG/GSA., Cited by: §4.2.1.
  • N. C. Williams (2008) Geologically-constrained ubc-gif gravity and magnetic inversions with examples from the agnew-wiluna greenstone belt, western australia. Ph.D. Thesis, University of British Columbia. External Links: Document, Link Cited by: §4.2.1.
  • P. Yan, T. Kalscheuer, P. Hedin, and M. A. Garcia Juanatey (2017) Two-dimensional magnetotelluric inversion using reflection seismic data as constraints and application in the COSC project. Geophysical Research Letters 44 (8), pp. 3554–3563. External Links: Document,, Link Cited by: §7.

Appendix A Working with nonlinear petrophysical relationships

Figure 11: Gaussian mixture with various polynomial relationships: one linear (no addition required), one quadratic and one cubic. In the main panel, thicker contour lines indicate higher probabilities. The 1D probability distribution for each physical property, and the respective histogram of each unit, are projected on the left and bottom panels.

While linear trends can be accounted through the covariance matrix to obtain elongated clusters, it is also possible to account for nonlinear relationships between physical properties in our framework (Fig. 11). This is achieved by composing the Gaussian function with the nonlinear relationship between the physical properties of the particular rock unit . This corresponds to using for each rock unit in the GMM in equation 9.

In this section, we present a joint inversion of two linear problems whose respective physical properties have nonlinear petrophysical relationships between each other. For simplicity, the physics of the two problems are the same and is based on the examples previously used in Li and Oldenburg (1996). The models are discretized on a 1D mesh defined on the interval and divided into 100 cells. Both datasets consist of data points at various frequencies equally distributed from to evaluated according to:


Each model presents three distinct units. The background for both models is set to while the two other units are linked through a quadratic and a cubic relationship respectively.

We inverted the two datasets using the multi-physics PGI framework, including a priori knowledge of the petrophysical distributions and relationships. The result can be seen in the first row of Figs 12(a) to (c). Note that the petrophysical relationships are well reproduced. This allows the recovery of details of the two models.

For comparison, we also jointly inverted the datasets but without the polynomial relationships (by merely fitting a Gaussian distribution to each unit). The result can be seen in the second row of Figs 12(d) to (f). While the overall structures are well recovered, we miss some details of the models. The background is not as flat as with the full information, the lower tip of the model of problem is completely missed.

We also inverted both datasets individually using the classic Tikhonov inversion. The result is shown in the last row of Figs 12(g) to (i). Results are smoother, as expected. The background presents even more variations, but the overall structures are recovered. Same as for the joint inversion without the polynomial relationships, the details are missing.

Figure 12: Linear joint inversions with various types of physical properties relationships (no trend, quadratic, cubic). Panels (a) to (c) show the inversion result with PGI using the known nonlinear relationships. The first panel (a) shows the result for the first problem, (b) for the second problem. The panel (c) shows the cross-plot of the models over the contour of the GMM with nonlinear relationships. Panels (d) to (f) show the result with PGI without nonlinear relationships. In panel (f), we show the used GMM without nonlinear relationships. Panels (g) to (i) show the Tikhonov inversions result.

Appendix B Updating the Gaussian mixture model in multi-dimensions

In this section, we generalize the learning of the GMM parameters presented in Astic and Oldenburg (2019) to multivariate Gaussian distributions, representing multiple physical properties. The main difference comes from the fact that the means are now vectors (they are scalars in 1-D) and that the scalar variances in 1D become covariance matrices. We define the problem as a Maximum A Posteriori (MAP) estimate of the GMM means, proportions and covariance matrices, using the MAP Expectation-Maximization (MAP-EM) algorithm (Dempster et al., 1977).


We choose to follow a semi-conjugate prior approach for the choice of the prior distributions

(Astic and Oldenburg, 2019; Murphy, 2012).

The computation of the responsibilities for the E-step of the MAP-EM algorithm stay the same as in Astic and Oldenburg (2019), except for the dimension of the parameters:


The update to the proportions stays similar to the univariate case (Astic and Oldenburg, 2019):


where is the confidence in the prior proportion of the rock unit , is the volume of the th cell and is the volume of the active mesh. This allows the estimates to be mesh-independent by using volumetric proportions instead of cell counts.

We generalize the update to the means proposed in Astic and Oldenburg (2019) to each physical property as:


where is the confidence in , which is the prior mean of the physical property of the rock unit . The confidences are thus consider as vectors. This is an important tool that we use in section 6.5.1 to formulate a geologic assumption about the model.

The update to the covariance matrices takes the form:


where is the confidence in the prior covariance matrix of the rock unit .

Appendix C Algorithm

Initialization :

  • Input:

    • Initial geophysical model , GMM parameters and geological model .

  • Parameters:

    • Objective function: data’s noise , and .parameters.

    • Localized prior: specific for available locations , weights .

    • GMM prior weights: for the means, variances and proportions.

    • Optimization: -cooling factor (), sufficient decrease rate (), tolerance on target misfit .

  • Output:

    • , , .

while  and  do

6       543 Objective Function Descent Step:
  • Compute a model perturbation using an inexact Gauss-Newton.

  • Line search with Wolfe condition to find an that satisfy a sufficient decrease of .

  • Return .

Update Petrophysical Distribution
  • Fit a new GMM on such as in equations 27, 30 and 32 until no sufficient increase of the posterior is observed.

  • Compute the membership of the current model as in equation 11 using .

  • Update and according to equations 12 and 13 respectively using .

Update regularization weights:
if  then
7               Decrease : .
8        else if  and  then
9               Increase : (equation 20).
10        if (optional) and and  then
11               Include in Smoothness.
13       12 Update geophysical data misfit weights:
if  then
14               update according to equations 21 to 23.
16 end while
Algorithm 1 PGI extended from Astic and Oldenburg (2019) for joint inversion

Appendix D Pseudo-code for the implementation in SimPEG

Here, we provide an overview for the use of our implementation of the multi-physics PGI framework by other scientists. We outline the main components of the multi-physics PGI implementation using SimPEG and other core tools in the Python ecosystem. As in the paper, we consider the multi-physics inversion of gravity and magnetic data.

We begin by creating a simulation mesh and defining mappings that translate the model vector, which contains all of the parameters we will invert for, to physical properties on the mesh to be used in each forward simulation. The inversion model is a single vector; for an inversion with multiple physical parameters, we stack them and use the Wires map to keep track of which indices in the vector correspond to each physical parameter.


Note that mappings can be composed; for example if we wanted to invert for log-susceptibility rather than susceptibility, then we would compose the wires.susceptibility mapping with an ExpMap instance. For examples, see Kang et al. (2015).

Next, we construct the forward simulations for both the gravity and magnetic data. The survey objects contain the locations of the receivers as well as parameters defining the source field for the magnetics simulation (magnitude, inclination, and declination).


With both gravity and magnetic simulations defined, we now have the ability to compute predicted data given a model. The next step is to construct the data misfit term as described in equation 18.


The mag_data and grav_data objects contain the observed data as well as user-specified uncertainties, and the scalar chi-values are initialized according to equation 23.

Next, we construct the GMM and regularization which consists of the petrophysical smallness term described in equation 10 and smoothness terms for both the density and susceptibility. To construct the GMM, we use the WeightedGaussianMixtureModel object, which inherits from and extends the sklearn.mixture.GaussianMixture object in Scikit-Learn (Pedregosa et al., 2011). The instantiated gmm object has methods for fitting a GMM on a model and computing membership of a model as needed in steps 4 and 5 in Algorithm 1.


Having defined the components of the objective function, we now specify the optimization routine, in this case, inexact Gauss Newton with projections for bound constraints.


Finally, we assemble the inversion using directives in SimPEG to orchestrate updates throughout the inversion. Each directive has an initialize method which is called at the beginning of the inversion and can be used to set initial values of parameters, for example to initialize and the -values, as well as an end_iteration method which is called after a model update (at the end of step 3 in Algorithm 1) and can be used to update parameters in the inversion. The updates in steps 4 through 7 in Algorithm 1 make use of the directives functionality.


Appendix E Additional models obtained by individual PGI for the DO- synthetic case study

Figure 13: (a) magnetic PGI using PK/VK magnetic signature. Note that the recovered volume is much larger than the volume of PK/VK recovered from the gravity inversion with this unit density contrast (Fig. 6c); (b) gravity PGI using the HK density contrast signature. Note again the larger volume compared to Fig. 6(d).

Appendix F Multi-physics PGI with no petrophysical information and no assumption

Figure 14: Results of the multi-physics PGI without providing the means of the physical properties for the kimberlite facies; a single kimberlite body is enough to meet the petrophysical requirements (the spread set by the covariances matrices) and reproduce the geophysical datasets; (a) Plan map, East-West and North-South cross-sections through the recovered density contrast model; (b) Plan map, East-West and North-South cross-sections through the magnetic susceptibility contrast model; (c) Cross-plot of the inverted models. The colour has been determined by the clustering obtained from this framework joint inversion process. In the background and side panels, we show the learned petrophysical GMM distribution; (d) Plan map, East-West and North-South cross-sections through the resulting quasi-geologic model.