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.

**Surface Geometry Inversion (SGI):**
We are continuing to develop our "surface geometry inversion" (SGI) approach. (Earlier work is also mentioned below.)
This approach works with a wireframe model of the Earth's subsurface.
The wireframe surfacess (typically surfaces made up of tessellations of triangular facets, as that's the most flexible
and general) are the contacts between the different geological units.
These wireframe surfaces are moved around during the inversion process.
We are using Genetic Algoritms as the global optimization method, with a follow-on Markov chain Monte Carlo
sampling procedure to get uncertainties in the optimal locations of the wireframe interfaces.
So far Peter and Chris Galley have implemented and used SGI for gravity and magnetics, and Xushan Lu is implementing
SGI for time-domain EM and DC resistivity.
More on this below.

**More joint inversion, and clustering inversion:**
Mehrdad Darijani, for his PhD project, applied Peter's joint and clustering inversion codes (gravity, magnetics,
seismic travel-time) to data from the Athabasca Basin as part of the CMIC "Footprints" project.
The challenge was to determine the thickness of the quaternary sediments, and to take this into account when modelling
and inverting the gravity data, with the ultimate aim being to get information about possible alteration zones
within the sandstones (which could be indicators of uranium deposits).
The variable thickness of the quaternary sediments nicely obfuscates any gravity signal coming from alteration deeper
in the sandstones.
Mehrdad tried various combinations of data---seismic refraction, magnetic, gravity---and various inversion
approaches---gravity inversion constrained by results of seismic inversion, joint gravity and magnetic inversion, joint
gravity and seismics, regular sum-of-squares and clustering---to see if any alteration in the sanstones could be made out.
More on this below.

**Meshfree forward modelling:**
Unstructured tetrahedral meshes are very flexible and powerful for discretizing volumes between arbitrary, complicated
surfaces (topography!).
However, for certain forward-modelling problems such as finite-element solutions to the EM equations, the "quality" of the
mesh, i.e., how *not* long and pointy the tetrahedral cells are, is important.
However, for models with lots of complicated surfaces, it can be hard to generate a high quality tetrahedral mesh.
Wouldn't it be great if we could have one mesh (unstructured tetrahedral mesh) for parameterizing the Earth model, and
a separate, distinct mesh (or simply a cloud of nodes!) for the numerical solution of the partial differential equation?
Jianbo Long investigated this for his PhD, implementing meshfree meshes for gravity and for EM.
More on this below.

**3-D EM forward modelling for real-life Earth models:**
We have successfully applied our various 3-D EM forward-modelling codes to a number of real-life
Earth models, including high conductivity contrast mineral exploration scenarios and marine CSEM
hydrocarbon exploration scenarios.
Between my two post-docs, Masoud Ansari and Hormoz Jahandari, and Dr. Jianhui Li who visited us last
year, we have finite-element and finite-volume codes both for E-field and vector-scalar potential
formulations.
More on this below.

**Surface-based inversion:**
We (well Peter Lelièvre is!) are investigating inversion methods that can work with the surfaces
in a wireframe Earth model.
The goal is to have the geophysical inversion work on exactly the same Earth model that the geologists
are working with.
One approach is to parameterize the surfaces using a small number of control nodes with subsequent
surface subdivision, and then to invert for the locations of these nodes using a global optimization
technique.
For a little more on this please see below.

**3-D EM forward modelling on unstructured meshes:**
Two of my PhD students, Masoud Ansari and Hormoz Jahandari, have developed 3-D controlled-source
EM forward modelling methods for Earth models parameterized in terms of unstructured tetrahedral
meshes.
Masoud has implemented a finite-element approach with the electric field decomposed into
vector and scalar potentials.
This enables the relative importance of the inductive and galvanic contributions to the electromagnetic
fields generated in the Earth to be investigated.
Hormoz has developed a finite-volume approach that essentially corresponds to a staggered grid
formulation for tetrahedral meshes.
Both Masoud and Hormoz use total field formulations so that sources on arbitrary topography can
be considered.
More on this below.

**Different discretizations of the subsurface, continued; computational geophysics
on unstructured grids:**
We have been working on computational methods for unstructured grids for a couple of years
now. (Initial thoughts are preserved below.)
We've been working through most of the data types one would want.
So far we've developed code for modelling gravity, gravity gradient, magnetic, seismic
travel-time, DC resistivity, and electromagnetic data.
This is using finite element and finite volume methods on tetrahedral and Voronoï grids for
the potential field methods and EM, and an implementation of a fast marching method on tetrahedral
grids for seismic travel times.
We currently have inversion code for gravity, gravity gradiometry, magnetic and seismic travel-time
data (including joint inversion - see the next topic); inversion code for EM may follow.
A few more specifics on all this below.

**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.

As an alternative to fuzzy blob minimum-structure inversions (a second, complementary tool, not a replacement tool), and as an alternative to attempting to generate blocky, piecewise-constant models using non-L2 measures in minimum-structure inversion, we're investigating what we're calling "surface geometry inversion" (SGI). In this approach, the Earth model is parameteried in terms of wirefrace surfaces (tessellations of triangular facets) between otherwise uniform volumes. These are the geological contacts between the different geological units. We then use a global optimization process (Genetic Algorithms at the moment) to minimize the data misfit (the only thing being minimized in this form of the inverse problem). We also have a go with a Markov chain Monte Carlo sampling process to get estimates of the uncertainties on the locations of the wireframe surfaces in the best-fitting model.

The figures below show the results of Chris Galley's work applying the SGI approach to magnetic and gravity data from over the active mound at the Trans-Atlantic Geotraverse (TAG) hydrothermal vent field on the Mid-Atlantic Ridge. The left panels show the starting model and final model for the magnetic inversion. The orange-y colours show the uncertainty of the constructed surface. The images on the right show Chris' new and improved geological model based on the magnetic and gravity SGI results. (The vertical lines indicate the holes that have been drilled into this feature.)

Xushan Lu has also been implementing SGI for time-domain EM data, in the context of finding and delineating graphitic fault zones in the Athabasca Basin, which can be tell-tale indicators of uranium deposits in the area.

References:

Galley, C., P. Lelièvre, A. Haroon, S. Graber, J. Jamieson, F. Szitkar, I. Yeo, C. Farquharson, S. Petersen, R. Evans, 2021, Magnetic and gravity surface geometry inverse modeling of the TAG active mound, Journal of Geophysical Research: Solid Earth, 126, e2021JB022228. (10.1029/2021JB022228)

Lu, X., C.G. Farquharson, P. Lelièvre, G. Harrison and J.-M. Miehé, Surface geometry inversion of Preston Lake transient electromagnetic data, First International Meeting for Applied Geoscience and Energy (i.e., Annual Meeting of SEG, joint with AAPG), Denver, 26 Sept - 1 Oct, 2021.

Galley, C.G., J.W. Jamieson, P.G. Lelièvre, C.G. Farquharson and J.M. Parianos, 2020, Magnetic imaging of subseafloor hydrothermal fluid circulation pathways, Science Advances, 6 (44), eabc6844. (DOI:10.1126/sciadv.abc6844)

Galley, C.G., P.G. Lelièvre and C.G. Farquharson, 2020, Geophysical inversion for 3D contact surface geometry, Geophysics, 85, K27-K45. (DOI:10.1190/geo2019-0614.1)

Lu, X., C. Galley, P.G. Lelièvre and Colin Farquharson, Surface geometry inversion of potential field and electromagnetic geophysical data, AGU Fall Meeting, San Francisco, 1-17 December, 2020.

Galley, C., J.W. Jamieson, P.G.Lelièvre, C. Farquharson and J. Parianos, Imaging Subseafloor Hydrothermal Systems, AGU Fall Meeting, San Francisco, 1-17 December, 2020.

Lu, X., P. Lelièvre and C.G. Farquharson, Surface geometry inversion of time-domain EM data, SEG Annual Meeting, Houston, 11-16 October, 2020.

Galley, C., P. Lelièvre, C. Farquharson, J. Jamieson, A. Haroon, S. Graber, S. Petersen and F. Szitkar, Modelling the geometry of the Trans-Atlantic Geotraverse seafloor massive sulphide deposit through magnetic surface geometry inversion, SEG Annual Meeting, Houston, 11-16 October, 2020.

The figures below show the results of constrained inversion of real gravity data from the Athabasca Basin, with the overburden layer derived from seismic refraction inversion used as the main constraint in the gravity inversion. The top two panels show the refraction inversion result and the fits to the observed travel-time data. This was a clustering inversion. The readily identified base of overburden in this clustering inversion result was used to define an overburden layer of uniform density that was then used to provide bounds on the densities in the upper parts of the density model during the gravity inversion. The lower panels show the gravity inversion results. The black outline indicates a zone of suspected alteration, wich would be responsible for a density low, which is indeed what is apparent in this model constructed by the constrained seismic-then-gravity inversion.

Mehrdad also did joint gravity and magnetic inversion, using a fuzzy-clustering, physical-property-based joint inversion approach. The image up with the summary teaser paragraph is from the joint inversion of real gravity and magnetic data from the Athabasca Basin. This joint inversion was also heavily constrained as we know the sandstones are effectively non-magnetic, with the only magnetic units being the glacial sediments (sort of) and the crystalline basement.

References:

Darijani, M., C.G. Farquharson and P.G. Lelièvre, 2021, Joint and constrained inversion of magnetic and gravity data: A case history from the McArthur River area, Canada, Geophysics, 86, B79-B25. (DOI:10.1190/GEO2019-0818.1)

Darijani, M., C.G. Farquharson and P.G. Lelièvre, 2020, Clustering and constrained inversion of seismic refraction and gravity data for overburden stripping: Application to uranium exploration in the Athabasca Basin, Canada, Geophysics, 85, B133-B146. (DOI: 10.1190/geo2019-0525.1)

Lelièvre, P., M. Darijani, C. Galley, P. Zheglova and C. Farquharson, No Magic Bullet: Three Integrated Imaging Problems, Three Solutions, IUGG General Assembly, Montréal, 8-18 July, 2019.

Darijani, M., and C.G. Farquharson, Synthetic modeling and joint inversion of gravity and seismic refraction data for overburden stripping in the Athabasca Basin, Canada, SEG Annual Meeting, Dallas, 16-21 October 2016.

For finite-element and finite-volume solutions of the geophysical EM forward-modelling problem, the quality of the mesh used to discretize the computational domain can impact the efficiency and effectiveness of the solution process. This is not so much the case if a direct solver is used to solve the system of equations resulting from the finite element or finite volume discretization, but can be very much the case if an iterative solver is used. Generating a quality mesh can be challenging for complicated (real-life) Earth models. Meshfree approaches for solving partial differential equations only need a cloud of points, not a celluralized mesh. Jianbo Long has investigated thes meshfree approaches, trying them out for first gravity (solving Poisson's equation) then for MT and controlled-source EM.

The figures below show Jianbo's MT results for the COMMEMI 3D-1A model. The top row of images shows the model, an unstructured tetrahedral mesh discretization of the model, and the cloud of nodes used for the meshfree forward-modelling approach. The connectivity between the nodes in the tetrahedral mesh (i.e., the edges connecting nodes into tetrahedra) introduces an extra level of compexity into a cellularized, voxellated mesh that's not present in the cloud of nodes used in a meshfree approach. The lower panels are the MT results along two orthogonal profiles crossing over the conductive block. Jianbo's meshfree results are compared with the COMMEMI results and results computed by Hormoz.

See also Xin Huang's work on the spectral element method in which the discretization for the spectral element basis functions and numerical quadrature is decoupled from the mesh describing the Earth model.

References:

Huang, X., C.G. Farquharson, C. Yin, L. Yan, X. Cao and B. Zhang, 2021, A 3D forward-modeling approach for airborne electromagnetic data using a modified spectral-element method, Geophysics, 86, E343-E354. (DOI:10.1190/GEO2020-0004.1)

Long, J., and C.G. Farquharson, Meshfree modelling of 2D MT data with RBF-FD and unstructured point, SEG Annual Meeting, Houston, 11-16 October, 2020.

Long, J., and C.G. Farquharson, 2019, On the forward modelling of three-dimensional magnetotelluric data using a radial-basis-function-based mesh-free method, Geophysical Journal International, 219, 394-416. (DOI: 10.1093/gji/ggz306)

Long, J., and C.G. Farquharson, 2019, Three-dimensional forward modelling of gravity data using mesh-free methods with radial basis functions and unstructured nodes, Geophysical Journal International, 217, 1577-1601. (DOI: 10.1093/gji/ggz115)

Long, J., and C.G. Farquharson, Meshfree modelling of 3-D controlled-source EM data: A new method to treat the singular source terms, SEG Annual Meeting, San Antonio, 15-20 September, 2019.

Long, J., and C.G. Farquharson, Three-dimensional magnetotelluric modelling with an RBF-based meshfree method and unstructured points, 24th EM Induction Workshop, Helsingør, Denmark, 13-20 August 2018.

Farquharson, C., J. Long, X. Lu, and P.G. Lelièvre, Electromagnetic forward modelling for realistic Earth models using unstructured tetrahedral meshes and a meshfree approach, AGU Fall Meeting, New Orleans, 11-15 December 2017.

Long, J., and C. Farquharson, Three-dimensional controlled-source EM modeling with radial basis function-generated finite differences: A meshless approach, SEG Annual Meeting, Houston, 24-29 September 2017.

Long, J., and C.G. Farquharson, Geophysical electromagnetic data modelling with radial basis function generated finite differences, 6th International Symposium in Three-Dimensional Electromagnetics, Berkeley, 28-30 March 2017.

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)

Dunham, M.W., S. Ansari, and C.G. Farquharson, Three-dimensional computer modeling of realistic marine CSEM earth models in the Flemish Pass Basin, AAPG Convention & Exhibition, Calgary, 19-22 June, 2016. [PDF file of poster (3.2M).]

Dunham, M.W., S. Ansari, and C.G. Farquharson, Three-dimensional Marine CSEM Forward Modelling in the Flemish Pass Basin Using Realistic Unstructured Meshes, EAGE Conference and Exhibition, Vienna, 30 May-2 June, 2016.

Nalepa, M., S. Ansari, and C.G. Farquharson, Synthetic modelling of marine controlled-source electromagnetic data for hydrocarbon exploration, GeoConvention, Calgary, 7-11 March, 2016. [PDF file of poster (9M).]

Jones, D., S. Ansari and C.G. Farquharson, Synthesizing 3D Time-Domain Electromagnetic Geophysical Forward Models for Uranium Exploration in the Athabasca Basin, PDAC-SEG Student Minerals Colloquium, Toronto, 5 March 2016. [PDF file of poster (2M).]

Jahandari, H., and C.G. Farquharson, 2015. Finite-volume modelling of geophysical electromagnetic data on unstructured grids using potentials, Geophysical Journal International, 202, 1859-1876. (DOI: 10.1093/gji/ggv257)

Ansari, S., C.G. Farquharson, and S. MacLachlan, Gauged vector finite-element schemes for the geophysical EM problem for unique potentials and fields, SEG Annual Meeting, New Orleans, 18-23 October 2015. [PDF file of talk (71M).]

Jahandari, H., and C.G. Farquharson, Finite-volume modeling of geophysical electromagnetic data using potentials on unstructured staggered grids, SEG Annual Meeting, New Orleans, 18-23 October 2015. [PDF file of talk (9M).]

Back to top.

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.

Lelièvre, P.G., R. Bijani, and C.G. Farquharson, Lithological and Surface Geometry Joint Inversions Using Multi-Objective Global Optimization Methods, EGU General Assembly, Vienna, 17-22 April, 2016. [PDF file of presentation (3.3M).]

Lelièvre, P.G., R. Bijani, and C.G. Farquharson, Geophysical Inversion With Multi-Objective Global Optimization Methods, EGU General Assembly, Vienna, 17-22 April, 2016.

Lelièvre, P.G., C.G. Farquharson, and R. Bijani, 3D potential field inversion for wireframe surface geometry, SEG Annual Meeting, New Orleans, 18-23 October 2015. [PDF file of talk (7M).]

Lelièvre, P., C. Farquharson, and R. Bijani, 3D stochastic geophysical inversion for contact surface geometry, EGU General Assembly, Vienna, 12-17 April 2015. [PDF file of presentation (10M).]

Back to top.

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.

The figures below show a comparison of the fields computed by Masoud's and Hormoz's codes for a horizontal coplanar source-receiver configuration "flying" over a graphite cube immersed in brine. This is Ken Duckworth's physical scale modelling example that I used to test my integral equation code on a number of years ago. (See below for more information on this.) The grey shaded diagram below (to scale) shows the geometry of the situation. The conductivity of the graphite cube was 63,000S/m and the conductivity of the brine was 7.3S/m! The coloured panels on the left show vertical sections through the Delauny tetrahedral grid and its dual Voronoï grid. There is lots of refinement in the mesh along where the transmitters and receivers are located, and a bit of refinement in and around the cube. The ability of the unstructured meshes to be refined considerably in certain locations, to transition into regions of coarse discretization, and to do both of these things while maintaining a quality discretization, is very useful for EM numerical modelling. The panels on the right show the real and imaginary parts of the secondary H-field as the transmitter-receiver pair moves over the cube for a number of different frequencies. The circles indicate the values computed by Hormoz, the black lines indicate the values computed by Masoud, the red lines show the physical scale modelling measurements, and the gold lines show the integral equation results. We're definitely satisfied with this comparison.

References:

Jahandari, H., and C.G. Farquharson, 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids, Geophysics, 79, E287-E302. (DOI: 10.1190/geo2013-0312.1)

Ansari, S., and C.G. Farquharson, 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids, Geophysics, 79, E149-E165. (DOI: 10.1190/geo2013-0172.1)

Ansari, S., and C.G. Farquharson, A Potential Method for Three-dimensional Numerical Modeling of Geophysical Electromagnetic Problems, EAGE Conference and Exhibition, Amsterdam, 16-19 June, 2014. [PDF file of talk (75M).]

Jahandari, H., and C.G. Farquharson, Forward Modelling of Geophysical Electromagnetic Data on Unstructured Grids Using a Finite-volume Approach, EAGE Conference and Exhibition, Amsterdam, 16-19 June, 2014. [PDF file of talk (16M).]

Farquharson, C.G., P.G. Lelièvre, S. Ansari, and H. Jahandari, Towards Real Earth Models -- Computational Geophysics on Unstructured Tetrahedral Meshes?, EAGE Conference and Exhibition, Amsterdam, 16-19 June, 2014. [PDF file of talk (18M).]

Jahandari, H., and C.G. Farquharson, A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids, SEG Annual Meeting, Houston, 22-27 September 2013. [PDF file of talk (15M).]

Ansari, S., and C.G. Farquharson, Numerical modeling of geophysical electromagnetic inductive and galvanic phenomena, SEG Annual Meeting, Houston, 22-27 September 2013. [PDF file of talk (44M).]

Jahandari, H., and C.G. Farquharson, Finite volume modelling of electromagnetic data using unstructured staggered grids, 3DEM-5 5th International Symposium on Three-Dimensional Electromagnetics, Sapporo, 7-9 May 2013. [PDF file of talk (11M).]

Ansari, S., and C.G. Farquharson, Three dimensional modeling of controlled-source electromagnetic response for inductive and galvanic components, 3DEM-5 5th International Symposium on Three-Dimensional Electromagnetics, Sapporo, 7-9 May 2013. [PDF file of poster (29M).]

Back to top.

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:

Farquharson, C.G., P.G. Lelièvre, S. Ansari, and H. Jahandari, Towards Real Earth Models -- Computational Geophysics on Unstructured Tetrahedral Meshes?, EAGE Conference and Exhibition, Amsterdam, 16-19 June, 2014. [PDF file of talk (18M).]

Lelièvre, P.G., and C.G. Farquharson, 2013. Gradient and smoothness regularization operators for geophysical inversion on unstructured meshes, Geophysical Journal International, 195, 330-341. (DOI: 10.1093/gji/ggt255)

Jahandari, H., and C.G. Farquharson, 2013. Forward modeling of gravity data using finite-volume and finite-element methods on unstructured grids, Geophysics, 78, G69-G80. (DOI:10.1190/geo2012-0246.1)

Lelièvre, P.G., C.J. Tycholiz, and C.G. Farquharson, Constrained geophysical inversion on unstructured meshes, SEG Annual Meeting, Las Vegas, 4-9 November, 2012.

Lelièvre, P.G., A. Carter-McAuslan, C. Tycholiz, C.G. Farquharson, and C.A. Hurich, Unstructured grid modelling to create 3D Earth models that unify geological and geophysical information, GAC-MAC Joint Annual Meeting, St. John's, 27-29 May 2012. [PDF file of talk (2.7M).]

Tycholiz, C., P.G. Lelièvre, and C.G. Farquharson, Geologically constrained inversion of airborne gravity gradiometer data using unstructured tetrahedral meshes, GAC-MAC Joint Annual Meeting, St. John's, 27-29 May 2012.

Carter-McAuslan, A., P.G. Lelièvre, G. Blades, and C.G. Farquharson, Testing joint inversion code with geologically realistic models, GAC-MAC Joint Annual Meeting, St. John's, 27-29 May 2012. [PDF file of talk (8.5M).]

Lelièvre, P.G., A. Carter-McAuslan, C. Tycholiz, C.G. Farquharson, and C.A. Hurich, Refining 3D Earth models through constrained joint inversion on flexible unstructured meshes, KEGS-DMEC Geophysical Symposium, Toronto, 3 March, 2012. [PDF file of talk (4.2M).]

Lelièvre, P.G., A. Carter-McAuslan, C. Farquharson, and C. Hurich, 2012. Unified geophysical and geological 3D Earth models, The Leading Edge, 31, 322-328.

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-K15.

Lelièvre, P.G., C.G. Farquharson, and C.A. Hurich, 2011. Inversion of first-arrival seismic traveltimes without rays, implemented on unstructured grids, Geophysical Journal International, 185, 749-763.

Lelièvre, P.G., C.G. Farquharson, and C.A. Hurich, 2011. Computing first-arrival seismic traveltimes on unstructured 3-D tetrahedral grids using the Fast Marching Method, Geophysical Journal International, 184, 885-896.

Ansari, S., and C.G. Farquharson, 3D finite-element simulation of electromagnetic data for inductive and galvanic components, SEG Annual Meeting, San Antonio, 18-23 September, 2011.

Jahandari, H., and C.G. Farquharson, Forward modelling of gravity data for unstructured grids using the finite-volume method, SEG Annual Meeting, San Antonio, 18-23 September, 2011.

Back to top.

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.

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.

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.

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.

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.

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).

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: 2nd January 2022.