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.

3-D EM forward modelling for real-life Earth models

We have now progressed to running our 3-D EM forward modelling codes rather routinely on real-life Earth models. These Earth models are parameterized in terms of unstructured tetraherdal meshes (of course!), which are constructed from wireframe surfaces defining the geological contacts. We are doing this using the finite-element and finite-volume codes of my now post-docs Masoud Ansari and Hormoz Jahandari, as well as the finite-element code of Jianhui Li (who we were lucky enough to have visit us last year).

The top figure below shows a zoomed-in view of the model of the massive sulphide Ovoid ore deposit at Voisey's Bay, Labrador, and the refinement indication the source and receiver locations for a DIGHEM survey that we modelled. The figure on the bottom left shows the real (lines) and synthesized (crosses) DIGHEM data along this survey line over the Ovoid. The figure on the bottom right shows the strength of the current density within the Ovoid for the transmitter located above the centre of the Ovoid.

We have also synthesized marine CSEM data for a number of scenarios off-shore Newfoundland. For more information, and pictures, for these examples, please see the abstracts and links below.

References:

LI Jian-Hui, Colin G. Farquharson, HU Xiang-Yun, ZENG Si-Hong, 2016. A vector finite element solver of three-dimensional modelling for a long grounded wire source based on total electric field, Chinese Journal of Geophysics, 59(4): 1521-1534. (DOI: 10.6038/cjg20160432)

Surface-based inversion

Peter Lelièvre, with me trying to keep up, is developing methodology and software that inverts for the locations of surfaces in an Earth model, i.e., a surface-based inversion rather than the current usual pixellated, cell-based inversion. This is motivated by working with geolgical Earth models parameterized in terms of tessellated wireframe surfaces -- why can't we just move those surfaces around and morph those surfaces in our geophysical inversion -- and a continuing search for an effective way of constructing piecewise-constant models.

The top two figures below show the result of inverting a real data-set for an amorphous blob that could be an ore body. The green nodes and edges indicate the control frame, and the coloured wire-frame mesh (with black outlines for the triangular facets) indicates the actual surface constructed by the inversion. The colouring indicates the uncertainty in the locations of the varioius parts of the surface: colder colours corrsepond to smaller uncertainty. The bottom two images are for a synthetic example in which a surface between two units of different density is reconstructed during the inversion.

References:

Lelièvre, P.G., and C.G. Farquharson, Inversion of Potential Field Data for Contact Surface Geometry, EAGE Conference and Exhibition, Vienna, 30 May-2 June, 2016.

3-D EM forward modelling on unstructured meshes

We have developed methodology and software for synthesizing controlled-source EM data for 3D Earth models specified in terms of unstructured tetrahedral meshes. In particular, Masoud Ansari has derived, implemented and coded-up a finite-element solution for vector and scalar potentials, and Hormoz Jahandari has derived, implemented and coded a finite-volume solution that amounts to a staggered grid scheme for tetrahedral meshes. Both Masoud's and Hormoz's formulations are for total fields (so that we can handle general topography) and for the frequency domain.

Different discretizations of the subsurface, continued; computational geophysics on unstructured grids

Unstructured tetrahedral grids are an obvious means of parameterizing Earth models because of the flexibility they offer. However, the numerical methods needed to compute synthetic geophysical data for such grids are more complicated, sometimes significantly so, than the methods needed to compute things using rectilinear grids: more challenging, more intellectually intriguing, and ultimately satisfying. Dr. Peter Lelièvre, a research scientist in my group, has written a seismic first-arrival travel-time solver for unstructured tetrahedral grids. This uses a fast marching algorithm derived for unstructured tetrahedral grids. Peter has also written a travel-time inversion code, and has developed joint inversion methodology and software - all for tetrahedral unstructured grids - for gravity and travel-time data (see next section). Angela Carter-McAuslan, an MSc student with me, has been applying Peter's joint inversion code to data-sets generated from realistic models of the Voisey's Bay deposits. Cassandra Tycholiz, also an MSc student with me, is using a gravity gradiometry variant of Peter's code and is investigating its capabilities, and that of airborne gradiometry data, on realistic synthetic data-sets from Voisey's Bay. A PhD student with me, Hormoz Jahandari, has written finite-element and finite-volume codes for computing gravity data for unstructured grid Earth models. (Finite element and finite volume methods for gravity because they can be used for memory efficient implicit computation of the Jacobian matrix in inversions, just as is done for MT.) Hormoz is now turning his attentions to EM data. Another PhD student with me, Seyedmasoud Ansari, is also investigating forward modelling methods for EM data for unstructured grids, in particular, finite elements and vector-scalar potential decompositions of the electric field. All very interesting and stimulating, and hopefully useful. There are more details about the current state of the work mentioned above in the publications and presentations listed at the end of this section.

Below are a couple of views of three travel-time iso-surfaces for first-arrival travel-time computations through a model of the Ovoid massive sulphide ore body at Voisey's Bay. The sulphide is significantly slower than the host rock. The model was specified using an unstructured tetrahedral grid that filled the volumes within and outside the wireframe surface of the massive sulphide unit that had been derived from drilling. (The troctolite was not included in this model.) The computations were done using an implementation of a fast marching method derived specifically for unstructured tetrahedral grids.

References:

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:

Carter-McAuslan, A., P.G. Lelièvre, and C.G. Farquharson, 2015. A study of fuzzy c-means coupling for joint inversion, using seismic tomography and gravity data test scenarios, Geophysics, 80, W1-W15. (DOI: 10.1190/geo2014-0056.1)

Carter-McAuslan, A., P.G. Lelièvre, and C.G. Farquharson, Joint inversion of seismic traveltime and gravity data: A synthetic study using geologically realistic models from the Voisey's Bay deposit, AGU Fall Meeting, San Francisco, 3-7 December, 2012.

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.

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.

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.

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.

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.