Colin G. Farquharson:  Research

I'm interested in all aspects of forward-modelling and inversion of geophysical data. Below are some specific topics I'm working on now or have done so in the past.

Joint inversion: Inverting two or more data-sets for some kind of common Earth model is arguably the most active research area in geophysical inversion at the moment. I am very lucky to be working with Dr. Peter Lelièvre, an NSERC Post-Doctoral Fellow in my group. (Peter's home page is here.) Peter has been working on the joint inversion of seismic travel-time and gravity data. Our current strategy is to optimize an objective function that comprises the usual measures of data misfit and model structure for the individual data-sets and physical properties plus a measure of model similarity. We are using a combination of one or more of the following for the measure of similarity: explicit analytic relationship, implicit analytic relationship, statistical relationship. More below.

Different discretizations of the subsurface: The usual way at the moment for we geophysicists to discretize the Earth, and thus obtain a model of the spatial variation of some physical property in the subsurface, is to use a regular rectangular grid. This is by far the easiest kind of mesh to deal with, especially in terms of the information required to keep track of where all the cells and nodes are in the model. However, this is perhaps not the best kind of mesh for representing the possible geological features in the Earth. Wouldn't it be better to have a mesh that could follow known (or projected) interfaces such as faults and contacts between rock units? More below. (No, the picture to the left, which is of the wire-frame geological model of the Ovoid deposit at Voisey's Bay, Labrador, is not squint.)

Blocky models in minimum-structure inversions: Minimum-structure inversions (in which the model is finely parameterized, and a measure of model complicatedness is minimized along with a measure of data-misfit) have many good points, notably reliability and robustness. Typical, traditional implementations use an l2, or sum-of-squares, measure for the complicatedness. This is why the models which are produced have their characteristic fuzzy, smeared-out appearance. If an l1 measure is used instead, models with sharp interfaces between zones of uniform physical parameter can be constructed. Also, if diagonal finite-differences are explicitly included in the measure of model structure, angled interfaces can be generated in the model. (If only finite differences in the x-, y- and z-directions are included, only horizontal and vertical interfaces are generated.) More below.

3-D magnetotelluric inversion: While at UBC-GIF, I wrote the MT-relevant bits of the 3-D MT inversion program MT3Dinv. This program shares components with the controlled-source EM forward-modelling and inversion programs EH3D and EH3Dinv that were developed by Eldad Haber, Roman Shekhtman and Doug Oldenburg. All these programs came out of the IMAGE consortium research project at UBC-GIF. I have since been testing MT3Dinv on various data-sets. For example, in conjunction with Jim Craven (GSC, Ottawa), I have inverted the AMT data-set from over the McArthur River uranium mine in Saskatchewan. The results are (quite?) good. But 3-D inversion of multi-frequency EM data is still so computationally expensive that suitably discretized models, and appropriately converged iterative linearized inversion procedures, are still elusive. More below.

Integral equation solution for large contrasts: Traditional implementations of the integral equation solution to the geophysical electromagnetic forward-modelling problem fail for large conductivity contrasts. I have developed an implementation that does not have this limitation. It has been tested against results from physical scale modelling performed by Ken Duckworth formerly at the University of Calgary. More below.

Program EM1DTM: This is the sister program to EM1DFM, which I wrote a few years ago while at UBC-GIF. (EM1DFM is available for commercial and academic licensing: see the relevant UBC-GIF page.) Program EM1DTM is for time-domain EM data, particularly lines or grids of many soundings (that is, what one gets from airborne surveys). Like EM1DFM, nifty things such as the GCV or L-curve criteria can be used to automatically estimate the trade-off parameter. Unlike EM1DFM though, EM1DTM includes general measures, and so can produce blocky, piecewise-constant models if desired. More below.

Joint inversion

There is just one subsurface beneath our feet. Shouldn't we be able to construct an Earth model via inversion that has values of multiple physical properties in each cell that both reproduce whatever geophysical survey data we may have and make geological sense even on a cell by cell basis? If there is an interface between two units that have different densities and magnetic susceptibilities, for example, that interface should appear in our model and the cell densities and susceptibilities on either side of the interface should be appropriate for the two units. But of course an interface between two units that differ only in, say, susceptibility and not density should also be present in the model with appropriate values of cell density and susceptibilities on either side of the interface.

Dr. Peter Lelièvre and I are investigating the use of "compositional" measures of model similarity within an otherwise standard minimum-structure inversion procedure. The objective function being minimised contains a data-misfit term for each geophysical data-set being considered, a measure of model structure for each physical property in the model, and the measure of model similarity. Our main application at the moment is to mineral exploration. In typical exploration scenarios, whether truly green-field prospecting or resource delineation, the number of possible (hoped for?) geological units and zones of mineralization within the volume of subsurface under investigation are limited. Also, the possible units and their physical properties can be anticipated reasonably well. Therefore, as measures of model similarity, Peter and I are trying explicit analytic relationships, implicit relationships such as measures of correlation, and statistical relationships such as clustering.

Below are images for inversions of synthetic gravity and seismic travel-time data for the Voisey's Bay Ovoid deposit. The images on the left show, using iso-surfaces, the anomalous density and slowness models constructed from independent inversions. The anomalous slowness is nicely concentrated where it should be, but the density anomaly is far too spread out. The images on the right show the results of a joint inversion. The density and slowness anomalies are now much more coincident.

References:

Lelièvre, P.G., C.G. Farquharson, and C.A. Hurich, 2012. Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration, Geophysics, 77, K1.

Farquharson, C.G., P.G. Lelièvre, and C.A. Hurich, Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration, Fall Meeting, AGU, San Francisco, Calif., 13-17 December, 2010.

Lelièvre, P.G., C.G. Farquharson, and C.A. Hurich, Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration, SEG Annual Meeting, Denver, 17-22 October, 2010. [PDF file of talk (5.9M).]

Lelièvre, P.G., C.G. Farquharson, and C.A. Hurich, Joint inversion of seismic travel times and gravity data on 3D unstructured grids with application to mineral exploration, EGU General Assembly, Vienna, 2-7 May, 2010. [PDF file of talk (4.0M).]

Back to top.

Different discretizations of the subsurface

The standard way at the present of generating a 3-D model of the subsurface that can be used for computing the response of a geophysical exploration technique is to divide the subsurface into a mesh of regular rectangular prisms. If these prisms are made small enough, then the model can, in principle, represent any spatial variation of a physical property. (Just like dots-per-inch for printing.) However, it would be better if the meshes we were working with could mimic the variability of real geological boundaries, such as faults and contacts between rock units, without having to resort to a fine stair-case of rectangular cells. This issue is particularly important as the 3-D models that geologist typically work with are wire-frame models that can represent complicated interfaces. Currently, either the geophysical model or the geological model has to be "transformed" so that it can be compared directly with the other. It would be much better, both in concept and in practice, if the meshing of our geophysical models was sophisticated and flexible enough to also be acceptable meshing for a geological model.

I've started work on this, in particular, investigating finite-element methods of forward modelling on unstructured tetrahedral meshes. Below is the wire-frame geological model of the Ovoid deposit at Voisey's Bay - an example of the kind of model that it would be good to be able to deal with in our geophysical forward-modelling and inversion programs - and an unstructured tetrahedral mesh of the good old block-in-a-halfspace model (generated by TetGen).

References:

Farquharson, C.G., P.G. Lelièvre, and C.A. Hurich, Unified geophysical and geological 3-D Earth models, EGU General Assembly, Vienna, 3-8 April, 2011. [PDF file of talk (1.2M).]

Back to top.

Blocky models in multi-dimensional minimum-structure inversions

Minimum-structure inversions generally involve discretizing the Earth using as fine a mesh as possible. The mesh remains fixed during the course of an inversion, with the values of the physical parameter in the cells being the quantities that are sought. The inversion is implemented by minimizing an objective function comprising both a measure of the amount of structure in the model as well as a measure of fit to the observations. Because the complexity of the model is kept in check by minimizing this composite objective function, minimum-structure inversions are usually reliable and robust.

The normal choice for the measure of the amount of structure in the model is an l2-measure of spatial derivatives of the model (because it results in a linear system of equations to solve). It is a consequence of the fundamental nature of the l2 measure that the resulting models have the familiar fuzzy, smeared-out appearance. If an l1 measure is used instead, models will be generated that have sharp interfaces between uniform regions. I have previously implemented l1 (and other) measures in 1-D inversion of controlled-source EM data. I have now implemented these measures in the 2-D magnetotelluric inverse problem, and the 3-D gravity inverse problem. Below are three different models produced by the inversion of the surface gravity data-set over the Ovoid deposit at Voisey's Bay, Labrador: the typical l2-measure inversion, an l1-measure inversion involving only finite differences in the x-, y- and z-directions, and an l1-measure inversion including also diagonal finite differences. (For technical details, see references below.)

References:

Farquharson, C.G, 2008. Constructing piecewise-constant models in multidimensional minimum-structure inversions, Geophysics, 73, K1-K9.

Farquharson, C.G., Constructing piecewise-constant models in multi-dimensional minimum-structure inversions, 76th Annual Meeting of the Society of Exploration Geophysicists, New Orleans, Louisiana, 1-6 October 2006. [PDF file of talk (2.2M).]

Farquharson, C.G., Blocky models in minimum-structure inversions, 18th International Workshop on Electromagnetic Induction in the Earth, El Vendrell, Spain, 17-23 September 2006. [PDF file of poster (0.5M).]

Back to top.

3-D magnetotelluric inversion

One of the products of the UBC-GIF IMAGE consortium research project was the 3-D magnetotelluric inversion program MT3Dinv. This was a combined effort of Eldad Haber, Roman Shekhtman and Doug Oldenburg and myself. (Other products of the IMAGE project were the related frequency-domain controlled-source EM forward-modelling and inversion programs EH3D and EH3Dinv.) The forward modelling is done using a finite-volume discretization, a primary-secondary field separation, and the relevant 2-D boundary conditions. The inversion algorithm is a basic minimum-structure, Gauss-Newton one, with the system of equations at each iteration solved (to some extent!) using conjugate gradients. The requisite action of the Jacobian matrix, or its transpose, on a vector is done by solving essentially the same system of equations as for the forward problem, but with the appropriate right-hand side. All these solutions of course try to exploit the sparseness of the forward-modelling matrix as much as possible.

One example data-set that I've had a go inverting with MT3Dinv is the data-set, collected as part of the EXTECH IV project, over the McArthur River uranium mine in Saskatchewan. This is in conjunction with Jim Craven at the GSC. Features of interest are the graphitic fault zones in the basement, and alteration halos in the sedimentary rocks above these fault zones. The graphitic faults have decent along-strike extent, and are conductive. The alteration halos are somewhat resistive relative to the sedimentary rocks, and are intermittent along the strike of the faults. Line-by-line 2-D inversions of the data show nice consistent conductors from one section to the next. The 3-D inversions (see images below) do show linear conductors, but with quite a bit of along-strike variability. (There's also the issue of the conductive features in my model being too deep, and the large amount of near surface structure (noise?).) This along-strike variability is actually consistent with the data. Despite the large-scale geology being 2-D, the finer-scale geology is probably more 3-D than the line-by-line 2-D inversions imply.

References:

Farquharson, C.G, and J.A. Craven. Three-dimensional inversion of magnetotelluric data for mineral exploration: An example from the McArthur River uranium deposit, Saskatchewan, Canada, Journal of Applied Geophysics, in press.

Craven, J.A., C.G. Farquharson, R.L. Mackie, W. Siripunvaraporn, V. Tuncer, and M.J. Unsworth. A comparison of two- and three-dimensional modelling of audiomagnetotelluric data collected at the world's richest uranium mine, Saskatchewan, Canada, 1ath International Workshop on Electromagnetic Induction in the Earth, El Vendrell, Spain, 17-23 September 2006. [PDF file of poster (1.9M).]

Farquharson, C.G., D.W. Oldenburg, and P. Kowalczyk. Three-dimensional inversion of magnetotelluric data from the Turquoise Ridge mine, Nevada, 74th Meeting of the SEG, Denver, Colorado, 10-15 October 2004.

Farquharson, C.G., D.W. Oldenburg, E. Haber and R. Shekhtman. An algorithm for the three-dimensional inversion of magnetotelluric data, 72nd Meeting of the SEG, Salt Lake City, Utah, 6-11 October 2002.

Back to top.

Integral equation solution for large contrasts

Traditional numerical solutions of the electric field integral equation for geophysical applications used pulse basis functions to represent the fields within the anomalous region, and used integrations over equivalent-volume spheres to approximate the integrations of both the current density within cells and surface charge density between cells. Such solutions give incorrect results for contrasts in conductivity between the anomalous region and the background that are greater than roughly 300:1. If the electric field within the anomalous region is instead represented by edge-element basis functions, and if the integration of surface charge density on interfaces across which there is a change in conductivity are explicitly included, a numerical implementation is possible which can successfully model the fields even when the conductivity contrast is high.

Results computed by this new integral equation formulation have been compared with the results of physical scale modelling experiments that were carried out by Ken Duckworth at his Electromagnetic Modelling Laboratory in the Department of Geology and Geophysics, University of Calgary. The figure below illustrates this for a graphite cube immersed in brine. The surface of the brine was 2cm above the top of the cube, and the transmitter-receiver pair was 2cm above the surface of the brine. (The conductivity of the graphite cube was 6.3E+04S/m, and that of the brine was 7.3S/m. The side-length of the cube was 14cm. The separation between transmitter and receiver was 20cm.) The heavy lines represent the scale modelling measurements; the narrow lines the numerical modelling values.

References:

Farquharson, C.G., K. Duckworth and D.W. Oldenburg, 2006. Comparison of integral equation and physical scale modeling of the electromagnetic responses of models with large conductivity contrasts, Geophysics, 71, G169-G177.

Farquharson, C.G., and D.W. Oldenburg, 2002. An integral-equation solution to the geophysical electromagnetic forward-modelling problem, in Three-Dimensional Electromagnetics: Proceedings of the Second International Symposium, M.S. Zhdanov and P.E. Wannamaker (eds).

Back to top.

Program EM1DTM

This program is designed specifically to handle lots and lots of time-domain EM soundings, which is what results from airborne surveys. The voltage (or magnetic field) measurements for the set of time channels for a particular spatial location are inverted to give a horizontally layered conductivity model. All the models can then be pasted side-by-side to give an image of the subsurface (see picture below).

Data from an airborne time-domain EM survey (top panel; dots indicate observations, lines show the predicted data), and the 2-D conductivity image of the subsurface constructed from the 1-D inversions of the soundings.

EM1DTM is the sister program to EM1DFM, and has many of the same features (under-determined, minimum structure inversions; four possible ways of determining the trade-off parameter: fixed or cooling schedule, discrepancy principle, GCV, L-curve), doesn't have others (just does conductivity, not susceptibility), and has some features EM1DFM doesn't (general measures: for constructing blocky, piecewise-constant models and using robust measures of data misfit).

About 3km's worth of the conductivity image produced by EM1DTM from a long line of GEOTEM data acquired during exploration for oil sands in Alberta. An L1-type norm was used for the measure of model structure, which produces this character of model with relatively sharp interfaces between units of different conductivities. Jamin Cristall did the inversions as part of his BASc thesis at UBC. (More information here.)

References:

Cristall, J., C.G. Farquharson and D.W. Oldenburg, Inversion of Airborne Electromagnetic Data: Application to Oil Sands Exploration, The 2004 Joint Assembly of the CGU, AGU, SEG and EEGS, Montréal, PQ, 17-21 May 2004.

Cristall, J., C.G. Farquharson and D.W. Oldenburg, Airborne Electromagnetic Inversion Applied to Oil Sands Exploration, 2004 CSEG National Convention, Calgary, AB, 10-13 May 2004.

Also, mainly for EM1DFM:

Farquharson, C.G., and D.W. Oldenburg, 2004. A comparison of automatic techniques for estimating the regularization parameter in nonlinear inverse problems, Geophysical Journal International, 156, 411-425.

Farquharson, C.G., D.W. Oldenburg and P.S. Routh, 2003. Simultaneous one-dimensional inversion of loop-loop electromagnetic data for magnetic susceptibility and electrical conductivity, Geophysics, 68, 1857-1869.

Back to top.

Return to my home page.

Last update: 4th February 2012.