MAGNETIC INVERSION TO RECOVER COMPACT STRUCTURES BY AN ARCTANGENT STABILIZER

2.1. Forward modelling

The magnetic data observed at Earth's surface are primarily caused by induced magnetization effects. Forward modelling aims to compute this magnetic field, which arises from the distribution of magnetic susceptibility (κ) in subsurface materials. The total magnetic field component () at an observation point r₀ outside the magnetic source region can be expressed as (Blakely, 1995):

[TeX:] \overset{⃑}{T}(r_{0}) = - C\nabla\int_{R}^{}{J\left( r \right).\nabla\frac{1}{\left| r_{0} - r \right|}\text{dv}}$\stackrel{⃑}{T}(r_0)=-C\nabla {\int }_R^{}J\left(r\right).\nabla \frac{1}{|r_0-r|}\mathrm{\text{dv}}$ (1)

where C=1×10^-7 (H.m^-1), R is the volume occupied by the causative body, r₀ denotes the observation point, r represents the position of a volume element dv of the source body and J=J_i+J_r is the total magnetization vector that is the sum of induced and remanent magnetizations. Often, we ignore the remanent component and the total magnetization vector can be obtained by the following formula:

[TeX:] J_{i} = \kappa H$J_i=\kappa H$ (2)

where H is the Earth's magnetic field and κ represents magnetic susceptibility. To compute the total magnetic field component in Equation 1, the subsurface beneath the survey area is discretized into a series of rectangular prismatic cells. Each prism has predefined dimensions and location, with an assigned uniform magnetic susceptibility value. The foundational method for calculating the magnetic response of such rectangular prisms was developed by Bhattacharyya (1964), with later refinements by Rao and Babu (1991) improving computational efficiency for modern computer-based implementations. In this study, we compute the magnetic signature of individual prismatic bodies using the Rao and Babu (1991) formulation. Blakely (1995) further contributed a concise mathematical representation of surface magnetic field data, compactly expressing the relationship between subsurface magnetic properties and surface observations as follows:

[TeX:] d_{N \times 1} = G_{N \times M} \times \kappa_{M \times 1}$d_{N\times 1}=G_{N\times M}\times {\kappa }_{M\times 1}$ (3)

where N represents the total number of measurement locations at the surface, and M denotes the number of rectangular prisms used to discretize the subsurface in the model. The vector κ contains the unknown model parameters, which in this context represent the magnetic susceptibilities of each subsurface prism. The vector d comprises the measured magnetic field data from the survey, while matrix G serves as the forward operator that describes how each subsurface prism contributes to the magnetic field at each observation point.

2.2. Inversion method

To perform quantitative magnetic data inversion, we construct an objective function comprising two key components: (1) a data misfit term measuring the discrepancy between observed and predicted data, and (2) a regularization term controlling model characteristics. The data misfit quantifies how well the inversion results match field measurements, while the regularization term constrains the solution to exhibit desired physical properties. Following standard practice in geophysical inversion, we employ Tikhonov regularization (Tikhonov & Arsenin, 1977) to ensure stable solutions:

[TeX:] \Phi\left( \kappa \right) = \left\| W_{d}\left( \text{Gκ} - d \right) \right\|_{2}^{2} + \lambda S\left( \kappa \right)$\Phi \left(\kappa \right)={\parallel W_d(\mathrm{\text{Gκ}}-d)\parallel }_2^2+\lambda S\left(\kappa \right)$ (4)

where 𝜆 is the regularization parameter. S(κ) represents the stabilizing functional and the first term is weighted data misfit, W_d represents the weighting matrix associated with the data:

[TeX:] W_{d} = diag\left( \frac{1}{\sigma_{1}},....,\frac{1}{\text{σi}} \right)$W_d=diag(\frac{1}{{\sigma }_1},....,\frac{1}{\mathrm{\text{σi}}})$ (5)

where σᵢ represents the standard deviation of the i^th measurement. This study investigates a compactness constraint based on an arctangent function, designed to improve the delineation of sharp boundaries between subsurface structures. The stabilizing functional can be expressed as:

[TeX:] S\left( \kappa \right) = \left\| W_{e}\kappa \right\|_{2}^{2}$S\left(\kappa \right)={\parallel W_e\kappa \parallel }_2^2$ (6)

where W_e is a diagonal model weighting matrix which is dependent on model parameters and can be given by (Hu et al., 2017):

[TeX:] W_{e} = diag\left( \sqrt{\frac{\arctan\sqrt{\kappa^{2} + \mu^{2}}}{\kappa^{2} + \varepsilon^{2}}} \right)$W_e=diag\left(\sqrt{\frac{arctan\sqrt{{\kappa }^2+{\mu }^2}}{{\kappa }^2+{\epsilon }^2}}\right)$ (7)

where ε and μ are small, non-zero constants introduced to prevent mathematical singularities and ensure numerical stability. While these parameters serve this essential computational purpose, they have negligible impact on the model parameters themselves. The arctangent functional provides an efficient focusing inversion method for magnetic data that offers two key advantages: (1) it eliminates the need for a focusing parameter, significantly simplifying implementation, and (2) it maintains robust performance in recovering subsurface structures.

Li and Oldenburg (1998) demonstrated that magnetic data inherently possesses limited depth resolution due to the fundamental nature of potential field measurements. This characteristic makes precise depth determination of magnetic sources particularly challenging. To mitigate the natural tendency of susceptibility models to concentrate near the surface, we incorporate a depth weighting matrix defined as:

[TeX:] W_{z} = diag\left\lbrack \frac{1}{z^{\frac{\eta}{2}}} \right\rbrack$W_z=diag\left[\frac{1}{z^{\frac{\eta }{2}}}\right]$ (8)

where z is the depth of the model parameter and η­ is set to 3 in magnetic inversion (Li & Oldenburg, 1998). Moreover, the algorithm must ensure that the model adheres to the specified upper and lower susceptibility limits until a satisfactory model is obtained (Portniaguine & Zhdanov, 1999, 2002).

[TeX:] \kappa = \left\{ \begin{matrix} b\quad\text{if}\mspace{6mu}\kappa < b \\ a\quad\text{if}\mspace{6mu}\kappa > a \\ \end{matrix} \right.\ $\kappa =\{\begin{array}{c}\hfill b\mathrm{\text{if}}\kappa <b\hfill \\ \hfill a\mathrm{\text{if}}\kappa >a\hfill \end{array}$ (9)

where b and a represent the lower and upper bounds of magnetic susceptibility, respectively. We solve this magnetic focusing inversion problem through an iterative approach using the reweighted regularized conjugate gradient (RRCG) method (Zhdanov, 2015), implemented as follows:

[TeX:] r_{n} = W_{d}\left( G\left( \kappa_{n} \right) - d, \right),\quad s_{n} = W_{\text{en}}W_{z}\left( \kappa_{n} \right),$r_n=W_d(G\left({\kappa }_n\right)-d,),s_n=W_{\mathrm{\text{en}}}{W}_z\left({\kappa }_n\right),$

[TeX:] I_{(k_{n})}^{\lambda} = I_{(k_{n})}^{\lambda_{n}} = G^{T}W_{d}^{2}r_{n} + \lambda_{n}W_{\text{en}}W_{z}s_{n},$I_{(k_n)}^{\lambda }=I_{(k_n)}^{{\lambda }_n}=G^T{W}_d^2{r}_n+{\lambda }_n{W}_{\mathrm{\text{en}}}{W}_z{s}_n,$

[TeX:] \beta_{n + 1}^{\lambda} = \frac{\left\| I_{\left( k_{n + 1} \right)}^{\lambda} \right\|^{2}}{\left\| I_{\left( k_{n} \right)}^{\lambda} \right\|^{2}},\quad{\widetilde{I}}_{(k_{n})}^{\lambda_{n}} = I_{(k_{n})}^{\lambda_{n}} + \beta_{n}^{\lambda_{n}}{\widetilde{I}}_{(k_{n - 1})}^{\lambda_{n - 1}},\quad{\widetilde{I}}_{0}^{\lambda_{0}} = I_{0}^{\lambda_{0}},${\beta }_{n+1}^{\lambda }=\frac{{\parallel I_{\left(k_{n+1}\right)}^{\lambda }\parallel }^2}{{\parallel I_{\left(k_n\right)}^{\lambda }\parallel }^2},{\tilde{I}}_{(k_n)}^{{\lambda }_n}=I_{(k_n)}^{{\lambda }_n}+{\beta }_n^{{\lambda }_n}{\tilde{I}}_{(k_{n-1})}^{{\lambda }_{n-1}},{\tilde{I}}_0^{{\lambda }_0}=I_0^{{\lambda }_0},$

[TeX:] \omega_{n}^{\lambda_{n}} = \frac{\left( {\widetilde{I}}_{(k_{n})}^{{\lambda_{n}}^{T}}I_{(k_{n})}^{\lambda_{n}} \right)}{{\widetilde{I}}_{(k_{n})}^{{\lambda_{n}}^{T}}\left( G^{T}W_{d}^{2}G + \lambda W_{\text{en}}^{2}W_{k}^{2} \right)I_{(k_{n})}^{\lambda_{n}}},${\omega }_n^{{\lambda }_n}=\frac{\left({\tilde{I}}_{(k_n)}^{{{\lambda }_n}^T}{I}_{(k_n)}^{{\lambda }_n}\right)}{{\tilde{I}}_{(k_n)}^{{{\lambda }_n}^T}(G^T{W}_d^{2}G+\lambda W_{\mathrm{\text{en}}}^2{W}_k^2)I_{(k_n)}^{{\lambda }_n}},$

[TeX:] \kappa_{n + 1} = \kappa_{n} - \kappa_{n}^{\lambda_{n}}{\widetilde{I}}_{(k_{n})}^{\lambda_{n}},\quad\gamma = \frac{\left\| s_{n + 1} \right\|^{2}}{\left\| s_{n} \right\|^{2}},${\kappa }_{n+1}={\kappa }_n-{\kappa }_n^{{\lambda }_n}{\tilde{I}}_{(k_n)}^{{\lambda }_n},\gamma =\frac{{\parallel s_{n+1}\parallel }^2}{{\parallel s_n\parallel }^2},$

[TeX:] \lambda_{1} = \frac{\left\| W_{d}\left( G\kappa_{1} - d \right) \right\|_{2}^{2}}{\left\| \kappa_{1} \right\|_{2}^{2}}${\lambda }_1=\frac{{\parallel W_d(G{\kappa }_1-d)\parallel }_2^2}{{\parallel {\kappa }_1\parallel }_2^2}$

[TeX:] \lambda_{n + 1} = \lambda_{n},\quad\text{if}\mspace{6mu}\gamma \leq 1,\quad\text{and}\quad\lambda_{n + 1} = \frac{\lambda_{n}}{\gamma},\quad\text{if}\text{\!\;}\gamma > 1

[TeX:] \lambda_{n + 1}^{'} = q\lambda_{n + 1},\quad q < 1,\mspace{6mu}\text{if}\mspace{6mu}\left\| W_{d}r_{n} \right\|^{2} - \left\| W_{d}r_{n + 1} \right\|^{2} < 0.01\left\| W_{d}r_{n} \right\|^{2}${\lambda }_{n+1}^{\prime }=q{\lambda }_{n+1},q<1,\mathrm{\text{if}}{\parallel W_d{r}_n\parallel }^2-{\parallel W_d{r}_{n+1}\parallel }^2<0.01{\parallel W_d{r}_n\parallel }^2$ (10)

The algorithm converges to the solution when the normalized misfit arrives at the given level δ₀ (Rezaie, 2019, 2020):

[TeX:] \frac{\left\| W_{d}\left( G\kappa_{n} - d \right) \right\|_{2}^{2}}{\left\| W_{d}d \right\|_{2}^{2}} \leq \delta_{0}$\frac{{\parallel W_d(G{\kappa }_n-d)\parallel }_2^2}{{\parallel W_{d}d\parallel }_2^2}\le {\delta }_0$ (11)

The RRCG is an efficient iterative solver that does not require constructing large M × M matrices or computing their inverses, unlike methods such as Gauss-Newton, Levenberg-Marquardt, and other IRLS algorithms. Another advantage of the RRCG algorithm is its ability to automatically select the regularization parameter through an adaptive process during iterations. This avoids the need to solve the inverse problem multiple times in each iteration.

2.3. Synthetic model

To illustrate the application of our proposed method, we examine a synthetic data scenario featuring distinct anomalous structures. According to Equation 7, we need to select two parameters μ and ε, that we set both equal to 1 × 10^-16 for all models in this paper.

The synthetic model consisted of four blocks with a susceptibility of 0.05 SI and the background susceptibility is uniform that is equals to zero. The total magnetic field anomaly was calculated at ground level, based on an inducing magnetic field which has the inclination of 65°and the declination of -25°. The inducing field strength is 50,000 nT. The model covers an area extending 90 m to the north and 90 m to the east. The subsurface was divided into 23328 equal-sized rectangular prisms for inversion. Each prism measured 12.5 m on each side, arranged in a 36 × 36 × 18 grid. There are four buried blocks in the area that the size of the buried blocks are expressed in Table 1.

Table 1. Size parameters of the buried prisms in the synthetic model

Figure 1a presents the calculated magnetic anomaly for the synthetic model, incorporating 5% random noise. The inversion process converged after 23 iterations in 6.3 s, achieving a normalized misfit of 0.05 (see Figure 1b).

Figure 1. The magnetic anomaly from the synthetic model with 5% Gaussian noise of each datum magnitude (a). The behaviour of the normalized misfit (b).

Figure 2 displays horizontal and vertical cross-sections of the true and inverted 3D susceptibility models. The horizontal slice is extracted at 75 m depth, while the vertical section is taken along the x-axis at y = 450 m. Together, these orthogonal views provide complementary visualization of the subsurface structure. The results demonstrate that the suggested inversion technique effectively reconstructs the synthetic model. In comparison with the true models in panels a and b, the recovered model accurately locates the four blocks in their intended positions, with distinct boundaries. Additionally, the estimated susceptibility values closely align with the actual values of the true blocks.

The model's computed magnetic response is presented in Figure 3a. A comparison between this calculated data and the measured data is illustrated in Figure 3b. This comparison demonstrates that the calculated data achieves an acceptable match to the observed data, indicating a successful data fitting process.

Figure 2. (a) Horizontal plan section of the true susceptibility distribution at a depth of 75 m. (b) A cross-sectional slice of the true susceptibility model at Y = 450 m. (c) Plan section through recovered model by inversion. (d) The vertical slice in the inversion result.

Figure 3. (a) The map of the magnetic anomaly from the focused model. (b) The residual map between computed and original data.

2.4. Application to the field data

The McFaulds Lake area in Canada is well-known for deposits, rich in nickel, copper, and platinum group elements (Ni-Cu-PGE) such as the Eagle's Nest deposit. This site lies near Webequie, northwestern Ontario,. The area is characterized by a highly magnetic ultramafic intrusion of mantle origin, recognized as the 'Ring of Fire'. This intrusion has been emplaced along the edge of a significant granodiorite pluton, within the rock formations of the Sachigo greenstone belt. The area has been the focus of exploration for various economically significant mineral deposits including magmatic deposits of Ni, Cu, and PGE, Chromite mineralization of magmatic origin, Volcanic-hosted massive sulfide deposits as well as Diamond-bearing kimberlite formations (Balch et al. 2010; Lin and Zhdanov, 2018). The Ring of Fire's rocks include mafic metavolcanic flows, felsic metavolcanic flows, pyroclastic rocks and layered mafic to ultramafic intrusions. These intrusions run roughly parallel to the belt's westernmost section, occasionally cutting across it at an angle. They are situated near a large granitoid batholith west of the belt. The main layered intrusion in this area is notable for two types of deposits, at its deeper parts there are high-grade Ni-Cu-PGE deposits and the shalower parts and further east contains stratiform chromite deposits in the intrusion's stratigraphy (Geological Survey and Geological Survey of Canada, 2011).

In the study area, chromite deposits are typically found within layered ultramafic intrusive rocks, often in association with magnetite and serpentine. This association allows for the magnetic tracing of chromite through its potential host rocks. Therefore, airborne geophysical surveys were carried out in the McFaulds Lake area in 2010, gathering magnetic data from the air (Geological Survey and Geological Survey of Canada, 2011). Figure 4 displays a geological map of the study area with known mineralizations (Laudadio et al. 2022).

Figure 4. Geological map of the study area (adapted from Laudadio et al. 2022)

The Ring of Fire intrusive body in the study area contains two main ultramafic intrusions, the Black Thor intrusion (BTI) and the Double Eagle intrusion (DEI), which are associated with significant chromite mineralization. These intrusions host several major deposits, including Black Thor, Black Label, Big Daddy, Black Horse, Blackbird, and Black Creek (see Figure 4). Collectively, the deposits in this region hold an estimated 343 million tonnes of geological resources averaging 32% chromite. Research suggests that the BTI and DEI were supplied by separate magma conduits. The BTI and DEI intrusions are primarily made up of a thick ultramafic lower section and a thinner mafic upper section. The ultramafic portion consists of dunite, peridotite (mainly harzburgite), pyroxenite, and chromite, all metamorphosed to lower greenschist facies. Within these rocks, dunite and peridotite are largely composed of serpentinized or talc-altered olivine, making up 75–95% of their mineral content. Magnetite is widespread, appearing both as fine disseminated particles and as rims around chromite grains, a result of serpentinization processes. Chromitite mineralization appears in diverse styles, including finely laminated to thickly bedded layers, variably disseminated grains, and interlayered massive to semi-massive chromite deposits. These occurrences are typically hosted within dunite and peridotite zones. (Laudadio et al. 2022).

To demonstrate how well our proposed method works with field data, we've applied it to a portion of the magnetic data that covers the Black Creek area which is known for its chromite mineralization. Figure 5a illustrates the magnetic anomaly of the study area in which the C1 and C2 lines intersect at the location of the Black Creek chromite deposit.

Figure 5. The map of magnetic anomaly over the McFaulds Lake area. The conjunction of the C1 and C2 lines represents the position of the Black Creek deposit (a). The behaviour of the normalized misfit (b).

The dimensions of the model domain was set to 6 km × 4.5 km × 2 km in x, y and z directions and the subsurface region was divided into cube-shaped segments, each measuring 100 m on all sides leading to 60 × 45 × 20 = 54000 cells. We ran the magnetic inversion by the arctangent stabilizing functional and the algorithm has been converged to the solution after 89 iterations in 34.38 s (see Figure 5b). Figure 6 presents a comparison of actual measurements (see Figure 5a) and calculated data from the obtained model (see Figure 6a) for the total magnetic intensity field. The visual representation demonstrates a high degree of agreement between the observed and predicted data. This close correspondence is further highlighted by the inclusion of a difference map (see Figure 6b), which visually depicts the discrepancies between the measured and predicted magnetic field.

Figure 7 illustrates cross-sections in the inverted model of the study area. The panels a and b display vertical slices through the reconstructed magnetic susceptibility distribution along two specific lines, labeled C1 and C2 in Figure 5a.

In the C1 cross-section, the inversion results clearly reveal the top boundaries of the identified mineral deposit. However, the deeper portions of this body exhibit lower magnetic properties. Examining the C2 cross-section, the magnetic inversion suggests that the dike structure on the right side contains a higher concentration of magnetic minerals. This increased magnetism might be indicative of chromite deposits in the Black Creek region.

The panel c presents a horizontal cross-section of anomalous magnetic susceptibility at a depth of 300 m. The horizontal slice, along with the previously mentioned vertical cross-sections, reveals anomalous targets characterized by distinct and partially continuous boundaries.

This representation offers valuable insight into the nature and distribution of mineral deposits in the area. This approach enhances our understanding of the subsurface geological structures and potential mineral resources. The results agree with those obtained by Lin and Zhdanov (2018) using the multinary joint inversion of gravity and magnetic data over the study area.

Figure 6. The map of the predicted magnetic data from the model (a). The map of the difference between the observed and predicted data (b).

Figure 7. Vertical cross-sections of the recovered anomalous magnetic susceptibility distribution along C1 line (a) and C2 line (b). The horizontal section at the depth of 300 m (c).