# Geophysics

- © 2006 Society of Exploration Geophysicists. All rights reserved.

## Abstract

Gravity and magnetic data are inverted jointly in terms of a model consisting of an interface separating two layers having a constant density and magnetization contrast. A damped least-squares inversion is used to determine the topography of the interface. The inversion requires knowledge of the physical property contrasts across the interface and its average depth. Since the relationship between model parameters and data is weakly nonlinear, a constant damped least-squares inverse is used during the iterative solution search. The effect of this inverse is closely related to a downward continuation of the field to the average interface depth. The method is used to map the base of the Sept-Iles mafic intrusion, Quebec, Canada, and the shape of the central uplift at the Chicxulub impact crater, Yucatan, Mexico. At Sept-Iles, the intrusion reaches a thickness of

## INTRODUCTION

Multiple data sets offer the possibility of reducing ambiguity and nonuniqueness in geophysical data interpretation. The joint inversion of two data sets requires that each set of observations is caused by the same source, i.e., the spatial variation of each physical property is the same or, as commonly modeled, sources with constant physical properties have the same boundaries. If a magnetic field and gravity field are caused by the same source in which the magnetization direction and ratio of magnetization intensity to density are both constant, they are related by Poisson's relation. Hence, for noise-free data, one field can be calculated from the other. Thus, the utility of jointly inverting gravity and magnetic data would appear to be limited. However, where coincident gravity and magnetic anomalies appear to be caused by the same source, joint inversion provides a means to quantitatively assess whether this is indeed the case. If inversion results in good fits to both data sets, then a common source is indicated. However, poor fits to portions of either data set will pinpoint those parts of the source where physical property variations do not vary in concert and/or the boundaries of the sources causing the gravity and magnetic fields do not coincide. In these areas, the geological model must be adjusted or alternative interpretation methods could be used to satisfactorily model each data set independently. Bott and Ingles (1972) use an equivalent-layer approach to investigate the variation in magnetization/density ratio of sources without direct inversion of the data. Their method is independent of source type and is based on correlating coincident pseudogravity (calculated from the magnetic field) and gravity profiles. This approach does highlight variations in magnetization/density ratios and, consequently, implied variations in source boundaries; however, joint inversion provides additional information in the form of an acceptable geological model for the common source.

For many geological bodies, it is unlikely that magnetization (induced or remanent) will vary in exactly the same fashion as density. Density is primarily controlled by lithology, i.e., silicate mineralogy. In contrast, magnetization, which is predominantly determined by the occurrence of magnetite, is less directly related to lithology but is strongly influenced by metamorphism, alteration, and emplacement conditions. Variations in density are also likely to be much smaller in magnitude than those in magnetization, so the ratio between the two is expected to vary, e.g., between different phases of a single pluton. Consequently, a joint inversion of the gravity and magnetic data can delineate areas where the common source assumption is appropriate and where it is not. This provides useful geological information in terms of physical-property variation and source structure. The difference in resolution power between magnetic and gravity data also means that each data set responds differently to sources at various depths and/or levels of structural complexity. Magnetic data contain more high-frequency content than gravity data, so they are often more effective at resolving shallow or complex structures.

Several approaches to inverting gravity and magnetic data jointly are described in the literature. Menichetti and Guillen (1983) use a generalized inversion technique to determine body shapes for the 2.5D case with specified densities and magnetizations. A stable inversion is achieved by removing the effects of small eigenvalues, providing a suggested improvement over separate inversion of either data set. Serpa and Cook (1984) use a least-squares inversion method to solve for the vertices and physical properties of a 2.5D model. Their method does not include damping or regularization, so stabilizing the solution is effected by holding some parameters fixed while inverting for others. Zeyen and Pous (1993) use a 3D subsurface model comprising of an ensemble of vertical prisms and solved for physical properties and the tops and bottoms of each prism. They control solution stability by specifying a priori parameter and data covariances. Gallardo-Delgado et al. (2003) extend the 3D approach to include a continuous density variation with depth and the direction of magnetization as unknown parameters. Within a statistical framework, Monte Carlo methods have been used for joint gravity and magnetic inversion producing an a posteriori probability density function describing a suite of acceptable models rather than a single inverted model (Bosch, 1999; Bosch et al., 2001).

All of these inversion methods are based on solving a combined set of equations relating the model to both data sets. Depending on the model's definition, the problem can be posed with the number of data points either greater or smaller than the number of model parameters. In either case, the problem can be solved by minimizing (1) the difference between both observed and predicted data sets plus, for solution stability, (2) a measure of model length or roughness. This leads to methods that find a minimum norm solution (Menke, 1989). Tighter constraints on linking the different physical properties being solved for can be incorporated into the inversion by introducing some measure of structure into the model. Haber and Oldenburg (1997) define a measure of the model structure based on its curvature and then minimize the difference in structure between the two models subject to the constraint that both data sets are fit to within their expected errors. A similar approach is used by Gallardo and Meju (2003), who use a structure measure based on the vertical and horizontal gradients of the physical properties being solved for. Both of these approaches are again posed as minimization problems, with extra terms for structure added to the standard terms of data misfit and model length.

In this study, gravity and magnetic data are inverted jointly using a damped least-squares approach. The standard way of dealing with nonlinearity is by reevaluating matrix quantities at each model estimate during an iterative search for an acceptable solution. Since matrix evaluation and/or inversion is computationally expensive for large data sets and models, a fixed inverse with a simple form is used throughout the search. This approach is demonstrated on synthetic data and applied to coincident gravity and magnetic anomalies over the Sept-Iles mafic intrusion in Quebec, eastern Canada, and the Chicxulub impact crater, Mexico. For the two real data examples, joint inversion is used to determine how well the anomalies are described in terms of a common source and to focus on areas where this assumption breaks down.

## METHOD

The model used in the following consists of two homogeneous layers, each having a constant density and magnetization, separated by an interface whose topography is to be determined. The nonuniqueness associated with this kind of model is predominantly related to the depth and physical-property contrast of the interface (see Skeels, 1947). For a constant interface depth, varying the physical property contrasts between the two layers results in a range of models with an inverse relationship between the interface topography and the density or magnetization, all fitting the data with the same precision. For the gravity case, in general, the average interface depth need not be known a priori; but for the algorithm used below, it must be specified, since this parameter is not solved for. Similarly, for the magnetic case, the average interface depth must be known since the magnetic data do not provide information on its value. Consequently, the average interface depth and the density and magnetization contrast between the two layers are specified prior to inversion.

The expressions relating to gravity data **g** and magnetic data **m** to the model **x** are given by Parker (1973):
**j** and **t** are unit vectors in the direction of the magnetization and the earth's magnetic field, *ρ* is the density contrast, *J* is magnetization contrast, *G* is the gravitational constant, F[ ] denotes Fourier transform, and **k** is the wavevector. Vector **K** is constructed from wavevector **k**, whose components in the *z* is measured positive downward, and **x** is the interface topography measured (positive upward) with respect to *z*. Since equations 1 and 2 consist of a series of Fourier transforms, **g**, **m**, and **x** must be located on coincident grids and the appropriate treatment applied to minimize edge effects. Linearized versions of the above relations are, in partitioned matrix form,
**G** and **M** can be written as (Pilkington and Crossley, 1986)
**E** and **Γ** and **Λ** are diagonal with elements
**G** and **M** in equations 5 and 6 results from taking the first term in the series in equations 1 and 2 and setting **x** to zero. Since **E** and **Γ** and **Λ** are diagonal, the factorizations in equations 5 and 6 are equivalent to the singular value decompositions (SVD) of **G** and **M**. Hence, **Γ** and **Λ** contain the eigenvalues of **G** and **M**, respectively.

Because the data sets have different units and magnitudes, one or both of them must be weighted to balance their contribution in the inversion. If, for example, numerical values of **M** and **m** are much greater than **G** and **g** in equation 3, the inversion will be dominated by the magnetic data and information from the gravity data will be lost. If both data sets are assumed to be contaminated by independent Gaussian errors, they can be transformed into sets of dimensionless values by dividing each observation by its expected standard deviation (Vozoff and Jupp, 1975; Serpa and Cook, 1984). If errors are dependent, the data transformation requires knowledge of the data covariances but is still implemented straightforwardly (Jackson, 1973). A similar weighting method consists of normalizing parameter and observed data sets by a priori model values (Hering et al., 1995). Finally, empirical weighting guided by estimates of data variances can be used. This is a flexible approach in that the weighting can be adjusted even if it requires that weights are not in agreement with estimated data variances (Lees and VanDecar, 1991; Zeyen and Pous, 1993).

In the statistical approach to inversion, weighting is quantified by the ratio of the standard deviation of the magnetic data errors to that of the gravity data. Determining reliable and usable estimates of data errors is difficult, to say the least, since the data used in inversion contain not only raw measurement errors but also the combined effects of error introduced by various data-reduction steps and, most serious of all, gridding. Consequently, estimates of data errors can be used to determine the weighting; but if it is clear that one data set is dominating the inversion solution, then the weighting should be adjusted accordingly. The empirical approach adopted here is implemented by weighting just one set of equations, so that equation 3 becomes
*w* is the scalar weighting factor. An alternative approach that avoids the weighting issue is sequential inversion (Lines et al., 1988), in which one data set is inverted and the results are used as the starting model for inversion of the second data set. Coupling the two data sets may, however, prove weaker for this method.

Equation 9 is an overconstrained system solved by forming the damped least-squares solution (Menke, 1989) for the equivalent equation 4:
**I** is the identity matrix, and *θ* is the damping factor. Substituting for **G** and **M** using equations 5 and 6
**A** and forming its generalized inverse. Because **Λ**, and **x** can be carried out efficiently in the wavenumber domain using gridded data without any need for explicit matrix inversion. Nonlinearity is taken care of by using the iterative scheme
*k* is iteration,

Use of a constant inverse in equation 14 is possible and effective because of the weak nonlinearity of the expressions in equations 1 and 2. Updating *θ*, which is kept constant throughout the iterative procedure, also contributes to the stability of the iterations by suppressing the effects of small eigenvalues to the current model estimate. Any small eigenvalues in **G** and **M** lead to large numerical values of the terms in parentheses in equations 12 and 13, containing the inverse eigenvalues. If the latter become large, this can cause unacceptably high-amplitude, oscillatory behavior in the model and equation 14 diverges. Damping is always necessary because of the exponential terms in the eigenvalues (equations 7 and 8). The inversion is equivalent to a downward continuation of the field and conversion of the field at the average interface depth into topography through division by the nonexponential terms in equations 7 and 8. In the gravity case, this is equal to converting the continued field into topography using the Bouguer slab formula.

The effect of introducing damping through addition of *θ* is clarified when the damped inverse eigenvalues are compared with the undamped versions. Since the index *i* in equations 7 and 8 corresponds to wave vector, forming the ratio of damped to undamped inverse eigenvalues quantifies the frequency response of the smoothing introduced by the damping factor. Figure 1 shows an example case where *θ* is numerically equal to one of the eigenvalues (*θ* from the range of

## SYNTHETIC DATA EXAMPLE

The application of the iterative inversion scheme (equation 14) is demonstrated on a synthetic data set. Figure 2a shows a model topography with an average depth of *w*. The inversion in Figure 2 uses

## SEPT-ILES MAFIC INTRUSION, QUEBEC

The Sept-Iles mafic intrusion is located on the north shore of the St. Lawrence River in Quebec, Canada. Cambrian in age, it has a diameter of

So far, interpretation of geophysical data over Sept-Iles has consisted of 2D modeling of gravity data only. Goodacre and Nyland (1966) model the gravity anomaly as a vertical cylinder with a radius of

The likelihood that both the gravity and magnetic fields over the Sept-Iles intrusion are predominantly the result of the same structure is demonstrated by the good coincidence of the two in Figure 5. Joint inversion would appear appropriate in this case. However, the observed layering of the intrusion and consequent variation of physical properties throughout the body means that the assumption of a single density/magnetization contrast for the inversion is probably too simplistic. Nevertheless, joint modeling of the observations using the two-layer model with just a single and constant physical property contrast can still provide useful information in terms of pinpointing areas where sources have more complex property variations.

Paleomagnetic studies of the Sept-Iles intrusive rocks indicate a stable remanent magnetization component (inclination

Compared to the Loncarevic et al. (1990) 2D model, which shows a simple increase in depth toward the center of the intrusion, the model in Figure 6a has a more complex form, with its deepest point coinciding with the maximum gravity anomaly south-southeast of the intrusion center. Secondary features are the two crescent-shaped troughs to the north of this point, the inner coinciding with a magnetic and gravity high and the outer associated with only a magnetic high, on the northern edge of the body (see Figure 5). Both troughs are the result of the more dense and more magnetic late gabbro phase that forms ring dikes of varying thickness (Loncarevic et al., 1990). They show positive magnetic residuals but are well fit in the gravity case (Figure 6b and c), suggesting a higher magnetization value is appropriate for these parts of the intrusion. Over the gravity maximum (

In the case of the Sept-Iles intrusion, assuming a constant source density and susceptibility leads to artifacts in the solution that would not appear in a model where physical properties are variable. Nevertheless, the joint inversion demonstrates that most of the observed gravity and magnetic data result from the same source. Inspection of the residuals for both data types delineates areas where property variations and/or source boundaries do not vary in concert and indicates how they deviate. This kind of information can form the basis for subsequent detailed modeling.

## CHICXULUB IMPACT STRUCTURE, YUCATAN

The second data example is taken from the central portion of the Chicxulub impact structure located on the northwestern Yucatan coastline. This

Figure 7a and b shows the residual gravity and magnetic fields used in inversion with the corresponding calculated fields from the topography solution in Figure 8a, given in Figure 7c and d. Initial density and susceptibility values were set to values previously used in 2D modeling. These were adjusted to

Figure 9 shows the behavior of the gravity and magnetic data fits as a function of varying density and susceptibility values and carrying out inversions. The particular solution depicted in Figure 8a represents a compromise between the two data fits, with rms values of

Relaxing some of the prior assumptions, e.g., the constant density and magnetization, would likely produce a better fit to the data. Since there is no strong observational or theoretical evidence to support the existence of more than one peak at central uplifts in impact craters, the gravity data could be better modeled with a variable density structure. Preexisting basement density contrasts will persist through uplift formation, so a single constant density as assumed here may not be the optimal choice. Conversely, the magnetic data appear compatible with a constant magnetization throughout the uplift.

## CONCLUSIONS

For simple two-layer models, gravity and magnetic data can be inverted jointly using a frequency-domain iterative scheme requiring no matrix inversion. The frequency response of the inverse eigenvalues used in the inversion algorithm (equation 14) show that the magnetic data contribute more short-wavelength information to the solution than the gravity data. The synthetic example also demonstrates this point.

The inversion of data sets from an intrusion and an impact crater central uplift shows that a simple common source assumption for these two bodies could explain most of the observed anomalies. The presence of residuals for either data type indicates areas that could not be satisfactorily explained in terms of the assumptions adopted. These regions are most likely caused by spatial variations in either or both density and magnetization. If such a subarea can be isolated, in which physical properties are constant (but not equal to the initial inversion values), then these property values can be adjusted and the data reinverted. Where the physical-property variations exist over small scales, alternative modeling methods must be used. The two data examples show that even the restrictive assumption of constant physical properties for the models can result in meaningful structural interpretations and provide a solid starting point for more detailed subsequent modeling if needed.

## ACKNOWLEDGMENTS

John Millar (University of Calgary) is thanked for providing the residual gravity data over Chicxulub. Two reviewers and the associate editor (J. Silva) provided detailed comments and suggestions that improved an earlier version of the manuscript. This is Geological Survey of Canada contribution 2004016.

- Received May 14, 2004.
- Revision received June 15, 2005.