Mineral exploration under a thick sedimentary cover naturally relies on geophysical methods. We have used high-resolution airborne magnetic and gravity gradient data over northeast Iowa to characterize the geology of the concealed Precambrian rocks and evaluate the prospectivity of mineral deposits. Previous researchers have interpreted the magnetic and gravity gradient data in the form of a 2D geologic map of the Precambrian basement rocks, which provides important geophysical constraints on the geologic history and mineral potentials over the Decorah area located in the northeast of Iowa. However, their interpretations are based on 2D data maps and are limited to the two horizontal dimensions. To fully tap into the rich information contained in the high-resolution airborne geophysical data, and to further our understanding of the undercover geology, we have performed separate and joint inversions of magnetic and gravity gradient data to obtain 3D density contrast models and 3D susceptibility models, based on which we carried out geology differentiation. Based on separately inverted physical property values, we have identified 10 geologic units and their spatial distributions in 3D which are all summarized in a 3D quasi-geology model. The extension of 2D geologic interpretation to 3D allows for the discovery of four previously unidentified geologic units, a more detailed classification of the Yavapai country rock, and the identification of the highly anomalous core of the mafic intrusions. Joint inversion allows for the classification of a few geologic units further into several subclasses. We have demonstrated the added value of the construction of a 3D quasi-geology model based on 3D separate and joint inversions.


The Midcontinent Rift system (MRS) and the surrounding Precambrian rocks are known to host highly significant mineral resources. One such example is the copper-nickel-platinum group element deposits in the Lake Superior region (Miller and Nicholson, 2013) where the intrusions related to the MRS, such as the Duluth Complex (Miller et al., 2002), crop out. However, more than 90% of the region affected by the MRS is covered by Phanerozoic sedimentary rocks and sediments with very limited drillholes, such as the Decorah Complex region. A direct consequence is that much of the geology of the Precambrian basement under the Phanerozoic cover is poorly understood and the mineral resource potential is largely unknown. Geophysical methods are useful for mapping the lithology and structure of concealed rocks under cover because they would produce geophysical responses that can be readily measured on or above the surface of the earth, if lateral or vertical variations in physical properties exist in and among these rocks.

To evaluate the mineral resource potentials over the northeast Iowa region, high-resolution aeromagnetic and airborne gravity gradient (AGG) data were acquired. Figure 1 displays the Gzz component of the AGG data after terrain correction and the total magnetic intensity (TMI) data over the Decorah area. Large-amplitude gravity and magnetic highs are observed and interpreted to result from a buried intrusive complex of mafic-ultramafic rocks (Drenth et al., 2015). Drenth et al. (2015) interpret the measured AGG and magnetic data to construct an interpreted geologic map of the Precambrian rocks under the sedimentary cover. They mainly identify two groups. One group consists of older units including the host rocks, a northern mafic intrusion, and a silicic pluton. The other group, with younger rocks, comprises the southern mafic intrusion (the Decorah Complex), mafic-ultramafic and intermediate/silicic rocks that intrude the host rocks, and a few northeast-trending dikes. However, their interpretation is largely based on 2D data maps.

The motivation for our work is to extend the 2D interpretations in Drenth et al. (2015) into 3D by creating a 3D quasi-geology model. We began with independent 3D inversions of gravity gradient and magnetic data and reconstructed a 3D density contrast model and a 3D susceptibility model. The estimated physical property values were visualized in a scatterplot. Next, we classified the physical property values in the scatterplot into 10 different classes (or clusters) with each class corresponding to a distinct geologic unit. The classification was achieved based upon the natural grouping patterns observed in the physical property values (including the trends, grouping, etc) and the expected ranges of physical property values for each geologic unit. The 10 distinct geologic units identified on the scatterplot were subsequently mapped into 3D space. Following the work of Li et al. (2019), we refer to the differentiated units in 3D as the 3D quasi-geology model, and the above-mentioned process of arriving at such a model as geology differentiation henceforth. The 3D quasi-geology model is an approximation of the geology in subsurface based on the integration of physical properties and subject to the resolution of the geophysical data. Our 3D quasi-geology model shows that the northern intrusion and Decorah Complex have a core of higher density and susceptibility that could be associated with mineralization. Additionally, the undefined mafic intrusive rocks have similar physical properties to the northern intrusion and Decorah Complex and should be investigated further when more data are available. Our geology differentiation work also shows a dichotomy in the host rocks with two distinct domains and a partitioning of the central silicic intrusion into two parts. The 3D quasi-geology model complements the 2D geologic map created by Drenth et al. (2015) and allows for a better understanding of the structures and compositions of the igneous intrusions in all three dimensions.

We also carried out joint cross-gradient inversion (Gallardo and Meju, 2003, 2004) and obtained a 3D density contrast model and a 3D susceptibility model that are structurally similar. The same methodology for geology differentiation described above was then applied to the jointly inverted density and susceptibility values. We have found that the enhanced structural similarity between density and susceptibility models allows for finer classification of classes 1, 3, 4, and 10, previously identified based on separate inversion results, into a few subclasses. Each subclass is characterized by a distinct density-susceptibility relationship that was previously obscured by the strong scattering among the separately inverted values, but revealed by joint inversion. The results show that the joint inversion can individualize specific bodies or its parts and, therefore, better detail the geologic units undercover. The new discovery based on jointly inverted models includes the segmentation of the core of the northern intrusion and Decorah Complex into several parts with different susceptibility ranges, which can be an important guideline for planning future drilling.

Below, we begin with a summary of the geologic setting of our study area, followed by a description of geophysical data collection and processing. We also detail the procedures that we have developed to estimate the data errors, a tricky yet important task for practical inversions of field data. Then, we describe the inversion methods and report the tuning parameters that we used to convert the measured data into 3D physical property models. The focus of our work is on geology differentiation, in which we differentiate and characterize 10 geologic units in 3D based on the independently and jointly inverted 3D models and prior interpretations of Drenth et al. (2015).

Geologic setting

The study area focuses on magnetic and gravity anomalies located in the northeast Iowa and southeast Minnesota border. Following Drenth et al. (2015), we will refer to this region as the Northeast Iowa Intrusive Complex (NEIIC) henceforth. Although the age of the NEIIC is still unproven, the complex is located near the MRS, which extends across the central United States and is responsible for the formation of extensive volumes of magma, including the Duluth Complex, which hosts large deposits of copper-nickel-platinum group elements (Miller et al., 2002). The mafic intrusions of the Duluth Complex are associated with gravity and magnetic anomalies (Chandler, 1990; Miller et al., 2002). Similarly, significant magnetic and gravity anomalies are also observed in the NEIIC; however, the full potential of the area is still underexplored because of the Phanerozoic sediments (200–700 m) that cover the region (Drenth et al., 2015).

In the NEIIC region, Drenth et al. (2015) interpret magnetic and gravity gradient data of the Decorah Complex and integrate its geophysical signature with information from two boreholes and the regional geology. The authors present an interpreted geologic map of the crystalline Precambrian rocks (Figure 2), which shows five main features: (1) a southern mafic intrusion (the Decorah Complex), (2) a central silicic intrusion, (3) a northern mafic intrusion, (4) the host rocks, which are intruded by small (5) mafic and intermediate or silicic bodies. Drenth et al. (2015) further detail the above five main groups into 12 geologic units, which belong to two age periods: older units that show signs of deformation (1.8–1.72 Ga rocks of the Yavapai province) and younger units that are largely undeformed (approximately 1.1 Ga and possibly Keweenawan). In the group of older units, they detailed the host rocks into undifferentiated metavolcanics, plutons, and metasediments (unit Y), the central silicic intrusion into a possibly S-type silicic pluton (unit Ysp), and the northern mafic intrusion (Van Schmus et al., 1989) into undifferentiated mafic rocks (unit Ym), a metagabbro (unit Ymg), a weakly magnetized part of the intrusion (unit Ywm), and a strongly magnetized part of the intrusion (unit Ysm). In the group of younger units, they detail the southern mafic intrusion into the gabbro of Decorah Complex (dg) and the central weakly magnetized rocks of the Decorah Complex (dwm). Also in the younger units, there are the mafic intrusive (unit mi), the intermediate or silicic intrusions (unit ii), and dikes. Only the metagabbro (unit Ymg) and the gabbro of the Decorah Complex (dg) were confirmed by drilling.

Although the lack of drilling does not allow confirmation of most of the interpreted units, the integration and interpretation performed by Drenth et al. (2015) permit a better understanding of the geologic setting under cover. However, this first step to unveil the geology under cover is still limited to 2D data interpretation. For this reason, we further explore the magnetic and gravity gradient data to build physical property models in 3D aiming to understand the distribution of the anomalous bodies at depth.

Geophysical data

From December 2012 to January 2013, a high-resolution airborne magnetic, gravity gradient, and time-domain electromagnetic survey was conducted over part of the NEIIC region. All airborne geophysical data are publicly accessible at http://mrdata.usgs.gov/. For the study reported in this manuscript, we used only the magnetic and gravity gradient data. A total of 94 east–west traverse lines were flown with 400 m spacing and a nominal clearance of 80–100 m over largely flat terrain, in addition to nine tie lines with 4000 m spacing.

The full gradient tensor data are measured using a full-tensor-gradiometry (FTG) system onboard an airplane. The FTG system contains three gravity gradient instruments (GGIs), each consisting of two opposing pairs of accelerometers on a rotating disc. The acquired data underwent a series of proprietary standard processing steps performed by the contractor. Figure 1a shows the terrain-corrected Gzz component of the measured AGG data. A density value of 2.4  g/cm3 found to locally minimize the terrain effects (Drenth et al., 2015) and is therefore used in our study. Regarding the use of full tensor components versus a single component in the inversion, Paoletti et al. (2016) find that, when the algebraic ambiguity (the difference between the number of unknowns and the number of data) and signal-to-noise ratio (S/N) are comparable, the overall quality of the inversion results does not depend upon whether multiple tensor components or only a single tensor component are used in the inversion. Pilkington (2012) shows that, at larger measurement-source distances, all tensor components and combinations provide similar information content concerning the source density distribution; therefore, including more components to the inversion does not provide any advantages. Based on the work of Pilkington (2012) and Paoletti et al. (2016), we argue that whether including more components in the inversion would result in better density models should be assessed on a case-by-case basis by taking algebraic ambiguity, S/Ns, and measurement-source distances all into account. We also argue that, if the S/Ns of the other components are significantly higher than that of Gzz, it is very likely that these components would not contribute any new information to the construction of a density model. Worse still, including these components into an inversion could deteriorate the inversion results because the inversion needs to fit not only the high signal-to-noise components but also the low signal-to-noise components. In our study, we chose to use only Gzz data instead of the full tensor components for two reasons. First, a thorough investigation of whether more tensor components would result in a better density model for our study area is nontrivial and deserves a full paper itself. Second, given the very limited geologic and geophysical work performed in our study area as well as the lack of any ground-truth data to compare with, it is nearly impossible to compare the density model obtained from multiple tensor components and that from only Gzz, and to tell which one is better.

The magnetic data were collected using a cesium vapor magnetometer. A lag correction was applied to correct for the distance between the magnetometer and the GPS antennas, followed by removal of the earth’s magnetic field and leveling. The ambient magnetic field has an inclination of 70 and a declination of zero. Figure 1b summarizes the processed magnetic data.

Data uncertainty estimation

Geophysical measurements always come with noise. It is important to understand the noise level before any inversion is carried out. Otherwise, we risk letting noise manifest itself in the recovered models without realizing its presence. Data-error estimates (or equivalently, uncertainties) are a required component of many inversion algorithms. However, for the data in this survey, data uncertainties were not provided and, therefore, must be estimated. We took two heuristic approaches to estimating data errors. In the first approach, we intended to use the L-curve method (Hansen and O’Leary, 1993) as an error estimator (Kindermann and Raik, 2019). First, we assigned 5% of the absolute values of the measurements (plus a small positive threshold value to prevent the occurrence of zero-error estimates) as an initial guess. Because it is unlikely that the random initial guess is a reasonable estimate of the data errors, we then adjusted the initial guess by constructing a Tikhonov curve. We achieved this by constructing a series of equivalent source layers that fit the measurements to different degrees (by using different regularization parameters spanning several orders of magnitude). To speed up the computations, we decimated the measured data along the flight lines by a factor of three and end up with a total of 17,548 measurements for the Gzz data and magnetic data. For each equivalent source, we calculated the data misfit value and the regularization term value, based on which we constructed two Tikhonov curves, one for the Gzz inversion (3a) and the other for the magnetic inversion (3b). We expected a well-defined “elbow” on each of the Tikhonov curves because the data-misfit value corresponding to the elbow would then represent a reasonable target data misfit. However, when plotted on the log-log scale as shown in Figure 3, the Tikhonov curves do not show such elbows and, consequently, are not amendable to being used with the L-curve criterion.

Then, we took a second approach by manually examining the spatial patterns of data residuals resulting from the previous equivalent source layers constructed using vastly different regularization parameters. When the regularization parameter is large, the data residuals show coherent and strong spatial patterns, indicating that the measured data were underfitted. We visually inspected the data residuals resulting progressively from larger regularization parameters to smaller ones. As the regularization parameters become smaller and smaller, the corresponding data residual maps start to show fewer and fewer coherent features. At some point, the data-residual maps are dominated by mostly random patterns. We then used the corresponding data misfit values as the appropriate data misfit values, based on which we adjusted our initial estimates of data errors. Our final error estimates for Gzz and TMI data are 17.4% and 0.25%, respectively, of the absolute values of the measured data.

Geophysical inversions

To make informed inferences about the structure and composition of the Precambrian basement, we performed separate inversions of the airborne geophysical data displayed in Figure 1. We define separate inversions as inversions that are carried out independently, and without any exchange of information between inversions. For our study, we performed two inversions, one inversion of gravity gradient Gzz data and a second inversion of magnetic data. These two inversions were performed separately and independently of each other, using two sets of inversion codes. This is in contrast to joint inversion in which only one inversion is performed to simultaneously recover at least two different physical property models. During the joint inversion, the multiple physical property models are constantly exchanging structural and/or physical property information, and information from one physical property is used to constrain the construction of another physical property model and vice versa. In our study, we implement the cross-gradient joint inversion (Gallardo and Meju, 2003, 2004) to seek a density and a susceptibility model that are structurally similar.

Below, we describe the inversion methodologies and report the adjustable parameters that we used when implementing 3D inversions.

Separate inversions

We used the standard smoothness-based Tikhonov regularized inversion method, as described in Li and Oldenburg (1996, 1998), to perform separate inversions of Gzz and magnetic data. With this method, the goal is to seek one model that can adequately reproduce the geophysical measurements while at the same time having the maximum degree of smoothness in all three spatial directions. Mathematically, an inversion is carried out by minimizing the following objective function:  
ϕ(m)=Wd(dobsGm)||22+β(αsWs(mmref)22+αxWx(mmref)22+αyWy(mmref)22+αzWz(mmref)22),s.t.  mminmmmax,
where m is the physical property model to be recovered; mmin and mmax are, respectively, the lower and upper bounds for the recovered physical property values; Wd is a data weight matrix consisting of data uncertainties and their correlations; dobsis the geophysical measurements to be inverted; G is a symbolic representation of the forward modeling operator that maps any given physical property model to its corresponding geophysical response; Wi (i=s,x,y, and z) are discretized spatial weighting matrices; and mref is a reference model that summarizes any prior knowledge or assumptions about the geology. In our study, we assume no prior knowledge and use a zero reference, that is, mref=0. We also incorporate the sensitivity-based spatial weighting (Li and Oldenburg, 2000) into Wi (i=s,xy, and z) matrices to account for the decay of potential field data with distance (Li and Oldenburg, 1996, 1998). The parameter β is commonly termed the regularization parameter in the literature. It determines the goodness of fit between observed and predicted data and the structural complexity of the final recovered model. In our work, the value of the regularization parameters is determined by a target data misfit value. The weighting parameters αi (i=s,x,y,, and z) are directly related to the correlation lengths in the northing, easting, and depth directions.

For inversion of Gzz data, we set unrealistically large lower and upper bounds, 10 and 10  g/cm3, which caused the bound constraints to have very little effect on the recovered density model. This was done because we expect positive and negative density contrasts to explain the features observed in Figure 1a. We used a regularization parameter of 120. The subsurface area of interest was divided into 720,000 rectangular prisms. The cells in the core region are 250 × 250 × 250 m3. We also added padding cells with progressively increasing widths outside the core region. The entire inversion domain spans 35 km in the west–east direction, 37.5 km in the north–south direction, and 10 km in depth. We set the correlation length in the northing, easting, and depth directions to five times, five times, and three times, respectively, the cell widths in each direction. The topography is fully accounted for in the Gzz and magnetic inversions.

For inversion of magnetic data, we used the same lower and upper bounds, the same discretization strategy (i.e., the same mesh), and the same correlation lengths to maximize the compatibility of the resulting susceptibility and density models. This is especially important for the subsequent joint interpretation of the inverted models because we would like any notable differences between the two models to result from geology as opposed to the different setup parameters used in the two separate inversions. The regularization parameter we used for magnetic inversion is 300. We note that we took remanence into account by allowing the magnetization directions to be in either the direction of the present-day earth’s magnetic field or the opposite. We acknowledge the potential limitation of this approach, but we also point out that we had no trouble fitting the observed magnetic data with this strategy. Indeed, we can fit the magnetic observations to any desired degree, as evidenced by the wide range of data misfit values that we achieved in Figure 3b. This is a strong indication that introducing more possibilities of magnetization directions is not required by the magnetic measurements. However, we emphasize that this strategy might not work well in other intrusive complexes, depending on the age and nature of the underlying rocks’ formation and the subsequent geologic history unique to the area of study under investigation. In such areas, full-3D magnetization inversions (Lelièvre and Oldenburg, 2009; Fournier, 2015; Li and Sun, 2016; Fournier, 2019) might be necessary.

The separately inverted 3D density and susceptibility models as well as their respective predicted data are summarized in Figure 4.

Joint inversion

The objective function in equation 1 can be extended to joint inversion of two sets of geophysical data, as shown below:  
where β1 and β2 are the two regularization parameters, one for each data set, and λ is the weight on the cross-gradient term. The last term measures the structural differences between a density and a susceptibility model, and it also determines how structural information from one model is used to constrain the other model. The matrix Wm1 results from a consolidation of the spatial weighting matrices Wi (i=s,x,y, and z) and the associated weighting parameters αi (i=s,x,y, and z) previously described for separate inversions. The same is true for the matrix Wm2. In our implementation, we set β1 to 108, β2 to 270, and λ to 3e+15. The rest of the parameters are the same as those described above for the separate inversions.

Geology differentiation

The inverted density and susceptibility models show the behavior at depth of individual bodies and structures in the study area. However, the physical property models are only proxies of geology, and as pointed out recently by Li et al. (2019), an extra step of integrating the inverted physical property values into a quasi-geology model is needed. The quasi-geology model (Kowalczyk et al., 2010; Martinez and Li, 2015; Sun and Li, 2015; Devriese et al., 2017; Fournier et al., 2017; Kang et al., 2017; Melo et al., 2017; Li et al., 2019) shows the spatial distributions of various geologic units that can be identified by geophysical methods, and it allows geologists to develop useful insights that guide further tasks such as drilling for mineral exploration and maximize the value of geophysical inversions.

As the first step to geology differentiation, we integrated the inverted density contrast and susceptibility values by constructing a scatterplot of both properties, in which each point in the plot corresponds to one cell in the model. Figure 5a and 5b shows the crossplot of separately and jointly inverted density and susceptibility values, respectively. Both crossplots summarize the variations of the estimated physical property values, and they exhibit some interesting trends and natural grouping. We observe that, to the first order, the separately inverted values in Figure 5a and the jointly inverted values in 5b are distributed in a similar way. But the former exhibits more scattering, whereas the latter displays several well-defined linear features that stick out in the upper-right portion of the crossplot. We note that such linear features have also been reported by previous researchers (Fregoso and Gallardo, 2009; Doetsch et al., 2010; Infante et al., 2010) and are used for differentiating geologic units (Linde et al., 2006, 2008).

Our next step is to classify the inverted values in the scatterplots into different classes (or, groups), with each class characterized by a distinct range of density and susceptibility values. Hereafter, for the convenience of explanation, we refer to a scatterplot of the inverted physical property values as an actual scatterplot to be contrasted with the conceptual scatterplot introduced below.

To differentiate the inverted values into different groups, we need to first determine the expected ranges for all geologic units. We did this for the 10 geologic units identified in Drenth et al. (2015). Here, we exclude the normally and reversely magnetized diabase dikes from this step because the identification of these dikes were based on magnetic lineaments and we do not expect 3D smoothness-based inversions to recover these dikes. The expected physical property value ranges are determined by correlating the inverted physical property anomalies with the characteristics of the geophysical data anomalies that Drenth et al. (2015) identify for each unit. In our work, this was done by visually examining and linking the inverted density and susceptibility anomalous features with the different units described in Drenth et al. (2015). We acknowledge the subjectivity involved in this step and note that automated tools are still lacking and machine learning might be useful to automate the process. The final results of the expected physical property value ranges for different units are summarized in a conceptual scatterplot (Figure 6), which contrasts with the actual scatterplot shown in Figure 5a. The conceptual scatterplot is critical for geology differentiation because it guides our efforts of differentiating (or grouping) the actually obtained physical property values in Figure 5a into different geologic units. We point out that the final grouping in an actual scatterplot such as those in Figure 5 might very well be different from the expected ranges of physical property values in a conceptual scatterplot. Such differences represent adjustments made to the expected behavior after inversions and new information contents from inversions.

In the conceptual scatterplot, shown in Figure 6, we were able to identify seven geologic units:

  1. 1)

    one unit comprising the metagabbro (Ymg) and the gabbro of the Decorah Complex (dg) with the highest density contrast and susceptibility values

  2. 2)

    the mafic intrusive (mi) with moderate to high density and susceptibility

  3. 3)

    the strongly magnetized part of the layered intrusion (Ysm) with moderate density contrast and high susceptibility

  4. 4)

    one unit comprising the undifferentiated mafic (Ym), the weakly magnetized part of the layered intrusion (Ywm), and the weakly magnetized part of the Decorah Complex (dwm) with moderate density contrast and susceptibility values

  5. 5)

    the metavolcanics, plutons, and metassediments (Y) with low density contrast and low susceptibility

  6. 6)

    the intermediate or silicic intrusive (ii) with very low density contrast and moderate susceptibility

  7. 7)

    the silicic pluton (Ysp) with the largest negative density contrast and the lowest susceptibility values.

In the next step, we applied the conceptualization summarized in Figure 6 to the inverted physical properties scatterplots (Figure 5) with the objective of identifying the same units by empirically testing the ranges for each unit. The conceptual plot guided the identification of some of the main units. Although it did not allow the identification of all of them, it helped reveal some units that were not previously described, contributing to new understanding of the area at depth. Below, we first report our differentiation results based on separate inversions, followed by results based on joint inversion.

Geology differentiation based on separate inversions

Our geology differentiation results based on separate inversions are summarized in three forms: classification on the scatterplot in Figure 7a, classification on a 2D spatial map in Figure 8b, and a 3D view consisting of several depth slices in Figure 9. With the conceptualization in Figure 6, we successfully identified class 1 as the metagabbro (Ymg) of the northern intrusion and the core of the gabbro of Decorah Complex (dg). Further, differentiation between these two units is possible because three linear trends are visible in class 1. However, they are obscured by the strong scattering from the smoothness-based geophysical inversions. Therefore, we chose to group them into one class. But as is shown in the next section, joint inversion is be able to refine this level of detail and allows for further differentiation within class 1. Additionally, class 1 includes the central part of the undifferentiated mafic (Ym), of the northern intrusion, in the same class. The concept also has good correspondence between class 2 and the mafic intrusive (mi). This class also includes the outer part of the northern intrusion (Ysm, Ywm, and part of Ym), showing that unit Ysm has lower susceptibility and higher density contrast and units Ym and Ywm have higher susceptibility and density contrast than initially predicted. In addition, class 2 includes the outer part of the Decorah Complex (dg), suggesting that a layering or chemical zonation of the mafic intrusion might be present.

Class 3 differentiates a specific part of the Decorah Complex (dg), showing that the southwestern part of the intrusion has higher susceptibility (the white ellipse in the susceptibility model in Figure 9b). Class 4 shows a specific mafic intrusion in the east, one intrusion with lower susceptibility values than the other mafic intrusions. Most of the intrusion is characterized by a reversed magnetic polarity, as shown in Figure 7a. Classes 3 and 4 were not predicted by the conceptual plot.

In Drenth et al. (2015), the Yavapai province host rock (unit Y) is interpreted as having nondescript geophysical signatures or geophysical anomalies that do not resemble those over other interpreted geologic units. This unit, shown in the interpreted basement geologic map in Drenth et al. (2015) as a homogeneous background, occupies the whole study area wherever no other geologic unit is expected. But, our joint interpretation of 3D density and susceptibility models allows for the identification of two domains within this unit, classes 5 and 6, with the same susceptibility values, but class 5 shows reversed magnetic polarity. Classes 5 and 6 were correctly predicted by the conceptual plot (unit Y in Figure 6), but one body with even higher susceptibility (the pink ellipse in the susceptibility model in Figure 9c) was identified in the northwestern part of the Y unit, which corresponds to class 7 and is new in our interpretation.

Classes 8 and 9 correspond to the silicic pluton (Ysp), and according to the prediction, they have the lowest density and susceptibility values. However, the new information from our interpretation is that class 8 corresponds to a deep eastern extension of the pluton with relatively high susceptibility (but the negative sign shows that it has reversed magnetic polarity), whereas class 9 shows a part of the pluton with lower susceptibility (also with reversed magnetic polarity) than class 8. Class 10 was successfully predicted by the conceptual plot, but it only shows two bodies of the intermediate or silicic intrusion (unit ii). One is in the north and shows good spatial correspondence with the interpreted map of Drenth et al. (2015), and the other was previously interpreted as host rock (Y). Interestingly, these two bodies are elongated trending to the northeast direction, the same trend of the dikes and faults in the region.

Our interpretation based on 3D inversion shows that the weakly magnetized part of the Decorah Complex (dwm) has a similar geophysical signature to the less magnetic host rock (class 5). Given no other data available to differentiate the former from the latter, we classified them into the same unit (class 5). Our 3D inversions also show that the center of the silicic pluton (Ysp) has a higher density contrast than its surroundings and similar susceptibility magnitudes but of an opposite sign. We classified this region into class 6. We also note that unit dg is represented as one single body by Drenth et al. (2015), but it is further detailed into three units (i.e., classes 1, 2, and 3) in our 3D-inversion-based geology differentiation results.

In our classification and interpretation, we attributed the recovered negative susceptibility values to a reversed magnetic polarity. To verify that these negative susceptibility values recovered from inversions are not artifacts from the smoothness regularization but are real geologic features, we set these negative susceptibilities to zero and used the resulting 3D nonnegative susceptibility model for a forward calculation of the magnetic total-field anomaly. We then calculated the difference between the measured magnetic data (shown in Figure 10a) and the simulated data due to a model without negative susceptibility values (shown in Figure 10b). We observe that the data-difference map contains a significant amount of features. These features cannot be explained by a model with all positive susceptibilities. We conclude that negative susceptibility values must be included into our 3D susceptibility model to explain all the features in the magnetic measurements in Figure 10a and, therefore, they are real features.

Geology differentiation based on joint inversion

We summarize our geology differentiation results based on joint inversion in a classified scatterplot in Figure 7b, a 2D spatial map in Figure 8c and a 3D volume-rendered visualization in Figure 11. In Figure 7b, we have classified the jointly inverted values into 10 classes, the same as the classification shown in Figure 7a, but with a few classes further divided into subclasses. For example, class 1 now consists of class 1a, 1b, 1c, and 1d, each of which corresponds to one well-defined linear feature in the scatterplot. Specifically, class 1a is characterized by higher susceptibilities among all four subclasses. It spatially corresponds to the border of the core of Decorah intrusion (dg) shown in Figure 8c. Class 1b shows moderate to high susceptibilities, and it manifests itself spatially in the core of the Decorah intrusion (in the south) and in the north Ymg and central Ym units. Class 1c has moderate susceptibilities, and it shows up in the northwest cores of unit dg. Class 1d has moderate-low susceptibilities located at the border of the northwest cores of dg. This zonation within class 1 can potentially represent petrological zonation of the mafic intrusion, but this can only be confirmed by drilling in the future.

Class 3, previously located in the southwest of the unit dg, now comprises two bodies: class 3a representing the core of a mafic intrusion in the northwest of the study area and class 3b in the south part of dg intrusion. Both units have the same susceptibility ranges, but class 3a has higher density contrast than class 3b. Class 4, a potential mafic intrusion (mi) in the east, was also divided into two subclasses, with class 4b having higher density contrast and slightly higher susceptibilities than class 4a. This finer classification might also be attributed to petrological zonation in the mafic intrusion. For class 10, the jointly inverted susceptibility values were used to separate the previously single body in the most north–northwest area into two bodies. Class 10b is located further to the north and has a higher susceptibility than class 10a. These classes correspond to the intermediate or silicic intrusive rocks (unit ii) interpreted by Drenth et al. (2015); however, they seem to be part of a different geologic unit.

Joint inversion allows for finer differentiation of the inverted values into more distinct geologic units. Due to the lack of drilling or any other ground-truth data, the differentiation results cannot be verified. But we point out that the identification of these subclasses as well as all of the other units can be used to guide future drilling activities and future geophysical data acquisition.


Data interpretation heavily depends on the interpreter’s experience and varies with different areas and available geophysical data. Nonetheless, it is a fundamental step for any decision making. For this reason, interpretation methods based on integration, which combines different yet complementary information from multiple data sets, are important to improve the quality of interpretation. The process becomes more challenging in undercover areas, such as our study area, where information from outcrops and boreholes is scarce to guide and confirm the interpretation. However, as mentioned above, interpretation helps to formulate hypotheses about the lithologies and structures underground and potentially guide exploration efforts with more specific information (as opposed to less specific information when such hypotheses are not available).

Geophysical data interpretation is performed as a continuum of building blocks, where the first stage is 2D data interpretation in maps. For the Decorah complex region, the work of Drenth et al. (2015), based on 2D geophysical data maps and limited geologic data, lays the foundation for an initial understanding of the area under the Phanerozoic cover. The second interpretation stage is the construction of separately and/or jointly inverted physical property models to explore the bodies and structures in 3D. In the third stage, geology differentiation is performed based on the expected ranges of physical property values for different geologic units and the natural groupings and features among the inverted physical properties. Finally, the geology differentiation results are summarized and visualized in a 3D quasi-geology model. We believe that geology differentiation and its subsequent visualization in a 3D quasi-geology model are a critical step toward getting closer to the “real” geology.

In summary, the geology differentiation method that we apply in this study identified the following:

  1. 1)

    For the northern intrusion, its central part is the same as what was identified in the borehole BO-1 (class 1b).

  2. 2)

    The outer border of the northern intrusion does not seem to have so many divisions as what was proposed by Drenth et al. (2015).

  3. 3)

    Classes 3, 4, and their subclasses show the possibility of differentiating specific mafic intrusions and their various parts in the study area.

  4. 4)

    The mafic intrusive (mi unit) has the same signature of the northern mafic intrusion and the gabbro of Decorah Complex in the geology differentiation. Therefore, it should be considered when evaluating the mineral potential of the region.

The segmentation of the scatterplots (Figure 7a and 7b) shows that density information allows for the differentiation of two main lithologic groups: the first, mafic intrusions (such as classes 1–4) and the second, silicic plutons, metasediments, and metavolcanic (such as classes 6–10). But the susceptibility information allows us to fine-tune different parts or bodies in the main classes. Examples include the differentiated parts of the mafic intrusion, the extension of the silicic pluton to the east, and the two parts of the host rock.


To develop a better understanding of the geology of the Precambrian basement rock over the Decorah region, we constructed a 3D quasi-geology model by performing 3D separate and joint inversions of gravity gradient and magnetic data, followed by a subsequent geology differentiation. A critical component of geology differentiation is to establish a conceptual scatterplot consisting of the expected physical property values ranges for all geologic units. We accomplished this by correlating the inverted density and susceptibility features with the geologic units identified by previous researchers. The established conceptual scatterplot was used to guide the differentiation of the inverted density and susceptibility values into 10 different units. The 10 distinct geologic units identified on the scatterplot were subsequently mapped into 3D space, resulting in a 3D quasi-geology model. We identified four new geologic units (classes 3, 4, 7, and 8), a dichotomy in the Yavapai host rocks with two distinct domains, a partitioning of the central silicic intrusion into two parts, and the highly anomalous cores of the mafic intrusions. Joint inversion allows for finer differentiation of the inverted physical property values into more distinct geologic units. Our work extends the previous 2D interpretations into 3D, allowing for a better understanding of the igneous intrusion in the study area in all three dimensions. Although we cannot verify our differentiation results, the differentiated geologic units represent testable hypotheses and provide guidance for future drilling activities and geophysical data acquisition. With the availability of high-resolution airborne geophysical data, the ever-increasing computation power and the constant developments in inversion methodologies, we expect 3D geology differentiation to be performed more often in the future, and our work demonstrates the feasibility and benefits of such efforts.


We thank H. Hu for his help with the MATLAB visualization that facilitated interpretation of the inverted models. We also would like to thank the reviewers B. Drenth, H. Ugalde, and Y. Li as well as AE G. Liu and ASE I. Filina for their constructive comments that helped improve the manuscript.

Data and materials availability

Data associated with this research are available and can be accessed via the following URL https://mrdata.usgs.gov/magnetic/show-survey.php?id=IA_10002.

Jiajia Sun recevied a B.S. (2008) in geophysics from China University of Geosciences (Wuhan), and a Ph.D. (2015) in geophysics from the Colorado School of Mines. He was postdoctoral fellow at the Colorado School of Mines from 2015 to 2017, and he has been an assistant professor of geophysics in the Department of Earth and Atmospheric Sciences at the University of Houston since 2017. His research interests include (1) tackling the magnetic remanence problem by 3D magnetization vector inversions, (2) developing joint inversion algorithms for integrated imaging of the earth based on multiphysical geoscience data sets, (3) quantifying uncertainties of geophysical inversions in deterministic and Bayesian inversion frameworks, and (4) automating geophysical data interpretation via machine learning. He received honorable mention for Best Paper in Geophysics in 2015, and Best Paper in the Mining sessions at the 2016 SEG Annual Meeting. He is an active member of SEG, a member of AGU, EAGE, and GSH, and serves on the SEG Gravity and Magnetics Committee and the Mining Committee.

Aline Tavares Melo received a B.S. (2008) in geology an M.S. (2012) in geophysics from the Universidade de Brasilia (Brazil), and a Ph.D. (2018) in geophysics at the Colorado School of Mines. She has worked as a mining geophysicist for Vale S.A. from 2008 to 2013. She is currently an assistant professor at the Universidade Federal de Minas Gerais. Her research interest is integrated quantitative interpretation of multiple geophysical data for geologic differentiation, inversion, mineral exploration, and machine learning. She is a cofounder of the Brazilian Association of Women in Geosciences.

Jae Deok Kim received a B.S. (2017) in energy resources engineering from Seoul National University, and an M.S. (2020) in geophysics from the University of Houston. He is currently a Ph.D. student at the MIT-WHOI Joint Program focusing on electromagnetic methods. His career in geophysics started in 2016 as a research intern at the Applied Geophysics Laboratory of Ecole Polytechnique in Montreal, Canada. More recently, he worked as a research assistant at the University of Houston, researching integrated interpretation methods in mineral exploration using joint inversions and geology differentiation. His research interests lie in quantitative integrated methods for the interpretation of multiple types of geophysical and geologic data using joint inversions, geology differentiation, and Bayesian methods; his areas of interest include resource exploration, environmental and hydrological problems, and the structure of the crust and mantle. He has been the recipient of numerous accolades, including the 2019–2020 SEG Robert E. and Margaret S. Sheriff Scholarship, the 2019 SEG/ExxonMobil Student Education Program, and the University of Houston Outstanding Academic Achievement Award. He is a member of SEG and AGU.

Xiaolong Wei received a B.S. (2015) in geophysics from the China University of Geosciences, Beijing, and an M.S. (2018) in mineral exploration from Northwest University, Xi’an, China. He is currently pursuing a Ph.D. in geophysics from the University of Houston. Since 2019, he has been a research assistant with the Department of Earth and Atmospheric Sciences, the University of Houston. His current research interests include geophysical inverse problems, joint inversion, and uncertainty analysis in both deterministic and Bayesian inversion frameworks. He is a member of SEG.

Freely available online through the SEG open-access option.