**The E**ffi**cient 3D Gravity Focusing Density Inversion Based on Preconditioned JFNK Method under Undulating Terrain: A Case Study from Huayangchuan, Shaanxi Province, China**

#### **Qingfa Meng 1,2, Guoqing Ma 1,2, Taihan Wang 1,2,\* and Shengqing Xiong 3,\***


Received: 2 July 2020; Accepted: 21 August 2020; Published: 22 August 2020

**Abstract:** Since polymetallic ores show higher anomalies in gravity exploration methods, we usually obtain the position and range of ore bodies by density inversion of gravity data. The three-dimensional (3D) gravity focusing density inversion is a common interpretation method in mineral exploration, which can directly and quantitatively obtain the density distribution of subsurface targets. However, in actual cases, it is computation inefficient. We proposed the preconditioned Jacobian-free Newton-Krylov (JFNK) method to accomplish the focusing inversion. The JFNK method is an efficient algorithm in solving large sparse systems of nonlinear equations, and we further accelerate the inversion process by the preconditioned technique. In the actual area, the gravity anomalies are distributed on the naturally undulating surface. Nowadays, the gravity inversion under undulating terrain was mainly achieved by discretizing the ground into unstructured meshes, but it is complicated and time-consuming. To improve the practicality, we presented an equivalent-dimensional method that incorporates unstructured meshes with structured meshes in gravity inversion, and the horizontal size is determined by the gradient of observed gravity and terrain data. The small size meshes are adopted at the position where the terrain or gravity gradient is large. We used synthetic data with undulating-terrain to test our new method. The results indicated that the recovered model obtained by this method was similar to the inversion method of unstructured meshes, and the new method computes faster. We also applied the method to field data in Huayangchuan, Shaanxi Province. The survey area has complicated terrain conditions and contains multiple polymetallic ores. Based on the high-density characteristics of polymetallic ore bodies in the area, we calculate the field data into 3D density models of the subsurface by the preconditioned JFNK method and infer six polymetallic ores.

**Keywords:** preconditioned jacobian-free Newton-Krylov (JFNK) method; undulating terrain; gravity focusing density inversion; adaptive equivalent-dimension; polymetallic minerals; unstructured mesh and structured mesh

#### **1. Introduction**

The gravity exploration method is used to detect polymetallic ores because of its higher density feature, and they can use the density inversion method of gravity anomaly to obtain the approximate horizontal position of the ores. The most commonly gravity inversion method is three-dimensional gravity inversion, which can directly obtain the subsurface density distribution to interpret the

three-dimensional area of ores [1–3]. The high-density borders of general gravity inversion are not sharp enough to match the ores. Last and Kubic presented compact gravity inversion, which is the precursor of gravity focusing inversion [4–6]. The focusing inversion uses the minimum support functional to increase the high-density value with the iteration of solution, so the volume of high-density results of inversion is minimal. In this way, the focusing inversion has a higher resolution and better convergence speed than the standard smooth inversion. However, the gravity method is suitable for the interpretation of mass data over a large area. Therefore, many geophysics make an effort to improve the computational efficiency of gravity inversion [7,8].

The JFNK method is designed for solving the large implicit nonlinear equations, and it is the combinations of Newton's method and Krylov subspace methods. This method uses the differential of vector instead of the jacobian matrix, so it does not need to form and store the elements of the true jacobian matrix. The jacobian matrix consists of the first-order partial derivatives of objective equation, which represents the optimal linear approximation of the equation to a given solution, so memory is intensive in solving large sparse nonlinear equations. Therefore, the JFNK method has a much faster convergence and smaller computation cost. This method has been successfully applied in fluid computing and other fields [9,10]. The Krylov subspace method is a kind of common iterative algorithm, and the conjugate gradient (CG) is a representative method of the Krylov subspace method which is commonly used in gravity inversion. Vandecar and Snieder used the preconditioned conjugate gradient method for solving the inversion problem of large data [11].

The study on the gravity inversion method considering the undulating terrain can be divided into two categories. The first type is to reduce the curved surface into a horizontal plane to offset the effects of undulating terrain, e.g., the equivalent source method, finite element method, boundary element method, Taylor series method and iterative method [12–16]. Yao also used the data of field and vertical gradient to reduce the curved surface into a horizontal plane to obtain more accurate and stable results [17]. To improve the direct solution method, Liang used the gradient data and boundary element method to obtain the potential field of the horizontal plane directly from the curved surface data [18,19]. Unfortunately, this method assumes that there is no field source between the undulating terrain and the converted horizontal plane. The data obtained by this way are only an approximation of the actual data, which reduces the precision of inversion results. The second type is to mesh the subsurface that corresponds to the undulating terrain. Currently, the common method is unstructured mesh technology. This method can fit the surface by dividing the subsurface into multiple triangles or tetrahedrons in two-dimensional (2D) or 3D inversions. Zhong used tetrahedrons to fit the terrain to realize terrain correction [20]. Geophysicists obtained the mesh with undulating terrain by using the triangulation technology to be applied on the inversion of the electrical prospecting with the complicated terrain [21,22]. Zhang applied triangulation technology on forward modeling and inversion of gravity and magnetic profiles [23]. In the gravity and magnetic inversion field, the triangulation was mainly used to fit geological models to obtain more accurate forward results [24]. The unstructured mesh method is complex and computationally intensive, which is not suitable for the inversion of gravity with a large amount of data. However, the common structured mesh method cannot fit the undulating terrain. Therefore, it is necessary to combine the structured and unstructured mesh methods in the gravity and magnetic inversion field.

In this paper, we presented the preconditioned JFNK gravity focusing density inversion method to obtain the density feature of the ores with accuracy and efficiency. To solve the inversion meshing under the undulating terrain, we also presented an adaptive equivalent-dimension method, which processed the subsurface by combining triangulation unstructured and rectangular structure mesh. Considering that the density variation and the terrain fitting are related to the anomaly and the change rate of the terrain, the horizontal size of each grid unit depends on the gradient of anomaly or terrain. We used the models with the undulating terrain and compared common unstructured meshing methods to gauge the performance of the algorithm and test its accuracy and efficiency. Then, we applied this approach to interpret gravity data collected in the Huayangchuan polymetallic mining area in Shaanxi Province.

#### **2. Methodology of Preconditioned JFNK Gravity Focusing Inversion**

The gravity inversion involved discretizing the subsurface into a finite number of rectangular units, and the discrete map is shown in Figure 1. There are *n* observation points and *m* rectangular units, the forward expression of the gravity anomaly could be expressed as:

$$\mathbf{g}\_{n \times 1} = \mathbf{G}\_{n \times m} \rho\_{m \times 1} \tag{1}$$

where *gn*×<sup>1</sup> is the gravity anomaly, *Gn*×*<sup>m</sup>* is the contribution to the *n*th datum of a unit density in the *m*th cell called kernel function matrix, and ρ*m*×<sup>1</sup> is the density matrix of rectangular units. 1 1

**Figure 1.** Discretization of underground space.

The kernel function matrix is [25]:

$$\mathbb{G}(n,m) = -\gamma \sum\_{i=1}^{2} \sum\_{j=1}^{2} \sum\_{k=1}^{2} \mu\_{ijk} \left[ \mathbf{x}\_i \ln(\mathbf{y}\_j + r\_{ijk}) + y\_j \ln(\mathbf{x}\_i + r\_{ijk}) - z\_k \arctan(\frac{\mathbf{x}\_i \mathbf{y}\_j}{z\_k r\_{ijk}}) \right] \tag{2}$$
  $\mathbf{x}\_i = \mathsf{S} - \xi\_i; y\_j = y$   $\mathbf{\mathcal{H}}$   $\eta\_j; z\_k = z$   $\mathsf{\mathcal{L}}\_{ki}$  $r\_{ijk} = \sqrt{\mathbf{x}\_i^2 + y\_j^2 + z\_k^2}; \mu\_{ijk} = (-1)^i (-1)^j (-1)^k.$ 

γ <sup>111</sup> (, , ) <sup>222</sup> (, , ) In the expression, γ is the gravitational constant, (*x*, *y*, *z*) are the coordinates of the observation points, and (ξ1, η1, ζ1) and (ξ2, η2, ζ2) are the coordinates of the minimum and maximum corner points of the cells.

1 Performing gravity inversion requires solving Equation (1) for ρ*m*×1. Since the number of underground cells was greater than the data, the Equation (1) was underdetermined.

We solved this underdetermined equation in the following manner [26]:

$$\boldsymbol{\phi} \cdot \boldsymbol{\phi} = \boldsymbol{\phi}\_{\mathcal{S}} + \boldsymbol{\mu} \boldsymbol{\phi}\_{\mathcal{P}} = \left(\boldsymbol{g}\_{\text{obs}} - \boldsymbol{G} \boldsymbol{W}\_{z}^{-1} \boldsymbol{W}\_{z} \boldsymbol{\rho}\right)^{T} \left(\boldsymbol{g}\_{\text{obs}} - \boldsymbol{G} \boldsymbol{W}\_{z}^{-1} \boldsymbol{W} \boldsymbol{\rho}\right) + \boldsymbol{\mu} \boldsymbol{\rho}^{T} \boldsymbol{W}\_{\varepsilon}^{T} \boldsymbol{W}\_{z}^{T} \boldsymbol{W}\_{\varepsilon} \boldsymbol{\mathcal{W}}\_{\varepsilon} \boldsymbol{\rho} \to \min \tag{3}$$

$$\mathcal{W}\_{\ell} = \left(\rho^2 + e^2\right)^{-1/2}\prime,\tag{4}$$

φ*g* ρ where φ*<sup>g</sup>* is the square norm of the difference between the observed anomaly(*gobs*) and the calculated anomaly(*G*ρ), and it represents the fitting functional of the data. In the expression, φ<sup>ρ</sup> is the stabilizing functional that constrained the inversion results to the real conditions. *u* is a regularization parameter. The value of the regularization parameter is usually 10*<sup>n</sup>* , and the *n* is determined by a "trial and error" method. *W<sup>z</sup>* is the depth-weighting function which was used to counteract the inherent decay of the kernel function [27]. *Wz*ρ is the weighted solution variable, we replaced it with ρ*w*. The Equation (4) is

the minimum support function [4]. In general, the e is 0.1, which is a focusing parameter to determine the sharpness of the inversion results.

The derivative formula of Equation (3) is:

$$\frac{\partial \phi}{\partial \rho} = \left[\mathcal{W}\_{\varepsilon}^{-2} \left(\mathcal{G}\mathcal{W}\_{z}^{-1}\right)^{T} \left(\mathcal{G}\mathcal{W}\_{z}^{-1}\right) + \mu\mathbb{I}\right] \times \rho\_{w} - \mathcal{W}\_{\varepsilon}^{-2} \left(\mathcal{G}\mathcal{W}\_{z}^{-1}\right)^{T} \mathcal{g}\_{\text{obs}} \tag{5}$$

We can substitute the solution of Equation (3) with the minimization of Equation (5). Finally, we used the JFNK method to obtain the ρ*<sup>w</sup>* of Equation (5), and the results ρ can be obtained by removing the depth-weighting function.

The JFNK method was introduced as follows. The *f*(*x*) = 0 can be processed as:

$$f(\mathbf{x}) \approx f(\mathbf{x}\_k) + f'(\mathbf{x}\_k)(\mathbf{x} - \mathbf{x}\_k) = \mathbf{0},\tag{6}$$

$$f'(\mathbf{x}\_k)(\mathbf{x} - \mathbf{x}\_k) = -f(\mathbf{x}\_k), k = 0, 1, \dots \tag{7}$$

The Krylov subspace method was created to solve the Equation (7), and we used the CG method which is a common solution technique with the Krylov subspace [26]. The matrix-vector product required by CG method was approximated by forward finite difference schemes:

$$f'(\mathbf{x})v \approx \frac{f(\mathbf{x} + \varepsilon v) - f(\mathbf{x})}{\varepsilon}, \varepsilon \neq 0,\tag{8}$$

where *v* is a vector of the Krylov subspace. The process of the preconditioned CG method used to solve Equation (7) is:

$$\begin{aligned} A &= f'(\mathbf{x}); \quad m = (\mathbf{x} - \mathbf{x}\_k); \; b &= -f(\mathbf{x}),\\ k &= 0; \qquad m\_0 = \text{initial solution estimate}; \; r\_0 = A^T(b - Am\_0),\\ \text{while } &\qquad r\_k \neq 0\\ z\_k &= Sr\\ k &= k + 1\\ \text{if } k &= 1\\ p\_1 &= z\_0,\\ \text{else} \\ \beta\_k &= r\_{k-1}^T z\_{k-1} / r\_{k-2}^T z\_{k-2} \\ \beta\_k &= z\_{k-1} + \beta\_k p\_{k-1} \\ \text{end}\\ q\_k &= Ap\_k\\ \alpha\_k &= r\_{k-1}^T z\_{k-1} / q\_k^T q\_k\\ m\_k &= m\_{k-1} + \alpha\_k p\_k\\ r\_k &= r\_{k-1} - \alpha\_k A^T q\_k\\ \text{end} \end{aligned} \tag{9}$$

where *k* is the number of iterations, and as the number of iteration increased, the residual *r* was minimized by moving a distance α in search direction *p*. The *S* is a preconditioner of the search direction, and we could obtain the new search direction. The optimum preconditioning is (*A* <sup>T</sup>*A*) -1. It produced the least-squares solution to Equation (7) in a single iteration. Since we were unable to obtain the optimum preconditioning directly, we substituted the optimum preconditioning with a depth-weighting function [10]. All the other parameters were intermediate variables with no significance.

Under the condition of the natural surface, the structured mesh method was not applicable, and the unstructured mesh method w triangulation (shown in Figure 2a). The structured method is computationally intensive. We proposed the adaptive equivalent-dimensional mesh method which combines the unstructured and structure mesh to achieve gravity inversion with the undulating terrain. Only the surface layer (the yellow area in Figure 2b) was discretized by the triangular unstructured mesh method, and the other subsurface was discretized into structure rectangular mesh. The three blue cells in Figure 2b have the same height, and it is shown that the structure rectangular mesh was the *n* dimension in each horizontal position. Therefore, many judgment processes were avoided in the calculation, and the calculation efficiency was improved; the memory was not increased in the case of fitting the surface as far as possible. This method could efficiently and accurately obtain the 3D density models, so it was more applicable to the gravity inversion with a large amount of data.

**Figure 2.** Mesh with triangulation and equivalent dimension in undulating terrain. (**a**)Triangulation; (**b**) equivalent dimension.

As the adaptive equivalent-dimensional mesh method, we derived the equations to determine intervals (*dx*, *dy*) of the mesh according to the terrain and anomaly:

$$\mathrm{d}x\_{l} = \frac{\mathrm{X}}{\mathrm{2}\eta\_{\mathrm{x}}} + \frac{\frac{\mathrm{u}\_{\mathrm{y}}}{\mathrm{\underline{\mathrm{x}}}} \mathrm{THD}\_{\mathcal{S}}(i,j)}{\mathrm{max}[\sum\_{j=1}^{\mathrm{N}} \mathrm{THD}\_{\mathcal{S}}(:,j)]} - \frac{\frac{\mathrm{u}\_{\mathrm{y}}}{\mathrm{\underline{\mathrm{x}}}} \mathrm{THD}\_{\mathrm{h}}(i,j)}{\mathrm{max}[\sum\_{j=1}^{\mathrm{N}} \mathrm{THD}\_{\mathcal{S}}(:,j)]}}{\sum\_{i=1}^{\mathrm{N}} \mathrm{2} - \frac{\sum\_{j=1}^{\mathrm{N}} \mathrm{THD}\_{\mathcal{S}}(i,j)}{\mathrm{max}[\sum\_{j=1}^{\mathrm{N}} \mathrm{THD}\_{\mathcal{S}}(:,j)]} - \frac{\frac{\mathrm{u}\_{\mathrm{y}}}{\mathrm{\underline{\mathrm{x}}}} \mathrm{THD}\_{\mathcal{S}}(i,j)}{\mathrm{max}[\sum\_{j=1}^{\mathrm{N}} \mathrm{THD}\_{\mathcal{S}}(:,j)]}}{\mathrm{max}[\sum\_{j=1}^{\mathrm{N}} \mathrm{THD}\_{\mathcal{S}}(:,j)]}$$

$$dy\_{j} = \frac{Y}{2n\_{y}} + \frac{Y}{\frac{\max\left[\sum\_{i=1}^{n\_{\mathcal{X}}} THD\_{\mathcal{S}}(i; \cdot)\right]}{\max\left[\sum\_{i=1}^{n\_{\mathcal{X}}} THD\_{\mathcal{S}}(i; \cdot)\right]} - \frac{\sum\_{i=1}^{n\_{\mathcal{X}}} THD\_{h}(i; \cdot)}{\max\left[\sum\_{i=1}^{n\_{\mathcal{X}}} THD\_{h}(i; \cdot)\right]}}{\sum\_{j}^{n\_{\mathcal{X}}} \left[2 - \frac{\sum\_{i=1}^{n\_{\mathcal{X}}} THD\_{\mathcal{S}}(i; \cdot)}{\max\left[\sum\_{i=1}^{n\_{\mathcal{X}}} THD\_{h}(i; \cdot)\right]} - \frac{\sum\_{i=1}^{n\_{\mathcal{X}}} THD\_{h}(i; \cdot)}{\max\left[\sum\_{i=1}^{n\_{\mathcal{X}}} THD\_{h}(i; \cdot)\right]}\right]} \tag{11}$$

$$THD\_f = \sqrt{\left(\frac{\partial f(\mathbf{x}, y)}{\partial \mathbf{x}}\right)^2 + \left(\frac{\partial f(\mathbf{x}, y)}{\partial y}\right)^2} \,. \tag{12}$$

$$\mathbf{x} = \mathbf{x}\_1 \dots \mathbf{x}\_N \qquad \mathbf{x}\_1 \dots \mathbf{x}\_{N-1} \dots \mathbf{x}\_1 \qquad \mathbf{x}\_{N-1} \dots \mathbf{x}\_1 \qquad \dots \qquad \mathbf{x}\_N \dots \mathbf{x}\_1 \dots \mathbf{x}\_N$$

1 1 1 1 [2 ] max[ ( ,:)] max[ ( ,:)] In the expression, *X* and *Y* represent the total length of the measuring area along x and y directions. *THD<sup>g</sup>* and *THD<sup>h</sup>* represent the total horizontal derivative (THD) of gravity anomaly *g* and elevation of measure points *h*, respectively. Its equation is shown in Equation (12). In this way, the *dx* and *dy* of the cells would be 0.5~1.5 times the average length. Therefore, the intervals of mesh were close in the area where the gravity anomaly or terrain varied greatly.

#### **3. Gravity Model Tests**

We tested our method in the inversion consisting of typical models of 150 × 150 × 100 m with undulating terrain, and the results are shown in Figure 3. The buried depth of the center of models was 100 m from the surface, and the central coordinates of two models were (225, 250) m and (725, 250) m. As we can see in Figure 3b, the extreme value position of the gravity effect was asymmetric because of the asymmetry of the undulating terrain. Therefore, the x-coordinates of the two extreme value positions of the gravity anomaly were 200 and 800 m, respectively. There was deviation from the x-coordinates of 225 and 725 m of the two actual models. When the terrain fitting accuracy was consistent, the terrain grid size was positively correlated with the fineness of the subsurface mesh. Figure 3c shows the run time of programs in different terrain grid sizes with twenty-five obverse points. It is shown in Figure 3c that the run time of programs with equivalent-dimensional mesh was obvious shorter than with triangulation mesh.

2 2 (, ) (, ) ( ) ( )

**Figure 3.** Gravity inversion by the preconditioned JFNK method. (**a**) Surface topography; (**b**) gravity anomalies under undulating terrain; (**c**) comparison of mesh computing time; (**d**) slice at y = 250 m of inversion result with equivalent-dimensional meshes; (**e**) slice at y = 250 m of inversion result with triangulation mesh.

The 3D gravity inversion in Figure 3d was performed based on equivalent-dimensional mesh. Figure 3d shows the slice of the inversion result at y = 250 m, and the red curve is consistent with the undulating terrain and represents the top of the subsurface, and the black rectangles represent real positions of typical models. It can be seen in Figure 3d that the high-density units in the inversion result were similar to the actual range of models. Therefore, this method could obtain the information of models directly by the inversion with equivalent-dimensional mesh. Figure 3e shows the 3D gravity inversion with triangulation mesh, and the result of the inversion with equivalent-dimensional mesh was similar to the inversion with triangulation mesh, comparing Figure 3d,e. Therefore, the new meshing method proposed by us can obtain gravity inversion results consistent with the traditional method in the regions with undulating terrain efficiently.

Figure 4a shows the gravity anomaly of Figure 3a with noise, and the anomaly showed some distortion after the addition of 30 signal noise ratio noise. Figure 4b is the slice of the inversion result, and Figure 4c is the 3D view of slices. In Figure 4d, blue lines represent the edge of the models, and the yellow and blue blocks represent densities greater than 0.28 and 0.16 g/cm<sup>3</sup> in the inversion results,

**Figure 4.** Gravity inversion containing noise by the preconditioned JFNK method. (**a**) Synthetic gravity anomalies under undulating terrain containing noise; (**b**) slice at y = 250 m of inversion results; (**c**) slices show of inversion results; (**d**) 3D view of inversion results.

We processed a gravity anomaly with complex terrain to verify the method of equivalentdimensional mesh. There were two peaks and one valley—shown in Figure 5b. In the 100 m depth from the surface, there were three prisms of 100 × 100 × 100 m; the central coordinates of which were (250, 250), (250, 750), and (750, 500) m, respectively. Figure 5a is the gravity anomaly of models, and the gravity containing 30 SNR noise is shown in Figure 5c. Because of the complex terrain's effect, the coordinates of the extreme value points shown in Figure 5a were (250, 200), (250, 800), and (750, 550) m, and the coordinates were not consistent with the models' coordinates. We achieved the density inversion based on the mesh with the equivalent dimension and obtained the slice diagram of the inversion results shown in Figure 5e–h by cutting along the white dashed lines numbered 1 and 2 shown in Figure 5a–c. Figure 5d shows the high-density blocks of a density greater than 0.3 g/cm<sup>3</sup> , and the blue lines represent the edge of the models.

Although the extreme values' coordinates of the gravity anomaly deviated from the model positions at lines 1 and 2, the central position of the inversion results obtained by the gravity inversion based on equivalent-dimensional mesh were still consistent with the model position. The results shown in Figure 5d–h indicate that this method is applicable and stabilized in complex terrain.

**Figure 5.** Gravity inversion under complex terrain. (**a**) Gravity anomalies; (**b**) surface topography; (**c**) gravity anomalies containing noise; (**d**) 3D view of inversion results; (**e**) slice at x = 250 m of inversion results; (**f**) slice at x = 250 m of inversion results computed with noise; (**g**) slice at x = 750 m of inversion results; (**h**) slice at x = 750 m of inversion results computed with noise.

#### **4. Actual Data Processing**

In order to validate the applicability of the adaptive equivalent-dimensional mesh method in the field data, we processed the gravity data of polymetallic mining areas in Shaanxi province of China. The area is located in Huayangchuan in Shaanxi province. Its tectonic location belongs to Xiaoqinling intracontinental orogenic belt in the southern margin of north China landmass. The main outcrops are the Neoarchean erathem Taihua group (deep metamorphic crystalline complex), the Mesoproterozoic erathem Xionger group (volcanic–sedimentary), Gaoshanhe formation (clasolite), and the Luonan group (carbonatite); some areas have the Sinian series and the Cambrian system clasolite. Paleogene,

Neogene, and Quaternary are deposited in Cenozoic basins, and the Ordovician, Silurian, and Paleozoic to Mesozoic stratums are missing [28,29]. The geological condition of the survey area is shown in Figure 6a.

**Figure 6.** Geological and gravity of survey area. (**a**) The geological map; (**b**) Bouguer gravity anomaly.

Figure 6b shows the Bouguer gravity anomaly of the survey area. The center area obviously had a high-value gravity anomaly in Figure 6b, so there was a high-density body in the subsurface. There were polymetallic mines in the vicinity of this area, and the geophysical characteristics of the metal minerals were high values of density. Therefore, we inferred that the high-value gravity in this area was a polymetallic deposit. Metallic minerals were produced by magma intruding into the shallow strata of the subsurface through the faults, so the location of faults had great value for us to infer the range of minerals.

The total horizontal derivative, of which the maximum value indicates the location of the faults, is commonly used to detect faults in gravity interpretation [30]. Therefore, by comprehensively analyzing the geological map and the gravity data of this region, we finally divided the faults of this region into the total horizontal derivative and local gravity field map shown in Figure 7a.

**Figure 7.** Fracture interpretation results on total horizontal derivative and local gravity field; the dashed lines F1–F9 represent the faults interpreted by the total horizontal derivative. (**a**) Total horizontal derivative; (**b**) field of local gravity anomaly.

− The Bouguer gravity anomaly is the comprehensive response of the density at all depths. The 3D density inversion depth range was an altitude of −1500 m from the surface, so we needed to perform separation of the potential field by using the matched filter method. The expression of the matched filter is [31,32]:

$$\text{LnE}(\omega) = 2\text{LnB} - 2\text{H}\omega \tag{13}$$

( ) where *LnE*(ω) is the average logarithmic power of the potential field, *H* is the depth of the layer corresponding to the anomaly, ω is the wave number, and *B* is a constant related to the deep or shallow source anomaly.

− In this way, we removed the regional field anomaly corresponding to the density below −1500 m and extracted the shallow response of the gravity anomaly, namely the local field. The result is shown in Figure 7b, indicating that the local anomalies in this area were separated by faults; this is consistent with the geological characteristics of the ores distributed around the faults, which shows the accuracy of the matched filter method. The terrain in this region was complex with an average altitude of 1480 m and a maximum height difference of 509 m, as shown in Figure 8a. Considering the regional terrain relief, our adaptive equivalent-dimensional mesh method was suitable for this real case.

.

**Figure 8.** 3D inversion results of actual area. (**a**) Elevation of the terrain; (**b**) field of local gravity anomaly; (**c**) 3D view of inversion results greater than 0.12 g/cm<sup>3</sup> ; (**d**) slices of inversion results.

− The depth of −1500 m extending to surface of this region was meshed by the proposed method in this paper, and 3D gravity inversion was performed to obtain the density in the subsurface, as shown in Figure 8c,d. Figure 8b is the local field of gravity, in which the white dotted lines represent the horizontal position of the slices, and Figure 8d shows the slices of the 3D density model obtained from the gravity inversion results shown in Figure 8c. The density model in each section show the high-density units, which represent the position of the explained minerals.

Through the recognition of high-density units in the 3D density inversion model, the range of metallic ore bodies in this region was interpreted in Figure 9 by a red range. The ore bodies distributed near the faults that conform to the principle of the metal minerals were generated in the magma which intruded along the faults. There were six ore bodies explained in this area, and the number shown in Figure 9 around the red range is the height above sea level of these ore bodies.

**Figure 9.** The range of ore bodies interpreted from 3D inversion results.

#### **5. Conclusions**

We proposed the preconditioned JFNK method to perform the intensive inversion computation; therefore, the 3D focusing inversion algorithm became more efficient. In addition, we combined unstructured and structured meshes to achieve gravity inversion with undulating terrain. Therefore, the computational complexity is less than that of the unstructured mesh method in the forward operation of the kernel function. In order to maintain the efficiency, we improved the terrain fitting and fineness of meshes by the adaptive equivalent-dimensional method. The fineness of meshes is positively correlated with the gradient of gravity and terrain. Synthetic tests showed that our method was suitable for gravity inversion under undulating terrain, and the algorithm was fast with accurate inversion results. Finally, we applied it to the field data processing in Huayangchuan, Shaanxi Province, and the distribution range and depth of the ore bodies were inferred from recovered density models, which further verified the stability and practicability of the new method.

**Author Contributions:** Conceptualization, Q.M. and G.M.; data curation, T.W.; formal analysis, G.M.; funding acquisition, G.M.; investigation, Q.M.; methodology, G.M.; project administration, S.X.; resources, S.X.; software, Q.M.; supervision, S.X.; validation, T.W.; visualization, Q.M.; writing—original draft, Q.M.; writing—review and editing, G.M. and T.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the national key research and development plan issue (Grant No.2017YFC0602203), the Excellent Young Talents Fund project of Jilin Province (20190103011JH), National Science and Technology Major Project task (No. 2016ZX05027-002-003), the Ministry of Education, and the China Postdoctoral Science Foundation (2019M651209).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Recovering Magnetization of Rock Formations by Jointly Inverting Airborne Gravity Gradiometry and Total Magnetic Intensity Data**

**Michael Jorgensen 1,2 and Michael S. Zhdanov 1,2,\***


**Abstract:** Conventional 3D magnetic inversion methods are based on the assumption that there is no remanent magnetization, and the inversion is run for magnetic susceptibility only. This approach is well-suited to targeting mineralization; however, it ignores the situation where the direction of magnetization of the rock formations is different from the direction of the induced magnetic field. We present a novel method of recovering a spatial distribution of magnetization vector within the rock formation based on joint inversion of airborne gravity gradiometry (AGG) and total magnetic intensity (TMI) data for a shared earth model. Increasing the number of inversion parameters (the scalar components of magnetization vector) results in a higher degree of non-uniqueness of the inverse problem. This increase of non-uniqueness rate can be remedied by joint inversion based on (1) Gramian constraints or (2) joint focusing stabilizers. The Gramian constraints enforce shared earth structure through a correlation of the model gradients. The joint focusing stabilizers also enforce the structural similarity and are implemented using minimum support or minimum gradient support approaches. Both novel approaches are applied to the interpretation of the airborne data collected over the Thunderbird V-Ti-Fe deposit in Ontario, Canada. By combining the complementary AGG and TMI data, we generate jointly inverted shared earth models that provide a congruent image of the rock formations hosting the mineral deposit.

**Keywords:** inversion; gravity; magnetics; magnetization vector; remanent magnetization; three-dimensional; joint inversion

#### **1. Introduction**

In mineral exploration, magnetic data have traditionally been inverted to produce magnetic susceptibility models, which represent magnetization induced by the current magnetic field. This does not take into account, however, the remanent magnetization of the rocks produced by the ancient magnetic field. More information about rock formations and geological processes can be obtained by inverting magnetic data for magnetization vector, as opposed to magnetic susceptibility only. At the same time, increasing the number of inversion parameters (i.e., inverting for the scalar components of magnetization vector) increases the non-uniqueness of the inverse problem. We overcome this problem by jointly inverting multiple geophysical data sets. The total magnetic intensity (TMI) and airborne gravity gradiometry (AGG) data are typically gathered in coincident surveys, making them a natural choice for joint interpretation.

Over the last decade, significant efforts were made by many researchers to develop the methods of magnetic data inversion in the presence of remanent magnetization. To include both induced and remanent magnetization, one needs to model the distribution of magnetization vector within the rock formation rather than the distribution of scalar susceptibility alone. This approach was used by [1,2] to invert TMI data. This enables explicit inversion

**Citation:** Jorgensen, M.; Zhdanov, M.S. Recovering Magnetization of Rock Formations by Jointly Inverting Airborne Gravity Gradiometry and Total Magnetic Intensity Data. *Minerals* **2021**, *11*, 366. https:// doi.org/10.3390/min11040366

Academic Editors: Stanisław Mazur and Amin Beiranvand Pour

Received: 5 March 2021 Accepted: 29 March 2021 Published: 31 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of the magnetization direction and amplitude, rather than just the magnetization amplitude only (e.g., [3]). Ellis et al. [4] reported further progress in the solution of this problem; they introduced a technique for regularized inversion for the magnetization vector. From the magnetization vector, one can recover information about both the remanent and induced magnetization. A detailed overview of the published papers can be found in the comprehensive review paper by Y. Li [5], where one can find a long list of publications on the subject. We reference the interested readers to this publication for more information about the recent advances in 3D inversion of magnetic data in the presence of significant remanent magnetization

The results presented in the cited papers, however, illustrate the practical difficulties of the inversion related to the fact that, in this case, one has to determine from the observed scalar TMI data three unknown components of the magnetization vector within each cell instead of one unknown value of susceptibility. The main problem is with the increased practical non-uniqueness associated with the inverse problem.

In order to reduce the non-uniqueness of the inverse problem, we use a method of direct inversion of the magnetic data for the magnetization vector based on applying a Gramian stabilizer in the framework of the regularized inversion [6]. To this end, we also present a method of joint inversion of AGG and TMI data in the presence of remanent magnetization. We consider two different approaches to the joint inversion. The first approach is joint inversion with Gramian constraints [7], where the structural correlation of the different model gradients is enforced. The Gramian constraint is analogous to the cross-gradients method [8]; however, whereas the cross-gradients approach typically relies on an approximation, the Gramian constraint utilizes an exact analytical formula for the gradient direction of the parametric functional, which guarantees rapid convergence of the inverse solution [9].

The second approach presented is joint inversion with a joint focusing stabilizer [10,11], which enforces null-space sparsity in both models via a modified minimum support constraint. It is well known that the focusing stabilizers minimize the areas with anomalous physical properties (in the case of the minimum support stabilizer) or the areas where major changes in physical properties occurs (in the case of the minimum gradient support stabilizer). The joint focusing stabilizers force the anomalies of different physical properties to either overlap or experience a rapid change in the same areas, thus enforcing the structural correlation.

To illustrate both of these approaches, we present a case study of the joint inversion of the AGG and TMI data, which were collected over the Thunderbird V-Ti-Fe deposit in the Ring of Fire area of Ontario, Canada. The Ontario Geological Survey (OGS) and the Geological Survey of Canada (GSC) collaboratively gathered the airborne data with the Fugro Airborne Surveys gravity gradiometer and magnetic system between 2010 and 2011. Appropriate filters are applied to the data, i.e., reduction-to-pole (RTP) and deep regional trend removal. The filtered RTP data were used to invert the susceptibility model. The original TMI data with the removed regional trend were inverted for the magnetization vector. Standalone inverse models are produced using a single datum, i.e., AGG or TMI only, yielding standalone models to compare and contrast to the jointly inverted models and to produce appropriate model weights used in the joint inversions.

We compare and contrast the 3D density, magnetic susceptibility, and magnetization vector models of the Thunderbird V-Ti-Fe deposit obtained from standalone, Gramian, and joint-focused inversions. Application of both joint inversion methodologies to these data has resulted in 3D models of subsurface formations that have sharper geospatial boundaries and stronger structural correlations than the standalone inverted models. The Gramian constraint approach yielded the highest level of numerical structural correlation. The joint focusing approach provided similar results but was significantly faster and easier to implement.

#### **2. Forward Modeling of the Gravity and Magnetic Data**

#### *2.1. Gravity Forward Modeling*

The gravity field can be computed from the gradient of the gravity potential *U*:

$$\mathbf{g}(\mathbf{r}) = \nabla \mathcal{U}(\mathbf{r}) = \gamma \int \int \int\_{-D} \rho(\mathbf{r}') \frac{\mathbf{r}' - \mathbf{r}}{\left|\mathbf{r}' - \mathbf{r}\right|^3} d\mathbf{v},\tag{1}$$

where γ is the universal gravity constant; *r* and *r* ′ are the observation and integration points, respectively; ρ is the density; and the gravity potential, *U*, has the following form:

$$
\mathcal{U}I(\mathbf{r}) = \gamma \int \int \int\_{-D} \frac{\rho(\mathbf{r}')}{|\mathbf{r}' - \mathbf{r}|} dv. \tag{2}
$$

The second spatial derivatives of the gravity potential *U* form a symmetric gravity gradient tensor,

*g*ˆ = *gxx gxy gxz gyx gyy gyz gzx gzy gzz* , (3)

where

$$\mathcal{g}\_{\alpha\beta}(r) = \frac{\partial^2}{\partial \alpha \partial \beta} \mathcal{U}(r); \alpha, \beta = \propto, y, z. \tag{4}$$

The expressions for the gravity tensor components can be written as follows:

$$\mathbf{g}\_{\alpha\beta}(\mathbf{r}) = \gamma \int \int \int\_{D} \frac{\rho(\mathbf{r}')}{\left|\mathbf{r}' - \mathbf{r}\right|^{3}} \mathbf{K}\_{\alpha\beta}\left(\mathbf{r}' - \mathbf{r}\right) dv\_{\prime} \tag{5}$$

where kernels *K*αβ are equal to:

$$\mathbb{K}\_{\alpha\beta} \left( r' - r \right) = \begin{cases} \frac{3(\alpha - \alpha')(\beta - \beta')}{|r' - r|^2}, \alpha \neq \beta \\\ \frac{3(\alpha - \alpha')^2}{|r' - r|^2} - 1, \alpha = \beta \end{cases}; \alpha, \beta = \text{x, y, z}. \tag{6}$$

One can use the point-mass approximation to calculate Formula (5) in discretized form by considering each cell as a point mass [12].

#### *2.2. Magnetic Forward Modeling and Magnetization Vector*

The earth's current magnetic field is the vector sum of contributions from two main sources: the background magnetic field due to the dynamo in the earth's liquid core, and the crustal field due to magnetic minerals. Conventional magnetic inversion methods produce a 3D magnetic susceptibility model from the magnetic vector field, **H**, or from the total magnetic intensity (TMI) data, *T*. This is an effective tool for targeting magnetic mineralization; however, the following assumptions must generally be made: (1) there is no remanent magnetization, (2) self-demagnetization effects are negligible, and (3) the magnetic susceptibility is isotropic (e.g., [13–16]).

Under these assumptions, the intensity of magnetization, **I**(*r*), is linearly proportional to the inducing magnetic field, **H**0(*r*),

$$\mathbf{I}(\mathbf{r}) = \chi(\mathbf{r})\mathbf{H}\_0(\mathbf{r}) = \chi(\mathbf{r})H\_0\mathbf{1}(\mathbf{r}),\tag{7}$$

where χ(*r*) is the magnetic susceptibility and **l**(*r*) = *lx*, *ly*, *lz* is the unit vector along the direction of the inducing field. Given the inclination (*I*), declination (*D*), and azimuth (*A*) from the International Geomagnetic Reference Field (IGRF), the direction of the inducing

magnetic field can be computed as follows (assuming that the *x* axis is directed eastward, the *y* axis has a positive direction northward, and the *z* axis is downward):

$$\begin{array}{l} l\_{\boldsymbol{X}} = \cos(I)\sin(D - A) \\ l\_{\boldsymbol{Y}} = \cos(I)\cos(D - A) \\ l\_{\boldsymbol{z}} = \sin(I) . \end{array} \tag{8}$$

Thus, the anomalous magnetic field can be presented in the following form:

$$\mathbf{H}(\mathbf{r}) = -H\_0 \int \int \int\_{-D} \frac{\chi(\mathbf{r}')}{\left|\mathbf{r}' - \mathbf{r}\right|^3} \left[1 - \frac{\Im(\mathbf{1} \cdot (\mathbf{r}' - \mathbf{r}))(\mathbf{r}' - \mathbf{r})}{\left|\mathbf{r}' - \mathbf{r}\right|^2}\right] dv. \tag{9}$$

In airborne magnetic surveys, the total magnetic intensity (TMI) field is measured, which can be computed approximately as follows:

$$T(\mathbf{r}) \approx \mathbf{1} \cdot \mathbf{H}(\mathbf{r}) = -H\_0 \int \int \int\_{-D} \frac{\chi(\mathbf{r}')}{\left|\mathbf{r}' - \mathbf{r}\right|^3} \left[1 - \frac{\Im(\mathbf{1} \cdot (\mathbf{r}' - \mathbf{r}))^2}{\left|\mathbf{r}' - \mathbf{r}\right|^2}\right] d\mathbf{v}.\tag{10}$$

In the general case, however, the earth's magnetic field is nonstationary over geological time, meaning that the direction of magnetization in rock formations may differ from the direction of today's magnetic field. This is due to several factors, i.e., geomagnetic reversals and wandering poles. The ancient magnetization can be particularly useful when exploring magnetic structures such as kimberlites, dykes, iron-rich ultramafic pegmatitoids (IRUP), platinum group element (PGE) reefs, and banded iron formations (BIF). Moreover, remanence can also be used to determine the age of intrusive or alteration events.

By inverting for magnetization vector, one can take into account the effects of selfdemagnetization, anisotropy, and remanent magnetization. To include both induced and remanent magnetization, we need to model on the magnetization vector rather than the scalar susceptibility. This modifies Equation (7) as follows:

$$\mathbf{I}(\mathbf{r}) = H\_0 \mathbf{M}(\mathbf{r}),\tag{11}$$

where **M** has two parts: induced, **M***ind*, and remnant, **M***rem*, magnetizations:

$$\mathbf{M}(r) = \mathbf{M}\_{ind}(r) + \mathbf{M}\_{rem}(r). \tag{12}$$

Therefore, Equation (10) should also be modified:

$$T(r) \approx \mathbf{1}(r) \cdot \mathbf{H}(r) = $$

$$-H\_0 \mathbf{1}(r) \cdot \int \int \int\_V \frac{1}{\left| \mathbf{r'} - \mathbf{r} \right|^3} \left[ \mathbf{M} - \frac{\mathbf{3} (M\_k \cdot (\mathbf{r'} - \mathbf{r})) (\mathbf{r'} - \mathbf{r})}{\left| \mathbf{r'} - \mathbf{r} \right|^2} \right] dv\_r \tag{13}$$

One can compute the volume integral in Equation (13) in closed form for magnetic susceptibility. We can also evaluate the volume integral numerically with sufficient accuracy using single-point Gaussian integration with pulse basis functions provided the depth to the center of the cell exceeds twice the dimension of the cell [12].

#### **3. Inversion for Magnetization Vector Using Gramian Constraints**

In this section, we present a summary of the principles of the robust inversion for the magnetization vector using Gramian constraints.

It is well-known that the regularized solution of the geophysical inverse problem can be formulated as minimization of the Tikhonov parametric functional,

$$P^{\mathfrak{X}}(\mathfrak{m}) = \mathfrak{q}(\mathfrak{m}) + \alpha \mathbb{S}\_{\text{MN}, \text{MS}, \text{MGS}}(\mathfrak{m}) \to \min,\tag{14}$$

where *m* = *Mx*, *My*, *M<sup>z</sup>* is the magnetization vector and ϕ(*m*) is a misfit functional defined as the squared *L*<sup>2</sup> norm of the difference between the predicted, *A*(*m*), and observed, *d*, data,

$$\mathfrak{sp}(\mathfrak{m}) = ||A(\mathfrak{m}) - \mathfrak{d}||\_{L\_2}^2. \tag{15}$$

In the last formula, *A* is the forward modeling operator for the magnetic problem described by a discrete form of Equation (13).

Notations *SMN*,*MS*,*MGS* in Equation (14) indicate in compact form that we can use any one of three types of stabilizing functionals, *SMN*, *SMS*, or *SMGS*, based on minimum norm, minimum support, and minimum gradient support constraints, respectively. These stabilizers are determined as follows [17]:

$$\mathcal{S}\_{\rm MN}(\mathfrak{m}) = \left| \left| \mathfrak{m} - \mathfrak{m}\_{\rm apr} \right| \right|\_{L\_2}^2 = \int\_V \left( \mathfrak{m} - \mathfrak{m}\_{\rm apr} \right)^2 dv\_\prime \tag{16}$$

$$\mathcal{S}\_{\rm MS}(\mathfrak{m}) = \int\_{V} \frac{\left(\mathfrak{m} - \mathfrak{m}\_{apr}\right)^{2}}{\left(\mathfrak{m} - \mathfrak{m}\_{apr}\right)^{2} + \varepsilon^{2}} dv,\tag{17}$$

and

$$\mathcal{S}\_{\rm MGS}(m) = \int\_{V} \frac{\nabla m \cdot \nabla m}{\nabla m \cdot \nabla m + e^2} dv\_{\prime} \tag{18}$$

where *e* is a focusing parameter, which can be selected using an L-curve method [9].

The minimum norm stabilizer, *SMN*, minimizes the variations of the solutions, *m*, from the a priori model, *mapr*. The minimum support stabilizer, *SMS*, minimizes the volume occupied by the anomalous magnetization, while the minimum gradient support stabilizer, *SMGS*, selects the inverse models with sharp boundaries between the formations with different magnetic properties. Thus, by selecting the proper stabilizing functionals, the user may emphasize different properties of the inverse models. The minimum norm stabilizer usually results in relatively smooth distributions of magnetization, while the focusing stabilizers, *SMS* and *SMGS*, generate models with sharp boundaries. We use regularized inversion with focusing stabilization, as this recovers models with sharper boundaries and higher contrasts than the regularized inversion with smooth (e.g., minimum norm) stabilization.

The minimization problem (14) can be solved using a variety of optimization methods. For improved convergence and to avoid any matrix inversions, we minimize Equation (14) using the re-weighted regularized conjugate gradient (RRCG) method. We refer the reader to [17] for further details of this method.

Inverting for the magnetization vector is a more challenging problem than inverting for scalar magnetic susceptibility, because we have three unknown scalar components of the magnetization vector for every cell. At the same time, in the case of induced magnetization, there is an inherent correlation between the different components of the magnetization vector. In order to ensure a smooth change in the direction of the magnetization vector, we impose a similar condition by requiring that the different components of the magnetization vector should be mutually correlated as well. The results of numerical model experiments and case studies show that this approach to regularization ensures a robust inversion for magnetization vector.

It was demonstrated in [7] that one can enforce the correlation between the different model parameters by using the Gramian constraints. Following the cited paper, we have included the Gramian constraint in Equation (14) as follows:

$$P^{\mathfrak{\alpha}}(\mathfrak{m}) = \mathfrak{q}(\mathfrak{m}) + \lambda c\_1 \mathbf{S}\_{\mathbf{M}\mathbf{N},\mathbf{M}\mathbf{S},\mathbf{M}\mathbf{G}}(\mathfrak{m}) + \lambda c\_2 \sum\_{\mathfrak{\alpha} = \mathbf{x}, \mathbf{y}, \mathbf{z}} \mathbf{S}\_{\mathbf{G}} \left(\mathfrak{m}\_{\mathbf{x}\mathbf{y}} \mathbf{x}\_{\mathbf{e}ff}\right),\tag{19}$$

where *m* is the **3***N<sup>m</sup>* length vector of magnetization vector components; *m*<sup>α</sup> is the *N<sup>m</sup>* length vector of the α component of magnetization vector, α = *x*, *y*, *z*; χ*e f f* is the *Nm*-length vector of the effective magnetic susceptibility, defined as the magnitude of the magnetization vector,

$$\chi\_{eff} = \sqrt{M\_{\text{x}}^2 + M\_{\text{y}}^2 + M\_{\text{z}}^2}. \tag{20}$$

Functional *S<sup>G</sup>* is the Gramian constraint,

$$\mathbf{S}\_{G}\left(\mathfrak{m}\_{\alpha}\mathbf{x}\_{\mathcal{eff}}\right) = \left| \begin{pmatrix} (\mathfrak{m}\_{\alpha}, \mathfrak{m}\_{\alpha}) & \begin{pmatrix} \mathfrak{m}\_{\alpha}\,\chi\_{\mathcal{eff}} \\ \chi\_{\mathcal{eff}}, \mathfrak{m}\_{\alpha} \end{pmatrix} \\ \begin{pmatrix} \mathfrak{x}\_{\mathcal{eff}}, \mathfrak{m}\_{\alpha} \end{pmatrix} & \begin{pmatrix} \mathfrak{x}\_{\mathcal{eff}}, \chi\_{\mathcal{eff}} \end{pmatrix} \end{pmatrix} \right|, \tag{21}$$

where (∗, ∗)*L*<sup>2</sup> denotes the *L*<sup>2</sup> inner product operation [9].

Using the Gramian constraint (21), we enhance a direct correlation between the scalar components of the magnetization vector with χ*e f f* , which is computed at the previous iteration of an inversion and is updated on every iteration. The advantage of using the Gramian constraint in the form of Equation (21) is that it does not require any a priori information about the magnetization vector (e.g., direction and the relationship between different components) since the amplitude, χ*e f f* is computed at the previous iteration. On the first iteration, the scalar components are determined independently

The minimization problem (19) is solved using the re-weighted regularized conjugate gradient (RRCG) method. Details of the RRCG method and conjugate gradient derivations for the parametric functional (19) can be found in [9,17].

#### **4. Joint Inversion of Gravity and Magnetic Data**

#### *4.1. Gramian Joint Inversion*

In a case of joint inversion of gravity and magnetic data, the geophysical inverse problem is given by the following set of operator Equations:

$$\mathbf{d}^{(i)} = \mathbf{A}^{(i)} \left( \mathfrak{m}^{(i)} \right), (i = 1, 2), \tag{22}$$

where *m*(1) = *ρ* is the density model, *m*(2) = *Mx*, *My*, *Mz* is the magnetization vector model, *A* (*i*) are the forward modeling operators, *d* (*i*) are the data, and the superscript (*i* = 1, 2) indicates the gravity and magnetic problems, respectively. Operators *A* (*i*) are usually represented as discrete forms of the corresponding integral expressions (5) and (13) for gravity gradient and TMI fields, respectively.

We should note that the units and scales of the density and magnetization vector are very different. Therefore, in practical applications, we should operate with the normalized, dimensionless model parameters [9]:

$$
\widehat{\mathfrak{m}^{(i)}} = \mathfrak{W}\_{\mathfrak{m}}^{(i)} \mathfrak{m}^{(i)}, i = 1, 2; \tag{23}
$$

where *W* (*i*) *<sup>m</sup>* , *i* = 1, 2, are the corresponding linear operators of the model weighting. A trivial example of these operators could be just division of the model parameters by their upper or lower bounds, *m* (*i*) *max* or *m* (*i*) *min*, or by the corresponding interval of the model parameters distribution, *m* (*i*) *max* − *m* (*i*) *min* . This simple normalization makes the normalized parameters vary on the same scale. In a general case, the model weights can be defined based on integrated sensitivities, which will be discussed in details below.

Similarly, different data sets, as a rule, have different physical dimensions as well. Therefore, it is convenient to consider dimensionless weighted data, *d*g(*i*) , defined as follows:

$$
\overline{\mathbf{d}^{(i)}} = \mathbf{W}\_d^{(i)} \mathbf{d}^{(i)},\tag{24}
$$

where *W* (*i*) *d* , *i* = 1, 2, are the corresponding linear operators of the data weighting.

The solutions of geophysical inverse problems (22) are typically ill-posed [9]. A corresponding well-posed problem can be constructed by applying regularization and minimizing a parametric functional using the conjugate gradient method. In this case, misfit and stabilizing terms, each corresponding to the gravity gradiometry and magnetic data, are incorporated into a joint parametric functional, which is subject to the Gramian structural constraint:

$$P\left(\widehat{\mathfrak{m}^{(1)}},\widehat{\mathfrak{m}^{(2)}}\right) = \sum\_{i=1}^{2} \varphi\left(\widehat{\mathfrak{m}^{(i)}}\right) + \alpha \sum\_{i=1}^{2} \mathcal{S}\_{\text{MN,MS,MGS}}\left(\widehat{\mathfrak{m}^{(i)}}\right) + \beta G\left(\nabla \widehat{\mathfrak{m}^{(1)}}, \nabla \widehat{\mathfrak{m}^{(2)}}\right), \tag{25}$$

where *α* and *β* are the regularization parameters, responsible for degrees of smoothing/focusing (*α*) of the models and enforcing similarities (*β*) between the density and magnetization distributions, produced by a joint inversion. We use an adaptive algorithm of selecting these parameters, described in [9].

The misfit terms are defined as follows,

$$\mathfrak{sp}\left(\widehat{\mathfrak{m}^{(\overline{i})}}\right) = \left| \left| \widehat{\overline{A^{(\overline{i})}}} \left( \widehat{\overline{m^{(\overline{i})}}} \right) - \widehat{\overline{d^{(\overline{i})}}} \right| \right|\_{L\_2}^2 i = 1, 2; \tag{26}$$

where *A*g(*i*) *m*g(*i*) are the weighted predicted data:

$$
\widehat{A^{(i)}}\left(\widehat{m^{(i)}}\right) = \mathcal{W}\_d^{(i)} A^{(i)}\left(\widehat{m^{(i)}}\right). \tag{27}
$$

The stabilizing terms, *SMN*, *SMS*, *SMGS*, are based on minimum norm, minimum support, and minimum gradient support constraints, respectively, as defined by Formulas (16)–(18).

The Gramian term is defined by the following formula,

$$G\left(\widehat{\nabla \mathfrak{m}^{(1)}}, \widehat{\nabla \mathfrak{m}^{(2)}}\right) = G\left(\widehat{\nabla \mathfrak{m}^{(1)}}, \widehat{\nabla \widetilde{\mathbf{M}\_{\gamma}}}\right) = \sum\_{\gamma = x, y, z} \left| \begin{array}{c} \left(\widehat{\nabla \mathfrak{m}^{(1)}}, \widehat{\nabla \mathfrak{m}^{(1)}}\right) \\ \left(\widehat{\nabla \mathbf{M}\_{\gamma}}, \widehat{\nabla \mathfrak{m}^{(1)}}\right) \end{array} \begin{array}{c} \left(\widehat{\nabla \mathfrak{m}^{(1)}}, \nabla \widehat{\mathbf{M}\_{\gamma}}\right) \\ \left(\widehat{\nabla \mathbf{M}\_{\gamma}}, \nabla \widehat{\mathbf{M}\_{\gamma}}\right) \end{array} \right|,\tag{28}$$

where ∇*m*(1) is the gradient of the density model, ∇*M*g*<sup>γ</sup>* are the gradients of the scalar components of magnetization vector, and (∗, ∗)*L*<sup>2</sup> denotes the *L*<sup>2</sup> inner product operation [9]. Minimization of the Gramian aligns the model gradients, which in turn enforces the structural similarity of the shared earth model.

In a case of structural Gramian constraints defined by Formula (28), one can normalize the gradient vectors by their length:

$$
\widetilde{\nabla \mathfrak{m}^{(1)}} = \frac{\nabla \mathfrak{m}^{(1)}}{|\nabla \mathfrak{m}^{(1)}|}, \widetilde{\nabla \mathbf{M}\_{\Upsilon}} = \frac{\nabla \mathbf{M}\_{\Upsilon}}{|\nabla \mathbf{M}\_{\Upsilon}|}, \gamma = \mathbf{x}, y, z; \tag{29}
$$

with expression for Gramian (28) taking the form:

$$G\left(\widetilde{\nabla \mathfrak{m}^{(1)}}, \widetilde{\nabla \mathbf{M}\_{Y}}\right) = \sum\_{\gamma = \mathbf{x}, \mathbf{y}, \mathbf{z}} \left| \begin{pmatrix} \frac{\nabla \mathfrak{m}^{(1)}}{|\nabla \mathfrak{m}^{(1)}|} \,' \frac{\nabla \mathfrak{m}^{(1)}}{|\nabla \mathfrak{m}^{(1)}|} \\ \frac{\nabla \mathbf{M}\_{\gamma}}{|\nabla \mathbf{M}\_{\gamma}|} \,' \frac{\nabla \mathfrak{m}^{(1)}}{|\nabla \mathfrak{m}^{(1)}|} \end{pmatrix} \begin{pmatrix} \frac{\nabla \mathfrak{m}^{(1)}}{|\nabla \mathfrak{m}^{(1)}|} \,' \frac{\nabla \mathbf{M}\_{\gamma}}{|\nabla \mathbf{M}\_{\gamma}|} \\ \left| \frac{\nabla \mathbf{M}\_{\gamma}}{|\nabla \mathbf{M}\_{\gamma}|} \,' \frac{\nabla \mathbf{M}\_{\gamma}}{|\nabla \mathbf{M}\_{\gamma}|} \right| \end{pmatrix} \right|. \tag{30}$$

This normalization enforces the correlations of the unit vectors in the gradient directions, thus resulting in the structural similarity between the models representing density and magnetizations of the rock formations.

#### *4.2. Joint Focusing Inversion*

Another approach to joint inversion can be based on a joint total variation functional [18], or on joint focusing stabilizers, e.g., minimum support and minimum gradient support constraints [10,11]. It is well known that the focusing stabilizers minimize the areas with anomalous physical properties (in the case of the minimum support stabilizer) or the areas where major changes in physical properties occur (in the case of the minimum gradient support stabilizer). The joint focusing stabilizers force the anomalies of different physical properties to either overlap or experience a rapid change in the same areas, thus enforcing the structural correlation. We will demonstrate these properties in the model study presented in our paper.

The AGG and TMI data misfit terms are incorporated in the joint parametric functional and subject to the joint focusing stabilizers:

$$P\left(\widehat{\mathfrak{m}^{(1)}}, \widehat{\mathfrak{m}^{(2)}}\right) = \sum\_{i=1}^{2} \varphi\left(\widehat{\mathfrak{m}^{(i)}}\right) + \mathfrak{a}S\_{IMS,IMGS}\left(\widehat{\mathfrak{m}^{(1)}}, \widehat{\mathfrak{m}^{(2)}}\right). \tag{31}$$

The data misfit functionals are defined as above in Equation (26). The terms *SJMS*, and *SJMGS* are the joint stabilizing functionals, based on minimum support and minimum gradient support constraints, respectively [10,11].

For example, a joint minimum support stabilizer can be introduced as follows:

$$S\_{JMS} = \int\_{V} \frac{\sum\_{i=1}^{2} \left( \overbrace{\mathfrak{m}^{(i)}}^{-} - \overbrace{\mathfrak{m}\_{apr}^{(i)}}^{-} \right)^{2}}{\sum\_{i=1}^{2} \left( \overbrace{\mathfrak{m}^{(i)}}^{-} - \overbrace{\mathfrak{m}\_{apr}^{(i)}}^{-} \right)^{2} + \varepsilon^{2}} dv\_{\prime} \tag{32}$$

where *e* is the focusing parameter, which can be selected using an L-curve method [19]. In a similar way, we can introduce a joint minimum gradient support functional (JMGS):

$$S\_{IMGS} = \int\_{V} \frac{\sum\_{i=1}^{2} \left( \nabla \widehat{\mathfrak{m}^{(i)}} \cdot \nabla \widehat{\mathfrak{m}^{(i)}} \right)}{\sum\_{i=1}^{2} \left( \nabla \widehat{\mathfrak{m}^{(i)}} \cdot \nabla \widehat{\mathfrak{m}^{(i)}} \right) + e^{2}} dv. \tag{33}$$

The joint focusing constraint is not applied until the data misfit functionals have been sufficiently minimized (e.g., *χ* <sup>2</sup> = 2). This avoids the introduction of inversion artifacts arising from the focusing constraint.

#### *4.3. Representation of a Stabilizing Functional in the Form of a Pseudo-Quadratic Functional*

In this and the following sections, for simplicity, we will drop the tilde sign above the dimensionless weighted data and model parameters, assuming that they have already been normalized according to Formulas (23) and (24).

The focusing stabilizing functionals introduced above can be represented as pseudoquadratic functionals of the model parameters [9,17]:

$$\begin{split} S(\mathfrak{m}) &= \left( \mathcal{W}\_{\varepsilon}(\mathfrak{m} - \mathfrak{m}\_{\mathfrak{a}pr}), \mathcal{W}\_{\varepsilon}(\mathfrak{m} - \mathfrak{m}\_{\mathfrak{a}pr}) \right)\_{L\_{2}} \\ &= \int\_{V} \left| \mathcal{W}\_{\varepsilon}(\mathfrak{r}) \left( \mathfrak{m}(\mathfrak{r}) - \mathfrak{m}\_{\mathfrak{a}pr}(\mathfrak{r}) \right) \right|^{2} dv\_{\prime} \end{split} \tag{34}$$

where *W<sup>e</sup>* is a linear operator of the product of the model parameters function, *m*(*r*), and the weighting function, *we*(*r*), which may depend on *m*. If operator *W<sup>e</sup>* does not depend on *m*(*r*), we obtain a quadratic functional, i.e., a minimum norm-stabilizing functional. In general cases, *w<sup>e</sup>* may even be a nonlinear function of *m*, like the minimum support (32) or minimum gradient support (33) functionals. In these cases, the functional *S*(*m*), determined by Formula (34), is not quadratic. That is why we call it a "pseudo-quadratic" functional. However, it was shown in [9,17] that presenting a stabilizing functional in a

pseudo-quadratic form enables the development of a unified approach to regularization with different stabilizers and simplifies the solution of the regularization problem.

For example, the joint minimum support functional, *sJMS*(*m*), can be written as follows:

$$\begin{split} \mathbf{S}\_{IMS} &= \sum\_{i=1}^{2} \int\_{V} \frac{\left(\mathfrak{m}^{(i)} - \mathfrak{m}\_{apr}^{(i)}\right)^{2}}{\sum\_{j=1}^{2} \left(\mathfrak{m}^{(j)} - \mathfrak{m}\_{apr}^{(j)}\right)^{2} + \varepsilon^{2}} d\upsilon \\ &= \sum\_{i=1}^{2} \left(\mathcal{W}\_{\varepsilon} \left(\mathfrak{m}^{(i)} - \mathfrak{m}\_{apr}^{(i)}\right), \mathcal{W}\_{\varepsilon} \left(\mathfrak{m}^{(i)} - \mathfrak{m}\_{apr}^{(i)}\right)\right)\_{\mathcal{L}\_{2}} \end{split} \tag{35}$$

where *W<sup>e</sup>* is a linear operator of the multiplication of the model parameter function, *m*(*r*), by the following weighting function, *we*(*r*):

$$w\_{\varepsilon}(\mathbf{r}) = w\_{\varepsilon}^{\text{JMS}}(\mathbf{r}) = \frac{1}{\left[\sum\_{j=1}^{2} \left(m^{(j)}(\mathbf{r}) - m\_{apr}^{(j)}(\mathbf{r})\right)^{2} + \varepsilon^{2}\right]^{1/2}}.\tag{36}$$

For the joint minimum gradient support functional, *sJMGS*(*m*), we assume *mapr* = 0, and find

$$\begin{split} S\_{IMGS} = \sum\_{i=1}^{2} \int\_{V} \frac{\left(\nabla \mathfrak{m}^{(i)} \cdot \nabla \mathfrak{m}^{(i)}\right)}{N\left(\nabla \mathfrak{m}^{(i)} \cdot \nabla \mathfrak{m}^{(i)}\right) + \varepsilon^{2}} d\upsilon \approx \sum\_{i=1}^{2} \int\_{V} \left| w\_{\varepsilon}^{(i)}(r) \mathfrak{m}^{(i)}(r) \right|^{2} dr \\ = \sum\_{i=1}^{N} \left( \mathbf{W}\_{\varepsilon}^{(i)} \mathfrak{m}^{(i)} \prime \mathbf{W}\_{\varepsilon}^{(i)} \mathfrak{m}^{(i)} \right)\_{L\_{2}} \,' \end{split} \tag{37}$$

where *W* (*i*) *e* is a linear operator of the multiplication of the model parameter function, *m*(*r*), by the following weighting function, *w* (*i*) *<sup>e</sup>* (*r*) :

$$w\_{\varepsilon}^{(i)}(r) = w\_{\varepsilon}^{(i)IMGS}(r) = \frac{\nabla \mathfrak{m}^{(i)}(r)}{\left(\sum\_{j=1}^{2} \left(\nabla \mathfrak{m}^{(j)} \cdot \nabla \mathfrak{m}^{(j)}\right) + \varepsilon^{2}\right)^{1/2} \left[\mathfrak{m}^{(i)2} + \varepsilon^{2}\right]^{1/2}},\tag{38}$$

and *ε* is a small number similar to a focusing parameter *e*.

Using the pseudo-quadratic form (37) of stabilizing functional, we can present the corresponding parametric functional (31) as follows:

$$P^{\alpha} \left( \mathfrak{m}^{(1)}, \mathfrak{m}^{(2)} \right) = \sum\_{i=1}^{2} \varphi \left( \mathfrak{m}^{(i)} \right) + \sum\_{i=1}^{2} \alpha^{(i)} \left( \mathbf{W}\_{\varepsilon}^{(i)} \mathfrak{m}^{(i)}, \mathbf{W}\_{\varepsilon}^{(i)} \mathfrak{m}^{(i)} \right)\_{L\_{2}}.\tag{39}$$

Therefore, the parametric functional introduced by Equation (31) can be minimized in the same manner as a conventional Tikhonov parametric functional. The only difference is the introduction of some variable weighting operator, *W* (*i*) *e* , which depends on the model parameters. A practical technique of minimizing the parametric functional (39) is based on application of the re-weighted regularized conjugate gradient (RRCG) method. Interested readers can find a detailed description of the RRCG algorithm in [9].

#### *4.4. Data and Model Weights*

Data and model weights are critical in both joint inversion scenarios presented, because they ensure that both inverse problems (gravity and magnetic) converge towards a satisfactory solution at roughly the same rate. Appropriate data and model weights also ensure that the models do not unduly bias one another where no structural corellation exists. As such, in both joint inversion scenarios, two levels of data and model weighting are employed. Initial data weights are given by the following function of the errors:

$$\mathbf{W}\_d^{(i)} = \frac{1}{\left(err\_f^{(i)}\mathbf{d}^{(i)} + err\_{\text{abs}}^{(i)}\right)}\,\tag{40}$$

where *err* (*i*) *f* are the relative errors, and *err* (*i*) *abs* are the absolute error floors. These intial data weights are used in standalone and joint inversion; however, in the joint inversions, data weights are further scaled by the initial data misfit for each data set:

$$\mathbf{W}\_{d,j}^{(i)} = \frac{\mathbf{W}\_d^{(i)}}{\boldsymbol{\varphi}\_{\text{ini}} \left(\mathbf{m}^{(i)}\right)} \,\prime \tag{41}$$

where <sup>ϕ</sup>*ini m*(*i*) are the initial misfits. The matrices of initial model weights are given by the integrated sensitivity [9,17]:

$$\mathbf{W}\_{m}^{(1)} = \operatorname{diag} \left( \mathbf{F}^{(1) \text{T}} \mathbf{F}^{(1)} \right)^{\frac{1}{4}} , \tag{42}$$

$$\mathbf{W}\_{M\_{\mathbf{Y}}}^{(2)} = \left\{ \operatorname{diag} \left( \mathbf{F}\_{\mathbf{Y}}^{(2) \mathbf{T}} \mathbf{F}\_{\mathbf{Y}}^{(2)} \right)^{\frac{1}{4}} \right\} , (\boldsymbol{\gamma} = \boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}), \tag{43}$$

where *F* (*i*) is the Fréchet derivative of *A* (*i*) *m*(*i*) and *F* (*i*)*T* is the transposed matrix. In the joint inversions, the model weights are further scaled by normalization with the maximum model contrast obtained from standalone inversions:

$$\mathbf{W}\_{m,j}^{(1)} = \frac{\mathbf{W}\_m^{(1)}}{\max\left(m\_{sep}^{(1)} - m\_b^{(1)}\right)}\tag{44}$$

$$\mathbf{W}\_{m,j}^{(2)} = \left\{ \frac{\mathbf{W}\_{M\_\gamma}^{(2)}}{\max\left(\mathbf{M}\_{sep\_\gamma}^{(2)} - \mathbf{M}\_{b\_\gamma}^{(2)}\right)} \right\} \prime (\gamma = x, y, z) \,, \tag{45}$$

where *m* (1) *sep*, *M* (*i*) *sepγ* are the density and magnetization models obtained from standalone (separate) inversions, and *m* (1) *b* , *M* (*i*) *bγ* are the background density and magnetization values.

The constants *α*, *β* are adaptively and monotonically reduced to ensure stable convergence towards a satisfactory solution [9,12]. The inversion is halted when both misfit terms have been minimized to the noise level.

#### **5. Model Study**

To illustrate the developed methods, we present the results of their application to synthetic gravity gradiometry and TMI data generated for a two-body model shown in Figure 1. The bodies are both parallelepipeds of dimensions 200 × 50 m, situated at 50 m and 100 m depths. The left body has physical properties of 0.3 g/cm<sup>3</sup> , −0.06 A/m in the X direction, and 0.06 A/m in the Z direction. The right body has physical properties of 0.4 g/cm<sup>3</sup> , 0.06 A/m in the X direction, and 0.06 A/m in the Z direction. The magnetization vectors for both bodies have inclinations of −45◦ and opposing declinations of −90◦ and 90◦ , or east/west. The inducing field has inclination 75.5◦ and declination −6.2◦ . The Gxx, Gyy, and Gzz gravity gradiometry components and TMI data are simulated at equidistant stations laterally spaced at 25 m with a flight height of 30 m, also shown in Figure 1. All data are contaminated with 5% Gaussian noise.

90<sup>∘</sup> 75.5<sup>∘</sup> −6.2<sup>∘</sup>

() = ቊ൫<sup>ஓ</sup>

, (ଵ) =

(ଶ) = ቐ ெ<sup>ಋ</sup>

max ቀ௦<sup>ಋ</sup>

()

,

 (ଵ) , ം ()

௦ (ଵ)

, ௦ം ()

,

(ଵ) = ൫()

(ଶ) ஓ (ଶ) ൯ ଵ ସ

() ൫() (ଵ) ൯ ଵ ସ ,

൯

(ଵ) − (ଵ) ൯ ,

> (ଶ) ቁ

 (ଵ)

max൫௦

(ଶ) − <sup>ಋ</sup>

(ଶ)

ቋ , (γ = , , ),

()்

ቑ , (γ = , , ),

−

−45<sup>∘</sup> −90<sup>∘</sup>

**Figure 1.** Panels (**a**,**b**) show 3D views of the true density model and vertical component of magnetization vector model, respectively. The black dots represent the location of the observation stations.

We have inverted the synthetic data using the standalone and joint inversions, introduced in the previous sections. Note that, in this model study and in the following case study, we have inverted magnetic data for all three scalar components of the magnetization vector, which is the principal distinguished feature of the developed method in comparison to the conventional approach based on inversion for magnetic susceptibility only. At the same time, there are several possible approaches for how to calculate the joint focusing functionals required for joint inversion. For example, we can consider three joint focusing stabilizers—each involving a scalar component of magnetization vector and density—or jointly focus the amplitude of magnetization vector and density. We have found that jointly focusing the dominant scalar component and density works the best. Therefore, we only include the vertical component of magnetization vector, *Mz*, in the focusing stabilizers. Based on the results of the standalone magnetization vector inversion, the vertical component of magnetization vector, *Mz*, is dominant in this area. As such, the joint focusing term includes only the density, *m*(1) = *ρ*, and the vertical component of magnetization vector, *m*(2) = *Mz*, to reduce non-uniqueness.

The inversion domain was discretized into rectangular cells with a horizontal cell size of 25 × 25 m and with 20 logarithmically spaced vertical layers for a total cell count of ~45,000. All inversions were terminated when both data sets reach the simulated noise level *χ* <sup>2</sup> = 1 .

Vertical sections of the true model and inverse solutions from standalone and joint inversions are shown in Figure 2. They illustrate that the joint inversions produce a better image of the model than the standalone inversions while maintaining a similar level of data misfit. Observed and predicted data for the Gzz component of gravity gradiometry and the TMI data for the various inversion scenarios are shown in Figure 3. We can see that the joint-focused inversion slightly improves data fidelity compared to the other inversion types. Between the joint inversion scenarios, the Gramian and joint-focused inversions arrive at similar models and data fits; however, inversion with joint focusing is three times faster with respect to both iterations required and computational speed.

<sup>௭</sup> , (ଵ) = ,

(ଶ) = <sup>௭</sup>

( <sup>ଶ</sup> = 1) ,

<sup>௭</sup> ,

**Figure 2.** Panels (**a**,**b**) show the true density model and magnetization vector model, respectively. Panels (**c**,**d**) show vertical sections of the density and magnetic vector models inverted standalone, respectively. Panels (**e**,**f**) show vertical sections of the density and magnetic vector models inverted with Gramian constraints, respectively. Panels (**g**,**h**) show vertical sections of the density and magnetic vector models inverted with joint focusing constraint, respectively. Panels (**b**,**d**,**f**,**h**) show the vertical component of the magnetic vector in color; the full magnetic vector is shown by black and white arrows, and the direction of the inducing field is shown by red arrows.

**Figure 3.** Model study. Panels (**a**,**b**) show observed and predicted Gzz gravity gradient data inverted standalone, respectively. Panels (**c**,**d**) show observed and predicted total magnetic intensity (TMI) data inverted standalone, respectively. Panels (**e**,**f**) show observed and predicted Gzz gravity gradient data inverted with Gramian constraints, respectively. Panels (**g**,**h**) show observed and predicted TMI data inverted with Gramian constraints, respectively. Panels (**i**,**j**) show observed and predicted Gzz gravity gradient data inverted with joint focusing, respectively. Panels (**k**,**l**) show observed and predicted TMI data inverted with joint focusing, respectively.

#### **6. Case Study of Joint Inversions of Airborne Gravity Gradiometer (AGG) and Magnetic Data**

*6.1. Regional Geology in the Ring of Fire*

McFaulds Lake covers the Ring of Fire intrusive complex located in the James Bay lowlands of northwestern Ontario. Ring of Fire is a roughly north–south trending Archean green belt (Figure 4) consisting of mafic and ultramafic rocks. Due to flat topography and Paleozoic carbonate rocks cover, the area was not extensively studied until the early 2000s, when, as a part of kimberlite exploration campaign, McFaulds volcanic massive sulfide (VMS) deposits were discovered. The recognition that the area hosts a greenstone belt has triggered an exploration rush which resulted in three major types of economic mineral deposit discoveries, including magmatic Ni-Cu-PGE, magmatic chromite mineralization, and volcanic massive sulfide mineralization.

**Figure 4.** Geological map of the Ring of Fire area with marked known deposits (from [20]).

The Ring of Fire complex is primarily composed of mafic metavolcanic flows, felsic metavolcanic flows, and pyroclastic rocks. The Thunderbird deposit is associated with one of the numerous layered mafic-to-ultramafic intrusions, which trend subparallel with and obliquely cut the westernmost part of the belt, close to a large granitoid batholith lying west of the belt [21]. These layered intrusions are associated with high-grade Ni-Cu-PGE, Chromite, and V-Ti-Fe deposits, with Thunderbird being the latter.

#### *6.2. Airborne Gravity and Magnetic Surveys*

In the exploration rush of the mid-2000s, there were numerous smaller-scale exploration projects performed; however, the overall picture of the area was lacking. In order to map regional geology and locate further potential mineral resources, an airborne geophysical survey was carried out in the McFaulds Lake region between 2010 and 2011 using the Fugro Airborne Surveys system. Both airborne gravity gradiometer (AGG) and magnetic data were collected. This project was collaboratively operated between the Ontario Geological Survey (OGS) and the Geological Survey of Canada (GSC). Figure 5 shows the McFaulds Lake region survey block. The survey parameters included:


**Figure 5.** Airborne geophysical survey block in the McFaulds Lake region. airborne gravity gradiometry (AGG) and magnetic data were collected through the Fugro AGG and magnetic system (from [21]).

#### *6.3. Thunderbird Deposit*

The Thunderbird deposit consists of semi-massive vanadium- and titanium-enriched magnetite, which corresponds to a strong gravity and magnetic anomaly. Noront Resources, owner of the claim, estimates the ore body at 1.6 km long, 400 m wide, and 500 m deep, based on gravity and magnetic data and limited core drilling. This deposit has not been developed yet, and a more detailed analysis is lacking.

We have inverted the Gxy, Gxx, Gyy, and Gzz gravity gradient components and TMI data separately and jointly. As discussed above, the filtered RTP magnetic data were used for susceptibility inversion, while the inversion for magnetization vector was based on the TMI data with the removed regional trend. As an example, panel (a) of Figure 6 shows the Gzz component of the observed AGG data, while panel (b) presents the filtered TMI data map of the Thunderbird area. The deposit is assumed to be a semi-massive V-Ti-enriched magnetite with an approximate volume of 0.32 km, which is based on potential field data analysis and limited core drilling.

The inversion grid was discretized with a 50 × 50 m horizontal spacing and a logarithmic vertical discretization ranging from 25–150 m. The total number of cells inverted was ~250,000 cells. A 16-core Intel Xeon workstation with 128 GB memory was used to run all of the inversions. The total inversion runtime was ~10 min for the standalone AGG and TMI inversions, ~15 min for the joint inversion with joint focusing stabilization, and ~45 min for the joint inversion with Gramian constraints.

**Figure 6.** Panel (**a**) shows the Gzz component of the observed AGG data in UTM (16N) coordinates in panel (**a**). Panel (**b**) presents the filtered TMI data map with the removed regional trend. Profile AA' is shown by the black line. The location of drill hole NOT-09-2G25 is shown by the black asterisk.

A comparison of vertical sections of the 3D models obtained from the various inversion methodologies is shown in Figure 7. It is clear that the models obtained from a joint inversion with either Gramian or joint focusing constraints have sharper geospatial boundaries and a higher degree of structural correlation, when compared to the models obtained from standalone inversions. It is important to note that all inverted 3D models share the same level of data misfit.

**Figure 7.** Panels (**a**,**b**) show vertical sections of the 3D density and magnetization vector models obtained from standalone inversions, respectively. Panels (**c**,**d**) show vertical sections of the 3D density and magnetization vector models obtained from joint inversion with Gramian constraints, respectively. Panels (**e**,**f**) show vertical sections of the 3D density and magnetizations vector models obtained from joint inversion with joint focusing constraint, respectively. Panels (**b**,**d**,**f**) show the vertical component of magnetization vector in color; the full magnetization vector is shown by the black arrows, and the inducing field direction is shown by the red arrows.

We have included both the full magnetization vector and the inducing field direction in the magnetic panels for clarity. Drill hole NOT-09-2G25, drilled by Noront Resources Limited in 2009 [22] and situated in the center of these profiles, found a gabbro schist unit from 21–462 m containing magnetite and titanomagnetite with concentrations up to 10% and smaller pockets containing up to 50% magnetite. This matches quite well with the jointly inverted results. Moreover, an iron formation with up to 20% magnetite concentration was found at 166–170 m which may be what the focused inversion is resolving in the near surface.

A comparison of selected observed and predicted data is shown in Figure 8. All inversion methodologies reach an excellent data fit. Figure 9 shows model parameter cross plots for the different inversion scenarios. Additionally, the correlation coefficient was calculated for the density model versus the vertical component of the magnetization vector and is summarized in Table 1.

**Figure 8.** Case study. Panels (**a**,**b**) show observed and predicted Gzz gravity gradient data inverted standalone, respectively. Panels (**c**,**d**) show observed and predicted TMI data inverted standalone, respectively. Panels (**e**,**f**) show observed and predicted Gzz gravity gradient data inverted jointly with Gramian constraints, respectively. Panels (**g**,**h**) show observed and predicted TMI data inverted jointly with Gramian constraints, respectively. Panels (**i**,**j**) show observed and predicted Gzz gravity gradient data inverted with joint focusing, respectively. Panels (**k**,**l**) show observed and predicted TMI data inverted with joint focusing, respectively.

Figures 10 and 11 show a comparison between the results of the susceptibility-only inversion and magnetization-vector inversion. The left panel in Figure 10 shows the vertical section of the inverse model for susceptibility inversion, while the right panel shows the vertical section of the inverse model for magnetization vector inversion. The small black arrows in the right panel represent the magnetization vectors, while the color bar reflects the vertical component of magnetization.

**Figure 9.** Model parameter cross plots shown in Panels (**a**–**c**) correspond to the standalone inversions, joint Gramian inversion, and joint focusing inversion, respectively.


**Figure 10.** Left panel (**a**) shows the vertical section of the inverse model for susceptibility inversion. Right panel (**b**) shows a vertical section of the inverse model for magnetization vector inversion. The small black arrows in the right panel represent the magnetization vectors, while the color bar reflects the vertical component of magnetization.

**Figure 11.** Comparison of 3D inverse models with drilling results. Left panel (**a**) shows the volume **Figure 11.** Comparison of 3D inverse models with drilling results. Left panel (**a**) shows the volume image of the inverse susceptibility model. Right panel (**b**) presents the volume image of the vertical component of the magnetization vector. The black–yellow–black solid line shows the location of the borehole drilled in the survey area. The yellow color indicates the mineralization zone confirmed by drilling.

Figure 11 shows a comparison of 3D inverse models with drilling results. The left panel (a) presents the volume image of the inverse susceptibility, while the right panel (b) shows a similar image of the vertical component on the magnetization vector. The black– yellow–black solid line is the borehole drilled in the survey area. The yellow color indicates the mineralization zone confirmed by drilling. One can see that the image produced by inversion for the magnetization vector correctly represents the location of the mineralization zone. In other words, the volume images also show that you can better pinpoint anomalous bounds with magnetization vector in the presence of remanent magnetization for the Thunderbird Vanadium-Titanium target in Ontario, Canada.

Finally, Figure 12 presents a comparison of the inverse density (top panel) and vertical component of magnetization models (bottom panel) recovered by joint focusing inversion. It was shown in [11] that it can be challenging to resolve a steeply dipping magnetic dike and recover a good data fit when inverting for magnetic susceptibility only. By inverting for the full magnetization vector, it is much easier to fit these structures and still retain a good data fit. Including these additional inversion parameters increases the non-uniqueness of the potential field inverse problem; however, we have remedied this with the addition of the Gramian and joint focusing constraints.

**Figure 12.** Top panel (**a**) shows the cut-off 3D image of the inverse density model. Bottom panel (**b**) presents the volume image of the vertical component of the magnetization vector. The black– yellow–black solid line shows the location of the borehole drilled in the survey area. The yellow color indicates the mineralization zone confirmed by drilling.

#### **7. Conclusions**

Magnetic data contain information about the remanent magnetization of rock formations that are often ignored when inverting for magnetic susceptibility only. We have introduced two novel methods of joint inversion of potential field data in the presence of remanent magnetization: (1) joint inversion subject to the Gramian constraint and (2) joint inversion subject to the joint focusing constraint. We have demonstrated the methods in a model study, and in a case study using airborne data acquired over the Thunderbird V-Ti-Fe deposit in the Ring of Fire area of Ontario, Canada. We have presented a comparison of the 3D models obtained using conventional single-datum standalone inversion and joint inversion subject to either the Gramian constraint or the joint focusing constraint. Comparison of the different inversion methodologies illustrates that 3D models obtained from joint inversion with either the Gramian or joint focusing constraint were able to resolve a shared earth model that featured sharper anomalous geospatial boundaries and a higher degree of structural correlation, as opposed to the 3D models obtained from standalone inversions. The joint inversion subject to the Gramian constraint did yield the highest degree of structural correlation; while the joint inversion subject to the joint focusing constraint provided comparable results but was much faster. By jointly inverting AGG and TMI data for density and magnetization vector, we demonstrated the efficiency of the novel methods in defining rock formations with remanent magnetization which are otherwise difficult to recover.

**Author Contributions:** Conceptualization, M.S.Z. and M.J.; methodology, M.S.Z.; software, M.J.; validation, M.J., and M.S.Z.; formal analysis, M.S.Z.; investigation, M.J.; resources, M.S.Z.; writing—original draft preparation, M.J.; writing—review and editing, M.S.Z.; visualization, M.J.; supervision, M.S.Z.; project administration, M.S.Z.; funding acquisition, M.S.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by Consortium for Electromagnetic Modeling and Inversion (CEMI), and TechnoImaging and received no external funding.

**Data Availability Statement:** The data used in this study are publicly available from the Ontario Geologic Survey.

**Acknowledgments:** The authors acknowledge the Consortium for Electromagnetic Modeling and Inversion (CEMI) at the University of Utah and TechnoImaging for the support of this research. The AGG and TMI data were collected by Fugro and made available by the Ontario Geological Survey, Canada.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Minerals* Editorial Office E-mail: minerals@mdpi.com www.mdpi.com/journal/minerals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-1739-1