## ABSTRACT

We have developed a general 3D amplitude inversion algorithm for magnetic data in the presence of self-demagnetization and remanent magnetization. The algorithm uses a nonlinear conjugate gradient (NLCG) scheme to invert the amplitude of the magnetic anomaly vector within a partial differential equation framework. Three quantities— the amplitude of the anomalous magnetic field, the analytic signal, and the normalized source strength, defined as the amplitudes of magnetic data that are weakly dependent on the magnetization direction — are inverted to recover the 3D distribution of the subsurface magnetic susceptibility. Numerical experiments indicate that our NLCG amplitude inversion algorithm has a rapid convergence rate that provides a reasonable inversion solution in the absence of knowing the total magnetization direction. High-resolution aeromagnetic data collected from the Pea Ridge iron oxide-apatite-rare earth element deposit, southeast Missouri, USA, are used to illustrate the efficacy of our amplitude inversion algorithm. This algorithm is generally applicable for tackling the large-scale inversion problem in the presence of self-demagnetization and remanent magnetization.

## INTRODUCTION

The effects of self-demagnetization on magnetic data are associated with the high susceptibility of the magnetic source. A source with susceptibility greater than 0.1 SI can lead to significant self-demagnetization effects, which cannot be ignored in magnetic data interpretation. Many previous studies indicate that strong self-demagnetization can alter the magnitude and direction of source magnetization (Jackson, 1999; Lelièvre and Oldenburg, 2006; Krahenbuhl and Li, 2017). This demagnetization effect has been widely researched, particularly in mineral exploration (Guo et al., 2011; Hillan, 2013) and unexploded ordnance detection (Altshuler, 1996; Billings, 2004; Wallace, 2006), and various inverse data modeling methods have been proposed. Sharma (1966) suggests a numerical solution to the magnetostatic problem with high susceptibility that used an integral equation formulation. Purss et al. (2001) propose an iterative 3D numerical modeling method for high-susceptibility and complex targets. Lelièvre and Oldenburg (2006) develop a finite-volume-discretization (FVD) method for modeling the total magnetostatic field and propose using the Gauss-Newton inversion algorithm to recover 3D distributions of subsurface magnetic susceptibility. This method has the advantage of inverting magnetic data caused by a high-susceptibility source. Krahenbuhl and Li (2017) thoroughly investigate self-demagnetization effects with an integral equation formulation and suggest using amplitude data to interpret self-demagnetization problems for complex source geometries.

Previous research (Li et al., 2010; Liu et al., 2015) shows that remanent magnetization also cannot be neglected in the inversion procedure. Many collected iron ore samples show that high-susceptibility ore often contains abundant magnetite with various magnetism directions (Leão-Santos et al., 2015; Liu et al., 2018). Li et al. (2010) summarize that two strategies can be used to deal with the problem of remanent magnetization in data inversion. The conventional route estimates the magnetization direction using forward modeling (Li and Oldenburg, 1996; Pilkington, 1997; Clark, 2012). This approach is more applicable for data with a homogeneous magnetization direction within an estimated bias of 15° (Shearer, 2005; Liu et al., 2015).

The alternative method converts magnetic data into a quantity that is weakly dependent on the magnetization direction. Paine et al. (2001) invert the analytic signal (AS) amplitude of magnetic anomaly data with the general 3D inversion algorithm proposed by Li and Oldenburg (1996). Stavrev and Gerovska (2000) analyze the effects of reducing the sensitivity with different quantities of transforms. Stavrev (2006) uses amplitude transforms as a tool for magnetic data inversion of a 2.5D anomalous field. Shearer and Li (2004) and Shearer (2005) invert the amplitude and total gradient of the anomalous magnetic field in the presence of strong remanent magnetization. Li et al. (2010) propose using 3D amplitude inversion in the presence of significant remanence. Pilkington and Beiki (2013) use normalized source strength (NSS) data in a standard 3D inversion algorithm (Beiki et al., 2012); their research indicates that NSS data are weakly dependent on strong remanent magnetization effects and that NSS inversion can recover a reliable source distribution. Li and Li (2014) tackle the challenge of unknown magnetization direction in aeromagnetic data interpretation using 3D amplitude inversion. Leão-Santos et al. (2015) present a case study applying amplitude inversion to an iron oxide-copper-gold deposit with strong remanent magnetization. Furthermore, Krahenbuhl and Li (2017) report, a complex spatial distribution of source susceptibility will induce spatial heterogeneity in the magnetization direction. Thus, a single estimated average magnetization direction may not be applicable in complex data interpretation.

Considering the inversion challenge of complex remanent magnetization in a high-susceptibility source, we propose directly inverting the amplitude of the anomalous magnetic field using Maxwell’s partial differential equation (PDE) and using an efficient approach to solve such a large-scale 3D magnetic data inversion problem. Many previous geophysical inverse optimization studies have shown that the nonlinear conjugate gradient (NLCG) algorithm is more efficient than the Gauss-Newton algorithm regarding computational time and memory (Rodi and Mackie, 2001; Siripunvaraporn and Egbert, 2007). Here, we propose a 3D NLCG inversion algorithm for inverting magnetic amplitude data based on a PDE framework. We consider five different amplitude quantities of magnetic data inversion. We then numerically analyze the sensitivity of these quantity transformations due to the strong remanent magnetization and the self-demagnetization effect. In addition, we invert three transformed quantities to validate the efficacy of the proposed inversion algorithm. Finally, we verify the method using field data from an iron oxide-apatite rare-earth element deposit.

## PDE MODELING AND INVERSION OF MAGNETIC AMPLITUDE DATA

The goal of this study is to recover the 3D distribution of magnetic susceptibility for a high-susceptibility source with remanent magnetization using magnetic amplitude data. Li and Oldenburg (1996) propose the general 3D inversion framework with a weighting function and a smoothing function in the integral equation domain. Li et al. (2010) propose inverting amplitude data to solve the problem of significant remanence. Lelièvre and Oldenburg (2006) develop an inversion algorithm to recover a high-magnetic-susceptibility source in the PDE domain. Our work builds on their research. To tackle the problem of self-demagnetization, a PDE solution was used for source-free magnetostatics forward modeling, which can account for self-demagnetization caused by a highly susceptible source with complex geometry. The algorithm is realized in a general 3D inversion framework to deal with the nonuniqueness of inversion. The innovation of this study is in the formulation of the amplitude data inversion in the PDE domain, using the NLCG scheme.

The regularization parameter $\beta $ controls the trade-off between the data misfit and the complexity of inverted model. In this research, the regularization parameter $\beta $ is identified through a traditional L-curve algorithm. We discuss the choice of $\beta $ in the following sections.

## AMPLITUDE SENSITIVITY ANALYSIS WITH SELF-DEMAGNETIZATION

They conclude that model physical properties influence the sensitivity analysis. Here, we use the prism divided cell to represent the subsurface source. A high-susceptibility (1 SI) cuboid model is used in the numerical experiment (Figure 1). The source extends at a depth from 4 to 8 m. The induced field is set to 50,000 nT with $I=90\xb0$ and $D=0\xb0$. The Konigsberg ratio ($Q=|J\u21c0rem|/|J\u21c0ind|$) in this experiment is set to two.

The sensitivity values of these amplitude transformations are measured by the model using varying $I$ and $D$ of remanence (Figure 2).

The sensitivity coefficients $S(I,D)$ of all transformed quantities decrease monotonically with $I$ close to that of the inducing field (90°) (Figure 2a–2e). In contrast, $S(I,D)$ show an increasing tendency when $D$ approaches the declination of the inducing field (0°) (Figure 2a–2e). We consider this as deriving from the heterogeneous distribution of the synthetic model (Figure 1) in the horizontal directions. To prove this assumption, we rotated the synthetic model (Figure 1) 90° horizontally and remeasured the sensitivity coefficients $S(I,D)$ (Figure 2g and 2h). As shown in Figure 2g and 2h, $S(I,D)$ decrease with reduced declinations $D$ after the distribution of the model has been altered. In this experiment, the quantity NSS shows a higher level of sensitivity to the remanent magnetization compared with AS and E when $I$ is smaller than 30° and $D$ locates in the high sensitivity range, as shown in Figure 2i and 2j.

In total, the source distribution did not influence the differences of the sensitivity levels of these quantities. The term $Ba$ represents a lower level of sensitivity to the remanent magnetization as compared with other transformed quantities (AS, NSS, $E$, and $Q$) (Figure 2f and 2h). The quantity $Q$ (the black curve) is the most sensitive to remanence in this experiment, so we do not process it further. The performance of quantities AS and $E$ is the same, so we only discuss AS in the following section for brevity.

To validate the performance of these methods, we compare these data for a remanent magnetization inclination of 45° (Figure 3).

In Figure 3d, quantity $Ba$ exhibits a lower level of sensitivity to the remanent magnetization as compared with other amplitude transformations. The quantities NSS, AS, $Ba$, and $Q$ (Figure 3b–3f) show clearer anomalous shifts that derive from the remanence. This agrees with the static result of Figure 2. In fact, the behaviors of different quantities greatly affect the inversion procedure, which is analyzed in detail in the next section.

## VALIDATION USING NUMERICAL DATA

The model shown in Figure 1 is used in this amplitude data inversion experiment. The inversion region is divided into a $40\xd740\xd732$ grid (a total of 51,200 1 m cells). The grid region includes two parts: an inactive region and an active region. There are 8 free-space padding cells on both sides in the $x$- and $y$-directions and 11 padding cells above the surface in the $z$-direction. Thus, the active region (model region) contains $24\xd724\xd721$ cells with unknown susceptibility values to be determined. In this study, all of the reference models and the initial inversion model are set to 0.00 and 0.01, respectively. The regularization parameter $\beta $ is set to 10 for all synthetic experiments.

Three quantities NSS, AS, and $Ba$ are analyzed. First, the total field anomaly data (Figure 3a) are inverted assuming that the effective magnetization direction is known. The recovered model is shown in Figure 4a and 4b. With the prior information of remanence, this result shows an ideal recovered model without any bias induced by the remanence. This prior information is ignored in the following inversion tests. Figure 4c–4f shows the inversion results of amplitude data NSS (Figure 3b), AS (Figure 3c) and $Ba$ (Figure 3d), respectively, assuming that only the inclination (90°) and declination (0°) of the induced field are known. These three amplitude transformations are not completely independent of the remanence magnetization, so the recovered inversion models show some bias. The distribution of the recovered model (Figure 4g and 4h) from $Ba$ is closest to the true model, and the recovered model from NSS also exhibits good performance (Figure 4c and 4d). The inverted results of AS show a large anomalous source distribution range (Figure 4e and 4f), but the location is correct. Furthermore, all of the inverted results agree with the sensitivity analysis shown in Figure 2.

Essentially, the recovered models (Figure 4) are consistent with the true source, although the maximum values of the models are much higher than the true susceptibility values (1 SI) (Figure 4a–4h). This is considered to be due to the strong remanence that influences the amplitude data. As is common in most magnetic inversions, the susceptibility values of the inverted model on the fuzzy boundaries are much smaller than those of the actual susceptibility (Figure 4a–4h). In the second numerical example, we invert the data simulated for a 3D dipping slab shown in Figure 5a and 5b. A northerly directed inducing field ($D=0\xb0$, $I=75\xb0$) and an easterly directed remanent magnetization ($D=90\xb0$, $I=\u221215\xb0$) were used. The Konigsberg ratio ($Q=|J\u21c0rem|/|J\u21c0ind|$) of the remanent magnetization for this model is two (Figure 5a). Figure 6a–6d shows that, through the amplitude transformations, the influence of the inclination and remanence has been somewhat reduced.

First, the data were directly inverted, assuming that the directional information of the effective magnetization is known. Figure 7a and 7b illustrates that the inverted model shows significant deviations from the spatial distribution of the true model. The inversion is severely affected by the strong remanence component, which dramatically changed the anomalous distribution of the magnetic field. The data cannot derive a geologically meaningful inversion with an effective magnetization direction. The presence of remanence may complicate the magnetic inversion (Morris et al., 2007). To demonstrate this, the direction of remanent magnetization from Figure 8 was altered to $D=60\xb0$ and $I=15\xb0$. The inverted model is shown in Figure 8c and 8d.

Numerically, the data misfit of the inverted model (Figure 8) was reduced (8.71%, 30 iterations) compared to the data misfit of the inversion in Figure 7 (13.88%, 30 iterations). Although a compact and reasonable model is recovered (Figure 8) with effective magnetization, significant errors remain. The experiments also illustrate that it is impossible to reach an ideal model that can provide a good fit with anomalous magnetic data containing strong remanent magnetization. Because the relationship between the remanent magnetization and the magnetic inversion is beyond the scope of this study, it is not discussed further.

Then, the amplitude data related with the model in Figure 5 were inverted ($D=0\xb0$, $I=75\xb0$ for the inducing field and $D=90\xb0$, $I=\u221215\xb0$ for the remanent magnetization), and it was assumed that the direction of effective magnetization was unknown in all of the amplitude inversions. As Figures 9 and 10 show, the recovered models of NSS and AS are not significantly affected by the remanent component. Although the recovered models do not clearly indicate the dipping angle of the slab, they successfully reveal the correct location of the true model.

The recovered $Ba$ model agrees with the dip direction (Figure 11a and 11b), depth, and location of the true model, and it is closest to the true model in this amplitude inversion test.

As Figures 9–11 show, similar to the common magnetic inversion, the recovered models have a broader extent of lower susceptibility than the true model. There are significant deviations among the recovered results of NSS, AS, and $Ba$ in the model spatial distribution. The inversion of AS represents the excessive structure distribution across the depth ranges (Figure 10). The $Ba$ inverted model is the most compact (Figure 11), and it delineates the critical features of the actual model, such as the depth and dip direction. Although the inverted model of NSS does not indicate the dip direction of the slab, it distributes at a reasonable depth and position. In summary, these experiments illustrate that the amplitude data $Ba$ are efficient for recovering a source with a strong remanence.

Finally, we compare the convergence of the GN and NLCG algorithms using these two numerical examples (Figure 12). The GN and NLCG converge effectively in 100 s. The number of iterations until convergence and the time for each iteration differ for NLCG and GN. In these two synthetic experiments, the inversion of the NLCG algorithm shows a faster rate of data misfit minimization than that of the GN algorithm. Moreover, by comparing the convergence rates of the two synthetic model inversions (Figure 12), we conclude that additional iterations and time are needed for the recovery of a spatially complex source with high susceptibility. Overall, NLCG requires less time to converge to a lower data misfit than GN in both inversion examples.

## VALIDATION USING FIELD DATA

High-resolution magnetic data were collected from the Pea Ridge iron oxide-apatite-rare earth element (IOA ± REE) deposit by the U.S. Geology Survey in southeast Missouri, USA (Figure 13). This deposit is one of the best-characterized large Fe oxide deposits. The Pea Ridge deposit is hosted by volcanic rocks of early Mesoproterozoic age (approximately 1500–1440 Ma) (Aleinikoff et al., 2016; Day et al., 2016) and is covered by thick Paleozoic sedimentary rocks. The magnetite and hematite orebody have apparent high magnetic susceptibility ($0.2<\kappa <0.7\u2009\u2009SI$). According to Harlov et al. (2016) and McCafferty et al. (2016), the host rock is a high-silica rhyolite with strong remanence ($Q\u22481$). Allingham (1976) reports that the direction of the remanent vectors follows a random distribution. The geology of Pea Ridge has been fully depicted by underground mining and mapping expeditions (Seeger, 2001). For these reasons, this data set was used as field data for validating our proposed inversion algorithm.

The data were acquired by an airborne Cesium stinger-mounted magnetometer at a nominal height of 80 m in 2014. The survey contained 90 north–south flight lines spaced every 400 m. We selected 8509 observation points from the data over an area of $8.0\xd76.4\u2009\u2009km$. The geomagnetic flux of the study region is 53,036 nT, and the inclination and declination are 66° and 359.6°, respectively. We calculated the quantities of $Ba$, AS, and NSS from the total anomaly data (Figures 14 and 15). The positions of two cross sections encompassing the data, A-A′ and B-B′, are marked in Figures 14b and 15.

Although the orebody is distributed across the depth of 400–815 m, which is covered by 11 geologic maps (Seeger et al., 2001), McCafferty et al. (2016) propose that the magnetic orebody of Pea Ridge is vertically extensive and traceable to depths of approximately 2 km. Therefore, we set the inversion model region depth to 2400 m. In total, 5,148,000 cells were involved in the inversion calculation (220 cells east–west × 180 cells north–south $\xd7$ 130 cells top-bottom), and the mesh for the active model region contains 1,920,000 cells (200 cells east–west $\xd7$ 160 cells north–south $\xd7$ 60 cells surface-bottom). Although the nonorthogonal grid can accelerate the inversion procedure and works well with the relatively small-scale grid used in our experiment, we found that the proposed algorithm is more stable with an orthogonal grid for solving large-scale problem. Therefore, in this field data example, we used uniform 40 m cubic cells to form an orthogonal grid for inversion.

The methods suggested by Hansen (1992) and Leão-Santos et al. (2015) were used to identify the regularization parameter $\beta $, as shown in Figure 16. In this regularization parameter identification experiment (Figure 16), the curve is well-fitted to the monotonic functions related to the regularization parameter $\beta $. There is an obvious trade-off point near the regularization parameter value of one. The extensive inversion experiments show that an over-smoothing and less structure-inverted model will exhibit an inversion when the value of $\beta $ is greater than one. Conversely, a loose recovered model, with excessive structures, will be inverted with a relatively small $\beta $ value.

The effective susceptibility of the orebody was recovered using the proposed NLCG inversion algorithm with the amplitude quantities, as shown in Figures 17, 18, and 19. For these recovered models, the objective functions converge to approximately 5% after 20 iterations. Although these objective function misfits are not very small, they do provide reasonable estimations of the subsurface orebody property distributions that are consistent with the calculated quantities. The predicted data from the inversion are displayed in Figure 20, which is well-matched with the input amplitude data (Figures 14b and 15).

The experiments suggest that the grid design is important for the proposed algorithm. Insufficient padding cells or cells that are too large may lead to a highly incorrect solution. By using appropriate choices of the padding cell number and boundary surface sizes, the inversion should generate a reasonable inversion model with a rapid convergence rate. Based on the experiments, we suggest padding at least one-tenth of the number of active cells in each direction outside of the grid as the inactive free space cells to increase the accuracy of the approximated boundary conditions.

According to McCafferty et al. (2016), 143 apparent magnetic susceptibility samples of the iron orebody have been measured; these came from several drill cores mined from the Pea Ridge deposit (Figure 21). The average apparent magnetic susceptibility of the orebody (magnetite and hematite) is approximately 0.2–0.5 SI. The abundant magnetite ore has a significantly higher apparent magnetic susceptibility compared to the host rock. Based on this information and considering that the common recovered effective susceptibility estimates are lower than the actual apparent magnetic susceptibility measurements, the cutoff display value for all inverted models was set as 0.25 SI in this field data experiment.

According to the vertical geologic cross sections of different mine levels (400–815 m) at the Pea Ridge iron mine (Seeger et al., 2001), the orebody is a single dipping body (Figure 22). We compare our inverted results along the cross sections shown in Figure 14b (Figure 22).

As Figure 22a and 22b shows, a hematite zone occurs at the top of the deposit. Although pure hematite is essentially nonmagnetic, McCafferty et al. (2016) conclude that three discovered iron oxide deposits (Pilot Knob, Iron Mountain, and Pea Ridge) in southeast Missouri, regardless of the iron mineralization style, have a high apparent magnetic susceptibility. The samples of the hematite zone at Pea Ridge also have high apparent magnetic susceptibility (Figure 21). This indicates that this shallow hematite-rich zone surrounds all of the inferred magnetite-rich bodies and significant magnetite alteration exists in the hematite zone. Harlov et al. (2016) reveal that this hematite zone forms a gradational contact with the magnetite ore and contains patches of magnetite. This also provides another interpretation for the high magnetic susceptibility of the shallow hematite zone in Pea Ridge. In addition, primary, low-volume, REE-rich breccia pipes occur along the boundary of the magnetite ore, and cannot be distinguished geophysically (Seeger, 1992).

Here, we compare the inversion results (Figures 17–19) with these geologic interpretations. First, all inversion models (Figures 17–19) show a compact magnetic body, and the value of magnetic susceptibility agrees with prior information (0.2–0.5 SI). We note that the inversion results of $Ba$ (Figure 17) show the correct dip compared to the vertical geologic cross-section map (Figure 22). Almost all inversion models show that the depth of the top of the dipping orebody is approximately 400 m at the correct horizontal position. However, the recovered models of $Ba$ and NSS indicate that the depth of the orebody bottom is 2150 m (Figures 17b and 19b), whereas the inversion result of AS shows that the orebody extends to depths beyond the inversion region. Considering that the recovered models of AS also have a broad susceptibility distribution and extents with depth in our synthetic experiment, it is likely that the inverted results, which estimate that the orebody is located at a depth between 400 and 2150 m, are correct. These inverted models (Figures 17b and 19b) partly validate the Seeger (1992) hypothesis, suggesting a greater volume of the magnetic orebody at depth. We also compare the geologic cross sections provided by Seeger et al. (2001) (Figure 23a), the estimated strike directions recovered models of $Ba$, AS, and NSS (Figure 22b–22d) exhibit the correct strike N60°E at a depth of 510 m (Figure 23a). As our inverted model shows (Figure 17b), the deepest mine level reached during exploration (815 m) is a relatively shallow part of the deposit.

The positions of the vertical geologic cross sections A-A′ and B-B′ are marked in Figures 14 and 15. In the inversion results (Figure 24), only $Ba$ shows a dipping angle that is consistent with the magnetite ore described by the geologic cross section maps.

Next, we compare the convergence of algorithms NLCG and GN in this field data inversion, as shown in Figure 25. NLCG finishes 25 iterations and converges in 6 h, whereas GN converges in 12 h and finishes 25 iterations in 34 h. Obviously, the number of iterations and convergence times differs between NLCG and GN algorithms. NLCG requires more steps but is faster at each step, whereas GN requires fewer steps but takes a long time for each step. In addition, NLCG reaches a larger initial reduction than GN. NLCG is more effective for this field data inversion example.

The inverted model of the NLCG and GN algorithms was compared in this field data inversion, as shown in Figure 26. The inverted model of NLCG (Figures 17, 24a, and 24b) and the inverted result of GN also reveal a dipping model that is consistent with the previously obtained geologic information (Figure 22). However, there are some locations where the result differs from that of the inverted model of NLCG. First, as Figures 17 and 22a show, the orebody has a small syncline structure at a depth of approximately 400–700 m. This detailed structure is not present in the inverted model of GN. Second, as shown in Figure 26b, the inverted model of GN only extends to a depth of 1500 m, which is slightly shallower than the previous geologic interpretation (2000 m). However, there are large volumes of anomalous sources that have relatively low susceptibility ($0.1<\kappa <0.25$) below the depth of 1500 m (Figure 26c and 26d). The inverted model of GN has more fuzzy boundaries and a broader extent than the NLCG inverted model at depth. In general, the shapes of the GN and NLCG inverted models are very similar. In summary, it is reasonable to state that the proposed algorithm can be used to efficiently invert a high-susceptibility magnetic source exhibiting strong remanence.

## CONCLUSION

We solve the PDE solution of magnetic data amplitude using a conjugate gradient algorithm. We tested three primary amplitude transformation methods of anomalous magnetic data using inversion tests. The numerical and field data results indicated that the amplitude of the anomalous magnetic field is an efficient amplitude transformation that can be used in the PDE inversion algorithm. The proposed NLCG algorithm showed a fast convergence rate regarding the computational time required to compute accurate solutions in our experiments. Using an FVD to solve the PDE problem is a discretized integral weak form method, in which the error comes from finite discretization. One of the most efficient methods to reduce this error component is to use a fine mesh. It is important to consider the speed of convergence of an inversion algorithm with a fine mesh, and this will be helpful to improve the accuracy of the inversion from other points of view.

The results also showed that the NLCG algorithm is a feasible approach to the large-scale inversion problem. Possible extensions to the method involve using preconditioning to enhance the performance of the algorithm. This case study demonstrated the efficacy of using amplitude inversion within the PDE framework in the presence of strong self-demagnetization and remanent magnetization. The proposed algorithm does not depend on site-specific prior geologic information, and it is generally applicable for data interpretation in other areas with analogous exploration targets and challenges.

## ACKNOWLEDGMENTS

This work was supported by grants from the National Natural Science Foundation of China (grant nos. 41674110 and 41630317). The authors thank the assistant editor J. Shragge and another anonymous assistant editor, the reviewers M. Leão-Santos, B. Morris, and another anonymous reviewer for their constructive comments and suggestions that significantly improved the manuscript. We also thank M. A. Kass for his valuable suggestion.

## DATA AND MATERIALS AVAILABILITY

Data associated with this research are confidential and cannot be released.