# Joint Inversion of Electromagnetic and Direct Current Resistivity Data Using Trust Regions. Application to Uranium Exploration in the Athabasca Basin

Ali Mohand-Said, Guy Marquis, Serge Sambolian, Jean-françois Girard, Grant Harrison, Elodie Williard

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

## To cite this version:

Ali Mohand-Said, Guy Marquis, Serge Sambolian, Jean-françois Girard, Grant Harrison, et al., Joint Inversion of Electromagnetic and Direct Current Resistivity Data Using Trust Regions. Application to Uranium Exploration in the Athabasca Basin. Geophysical Prospecting, 2025, 73 (8), pp.e70079. 10.1111/1365-2478.70079. hal-05484286

HAL Id: hal-05484286

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

Submitted on 2 Mar 2026

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

ORIGINAL ARTICLE

OPEN ACCESS

# Joint Inversion of Electromagnetic and Direct Current Resistivity Data Using Trust Regions. Application to Uranium Exploration in the Athabasca Basin

Ali Mohand-Said¹ | Guy Marquis¹ | Serge Sambolian¹ | Jean-François Girard¹ | Grant Harrison² | Elodie Williard³

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

¹ITES, Université de Strasbourg/CNRS/ENGEES, Strasbourg, France | ²Orano Canada Inc., Saskatoon, Canada | ³Orano Mining, Châtillon, France

**Correspondence:** Ali Mohand-Said (amohandsaid@unistra.fr)

**Received:** 24 December 2024 | **Accepted:** 28 August 2025

**Funding:** This study is supported by the GeoMin3D Industrial Chair (PI Dr. Julien Mercadier, Géoressources Nancy), funded by Orano and ANR.

## ABSTRACT

The Athabasca Basin (Saskatchewan, Canada) is a world-class uranium mining province hosting high-grade high-tonnage deposits. Electromagnetic data and direct current resistivity data are essential tools to detect deep geoelectric structures associated with mineralization. Both methods are sensitive to electrical resistivity but highlight different structures. On the one side, electromagnetic methods reveal deeply buried, highly conductive graphitic structures. On the other side, direct current resistivity methods reveal milder contrasts of resistivity at shallower depths. We are here exploring the benefits that can be expected from two-dimensional joint inversion of electromagnetic and direct current resistivity data for the exploration of unconformity-related uranium deposits of the Athabasca Basin. Our methodology is recovering a single resistivity model to fit both datasets. We used a trust-region globalization approach to regularize the local minimization sub-problems, thus avoiding the task of regularization parameter tuning. Several tests are first conducted on synthetic models. These tests show that stand-alone electromagnetic inversions are able to recover the position of conductive plates, but their geometry remains uncertain. On stand-alone direct current resistivity inversions, the layered background is recovered, as well as smeared anomalies of resistivity associated with graphitic conductors. Whenever a conductive halo overlies a conductive plate, the wide anomaly associated with the plate appears more conductive and slightly shallower but the conductive plate and the halo cannot be distinguished. In the presence of two closely spaced conductors, the conductive anomaly appears more conductive and wider, so that they cannot be distinguished. On stand-alone electromagnetic inversions, however, their separation is clear, but electromagnetic measurements are blind to alteration halos. Joint inversions give the most reliable models. Both the resistive background and the conductive plates are recovered. A better constrained background allows to recover more contrasted plates. Synthetic tests allowed us to confirm the potential to recover the footprint of a hectometric-scale conductor overlying a plate using joint inversion, where both stand-alone inversions failed. Following these synthetic tests, we present an application of our methodology to a dataset from the Waterbury–Cigar Lake area. Joint inversion allows to recover a geoelectric model reconciling both datasets. The model shows conductors better constrained below the depth of unconformity, allowing for interpretations of resistivity variations above them.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly
cited.

© 2025 The Author(s). *Geophysical Prospecting* published by John Wiley & Sons Ltd on behalf of European Association of Geoscientists & Engineers.
The Athabasca Basin (Saskatchewan, Canada) is a world-class uranium mining province hosting high-grade high-tonnage deposits. Uranium mineralization is found at the unconformity (UC), where Proterozoic basin sediments overlie Archean to Paleoproterozoic metamorphic basement rocks. Many exploration projects in the Athabasca Basin now target mineralization at depths of several hundred metres. Apart from the odd radioactive glacial boulders, there is no evidence at the surface of mineralization at depth (IAEA 2018). Thus, geophysical methods are first-order tools to define drilling targets (Powell et al. 2007).

Deposits are usually overlying deep-rooted, steeply-dipping, graphitic crustal structures. The current metallogenic model attributes ore formation to massive hydrothermal fluid flow through these structures and the UC surface, which constitute favourable hydrothermal fluid pathways. Hydrothermal events led to pervasive alteration of the sandstone layers hosting the mineralization. This alteration typically forms a halo encapsulating the deposits (Jefferson et al. 2007).

The association of the UC, graphitic structures and clay alteration in the sandstone is the main geophysical exploration target. Graphitic structures and alteration halos are both associated with electrical resistivity anomalies. Thus, localization of these objects is achieved using geophysical methods sensitive to electrical resistivity (Powell et al. 2007; Smith et al. 2011; Nimeck and Koch 2008). Graphitic structures are several orders of magnitude more conductive than basement gneissic and granitic rocks. Inductive methods are efficient in this context as the intensity of the eddy currents increases with conductivity. Airborne time-domain electromagnetic (TEM) methods are used in a first pass to identify conductive trends. Ground surveys then follow to focus on the most prospective targets (Nimeck and Koch 2008). Direct current (DC) resistivity surveys are also sensitive to graphitic conductors, but their main interest is to investigate resistivity variations in the sandstones, particularly those that are attributed to alteration halos (Powell et al. 2007; Legault et al. 2008). Indeed, alteration halos are usually not conductive enough to generate a strong secondary magnetic field while resistivity contrasts, even between resistive rocks, affect the flow of the current injected into the ground in a DC survey, and thus the electric potential at the surface.

Inversion is essential for the interpretation of geophysical data. When doing an inversion, we aim to obtain a model of a physical parameter (here, electrical resistivity), allowing us to reproduce a set of field observations, within a certain margin of error. In the Athabasca Basin, DC resistivity data inversion is routinely applied. Multiple case histories are reported in the literature (Powell et al. 2007; Legault et al. 2008; Long et al. 2017). However, interpretation of conductive structures on inversions is complex: the response of the graphitic conductors overprints that of less-contrasted objects above or at the UC. Published case studies of three-dimensional (3D) inversions of central-loop airborne TEM (Mir et al. 2022) and ground moving-loop TEM data (Keller 2019), recovered conductive bodies well constrained below the expected depth of unconformity. However, plate-modelling software such as Maxwell Plate (EMIT) and MultiLoop (Lamontagne Geophysics), allowing to approximate plate-like

conductor geometry by trial-and-error, remain an important tool for TEM data interpretation among practitioners, thanks to their small computational requirements.

Smith et al. (2011) tried to identify alteration markers, comparing DC resistivity and airborne electromagnetic (EM) inversion results in the context of exploration of UC-associated uranium deposits. Their results showed that airborne EM methods have very limited sensitivity to alteration in the presence of conductive overburden and lakes. However, DC resistivity inversion shows a low-resistivity zone connecting a shallow conductive structure, attributable to a lake or overburden, and a deep conductive structure, attributable to a graphitic conductor. With strong prior knowledge, that is drilling, it is reasonable to interpret it as alteration. However, without prior evidence of alteration, one can wonder if it is not a result of the smearing out of the conductors' footprints.

Moorkamp (2017) reviewed joint inversion involving EM measurements. In a joint inversion, instead of trying to fit one set of observations, many sets of observations are fitted simultaneously. In our case, TEM and DC resistivity measurements are related to the same petrophysical property, but have different sensitivities to it. Hence, the joint inversion results in a single model that will encompass the specific sensitivities of the different observations. Vozoff and Jupp (1975) first presented the outcome of a joint inversion of DC resistivity and magnetotelluric data to recover a resistivity model showing layer interfaces not seen on stand-alone inversions. Recent examples of joint inversions of data presenting different sensitivities to subsurface electrical resistivity showed more representative geoelectric models (Amatyakul et al. 2017; Christensen 2022; Ozyıldırım et al. 2020; Yi and Sasaki 2015).

The main thread of the research work we are presenting in this paper is the evaluation of the complementarity of TEM and DC resistivity methods, through joint inversion, in the context of uranium exploration in the Athabasca Basin. We are expecting, by inverting TEM and DC jointly, to alleviate the limitations identified in separate inversions and to get a more accurate geoelectric model. Modelling conductive anomalies associated with graphitic conductors could enhance the sensitivity of DC data to resistivity contrasts in sandstones and close to the UC, potentially related to fluid flow and recorded in the signal of DC data. And vice versa, shallow resistivity structures highlighted by DC measurements could lead to a refinement of the geometry of the conductive plates already identified with EM inversions. An assessment of the outcomes of joint inversion is made first by running two-dimensional (2D) stand-alone and joint inversions on synthetic datasets of ground moving-loop TEM and pole-dipole DC resistivity profiles on several canonical models. This first step allows us to evaluate our ability to recover resistivity structures by comparing the inversion results with the canonical models. The next step is to apply our joint inversion methodology to data acquired in the Waterbury–Cigar Lake area in the Athabasca Basin.

Along with the joint inversion study, we introduce an inversion strategy based on trust regions (Conn et al. 2000; Nocedal and Wright 2006, chap. 4). Non-linear inverse problems in geophysics are, in most of the cases, solved iteratively by minimizing local quadratic approximations of the cost function, made around suc-
cessive models closer and closer to an optimal model. Local sub-problems are often ill-conditioned, leading to solutions driving the inverse problem towards unrealistic solutions when solved without regularization. In its most common form, regularization is included in the inverse problem as a weighted functional added to the data misfit. Such an approach is referred to as Tikhonov regularization (Tikhonov and Arsenin 1977). Choosing the optimal weight for regularization can be a tedious task and may have a significant impact on the results of the inversion. Indeed, excessive model regularization can prevent recovering the structures needed to reproduce the data, while too little regularization can lead to an unrealistically rough recovered model, for marginal gains in data fitting. Trust-region methods are globalization algorithms that work iteratively by minimizing quadratic approximations of a cost function, say the L²-norm of the data discrepancy vector. Their specificity comes from the fact that the solutions of the local sub-problems should lie within trusted regions, defined by their radius size, where the quadratic models are considered to accurately represent the cost function. The constraints on the solution's size act as an implicit form of regularization, avoiding excessive perturbations where the model parameters' contribution to the data fitting is low. Exhaustive descriptions of the trust-region methods are given in Conn et al. (2000) and Nocedal and Wright (2006, chap. 4). Our implementation of the algorithm is developed further in this paper. Although trust-region methods are not common in geophysical inverse problems, cases are reported in the literature, mainly applied to seismic problems. In Adriaens et al. (2023) and He and Wang (2022), trust-region and line-search globalization strategies are compared on full-waveform inversion problems. To our knowledge, the use of trust-region strategies in EM inverse problems has only been reported in Malecki et al. (2020). Our aim here is not to compare trust-region methods with other algorithms in terms of computational cost or results obtained; therefore, all the results presented here were computed with the same algorithm.

# 2 | Modelling and Inversion Methodology

## 2.1 | Modelling

Maxwell's equations describe the relationship between electric and magnetic fields and constitute the basis for formulating resistivity problems in TEM (Nabighian and Macnae 1991) and DC (Binley and Kemna 2005). For our study, TEM and DC problems are discretized on tensor meshes and solved using the finite volume approach implemented in SimPEG open-source codes (Cockett et al. 2015; Heagy et al. 2017).

Modelling of electric potential for the DC resistivity problem is made through a 2.5-D approximation, considering a 2D resistivity model and a 3D source. For the TEM problem, the resistivity model was defined on a 2D grid and mapped on a 3D grid for computation of the B-field, considering uniform resistivity in the direction perpendicular to the model.

The equations defining TEM and DC problems are both dependent on electrical conductivity (reciprocal of resistivity). Thus, the inversion is run on a geoelectric model defined on a 2D grid common to both forward problems.

## 2.2 | Inversion Methodology

The inverse problem is addressed as a data-fitting problem minimizing the squared L²-norm of the difference between the observed $\mathbf{d}_{obs}$ and predicted data $\mathbf{d}_{pred}$, weighted by $\mathbf{W}_d$, a diagonal matrix containing the reciprocal of data uncertainties,

$$ \Phi = \frac{1}{2} ||\mathbf{W}_d(\mathbf{d}_{pred} - \mathbf{d}_{obs})||_2^2, \quad (1) $$

$$ \mathbf{d} = \begin{bmatrix} \mathbf{d}_{DC} \\ \mathbf{d}_{EM} \end{bmatrix}, \quad (2) $$

where the DC and EM subscripts describe the type of data. Data vectors are thus a simple concatenation of DC resistivity and TEM data. DC and EM data amplitudes span different orders of magnitude, and EM data amplitudes themselves can span several orders of magnitude. Considering this issue, we set the data uncertainties as percentages of the absolute values of the observed data to get an adimensional measure of the misfit, independent of the data type. An additional, constant noise floor was added to the uncertainties of EM data, thus increasing the relative noise of the lowest amplitude data. For the joint inversion case, we added factors weighting DC ($w_{dc}$) and EM ($w_{em}$) data contributions to the cost function, following Aquino et al. (2022). These weights are adjusted to avoid the dataset with the most data that dominates in the inversion. Formulae to compute the weights are given in (3a). $N_{em}$ and $N_{dc}$ correspond to the number of EM and DC data, respectively. $N_{data}$ corresponds to the total number of data. $A_{em}$ and $A_{dc}$ are the factors that are tuned, keeping their sum equal to 1, to balance dataset's relative contributions to the misfit.

$$ w_{em} = A_{em} \frac{N_{data}}{N_{em}}, \quad (3a) $$

$$ w_{dc} = A_{dc} \frac{N_{data}}{N_{dc}}. \quad (3b) $$

Thus, for the joint inversion case, we have:

$$ \mathbf{W}_d = \begin{bmatrix} \frac{w_{dc}}{\sigma_{DC}} & \mathbf{0} \\ \mathbf{0} & \frac{w_{em}}{\sigma_{EM}} \end{bmatrix}, \quad (4) $$

with $\sigma_{EM}$ and $\sigma_{DC}$ corresponding to the uncertainties associated with EM and DC data, respectively.

The non-linear inverse problem is addressed by solving a succession of linear sub-problems. Each sub-problem aims to minimize a local quadratic approximation of the cost function. The quadratic approximation, $\tilde{\Phi}_i$, at iteration $i$, is built through a second-order Taylor approximation of the cost function around a model $\mathbf{m}_i$:

$$ \begin{aligned} \Phi(\mathbf{m}_i + \Delta\mathbf{m}) &\approx \tilde{\Phi}_i(\Delta\mathbf{m}) \\ &\approx \Phi(\mathbf{m}_i) + \Delta\mathbf{m}^T \Phi'(\mathbf{m}_i) + \frac{\Delta\mathbf{m}^T \Phi''(\mathbf{m}_i) \Delta\mathbf{m}}{2} \end{aligned} \quad (5) $$
Minimizing $\tilde{\Phi}_i(\Delta\mathbf{m})$ consists of finding the model correction vector $\Delta\mathbf{m}$ for which its gradient with respect to $\Delta\mathbf{m}$, $\partial\tilde{\Phi}_i(\Delta\mathbf{m})/\partial\Delta\mathbf{m}$, is equal to 0. The following linear system of equations gives the solution to this sub-problem:

$$ \Phi''(\mathbf{m}_i)\Delta\mathbf{m}_i = -\Phi'(\mathbf{m}_i). \qquad (6) $$

The linear sub-problems thus involve the computation of first-order and second-order derivatives with respect to the model parameters. The derivatives are, respectively, the gradient $\Phi'(\mathbf{m}_i)$ and the Hessian $\Phi''(\mathbf{m}_i)$ of the cost function. Computation of these derivatives represents a significant computational effort. A Gauss–Newton or a quasi–Newton approximation of the Hessian is often used to reduce this effort. For large-scale problems, implicit computations of gradient and Hessian-vector products are used to limit storage of large, dense matrices. An alternative is to use limited-memory versions of quasi-Newton algorithms, storing only a few vectors to build an approximation of the Hessian or its inverse (Liu and Nocedal 1989). In our case, the inverse model is defined on a 2D grid with a relatively small number of parameters, allowing us to easily handle full sensitivity matrices. As SimPEG codes allow us to get implicit computation of sensitivity-vector products $\mathbf{J}^T\mathbf{v}$, using the parallelized mira-simpeg (https://pypi.org/project/Mira-SimPEG/) fork of the SimPEG package developed by Mira Geoscience, the full Jacobian matrix is efficiently constructed row by row in parallel. Each row corresponds to the Fréchet derivatives of recorded data at a receiver, with respect to the model parameters. For the joint inversion case, as the Fréchet derivatives of the DC and TEM problems correspond to the same resistivity model parameters, we built the Jacobian matrix upon a concatenation of separately computed Jacobian matrices of DC and TEM problems:

$$ \mathbf{J} = \begin{bmatrix} \mathbf{J}_{DC} \\ \mathbf{J}_{EM} \end{bmatrix}. \qquad (7) $$

From there, a Gauss–Newton approximation can straightforwardly be built with $\mathbf{J}(\mathbf{m}_i)$ being the Jacobian matrix of the current conductivity model $\mathbf{m}_i$:

$$ \Phi''(\mathbf{m}_i) \approx \mathbf{J}(\mathbf{m}_i)^T \mathbf{W}_d^T \mathbf{W}_d \mathbf{J}(\mathbf{m}_i). \qquad (8) $$

However, the Gauss–Newton system can be ill-conditioned, and some sort of regularization needs to be added to compute a consistent solution. Trust-region globalization (Conn et al. 2000; Nocedal and Wright 2006, chap. 4) methods allow us to handle ill-conditioning of the Hessian matrix by enforcing the solution to lie within a trust region defined by a radius $\Delta$. Meaning that the norm, further referred to as the size, of the model correction vector should be less than or equal to the radius of the trust region.

At each iteration $i$, we need to solve

$$ \min_{s.t. \quad ||\Delta\mathbf{m}||_2 \le \Delta_i} \quad \tilde{\Phi}_i(\Delta\mathbf{m}) = \Phi(\mathbf{m}_i) + \Delta\mathbf{m}^T \Phi'(\mathbf{m}_i) \\ + \frac{\Delta\mathbf{m}^T \mathbf{J}(\mathbf{m}_i)^T \mathbf{W}_d^T \mathbf{W}_d \mathbf{J}(\mathbf{m}_i) \Delta\mathbf{m}}{2}. \qquad (9) $$

Once a solution $\Delta\mathbf{m}_i$ is obtained, we evaluate the quality of the cost function approximation via the ratio of the effective misfit decrease and the expected misfit decrease based on the local quadratic approximation of the cost function. The effective misfit

decrease corresponds to the difference between the current misfit, $\Phi(\mathbf{m}_i)$, and the misfit obtained with the model perturbation $\Phi(\mathbf{m}_i + \Delta\mathbf{m}_i)$. The expected misfit decrease corresponds to the difference between the quadratic approximation of the cost function for a model perturbation being a null vector $\tilde{\Phi}_i(\mathbf{0})$ and for a model perturbation being $\Delta\mathbf{m}_i$, that is $\tilde{\Phi}_i(\Delta\mathbf{m}_i)$:

$$ r = \frac{\Phi(\mathbf{m}_i) - \Phi(\mathbf{m}_i + \Delta\mathbf{m}_i)}{\tilde{\Phi}_i(\mathbf{0}) - \tilde{\Phi}_i(\Delta\mathbf{m}_i)}. \qquad (10) $$

If the ratio $r$ is smaller than a given threshold, it means that the quadratic approximation of the cost function is not sufficiently accurate for the current radius $\Delta_i$. The step is then rejected, $\Delta_i$ is divided by 4, and a new step is evaluated. This operation is repeated in an inner loop until $r$ is greater than or equal to the threshold value, here chosen as 0.25. The inner loop iterations are numbered by the index $k$. The size of the solution vector $\Delta\mathbf{m}_i(\mu_k)$, measured by its norm, at each iteration, $k$ is effectively modulated by damping the approximate Hessian:

$$ \Delta\mathbf{m}_i(\mu_k) = -(\mathbf{J}(\mathbf{m}_i)^T \mathbf{W}_d^T \mathbf{W}_d \mathbf{J}(\mathbf{m}_i) + \mu_k \mathbf{I})^{-1} \Phi'(\mathbf{m}_i). \qquad (11) $$

The inner loop is initialized with $\mu_0$ equal to 0. If the sub-problem (11) cannot be solved, because it is ill-conditioned, with $\mu_0$, $\mu_k$ is estimated by minimizing an additional function $\Phi_\mu$:

$$ \Phi_\mu(\mu_k) = \frac{1}{||\Delta\mathbf{m}_i(\mu_k)||_2} - \frac{1}{\Delta_i}. \qquad (12) $$

The problem of minimizing (12) is solved using Newton's method in a sequence of iterations numbered by $l$. Starting with $\mu_k^0$ set as the absolute value of the lowest eigenvalue of the approximate Hessian plus the precision of the machine. Then $\mu_k^l$ is updated following (13) until a reasonably accurate solution is obtained, usually after two to three iterations. This process is described in Nocedal and Wright (2006, chap. 4) and Malecki et al. (2020).

$$ \mu_k^{l+1} = \mu_k^l - \frac{\Phi_\mu(\mu_k^l)}{\Phi_\mu'(\mu_k^l)} \\ = \mu_k^l + \frac{\frac{||\Delta\mathbf{m}_i(\mu_k^l)||_2^2}{\Delta\mathbf{m}_i(\mu_k^l)^T (\mathbf{J}(\mathbf{m}_i)^T \mathbf{W}_d^T \mathbf{W}_d \mathbf{J}(\mathbf{m}_i) + \mu_k^l \mathbf{I})^{-1} \Delta\mathbf{m}_i(\mu_k^l)}}{\frac{||\Delta\mathbf{m}_i(\mu_k^l)||_2 - \Delta_i}{\Delta_i}}. \qquad (13) $$

The less $\tilde{\Phi}_i$ is accurate, the smaller the radius is going to be, and the greater will be the damping factor $\mu_k$, leading to more regularization. Once a step is accepted, if $r$ is greater than a second given threshold, chosen as 0.75, the radius of the trust region of the next iteration $\Delta_{i+1}$ is set as the double of the current radius $\Delta_i$. Otherwise, the current radius $\Delta_i$ is set as $\Delta_{i+1}$.

Close to a minimum of $\Phi(\mathbf{m})$, the solution $\Delta\mathbf{m}$ from (6) becomes smaller, leading to a more accurate estimate of the cost function. Thus, $r$ will get closer and closer to 1, leading to successive increases of $\Delta$ and decreases of $\mu$, until $\Delta\mathbf{m}_i(\mathbf{0})$ lies within $\Delta_i$ and the Gauss–Newton step is accepted directly. As we get closer to a global solution, the trust-region constraint and the regularization both vanish. One of the advantages of trust-region methods is that as the radius evolves, so does the regularization. An illustration
of this process is presented in the Appendix. Instead of selecting a regularization parameter before starting an inversion, the user specifies an initial radius, for example a multiplier of the norm of the gradient, and the parameters for radius updates, that is thresholds of the ratio r, factors of the radius update and bounds of radius. Modification of these parameters could affect the convergence speed of the inversion, as well as its final result (Conn et al. 2000, chap. 17).

Using simple damping on the approximate Hessian, the solution of sub-problem (11) is stabilized by reducing the sensitivity to data noise. However, the effect of the inverse of the approximate Hessian on the gradient is mitigated, which leads to a loss of resolution where the model is poorly illuminated by the survey. Constraining the roughness of the solution instead of its size allows us to maintain a degree of sensitivity to poorly illuminated parameters of the model. This can be handled in trust-region methods by changing the inner product defining the norm that measures the size of the solution. The variability of results obtained with different inner products has been shown in Adriaens et al. (2023). The positive definite roughness operator **M**, introduced in Zuberi and Pratt (2018), that approximates the discrete Laplacian of model $\nabla^2$, is now used to define the trust region. Hence, the trust-region size controls the roughness of the model, $\beta$ corresponding to the smoothing length. In this framework, (9) becomes

$$ \min_{\text{s.t.} \ ||\Delta\mathbf{m}||_{M} \le \Delta i} \quad \Phi_i(\Delta\mathbf{m}) = \Phi(\mathbf{m}_i) + \Delta\mathbf{m}^T \Phi'(\mathbf{m}_i) \\ + \frac{\Delta\mathbf{m}^T \mathbf{J}(\mathbf{m}_i)^T \mathbf{W}_d^T \mathbf{W}_d \mathbf{J}(\mathbf{m}_i) \Delta\mathbf{m}}{2} \quad (14) $$

with

$$ \mathbf{M} = \mathbf{I} - \beta^2 \nabla^2. \quad (15) $$

The sub-problem then becomes equivalent to

$$ \Delta \mathbf{m}_i(\mu_k) = -(\mathbf{J}(\mathbf{m}_i)^T \mathbf{W}_d^T \mathbf{W}_d \mathbf{J}(\mathbf{m}_i)) + \mu_k \mathbf{M}^{-1} \Phi'(\mathbf{m}_i), \quad (16) $$

where the identity matrix is replaced by the positive definite operator **M**, giving a smoother solution and the size of the solution is evaluated with the M-norm. To solve (14), we used the iterative preconditioned Lanczos method described in Conn et al. (2000, chap. 7), using **M** as a preconditioner. As an extra constraint on the model parameters, we set limits on the logarithms of the conductivity values with the band-pass transform proposed by Key (2016). This condition prevents us from obtaining unrealistic conductivity values.

## 3 | Application to Synthetic Models

A set of four synthetic canonical models for the resistivity distribution around a UC-related deposit was built to assess the effectiveness of our joint inversion algorithm and the recoverability of the different objects. In all configurations, an identical layered background hosts different types of targets. The three layers are from the surface: a 50-m thick overburden of glacial tills, a 400-m thick resistive sandstone layer and a basement.

The scenarios considered are (A) a conductive plate in the basement (Figure 1A); (B) a conductive plate in the basement with an alteration halo overlying the plate with 200 m of vertical extent (Figure 1B); (C) a conductive plate with a flat lying alteration along the unconformity (Figure 1C); (D) two conductive plates and a vertical offset of the unconformity (Figure 1D).

Resistivity values of the basement (5000 Ω m), fresh sandstones (2000 Ω m) and altered sandstones (500 Ω m) are comparable to those observed from statistics on petrophysical measurements on rock samples made by Orano Canada Inc. We fixed the till resistivity at 1000 Ω m and the plate resistivity at 10 Ω m. The overall resistivity values range from 10–5000 Ω m, but, for ease of interpretation, in particular, in the range of values of the alteration, we narrowed the colour scale to 100–4000 Ω m for all the resistivity models shown in this section.

For each synthetic model, DC apparent resistivity data and the vertical component of the B-field were calculated on a grid finer than the one used for the inversions. We used the same survey parameters as the contractors in the Waterbury–Cigar Lake area: details are given in Section 4. Synthetic apparent resistivity data were contaminated with a 5% Gaussian noise. B-field measurements were contaminated with 2.5% Gaussian noise and a random white noise of 0.01 pT. We selected the noise levels based on those observed in the field data. Sample curves showing the added noise for three EM stations are shown in Figure 2. The absolute noise of low-amplitude data of station X = −700 m (Figure 2A) is lower than for high-amplitude data of station X = 100 m (Figure 2B), but for relative values (Figure 2D,E) the trend is reversed. Estimates of measurement errors during data processing can yield lower noise levels for the highest amplitude data or higher noise for the lowest amplitude data. Thus, for the inversion of field data, we conservatively chose to stop the inversion before the root-mean-square (RMS) reached 2.0. For the synthetic cases, of course, we do not have this issue. Profiles with added noise are shown in Figures 3 and 4.

Separate DC and TEM inversions, as well as a joint inversion, are run for all models. For each inversion, the discretization and algorithm remain the same. Starting models for the inversions are half-spaces of 1000 Ω m. Inversions are run for up to 59 iterations or until the radius of the trust region reaches a predefined lower limit, which was set at one-millionth of its initial value to prevent the inversions from being stuck in an infinite loop of radius reduction and to give enough margin in case the initial radius was set too high. In practice, this lower limit was, however, never reached. Results shown here correspond to the first iterations reaching an RMS threshold value of 1.1, except for the joint inversion of the two-plate model (9D) for which the limit was set at 1.2 to avoid overfitting of the EM data. Cumulated sensitivities, introduced in Christiansen and Auken (2012), give insight on the depth of investigation (DOI) and are calculated for all models. However, interpretation of the DOI based on cumulated sensitivities should be done carefully as the threshold value selection is user-dependent. For all the figures, we kept the same threshold of 0.001 to avoid biases between the different inversions. Thus, it might over- or underestimate the actual DOI. We use it mainly as an indication that the sensitivity values decrease significantly at some depth.
<figure><img src="image_5.png" /><figcaption>**FIGURE 1** | Synthetic resistivity models. (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 2** | Absolute noise added to the synthetic data of the single plate model (Figure 1A), normalized by the source current for three stations: (A) Station at -700 m, (B) station at 100 m and (C) station at 2500 m. Total percentage of noise added to the same synthetic data for (D) station at -700 m, (E) station at 100 m and (F) station at 2500 m.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**FIGURE 3** | Synthetic moving loop TEM noisy data: vertical component of magnetic flux normalized by source current. Each line corresponds to a single time channel. Red colour corresponding to the earliest time channel (0.1 ms) and blue colour corresponding to the latest time channel (2.06 ms). (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

## 3.1 | Stand-Alone Electromagnetic Inversions

The first set of inversions is done with synthetic EM data only. All four inversions reached an RMS value of 1.1. Profiles of calculated data are shown in Figure 5. On all models, conductive anomalies embedded within a resistive background are recovered

(Figure 6). The overburden layer is not recovered, but an increase in resistivity related to the unconformity is visible.

The three models obtained from data simulated on the single-plate models (i.e. A, B, C) are very similar (Figure 6A–C), showing a single conductive anomaly centred just below the
<figure><img src="image_8.png" /><figcaption>**FIGURE 4** | Synthetic pole-dipole DC resistivity noisy data. (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 5** | Calculated TEM data corresponding to final inverse models of synthetic TEM data inversions. Crosses correspond to observed data and plain lines to calculated data. Each line corresponds to a single time channel. Red colour corresponding to the earliest time channel (0.1 ms), and blue colour corresponding to the latest time channel (2.06 ms). (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

unconformity, at the expected position of the plate. Resistivity values of these anomalies are of the expected order of magnitude with a minimum of 10 Ω m. However, the conductive anomalies do not extend at depth. It is thus not possible to infer the dip or vertical extent of a plate from these inversion results.

Although an effect from the occurrence of an alteration body overlying the conductive plate can be seen on synthetic EM profiles (Figure 5), the models obtained by inversions of data from cases B and C (that include alteration bodies) are not significantly different from the one inverting case A (plate-only) data.

The model obtained from the inversion of data from model D (two plates) shows two conductive anomalies related to the top of the conductive plates. Their resistivity are also in the right range of resistivity, with minima at 10 Ω m. These two anomalies are merging at depth, which is misleading, as they are supposed to be parallel. Thus, as for the test with single-plate models, it is not possible, from this result, to make interpretations of plate dip and vertical extent.

## 3.2 | Stand-Alone Direct Current Resistivity Inversions

A set of four inversions was also run with the synthetic DC resistivity data only, and all reached the RMS threshold of 1.1. Pseudo-sections of apparent resistivity are shown in Figure 7. For all models (Figure 8), the resistive, layered structure of the background is recovered. A shallow conductive layer of 700 Ω m corresponding to the overburden is recovered, followed by a more resistive layer of sandstone and an even more resistive layer below the unconformity. On all four models, conductive anomalies (100–400 Ω m), matching with the original plate positions, are recovered. Widths and resistivity values vary between the different models. Comparing the single-plate models (Figure 8A–C), we see that the models with alterations (Figure 8B,C) show anomalies that are more conductive than the model without alteration (Figure 8A). However, there is no distinctive anomaly attributable to alteration. Although the recovered models differ, it could be assumed that it is related to conductive plates being more conductive. The model obtained from inversion of two-plate data
<figure><img src="image_10.png" /><figcaption>**FIGURE 6** | Inverse model of inversions of synthetic TEM data. White dashed lines correspond to cumulated sensitivity values of 0.001. These models correspond to the first models having an RMS lower than 1.1. (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**FIGURE 7** | Calculated DC resistivity data corresponding to final inverse models of synthetic DC resistivity data inversions. (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 8** | Inverse model of inversions of synthetic DC resistivity data. White dashed lines correspond to cumulated sensitivity values of 0.001. These models correspond to the first models having an RMS lower than 1.1, except for case C, which broke down at the tenth iteration. (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

(Figure 8D) shows a single conductive anomaly centred between the true positions of the two plates. Although the anomaly is wider and its conductivity is higher than on the three other models, it is not possible to assess whether there are several plates or only a single wide conductive structure.

### 3.3 | Joint Inversion

Results from stand-alone DC and EM inversions clearly show different structures. Shallow and resistive structures are well identified on DC inversions, while EM is almost blind to them
<figure><img src="image_13.png" /><figcaption>**FIGURE 9** | Inverse model of joint inversions of synthetic DC resistivity and TEM data. White dashed lines correspond to cumulated sensitivity values of 0.001. These models correspond to the first models having an RMS lower than 1.1, except for case D which stopped when the RMS reached 1.2 to avoid overfitting the EM dataset. (A) Single plate model; (B) Single plate with rectangle alteration model; (C) Single plate with flat alteration model; (D) Two plates model.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 10** | Calculated TEM data corresponding to final inverse models of joint inversions of synthetic DC resistivity and TEM data. Crosses correspond to observed data and plain lines to calculated data. Each line corresponds to a single time channel. Red colour corresponds to the earliest time channel (0.1 ms), and blue colour corresponds to the latest time channel (2.06 ms). (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

but shows significant increases in conductivity associated with conductive plates. Thus, we expect that a joint inversion with equivalent RMS of DC and EM measurements would produce models reflecting the sensitivities of the two methods.

Joint inversions were run on the four pairs of synthetic datasets. As stated in the inversion methodology section, an additional weighting was added to balance the contribution of EM and DC datasets (cf. Equation 3a). These weights were adjusted to ensure that the RMS of EM and DC datasets stay close when reaching the threshold value. For the single plate models (Figure 9A–C), the weights (3a) were kept constant, with $A_{em} = 0.2$ and $A_{dc} = 0.8$. But for the two-plate model (Figure 9D), we had to modify these weights to put more emphasis on the EM dataset $A_{em} = 0.35$ and $A_{dc} = 0.65$ in order to converge towards a solution close to the true model. Thus, to compensate for the higher contribution of the EM dataset and avoid any overfitting of it, the RMS threshold was set at 1.2. Looking at the pseudo-sections of apparent resistivity (Figure 11) and EM profiles (Figure 10), the models obtained fit equivalently DC and EM datasets. The model

with flat alteration (Figure 9C) along the basement shows an RMS of EM data (Figure 10C), suggesting overfitting but given that it is not significantly lower than 1, we kept it as the final model.

Looking at the inverse models (Figure 9), the layered structure of the background is recovered on all four models. Overburden, sandstones and basement are well differentiated. For the two-plate model (Figure 9D), the sandstone layer, although less resistive than the basement, is less well resolved than in the three other cases (Figure 9A–C). The conductive plates are also well identified on all models and seem to be better constrained spatially than on stand-alone inversion results (Figures 6 and 8).

The model obtained from flat alteration data appears slightly different from the model from single plate data (Figure 9C). A sub-horizontal anomaly of 800 Ω m lying above the unconformity between X = 100 m and X = 800 m matches with the contours of the true alteration body (Figure 1C) but does not recover it entirely and is not connected to the anomaly related to the plate. The identification of an alteration along the unconformity, associated
<figure><img src="image_15.png" /><figcaption>**FIGURE 11** | Calculated DC resistivity data corresponding to final inverse models of joint inversions of synthetic DC resistivity and TEM data. (A) Single plate model, (B) single plate with rectangle alteration model, (C) single plate with flat alteration model and (D) two plates model.</figcaption></figure>

with a conductive plate, does not seem to be straightforward given the results obtained. The model with the rectangle alteration (Figure 9B), however, is well resolved, and we can see the anomaly of the conductive plate spilling out into the contours of the original alteration halo. The resistivity inside the box of the alteration halo is higher (800 Ω m) than the expected value (500 Ω m), but is within the same order of magnitude. On the three other tests, the conductive plates are well constrained below the unconformity and their footprint does not spill out into the sandstones. Thus, it suggests that the occurrence of an alteration halo of relatively low resistivity could be imaged with joint inversion.

The model obtained from the inversion of two-plate data (Figure 9D) shows two distinct conductive bodies, both matching with the original plates' positions. The recovered anomalies seem misleadingly linked at depth. The sandstone layer is less homogeneous, and some misleading artefacts are seen below the depth of investigation defined by the threshold (0.01) of the cumulated sensitivity. However, these artefacts are not connected to the plates and could not be misinterpreted as an alteration halo lying above the plate as in the case of Figure 1B.

Joint inversion of DC and EM synthetic data showed significant improvement in the geoelectric models obtained. As expected, the models fit both datasets (Figures 10 and 11). Joint inversion led to a model with more contrasted conductive structures, well constrained below the unconformity. Assessment of the recoverability of alteration markers was of interest in these tests. Our results show that a blocky alteration halo of 500 Ω m overlying a more conductive body of 10 Ω m can be recovered. But a thin alteration halo extending along the unconformity could not be interpreted with our methodology. If multiple parallel plates, relatively close, are present, we can model their tops. But at greater depths, it seems that the plates are merging into a common root, which is not actually the case; thus, this interpretation of plate location should be confined to the very top part of the conductor.

**4 | Application to Data From the Athabasca Basin**

The area of interest (Figure 12) is in the vicinity of the Cigar Lake mine, in the eastern part of the Athabasca Basin (Saskatchewan,

Canada). The deposit is hosted in the sandstones and conglomerates of the Manitou Falls formation, directly overlying the UC with the basement pelitic gneisses, at a depth of around 450 m. The mineralization is expressed as elongated lenses of metre-scale height and width. The lenses spread over a cumulative distance of close to 2 km in the E-W direction. An alteration halo, extending up to 300 m above the mineralization, surrounds the deposit. The nature and intensity of the alteration vary with the distance to the deposit. Clay alteration of the sandstones is occurring close to the mineralization, then the sandstones are marked by quartz dissolution and finally bleaching. These alteration markers are associated with a decrease in resistivity. Silicification pods are locally identified in the bleached sandstones and are associated with an increase in resistivity. The mineralization is structurally associated with a basement graphitic conductor striking E-W and steeply dipping to the S. Several drillholes confirmed the extension of this conductor many kilometres to the east of the deposit. However, no mineralization was identified on the eastern part of the conductor (Bruneton 1993; Martz et al. 2019). The area of interest in the Waterbury–Cigar property is mainly covered by Waterbury Lake. Drumlins, eskers and muskegs are shaping the topography of the emerged lands. The Manitou Falls formations sediments are covered by an overburden of glacial sediments. The thickness of the overburden is variable, ranging from 4 to 76 m, averaging 30 m from drillhole data.

A TEM moving-loop survey was conducted in 2015. The system is composed of a Geonics EM-37 transmitter working with a single turn square loop of 400 m by 400 m. A bipolar square wave signal of 23.5 A is emitted in the loop at a base frequency of 5 Hz. An ElectroMagnetic Imaging Technology SMARTem24 receiver records the measurements of three components of secondary B-field from a Leibniz Institute of Photonic Technology/Supracon High Temperature Jessy Deep SQUID Transient Electromagnetic sensor over 29 geometrically spaced time gates ranging from 0.10 to 42.68 ms offtime. Several profiles oriented N-S, perpendicular to the main graphitic conductor and previously identified magnetic trends, were made in a Slingram configuration, that is the transmitter loop and the receiver move together along the profile between each measurement. Measurements are made every 100 m, with a constant offset of 800 m between the transmitter and the receiver, the receiver being located north of the transmitter.
<figure><img src="image_16.png" /><figcaption>**FIGURE 12** | Map of the area of interest. Lines highlighted in white correspond to the profiles used in the field data inversions.</figcaption></figure>

The next year, a DC resistivity survey was conducted along three profiles parallel to the aforementioned moving-loop profiles, two of them being coincident with moving-loop profiles. Data were acquired using a pole-dipole array with a dipole length of 100 m. The current electrodes (pole) are at the midpoints between potential electrodes (dipole). A maximum of 60 ‘n’ levels was recorded, ranging from 0.5 to 59.5.

For our work, we selected profile L9500E of each survey. Several drillholes have been made in the surroundings of these profiles, all along their length, so they are ground-truthed over a greater length. The drillholes document the depth of UC and the location of the conductive plates. Unfortunately, downhole resistivity measurements were not available there. The two profiles selected are coincident with almost all their common length. The northern part of the DC resistivity slightly deviates to the east. The moving-loop profile is longer, station positions covering 6.5 km, while the electrode positions of the DC resistivity profile cover 5.6 km. To reduce the computational cost of modelling, the moving-loop dataset was sampled down to 26 stations, with a spacing of 200

m from X = -2500 m to X = 2500 m. From the measurements recorded, only the vertical component of the 15 earliest time gates, ranging from 0.10 to 2.07 ms, was kept (Figure 13).

Given the lower signal-to-noise ratio when the receiver dipoles
are far away from the source, the DC resistivity survey measure-
ments were also down-sampled. Measurements with a geometric
factor greater than 100,000 were discarded. The pseudo-section of
apparent resistivity is shown in Figure 14.

4.1 | Inversions

The inversion parameters and noise levels used for the synthetic
case were kept for the real data inversions. As a reminder, 5% of
the DC data and 2.5% of the absolute value of the EM data plus a
noise floor of 0.01 pT. For the datasets weighting, we found that
the best balance is found for Aem = 0.3 and Adc = 0.7. Inversions
were run for 59 iterations. However, we decided to keep the first
models reaching an RMS value of 2.0 (Figure 15) to compensate
<figure><img src="image_17.png" /><figcaption>**FIGURE 13** | Moving Loop TEM profile L9500E: Vertical component of secondary B-field normalized by source current. Each line corresponds to a single time channel. Red colour corresponds to the earliest time channel (0.1 ms), and blue colour corresponds to the latest time channel (2.06 ms). The data shown correspond to the sample kept for the inversion.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**FIGURE 14** | Pseudo-section of apparent resistivity of profile L9500E. The data shown correspond to measurements having a geometric factor lower than 100,000.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**FIGURE 15** | RMS values through iterations for stand-alone DC resistivity (blue), TEM (red) and joint (green) inversions of Athabasca Basin field data. The dashed black line corresponds to the threshold RMS value of 2.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**FIGURE 16** | Inverse model of Athabasca basin field data. White dashed lines correspond to cumulated sensitivity values of 0.001. (A) Stand-alone DC resistivity inversion, (B) stand-alone TEM inversion and (C) DC resistivity and TEM joint inversion.</figcaption></figure>

for the uncertainty on the noise level and avoid overfitting to the noise as much as possible, although the data are of high quality. These inverse models are shown on Figure 16. DC and TEM predicted data calculated on for each inverse model are shown on Figures 17 and 18 respectively.

The stand-alone DC inversion shows a widespread conductive anomaly ranging from $X = -1000$ m to $X = 2700$ m laterally, with a top close to 400 m (Figure 16A). The core of this large body is located at 100 m, along the X axis, with a resistivity of 20 Ω m. This large anomaly underlies a resistive layer of few thousand Ω m, almost reaching the surface. A surface layer, of 50 m thickness, which corresponds to the minimum cell size, of few hundred to a thousand Ω m is modelled. Several lateral discontinuities cut this uppermost layer that could be attributed to lakes, land or topography variations. Synthetic moving-loop

TEM data were calculated for this model and are very different from the measured data, the RMS of EM data reaching 31.929 (Figure 18A). The low resolution of DC data at depth does not allow us to recover a conductor strong enough to generate the observed EM signals. Overall, the model is consistent with the drillhole data, the interfaces between the geological units being consistent with the layered structure recovered.

The stand-alone EM inversion result shows four distinct, nearly vertical conductive anomalies in a resistivity background of few thousand Ω m (Figure 16B). Two bodies centred at $X = -500$ m and $X = 0$ m of a few Ω m seem to be linked but are not showing a great continuation at depth nor a dipping trend. From the results obtained from the synthetic case (Figure 6D), it seems that the connection of the two plates is probably an artefact. Their top positions are located at a depth of 350 m. The following
<figure><img src="image_21.png" /><figcaption>**FIGURE 17** | Calculated DC resistivity data corresponding to final inverse models. (A) Stand-alone DC resistivity inversion, (B) stand-alone TEM inversion and (C) DC resistivity and TEM joint inversion.</figcaption></figure>

conductive body centred on X = 1200 m is slightly less conductive, with a resistivity of 20 Ω m but shows a clear, steeply dipping north direction. This body is also deeper than the previous ones, with a depth of 450–500 m. Another anomaly centred at X = 2400 m, not well constrained by the data as it is located at the end of the profile, can be seen. Its depth is well identifiable at 450 m with a low resistivity of 10 Ω m. The top of conductive plates recovered with stand-alone EM inversion is consistent with the expected depth of unconformity, and their conductivity contrasts sharply with the background. This result is consistent with the geological knowledge of the area. However, DC synthetic resistivity data calculated from this model show a large difference with observed data, with an RMS value of 40.419 (Figure 17B).

On the joint inversion results, the conductive plates identified on stand-alone EM inversions are recovered with slight changes (Figure 16C). The two central anomaly footprints are better constrained and do not extend towards the surface, except for the anomaly centred on X = 0 m. A low-conductivity zone of 300 Ω m extending from the top of the plate to 200 m of depth is

visible. This is echoed with the synthetic test on Figure 9B, which shows a low conductivity zone above a plate, unseen on both EM (Figure 6B) and DC (Figure 8B) stand-alone inversions. The plate centred on X = 1200 m is now more conductive, but its contrast is still not as sharp as for the two plates to its left. A low-conductivity zone of few hundred of Ω m is visible above the most conductive part of the anomaly, but it stays confined at a depth where the unconformity is expected. Thus, it would not be interpreted as an alteration marker. The right-most conductive body remains relatively unchanged.

To further validate the model obtained by joint inversion and evaluate the stability of our algorithm, inversions were run with the same data uncertainties and weighting but with starting models significantly less resistive (100 Ω m, Figure 19A) and more resistive (10,000 Ω m, Figure 19C). The results remain almost unchanged except for the deepest part of the model, below 1000 m, being more or less resistive. The sensitivities of DC and TEM data to that part of the model are, however, very low.
<figure><img src="image_22.png" /><figcaption>**FIGURE 18** | Calculated TEM data corresponding to final inverse models. Crosses correspond to observed data and plain lines to calculated data. Each line corresponds to a single time channel. Red colour corresponds to the earliest time channel (0.1 ms), and blue colour corresponds to the latest time channel (2.06 ms). (A) Stand-alone DC resistivity inversion, (B) stand-alone TEM inversion and (C) DC resistivity and TEM joint inversion.</figcaption></figure>

Other tests were run where we changed the initial radius of the trust region (Figure 20), selecting a smaller starting radius (Figure 20A). Again, the resulting model remains relatively unchanged, the main difference being that the inversion needed more iterations to converge, as it started slowly with very small steps, until the trust-region radius became large enough to bring significant model updates.

On all joint inversion results, a common conductive root at 1500 m depth seems to connect the different conductive plates. This feature is echoing the work of Tschirhart et al. (2022), showing an interconnection of UC-related uranium ore systems to common, crustal-scale, deep heat sources driving fluid flows through fault zones of the Athabasca Basin interpreted from magnetotelluric

and gravity inversions. However, at this depth, we are out of the limit marked by the cumulated sensitivity threshold, and the resistivity of this anomaly varies between the different runs. The reliability of this common conductive root on our inversion results, therefore, remains questionable. Inclusion of later time channels in the inversion could enhance the sensitivity to this part of the geoelectric model.

The geoelectric model obtained by joint inversion satisfies both EM and DC measurements; RMS values are, respectively, 1.970 (Figure 18C) and 1.902 (Figure 17C).

A closer look at the RMS of DC measurements (Figure 21) shows that the highest misfit values are attributed to data with the
<figure><img src="image_23.png" /><figcaption>**FIGURE 19** | Joint inversion models of Athabasca basin field data obtained from different starting models. White dashed lines correspond to cumulated sensitivity values of 0.001. (A) 100 Ω m starting model, (B) 1000 Ω m starting model (results are presented in Figure 16C) and (C) 10,000 Ω m starting model.</figcaption></figure>

shortest spacings, for both stand-alone DC inversion (Figure 21A) and joint inversion (Figure 21B), which likely results from a too coarse discretization at electrode positions. The inversion would highly benefit from a finer discretization of the shallowest tens of metres; however, this would radically increase the computation time. It can be noticed that the fit to the deepest data plotted around X = 0 m decreased with the joint inversion, compared to the stand-alone inversion. These data are associated with the two close conductive anomalies identified on EM and joint inversion results (Figure 16B,C).

RMS of EM data for each time channel are shown in Figure 22. The trends remain the same for the EM and joint inversion data. The mid-time channels are the closest to observed measurements. The amplitudes are slightly lower overall for joint inversion, as

the global EM RMS is also lower. It is worth noticing that joint inversion of EM and DC data seemed to affect the RMS of the earliest time channels, as the RMS of time channels from 0.1 to 0.191 ms are shows a greater RMS for joint inversion data than that of the EM inversion.

The same comparison about RMS values can be made station by station for EM data (Figure 23). The stations associated with the lowest RMS on the EM inversion model and the joint inversion model are not the same. Measurements from stations X = 300 m to X = 2100 m are, overall, better reproduced by the joint inversion model, while measurements from stations from X = -500 m to X = -100 m are better reproduced by the EM inversion model. These stations are those with the highest RMS on the joint inversion model and correspond to the location of important conductors.
<figure><img src="image_24.png" /><figcaption>**FIGURE 20** | Joint inversion models of Athabasca basin field data obtained from different starting trust-region radii. White dashed lines correspond to cumulated sensitivity values of 0.001. (A) Initial radius equal to 1, (B) initial radius equal to five times the norm of the gradient of the cost function (results are presented in Figure 16C) and (C) 10,000 Ω m starting model.</figcaption></figure>

Thus, it prevents excessive interpretation of conductivities at these locations.

Making the connection between the comparisons of misfits from stand-alone inversions and joint inversion, a common trend is emerging. The measurements related to areas of common sensitivity, around X = 0 m for DC data, early time channels and stations X = -500 m to X = -100 m for EM data, are better fitted on stand-alone inversions. The compromise to be found in these areas may slow down the inversion, requiring more iterations to achieve an equivalent misfit. Emphasising these data, either by discarding other data or by tuning their weight, would help to focus the inversion on areas of common sensitivity.

Just like for the synthetic cases, we see that DC and EM stand-alone inversions result in different resistivity models, each one lacking crucial information to explain both datasets. Joint inversion of these datasets gives a geoelectric model explaining both datasets and showing a reshaping of the conductive structures identified on the EM inversion. The shallowest features from 0 to 250 m depth identified on DC resistivity inversion remained unchanged. A relatively conductive object linked to the plate centred on X = 0 m shows similarities with the synthetic case with rectangle alteration. This suggests the potential presence of a conductive alteration zone in the sandstone above this conductive plate. Although the detailed study of the RMS shows that the data associated with structures that contribute to both EM and

DC data may be harder to recover, the stability of joint inversions with varying starting models and radius of trust region gives a good confidence in the structures modelled on the first 600 m of the model.

These results are consistent with the geology of UC-related uranium deposits of the Athabasca Basin, and bring significant insights on the subsurface geoelectric structures of the Waterbury–Cigar Lake area. However, results are dependent on the data uncertainties and the balance between fitting EM and DC datasets. In our case, the data uncertainties were set as a percentage of the absolute values of the data. A noise floor was also added to the EM data. This was not a major concern for the synthetic tests as it corresponds exactly to the noise added to the data. For the true data inversion, a proper evaluation of the data uncertainties may improve the reliability of the model. Indeed, the noise level might be higher at later times for the EM data as well as for DC resistivity data with the largest offsets between transmitter and receiver. In our study, we compensated for it by setting data uncertainties higher than expected and by stopping earlier the inversion at an RMS value of 2.

Another axis of improvement concerns the dimensions of the resistivity model. We made the assumption that the structures are two dimensional, but, in reality, the structures are three dimensional; they can be discontinuous or deviated laterally. Our results may benefit from upscaling of our methodology to three dimensions. At the same time, inclusions of all the components
<figure><img src="image_25.png" /><figcaption>**FIGURE 21** | DC data weighted residuals $W_d(d_{pred} - d_{obs})$. Each point corresponds to one measurement. (A) Stand-alone DC inversion and (B) joint inversion.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**FIGURE 22** | Joint inversion (green) and EM inversion (red) RMS value of TEM measurements for each time channel.</figcaption></figure>

of the B-field from the moving-loop survey would help to better constrain the 3D geometry of the conductive plates.

**5 | Conclusion**

This work aimed to assess the outcome of 2D joint inversion of DC resistivity and TEM data in the context of uranium exploration in the Athabasca Basin. We showed on synthetic tests that both

EM and DC stand-alone inversions, although they give insight on different parts of the resistivity structures, failed to recover geoelectric models reconciling both datasets. Models obtained by joint inversions are consistent with both EM and DC data. These models show the conductive plates identified on EM inversions and the layered resistive background identified on DC inversions. Inclusions of information from shallow structures led to improved structural orientations and contrasts of conductive plates. However, although the results are improved on joint
<figure><img src="image_27.png" /><figcaption>FIGURE 23 | Joint inversion (green) and EM inversion (red) RMS value of TEM measurements for each station.</figcaption></figure>

inversion in the case of multiple plates, compared to stand-alone
inversion, the geometry of the plates is not perfectly recovered
as the depth of investigation is limited to their shallowest part.
One of the objectives of the joint inversion of DC and EM data
was to highlight evidence of alteration markers that are invisible
on stand-alone inversions. Tests on occurrences of alterations
showed that a blocky alteration halo of hundreds of metres scale
overlying a conductor could be recovered with joint inversion.
We applied our methodology to real data from the eastern part
of the Athabasca Basin. Our results highlight four near-vertical
conductors at depths between 350 and 500 m. Information on
shallow parts of the model recorded in DC has a significant
effect on the footprint of the identified conductors. However,
more geological prior knowledge would be needed to strengthen
such an interpretation, but the conductive anomalies overlying
the conductor centred on X = 100 m suggest the occurrence
of alteration markers or conductive structure associated with a
graphitic conductor in the basement. Our results were obtained
using an algorithm of inversion based on a trust-region method.
Using it, we were able to bypass the cumbersome task of
tuning a regularization parameter. The stability of our inversions
with different starting parameters confirmed the effectiveness of
our methodology.

Acknowledgements

This research is supported by the GeoMin3D Industrial Chair (PI Dr. Julien Mercadier, Géoressources Nancy), funded by Orano and ANR. We thank Cameco and Orano for providing the Cigar Lake data. This work benefited from discussions with our colleague Dr. Alessandro Vinciguerra and from detailed reviews from Dr. Marion Jegen and an anonymous reviewer. The authors would also like to acknowledge the High Performance Computing Center of the University of Strasbourg for providing scientific support and access to computing resources. Part of the computing resources was funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data.

Data Availability Statement

Data are proprietary and are not shared.

References

Adriaens, X., L. Métivier, and C. Geuzaine. 2023. “Inner Product Preconditioned Trust-Region Methods for Frequency-Domain Full Waveform Inversion.” *Journal of Computational Physics* 493: 112469. https://doi.org/10.1016/j.jcp.2023.112469.

Amatyakul, P., C. Vachiratienchai, and W. Siripunvaraporn. 2017. “WSJointInv2D-MT-DCR: An Efficient Joint Two-Dimensional Magnetotelluric and Direct Current Resistivity Inversion.” *Computers and Geosciences* 102: 100–108. https://doi.org/10.1016/j.cageo.2017.02.010.

Aquino, M., G. Marquis, and J. Vergne. 2022. “Joint One-Dimensional Inversion of Magnetotelluric Data and Surface-Wave Dispersion Curves Using Correspondence Maps.” *Geophysical Prospecting* 70: 1455–1470. https://doi.org/10.1111/1365-2478.13239.

Binley, A., and A. Kemna. 2005. DC Resistivity and Induced Polarization Methods. In *Hydrogeophysics*, edited by Y. Rubin and S.S. Hubbard, 129–156. Water Science and Technology Library, vol. 50. Springer Netherlands. https://doi.org/10.1007/1-4020-3102-5_5.

Bruneton, P. 1993. "Geological Environment of the Cigar Lake Uranium Deposit." *Canadian Journal of Earth Sciences* 30: 653-673. https://doi.org/10.1139/e93-054.

Christensen, N. 2022. “Joint Inversion of Airborne TEM Data and Surface Geoelectrical Data. The Egebjerg Case.” *Journal of Applied Geophysics* 196: 104511. https://doi.org/10.1016/j.jappgeo.2021.104511.

Christiansen, A., and E. Auken. 2012. “A Global Measure for Depth of Investigation.” *Geophysics* 77: WB171–WB177. https://doi.org/10.1190/geo2011-0393.1.

Cockett, R., S. Kang, L. Heagy, A. Pidlisecky, and D. Oldenburg. 2015. “SimPEG: An Open Source Framework for Simulation and Gradient Based Parameter Estimation in Geophysical Applications.” *Computers and Geosciences* 85: 142–154. https://doi.org/10.1016/j.cageo.2015.09.015.

Conn, A., N. Gould, and P. Toint. 2000. *Trust Region Methods*. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898719857.

He, Q., and Y. Wang. 2022. “Trust-Region Versus Line Search Globalization Strategies for Inexact Newton Method and Application in Full Waveform Inversion.” *Journal of Applied Geophysics* 201: 104639. https://doi.org/10.1016/j.jappgeo.2022.104639.
Heagy, L., R. Cockett, S. Kang, G. Rosenkjaer, and D. Oldenburg. 2017. "A Framework for Simulation and Inversion in Electromagnetics." *Computers and Geosciences* 107: 1–19. https://doi.org/10.1016/j.cageo.2017.06.018.

IAEA. 2018. "Unconformity-Related Uranium Deposits." TECDOC Series, IAEA-TECDOC-1857. International Atomic Energy Agency.

Jefferson, C., D. Thomas, S. Gandhi, et al. 2007. "Unconformity-Associated Uranium Deposits of the Athabasca Basin, Saskatchewan and Alberta." In *Mineral Deposits of Canada: A Synthesis of Major Deposit-Types, District Metallogeny, the Evolution of Geological Provinces, and Exploration Methods*, 273–305. Mineral Deposits Division, Special Publication No. 5. Geological Association of Canada. https://doi.org/10.4095/223744.

Keller, C. 2019. "Advancements in 3D Time-Domain Electromagnetic Inversion in the Athabasca Basin, Canada: A Unique Case Study From the West McArthur Project." In *SEG Technical Program Expanded Abstracts 2019*, 2133–2137. Society of Exploration Geophysicists. https://doi.org/10.1190/segam2019-3215196.1.

Key, K. 2016. "MARE2DEM: A 2-D Inversion Code for Controlled-Source Electromagnetic and Magnetotelluric Data." *Geophysical Journal International* 207: 571–588. https://doi.org/10.1093/gji/ggw290.

Legault, J., D. Carriere, and L. Petrie. 2008. "Synthetic Model Testing and Distributed Acquisition DC Resistivity Results Over an Unconformity Uranium Target From the Athabasca Basin, Northern Saskatchewan." *The Leading Edge* 27: 46–51. https://doi.org/10.1190/1.2831679.

Liu, D., and J. Nocedal. 1989. "On the Limited Memory BFGS Method for Large Scale Optimization." *Mathematical Programming* 45: 503–528. https://doi.org/10.1007/BF01589116.

Long, S., R. Smith, and R. Hearst. 2017. "Ground Resistivity Method and DCIP2D Forward and Inversion Modelling to Identify Alteration at the Midwest Uranium Deposit, Northern Saskatchewan, Canada." *Exploration Geophysics* 48: 383–393. https://doi.org/10.1071/EG15059.

Malecki, S., R.-U. Börner, and K. Spitzer. 2020. "A New Method for Absolute Underground Positioning Based on Transient Electromagnetics." *Geophysical Journal International* 221: 87–96. https://doi.org/10.1093/gji/ggz540.

Martz, P., J. Mercadier, M. Cathelineau, et al. 2019. "Formation of U-Rich Mineralizing Fluids Through Basinal Brine Migration Within Basement-Hosted Shear Zones: A Large-Scale Study of the Fluid Chemistry Around the Unconformity-Related Cigar Lake U Deposit (Saskatchewan, Canada)." *Chemical Geology* 508: 116–143. https://doi.org/10.1016/j.chemgeo.2018.05.042.

Mir, R., P. Fullagar, M. Darijani, et al. 2022. "Forward Modeling and 3D Inversion of Electromagnetic Data Collected Over the McArthur River Uranium Deposit in the Athabasca Basin, Canada." *Geophysics* 87: B129–B143. https://doi.org/10.1190/geo2021-0199.1.

Moorkamp, M. 2017. "Integrating Electromagnetic Data With Other Geophysical Observations for Enhanced Imaging of the Earth: A Tutorial and Review." *Surveys in Geophysics* 38: 935–962. https://doi.org/10.1007/s10712-017-9413-7.

Nabighian, M. N., and J. C. Macnae. 1991. "Time-Domain Electromagnetic Prospecting Methods." In *Electromagnetic Methods in Applied Geophysics: Volume 2, Application, Parts A and B*, 427–520. Investigations in Geophysics, vol. 3. Society of Exploration Geophysicists. https://doi.org/10.1190/1.9781560802686.ch6.

Nimeck, G., and R. Koch. 2008. "A Progressive Geophysical Exploration Strategy at the Shea Creek Uranium Deposit." *The Leading Edge* 27: 52–63. https://doi.org/10.1190/1.2831680.

Nocedal, J., and S. Wright. 2006. *Numerical Optimization*. Springer Series in Operations Research and Financial Engineering. Springer New York. https://doi.org/10.1007/978-0-387-40065-5.

Ozyıldırım, O., I. Demircil, N. Gündoğdu, and M. Candansayar. 2020. "Two Dimensional Joint Inversion of Direct Current Resistivity and Radiomagnetotelluric Data Based on Unstructured Mesh." *Journal of Applied Geophysics* 172: 103885. https://doi.org/10.1016/j.jappgeo.2019.103885.

Powell, B., G. Wood, and L. Bzdel. 2007. "Advances in Geophysical Exploration for Uranium Deposits in the Athabasca Basin." In *Proceedings of Exploration 07: Fifth Decennial International Conference on Mineral Exploration*, edited by B. Milkereit, 771–790.

Smith, R., R. Koch, G. Hodges, and J. Lemieux. 2011. "A Comparison of Airborne Electromagnetic Data with Ground Resistivity Data Over the Midwest Deposit in the Athabasca Basin." *Near Surface Geophysics* 9: 319–330.

Tikhonov, A., and V. Arsenin. 1977. *Solutions of Ill-posed Problems*. Winston.

Tschirhart, V., E. Potter, J. Powell, E. Roots, and J. Craven. 2022. "Deep Geological Controls on Formation of the Highest-Grade Uranium Deposits in the World: Magnetotelluric Imaging of Unconformity-Related Systems From the Athabasca Basin, Canada." *Geophysical Research Letters* 49: e2022GL098208. https://doi.org/10.1029/2022GL098208.

Vozoff, K., and D. Jupp. 1975. "Joint Inversion of Geophysical Data." *Geophysical Journal of the Royal Astronomical Society* 42: 977–991. https://doi.org/10.1111/j.1365-246X.1975.tb06462.x.

Yi, M.-J., and Y. Sasaki. 2015. "2-D and 3-D joint inversion of loop-loop electromagnetic and electrical data for resistivity and magnetic susceptibility." *Geophysical Journal International* 203, no. 2: 1085–1095. https://doi.org/10.1093/gji/ggv354.

Zuberi, M., and R. Pratt. 2018. "Mitigating Nonlinearity in Full Waveform Inversion Using Scaled-Sobolev Pre-Conditioning." *Geophysical Journal International* 213: 706–725. https://doi.org/10.1093/gji/ggx549.

## APPENDIX

As demonstrated in Section 2, the step size is modulated by the $\mu$ parameter. Thus, the radius of the trust-region $\Delta$ and $\mu$ are closely related. Figure A.1A shows Rosenbrock's function, which represents our cost function $\Phi$ here. Its minimum is at coordinates (1, 1). Let us say we are at the position $m_i$ (-1, -1). Figure A.1B shows the quadratic approximation of our cost function around $m_i$: $\Phi_i(\Delta m)$. In Figure A.1C, the solution of the local sub-problem (cf. Equation 9) is shown for different radii. We can clearly see that as the radius decreases, $\mu$ increases. It is worth seeing that not only the size of the solution changes but also its direction, getting closer and closer to the gradient direction as the damping increases.

Figure A.2 shows the evolution of trust-region radius and $\mu$ parameter through iterations for the joint inversion case of the real dataset presented in Section 4. As the trust-region radius increased during the first two iterations, the regularization stayed inactive and the Gauss-Newton steps were accepted directly. At iteration 3, the calculation of a model perturbation allows us to decrease the misfit required for regularization. The trust-region radius got smaller, and the regularization became active. From there, the trust-region radius oscillated until iteration 21, with the regularization parameter decreasing overall but being inactive at some iterations. After iteration 21, the regularization became inactive and the Gauss-Newton steps were constantly accepted directly.
<figure><img src="image_28.png" /><figcaption>**FIGURE A.1** | An example of trust-region radius size and $\mu$ parameter around a vector $\mathbf{m}_i$ on Rosenbrock's function. (A) Contours of Rosenbrock's function. Blue colour corresponds to the minimal value; the global minimum is at (1, 1). (B) A quadratic model approximating Rosenbrock's function around $\mathbf{m}_i$ is built. (C) Minimum of the quadratic approximation is shown for three different trust-region radii.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**FIGURE A.2** | Plot of trust-region radius (red) and regularization parameter satisfying the step size constraint (blue) through iteration for the joint inversion of field data.</figcaption></figure>