# Disjoint interval bound constraints using the alternating direction method of multipliers for geologically constrained inversion: Application to gravity data

Vitaliy Ogarko, Jérémie Giraud, Roland Martin, Mark Jessell

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

## To cite this version:

Vitaliy Ogarko, Jérémie Giraud, Roland Martin, Mark Jessell. Disjoint interval bound constraints using the alternating direction method of multipliers for geologically constrained inversion: Application to gravity data. Geophysics, 2021, 86 (2), pp.G1-G11. 10.1190/geo2019-0633.1. hal-04904259

HAL Id: hal-04904259

https://ut3-toulouseinp.hal.science/hal-04904259v1

Submitted on 21 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Disjoint interval bound constraints using the alternating direction method of multipliers for geologically constrained inversion: Application to gravity data

Vitaliy Ogarko¹, Jérémie Giraud², Roland Martin³, and Mark Jessell²

## ABSTRACT

To reduce uncertainties in reconstructed images, geologic information must be introduced in a numerically robust and stable way during the geophysical data inversion procedure. In the context of potential (gravity) data inversion, it is important to bound the physical properties by providing probabilistic information on the number of lithologies and ranges of values of possibly existing related rock properties (densities). For this purpose, we have introduced a generalization of bounding constraints for geophysical inversion based on the alternating direction method of multipliers. The flexibility of the proposed technique enables us to take into account petrophysical information as well as probabilistic geologic modeling, when it is available. The algorithm introduces a priori knowledge in terms of physically acceptable bounds of model parameters based on the nature of the modeled lithofacies in the region under study.

Instead of introducing only one interval of geologically acceptable values for each parameter representing a set of rock properties, we define sets of disjoint intervals using the available geologic information. Different sets of intervals are tested, such as quasidiscrete (or narrow) intervals as well as wider intervals provided by geologic information obtained from probabilistic geologic modeling. Narrower intervals can be used as soft constraints encouraging quasidiscrete inversions. The algorithm is first applied to a synthetic 2D case for proof-of-concept validation and then to the 3D inversion of gravity data collected in the Yerrida Basin (Western Australia). Numerical convergence tests show the robustness and stability of the bound constraints that we apply, which is not always trivial for constrained inversions. This technique can be a more reliable uncertainty reduction method as well as an alternative to other petrophysically or geologically constrained inversions based on the more classic "clustering" or Gaussian-mixture approaches.

## INTRODUCTION

One of the main drivers for the usage of prior geologic information in inversion is that it may reduce the severity of nonuniqueness while increasing the geologic plausibility of the recovered model. An obvious and attractive proposal to reach such an end is to force inverted model petrophysical values to lie between geologically realistic bounds, usually using a lower and/or an upper value. Several tech-

niques for maintaining model property values within the given bounds throughout the inversion are available. For instance, Zhang (1997), Zhang and Oldenburg (1999), and Sosa et al. (2013) impose hard bounds on inverted values. Others invert for the logarithm or the square root of the considered property to ensure that the inverted values remain positive (Rijo et al., 1977; Li and Oldenburg, 2003; Newman and Commer, 2005; Lelièvre and Oldenburg, 2006). Alternatively, Kim et al. (1999), Commer (2011), and Zhou et al. (2015) apply model

Manuscript received by the Editor 7 October 2019; revised manuscript received 18 October 2020; published ahead of production 24 December 2020; published online 11 February 2021.

¹The University of Western Australia, The International Centre for Radio Astronomy Research, 7 Fairway, Crawley, Western Australia 6009, Australia and ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D). E-mail: vitaliy.ogarko@uwa.edu.au (corresponding author).

²The University of Western Australia, Centre for Exploration Targeting, School of Earth Sciences, 35 Stirling Highway, Crawley, Western Australia 6009, Australia. E-mail: jeremie.giraud@uwa.edu.au; mark.jessell@uwa.edu.au.

³CNRS, Géoscience Environnement Toulouse, UMR CNRS 5563, Université Paul Sabatier, Observatoire Midi-Pyrénées, 14 Avenue Edouard Belin, Toulouse 34000, France. E-mail: roland.martin@get.omp.eu.

© 2021 The Authors. Published by the Society of Exploration Geophysicists. All article content, except where otherwise noted (including republished material), is licensed under a Creative Commons Attribution 4.0 International License (CC BY). See https://creativecommons.org/licenses/by/4.0/. Distribution or reproduction of this work in whole or in part commercially or noncommercially requires full attribution of the original publication, including its digital object identifier (DOI).
mapping to strictly enforce the inequality constraints. Zhang et al. (2015) implement inequality constraints using Lagrangian multipliers in the inversion's cost function. In Zhdanov and Lin (2017), a multi- nary transformation of model parameters is introduced to approximate discrete inverted model values. More recently, Aghamiry et al. (2019) introduce lower and upper bound constraints using the alternating di- rection method of multipliers (ADMM) (Dykstra, 1983; Boyd et al., 2011; Bertsekas, 2016).

The efficiency of the above techniques is not questioned and has been demonstrated, but, similar to like all constraints applied to geophysical inversion, they have limitations. Although they effectively restrict the solution space, they also limit the space to be explored by inversion before convergence. Bound constraints restricting the model space to a single interval either globally (i.e., the same interval for all elements of the model space) or locally (i.e., a single interval whose bounds vary across the elements of the model space) fail to encompass the full ensemble of geologically plausible models. Consequently, bound constraints might prevent optimum (global) convergence of the algorithm because of the limited model space exploration during the inversion. In particular, when a single interval is used to bind the solution for every model cell, such constraints might prematurely precipitate deterministic inversions in a local minimum. In addition, as hinted above, the bounds within which the inverted values are forced to remain might contradict the true geology of the area, and they might generate numerical instabilities due to the nonlinearity that they introduce in the inverse problem

The objective of this work is to alleviate some of the limitations highlighted above. To this end, we propose an algorithm using multiple bound constraints that allows the inverted model to venture outside such bounds to fit geophysical measurements while maintaining the possibility to enforce them strictly after sufficient exploration of the model space. This is achieved by coupling the well-known least squares with QR-factorization (LSQR) algorithm of Paige and Saunders (1982) with ADMM. In our implementation, model updates are dominated by the contribution of geophysical data and regularization terms until their relative importance has decreased sufficiently for ADMM to progressively guide inversion toward a model honoring the bound constraints as well. We take advantage of the flexibility that ADMM offers and develop the possibility to allow multiple interval bounds determined either by expert knowledge or by geologic modeling for a given model cell. This enables the algorithm to treat the inversion as a supervised clustering problem where the number of clusters and the spread around each of them can be set arbitrarily in every model cell, or according to prior information. We implemented the proposed inversion framework into the open source 3D parallel code Tomofast-x (Martin et al., 2018; Giraud et al., 2019a, 2019b, 2020).

In the remainder of this paper, we first formulate the inverse problem. Then, we detail the use of ADMM to enforce bound constraints and show how results from geologic modeling can be used to set spatially varying bounds. In the ensuing section, we present an application example using data from the Mansfield area (Victoria, Australia). Using this model, we illustrate the capabilities of ADMM described in previous sections and perform a numerical analysis of the inversion. Finally, we place the proposed approach in the general context of geophysical exploration.

METHODOLOGY

Formulation of the inverse problem

We formulate the geophysical inverse problem in the least-squares framework (see chapter 3 in Tarantola, 2005). Let the objective function that we optimize during geophysical inversion be

$$
\theta(\mathbf{d}, \mathbf{m}) = \|\mathbf{d} - g(\mathbf{m})\|_{2}^{2} + \alpha^{2} \|\mathbf{W}(\mathbf{m} - \mathbf{m}^{\text{ref}})\|_{2}^{2}, \quad (1)
$$

where **d** is the observed data and **g**(*m*) is the geophysical response of model *m* = (*m*₁, *m*₂, ..., *m*ₙ) ∈ ℝⁿ, with *n* being the total num- ber of model cells. The second term corresponds to the conventional model damping, with damping weight *α*, depth weighting operator **W**, and reference model *m*^ref. The operator **W** is defined as the square root of the integrated sensitivity matrix (Portnai guine and Zhdanov, 2002). For more details, see Appendix A.

Our aim is to solve equation 1 while encouraging the model to honor geologic bound constraints. This problem is expressed in its generic form as

$$
\text{minimize} \quad \theta(\mathbf{d}, \mathbf{m}) \quad \text{subject to} \quad m_i \in \mathcal{B}_i, 1 \le i \le n, \quad (2)
$$

where $\mathcal{B}_i$ is the interval binding each of the $n$ model cells.

For the ith model cell, $B_i$ is defined as a union of disjointed in-
tervals:

$$
B_i = \bigcup_{l=1}^{L_i} [a_{i,l}, b_{i,l}], \quad \text{with } b_{i,l} > a_{i,l}, \quad \forall l \in [1, L_i], \quad (3)
$$

where $a_{i,l}$ and $b_{i,l}$ are the lower and upper bounds for the $i$th model cell, respectively, and $l$ is the lithology index. Here, $L_i$ is the total number of bounds allowed for the considered cell, corresponding to the number of lithologies.

ADMM in geophysical inversion

We reformulate the geophysical inverse problem in equations 1 and 2 in its ADMM form as

$$
\text{minimize} \quad \theta(\mathbf{d}, \mathbf{m}) + f(z) \quad \text{subject to} \quad \mathbf{m} - z = 0, \quad (4)
$$

where *f* is the indicator function of *B* (for more details, see section 5, equations 5.1 in Boyd et al., 2011). Following Boyd et al. (2011), this problem is solved iteratively by alternating the updates of **m** and *z* as

$$
\mathbf{m}^{k+1} = \arg \min_{\mathbf{m}} (\theta(\mathbf{d}, \mathbf{m}) + \tau^2 \|\mathbf{m} - \mathbf{z}^k + \mathbf{u}^k\|_2^2), \quad (5)
$$

$$
z^{k+1} = \pi_B(\mathbf{m}^{k+1} + \mathbf{u}^k), \quad (6)
$$

$$
\mathbf{u}^{k+1} = \mathbf{u}^k + \mathbf{m}^{k+1} - \mathbf{z}^{k+1}, \quad (7)
$$

where **u** is called the dual variable, $\tau \in \mathbb{R}^+$ is the penalty weight, and $k$ is the current iteration. The ADMM variable $z$ is calculated through the projection onto $\mathcal{B}$:

$$
\pi_B(\mathbf{x}) = [\pi_{B_1}(x_1), \pi_{B_2}(x_2), \dots, \pi_{B_n}(x_n)], \quad \text{with} \quad (8)
$$
$$ \pi_{B_i}(x_i) = \arg \min_{y \in B_i} \|x_i - y\|_2. \qquad (9) $$

The new model $m^{k+1}$ is obtained using the LSQR algorithm of Paige and Saunders (1982). For more details see Appendix A.

The consequence of the formulation of equations 5–7 is that the second term in equation 5 will be increasing at each iteration for the cells where $m \notin B$ and ultimately reaches zero for cells satisfying $m \in B$. This allows early iterations of the inversion to explore the model space outside of the bounds and to gradually encourage the inversion to find geophysically equivalent models honoring the conditions in equation 2. That is, the inversion will seek model updates from the null space (Deal and Nolet, 1996; Muñoz and Rath, 2006) in the vicinity of the current solution that respect the condition set in equation 2. We illustrate the application of projection $\pi_B$ in Figure 1 for three arbitrary intervals.

<figure><img src="image_4.png" /><figcaption>Figure 1. Example of projection using $\pi_B$ per equation 9. The example shown here uses three intervals for a given cell i characterized as follows: [2, 2.2], [4.5, 6], and [8, 9.5].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2. True density contrast model. The densities used are 0, 100, 200, and 300 (kg/m³) for lithologies 4, 3, 2, and 1, respectively.</figcaption></figure>

## Integration with geologic modeling

The set $B$ can be set globally (i.e., homogeneous bound constraints) or in accordance with some knowledge of the local geology (i.e., heterogeneous bound constraints). When prior geologic models are available, they can be used to determine the bounds to be used and their spatial variations. Here, we calculate spatially varying intervals for bound constraints using probabilistic geologic modeling. More specifically, we rely on probabilistic geologic modeling following Pakyuz-Charrier et al. (2018a, 2018b), who perform Monte Carlo sampling of the geologic model space by perturbation of the geologic input data (foliations, dip directions, and contacts). This allows the estimation of lithology observation probabilities. From these, we calculate the bounds of the ith model-cell $B_i$ as

$$ \mathcal{B}_i = \bigcup_{l=1}^{L_i} [a_{i,l}, b_{i,l}], \qquad (10) \\ \Psi_{i,l} > \Psi_{t,l} $$

where $\Psi_{i,l}$ is the observation probability for the ith lithology in the ith model cell. The value chosen for the probability threshold $\Psi_{t,l}$ is used to determine the use of the bound intervals. In the remainder of

<figure><img src="image_6.png" /><figcaption>Figure 3. Probabilistic geologic model (1 refers to 100%). The gray area symbolizes the part of the model where the observation probabilities are equal to zero ($\psi_{i,l} = 0$). The blue dots show the boundary of the region with $\psi_{i,l} > 0$.</figcaption></figure>
this paper, we use $\Psi_{i,l} = 0$ to account for all intervals compatible with geologic modeling.

Note that when the limits of the interval are close to each other, that is, $a_{i,l} = b_{i,l} + \varepsilon$ in equation 10, with $\varepsilon$ being arbitrarily small such that $\varepsilon \ll b_{i,l}$, the inversion algorithm remains stable. We can then perform what we further refer to as *quasidiscrete inversions* (see the example in the next section). In such cases, it can be useful, for instance, to perform geophysical validation or falsification of a series of geologic models produced by geologic modeling engines.

## Convergence and analysis of the constrained problem

The addition of the constraints enforced through ADMM to the LSQR approach does not prevent convergence of the inversion. In particular, the convergence properties of the algorithm are such that it converges for any value of weight $\tau$ used in equation 5 (Ghadimi et al., 2015). The constraints are enforced “only in the limit as the

<figure><img src="image_7.png" /><figcaption>Figure 4. (a) True density model and inversions results for the different uses of ADMM: (b) without bound constraints, (c) with global upper and lower bounds, (d) with global multiple-valued bound constraints without geologic information, and (e) with geologically derived bound constraints. The model consists of 32 × 80 cells.</figcaption></figure>

iterates converge to an optimal solution” (Aster et al., 2019). Convergence of the ADMM is monitored by the value $r$ of the second term of equation 5, unweighted:

$$r = \|\mathbf{m}^k - \mathbf{z}^k + \mathbf{u}^k\|_2^2. \quad (11)$$

We also monitor the geophysical data misfit through the normalized data misfit term in equation 1. Finally, we analyze the results in terms of their model misfit when inverting synthetic data. The model misfit is given as the root-mean-square (rms) of the difference between the true and inverted models.

## MANSFIELD SYNTHETIC SURVEY

### Synthetic model

This subsection introduces synthetic tests to validate the applicability and efficiency of the techniques introduced in the previous section. The synthetic geophysical model is built in the same fashion as in Giraud et al. (2017), of which it constitutes an extension. The geologic structural framework is the same as Pakyuz-Charrier (2018), which is based on field geologic observations from the

<table><caption>Table 1. Inversion metrics for the different cases tested.</caption><thead><tr><th>Bound type</th><th>rms model</th><th>Relative rms data misfit</th><th>ADMM cost</th></tr></thead><tbody><tr><td>(b) No bound constraints</td><td>119.6</td><td>1.27 × 10<sup>-3</sup></td><td>N/A</td></tr><tr><td>(c) Lower and upper</td><td>62.0</td><td>2.8 × 10<sup>-3</sup></td><td>0.41</td></tr><tr><td>(d) Global multiple</td><td>63.5</td><td>3.3 × 10<sup>-3</sup></td><td>0.06</td></tr><tr><td>(e) Geologically derived</td><td>44.5</td><td>4.4 × 10<sup>-3</sup></td><td>0.006</td></tr></tbody></table>

Note: In the fourth column, the cost $r$ is calculated after 50 iterations. The rms model misfit is given in kg/m$^3$.

<figure><img src="image_8.png" /><figcaption>Figure 5. Measured Bouguer gravity anomaly of the Yerrida Basin area. The inset provides a map showing the survey area within Australia.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6. Observation probabilities in the studied areas calculated from probabilistic geologic modeling. The gray dots cover the part of the model where the lithology observation probability is strictly greater than zero. These define parts of the model in which local bound constraints will encourage density contrast values to lie within the corresponding bounds. The model constitutes 100 × 100 × 42 cells.</figcaption></figure>
Mansfield area (Victoria, Australia). We use the work of Pakyuz-Charrier et al. (2018a, 2018b) who perform probabilistic geologic modeling of this area to obtain the probability $\psi_{i,l}$ of occurrence of the different lithologies for every cell of the model. The reference geologic model used to generate the geophysical measurements is shown in Figure 2, and the probabilistic model is shown in Figure 3.

We use the model available for the Mansfield area as follows. We first populate the reference structural model with the true simulated lithology's density contrasts (Figure 2). This density contrast model is used to generate synthetic geophysical measurements. Using the observation probabilities of the different rock types (Figure 3), we derive a set of spatially varying intervals used as bound constraints by applying equation 10.

## Application of different types of bound constraints

Using the model described above, we test several applications of the algorithm presented here. For benchmarking purposes, we first perform inversion without application of bound constraints; that is, $B = (-\infty, +\infty)$ everywhere in the model. We then perform a series of inversions in which we increase the amount of information used to derive the bounds. In the first case, a single interval is used (i.e., $B = B_i = \text{const.}$) across the entire model. This is similar to other bound constraint approaches in which only the upper and lower bounds are provided. Second, spatially invariant disjoint intervals are used (i.e., the same interval $B_i = \cup_{l=1}^r [a_l, b_l]$ is used everywhere in the model). Using this setting, we perform quasidiscrete inversion by choosing a narrow interval. Following this, we use geologic information to derive bounds that vary in space. More specifically, we use heterogeneous bound constraints using equation 10. The bounds are heterogeneous in that the number of intervals and their central value may vary spatially. Using this setting, we test broad and narrow intervals.

The intervals that we use here are defined using a fixed density contrast value and adding $\pm 0.01 \text{ kg/m}^3$. The inversion results of the different cases tested are shown in Figure 4. For quantitative comparison, the corresponding metrics are given in Table 1. For the sake of consistency, we use starting and reference models equal to zero for all of the inversions shown here.

The significant model rms misfit difference between cases (b) and (c) (Table 1) is due to values exceeding $350 \text{ kg/m}^{-3}$ in the bottom right part of the model when no bound constraint is applied (Figure 4b), whereas Figure 4b and 4c shows similar features. The application of multiple bounds (Figure 4d) does not reduce the model misfit (Table 1), but it allows easier geologic interpretation and presents structural features closer to the true model. From Table 1 and Figure 4, it is clear that optimum results are obtained when geologic information is used to adjust the bound constraints locally (Figure 4e). The corresponding inverted model shows the lowest model misfit, and, from a geologic point of view, it is the closest to the geologic model. These synthetic inversion tests and the summarized observations offer an illustration of the potential applicability of the methodology introduced in this paper. We also note that $r$ (equation 11) decreases from case (b) to case (e) (Table 1).

To confirm the robustness of the method to biased geologic information, we tested the presence of the fictitious lithology (for more details, see Appendix B). We also confirmed the method's robustness to using a geologically derived starting model (for more details, see Appendix C).

# YERRIDA BASIN APPLICATION EXAMPLE

## Survey setup and context

The Yerrida Basin has been the subject of recent studies from geologic as well as geophysical viewpoints. This is driven, in part, by the potential presence of buried economic mineral deposits. The Yerrida Basin presents a challenge to mineral explorers because it consists of a sedimentary basin deposited above mafic and felsic formations known to be potential hosts for mineralization pathways. Its formation has been studied by Hawke et al. (2015) and Occhipinti et al. (2017).

Describing the geologic context in detail is beyond the scope of this paper because we focus on the introduction and application of local bound constraints in inversion. A concise introduction to the geology of the area is provided by Giraud et al. (2019a) who study the Yerrida Basin from a geophysical point of view. A more detailed analysis of the area was performed by Lindsay et al. (2020). We therefore provide information essential to this study and leave it to interested readers to refer to relevant references provided in the cited works.

From a geophysical point of view, the complete Bouguer anomaly presents values ranging from approximately -3 to 63 mGal, making the area suitable for gravity modeling. The geophysical measurements inverted for in this study are shown in Figure 5.

From a petrophysical point of view, the area of interest of this work is composed of lithologies that can be differentiated by their density contrasts and grouped accordingly. Density contrasts range from 0 to $330 \text{ kg/m}^3$. As mentioned above, probabilistic geologic modeling was performed by Pakyuz-Charrier et al. (2018a, 2018b) through Monte Carlo sampling of the structural geologic model space associated with the geologic measurements of the area. The resulting probability model for the lithologies associated with the different density contrasts observed in the area is shown in Figure 6. This probability volume was obtained from the simulation of 500 geologic models representative of the geologic measurements as in Lindsay et al. (2020) and made publicly available by Lindsay (2019). We use these probability models to define the bound constraints in the same fashion as shown in the example provided in the previous section. The density contrast used for modeling is given in Table 2.

<table><caption>**Table 2.** Mean density contrast of the different rock formations. These values are used to define intervals with a width of $\pm 0.01 \text{ kg/m}^3$ around the density contrast of the lithologies.</caption><thead><tr><td>Formation (lithology)</td><td>Density contrast (kg/m<sup>3</sup>)</td></tr></thead><tbody><tr><td>Goodin inlier and background lithologies</td><td>0</td></tr><tr><td>Nonmafic greenstone</td><td>30</td></tr><tr><td>Mooloogool group</td><td>130</td></tr><tr><td>Juderina</td><td>180</td></tr><tr><td>Mafic greenstone</td><td>230</td></tr><tr><td>Killara</td><td>330</td></tr></tbody></table>
Inversion results

We follow the same series of inversions as in the previous section to explore some of the possibilities offered by the method that we introduce. We first perform inversion without bound constraints

(Figure 7a). We then perform inversion using a single (global) interval binding the inversion (Figure 7b). Following this, we apply spatially invariant disjoint intervals to quasidiscrete inversion where the values given in Table 2 are enforced across the 3D volume,

<figure><img src="image_10.png" /><figcaption>Figure 7. Inversion results for cases with (a) no bound constraints, (b) lower and upper bound constraints, (c) spatially invariant multiple bounds, and (d) geologically derived multiple bound constraints. The data misfit maps for these models are shown in Appendix D.</figcaption></figure>
assuming no knowledge of the geology. Finally, we use geologic information in the form of the probability volumes shown in Figure 6 to derive bounds that vary in space accordingly with the probability of observation of each of the lithologies. Note that the intervals we use have a narrow width of ± 0.01 kg/m³ around the mean value of the density contrast of the lithologies observed in the studied area. Such narrow intervals were used for the algorithm to encourage inversion to recover values clustered as close as possible to the assumed density contrast of the lithologies observed in the area. Note that all inversions shown here have a posteriori data misfit values with an rms error less than or equal to 1.8 mGal.

The inverted models presented in Figure 7a and 7b show results
that are visually similar except in the regions where the contrast is

close to zero and negative, where the application of the bound con-
straints prevents inversion to produce model cells with density less
than -30 kg/m³. This can also be seen by comparing the histograms
in Figure 8a and 8b. Comparison of Figure 7b and 7c allows a better
understanding of the effect of constraints encouraging inversion to
return a quasidiscrete valued model. That is, the application of
bound constraints with a narrow interval, which produces what
we have introduced as quasidiscrete inversion (Figures 7c, 7d,
8c, and 8d), has a focusing effect and sharpens the model recovered
from inversion while maintaining geophysical consistency. Interest-
ingly, the histogram shown in Figure 8d presents more variability
around the values toward which the geologically derived constraints
drive the model. This variability can be modeled using a mixture
model that may approximate the natural variability of the density

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

Figure 8. Histograms of the recovered density contrast models for the different usage of bound constraints: (a) no bound constraints, (b) lower and upper bound constraints, (c) spatially invariant multiple-value, and (d) geologically derived constraints for quasidiscrete inversion. In panel (d), the blue bold line represents the mixture model matching the histogram. The number of distributions making up the mixture corresponds to the number of units used during modeling L (L = 6) and is a sum of the Gaussian and log-normal distributions. Note that the peaks for (c) and (d) correspond to the same density contrast values.
of rock units in the area (see Figure 8d). This represents the statistical distribution of the petrophysical properties that allows probabilistic geologic modeling and geophysical data to be reconciled. This suggests that Figure 7d corresponds to the recovered model presenting the highest degrees of geologic and petrophysical consistency. It can therefore be used as a starting point for further modeling.

Regarding computational requirements, for the most costly cases (with multiple bounds), the inversion run using 50 major iterations and 100 minor iterations (inside the LSQR solver) took approximately 1000 s using 360 Intel Ivybridge 2.8 GHz CPUs (using message passing interface for parallelization). The additional memory requirements to add bound constraints $M_{ADMM} \approx 4n \times 8$ bytes are small compared to the memory needed to store the full-sensitivity kernel $M_{sens} = N_{data}n \times 8$ bytes (where we are reminded that $n$ is the total number of model cells).

## CONCLUSION

In this work, we have introduced bounded constraints using the ADMM in a 3D parallel inversion framework (using the Tomofast-x code) to guide and constrain the geophysical inversion. We have generalized previous studies by enabling the use of multiple intervals of arbitrary width. These intervals can be determined using petrophysical distributions and/or the available geologic information, thereby restricting the inverted property values to the set lithologies present in the studied area. This allowed us to reduce the uncertainties on the models during the inversion process. Furthermore, intervals can be arbitrarily narrow to allow what we called quasidiscrete inversion. We have trialed the method using synthetic and field data sets. The natural application case that we presented demonstrates the applicability of the technique that we introduced. In particular, it illustrates how we can effectively leverage existing geologic modeling in the area to apply spatially varying sets of the interval as the bound constraints.

This method can be seen as a hybrid method between clustering and hard bound constraints. This can provide information useful to constrain the geology of the area or to rule out the proposed geologic scenarios.

The implementation of the presented algorithm is straightforward because adding this type of constraint is computationally equivalent to adding a well-used model damping term $||m - m^{\text{ref}}||_2^2$ to the cost function. For these reasons, this algorithm also has good scalability properties, as confirmed by runs on a supercomputer using hundreds of CPUs.

The algorithm presents good convergence properties (stability and convergence speed) because it relies on piecewise linear projections instead of the commonly used nonlinear mappings involving the logarithm, exponential, or square root functions.

The ADMM-based bound constraints presented here also can be easily complemented by other constraints relying on clustering based on petrophysical sample distributions (e.g., Gaussian mixtures), or structural similarity constraints-based approaches, or even with geologic probability-based constraints.

Finally, the good scalability of the algorithm that we have developed offers the possibility of testing large domains and a high number of intervals. Furthermore, the method can deal with different geophysical data sets and properties to further reduce model uncertainties.

## ACKNOWLEDGMENTS

Vitaliy Ogarko acknowledges the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) for supporting some of his research efforts. Appreciation is expressed to the CALMIP supercomputing centre (Toulouse, France), for their support through Roland Martin's computing projects no. P1138_2017 and no. P1138_2018 and for the computing time provided on the EOS and Olympe machines. Jeremie Giraud and Mark Jessell are supported, in part, by Loop-3D Enabling Stochastic 3D Geological Modelling (LP170100985) and the Mineral Exploration Cooperative Research Centre (MinEx CRC) whose activities are funded by the Australian Government's Cooperative Research Centre Program. This is MinEx CRC Document 2021/8. The authors also thank Mark D. Lindsay for interesting discussion relating to the utilization of data and model from the Yerrida basin. They also express their appreciation to L. Ailleres and the rest of the Loop consortium's researchers for discussions around geology/geophysics integration.

## DATA AND MATERIALS AVAILABILITY

Data associated with this research are available and can be obtained by contacting the corresponding author.

## APPENDIX A

### INVERSE PROBLEM ALGORITHM

In this study, the least-squares algorithm is used to solve the matrix system. The solution of equation 5 is obtained by an iterative LSQR algorithm of Paige and Saunders (1982). At each kth inversion cycle, we solve the linearized version of the system of equation 5, in the least-squares sense, as follows:

$$ \begin{bmatrix} SW^{-1} \\ \alpha I_d \\ \tau W^{-1} \end{bmatrix} \Delta \bar{m}^{k+1} = \begin{bmatrix} d - g(m^k) \\ -\alpha W(m^k - m^{\text{ref}}) \\ -\tau(m^k - z^k + u^k) \end{bmatrix}, \quad (A-1) $$

$$ \Delta m^{k+1} = W^{-1} \Delta \bar{m}^{k+1}, \quad (A-2) $$

$$ m^{k+1} = m^k + \Delta m^{k+1}, \quad (A-3) $$

where $I_d$ is the identity matrix, $u^0 = 0$, $z^0 = 0$, and the starting model $m^0 = 0$. The updated ADMM variables $z^{k+1}$ and $u^{k+1}$ are calculated using equations 6–7. Here, $S$ is the gravity sensitivity kernel linking the density anomaly model $m$ to the gravity response $g(m)$ such that $g(m) = S$. We precondition $S$, using column scaling, by the depth-weighting operator $W$, to account for gravity signal attenuation with depth. The diagonal operator $W$ is defined for each model parameter $m_i$, over all observed data (Portniaguine and Zhadnov, 2002), as

$$ W_{ii} = \sqrt{\left( \sum_{j=1}^{N_{\text{data}}} S_{ji}^2 \right)^{\frac{1}{2}}}. \quad (A-4) $$
<figure><img src="image_12.png" /></figure>

Figure B-1. Usage of ADMM for geologically derived bound constraints: (a) reminder of results obtained without a fictitious interval (Figure 4e) and (b) inverted model using a fictitious interval centered on a density contrast of 400 kg/m³ added to lithology 2.

In this work, the typical number of iterations used in a major inversion loop is 50 and the number of minor iterations (used in the LSQR solver) is 100.

## APPENDIX B

### BIASED GEOLOGY EXAMPLE

Here, we test the influence of biased geologic information when using spatially varying intervals defined by geology. To this end, we add a fictitious interval around a density contrast equal to 400 kg/m³ everywhere that lithology 2 has a probability of occurrence greater than zero (see Figure 3b). The results of this experiment are shown in Figure B-1b.

As observed in Figure B-1b, only a small portion of the model presents a “wrong” density contrast corresponding to the fictitious interval. This portion coincides to model cells where the true model shows high-density contrasts, that is, 300 kg/m³. It therefore appears that our method is mostly robust to the biased geologic information tested here.

## APPENDIX C

### GEOLOGICALLY DERIVED STARTING MODEL EXAMPLE

For the sake of consistency in the comparison of inversion results using the different bound constraints introduce here, we have used only homogeneous starting models equal to 0 kg/m³. For completeness, below are the results obtained when the starting model is derived from probabilistic geologic modeling (Figure C-1a). It is obtained from the probabilities of observing the different geologic units (Figure 3) multiplied by the corresponding density contrast (Figure 2).

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

Figure C-1. Inversion results using a geologically derived starting model: (a) starting model and (b) inverted model.

Using this starting model, the inverted model, shown in Figure C-1b, has a model misfit of 37 kg/m³ and a data misfit of 3.09 × 10⁻³, which is in line with the values reported in Table 1. Although some geologic units seem to be better resolved, the model presents, overall, features that are similar to Figure 4e. It is worth noting that the data misfit does not evolve significantly during inversion because the starting model produces a value approximating 3 × 10⁻³. However, comparison of the inverted and starting model reveals the changes brought by inversion. Therefore, although such models are geophysically equivalent (similar geophysical data

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

Figure D-1. Comparison of data misfit maps for the inverted models, for the Yerrida Basin example, using different type of bound constraints: (a) no bound constraints, (b) lower and upper bound constraints, (c) spatially invariant multiple bounds, and (d) geologically derived multiple bound constraints.
misfit), the geologic models differ. This means that inversion has navigated the inverse problem's null space from an archetype representative of a category of smooth models (the starting model) to an archetype of rough models represented by the inverted model (for more information about null-space navigation in geophysical problems, see Muñoz and Rath, 2006; Vasco, 2007; De Wit et al., 2012; Fichtner and Zunino, 2019). Importantly, this illustrates the robustness of the method to the starting model, and it shows that it can be applied to field data.

## APPENDIX D
# YERRIDA BASIN DATA MISFIT MAPS

Here, we present the data misfit maps obtained for the cases with no bound constraints, lower and upper bound constraints, spatially invariant multiple bounds, and geologically derived and multiple bound constraints applied to data from the Yerrida Basin (Figure D-1a to D-1d, respectively). These correspond to the model shown in Figure 7.

## REFERENCES

Aghamiry, H. S., A. Gholami, and S. Operto, 2019, Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: Application to large contrast media: Geophysical Journal International, **218**, 855-872, doi: **10.1093/gji/ggz189**.

Aster, R. C., B. Borchers, and C. H. Thurber, 2019, Sparsity regularization and total variation techniques, in Parameter estimation and inverse problems: Elsevier, 181-209.

Bertsekas, D. P., 2016, Nonlinear programming, 3rd ed.: Athena Scientific.

Boyd, S., N. Parikh, E. Chu, B. Peleato, and J. Eckstein, 2011, Distributed optimization and statistical learning via the alternating direction method of multipliers: Foundations and Trends® in Machine Learning, **3**, 1-122, doi: **10.1561/2200000016**.

Commer, M., 2011, Three-dimensional gravity modelling and focusing inversion using rectangular meshes: Geophysical Prospecting, **59**, 966-979, doi: **10.1111/j.1365-2478.2011.00969.x**.

Deal, M. M., and G. Nolet, 1996, Nullspace shuttles: Geophysical Journal International, **124**, 372-380, doi: **10.1111/j.1365-246X.1996.tb07027.x**.

De Wit, R. W. L., J. Trampert, and R. D. Van Der Hilst, 2012, Toward quantifying uncertainty in travel time tomography using the null-space shuttle: Journal of Geophysical Research, Solid Earth, **117**, 1-20, doi: **10.1029/2011JB008754**.

Dykstra, R. L., 1983, An algorithm for restricted least squares regression: Journal of the American Statistical Association, **78**, 837-842, doi: **10.1080/01621459.1983.10477029**.

Fichtner, A., and A. Zunino, 2019, Hamiltonian nullspace shuttles: Geophysical Research Letters, **46**, 644-651, doi: **10.1029/2018GL080931**.

Ghadimi, E., A. Teixeira, I. Shames, and M. Johansson, 2015, Optimal parameter selection for the alternating direction method of multipliers (ADMM): Quadratic problems: IEEE Transactions on Automatic Control, **60**, 644-658, doi: **10.1109/TAC.2014.2354892**.

Giraud, J., M. Lindsay, V. Ogarko, M. Jessell, R. Martin, and E. Pakyuz-Charrier, 2019a, Integration of geoscientific uncertainty into geophysical inversion by means of local gradient regularization: Solid Earth, **10**, 193-210, doi: **10.5194/se-10-193-2019**.

Giraud, J., V. Ogarko, M. Lindsay, E. Pakyuz-Charrier, M. Jessell, and R. Martin, 2019b, Sensitivity of constrained joint inversions to geological and petrophysical input data uncertainties with posterior geological analysis: Geophysical Journal International, **218**, 666-688, doi: **10.1093/gji/ggz152**.

Giraud, J., M. Lindsay, M. Jessell, and V. Ogarko, 2020, Towards plausible lithological classification from geophysical inversion: Honouring geological principles in subsurface imaging: Solid Earth, **11**, 419-436, doi: http://doi.org/10.5194/se-11-419-2020.

Giraud, J., E. Pakyuz-Charrier, M. Jessell, M. Lindsay, R. Martin, and V. Ogarko, 2017, Uncertainty reduction through geologically conditioned petrophysical constraints in joint inversion: Geophysics, **82**, no. 6, ID19-ID34, doi: **10.1190/geo2016-0615.1**.

Hawke, M. L., S. Meffre, H. Stein, P. Hilliard, R. Large, and J. B. Gemmell, 2015, Geochronology of the DeGrussa volcanic-hosted massive sulphide deposit and associated mineralisation of the Yerrida, Bryah and Padbury Basins, Western Australia: Precambrian Research, **267**, 250-284, doi: **10.1016/j.precamres.2015.06.011**.

Kim, H. J., Y. Song, and K. H. Lee, 1999, Inequality constraint in least-squares inversion of geophysical data: Earth, Planets and Space, **51**, 255-259, doi: **10.1186/BF03352229**.

Lelievre, P. G., and D. W. Oldenburg, 2006, Magnetic forward modelling and inversion for high susceptibility: Geophysical Journal International, **166**, 76-90, doi: **10.1111/j.1365-246X.2006.02964.x**.

Li, Y., and D. W. Oldenburg, 2003, Fast inversion of large-scale magnetic data using wavelet transforms: Geophysical Journal International, **152**, 251-265, doi: **10.1046/j.1365-246X.2003.01766.x**.

Lindsay, M., S. Occhipinti, C. Laflamme, A. Aitken, and L. Ramos, 2020, Mapping undercover: Integrated geoscientific interpretation and 3D modelling of a Proterozoic basin: Solid Earth, **11**, 1053-1077, doi: **10.5194/se-11-1053-2020**.

Lindsay, M. D., 2019, Yerrida basin 3D geological model and gravity inversion results (Version 1) [Data set]: Zenodo.

Martin, R., V. Ogarko, D. Komatitsch, and M. Jessell, 2018, Parallel three-dimensional electrical capacitance data imaging using a nonlinear inversion algorithm and Lp norm-based model: Measurement, **128**, 428-445, doi: **10.1016/j.measurement.2018.05.099**.

Muñoz, G., and V. Rath, 2006, Beyond smooth inversion: The use of null-space projection for the exploration of non-uniqueness in MT: Geophysical Journal International, **164**, 301-311, doi: **10.1111/j.1365-246X.2005.02825.x**.

Newman, G. A., and M. Commer, 2005, New advances in three dimensional transient electromagnetic inversion: Geophysical Journal International, **160**, 5-32, doi: **10.1111/j.1365-246X.2004.02468.x**.

Occhipinti, S., R. Hocking, M. Lindsay, A. Aitken, I. Copp, J. Jones, S. Sheppard, F. Pirajno, and V. Metelka, 2017, Paleoproterozoic basin development on the northern Yilgarn Craton, Western Australia: Precambrian Research, **300**, 121-140, doi: **10.1016/j.precamres.2017.08.003**.

Paige, C. C., and M. A. Saunders, 1982, LSQR: An algorithm for sparse linear equations and sparse least squares: ACM Transactions on Mathematical Software, **8**, 43-71, doi: **10.1145/355984.355989**.

Pakyuz-Charrier, E., J. Giraud, V. Ogarko, M. Lindsay, and M. Jessell, 2018a, Drillhole uncertainty propagation for three-dimensional geological modeling using Monte Carlo: Tectonophysics, **747-748**, 16-39, doi: **10.1016/j.tecto.2018.09.005**.

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

Pakyuz-Charrier, E. I. G., 2018, Mansfield (Victoria, Australia) area original GeoModeller model and relevant MCUE outputs: Zenodo.

Portniaguine, O., and M. S. Zhdanov, 2002, 3-D magnetic inversion with data compression and image focusing: Geophysics, **67**, 1532-1541, doi: **10.1190/1.1512749**.

Rijo, L., W. H. Pelton, E. C. Feitosa, and S. H. Ward, 1977, Interpretation of apparent resistivity data from Apodi valley, Rio Grande do Norte, Brazil: Geophysics, **42**, 811-822, doi: **10.1190/1.1440749**.

Sosa, A., 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**, 1745-1762, doi: **10.1093/gji/ggv326**.

Tarantola, A., 2005, Inverse problem theory and methods for model parameter estimation: SIAM.

Vasco, D. W., 2007, Invariance, groups, and non-uniqueness: The discrete case: Geophysical Journal International, **168**, 473-490, doi: **10.1111/j.1365-246X.2006.03161.x**.

Zhang, Y., J. Yan, F. Li, C. Chen, B. Mei, S. Jin, and J. H. Dohm, 2015, A new bound constraints method for 3-D potential field data inversion using Lagrangian multipliers: Geophysical Journal International, **201**, 267-275, doi: **10.1093/gji/ggv016**.

Zhang, Z., 1997, Reconstruction of conductivity and susceptibility from the inversion of EM data: Ph.D. thesis, University of British Columbia.

Zhang, Z., and D. W. Oldenburg, 1999, Simultaneous reconstruction of 1-D susceptibility and conductivity from electromagnetic data: Geophysics, **64**, 33-47, doi: **10.1190/1.1444528**.

Zhdanov, M. S., and W. Lin, 2017, Adaptive multinary inversion of gravity and gravity gradiometry data: Geophysics, **82**, no. 6, G101-G114, doi: **10.1190/geo2016-0451.1**.

Zhou, J., X. Meng, L. Guo, and S. Zhang, 2015, Three-dimensional cross-gradient joint inversion of gravity and normalized magnetic source strength data in the presence of remanent magnetization: Journal of Applied Geophysics, **119**, 51-60, doi: **10.1016/j.jappgeo.2015.05.001**.

Biographies and photographs of the authors are not available.