1 Introduction
Mineral deposits, or other geologic features, are characterized by different physical properties, and hence multiple geophysical surveys are used to delineate them. For example, kimberlites have signatures that depend upon density, magnetic permeability, and electrical conductivity. They are often discovered through data collected in airborne surveys thanks to their signature of 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; Melo et al., 2017), multiple cases studies have shown that joint 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). However, those inversions require a coupling term that uses a priori information about the physical properties or geologic structures. Coupling methods for joint inversion schemes generally use one or a combination of those two types of information.
The firsts 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 crossgradient. This approach has become commonly used, and Gallardo and Meju (2011) provides an indepth review of the method. However, this strategy assumes that physical properties contrasts are colocated, 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. The framework proposed in this paper belongs to this family. Some of the earliest works used an experimental constitutive formula as their physical properties constraint (Afnimar et al., 2002; Onizawa et al., 2002). Stefano et al. (2011) combined this approach with the above mentioned crossgradient method for subsalt imaging. Some stochastic frameworks have also been proposed. Grana and Della Rossa (2010); Grana et al. (2017) proposed a Bayesian approach for the prediction of lithological facies and fluid properties from seismic attributes. Shamsipour et al. (2012) used geostatistical techniques of cokriging and conditional simulation for jointly inverting gravity and magnetic data assuming the density and magnetic susceptibility auto and crosscovariances follow a linear model of coregionalization. On the deterministic side, of which the framework we present belongs, recent frameworks used clustering techniques such as the fuzzy Cmeans 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. It allows to formulate less strict relationships between physical properties. An important contribution in the last few years is the work by Sun and Li (2015, 2016, 2017) which further developed the method. Besides the addition of the fuzzy Cmeans term, they developed an iterative update to the cluster centres throughout the algorithm; this starts introducing the notion of uncertainty for the petrophysical data. Giraud et al. (2017) focused on reducing uncertainties in stochastic geological modelling by linking geophysics and geological modelling through petrophysical information.
This study extends the framework for petrophysically and geologically guided inversion (PGI) developed in Astic and Oldenburg (2019) as a coupling term to perform joint inversions. This framework (Fig. 1) generalizes concepts presented in Grana et al. (2017), Giraud et al. (2017) and Sun and Li (2017)
. We show how geological and petrophysical information represented as a Gaussian mixture model (GMM) can be used to couple multiple voxelbased geophysical inversions. Each contrasting geological unit is represented by a multidimensional Gaussian distribution which summarizes its characteristics and petrophysical signature. The resulting objective function incorporates both petrophysical and geological information into a single smallness term in the regularization, and the iteration steps are decomposed into a suite of cyclic optimization problems across the geophysical, petrophysical and geological data.
In this paper, we start by reviewing the key concepts of PGI. We then generalize those concepts for the case of multiple physical properties and governing equations for performing joint inversions. We delineate the numerical implementation considerations that are necessary when dealing with multiple datasets and our strategy for reaching multiple target misfits. We also highlight our flexible and modular implementation within the SimPEG framework (Cockett et al., 2015; Heagy et al., 2017). Finally, we demonstrate the advantages of performing joint inversion by using a synthetic model of the DO Tli Kwi Cho kimberlite, Northwest Territories, Canada (Jansen and Witherly, 2004).
2 Key concepts for PGI joint inversion
In this section, we present our representation of petrophysical and geological information as a Gaussian mixture model (GMM) and its inclusion into a geophysical inverse problem.
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 tradeoff parameter .

Bold lowercase symbols designate vectors, such as the geophysical model
. 
Bold uppercase symbols designate matrices, such as the multidimensional Gaussian distribution covariance matrix .
For a joint inversion, the geophysical model is likely to contain multiple physical properties. Several surveys might point to the same physical property (e.g. gravity and gravity gradiometry) or one survey might take several physical properties (e.g. electromagnetic surveys). We thus adopt the following notations for the geophysical model indices:
(1)  
(2) 
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 mesh cells. For clarity, we stay consistent throughout this study with the index notation. The index always refers to the number of model’s parameters, from 1 to . We will then denote a single cell’s physical properties as:
(3) 
Likewise, we denote the vector model for a single physical property on the whole mesh with the index from 1 to with an exponent:
(4) 
Lastly, the geophysical model at iteration of an inversion will be denoted with parentheses .
2.2 Background: the Tikhonov inverse problem
Our inverse problem can be considered as an augmentation of the wellestablished Tikhonov inversion, first introduced in (Tikhonov and Arsenin, 1977). This is labeled process in our PGI framework (Fig. 1). In the Tikhonov approach, the geophysical inverse problem is cast as an optimization problem in which an objective function , such as shown in equation 5, is minimized. Using the same notation convention as Oldenburg and Li (2005), the goal is to find a solution that minimizes:
(5) 
In equation 5, the vector is the geophysical model, which represents physical properties on a mesh. The term is the data misfit, is the model regularization function and is a positive scalar that adjusts the relative weighting between the two terms. A value of is sought so that the data misfit is below an acceptable target misfit (Parker, 1977). The regularization term contains a priori information about the expected features of the geophysical model. It is usually composed of a smallness term that regulates the values the model can take and smoothness terms in all directions which control the variations between adjacent locations, weighted by scalar weights :
(6) 
The smallness term is essential for our framework. Next, we present our representation of petrophysical and geological information and how we include it into the regularization term.
2.3 Gaussian Mixture Model for petrophysical and geological information
To include physical properties information into the inversion, we need to define a way to model the distribution of the physical properties mathematically. Consider a petrophysical dataset with rock samples. For each sample, we have measured physical properties. We will denote those measurements for the various samples (indexed by ) by the vectors . We assume that for each rock sample of the dataset we know its current geological classification, denoted by where is the number of distinct geological units. A sample belonging to unit is noted .
We represent the petrophysical signature of each geological unit as a multidimensional Gaussian probability distribution. If the physical properties of the unit do not follow a Gaussian distribution, we can often project the physical properties into a transformed space where their distribution appears approximately Gaussian, such as a logspace for electrical conductivity or magnetic susceptibility (see the notion of mapping in the SimPEG framework as defined in Cockett et al. (2015); Kang et al. (2015)). The Gaussian distribution representing the physical properties distribution for each unit is defined by its mean (vector of size ), and its covariance matrix (matrix of size ), plus its proportion . Their respective definitions are given in equations 7 to 9, with labeled samples for each unit (total of samples):
(7) 
(8) 
(9) 
The expected probability function to observe an unlabeled physical property data point can be written as a Gaussian Mixture Model (GMM), which simply sums the probability of each known unit:
(10) 
The variable holds the GMM global variables . With some modifications, it can also represent nonlinear relationships (such as polynomial as presented in Onizawa et al. (2002)); see Appendix B for such an example.
We denote
as the categorical variable for the labels, i.e. the geological classification (or membership). Given an unlabeled rock sample with physical properties
, we define its membership as the most probable geological unit:(11)  
with:  
(12)  
(13) 
The categorical variable is key for building a quasigeological model (Li et al., 2019) from the physical properties model obtained by inversion. This corresponds to the process in our framework (Fig. 1).
In Fig. 2, we present an example of a GMM with three distinct rock units characterized by two physical properties (twodimensional distributions). The background is coloured according to the geological identification that would be made at each location using equation 11. On the bottom and the left of the Fig. 2 we represent the GMM probability distribution for each property individually. We notice that, while all three units are distinct in the twodimensional space, they can overlap significantly when only considering one property at the time. For property , the rock units and are distinct while rock unit is indistinguishable from rock unit . For property , rock units and are now distinct while unit is indistinguishable from either rock units or . This highlights the gains that can be made by jointly inverting for several physical properties. Units that might not be distinguishable in one survey may be in another one, and this distinction will guide the inversions towards separating them in both physical properties. 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 multiphysics inverse problem.
2.4 Gaussian mixture model definition
Suppose we are given petrophysical and geological information about the geophysical model that we want to recover in our inversion. Petrophysical information can include an expected number of distinct units, the mean physical property value for each geological unit, or their covariance. Geological information can consist of an expectation of encountering certain rock units in particular locations. We need to design a probability distribution that offers the maximum flexibility for representing that petrophysical and geological a priori knowledge.
In Astic and Oldenburg (2019), we developed a GMM prior to include petrophysical and geological information in inversions involving only one type of geophysical survey with only one physical property to recover. In the equation 14
below, we propose an extended version of that GMM prior that is designed to couple multiple physical properties as well as geological information. Note that now the means are vectors the size of the number of different physical properties and the scalar variance becomes a full positivedefinite matrix of the same size. The parameters are both spatially (index
) and lithologically (index ) dependent.(14) 
where:

is the number of distinct rock units.

is the number of active cells in the mesh.

represents the physical property values at location .

describes the membership to a rock unit.

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 (with different physical properties). 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 weigh each physical property accordingly to the survey on which it depends.

holds the GMM global variables .
This GMM probability distribution (equation 14), representing the current geological and petrophysical knowledge, is used to define the smallness term in the regularization. In the next subsection, we develop the modified objective function for the inverse problem.
2.5 PGI Joint Inversion
In Astic and Oldenburg (2019), we demonstrated how using the GMM defined in equation 14 modifies the smallness term of the Tikhonov inverse problem. The smallness term is still being defined as a leastsquares misfit between the current model, but the reference model is now updated at each iteration, as is the smallness matrix . This dynamic reference model and smallness matrix are determined based on the current geophysical model and the geological and petrophysical prior information. Generalizing the result obtained in Astic and Oldenburg (2019) to multiple physical properties, we obtain a single smallness term that encompasses and couple all of the model parameters:
(15)  
with:  
(16)  
(17)  
(18) 
2.6 Petrophysical and geological target misfit
The PGI smallness expresses a misfit between the petrophysical and geological information and the geophysical model. In Astic and Oldenburg (2019) we defined a measure , similar to the PGI smallness term but without the weights , for which we defined a target misfit . We can generalize the same approach here:
(19)  
with:  
(20)  
(21)  
(22) 
Each term in equation 19 is a multivariate Gaussian variable of dimension (the number of physical properties) of mean with an identity covariance. Thus,
follows a chisquared distribution with a sample size equal to the number of model parameters. Thus its target misfit value is defined as its expectation, such as defined in
Parker (1977) for a geophysical data misfit, and is given by:(23) 
with the number of active cells in the mesh and the number of physical properties. This generalizes the result obtained in Astic and Oldenburg (2019) ().
Our algorithm stops when all the target misfits, geophysical and petrophysical, are fitted. In the next section, we present our approach for handling multiple geophysical data misfits as well as an additional petrophysical misfit.
3 Numerical considerations for reaching multiple target misfits
We have multiple geophysical data misfits that we wish to fit. In addition, the inclusion of petrophysical and geological data with PGI adds another data misfit term; this also needs to reach its target misfit. In this section, we provide our strategies to choose and dynamically adjust 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 D.
3.1 Notations and objective function for multiple geophysical data misfits
We present in the equation below the objective function we are minimizing during the geophysical inversion process (Fig. 1, process 1). The terms contains multiple geophysical data misfits, each defined as a leastsquares norm. The smallness term is our coupling defined in equation 15. The smoothness terms, one for each direction and physical property, are contained within :
(24)  
with:  
(25)  
(26) 
where is the data misfit of the survey and the scalars are weighting terms for each data misfit. These become important to balance the various terms so that a solution is sought such that each is below its target misfit . This is developed later in this section. The forward operator generates the predicted data for the survey from . The notation denotes the subset of model parameters associated with the 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 magneticsusceptibility. The data measured by the survey is symbolized by and the uncertainty on those measurements by the matrix . In the smoothness terms (equation 26), the smoothness operator (usually first or second order difference) is represented by the matrix .
3.2 Review of previous strategies for balancing various geophysical data misfits and coupling terms
In this section, we present strategies that have been previously published with regards to balancing multiple geophysical data misfits along with coupling terms.
In the context of a single geophysical data misfit, we developed, in Astic and Oldenburg (2019), our empirical approach to fit a geophysical data misfit, as well as a petrophysical and geological data misfit. Our approach was based on updating the and values at each iteration. Similar work for linear problems has also been detailed in Sun and Li (2018) where the authors explore several parameters to optimize the clustering of their model, they do not define a target misfit for the petrophysical and geological information.
When several geophysical methods are involved, we need a strategy to scale the multiple geophysical data misfits so that each is properly fitted. This is where the weights in equation 25 become important. Several approaches have been published in the recent literature. Sun and Li (2016) and Sosa et al. (2013) normalize each geophysical survey data misfit by its number of data and keep those weights constant. In our experience, keeping the weights constant has led us to overfit some surveys while underfitting others. On the other hand, Lelièvre et al. (2012) devised a rigorously defined, but computationally expensive, strategy for balancing two geophysical data misfits. Their approach relies on adjusting the tradeoff parameter until the two data misfits are Paretooptimal. They then adjust the relative weights of the two surveys to fit both geophysical surveys. They next reinforce the importance of the coupling term before going into another round of adjustments of the tradeoff and surveys weights parameters. Their approach is somewhat related to what we developed in Astic and Oldenburg (2019), in the sense that they focus first on fitting the geophysical data misfit terms and then adjust the coupling term. On the contrary, the strategy presented in Gallardo and Meju (2004) favoured the coupling over the geophysical misfits.
Here our goal is to define a practical, computationally inexpensive, 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 to not require additional simulations at each parameter update. For the full algorithm, the reader can refer to Appendix D.
3.3 Scaling the regularization terms
Two types of weighting act on the regularization terms; they are the tradeoff parameter and the parameters. In Astic and Oldenburg (2019), we developed a strategy for cooling and warming to find a solution to the inverse problem that fits the geophysical and the petrophysical and geological information. This philosophy is still appropriate in a joint inversion framework and is what we use in this study.
For the smoothness terms, we use the weights to scale each physical property representation by its expected maximum amplitude (available through the GMM means if provided). This helps to equalize their contribution to the objective function value when parameters have widely different scales (like density contrast and magnetic susceptibility contrast). Note that the scaling of the physical properties in our extended smallness term (equation 15) is taken care of by the covariance matrices of the GMM.
3.4 Balancing the geophysical data misfits
We now define our strategy to scale the multiple geophysical data misfits to reach all the target values. We use the weights defined in equation 25. We dynamically update each weight based on its current misfit value as compared to the respective targets for the other surveys.
For initialization, a common normalization is to use the number of data as is done in previous strategies presented above. Here, we also propose a novel initialization. Our motivation is to not only account for the number of data but also for the relative noise levels, the scale of the physical properties, and the numerical operators, which all play major roles in the data misfit metrics. For example, density contrasts are often (not always) within while magneticsusceptibility is most often below SI. As for the associated operators, a given body with a density contrast of can create a gravity response of a few mGal while the same body with a magneticsusceptibility of
SI will produce a magnetic response in the hundreds of nT. There are two steps to our initialization. First, we account for the range of predicted data and noise levels of each geophysical survey, and second, we address the variability in scales of the physical properties. In the first part, for each geophysical data misfit term, we estimate the highest eigenvalue, denoted
, of , with the linearization (or sensitivity matrix, see McGillivray and Oldenburg (1990)) of the physical operator (equation 27). We choose to use a Rayleigh quotient with a power method (Ascher and Greif, 2011), which allows us to get a quick approximation of the . To account for the scale of physical properties, we multiply the highest eigenvalue estimate by the maximum expected value of the associated model parameters (available through the GMM means). The weight is defined as the inverse of the resulting product (equation 28). Because we will update those relative weights later, we prefer to always normalize the sum of the weights to unity (equation 30). This is to keep the importance of the total term relatively similar before and after adjusting the weights.(27)  
(28)  
then we normalize the sum:  (29)  
(30) 
This approach worked well on synthetic examples with linear operators (such as gravity and magnetic). This is what is used in the synthetic example in the present paper.
Good initial weights can help the inversion converge faster, but they have not appeared crucial in our experiments. What ensures the progress of all data misfits, while limiting the possibilities of overfitting, is to update the weights during the inversion. Our approach is philosophically similar to what we use for and to balance the geophysical data misfit and the petrophysical and geological misfit at each iteration in the inversion. 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 weights of the other geophysical data misfit terms and then normalize the sum of the weights to be equal to unity again. If several surveys are below their respective targets, we simply use the median of the ratios to warm the weights of the still unfitted surveys. Thus at any iteration of the geophysical inverse problem, if an ensemble of surveys has reached their respective targets, we warm the weights of the remaining surveys that are not yet fitted as:
(31) 
This strategy is what we used in the case study in Astic et al. (2020) where we invert together three field datasets, including complex petrophysical and geological relationships. The influence of the initial weights has proven minimal compared to the reweighting strategy for the fitting of all datasets.
3.5 A flexible interface for joint inversions in SimPEG
We implemented our framework as part of the opensource software SimPEG (Cockett et al., 2015; Heagy et al., 2017). SimPEG is designed to be a modular, extensible, objectoriented framework for simulations and inversions of geophysical data. In particular, two important features enabled us to focus our implementation efforts:

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

the use of directives which orchestrate updates to components of the inversion at each iteration.
The first point streamlines 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 simply sum them. We use operatoroverloading in Python so that when we express the addition of two objective functions in code, the evaluation of this creates a comboobjective function. This is an object that has the same evaluation and derivative methods as the individual data misfits, and thus readily interoperates with the rest of the simulation and inversion machinery in SimPEG. In pseudocode, the construction of the data misfit as expressed in equation 25 is as follows:
default combodatamisfit.py where the values are scalar weights.
The composability and interoperability extend to the construction of the regularization term. We can build each part of the regularization before combining them:
default comboregularization.py
Having built our framework within the SimPEG comboobjective function also means that our PGI approach can be readily incorporated into any of the geophysical survey types that are supported in the source code, including electrical and electromagnetic methods (Astic and Oldenburg, 2019) as well as magneticvector inversions and gravity gradiometry (Astic et al., 2020). This makes the construction of the objective function for any specific joint problem relatively easy for the user.
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 weights (equation 31) or smallness weights and reference model (equations 17 and 18), as well as for evaluating the target misfits and stopping criteria for the inversion.
Being opensource, working within the SimPEG framework allows us to be transparent about our implementation and to distribute the base code. This also offers us a unique opportunity to provide reproducible examples. Scripts and Jupyter notebooks to reproduce the examples presented in this study are available through Github at https://github.com/simpegresearch/Astic2020JointInversion (Astic, 2020).
4 Example: The DO kimberlite pipe
In this section, we present an illustration of how a joint PGI approach of several geophysical surveys, along with petrophysical data, can improve upon the Tikhonov inversions and the individual PGI results of each geophysical dataset. For that purpose, we compare 3D inversions results of synthetic gravity and magnetic data, based on the DO kimberlite pipe, composed of two different facies. We show that for this case study, the geophysical models obtained from Tikhonov inversions, or even from a PGI approach for each survey, only lead to a binary resolution: host rock or kimberlite pipe. We then demonstrate that using this full joint PGI framework allows distinguishing two clear, distinct kimberlite facies from the background host rock.
4.1 Setup
To illustrate the gains obtained with the joint PGI approach, we use simulated surface gravity and airborne magnetic data modelled from a synthetic model of the DO Kimberlite pipe (Northwest Territories, Canada), part of what was formerly known as the Tli Kwi Cho (TKC) kimberlite deposit (Jansen and Witherly, 2004) (Fig. 3). The pipe has two distinctive kimberlite units of interest that are embedded in a background consisting of a granitic basement covered by a thin layer of till (Fig. 3a). The first 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 an hypabyssal kimberlite (called HK) which has a strong magnetic susceptibility () and a weak density contrast () (Fig. 3b).
The inversion mesh is a tensor mesh with 375 442 active cells; each has a pair of densitymagnetic susceptibility values. The smallest cells are
m in all directions. We forward modelled the data over a grid of 961 receivers at the surface (Fig. 3c and d). We added Gaussian noise to the gravity and magnetic data with standard deviations of
mGal and nT everywhere. These values 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 only, and the magnetic susceptibility contrast values are null or positive. We used the sensitivity of each survey to define the weights. Sensitivitybased weighting is a common practice for potential fields inversions (Li and Oldenburg, 1996, 1998; Portniaguine and Zhdanov, 2002; Mehanee et al., 2005). Note that the sensitivity of the gravity survey is used to weigh the density contrast model, while the sensitivity of the magnetic survey is used to weigh the magneticsusceptibility model. The initial model is the background halfspace for all inversions.
4.2 Tikhonov inversions
We first run the individual inversions of the gravity and magnetic data using the wellestablished Tikhonov approach. The results, shown in Fig. 4, are relatively smooth. The gravity inversion (Fig. 4a) provides an approximate outline of the pipe. The magnetic inversion, on the other hand, shows a magnetic body centred on HK, but it is too diffuse to delineate a shape (Fig. 4b). Looking at the crossplot (Fig. 4c), we observe the expected continuous Gaussianlike distribution of the model parameters (in the region allowed by the bound constraints). Petrophysical signatures are not reproduced, with notably the highest density contrast values being associated with the highest magnetic susceptibility values. In both cases, the two anomalous units are indistinguishable from each other. To highlight this, we show an overlap of the two inversions in Fig. 4(d). We coloured each point relative to its density  magnetic susceptibility values (also shown in Fig. 4c); white for both null values, bluescale for only a density contrast, redscale for a magnetic susceptibility only and a purplescale for both). This highlights that combining both models does not show structures that seem closer to the ground truth. Thresholding would give highly variable results, depending on the values chosen to delineate units.
We then move to a PGI approach, by including petrophysical information, to invert each geophysical dataset individually. We assert what gains are made before moving to a full joint PGI.
4.3 Individual PGI
We apply the framework developed in Astic and Oldenburg (2019) to invert each geophysical dataset individually with the PGI approach. Results are shown in Fig. 5. For the prior petrophysical distribution, we use the true value for the means and proportions of each unit. For the petrophysical noise level, 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 contour plots on the left and the bottom for each physical property respectively in 5d). The GMM parameters are held fixed.
With that approach, we found that both magnetic and gravity data can be explained, individually, by assuming a single unit, either PK/VK or HK. For concision, we choose to show here the gravity result recovered using only the petrophysical signature of the PK/VK unit (Fig. 5a), which is the most responsible for the gravity anomaly; for 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. 5b). The additional models (gravity inversion with HK density signature and magnetic inversion with PK/VK magnetic signature) are shown in Appendix C where we also demonstrate the discrepancy in the recovered volumes of each unit between the magnetic and gravity inversions. The recovered volumes of each unit across inversions are not in agreement. For example, explaining the gravity anomaly with only a body with the same density as HK leads to a very large body (as PK/VK density contrast is four times more important), 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. 5a) gives useful information about the depth and delineation of the pipe, more than the Tikhonov inversion. The magnetic PGI using the HK petrophysical signature (Fig. 5b) places a body around the HK unit location, but miss its elongated shape. From the petrophysical perspective, both the gravity signature of PK/VK and the magnetic signature of HK are individually well reproduced, but not their combination. The densitymagnetic susceptibility crossplot for the individual PGI (Fig. 5c) is very far from the desired petrophysical distributions. Assuming even just one anomalous unit for each inversion as we did, there are still different combinations of physical properties. In this specific case, there is a clear anomalous body that has both the density signature of PK/VK and the magnetic signature of HK (Figs 5c and d), that we identify as ’undefined’ in the figure. It does not correspond to any unit signature and occupies a large volume. This hinders our ability to resolve two clear kimberlites from the inversions. Therefore we move to a joint inversion approach to take advantage of the densitymagnetic susceptibility relationship in the inversion and finally delineate two distinct kimberlite facies.
4.4 Joint PGI
We now apply our extended PGI framework to the joint inversion of the gravity, magnetic and petrophysical data allatonce. The parameters are again used to include the proper sensitivity weighting for each method. We use the same uncertainties that we used for the individual PGI. The nondiagonal elements of the covariance matrices are set to null, which just means we assume no correlations between the density and magneticsusceptibility variations within a single rock unit. The GMM parameters are held fixed. The results are shown in Fig. 6.
The improvement is significant. The joint inversion succeeds in recovering two distinct kimberlite facies that reproduce the provided petrophysical signatures. The surface outline is well recovered but now with two distinct facies, contrary to the individual PGI. The vertical extension is comparable to the truth. We also get a sense of the elongated shape and tilt of the HK unit (Fig. 6d). From the petrophysical perspective, the signature of each unit is well reproduced. The HK magneticsusceptibility of HK is slightly overestimated but still within the acceptable margins defined by the petrophysical noise levels we assigned (Fig. 6c). The discrepancy in the volume of each unit observed for the individual PGI has been resolved by providing the algorithm with the full 2dimensional petrophysical distributions.
4.5 Comment on joint inversions of potential fields data with limited petrophysical information
In Astic and Oldenburg (2019), we emphasized the benefits of learning a GMM during the inversion process to compensate for uncertain or unknown petrophysical information. This allows running our PGI framework without providing physical property mean values. This is performed by the petrophysical characterization process (Fig. 1), which updates the GMM parameters at each iteration based on the current geophysical model. While we generalize the learning process of the GMM parameters to a multidimensional case in Appendix A, we have not made use of it in the present study so far. In our current understanding, it is not as efficient for potential fields as it appeared for electromagnetics inversions. Our hypothesis to explain this is related to the multiple illuminations of the anomalous bodies, in addition to the nonlinearity of the problems, that electromagnetics methods offer compared to potential fields methods. The volume over contrast ratio is the quantity of importance in potential fields. The four individual PGI provide an excellent example of it. Each dataset, gravity or magnetic, could be fitted by either reproducing the signature of PK/VK or HK. The difference in physical property contrast was compensated by a difference in the volume of the recovered body.
We illustrate this point by running a joint PGI without providing physical property mean values nor proportions for the kimberlite units. The means and proportions of the kimberlite units are learned iteratively through the inversion. This is done by setting the and parameters to null for the kimberlite facies (see Appendix A). We fix the mean of the physical properties for the background to its true value (null). We keep the covariances similar to what we used previously. They define our petrophysical noise and how spread each petrophysical signature can be.
We find a volume made of a single kimberlite facies, with physical properties within the prescribed ranges, that fit both the gravity and magnetic data. The resulting geophysical model, learned GMM and geological model are shown in Fig. 7. It thus appears to us that, contrary to the electromagnetics examples provided in Astic and Oldenburg (2019), a good enough guess of the physical properties signature is necessary to perform reasonable PGI of potential fields data. We direct the reader to our DO field case study (Astic et al., 2020) for an example. This aspect will require more thorough research.
4.6 TKC example summary
In all of the inversions, we have shown here, the nearsurface features such as the east tip of the PK/VK unit are relatively well recovered. From the petrophysical perspective, the densitymagnetic susceptibility crossplots for the standard Tikhonov inversions are very removed from the expected distributions represented in contours. Individual PGI give us more information about the depth and delineation of the main PK/VK body. Individual datasets can both be reproduced using a single kimberlite facies. The geological models are, however, incompatible once put together and the petrophysical signature of each unit is only partially recovered. Finally, the joint inversion distinguishes the two units in the recovered models and even starts to give us information about the elongated shape and tilt of the HK unit. A correct estimation of the physical properties contrast appears to be essential for the success of the inversions of potential fields data with PGI.
5 Discussion
We have expanded our PGI framework to carry out joint inversions. We have focused on the methodology and demonstrating the capabilities of the PGI joint framework on an example whose groundtruth is known. A challenging case study that jointly inverts three geophysical surveys (gravity, gravity gradiometry and magneticvector inversion) with five model parameters can be found in Astic et al. (2020)). We have so far produced examples of joint linear problems. In our previous work (Astic and Oldenburg, 2019)
, we applied the PGI approach to multiple electromagnetics nonlinear problems (magnetotelluric, directcurrent resistivity and a field RESOLVE dataset (frequencydomain electromagnetics)). We plan to implement this approach for performing multiphysics inversions with electromagnetic methods.
Oldenburg et al. (2019) demonstrated a joint Tikhonovstyle inversion of several electromagnetics surveys also using the SimPEG package, but without physical property information and with no specific strategy for individually fitting the various data misfits.We 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 multiple surveys as well as petrophysical data. On the contrary, the initial misfit weights played a minor role in comparison. The robustness of our reweighting strategies for joint inversions remains to be tested with nonlinear geophysical problems.
The extension to several physical properties also offers the opportunity to include relationships that link physical properties together. Linear relationships are straightforward to implement with Gaussian distributions through the covariance matrix, which can define tilted, elongated probability distributions. An example of such a linear trend can be found in our case study (Astic et al., 2020). In Appendix B, 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 very modest, our example shown in Fig. 9
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
A). While we think it might be possible for nonlinear geophysical problems with multiple sources, our first assumption is that it might encounter similar issues like the ones mentioned above for the learning of the GMM parameters in potential fields inversion.Our PGI framework is composed of three regularized optimization problems (Fig. 1). In Astic and Oldenburg (2019), we have 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. The next step would be to add geological information, for example from a borehole. This could be used, for example, in our DO study to adjust the depth of the HK unit. The inclusion and dynamic extrapolation of borehole lithology information through the inversion is part of our current research.
6 Conclusion
We have expanded our framework to include petrophysical and geological information as a coupling term to perform multiphysics joint inversions. We have described our strategies for handling multiple geophysical target misfits as well as a petrophysical and geological target misfit. We have presented our efforts to make the implementation modulable, extensible and shareable. Finally, we have demonstrated through the DO 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 fields datasets allows us to delineate two kimberlite facies and to reproduce their petrophysical signature.
7 Acknowledgments
We sincerely thank the opensource 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 his work in the creation and maintenance of SimPEG. We also thank Condor Consulting Inc., Peregrine Diamonds Ltd. and Kennecott for making the DO example available for our research.
References
 Joint inversion of refraction and gravity data for the threedimensional topography of a sediment–basement interface. Geophysical Journal International 151 (1), pp. 243–254. External Links: Document Cited by: §1.
 A first course in numerical methods. edition, Society for Industrial and Applied Mathematics, Philadelphia, PA. External Links: Document, Link, https://epubs.siam.org/doi/pdf/10.1137/9780898719987 Cited by: §3.4.
 Joint inversion of potential fields data over the do27 kimberlite using a gaussian mixture model prior. Submitted to Interpretation. Cited by: §3.4, §3.5, §4.5, §5, §5.
 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: ISSN 0956540X, Document, Link, http://oup.prod.sis.lan/gji/articlepdf/219/3/1989/30144784/ggz389.pdf Cited by: Appendix A, Appendix A, §1, §2.4, §2.5, §2.5, §2.6, §2.6, §3.2, §3.2, §3.3, §3.5, §4.3, §4.5, §4.5, §5, §5, Algorithm 1.
 Collection of scripts for forward modeling and joint inversion of potential fields data. Note: https://doi.org/10.5281/zenodo.3571471Accessed: 20200131 External Links: Document, Link Cited by: §3.5.
 A new approach for kimberlite exploration using helicopterborne tdem data. In SEG Technical Program Expanded Abstracts 2018, pp. 1853–1857. External Links: Document, Link, https://library.seg.org/doi/pdf/10.1190/segam20182996206.1 Cited by: §1.
 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: §1, §2.3, §3.5.
 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 00359246, Link Cited by: Appendix A, Appendix A.
 Inversion of airborne geophysics over the do27/do18 kimberlites  part 1: potential fields. Interpretation 5 (3), pp. T299–T311. External Links: Document, Link, https://doi.org/10.1190/INT20160142.1 Cited by: §1.
 Zonation for 3d aquifer characterization based on joint inversions of multimethod crosshole geophysical data. Geophysics 75 (6), pp. G53–G64. External Links: Document, Link, https://doi.org/10.1190/1.3496476 Cited by: §1.
 Characterization of heterogeneous nearsurface materials by joint 2d inversion of dc resistivity and seismic data. Geophysical Research Letters 30 (13). External Links: Document Cited by: §1.
 Joint twodimensional dc resistivity and seismic travel time inversion with crossgradients constraints. Journal of Geophysical Research: Solid Earth 109 (B3), pp. . External Links: Document, Link, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2003JB002716 Cited by: §3.2.
 Structurecoupled multiphysics imaging in geophysical sciences. Reviews of Geophysics 49 (1). Note: doi: 10.1029/2010RG000330 External Links: Document, ISSN 87551209, Link Cited by: §1.
 Uncertainty reduction in joint inversion using geologically conditioned petrophysical constraints. Geophysics 82 (6), pp. 1–61. External Links: Document, ISSN 00168033, Link Cited by: §1, §1.
 Probabilistic petrophysicalproperties estimation integrating statistical rock physics with seismic inversion. Geophysics 75 (3), pp. O21. External Links: Document, ISBN 00168033, ISSN 00168033 Cited by: §1.
 Bayesian gaussian mixture linear inversion for geophysical inverse problems. Mathematical Geosciences 49 (4), pp. 493–515. External Links: ISSN 18748953, Document, Link Cited by: §1, §1.
 Joint inversion: a structural approach. Inverse Problems 13 (1), pp. 63. External Links: Link, Document Cited by: §1.
 A framework for simulation and inversion in electromagnetics. Computers & Geosciences 107, pp. 1–19. External Links: ISSN 00983004, Document, Link Cited by: §1, §3.5.
 The tli kwi cho kimberlite complex, northwest territories, canada: a geophysical case study. In 2004 SEG Annual Meeting, Cited by: §1, §4.1.
 Joint inversion of marine magnetotelluric and gravity data incorporating seismic constraints: preliminary results of subbasalt imaging off the faroe shelf. Earth and Planetary Science Letters 282 (14), pp. 47–55. External Links: Document Cited by: §1.
 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, https://doi.org/10.1190/geo20140122.1 Cited by: §1.
 Moving between dimensions in electromagnetic inversions. In SEG Technical Program Expanded Abstracts 2015, pp. 5000–5004. External Links: Document, Link, https://library.seg.org/doi/pdf/10.1190/segam20155930379.1 Cited by: §2.3, §7.
 Inversion of airborne geophysics over the do27/do18 kimberlites — part 3: induced polarization. Interpretation 5 (3), pp. T327–T340. External Links: Document, Link, https://doi.org/10.1190/INT20160141.1 Cited by: §1.
 Use of the analytic signal to identify magnetic anomalies due to kimberlite pipes. Geophysics 69 (1), pp. 180–190. External Links: Document, Link, https://doi.org/10.1190/1.1649386 Cited by: §1.
 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, https://doi.org/10.1190/geo20110154.1 Cited by: §1, §3.2.
 3D inversion of magnetic data. Geophysics 61 (2), pp. 394–408 (English). External Links: Document, ISBN 00168033, ISSN 1070485X Cited by: Appendix B, §4.1.
 Geology differentiation: a new frontier in quantitative geophysical interpretation in mineral exploration. The Leading Edge 38 (1), pp. 60–66. External Links: Document, Link, https://doi.org/10.1190/tle38010060.1 Cited by: §2.3.
 3d inversion of gravity data. Geophysics 63 (1), pp. 109–119. External Links: Document, Link, https://doi.org/10.1190/1.1444302 Cited by: §4.1.
 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: ISSN 03756742, Document, Link Cited by: §1.
 Methods for calculating fréchet derivatives and sensitivities for the nonlinear inverse problem: a comparative study. Geophysical Prospecting 38 (5), pp. 499–524. External Links: Document, Link, https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.13652478.1990.tb01859.x Cited by: §3.4.
 Weighted regularized inversion of magnetotelluric data. In SEG Technical Program Expanded Abstracts 1998, pp. 481–484. External Links: Document, Link, https://library.seg.org/doi/pdf/10.1190/1.1820468 Cited by: §4.1.
 Geophysical inversions applied to 3d geology characterization of an iron oxide coppergold deposit in brazil. Geophysics 82 (5), pp. K1–K13. External Links: Document, Link, https://doi.org/10.1190/geo20160490.1 Cited by: §1.
 3D electromagnetic modelling and inversion: a case for open source. Exploration Geophysics 0 (0), pp. 1–13. External Links: Document, Link, https://doi.org/10.1080/08123985.2019.1580118 Cited by: §5.
 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, https://doi.org/10.1190/1.1444246 Cited by: §1.
 Inversion for applied geophysics: a tutorial. In D. K. Butler, ed., NearSurface Geophysics, D. K. Butler (Ed.), pp. 89–150. External Links: Document, Link, https://library.seg.org/doi/pdf/10.1190/1.9781560801719.ch5 Cited by: §2.2.
 A method for simultaneous velocity and density inversion and its application to exploration of subsurface structure beneath izuoshima volcano, japan. Earth, Planets and Space 54 (8), pp. 803–817. External Links: ISSN 18805981, Document, Link Cited by: §1, §2.3, §5.

Cooperative inversion of 2d geophysical data sets: a zonal approach based on fuzzy cmeans cluster analysis
. Geophysics 72 (3), pp. A35–A39. External Links: Document, Link, https://doi.org/10.1190/1.2670341 Cited by: §1.  Understanding inverse theory. Annual Review of Earth and Planetary Sciences 5 (1), pp. 35–64. External Links: Document, Link, https://doi.org/10.1146/annurev.ea.05.050177.000343 Cited by: §2.2, §2.6.
 3‐d magnetic inversion with data compression and image focusing. Geophysics 67 (5), pp. 1532–1541. External Links: Document, Link, https://doi.org/10.1190/1.1512749 Cited by: §4.1.
 3D stochastic joint inversion of gravity and magnetic data. Journal of Applied Geophysics 79, pp. 27 – 37. External Links: ISSN 09269851, Document, Link Cited by: §1.
 Constrained optimization framework for joint inversion of geophysical data sets. Geophysical Journal International 195 (3), pp. 1745–1762. External Links: ISSN 0956540X, Document, Link, http://oup.prod.sis.lan/gji/articlepdf/195/3/1745/1722849/ggt326.pdf Cited by: §3.2.
 Multipledomain, simultaneous joint inversion of geophysical data with application to subsalt imaging. Geophysics 76 (3), pp. R69–R80. External Links: Document, Link, https://doi.org/10.1190/1.3554652 Cited by: §1.
 Multidomain petrophysically constrained inversion and geology differentiation using guided fuzzy cmeans clustering. Geophysics 80 (4), pp. ID1. External Links: Document, Link, /gsw/content_public/journal/geophysics/80/4/10.1190_geo20140049.1/2/geo20140049.pdf Cited by: §1.
 Joint inversion of multiple geophysical data using guided fuzzy cmeans clustering. Geophysics 81 (3), pp. ID37–ID57. External Links: Document, Link, https://doi.org/10.1190/geo20150457.1 Cited by: §1, §3.2.
 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.
 Magnetization clustering inversion — part 1: building an automated numerical optimization algorithm. Geophysics 83 (5), pp. J61–J73. External Links: Document, Link, https://doi.org/10.1190/geo20170844.1 Cited by: §3.2.
 Solutions of illposed 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.
Appendix A Updating the Gaussian mixture model in multidimensions
In this section, we generalize the learning of the GMM parameters presented in Astic and Oldenburg (2019)
to multidimensional Gaussian distributions, representing multiple physical properties at once. The main difference comes from the fact that the means are now vectors (scalar in 1D) and that the scalar variances in 1D become covariance matrices. We define the problem as a Maximum A Posteriori estimate of the GMM means, proportions and covariance matrices:
(32) 
We choose to follow a conjugate prior approach for the choice of the prior distributions
. We use the Maximum A Posteriori ExpectationMaximization (MAPEM) algorithm to estimate the parameters
(Dempster et al., 1977).The computation of the responsibilities for the Estep of the MAPEM algorithm (Dempster et al., 1977) stay the same as in Astic and Oldenburg (2019), except for the dimension of the parameters:
(33) 
The conjugate prior for the proportions follow a Dirichlet distribution:
(34)  
Thus giving the update:  
(35)  
with:  
(36)  
(37) 
where is the volume of the cell and is the volume of the active mesh. This allows the estimates to be meshindependent by using volumetric proportions instead of cell counts.
The semiconjugate for the means and covariance matrices follow respectively a Gaussian and InverseWishart distributions:
(38)  
and:  
(39) 
The update to the means become:
(40)  
with:  
(41) 
The update to the covariance matrices takes the form of:
(42)  
with:  
(43) 
Appendix B Working with nonlinear petrophysical relationships
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. 8). 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 14.
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:
(44) 
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 presented joint PGI framework, including a priori knowledge of the petrophysical distributions and relationships. The result can be seen in the first row of Figs 9(a) to (c). Note that the petrophysical relationships are well reproduced. This allows recovering in great detail 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 9(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 9(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.
Comments
There are no comments yet.