# The Leading Edge

- Copyright © 2001 Society of Exploration Geophysicists

Interpreting exploration data requires combining different types of information to solve the geologic puzzle. It implies bringing together all data components into an image that makes conceptual sense in terms of the geology of the exploration area. The identification of geologic objects and the inference of a spatial description of the lithology—consistent with all available information—are the objectives of the process.

The more information utilized, the more certain is the result of the inference. The interpretational process should serve to combine different types of geophysical data, petrophysical information on the rock properties, and information on the geology of the area. It is an interdisciplinary process, performed by an expert or a multidisciplinary team, the desired outcome being a common earth model, or class of models, known to satisfy both the conceptual geologic understanding and all available geoscientific data.

The object of the present work is to incorporate geophysical inversion methodologies as a tool in geologic interpretation, describing a process to jointly invert gravity and magnetic data that takes into account petrophysical and geologic constraints. Basically, the method seeks lithologic models that explain the overall data, helping with the task of quantitatively reconciling the available geologic and geophysical information. The process uses a geostatistical model to couple the lithology with realistic density and magnetic susceptibilities. Hence the values of density and magnetic susceptibilities are a priori conditioned by the lithology, avoiding unrealistic excursions allowed by common inversion approaches. On the other hand, prior information about the lithology, such as a geologic surface map or an interpreted drill hole, is incorporated in the inversion and satisfied by the resulting models.

## Joint model and multiple data sets

The first difference with common inversion approaches is that here we describe the media with a joint model—i.e., a model that simultaneously describes lithology, density, and magnetic susceptibility. The major lithologic categories in the area are identified from the surface geology or drill holes and grouped into lithotypes. Their subsurface geometry is postulated, initially on the basis of geologic interpretation alone, forming regions that should be consistent with the known surface geologic map and drill holes, if available.

Once the lithotypes have been defined, a statistical description of the density and the magnetic susceptibility is built for each region, based on laboratory measurements of rock samples, wireline petrophysical logging, or commonly reported data for the relevant rock types. The statistical description is independent for each lithologic region and provides control on mean values of the properties, their dispersion and correlation, and the spatial continuity of the property (covariance and cross-covariance functions) in each region. Briefly, the statistical model provides the link between the lithologic configuration of the subsurface and the physical properties (here, density and magnetic susceptibility) filling the subsurface regions.

As an illustration, we show a test model in Figure 1, together with the gravity and magnetic data calculated from it. The model describes both physical media properties and lithology. The density and magnetic susceptibility fields in the figure have been simulated from the lithologic field at the bottom of the figure, following a geostatistical model. The values of the magnetic susceptibility and the density obey different statistics (Figure 2).

Additionally, the geostatistical model is able to describe the texture of the property fields within each lithologic region. Figure 1 shows different ranges of the spatial continuity of the properties for each lithotype region. The density and the magnetic susceptibility change smoothly inside the anorthositic body, showing large areas of low susceptibility or high scatter from the mean value of the property, with smooth transitions between them. These spatial variations have a shorter range within the olivine gabbro body, because a shorter correlation distance (i.e., covariance function range) was used to describe the spatial continuity in this lithotype region than in the anorthositic one.

Simulating the gravity and magnetic data from the density and magnetic susceptibility subsurface configurations, correspondingly, makes the connection between the joint model and the geophysical data. Summarizing, we solve the forward problem, from the lithologic model to the geophysical data, as follows: Given a particular lithotype subsurface configuration, the geostatistical model provides reasonable simulations of what the density and the magnetic susceptibility should be. From them, we can calculate the gravity and magnetic fields with geophysical modeling.

## The inversion

A sufficiently general mathematical language is needed for formulating the inverse problem, as we are dealing with different types of information including geology, petrophysics, and geophysics. The mathematical language required is the language of probabilities. A given subsurface configuration (here, lithotype + density + magnetic susceptibility) has an associated probability density that will be larger as the calculated and true data are more closely fit and as the prior petrophysical and geologic information is better fulfilled. This provides a quantitative way of measuring how probable a particular model is. An arbitrary configuration of the joint model will be likely if the calculated geophysical data (forward model) fits the observed data within its associated error, and if the properties inside lithotype regions obey the prior distribution given for the type of rock in the area.

We suggest the review of previous work for a detailed description of the mathematics and algorithms behind the technique. Briefly, the probability density is built up as a product of three factors: (1) a likelihood function incorporating information about the misfit between the calculated and observed data; (2) a probability density describing the dependence of the physical properties on the lithology, according to the available petrophysical information; and (3) a probability density describing the prior information on the interpreted geologic structure of the area.

Following well-established statistical techniques, an algorithm is set up to generate joint models in proportion to the probability density function—i.e., models that fit the observed data and simultaneously satisfy the geologic and petrophysical requirements. This process forces consistency between the solution of the inverse problem and the conceptual realism of the solutions in terms of the geology and the physical property behavior—i.e., the resulting inverse models not only match the geophysical data but also, critically, “look right” to the geologist. The algorithm starts from an initial guess of the structure and provides modifications on an iterative basis. Each candidate model is then accepted or rejected according to the probability density computed from the data and prior constraints.

We illustrate the inversion method by taking the model in Figure 1 as a “true” model for a synthetic test. The goal is to infer the joint model from the “true” gravity and magnetic data shown in the same figure under the following prior constraints: (1) the set of lithotypes is fixed; (2) the lithotype map at the earth surface is known and fixed; and (3) the density and magnetic susceptibility should follow the distributions in Figure 2.

The method was implemented here in 2-D (depth versus horizontal direction) assuming constant properties along the third spatial direction and using a parameterization based on a triangulation of the 2-D section. The triangulation allows for a flexible representation of the boundaries between the lithologic regions. Properties are homogeneous within individual triangles but may be heterogeneous within each lithologic region, as already discussed. The method starts from an initial lithological model (Figure 3) and proceeds by iteratively modifying the model. Permitted modifications include changing the physical property of a triangle, changing the lithotype within a triangle, and moving the triangle vertices.

Figure 4 shows a joint model generated with the method, assigning a 5% uncertainty to the “true” data. The model fits the “true” data in Figure 1 within the uncertainties. It has the same lithology at the surface as the “true” model, and it fulfills the hypothesis on the distribution of the density and magnetic susceptibility, as shown in Figure 2. As more iterations are produced, the method generates many joint model simulations, with differences between them, but all consistent with the overall information.

After generating a large number of models consistent with the data, the method allows constructing plots of probabilities (Figure 5). The figure shows expected probabilities of where to find each of the three lithotypes in the model. Other probability plots could be produced for any other variable of interest (e.g., the depth to a lithologic region at some point or the relative volume of each of the lithotypes). Thus, interpretational decisions are firmly put on a quantitative statistical basis and the foundation is laid for quantitatively connecting exploration business decisions (“Should we drill here?”) to risk.

## Results for the Kiglapait data

Kiglapait is a region of mining exploration interest in Labrador, Canada, in the general region of the Voisey Bay nickel discovery, one of the world's most significant economic mineral discoveries of the past decade. An earth model was built using geologic mapping and ground geophysical data (all data shown here are courtesy of the Geological Survey of Canada), demonstrating that quantitatively constrained subsurface geologic modeling does not necessarily require drill holes or other subsurface data.

The large rectangular box containing the model, shown throughout Figure 6, is 40 km on a side. In Figure 6a, the two map contacts divide the area into three lithologic units: (1) an inner intrusive layer of olivine gabbro to ferrosyenite; (2) an outer intrusive layer of troctolite; and (3) an undifferentiated anorthositic host. The formation of these rock types is associated with magma differentiation and solidification. The three lithotypes have different density and magnetic susceptibility distributions, the olivine gabbro having the highest magnetic susceptibility and density of the three. Prior values to construct physical property distributions, shown in Figure 2, were taken from common reported values of the properties of these types of rocks.

The first step was construction of an initial model of the geologic structure in Figure 6. This model summarizes information on the geologic chart of the area, structural data on the layering orientation and the geologist's knowledge-based intuition. Figure 6b illustrates a cut section plane through the geologic model, used in the remainder of this paper as the section on which the 2-D modeling and inversion are described (this section is also illustrated in Figure 3). The goal of the geophysical inversion presented here is to define 2-D lithologic models based on the geologic structure shown along the section in Figure 6b—done initially without regard to the geophysical datA—and quantitatively constrained by the geophysical data and the physical rock properties.

The lithology at the surface is known from the geologic map of the area, and the inversion will use this information. Another interesting aspect of the geology is that as a result of the magna differentiation, the troctolite gabbro is an intermediary type of rock between the olivine gabbro and the anorthositic rocks. It is unlikely that the olivine gabbro could be in direct contact with the anorthositic rocks. This logic from the geology was also included in the inversion: lithologic models were not allowed to have olivine and anorthositic rocks in contact, without the intermediary troctolite gabbro between them.

Figure 7 shows two joint models generated by the inversion technique. The joint models fully satisfy the fixed lithology at the surface, the distributions of the density and magnetic susceptibility, and fit both types of data (gravimetric and magnetic) within uncertainties. Figure 8 shows the values of the properties within these models and Figure 9 the calculated probabilities for the lithotypes, according to position.

There are significant differences between the inferred configuration and the preliminary model (Figures 3 and 6b). The inferred configuration fully satisfies the observed data and the ensemble of petrophysical and geologic constraints, whereas the initial model satisfied only the known lithology at the surface, assumed continuity of down-plunge projection, and geologic intuition. An interesting feature in the result is the assymmetry in the location in depth of the bodies of troctolite and olivine gabbro to ferrosyenite. The major intrusive body extends to the right-hand side of the figure beneath its outcrop, in a direction consistent with the prolongation of the observed gravity anomaly in Figure 6c. There is a large probability of finding a buried body of gabbro to ferrosyenite, which is the lithotype exhibiting the highest density and magnetic susceptibility, below the indicated horizontal position of 47 km and within a depth range of 3–5 km. These features are a consequence of the asymmetry of the gravity and magnetic observed anomalies, together with the prior constraints on the property values, according to the lithotypes.

## Discussion and conclusions

Joint inversion of multiple types of geophysical data under lithologic constraints provides a method of integrating the multidisciplinary data upon which exploration teams base interpretation. It is a quantitative method that can summarize large amounts of data on probability-based assessments about the geologic structure. We believe that the results provided by the method have greater value than standard methods of “overlaying” multidisciplinary data, because the method exploits the structured relationships between data of different types. It is an explicit, quantitative method of integrating different types of information, in which the uncertainties and distributions provide the relative importance of each component in the inversion. In this manner, it helps in the understanding and reconciliation of information across the different disciplines involved in exploration.

As lithology is part of the model, the results depend on the starting hypothesis about the lithologic structure of the area, and in particular about the way in which lithology is described (e.g, the lithotypes considered in the model). Nevertheless, in cases in which competing initial hypotheses are justifiable, each one can be treated independently and the results compared.

The description of the physical rock properties, as a function of the lithology, is an important component of the method. Results naturally depend on the prior distributions providing the link between the lithologic and physical property model of the subsurface. A rough characterization of the distributions can be based on the average reported description for the types of rocks involved. Generally, however, we recommend acquisition of petrophysical data from the area under study as the most precise way to characterize the distributions, via in-situ logging or laboratory measurements of collected samples. Physical rock properties remain the link between geologic and geophysical data, and with present methodology in mind, petrophysical characterization of the area should be considered a critically important complement to geophysical exploration surveys.

## Suggested reading

Details on the mathematics and calculation behind the methodology are described in “Lithologic tomography: From plural geophysical data to lithology estimation” by Bosch (*Journal of Geophysical Research*, 1999). A description of the method with a detailed presentation of a field case can be found in “Lithologic tomography: An application to the Cote d'Armor region in French Brittany” by Bosch et al. (*Tectonophysics*, 2000). “Lithology discrimination from physical rock properties” by Bosch et al. (Geophysics, 2000) shows successful prediction of rock lithology from several physical properties measured in rock samples. An overview of the formulation of statistical techniques to solve geophysical inverse problems can be found in “Monte Carlo sampling of solutions to inverse problems” by Mosegaard and Tarantola (*Journal of Geophysical Research*, 1995). For data sources and previous interpretation regarding the Kiglapait intrusive, see “3-D visualization of structural field data and regional subsurface modeling for mineral exploration” by De Kemp and Desnoyers (Proceedings of Exploration 97: Fourth Decennial Conference on Mineral Exploration, 1997), and *The Kiglapait Layered Intrusion* by Morris (Geological Society of America, Memoir 112).

## Footnotes

Miguel Bosch is a visiting fellow of the Department of Earth Sciences at the University of Cambridge.