3D Modelling


Highly Conductive Massive Sulphide Bodies

Lamontagne Geophysics Ltd.

Owen Fernley


3D Forward Modelling Tool

  • Large loop and large volume EM problems
  • High contrast resistivity
  • Surface and downhole survey modelling
  • Ideal for massive sulphide exploration


Finite Difference 3D EM Modelling

  • Finite Difference Equation
  • The Staggered Grid
  • Boundary Conditions
  • Multigrid Solver


Results and Case Studies

  • Synthetic Sphere
  • Voisey's Bay Ovoid
  • Mystery Project

Part 1

Finite Difference 3D EM Modelling

Equation to Solve

$ (\vec{\nabla} \times \color{orange}{\rho} \vec{\nabla} \times \color{fuchsia}{\vec{\mathrm{H}}}) + \frac{\partial\mu_0\color{fuchsia}{\vec{\mathrm{H}}}}{\partial t} = -\frac{\partial\vec{\mathrm{B}}_p}{\partial t} $

Time Stepping

$ \vec{\mathrm{H}}_i = \vec{\mathrm{H}}_{i-1} + \vec{\mathrm{D}} $

Maxwell's Equations

$ \vec{\nabla} \bullet \vec{\mathrm{E}} = \frac{\rho}{\varepsilon_0} $

$ \vec{\nabla} \bullet \vec{\mathrm{B}} = 0 $

$ \vec{\nabla} \times \vec{\mathrm{E}} = -\frac{\partial\vec{\mathrm{B}}}{\partial t} $

$ \vec{\nabla} \times \vec{\mathrm{B}} = \mu_0 \left(\vec{\mathrm{J}} + \varepsilon_0\frac{\partial\vec{\mathrm{E}}}{\partial t}\right) $

Current Density

$ \color{red}{\vec{\mathrm{E}}} = \color{orange}{\rho}\vec{\mathrm{J}} $

Faraday's Law

$ \vec{\nabla} \times \color{red}{\vec{\mathrm{E}}} = -\frac{\partial\vec{\mathrm{B}}}{\partial t} $

Ampere's Law

$ \vec{\nabla} \times \vec{\mathrm{B}} = \mu_0 \vec{\mathrm{J}} $

Current Density

$ \color{grey}{\vec{\mathrm{E}} = \rho\vec{\mathrm{J}}} $

Faraday's Law

$ \vec{\nabla} \times \color{orange}{\rho}\vec{\mathrm{J}} = -\frac{\partial\vec{\mathrm{B}}}{\partial t} $

Ampere's Law

$ \vec{\nabla} \times \vec{\mathrm{B}} = \mu_0 \vec{\mathrm{J}} $

Current Density

$ \color{grey}{\vec{\mathrm{E}} = \rho\vec{\mathrm{J}}} $

Faraday's Law

$ \vec{\nabla} \times \color{orange}{\rho}\vec{\mathrm{J}} = -\frac{\partial\vec{\mathrm{B}}}{\partial t} $

Ampere's Law

$ \vec{\nabla} \times \color{fuchsia}{\vec{\mathrm{B}}} = \color{fuchsia}{\mu_0} \vec{\mathrm{J}} $

Current Density

$ \color{grey}{\vec{\mathrm{E}} = \rho\vec{\mathrm{J}}} $

Faraday's Law

$ \vec{\nabla} \times \color{orange}{\rho}\color{green}{\vec{\mathrm{J}}} = -\frac{\partial\vec{\mathrm{B}}}{\partial t} $

Ampere's Law

$ \vec{\nabla} \times \color{fuchsia}{\vec{\mathrm{H}}} = \color{green}{\vec{\mathrm{J}}} $

Current Density

$ \color{grey}{\vec{\mathrm{E}} = \rho\vec{\mathrm{J}}} $

Faraday's Law

$ \vec{\nabla} \times \color{orange}{\rho} \vec{\nabla} \times \color{fuchsia}{\vec{\mathrm{H}}} = -\frac{\partial\vec{\mathrm{B}}}{\partial t} $

Ampere's Law

$ \color{grey}{\vec{\nabla} \times \vec{\mathrm{H}} = \vec{\mathrm{J}}} $

Ampere's Law $\longrightarrow$ Faraday's Law

$ \vec{\nabla} \times \color{orange}{\rho} (\vec{\nabla} \times \color{fuchsia}{\vec{\mathrm{H}}}) = -\frac{\partial\vec{\mathrm{B}}}{\partial t} $

Secondary Primary Field

$ \vec{\nabla} \times \color{orange}{\rho} (\vec{\nabla} \times \color{fuchsia}{\vec{\mathrm{H}}}) = \left( \frac{\partial\mu_0\color{fuchsia}{\vec{\mathrm{H}}}}{\partial t} - \frac{\partial\vec{\mathrm{B}}_p}{\partial t} \right) $

Equation to Solve

$ (\vec{\nabla} \times \color{orange}{\rho} \vec{\nabla} \times \color{fuchsia}{\vec{\mathrm{H}}}) + \frac{\partial\mu_0\color{fuchsia}{\vec{\mathrm{H}}}}{\partial t} = -\frac{\partial\vec{\mathrm{B}}_p}{\partial t} $

Time Stepping

$ \color{fuchsia}{\vec{\mathrm{H}}}_i = \color{fuchsia}{\vec{\mathrm{H}}}_{i-1} + \color{fuchsia}{\vec{\mathrm{D}}} $

Finite Difference Equaton

$ \vec{\nabla} \times \color{orange}{\rho} \vec{\nabla} \times \vec{\alpha}\color{fuchsia}{\vec{\mathrm{D}}} + \frac{\mu_0\color{fuchsia}{\vec{\mathrm{D}}}}{h_t} = -\frac{\partial\vec{\mathrm{B}}_p}{\partial t} - \vec{\nabla} \times \color{orange}{\rho}\vec{\nabla} \times\color{fuchsia}{\vec{\mathrm{H}}} $

The Staggered Grid



Conductivity Subgrid

Conductivity Subgrid

Current Density

Current Density

The Staggered Grid

Boundary Conditions

Suitable Initial and Boundary Condition

  • $\vec{\mathrm{C}}$ and $\vec{\mathrm{J}}$ start at zero
  • $\vec{\mathrm{C}}$ steps on at $t=0$
  • $\vec{\mathrm{C}}$ and $\vec{\mathrm{J}}$ decay to zero at infinity

Zero Divergence: The $ \vec{\mathrm{C}} $ Vector

$ \vec{\mathrm{H}} = \vec{\mathrm{C}} - \vec{\nabla}\mathrm{g} $

$ \vec{\mathrm{H}} \longrightarrow \vec{\mathrm{C}} $ provided that $ \vec{\nabla} \bullet \vec{\mathrm{H}} = 0 $

Below Surface


C-Field Below Surface

in Free air

$ \vec{\mathrm{H}} = \vec{\nabla}\mathrm{g} $

at surface

$\vec{\mathrm{H}} = \vec{\mathrm{C}} - \vec{\nabla}\mathrm{g} $

below surface

$\vec{\mathrm{H}} = \vec{\mathrm{C}}$

Outer Edge

$\vec{\mathrm{H}_t} = 0$

Outer boundary can be pushed kilometers away

Multigrid Solver

  • Speeds up iteration time
  • Boundary conditions pushed far from area of interest
  • V and W cycle iterations: Interpolation and Sampling

Part 2

Results and Case Studies

Sphere Test (r = 100 m)


MGEM Software Components

  • Gridplot
  • 3D Viewer
  • MGEM Engine

Model Specifications

follow the law of geologic superposition

Grid Plot





sphere compare

3D Viewer

MGEM Engine


Four runs at once

Composing Models

  • background resistivity
  • cuboids
  • ellipsoids
  • polygonal prisms


  • dipping layered earths
  • quarter space
  • overburden
  • conductive or resistive dykes


ideal for lenses

Polygonal Prisms

detailing known bodies

All together now

Ovoid / OB / Hz
Ovoid / OB / Hz

Voisey's Bay Ovoid

Ovoid Photos
VB Section
VB field data

Iterative Forward Modelling

  • Ovoid Only
  • Overburden
  • Troctolyte Dykes
  • Conductive Conduit

Ovoid Image Source

VB field data

Gridplot Overlay

VB field data

Ovoid X-Section

VB field data

Ovoid Z-Section


Adding OverBurden


Troctolyte Dykes - Level 1


Troctolyte Dykes - Level 3

Level 3

Troctolyte Dykes - Level 6

Level 6

Adding a Conduit at Depth

Level 6 Conduit

Ovoid + OverBurden Hz

Ovoid / OB / Hz

Ovoid + OverBurden + Dykes Hz

Ovoid / OB / dykes / Hz

Ovoid + OverBurden + Dykes + Conduit Hz

Ovoid / OB / dykes / conduit / Hz
Voiseys Bay Compare Hz

Ovoid + OverBurden Hx

Ovoid / OB / Hz

Ovoid + OverBurden + Dykes Hx

Ovoid / OB / dykes / Hx

Ovoid + OverBurden + Dykes + Conduit Hx

Ovoid / OB / dykes / conduit / Hx
Voiseys Bay Compare HL

A More Realistic Exploration Example

plot B31A
plot B31A

Thank You!


Balch, S., 1999, Geophysical methods for nickel deposits with examples from Voisey’s Bay: GAC-MAC Meeting, St John’s Newfoundland.

Bengert, B., 2006, Successful application of electro-magnetics in the Reid Brook zone: targeting economic needles in a highly conductive haystack: ASEG Extended Abstracts 2006, 1–4.

Bochev, P. B., J. J. Hu, C. M. Siefert, and R. S. Tuminaro, 2008, An algebraic multigrid approach based on a compatible gauge reformulation of Maxwell’s equations: SIAM Journal on Scientific Computing, 31, 557–583.

Haber, E., D. Oldenburg, and R. Shektman, 2007, Inversion of time domain 3D electromagnetic data: Geophysical Journal International, 132, 1324–1335.

Hiptmair., R., 1997, Multigrid method for H(div) in three dimensions: Electronic Transactions on Numerical Analysis, 6, 133–152.

Hiptmair, R., 1999, Multigrid method for Maxwell’s equations: SIAM Journal on Numerical Analysis, 36, 204–225.

Hiptmair, R., and J. Xu, 2007, Nodal auxiliary space preconditioning in H(curl) and H(div) spaces: SIAM Journal on Numerical Analysis, 45, 2483–2509.

Jahandari, H., and C. G. Farquharson, 2014, A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids: Geophysics, 79, no. 6, E287-E302.

King, A., 2007, Review of geophysical technology for Ni-Cu-PGE deposits: Proceedings of the 5th Decennial International Conference on Mineral Exploration, 647–665.

Lamontagne, Y., and R. Langridge, 2013, UTEM ISR - induced source resistivity - Sierra Gorda project: ISR across a development-stage copper-molybdenum: KEGS PDAC Symposium 2013, 1-4.

Li, C., and A. J. Naldrett, 1999, Geology and petrology of the Voisey’s Bay intrusion: reaction of olivine with sulfide and silicate liquids: Lithos, 47, 1–31.

Lu, J. J., X. P. Wu, and K. Spitzer, 2010, Algebraic multigrid methods for 3D DC resistivity modeling: Chinese Journal of Geophysics, 53, 700–707.

Moucha, R., and R. C. Bailey, 2004, An accurate and robust multigrid algorithm for 2D forward resistivity modelling: Geophysical Prospecting, 52, 197-212.

Newman, G. A., 2014, A review of high-performance computational strategies for modeling and imaging of electromagnetic induction data: Surveys in Geophysics, 35, 85-100.

Polzer, B., 2000, The role of borehole EM in the discovery and definition of the Kelly Lake Ni- Cu deposit, Sudbury, Canada: 70th Annual International Meeting, SEG, Expanded Abstracts, 1063-1066.

Press, W. H., B. P. Flannery, S. A. Teukolksky, and W. T. Vetterling, 1988, Numerical recipes in C, the art of scientific computing: Cambridge University Press.

Yang, D., and D. W. Oldenburg, 2010, 3D forward modelling and inversion of inductive source resistivity data: 80th Annual International Meeting, SEG, Expanded Abstracts, 588-592.