Search
Close
Search
Advanced Search
Search Menu
AI Discovery Assistant

SUMMARY

We implement a spectral-infinite-element method (SIEM) to compute magnetic anomalies by solving a discretized form of the Poisson/Laplace equation. The SIEM combines the highly accurate spectral-element method with the mapped-infinite element method, which reproduces an unbounded domain accurately and efficiently. This combination is made possible by coupling Gauss–Legendre–Lobatto quadrature in spectral elements with Gauss–Radau quadrature in infinite elements along the infinite directions. Our method has two distinct advantages over traditional methods. First, the higher-order discretization accurately renders complex magnetized heterogeneities. Second, since the computation time is independent of the number of observation points, the method is efficient for very large models. We illustrate the accuracy and efficiency of our method by comparing calculated magnetic anomalies for various magnetized heterogeneities with corresponding analytical and commonly used computational solutions. We conclude with a practical example involving a complex 3-D model of an ore mine.

Magnetic anomalies: modelling and interpretation, Magnetic properties, Numerical modelling, Rock and mineral magnetism

1 INTRODUCTION

Analogous to gravity anomalies induced by density heterogeneities, magnetized objects produce magnetic anomalies, which can be measured on Earth’s surface or in space (Vine & Matthews 1963; Heirtzler et al.1968; Briais et al.1993). Similar to the gravitational potential, the magnetic potential is governed by Poisson’s equation inside a magnetization distribution and by Laplace’s equation elsewhere. Most existing methods for computing magnetic anomalies are based on a direct integral approach. Unlike mass density, magnetization is a vector quantity. Therefore, expressions for magnetic anomalies are often more complex than those for gravity anomalies. The procedure may be simplified by using the so-called Poisson’s relation, which takes advantage of analogies between magnetic and gravity fields to derive the magnetic field from the derivative of the gravitational potential (Blakely 1995). This method is often limited to simple isolated objects. There are several studies on computing magnetic anomalies based on the direct integral approach, for example, 2-D objects with a uniform magnetization (Talwani 1965; Won & Bevis 1987), infinitely thin laminae (Talwani 1965), laminae of finite thickness (Plouff 1976), rectangular prisms (Bhattacharyya 1964), polyhedrons (Bott 1963; Barnett 1976; Okabe 1979; Hansen & Wang 1988; Wang & Hansen 1990) and 3-D objects with arbitrary shapes (Talwani 1965). The software package MBOX is a widely used tool for computing magnetic anomalies due to uniform rectangular prisms (Plouff 1976; Blakely 1995). It can be used for bodies of arbitrary shape by dividing the body of interest into a number of small subprisms and superposing their individual contributions.

Although direct integral approaches are fast for small objects and a small number of observation points, the computational cost can quickly increase as the number of observation points grows. Furthermore, the volume integral approximation for complex objects or objects with variable properties may be inaccurate (Cai & Wang 2005; Martin et al.2017). The accuracy and efficiency of these methods may be significantly improved with numerical quadrature. For example, Martin et al. (2017) used spectral-elements to discretize the model and GLL quadrature for volume integration to calculate gravity anomalies.

We are not aware of any method for computing magnetic anomalies based on solving the weak form of Poisson’s equation. Cai & Wang (2005) used a finite-element method to compute gravity anomalies by imposing approximate boundary conditions on the outer surface of an extended model, under the assumption that the gravitational potential is of the order of r−4 at a distance r from the centre of mass. An alternative approach is to consider a very large domain that includes portions of outer space. This strategy requires large computational resources and is often inaccurate (Tsynkov 1998). A higher-order solver based on a convolution integral was also proposed to solve the unbounded Poisson equation (Hejlesen et al.2013). Alternatively, one may use spherical harmonics for spherically symmetric models (Dahlen 1974; Tromp & Mitrovica 1999).

In solid and fluid mechanics, infinite or semi-infinite domain problems are frequently solved using the displacement descent approach (Bettess 1977; Medina & Taylor 1983; El-Esnawy et al.1995) and the coordinate ascent approach (Beer & Meek 1981; Zienkiewicz et al.1983; Kumar 1985; Angelov 1991). In the displacement descent approach an element is mapped to a natural domain of interval [0, ∞]. This is realized by multiplying the standard interpolation functions by suitable decay functions. Since the integration interval is [0, ∞], classical Gauss–Legendre quadrature cannot be employed. Either Gauss–Legendre quadrature has to be modified to accommodate the [0, ∞] interval, or Gauss–Laguerre quadrature can be used (Mavriplis 1989). In terms of programming, both the Jacobian of the mapping and the numerical quadrature have to be modified from the classical finite-element method. In the coordinate ascent approach—also referred to as the ‘mapped infinite element’ approach—an element that extends to infinity in the physical domain is mapped to a standard natural element with interval [−1, 1]. This is achieved by defining shape functions using a pole—a point located outside the element opposite to the infinite direction. The corresponding shape function possesses a singularity at the far end of the infinite element. Unlike the displacement descent approach, standard Gauss–Legendre quadrature can be used, and only the Jacobian of the mapping needs to be modified from the classical finite-element method.

The spectral-element method (SEM) is a higher-order finite-element method that uses nodal quadrature, specifically, Gauss–Legendre–Lobatto (GLL) quadrature. For dynamic problems, the mass matrix is diagonal by construction. The method is highly accurate and efficient and is widely used in applications involving wave propagation (Faccioli et al.1997; Seriani & Oliveira 2008; Tromp et al.2008; Peter et al.2011), fluid dynamics (Patera 1984; Canuto et al.1988; Deville et al.2002), as well as for quasi-static problems (Gharti et al.2012a,b).

The spectral-infinite-element method (SIEM) combines the infinite-element approach based on coordinate ascent with the SEM. Gharti & Tromp (2017) originally developed the SIEM to calculate the background gravity field of the Earth, and subsequently applied the technique to several other problems, for example, gravity anomalies (Gharti et al.2018), coseismic and post-earthquake deformation (Gharti et al.2019a) and earthquake-induced gravity perturbations (Gharti et al.2019b). In this paper, we present an implementation of the SIEM for magnetic anomalies. The method is validated based on calculations of magnetic fields and total-field anomalies for a range of problems. Finally, we demonstrate a 3-D application of the method to an ore mine in Finland.

2 SPECTRAL-INFINITE-ELEMENT FORMULATION

2.1 Governing equation

Let Ω denote a magnetized region of interest with boundary Γ, as illustrated in Fig. 1. All of space is denoted by ◯, and the ‘boundary’ at infinity is denoted by Γ∞. The scalar magnetic potential, Φ, induced by a magnetization distribution, |$\mathbf {M}$|⁠, is governed by the Poisson’s equation
(1)
subject to the essential boundary conditions
(2)
and the natural boundary condition
(3)
Here |${}[\, \cdot \, ]_-^+$| denotes the jump in the enclosed quantity when going from the −ve side of boundary Γ to the +ve side, and |$\hat{\mathbf {n}}$| denotes the unit outward normal to the boundary, pointing from the −ve side to the +ve side. In eq. (1), the term |$\boldsymbol{\nabla }\cdot \mathbf {M}$| is defined throughout Ω , whereas the term |$\hat{\mathbf {n}}\cdot \mathbf {M}$| is restricted to the surface Γ . For a uniform magnetization, |$\boldsymbol{\nabla }\cdot \mathbf {M}$| vanishes, and only |$\hat{\mathbf {n}}\cdot \mathbf {M}$| contributes. Alternatively, the magnetization may also be expressed in terms of a volume and surface distribution of magnetostatic charges, Qv and Qs , respectively, such that
(4)
and
(5)
Outside the magnetization domain Ω the magnetization vanishes and the governing eq. (1) reduces to Laplace’s equation: ∇2Φ = 0 .
Schematic diagram of an unbounded domain with a fictitious boundary Γ∞ on which the magnetic potential vanishes, Φ = 0. The unbounded domain contains an unmagnetized region where $\mathbf {M}=\mathbf {0}$ , and a magnetized region, Ω, where $\mathbf {M}\ne \mathbf {0}$. The magnetized region Ω has a boundary Γ with unit outward normal $\hat{\mathbf {n}}$, pointing from the inside of the region (−) to the outside of the region (+).
Figure 1.

Schematic diagram of an unbounded domain with a fictitious boundary Γ∞ on which the magnetic potential vanishes, Φ = 0. The unbounded domain contains an unmagnetized region where |$\mathbf {M}=\mathbf {0}$| , and a magnetized region, Ω, where |$\mathbf {M}\ne \mathbf {0}$|⁠. The magnetized region Ω has a boundary Γ with unit outward normal |$\hat{\mathbf {n}}$|⁠, pointing from the inside of the region (−) to the outside of the region (+).

Open in new tabDownload slide

From hereon the approach is quite analogous to the approach of Gharti et al. (2018) for gravity anomalies problems, so we will be relatively brief, referring the reader to Gharti et al. (2018) for more detailed explanations of the method and its implementation.

2.2 Discretization

The weak form of the governing eq. (1) subject to the boundary conditions (2) and (3) is
(6)
where w denotes a test function and where we have used the fact that the magnetization vanishes outside of Ω.
The domain Ω is meshed using spectral elements (Fig. 2). A single layer of infinite elements is added outside the domain in order to reproduce the behaviour in outer space. Spectral and infinite elements share the same interpolation functions, namely Lagrange polynomials, but use different quadratures. Thus, the magnetic potential Φ is discretized in natural coordinates |$\boldsymbol{\xi }$| as
(7)
where Φα denotes the value of magnetic potential at quadrature point |$\boldsymbol{\xi }_\alpha$|⁠, and Nα an interpolation function. The total number of quadrature points in an element is denoted by n and is given by the product of the number of quadrature points in each dimension, nj, j = 1, 2, 3, that is, |$n=\prod _{j=1}^3 n^j$|⁠. The interpolation functions Nα in natural coordinates are determined by the tensor product of 1-D Lagrange polynomials, that is
(8)
such that
(9)
Here α denotes the index of quadrature point |$\boldsymbol{\xi }_\alpha =\lbrace \xi _{\alpha ^1},\xi _{\alpha ^2},\xi _{\alpha ^3}\rbrace$|⁠.
The unmagnetized region (light grey elements) and the magnetized region (interior dark grey elements) are discretized using spectral elements. A single layer of infinite elements is added outside the domain of interest to capture outer space (outer dark grey elements).
Figure 2.

The unmagnetized region (light grey elements) and the magnetized region (interior dark grey elements) are discretized using spectral elements. A single layer of infinite elements is added outside the domain of interest to capture outer space (outer dark grey elements).

Open in new tabDownload slide
The test function w is taken to be an interpolation function Nα, making the approach a Galerkin method. Upon substituting w = Nα and Φ given by eq. (7), in eq. (6), we obtain a set of elemental linear equations that may be written conveniently in matrix–vector form:
(10)
The quantities |$\mathsf {K}_\mathrm{ e}$| and |$\mathsf {F}_\mathrm{ e}$| are known, respectively, as the stiffness matrix and force vector of an element. Similarly, |$\mathsf {\Phi }_\mathrm{ e}$| is the magnetic potential vector. Symbolically, we write
(11)
Here
(12)
(13)
and
(14)
where
(15)
After assembling the elemental matrices and vectors, we obtain a set of global linear equations
(16)
where |$\mathsf {K}$| and |$\mathsf {F}$| are known, respectively, as the global stiffness matrix and global force vector. Similarly, |$\mathsf {\Phi }$| is the global magnetic potential vector.
After determining the magnetic potential Φ, we compute the perturbed magnetostatic field or the magnetic field intensity, |$\mathbf {H}$|⁠, by simply calculating a numerical derivative within a spectral element:
(17)
where |$\mathbf {J}^{-1}_\alpha$| denotes the inverse of the Jacobian of the mapping at node α.

One key advantage of the SEM is that derivatives are directly computed on quadrature points, rendering this operation highly accurate. However, derivatives may be discontinuous across elemental boundaries, and therefore on such boundaries we average over all elements that touch a given node.

Finally, the magnetic induction or the magnetic field, |$\mathbf {B}$| , is given by
(18)
where, μ0 denotes the magnetic permeability of free space and has a value |$4\, \pi \, \times \, 10^{-7}$| H m−1.
For a given regional magnetic field, |$\mathbf {F}$|⁠, the total-field anomaly can be approximated as (Blakely 1995)
(19)
where, |$\hat{\mathbf {F}}$| denotes the direction of the regional field |$\mathbf {F}$|⁠.

As described in detail in Gharti et al. (2018) within spectral elements we use GLL quadrature to construct the stiffness matrix, whereas in infinite elements we use Gauss–Radau quadrature in the infinite direction and GLL quadrature in other directions. Quadrature points on a spectral-infinite element interface coincide, thereby naturally coupling the two domains.

3 NUMERICAL EXAMPLES

We rely on MeshAssist (Gharti et al.2017) and Trelis/CUBIT (CUBIT 2017) for model preparation and meshing for all examples included in this paper. Each spectral element involves three GLL points in each direction, resulting in a total of 27 nodes per element. Our method is fully parallelized based on domain decomposition and the message passing interface. We use the mesh partitioning tool SCOTCH (Pellegrini & Roman 1996) to partition the mesh. We implemented parallel iterative Krylov solvers using PETSc, a portable and extensible toolkit for scientific computation (Balay et al.2015). We use a conjugate gradient solver with a geometric algebraic multigrid preconditioner (e.g. Stüben 2001) with a relative tolerance of 10−7. We have previously performed extensive convergence and strong scaling performance tests of our solver (Gharti et al.2018, 2019a).

3.1 Spherical magnetization

In this section we consider a spherical magnetization of radius 2 m buried at a depth of z = 5 m, as shown in Fig. 3. The sphere has a magnetization strength, |$M = \Vert \mathbf {M}\Vert = 10$| A m−1. We assume that the sphere is located at four different latitudes, namely, λ = 0°, 15°, 30° and 90° , such that the inclination of the magnetization, I  varies according to the geocentric dipole approximation
(20)
and we impose zero magnetic declination for all locations.
Model geometry for a spherical magnetization.
Figure 3.

Model geometry for a spherical magnetization.

Open in new tabDownload slide

We created a model of size 40 m × 20 m × 12 m that covers the sphere and the observation points. For reasons of accuracy, we slightly extended the model above the ground surface to avoid placing the infinite element layer adjacent to observation points. We meshed the model using hexahedral elements. Initially, we selected an element size of 0.5 m inside the sphere and 1 m outside the sphere, as shown in Fig. 4 (left-hand panel). The mesh consists of 8576 elements and honours the sphere. We added an infinite-element layer outside the model surface based on a pole positioned at the centre of the sphere, as shown in Fig. 4 (right-hand panel). Opposite faces of infinite elements must be parallel or diverging with respect to the pole so that they do not converge at infinity. We performed four different simulations for the magnetized sphere located at four different latitudes. Simulations were run on a single processor and took about 40 s for each simulation. The analytical expression for the magnetic induction due to a sphere is given by eq. (A1).

Left-hand panel: spectral-element mesh for a spherical magnetization. Right-hand panel: infinite elements radiating from the outer surface of the model (outer grey). The model is sectioned to visualize one quadrant of the magnetized sphere (dark grey). Only infinite elements from the visible outer surfaces are shown for clarity.
Figure 4.

Left-hand panel: spectral-element mesh for a spherical magnetization. Right-hand panel: infinite elements radiating from the outer surface of the model (outer grey). The model is sectioned to visualize one quadrant of the magnetized sphere (dark grey). Only infinite elements from the visible outer surfaces are shown for clarity.

Open in new tabDownload slide

Fig. 5 shows total-field anomaly profiles calculated at the ground surface, z = 0, along the x-axis from −20 to 20 m. The computed result is in good agreement with the analytical solution, although we observe some discrepancies near the largest magnitudes. Both magnitude and direction of the total-field anomaly vary with the location of the magnetization. For the sphere located at latitudes, 0° and 90°, the total-field anomaly is symmetric about the centre of the sphere. The total-field anomaly is asymmetric about the centre of the sphere for other locations.

Profiles for the total field anomaly at z = 0 along the x-axis (North) for a spherical magnetization.
Figure 5.

Profiles for the total field anomaly at z = 0 along the x-axis (North) for a spherical magnetization.

Open in new tabDownload slide

We refined the mesh by decreasing the element size to 0.25 and 0.5 m, respectively, inside and outside the dome, which increases the number of elements to 53 615. As shown in Fig. 5, the resulting total-field anomaly profile is very close to the analytical result, but we still observe some small discrepancies near the largest magnitudes. We further refined the mesh with element sizes of 0.125 and 0.25 m, respectively, inside and outside the dome, thereby increasing the total number of elements to 156 538. This refinement gives a nearly perfect match with the analytical solution. The spherical geometry poses a challenge for the mesher, because the coarser meshes fail to properly capture the spherical shape of the magnetization.

In Fig. 6, we plot the components of the magnetic induction computed on the ground surface of the model. As expected, we observe different patterns for each component depending on the location of the magnetized spheres. The figure also illustrates how the magnetic field gradually changes according to the angle of magnetic inclination. For magnetic inclinations of 0° and 90°, all three components are either symmetric or antisymmetric about the centre of the sphere. This is not the case for other locations. The magnetic field is stronger near the magnetized sphere and decays away from it.

Magnetic field components (left-hand panel: Bx; middle: By; right-hand panel: Bz) plotted on the ground surface for a magnetized sphere located at latitudes 0°, 15°, 30° and 90° (top to bottom).
Figure 6.

Magnetic field components (left-hand panel: Bx; middle: By; right-hand panel: Bz) plotted on the ground surface for a magnetized sphere located at latitudes 0°, 15°, 30° and 90° (top to bottom).

Open in new tabDownload slide

Finally, we plot the magnetic field lines on the vertical mid-plane, as shown in Fig. 7. It is interesting to see how the magnetic field lines emanate differently from the spheres for each location due to the different magnetic inclinations. These field lines explain the symmetry or asymmetry we observed for the total-field anomaly (Fig. 5).

Magnetic field lines plotted on a mid-vertical plane for a magnetized sphere.
Figure 7.

Magnetic field lines plotted on a mid-vertical plane for a magnetized sphere.

Open in new tabDownload slide

3.2 Rectangular prism

Next, we consider a rectangular prism of size 3 km × 4 km × 1 km with a magnetization strength M = 10 A m−1 located at a depth of 2.5 km from the ground surface, as shown in Fig. 8. The magnetization has an inclination of 30° and a declination of 0°. We created a model of size 28 km × 16 km × 8 km, which sufficiently covers the prism and the observation points. We extended the model slightly above the ground surface to avoid placing the infinite-element layer adjacent to observation points.

Model geometry for a magnetized rectangular prism.
Figure 8.

Model geometry for a magnetized rectangular prism.

Open in new tabDownload slide

Initially, we selected a coarse element size of 1 km and meshed the model using hexahedral elements. The mesh consists of only 4176 elements and honours the prism, as shown in Fig. 9. A single layer of infinite elements was added outside the mesh domain. We positioned the pole for the decay function required by the infinite elements at the centre of the prism.

Left-hand panel: spectral-element mesh of a model with a buried rectangular prism. The model is sectioned to visualize the prism (dark grey).
Figure 9.

Left-hand panel: spectral-element mesh of a model with a buried rectangular prism. The model is sectioned to visualize the prism (dark grey).

Open in new tabDownload slide

An analytical expression for the total-field anomaly of a prism is given by eq. (A5). Fig. 10 shows a total-field anomaly profile computed at the ground surface, z = 0. We observe good agreement with the analytical result. Although we used a very coarse mesh, we obtain acceptable results because the mesh perfectly honours the prism.

Profile for the total field anomaly at z = 0 for a rectangular prism.
Figure 10.

Profile for the total field anomaly at z = 0 for a rectangular prism.

Open in new tabDownload slide

Next, we refined the mesh by decreasing the element size from 1 to 0.5 km. The refinement increases the number of elements to 33 408. The resulting profile for the total-field anomaly is in excellent agreement with the analytical solution, demonstrating rapid convergence in terms of mesh refinement.

In Fig. 11, we plot the components of the magnetic induction computed on the ground surface. Due to the magnetic inclination, the magnetic field is not symmetric. It is stronger near the magnetized prism and decays ways from the prism.

Magnetic field components (Top: Bx; Middle: By; Bottom: Bz) plotted on the ground surface for a rectangular prism.
Figure 11.

Magnetic field components (Top: Bx; Middle: By; Bottom: Bz) plotted on the ground surface for a rectangular prism.

Open in new tabDownload slide

Finally, we plot the magnetic field lines on the mid-vertical plane as shown in Fig. 12. The magnetic field lines emanate from and converge at the magnetic source at 30°, that is, the angle of inclination.

Magnetic field lines plotted on the mid-vertical plane for a rectangular prism.
Figure 12.

Magnetic field lines plotted on the mid-vertical plane for a rectangular prism.

Open in new tabDownload slide

3.3 Multiple spheres with varying size and magnetization

In this section, we design a problem similar to the example discussed in Gharti et al. (2018) and Martin et al. (2017) for gravity anomalies. We consider 100 spheres with varying radii and magnetization strengths organized in a vertical plane, as shown in Fig. 13. The model consists of 10 horizontally equidistant lines in the vertical plane with a vertical spacing of 10 km. Each horizontal line consists of 10 equidistant spheres with a horizontal spacing of 14 km. The lth horizontal line has a constant sphere radius, rl, and magnetization strength, M, determined by the relations, |$r_l=2+0.2\, (l-1)$| and |$M=1.5\, l$| A m−1, respectively, such that lines |$l=1\, \mathrm{to}\, 10$| are numbered top to bottom. The first sphere in the first horizontal line is located at x = 40 km, y = 0 km and z = − 5 km from the origin. We set the magnetic inclination and declination to 35° and 10°, respectively, for all spheres.

Layout of an array of 100 spheres of varying size.
Figure 13.

Layout of an array of 100 spheres of varying size.

Open in new tabDownload slide

We created a model of size 206 km × 80 km × 170 km covering all spheres and observation points. The outer surfaces of the model lie 40 km in each direction from the nearest sphere centre. We selected an element size of |$0.2\, r_\mathrm{ s}$| inside a sphere and |$0.4\, r_\mathrm{ s}$| outside a sphere, where rs is the sphere’s radius. This approach gradually coarsens the mesh from top to bottom, varying the element size inside the sphere from 0.4 to 0.76 km and outside the sphere from 0.8 to 1.52 km. This strategy ensures the same number of elements in all spheres (Fig. 14). Although not strictly necessary, this maintains a relatively uniform resolution for all spheres. The mesh consists of 2 886 000 spectral elements for a total of 22 808 007 degrees of freedom. We used 160 processors, on which it took about 7 min to compute the magnetic field. We set the magnetic inclination and declination of the regional field to 35° and 10°, respectively.

Left-hand panel: spectral-element mesh for an ensemble of 100 spheres. Mesh elements on the spheres are not shown for clarity. Right-hand panel: zoom in on a frontal view of the mesh. Spheres with the same colour have the same magnetic properties.
Figure 14.

Left-hand panel: spectral-element mesh for an ensemble of 100 spheres. Mesh elements on the spheres are not shown for clarity. Right-hand panel: zoom in on a frontal view of the mesh. Spheres with the same colour have the same magnetic properties.

Open in new tabDownload slide

As a benchmark, we computed the total-field anomaly using MBOX (Blakely 1995). MBOX divides the model into a number of rectangular prisms. Each subprism uses the analytical solution given in eq. (A5) and the total-field anomaly is given by the superposition of the contributions of all the subprisms. We discretized the model using uniform cubic cells. Since the smallest spectral element has a size of 0.4 km and each spectral element has three GLL points in each direction, we chose a cubic cell of size 0.2 km. Consequently, the MBOX grid resolution is similar to the finest region of the spectral-element mesh, and hence the MBOX model is sampled more finely than the spectral-element model. There are a total of 1 611 488 cells with non-zero magnetization. We set 293 observation points each at two elevations, namely, 0 and 25 km, along the x-axis, resulting in a total of 586 observation points. The calculations took about 45 min on a single processor. Since the computational cost for the MBOX solution is proportional to the number of observation points for a given discretization, the total computational cost for the 22 808 007 points, that is, the total number of points in the spectral-element mesh, would be huge. On the other hand, the total computational cost for the SIEM is independent of the number of observation points for a given discretization.

We plot total-field anomaly profiles at two elevations along the x-axis, specifically at 0 and 25 km, as shown in Fig. 15. For both profiles, we observe excellent agreement of the SIEM and the MBOX solutions, with the analytical solution given by eq. (A3). At zero elevation, we clearly observe the influence of individual spheres as indicated by peaks and troughs. Magnitude of the peaks decreases and that of the troughs increases along x, that is, North direction due to the magnetic inclination. At far away from the sphere ensemble, that is, elevation 25 km, individual influences diminish and all spheres behave as a single entity as indicated by a smooth curve.

Total-field anomaly profiles at elevations of z = 0 km and z = −25 km for an ensemble of spheres.
Figure 15.

Total-field anomaly profiles at elevations of z = 0 km and z = −25 km for an ensemble of spheres.

Open in new tabDownload slide

In Fig. 16, we plot the components of the magnetic field on the ground surface. We clearly observe the influence of each column of magnetized spheres. Although the spheres are symmetrically placed about the x-axis, the magnetic field is not symmetric due to the inclination and declination.

Magnetic field components (Top: Bx; Middle: By; Bottom: Bz) plotted on the ground surface for an array of magnetic spheres.
Figure 16.

Magnetic field components (Top: Bx; Middle: By; Bottom: Bz) plotted on the ground surface for an array of magnetic spheres.

Open in new tabDownload slide

Finally, we plot magnetic filed lines on the mid-vertical planes through the sphere centres, as shown in Fig. 17. We observe the complicated interactions of the magnetic fields contributed by the individual spheres. Outside the spheres ensemble, magnetic field lines are regular as the individual influence of the spheres diminishes.

Magnetic field lines plotted on the mid-vertical plane for an array of magnetic spheres.
Figure 17.

Magnetic field lines plotted on the mid-vertical plane for an array of magnetic spheres.

Open in new tabDownload slide

3.4 Application to an underground ore mine

To demonstrate our method for realistic 3-D problems, we compute the magnetic field for an existing underground ore mine, namely, the Pyhäsalmi mine in Finland. This active, deep mine contains massive sulphide deposits with copper and zinc ore bodies (Puustjärvi 1999). The complete mine layout and surrounding infrastructure are shown in Fig. 18 (left-hand panel). For more details on the mine and its microseismic event characteristics, see Oye et al. (2005). Fig. 18 (right-hand panel) shows the 622 m × 622 m × 622 m 3-D mine model, which consists of host rock, ore body and stopes, that is, mined-out cavities. Due to the mine’s complex geometry, it is a very challenging problem to calculate the magnetic field. We assume uniform magnetization within the ore body with a magnetization strength of 20 A m−1, a magnetic inclination of 75° and an declination of 10°.

Left-hand panel: Pyhäsalmi mine with surrounding infrastructure: copper/zinc ore body (brown/pink), access tunnels (yellow), elevator shaft (dark blue) and seismic stations (numbered). The passage for quarried ore is labelled KN1. Right-hand panel: 3-D model of the Pyhäsalmi mine: stopes (blue) and ore body (brown). The remainder is host rock. Geophones are indicated by black dots and numbered.
Figure 18.

Left-hand panel: Pyhäsalmi mine with surrounding infrastructure: copper/zinc ore body (brown/pink), access tunnels (yellow), elevator shaft (dark blue) and seismic stations (numbered). The passage for quarried ore is labelled KN1. Right-hand panel: 3-D model of the Pyhäsalmi mine: stopes (blue) and ore body (brown). The remainder is host rock. Geophones are indicated by black dots and numbered.

Open in new tabDownload slide

We meshed the model using an average element size of 10 m. Due to the complex geometry, actual element sizes may differ in some areas, as shown in Fig. 19. The mesh involves 294 500 spectral elements and 2 517 769 degrees of freedom. We used 80 processors for the parallel simulation which took about 1 min.

Left-hand panel: spectral-element mesh for the Pyhäsalmi mine model. Right-hand panel: interior section of the mesh visualizing the ore body (brown) and stopes (cyan).
Figure 19.

Left-hand panel: spectral-element mesh for the Pyhäsalmi mine model. Right-hand panel: interior section of the mesh visualizing the ore body (brown) and stopes (cyan).

Open in new tabDownload slide

Fig. 20 shows the components of the magnetic field computed on the top surface, that is, at z = 0 m.

Magnetic field components (Top: Bx; Middle: By; Bottom: Bz) plotted on the ground surface for an ore mine.
Figure 20.

Magnetic field components (Top: Bx; Middle: By; Bottom: Bz) plotted on the ground surface for an ore mine.

Open in new tabDownload slide

Additionally, we plot magnetic field lines on the mid-vertical plane, as shown in Fig. 21. Imprints of cavities are observed in the field lines.

Magnetic field lines plotted on the mid-vertical plane for an ore mine.
Figure 21.

Magnetic field lines plotted on the mid-vertical plane for an ore mine.

Open in new tabDownload slide

Finally, we also compute the total-field anomaly on the top surface and compare with the MBOX solution, as shown in Fig. 22. For the MBOX simulation we use a cell size of 3.8875 m. There are a total of 25 921 observation points on the top surface. Since only the ore body is magnetized, there are only 149 739 non-zero cells. The MBOX calculation takes about 2 hr 42 min to run on a single processor to compute the total-field anomaly at 25 921 observation points. This is significantly larger than the SIEM, which takes only 1 min to run on 80 processors for 2 517 769 observation points. Since MBOX’s computational cost is proportional to the number of observation points, the total computational time would be huge if we had to compute the solution for all 2 517 769 observation points

Total-field anomaly plotted on the ground surface for an ore mine. (Left-hand panel) MBOX. (Right-hand panel) SIEM.
Figure 22.

Total-field anomaly plotted on the ground surface for an ore mine. (Left-hand panel) MBOX. (Right-hand panel) SIEM.

Open in new tabDownload slide

The total-field anomaly computed by the SIEM and MBOX are very similar in pattern, but the magnitude of the MBOX total-field anomaly is less than the SIEM. The cubic cells used for the MBOX solution are unable to accurately capture the complex shape of the ore body and generally underestimate the surface area or the volume of complicated 3-D objects.

4 CONCLUSIONS

We have successfully implemented an SIEM to solve the unbounded Poisson/Laplace equation for magnetic anomalies induced by magnetized objects. We benchmarked our results with analytical solutions and the MBOX package (Blakely 1995) for a range of problems and used the Pyhäsalmi mine as a practical 3-D example.

Our tests show that the method is highly accurate and efficient. Since there is only a single layer of infinite elements, the additional computational cost associated with ‘outer space’ is insignificant. The computational cost is independent of the number of observation points. Consequently, the method is very efficient for large-scale problems. For additional performance improvement, GPU acceleration will be of future interest.

Inversion of magnetic data is also a future interest, in particular, the planetary magnetic fields (Holme & Bloxham 1996; Lewis & Simons 2012; Plattner & Simons 2017). The main advantage of our numerical discretization of the unbounded Poisson equation is that it can be naturally coupled with the laws of continuum mechanics that govern geostatic and geodynamic deformation.

ACKNOWLEDGEMENTS

We thank Volker Oye, Katja Sahala and ISS for access to the mine model. We thank Frederik J. Simons for helpful discussions which inspired this work and important feedback on this manuscript. Parallel programs were run on computers provided by the Princeton Institute for Computational Science and Engineering (PICSciE). 3-D data were visualized using the open-source parallel visualization software ParaView/VTK (www.paraview.org). 2-D snapshots were created using the open-source Python plotting library Matplotlib (www.matplotlib.org). Most of the model sketches were created using the open-source vector graphics editor Inkscape (inkscape.org). This research was partially supported by NSF grants 1644826 and 1550901. Our software is open source and freely available via the Computational Infrastructure for Geodynamics (CIG; geodynamics.org). We thank Editor Prof. Richard Holme, Alain Plattner and an anonymous reviewer for their insightful comments which helped to improve the manuscript.

REFERENCES

Angelov
T.A.
,
1991
.
Infinite elements–theory and applications
,
Comput. Struct.
,
41
(
5
),
959
–
962
..
10.1016/0045-7949(91)90288-W

Google Scholar

OpenURL Placeholder Text

Crossref

 

Balay
S.
et al. .,
2015
.
PETSc Users Manual, Tech. Rep. ANL-95/11 - Revision 3.6
,
Argonne National Laboratory
.

Google Scholar

OpenURL Placeholder Text

 

Barnett
C.T.
,
1976
.
Theoretical modeling of the magnetic and gravitational fields of an arbitrarily shaped threedimensional body
,
Geophysics
,
41
(
6
),
1353
–
1364
..
10.1190/1.1440685

Google Scholar

OpenURL Placeholder Text

Crossref

 

Beer
G.
,
Meek
J.L.
,
1981
.
‘Infinite domain’ elements
,
Int. J. Numer. Methods Eng.
,
17
(
1
),
43
–
52
..
10.1002/nme.1620170104

Google Scholar

OpenURL Placeholder Text

Crossref

 

Bettess
P.
,
1977
.
Infinite elements
,
Int. J. Numer. Methods Eng.
,
11
(
1
),
53
–
64
..
10.1002/nme.1620110107

Google Scholar

OpenURL Placeholder Text

Crossref

 

Bhattacharyya
B.K.
,
1964
.
Magnetic anomalies due to prismshaped bodies with arbitrary polarization
,
Geophysics
,
29
(
4
),
517
–
531
..
10.1190/1.1439386

Google Scholar

OpenURL Placeholder Text

Crossref

 

Blakely
R.J.
,
1995
.
Potential Theory in Gravity and Magnetic Applications
,
Cambridge Univ. Press
.

Google Scholar

OpenURL Placeholder Text

 

Bott
M.H.P.
,
1963
.
Two methods applicable to computers for evaluating magnetic anomalies due to finite three dimensional bodies
,
Geophys. Prospect.
,
11
(
3
),
292
–
299
..
10.1111/j.1365-2478.1963.tb02039.x

Google Scholar

OpenURL Placeholder Text

Crossref

 

Briais
A.
,
Patriat
P.
,
Tapponnier
P.
,
1993
.
Updated interpretation of magnetic anomalies and seafloor spreading stages in the south China Sea: implications for the tertiary tectonics of Southeast Asia
,
J. geophys. Res.
,
98
(
B4
),
6299
–
6328
..
10.1029/92JB02280

Google Scholar

OpenURL Placeholder Text

Crossref

 

Cai
Y.
,
Wang
C.Y.
,
2005
.
Fast finite-element calculation of gravity anomaly in complex geological regions
,
J. geophys. Int.
,
162
(
3
),
696
–
708
..
10.1111/j.1365-246X.2005.02711.x

Google Scholar

OpenURL Placeholder Text

Crossref

 

Canuto
C.
,
Hussaini
M.Y.
,
Quarteroni
A.
,
Zang
T.A.
,
1988
.
Spectral Methods in Fluid Dynamics
,
Springer
.

Google Scholar

OpenURL Placeholder Text

 

CUBIT
,
2017
.
CUBIT 15.3 User Documentation
,
Sandia National Laboratories
,
Online accessed 2017 January 27
.

Google Scholar

OpenURL Placeholder Text

Dahlen
F.A.
,
1974
.
On the static deformation of an earth model with a fluid core
,
J. geophys. Int.
,
36
(
2
),
461
–
485
..
10.1111/j.1365-246X.1974.tb03649.x

Google Scholar

OpenURL Placeholder Text

Crossref

 

Deville
M.O.
,
Fischer
P.F.
,
Mund
E.H.
,
2002
.
High-Order Methods for Incompressible Fluid Flow
,
Cambridge Univ. Press
.

Google Scholar

OpenURL Placeholder Text

 

El-Esnawy
N.A.
,
Akl
A.Y.
,
Bazaraa
A.S.
,
1995
.
A new parametric infinite domain element
,
Finite Elem. Anal. Des.
,
19
(
1-2
),
103
–
114
..
10.1016/0168-874X(94)00060-S

Google Scholar

OpenURL Placeholder Text

Crossref

 

Faccioli
E.
,
Maggio
F.
,
Paolucci
R.
,
Quarteroni
A.
,
1997
.
2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method
,
J. Seismol.
,
1
,
237
–
251
..
10.1023/A:1009758820546

Google Scholar

OpenURL Placeholder Text

Crossref

 

Gharti
H.N.
,
Tromp
J.
,
2017
.
A spectral-infinite-element solution of Poisson’s equation: an application to self gravity
, preprint (arXiv:1706.00855).

 

Gharti
H.N.
,
Komatitsch
D.
,
Oye
V.
,
Martin
R.
,
Tromp
J.
,
2012a
.
Application of an elastoplastic spectral-element method to 3D slope stability analysis
,
Int. J. Numer. Methods Eng.
,
91
,
1
–
26
.

Google Scholar

OpenURL Placeholder Text

 

Gharti
H.N.
,
Oye
V.
,
Komatitsch
D.
,
Tromp
J.
,
2012b
.
Simulation of multistage excavation based on a 3D spectral-element method
,
Comput. Struct.
,
100–101
,
54
–
69
..
10.1016/j.compstruc.2012.03.005

Google Scholar

OpenURL Placeholder Text

Crossref

 

Gharti
H.N.
,
Langer
L.
,
Roth
M.
,
Tromp
J.
,
Vaaland
U.
,
Yan
Z.
,
2017
.
MeshAssist: An Open-Source and Cross-Platform Meshing Assistant Tool.
,
Zenodo
.

Google Scholar

OpenURL Placeholder Text

 

Gharti
H.N.
,
Tromp
J.
,
Zampini
S.
,
2018
.
Spectral-infinite-element simulations of gravity anomalies
,
J. geophys. Int.
,
215
(
2
),
1098
–
1117
..
10.1093/gji/ggy324

Google Scholar

OpenURL Placeholder Text

Crossref

 

Gharti
H.N.
,
Langer
L.
,
Tromp
J.
,
2019a
.
Spectral-infinite-element simulations of coseismic and post-earthquake deformation
,
J. geophys. Int.
,
216
(
2
),
1364
–
1393
..
10.1093/gji/ggy495

Google Scholar

OpenURL Placeholder Text

Crossref

 

Gharti
H.N.
,
Langer
L.
,
Tromp
J.
,
2019b
.
Spectral-infinite-element simulations of earthquake-induced gravity perturbations
,
J. geophys. Int.
,
217
(
1
),
451
–
468
..
10.1093/gji/ggz028

Google Scholar

OpenURL Placeholder Text

Crossref

 

Hansen
R.O.
,
Wang
X.
,
1988
.
Simplified frequencydomain expressions for potential fields of arbitrary threedimensional bodies
,
Geophysics
,
53
(
3
),
365
–
374
..
10.1190/1.1442470

Google Scholar

OpenURL Placeholder Text

Crossref

 

Heirtzler
J.R.
,
Dickson
G.O.
,
Herron
E.M.
,
Pitman
W.C.
,
Le Pichon
X.
,
1968
.
Marine magnetic anomalies, geomagnetic field reversals, and motions of the ocean floor and continents
,
J. geophys. Res.
,
73
(
6
),
2119
–
2136
..
10.1029/JB073i006p02119

Google Scholar

OpenURL Placeholder Text

Crossref

 

Hejlesen
M.M.
,
Rasmussen
J.T.
,
Chatelain
P.
,
Walther
J.H.
,
2013
.
A high order solver for the unbounded poisson equation
,
J. Comput. Phys.
,
252
(
0
),
458
–
467
..
10.1016/j.jcp.2013.05.050

Google Scholar

OpenURL Placeholder Text

Crossref

 

Holme
R.
,
Bloxham
J.
,
1996
.
The magnetic fields of Uranus and Neptune: methods and models
,
J. geophys. Res.
,
101
(
E1
),
2177
–
2200
..
10.1029/95JE03437

Google Scholar

OpenURL Placeholder Text

Crossref

 

Kumar
P.
,
1985
.
Static infinite element formulation
,
J. Struct. Eng.
,
111
(
11
),
doi:10.1061/(ASCE)0733-9445(1985)111:11(2355))
.

Google Scholar

OpenURL Placeholder Text

 

Lewis
K.W.
,
Simons
F.J.
,
2012
.
Local spectral variability and the origin of the martian crustal magnetic field
,
Geophys. Res. Lett.
,
39
(
18
),
L18201
,
doi:10.1029/2012GL052708
.

Google Scholar

OpenURL Placeholder Text

 

Martin
R.
et al. .,
2017
.
A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data-application to the western Pyrenees
,
J. geophys. Int.
,
209
(
1
),
406
–
424
.

Google Scholar

OpenURL Placeholder Text

 

Mavriplis
C.
,
1989
.
Laguerre polynomials for infinite-domain spectral elements
,
J. Comput. Phys.
,
80
(
2
),
480
–
488
..
10.1016/0021-9991(89)90112-5

Google Scholar

OpenURL Placeholder Text

Crossref

 

Medina
F.
,
Taylor
R.L.
,
1983
.
Finite element techniques for problems of unbounded domains
,
Int. J. Numer. Methods Eng.
,
19
(
8
),
1209
–
1226
..
10.1002/nme.1620190808

Google Scholar

OpenURL Placeholder Text

Crossref

 

Okabe
M.
,
1979
.
Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies
,
Geophysics
,
44
(
4
),
730
–
741
..
10.1190/1.1440973

Google Scholar

OpenURL Placeholder Text

Crossref

 

Oye
V.
,
Bungum
H.
,
Roth
M.
,
2005
.
Source parameters and scaling relations for mining related seismicity within the Pyhaesalmi ore mine, Finland
,
Bull. seism. Soc. Am.
,
95
,
1011
–
1026
..
10.1785/0120040170

Google Scholar

OpenURL Placeholder Text

Crossref

 

Patera
A.T.
,
1984
.
A spectral element method for fluid dynamics: laminar flow in a channel expansion
,
J. Comput. Phys.
,
54
,
468
–
488
..
10.1016/0021-9991(84)90128-1

Google Scholar

OpenURL Placeholder Text

Crossref

 

Pellegrini
F.
,
Roman
J.
,
1996
.
SCOTCH: a software package for static mapping by dual recursive bipartitioning of process and architecture graphs
,
Lecture Notes Comput. Sci.
,
1067
,
493
–
498
.

Google Scholar

OpenURL Placeholder Text

 

Peter
D.
et al. .,
2011
.
Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes
,
J. geophys. Int.
,
186
(
2
),
721
–
739
..
10.1111/j.1365-246X.2011.05044.x

Google Scholar

OpenURL Placeholder Text

Crossref

 

Plattner
A.
,
Simons
F.J.
,
2017
.
Internal and external potential-field estimation from regional vector data at varying satellite altitude
,
J. geophys. Int.
,
211
(
1
),
207
–
238
.

Google Scholar

OpenURL Placeholder Text

 

Plouff
D.
,
1976
.
Gravity and magnetic fields of polygonal prisms and application to magnetic terrain corrections
,
Geophysics
,
41
(
4
),
727
–
741
..
10.1190/1.1440645

Google Scholar

OpenURL Placeholder Text

Crossref

 

Puustjärvi
H.
,
1999
.
Pyhaesalmi modeling project, section B. Geology
,
Tech. rep., Geological Survey of Finland, and Outokumpu Mining Oy
.

Google Scholar

OpenURL Placeholder Text

 

Seriani
G.
,
Oliveira
S.P.
,
2008
.
Dispersion analysis of spectral-element methods for elastic wave propagation
,
Wave Motion
,
45
,
729
–
744
..
10.1016/j.wavemoti.2007.11.007

Google Scholar

OpenURL Placeholder Text

Crossref

 

Stüben
K.
,
2001
.
A review of algebraic multigrid
,
J. Comput. Appl. Math.
,
128
(
1
),
281
–
309
..
10.1016/S0377-0427(00)00516-1

Google Scholar

OpenURL Placeholder Text

Crossref

 

Talwani
M.
,
1965
.
Computation with the help of a digital computer of magnetic anomalies caused by bodies of arbitrary shape
,
Geophysics
,
30
(
5
),
797
–
817
..
10.1190/1.1439654

Google Scholar

OpenURL Placeholder Text

Crossref

 

Tromp
J.
,
Mitrovica
J.X.
,
1999
.
Surface loading of a viscoelastic Earth–I. General theory
,
J. geophys. Int.
,
137
(
3
),
847
–
855
..
10.1046/j.1365-246x.1999.00838.x

Google Scholar

OpenURL Placeholder Text

Crossref

 

Tromp
J.
,
Komatitsch
D.
,
Liu
Q.
,
2008
.
Spectral-element and adjoint methods in seismology
,
Commun. Comput. Phys.
,
3
(
1
),
1
–
32
.

Google Scholar

OpenURL Placeholder Text

 

Tsynkov
S.V.
,
1998
.
Numerical solution of problems on unbounded domains: a review
,
Appl. Numer. Math.
,
27
(
4
),
465
–
532
..
10.1016/S0168-9274(98)00025-7

Google Scholar

OpenURL Placeholder Text

Crossref

 

Vine
F.J.
,
Matthews
D.H.
,
1963
.
Magnetic anomalies over oceanic ridges
,
Nature
,
199
,
947
–
949
.

Google Scholar

OpenURL Placeholder Text

 

Wang
X.
,
Hansen
R.O.
,
1990
.
Inversion for magnetic anomalies of arbitrary three-dimensional bodies
,
Geophysics
,
55
(
10
),
1321
–
1326
..
10.1190/1.1442779

Google Scholar

OpenURL Placeholder Text

Crossref

 

Won
I.J.
,
Bevis
M.
,
1987
.
Computing the gravitational and magnetic anomalies due to a polygon: algorithms and fortran subroutines
,
Geophysics
,
52
(
2
),
232
–
238
..
10.1190/1.1442298

Google Scholar

OpenURL Placeholder Text

Crossref

 

Zienkiewicz
O.C.
,
Emson
C.
,
Bettess
P.
,
1983
.
A novel boundary infinite element
,
Int. J. Numer. Methods Eng.
,
19
(
3
),
393
–
404
..
10.1002/nme.1620190307

Google Scholar

OpenURL Placeholder Text

Crossref

 

APPENDIX A: ANALYTICAL SOLUTIONS

In this section, we list all analytical solutions used in our benchmarks. In what follows, depth z is positive downwards.

A1 Dipole or sphere

The magnetic field at a point |$\mathbf {x}_p$| due to a sphere of radius a with a magnetization |$\mathbf {M}$| and its centre located at |$\mathbf {x}_s$| is given by (Blakely 1995)
(A1)
where
(A2)
and the position vector |$\mathbf {r}$| of the observation point is measured from the centre of the sphere, that is, |$\mathbf {r}=\mathbf {x}_p-\mathbf {x}_s$|⁠. The point |$\mathbf {x}_p$| lies outside the sphere, and |$r=\Vert \mathbf {r}\Vert$|⁠.

A2 Multiple spheres

The magnetic field at a point |$\mathbf {x}_p$| due to ns spheres of radii as and magnetizations |$\mathbf {M}_s$| located at points |$\mathbf {x}_s$| is obtained by superposing the magnetic fields of individual spheres:
(A3)
where
(A4)
and the position vector |$\mathbf {r}_s$| of the observation point is measured from the centre of the individual spheres, that is, |$\mathbf {r}_s=\mathbf {x}_p-\mathbf {x}_s$|⁠. The point |$\mathbf {x}_p$| lies outside the spheres and |$r_s=\Vert \mathbf {r}_s\Vert$|

A3 Semi-infinite rectangular prism

The total-field anomaly observed at the origin due to a buried semi-infinite rectangular prism with a magnetization |$\mathbf {M}=M_x\, \mathbf {i}+M_y\, \mathbf {j}+M_z\, \mathbf {k}$| and dimensions x1 ≤ x ≤ x2, y1 ≤ y ≤ y2 and z1 ≤ z ≤ ∞ under a regional field directed parallel to |$\hat{\mathbf {F}}=\hat{F}_x\, \mathbf {i}+\hat{F}_y\, \mathbf {j}+\hat{F}_z\, \mathbf {k}$| is given by (Bhattacharyya 1964; Blakely 1995):
(A5)
where,
(A6)
and
(A7)

We may compute the total-field anomaly for a rectangular prism of thickness h by evaluating eq. (A5) twice. The total-field anomaly for this case is equal to the total-field anomaly evaluated for z1 minus the total-field anomaly evaluated for z1 = z1 + h.

© The Author(s) 2019. Published by Oxford University Press on behalf of The Royal Astronomical Society.
This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
Subject
Geomagnetism and Electromagnetism
Issue Section:
Geomagnetism, rock magnetism and paleomagnetism
Download all slides