Abstract
Geophysical electromagnetic survey methods are particularly effective in locating conductive mineral deposits or mineralization zones in a mineral resource exploration. The forward modelling of the electromagnetic responses over such targets is a fundamental task in quantitatively interpreting the geophysical data into a geological model. Due to the ubiquitous irregular and complex geometries associated with the mineral rock units, it is critical that the numerical modelling approach being used is able to adequately and efficiently incorporate any necessary geometries of the Earth model. To circumvent the difficulties in representing complex but necessary geometry features in an Earth model for the existing mesh-based numerical modelling approaches (e.g., finite element and finite difference methods), I present a meshfree modelling approach that does not require a mesh to solve the Maxwell’s equations. The meshfree approach utilizes a set of unconnected points to represent any geometries in the Earth model, allowing for the maximal flexibility to account for irregular surface geometries and topography. In each meshfree subdomain, radial basis functions are used to construct meshfree function approximation in transforming the differential equations in the modelling problem into linear systems of equations. The method solves the potential function equations of the Maxwell’s equations in the modelling. The modelling accuracy using the meshfree method is examined and verified using one magnetotelluric model and two frequency-domain controlled-source models. The magnetotelluric model is the well-known Dublin Test Model 2 in which the spherical geometry of the conductor in the shallow subsurface may pose as a challenge for many numerical modelling methods. The first controlled-source model is a simple half-space model with the electric dipole source for which analytical solutions exist for the modelling responses. The second controlled-source model is the volcanic massive sulphide mineral deposit from Voisey’s Bay, Labrador, Canada in which the deposit’s surface is highly irregular. For all modellings, the calculated electromagnetic responses are found to agree with other independent numerical solutions and the analytical solutions. The advantages of the meshfree method in discretizing the Earth models with complex geometries in the forward modelling of geophysical electromagnetic data is clearly demonstrated.
1 Introduction
Geophysical electromagnetic (EM) survey methods continue to be important in the exploration of mineral resources, particularly those with a high conductivity contrast from their host rocks (e.g., copper, zinc, iron, nickle) (; ; ; ; ). In recent years, due to the increasing acknowledgement of the important role of mineral resources in energy transition, various “critical mineral resource” initiatives have been proposed (e.g., ) and how we as a society can meet the demands has sparkled much discussion ().
A geophysical EM survey directly produces a map of the distribution of electrical conductivity of the subsurface. Naturally, the interpretation of any EM data collected over a region of interest becomes vital in determining the parameters of potential mineral deposits that may host economic resources. At the centre of a quantitative interpretation of EM data is the numerical modelling of EM responses (including controlled-source EM, magnetotelluric, transient EM data) which plays a critical role in the development of theories and methods of various EM survey techniques. Of the various advancements made over the last few decades (e.g., ; ; ), the numerical modelling of EM data has steadily evolved from closed-form, analytical computations of EM responses over relatively simple conductivity models that started around 1960s (; ) to fully numerical simulations of Maxwell’s equations over Earth models with complex geometries and nonlinear, anisotropic conductivity distributions nowadays (e.g., ; ).
The importance of representing realistically complex geometries of conductive mineral deposits or mineralization zones in the EM data modelling becomes obvious since mineral deposits or mineralization zones are naturally of irregular shapes of geometry, with some presenting quite extreme geometries (e.g., uranium deposits associated with thin graphite ; ). Despite the ubiquitous existence and importance of such realistic geometries, there are still significant challenges faced by numerical modelling techniques in terms of efficiently incorporating the geometries. These challenges are precisely what this study is trying to solve and in order to do so, a new type of modelling techniques called meshfree methods is used which I will present in detail in the following sections.
Numerical methods of forward modelling EM responses over a general three-dimensional (3-D) conductivity Earth model are often focused on mesh-based modelling methods in the applied geophysics which include finite difference (e.g., ; ; ; ; ; ), finite volume (e.g., ; ), integral equation (e.g., ; ; ; ; ) and finite element methods (e.g., ; ; ; ; ; ; ). They are termed mesh-based modelling methods in this study because they have the common feature of relying on a mesh-based discretization (e.g., rectilinear, triangular and tetrahedral meshes; see Figure 1A,B) of the conductivity model. Among these mesh-based methods, finite difference approaches may face more challenges than others in accurately representing complex topography surfaces and irregular surface geometries of a conductor since they require tensor-grid function approximation of differential equations. In contrast, finite volume, integral equation and finite element methods do not face such limitation. It may be argued that finite element modelling techniques, if combined with unstructured meshes whose automatic generations are facilitated by modern mesh generation software (; ), are the most flexible mesh-based approaches in the modelling of EM data over complex Earth models (; ; ; ).
FIGURE 1
For real-life geometries of exploration targets, unstructured meshes (e.g., triangular and tetrahedral meshes) possess a unique advantage in efficiently and accurately representing complex geometries that are important characteristics of potential exploration targets (; ). However, accurate numerical solution of EM responses using mesh-based modelling techniques, including finite element and finite volume methods, also require a certain degree of regularity of the mesh cells. In the finite element case, for example, the effect of the mesh quality (e.g., the ratio of the largest to smallest cell sizes, elongation, dihedral angles and radius-edge ratio of cells for tetrahedral meshes) on the computational accuracy is demonstrated to be significant (). Poor mesh quality can lead to very slow convergence or divergence in iteratively solving the resulting linear system of equations in modelling controlled-source EM data using a vector finite element implementation (). On the other hand, ensuring the quality of the mesh can lead to overwhelmingly excessive number of elements in the generated mesh, therefore intractable computational resources, in order to sufficiently conform to the real geometries in the model (; ).
The dilemma in balancing the quality of unstructured meshes and the number of mesh cells is often addressed using adaptive mesh refinement techniques (; ; ; ; ; ). In an adaptive mesh refinement approach, the current mesh used for the modelling of EM responses can be further refined or coarsened based on an estimate of the current numerical modelling error. The ideal result is that only the part of the mesh with the largest numerical errors is refined. In the unstructured mesh scenario, the adaptive refining process is often carried out by locally modifying the mesh, rather than re-meshing the whole Earth model, due to the concern of efficiency and robustness of the process (; ). However, complications may arise during the adaptive mesh refining. First, as the topology of the mesh changes at each refining step the mapping between the old mesh and the new mesh needs to be calculated in order to update the degrees of freedom, which can be an expensive and rather complicated process. Second, further dividing the cells in the current mesh may produce “hanging nodes” due to the non-conforming new cells within the parent cells (). The nonconformity of the new mesh may be eliminated at the cost of further refining neighbouring cells, often with a lower quality of the generated new cells.
Alternatively, the Earth model can be discretized using
meshfree points(see
Figure 1C) and the corresponding numerical modelling techniques are called
meshfree methods(
;
). A set of unconnected points, or meshfree points, serve for the same purpose as that of a quality mesh in obtaining an accurate numerical solution in forward modelling the EM data. Because of the lack of connectivity among the points, the physical property distribution (i.e., the conductivity distribution for EM data modelling) will be
sampledon the points when discretizing the partial differential equations. With a meshfree point discretization, the density and regularity of the point distribution are still important for accurate numerical modellings; however, the advantages of manipulating points over mesh generations are:
• With comparable regularity of a quality mesh, the generation of points requires much less effort in computer programming and is more straightforward. Also, the development of dedicated point generation software is also significantly easier ().
• Since there is no topology requirement, the generation of quality, unstructured meshfree points can be more robust than generating a quality unstructured mesh (; ).
• Adaptive point refining and/or coarsening is more computationally efficient than the same process when using meshes, since any addition or deletion of local points does not need to affect the rest of the points. The nonconformity issue and its complications in a mesh-based adaptive refining are completely removed ().
Based on a distribution of points, many meshfree methods for solving partial differential equations have been proposed (). For geophysical data modelling, however, only a few different meshfree methods have been proposed for seismic wave field modelling (e.g., ; ; ), gravity data modelling (e.g., ) and EM data modelling (e.g., ; ; ; ). In general, different meshfree methods differ in the choice of basis functions, the types of meshfree points (i.e., uniform or unstructured) and in that whether numerical integration is required in transforming the partial differential equations into the linear system of equations. The meshfree method demonstrated here, mostly known as radial-basis-function based finite difference (RBF-FD), does not need the potentially expensive step of numerical integration. It also naturally supports unstructured point distributions allowing for an efficient discretization of complex-geometry conductivity models. Meshfree modelling of EM data is considered to be more challenging than those of seismic and gravity data since the EM fields are discontinuous across conductivity discontinuities ().
The rest of the study is organized as follows. The details of the meshfree modelling method in the context of numerically solving Maxwell’s equations will be first presented, which is followed by the demonstration of the numerical accuracy of the method using a magnetotelluric example and two controlled-source EM examples. Further discussions for the applicability for other types of geophysical data of the modelling method are also presented before I conclude the study.
2 Methods
2.1 Maxwell’s equations for meshfree modelling
The frequency-domain Maxwell’s equations for the electromagnetic field in the quasi-static limit are expressed as ()for Faraday’s law and Ampère’s law, respectively. Here, and are the electric field and magnetic induction vector, respectively. with as the magnetic field and the magnetic permeability. is the conductivity distribution of the Earth model. with as the ordinary frequency in Hz, is the imaginary unit, and the convention of the time dependence is used here. represents any external EM sources as a current density vector; for example, the current density of an induction loop or of an electric dipole grounded into the Earth.
Eliminating in Eqs 1, 2 through simple substitutions leads to the second-order Helmholtz equation for :Here, the Earth materials are assumed to be non-ferromagnetic so that the magnetic permeability is just that of free space . As a result, Eq. 3 is further simplified as:A naive solving of Eq. 4 using numerical methods may lead to spurious or incorrect numerical solutions of EM responses if the discontinuous nature of at conductivity jumps is not considered. The ability of handling such discontinuities is one of the reasons behind the popularity of Yee-scheme finite difference methods () and vector finite element methods () when numerically solving Eq. 4.
In the meshfree modelling of EM responses, the degrees of freedom of the unknown function (e.g., in Eq. 4) are coincident with the point locations in a point discretization of the Earth model, a scenario similar to scalar finite element methods (). As demonstrated in detail by , the RBF-FD meshfree method using scalar meshfree basis functions, like scalar finite element methods, will force the numerical solution of the unknown function to be continuous everywhere. In this scenario, EM potential function equations can be used instead of the Helmholtz equation for the electric field. Using the vector magnetic potential and electric scalar potential defined via the relations ():we have the Helmholtz equation for the potential functions:which is obtained by substituting Eq. 5 into Eq. 4. Also, taking the divergence of Eq. 7 gives us the conservation of charge equation for the potential functions:In Eq. 8, the distribution of electric charges resulting from EM sources such as grounded electric dipoles is represented by the term . It is well known that the ungauged potential equations, Eqs 7, 8, does not provide a unique solution of the pair , despite that the electric and magnetic fields, which are calculated using Eqs 5 and 6, from solving the potential equations may still be unique (). Here, the Coulomb gauge condition is applied to Eq. 7 to stabilize the numerical solution (; ). Taking advantage of the vector identity , the Coulomb-gauged Helmholtz equation for the potential functions becomesBoth potential functions, and , are continuous across any conductivity jumps. In fact, the vector potential is also smooth at conductivity jumps (). The component-wise form of the pair of Eqs 9, 8 which are discretized here using the RBF-FD meshfree method isfor EM data modelling with a general source where represents in previous equations.
2.2 RBF-FD
The description of the RBF-FD meshfree numerical method here follows , ,, , and references therein. Like mesh-based numerical methods, the first step of the RBF-FD is to locally approximate an unknown function, , as some simple, rationale functions. In mesh-based finite difference methods, this is often done using Taylor expansions at a point using linear or quadratic functions depending on the finite difference scheme and order (e.g., first-order backward). In mesh-based finite element methods, this is typically done by using low-order polynomial basis functions within an element (i.e., a cell in the mesh). In the meshfree RBF-FD, the unknown function at any point (called support node, see Figure 2) is locally approximated as a linear combination of translations of a single radial basis function (RBF) using the neighbouring points in the subdomain of that point (Figure 2). Such interpolant can be written aswhere , is the norm, are the interpolation coefficients, and is the position of the th point which is also the center of the corresponding RBF . Note that unlike polynomial functions, a RBF is always radially symmetric around its center (). To determine the interpolation coefficients in Eq. 14, a local linear system of equations resulting from the Lagrange interpolation conditions (, , with as the sampled function values at the points), which can be written as or in a compact matrix formneeds to be solved. The 3-D RBFs , are used in the RBF-FD method for its computational efficiency and robustness in solving the local linear system in Eq. 16 (see detailed discussions in ). It can be proved that using the RBF, the symmetric matrix is always invertible as long as the local points are distinct (). This flexibility of point locations allows for arbitrary point distributions to be used which will be important in representing complex geometries in an Earth model. In practice, Eq. 15 is enriched with additional low-order polynomials to avoid numerical singularity in case the positions of meshfree points in a subdomain become too extreme (e.g., colinear, see ).
FIGURE 2
Using the above meshfree interpolant, any differential operator (e.g., in Eq. 10) can be discretized over the meshfree subdomains in the form of a linear combination of local function values, a process that is similar to the traditional mesh-based finite difference approximation but in a more general treatment: . In RBF-FDs, the discretization of the operator is multi-dimensional, while in the classical mesh-based finite differences the discretization of is restricted to be directional approximation (i.e., only 1-D). The weights, , are then obtained by solving the following local linear system
and form as the nonzeros in the corresponding rows of the coefficient matrix in the resultant global linear system ( is the total number of points). Here, denotes the value of at the location and can be readily calculated using the chain rules of the derivatives for the chosen RBF. is the vector of the unknown function values at the degrees of freedom (i.e., meshfree points in the RBF-FD). The proof of Eq. 17 is thoroughly presented in . The right-hand-side vector is formed from proper boundary conditions and the discretizations of EM source terms. Note that only Eq. 17 needs to be solved in transforming the differential equations (Eq. 10 to Eq. 13) into linear systems of equations. The well-known numerical analysis package LAPACK subroutines were used to numerically solve Eq. 17. In this study, all meshfree points are unstrutured and the size of meshfree subdomains is fixed as (the number of points in a meshfree subdomain) and the selection of the closest points for each subdomain is carried out using a kd-tree point selection method which is the same as in . The fixed number of points in subdomains means that the relative distances among the points in a meshfree subdomain can be smaller (e.g., near CSEM sources) when high numerical accuracies are needed.
3 Numerical results
In this section, the modellings of different EM data using the meshfree method are demonstrated. The global linear system from discretizing the - potential function equations is asymmetric, complex-valued and can be solved using either iterative solvers or direct solvers. In this study, all global linear systems are solved using the MUMPS direct solver (, version 5.3.3).
3.1 Magnetotelluric data
The magnetotelluric (MT) conductivity model for the demonstration of the meshfree modelling is the Dublin Test Model 2 (DTM2, ) in which a hemispherical conductor is buried at the top of the subsurface (Figure 3). The hemispherical conductor has the resistivity of 10 m and the background earth’s resistivity is 300 m. The radius of the hemisphere is km. Because of the perfect symmetry of the conductivity structure, there exists analytical solutions of MT responses at the galvanic limit (i.e., zero-frequency limit). For the same reason, the observed MT responses will be symmetric. These features of this model make it a good example for the comparison of different numerical modelling algorithms (). However, also due to the spherical surface geometry of the conductor, mesh-based modelling techniques, particularly those relying on tensor-grid or rectilinear meshes, will face challenges in accurately representing the geometry, and therefore in studying the effects of shallow inhomogeneities of the conductivity distribution on the MT data at sites near the edge of the conductor.
FIGURE 3
In MT data modellings, the EM source terms in Eqs 10–13, i.e., and its divergence, vanish as the actual EM sources are far away from the surface of the Earth (). In this model, the boundary conductivity distribution is that of the uniform subsurface and 1-D boundary conditions were used to compute boundary values (see details in ). The EM responses in the MT scenario are typically represented using apparent resistivity and phase data which are calculated from the electric and magnetic fields at the measurement locations. Although the previous study using the RBF-FD method () has demonstrated the effectiveness of the modelling capability, particularly how the discontinuous electric field can be correctly modelled using the developed RBF-FD method here, it does not demonstrate the flexible meshfree discretization of highly irregular surface geometries as we see in the DTM2 here. The unstructured point discretization with local refinements for this model is shown in Figure 4. The total number of points in the discretization is for a computational domain of .
FIGURE 4
Three MT sites (Figure 3) are designed here to examine the modelling accuracy of the meshfree method. Site 1 is at the origin of the coordinate system and at the center of the hemisphere conductor. Site 2 ( m; m) is 500 m away from the edge of the hemisphere and is inside the hemisphere (same as Site 10 in ). Site 3 ( m; m) is 100 m away from the edge of and outside the hemisphere (same as Site 18 in ). MT responses at Site 2 and 3 are expected to be significantly affected by the irregular shape of the hemisphere for long periods. Same as , the frequency range of Hz to 100 Hz (periods from 0.01 s to 10,000 s) were used for the examination. The calculated MT responses at the three sites using the meshfree RBF-FD method were compared with other independent solutions (all using mesh-based modelling methods) that are documented from and are shown in Figures 5–7 for the three sites. At each site, the apparent resistivity and phase data for the four components of the impedance tensor (i.e., , , and ) are plotted from top to bottom.
FIGURE 5
FIGURE 6
FIGURE 7
At Site 1 (Figure 5), which is at the origin of the model, the theoretical apparent resistivity of the MT responses for and are zero, which explains the extremely small and random numerical values of the apparent resistivity and phase data observed for all numerical solutions. For the off-diagonal components and , almost all numerical solutions agree with each other. At Site 2 (Figure 6) and Site 3 (Figure 7), all four components of the impedance tensor will be non-zero due to the edge effect of the conductor. For the phase data in the and components at these two sites, the meshfree solution appears to deviate from other solutions; this is because of the difference in the Coordinate systems being used and the phase angle calculation methods.1 The meshfree numerical solutions are validated by the following two observations: the symmetry in the solution (MT responses for and are the same, so are the and ) and the good agreement with the majority of other independent solutions. Note that among those independent solutions, a few solutions (e.g., Kiyan and Khoza) have a clear deviation from the main MT response curves in long periods (after 10 s for Site 2) due to insufficient mesh discretizations around the edge of the hemisphere conductor.
3.2 Frequency-domain controlled-source EM data
To compute the controlled-source EM (CSEM) responses, the external source terms in Eqs 10–13 ( and its divergence) will be non-zero at the locations of the source. Here, in the context of the meshfree RBF-FD method, the source handling method of is used. As shown in Figure 8, any CSEM source wire is initially represented by meshfree points in space with proper distances and regularity of distribution among them. For each point representing the source wire, a local unstructured mesh is contructed by connecting the points found in the subdomain of the point. Then a finite element-like weak-form treatment (e.g., ; ) is used to discretize the equations at the source point in which the electrical current of the wire will be coincident with the edges of the local mesh which is of tetrahedral type in this study. The number of nodes in this local mesh is very small and the mesh connectivity is generated automatically using common mesh generation software (e.g., Tetgen and Gmsh). For grounded wires, the divergence of the current density is only non-zero at the beginning and the ending source points.
FIGURE 8
The above method is capable of treating an arbitrarily shaped controlled source (grounded wires or current loops) in the point discretization of an Earth model. Because of this capability, the total-field approach of modelling the EM data, as described in Eqs 10–13 in the case of potential functions, is being used here and will provide more flexibility in handling complex topography and surface geometries. Under the total-field approach, the boundary values of the EM field on the computational domain is zero.
3.2.1 Half-space model
The first CSEM test model is that of a uniform subsurface with an electric dipole source at the Earth’s surface (Figure 9). Although the model is relatively ideal, analytical solutions exist for the EM fields at the surface which allows for a first-step examination of the accuracy of the developed RBF-FD meshfree method for CSEM data modellings. An -directed electric dipole, for example, will have the current density as ()where is the current intensity, is Heaviside function, and are the two ends of the grounded wire, and is the Dirac delta function. The closed-form, analytical formula (eq. 4.159 in ) for computing the inline electric field due to the above CSEM transmitter (Eq. 18) for any measurement locations at the surface of the subsurface (i.e., at ) is:where , is the wavenumber with as the conductivity of the subsurface. In the case of dipole sources where the length of the grounded wire approaches infinitesimally small in relative to the distance from the dipole to measurement locations, in Eq. 19.
FIGURE 9
The current density of the electric dipole source is set to be 1 A. For the meshfree solution, a set of unstructured points with local refinements around the dipole source was used. The actual CSEM source in the meshfree modelling is 1 m in length along the direction and is represented by six points with the equal spacing of 0.2 m from m to m. Four frequencies at and Hz were used for the accuracy examination. The uniform Earth’s subsurface has the conductivity of 0.02 S/m. With this conductivity value, the very low frequency of Hz will approach the direct current limit and the electric field responses will approach those of a direct current resistivity survey. The total number of the points for the model discretization is 81,951 which are distributed within the computational domain where the dipole is located at the center. The computed responses at the surface are shown in Figure 10. As evident from the comparison, the two solutions have an excellent agreement with each other, demonstrating the computational accuracy of the meshfree method. At the highest frequency Hz (Figure 10A), the imaginary from the meshfree numerical solution starts to deviate from the theoretical solution when km. This deviation is due to the decreased point density after the distance. At lower frequencies, such deviation vanishes due to less rapid changes of the EM field at the same locations. In addition, when the frequency approaches the direct current limit (Figure 10C,D), the frequency-domain response has the real part almost unchanged but the imaginary part continuously decreased until zero; that is, response is approaching the direct current resistivity response, which is another evidence of the correct modelling of the CSEM responses.
FIGURE 10
3.2.2 Ovoid mineral deposit model
The second CSEM test example is from the volcanic massive sulphide mineral deposit (termed as Ovoid deposit here) from Voisey’s Bay, Labrador, Canada, which has been under extensive studies for both geology and geophysical data modelling studies (e.g., ; ; ; ). The Ovoid deposit is a highly conductive iron-dominant deposit with complex surface geometry which serves as a perfect testing example for geophysical data modelling software. To test the meshfree modelling method developed here, the real topography of the Earth’s surface (also see ; ) has been used here. The nearest point of the deposit to the surface is about 70 m below the uneven surface (see Figure 11). Following , a grounded wire of 400 m long with a current density of 1 A along the easting direction (Figure 11A,B) is used as the source for a ground EM survey. The offset of the wire source is roughly 600 m away from the central part of the deposit. There are 143 measurement sites at the surface designed as receiver locations which are distributed evenly along 11 South-to-North profiles (Figure 11C). The profile spacing is 100 m. The site spacing along each profile is approximately 50 m. The highly conductive deposit is assumed to have a uniform conductivity of 1 S/m here and the relatively resistive background earth is assigned the conductivity of 0.001 S/m. Both conductivity values are taken from the previous test model of for facilitating a direct comparison of the meshfree modelling results with theirs.
FIGURE 11
For the model discretization, the meshfree points used here are also directly taken from the point distribution from the tetrahedral mesh used by which are visualized in Figure 11D. In their modelling, they have used a vector finite element method for computing the CSEM responses. For both finite element methods and meshfree methods, the local refinements around the transmitter and the receiver locations are necessary to improve the numerical modelling accuracy. The 400-m long grounded wire is represented in the meshfree method using 80 points located at the topographic surface with in average 5 m of distance apart from each other (see Figure 8). The total number of points for this discretization is 44,230 within the computational domain of . For long grounded wire sources, the inline electric field and vertical magnetic field components are often the main measurements. The computed and responses using the meshfree method are compared with the finite element results of which are shown in Figures 12–14 at the measurement site profiles m (also see Figure 11C), respectively, for the frequency of 500 Hz. It is seen that the two independent solutions, calculated using the same model discretization (i.e., meshfree points and tetrahedral mesh), have a very good agreement with each other for all sites (see supplemental materials for the comparison for the remaining profiles). The strong EM induction caused by the deposit is well reflected at the sites that are more closely above the deposit (see the profiles m and m). The higher frequency (i.e., Hz) responses for the profile (Profile 5) is also shown in Figure 15. Comparing Figure 15 with Figure 14, it is seen that at the higher frequency, the EM responses, particularly responses, attenuate faster in space.
FIGURE 12
FIGURE 13
FIGURE 14
FIGURE 15
4 Discussions
Direct current resistivity (DCR) survey methods are also frequently used for mineral resource exploration. Although I have not directly shown the modelling capability of the meshfree method for DCR data modelling, the first CSEM modelling example (in the case of ) is essentially a simple demonstration of how the developed meshfree modelling approach can be directly applied to compute DCR data for a general conductivity model. When the frequency is zero, the potential function equations described in Eqs 10–13 will be reduced to the exact potential function equation (i.e., Poisson’s equation for ) used for DCR data computation.
For model discretization, the unstructured meshfree points were generated using a combination of existing open-source software tools (including Tetgen and Paraview). The geophysical community is likely to continue to benefit from research in other fields. In the situation of dedicated software development for meshfree point generation, there have been significant research and development in the past decade (). It is anticipated that, like the history of the classical finite element methods, more open-source point generation tools will be available once the meshfree methods along with their capability of incorporating complex geometries become more widely known.
5 Conclusion
Earth models in the context of mineral resource exploration using geophysical survey methods often have rather complex surface geometries. It is important that the numerical forward modelling of geophysical EM data for such models is capable of efficiently handling these geometries. A meshfree modelling method that uses only unconnected points, instead of the traditional pixel cell-based meshes, to represent geometries has been developed and presented here. The - potential equations instead of the Helmholtz equation for the electric field are used for the continuity property of the potential functions. The meshfree method supports both uniform and non-uniform, unstructured point distributions with the latter being of particularly advantageous in discretizing Earth models with complex geometries with a minimal amount of points.
The modelling accuracy and the capability of handling highly irregular surface geometries of the method are demonstrated using three EM modelling examples. The first example is a magnetotelluric model in which the magnetotelluric impedance responses of a hemisphere-shaped near-surface conductor were modelled. The second example is an idealized half-space conductivity model excited with a grounded electric dipole source for which closed-form analytical solutions exist. The third example, which is also a controlled-source example, is the real-life highly conductive Ovoid mineral deposit model in which the realistic surface geometry of the deposit and the topography were used. The EM transmitter for this example is a 400 m long grounded wire. Through these examples, the feasibility of easily representing irregular surface geometries of the Earth models is clearly demonstrated. The point discretizations are considered to be more advantageous over traditional mesh discretizations for complex Earth models as they are easier to generate and manipulate. For all examples, the modelling accuracies of the meshfree method are verified using other independent numerical solutions or analytical solutions.
The demonstrated meshfree modelling method is also applicable to other geophysical data modellings in which numerical solutions and complex geometries of the model are important. The developed meshfree method for geophysical EM data modellings is shown to be effective for both natural-source magnetotelluric surveys and controlled-source EM surveys. For the magnetotelluric example, the meshing of the spherical surface geometry is a non-trivial task, especially for numerical methods that are restricted to the use of rectilinear meshes, as evident from the large differences of modelled magnetotelluric responses among some independent solutions (). The representation of such geometry is however quite straightforward and easy in the meshfree point discretization.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
JL: Conceptualization, Data curation, Investigation, Methodology, Software, Writing–original draft, Writing–review and editing.
Funding
The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.
Acknowledgments
Kara and Farquharson are thanked for graciously providing the mesh data to build the Ovoid deposit model discretization and for sharing their finite element solutions to compare with in this study.
Conflict of interest
The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Footnotes
1.^The solutions provided from only have the computed phase angles, rather than the impedance values themselves
References
1
AmestoyP.DuffI. S.KosterJ.L’ExcellentJ.-Y. (2001). A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM J. Matrix Analysis Appl.23, 15–41. 10.1137/s0895479899358194
2
AnsariS.FarquharsonC. G. (2014). 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids. Geophysics79, E149–E165. 10.1190/geo2013-0172.1
3
BadeaE. a.EverettM. E.NewmanG. a.BiroO. (2001). Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials. Geophysics66, 786–799. 10.1190/1.1444968
4
BuhmannM. D. (2003). Radial basis functions. Cambridge: Cambridge University Press. 10.1017/CBO9780511543241
5
ChenC.KruglyakovM.KuvshinovA. (2021). Advanced three-dimensional electromagnetic modelling using a nested integral equation approach. Geophys. J. Int.226, 114–130. 10.1093/gji/ggab072
6
ChenJ.-S.HillmanM.ChiS.-W. (2017). Meshfree methods: progress made after 20 years. J. Eng. Mech.143, 04017001. 10.1061/(asce)em.1943-7889.0001176
7
CoggonJ. H. (1971). Electromagnetic and electrical modeling by the finite element method. Geophysics36, 132–155. 10.1190/1.1440151
8
DuQ.GunzburgerM.JuL. (2002). Meshfree, probabilistic determination of point sets and support regions for meshless computing. Comput. methods Appl. Mech. Eng.191, 1349–1366. 10.1016/s0045-7825(01)00327-9
9
DuQ.WangD.ZhuL. (2009). On mesh geometry and stiffness matrix conditioning for general finite element spaces. SIAM J. Numer. Analysis47, 1421–1444. 10.1137/080718486
10
DyckA. V.WestG. F. (1984). The role of simple computer models in interpretations of wide-band, drill-hole electromagnetic surveys in mineral exploration. Geophysics49, 957–980. 10.1190/1.1441741
11
FabriA.GiezemanG.-J.KettnerL.SchirraS.SchönherrS. (2000). On the design of CGAL a computational geometry algorithms library. Softw. Pract. Exp.30, 1167–1202. 10.1002/1097-024x(200009)30:11<1167::aid-spe337>3.0.co;2-b
12
FarquharsonC. G.CravenJ. A. (2009). Three-dimensional inversion of magnetotelluric data for mineral exploration: an example from the McArthur River uranium deposit, Saskatchewan, Canada. J. Appl. Geophys.68, 450–458. 10.1016/j.jappgeo.2008.02.002
13
FarquharsonC. G.OldenburgD. W. (2002). An integral equation solution to the geophysical electromagnetic forward-modelling problem. Methods Geochem. Geophys. (Elsevier)35, 3–19. 10.1016/S0076-6895(02)80083-X
14
FasshauerG. E. (2007). Meshfree approximation methods with matlab, vol. 6 of Interdiscip. Math. Sci. World Sci. 10.1142/6437
15
FornbergB.FlyerN. (2015). Fast generation of 2-D node distributions for mesh-free PDE discretizations. Comput. Math. Appl.69, 531–544. 10.1016/j.camwa.2015.01.009
16
GehrmannR.NorthL. J.GraberS.SzitkarF.PetersenS.MinshullT.et al (2019). Marine mineral exploration with controlled source electromagnetics at the TAG hydrothermal field, 26°N mid‐atlantic ridge. Geophys. Res. Lett.46, 5808–5816. 10.1029/2019gl082928
17
GüntherT.RückerC.SpitzerK. (2006). Three-dimensional modelling and inversion of DC resistivity data incorporating topography-II. Inversion. Geophys. J. Int.166, 506–517. 10.1111/j.1365-246x.2006.03011.x
18
HanB.LiY.LiG. (2018). 3D forward modeling of magnetotelluric fields in general anisotropic media and its numerical implementation in Julia. Geophysics83, F29–F40. 10.1190/geo2017-0515.1
19
HohmannG. W. (1975). Three-dimensional induced polarization and electromagnetic modeling. Geophysics40, 309–324. 10.1190/1.1440527
20
JahandariH.AnsariS.FarquharsonC. G. (2017). Comparison between staggered grid finite–volume and edge–based finite–element modelling of geophysical electromagnetic data on unstructured grids. J. Appl. Geophys.138, 185–197. 10.1016/j.jappgeo.2017.01.016
21
JahandariH.BihloA.DonzelliF. (2021). Forward modelling of gravity data on unstructured grids using an adaptive mimetic finite-difference method. J. Appl. Geophys.190, 104340. 10.1016/j.jappgeo.2021.104340
22
JahandariH.FarquharsonC. G. (2014). A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids. Geophysics79, E287–E302. 10.1190/geo2013-0312.1
23
JiaX.HuT. (2006). Element-free precise integration method and its applications in seismic modelling and imaging. Geophys. J. Int.166, 349–372. 10.1111/j.1365-246X.2006.03024.x
24
JinJ.-M. (2014). The finite element method in electromagnetics. 3 edn. USA: John Wiley & Sons.
25
JonesA. G. (2023). Mining for net zero: the impossible task. Lead. Edge42, 266–276. 10.1190/tle42040266.1
26
JonesF. W.PascoeL. J. (1972). The perturbation of alternating geomagnetic fields by three-dimensional conductivity inhomogeneities. Geophys. J. Int.27, 479–485. 10.1111/j.1365-246X.1972.tb06103.x
27
KaraK. B.FarquharsonC. G. (2023). 3D minimum-structure inversion of controlled-source EM data using unstructured grids. J. Appl. Geophys.209, 104897. 10.1016/j.jappgeo.2022.104897
28
KeyK.OvallJ. (2011). A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling. Geophys. J. Int.186, 137–154. 10.1111/j.1365-246x.2011.05025.x
29
LelièvreP.Carter-McAuslanA.FarquharsonC.HurichC. (2012). Unified geophysical and geological 3D Earth models. Lead. Edge31, 322–328. 10.1190/1.3694900
30
LiB.LiuY.SenM. K.RenZ. (2017a). Time-space-domain mesh-free finite difference based on least squares for 2D acoustic-wave modeling. Geophysics82, T143–T157. 10.1190/geo2016-0464.1
31
LiJ.FarquharsonC. G.HuX. (2017b). 3D vector finite-element electromagnetic forward modeling for large loop sources using a total-field algorithm and unstructured tetrahedral grids. Geophysics82, E1–E16. 10.1190/geo2016-0004.1
32
LiuZ.RenZ.YaoH.TangJ.LuX.FarquharsonC. (2023). A parallel adaptive finite-element approach for 3-D realistic controlled-source electromagnetic problems using hierarchical tetrahedral grids. Geophys. J. Int.232, 1866–1885. 10.1093/gji/ggac419
33
LongJ.FarquharsonC. G. (2017). “Three-dimensional controlled-source EM modeling with radial basis function-generated finite differences: a meshless approach,” in SEG technical program expanded abstracts 2017 (China: Society of Exploration Geophysicists), 1209–1213.
34
LongJ.FarquharsonC. G. (2019a). “Meshfree modelling of 3-D controlled-source EM data: a new method to treat the singular source terms,” in SEG technical program expanded abstracts 2019 (China: Society of Exploration Geophysicists), 1050–1054.
35
LongJ.FarquharsonC. G. (2019b). On the forward modelling of three-dimensional magnetotelluric data using a radial-basis-function-based mesh-free method. Geophys. J. Int.219, 394–416. 10.1093/gji/ggz306
36
LongJ.FarquharsonC. G. (2019c). Three-dimensional forward modelling of gravity data using mesh-free methods with radial basis functions and unstructured nodes. Geophys. J. Int.217, 1577–1601. 10.1093/gji/ggz115
37
LongJ.FarquharsonC. G. (2020). “Meshfree modelling of 2D MT data with RBF-FD and unstructured points,” in SEG international exposition and annual meeting (SEG). China, SEG.
38
LongJ.FarquharsonC. G. (2024). Three-dimensional controlled-source electromagnetic data modelling with a hybrid meshfree-finite element approach. Germany: preparation.
39
LuX.FarquharsonC. G.MiehéJ.-M.HarrisonG. (2021). 3D electromagnetic modeling of graphitic faults in the Athabasca Basin using a finite-volume time-domain approach with unstructured grids. Geophysics86, B349–B367. 10.1190/geo2020-0657.1
40
MackieR. L.MaddenT. R.WannamakerP. E. (1993). Three-dimensional magnetotelluric modeling using difference equations—theory and comparisons to integral equation solutions. Geophysics58, 215–226. 10.1190/1.1443407
41
MiensopustM. P.QueraltP.JonesA. G.modellersD. M. (2013). Magnetotelluric 3-D inversion—a review of two successful workshops on forward and inversion code testing and comparison. Geophys. J. Int.193, 1216–1238. 10.1093/gji/ggt066
42
NabighianM. N. (1988). Electromagnetic methods in applied geophysics: voume 1, theory. Germany: Society of Exploration Geophysicists.
43
NalepaM.AnsariS.FarquharsonC. (2016). “Finite-element simulation of 3D CSEM data on unstructured meshes: an example from the East Coast of Canada,” in SEG technical program expanded abstracts 2016 (Germany: Society of Exploration Geophysicists), 1048–1052.
44
NamM. J.KimH. J.SongY.LeeT. J.SonJ.-S.SuhJ. H. (2007). 3D magnetotelluric modelling including surface topography. Geophys. Prospect.55, 277–287. 10.1111/j.1365-2478.2007.00614.x
45
NewmanG. A. (2014). A review of high-performance computational strategies for modeling and imaging of electromagnetic induction data. Surv. Geophys.35, 85–100. 10.1007/s10712-013-9260-0
46
NewmanG. A.AlumbaughD. L. (1995). Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophys. Prospect.43, 1021–1042. 10.1111/j.1365-2478.1995.tb00294.x
47
NewmanG. A.HohmannG. W.AndersonW. L. (1986). Transient electromagnetic response of a three-dimensional body in a layered earth. Geophysics51, 1608–1627. 10.1190/1.1442212
48
NguyenV. P.RabczukT.BordasS.DuflotM. (2008). Meshless methods: a review and computer implementation aspects. Math. Comput. Simul.79, 763–813. 10.1016/j.matcom.2008.01.003
49
OdenJ. T.PrudhommeS. (2001). Goal-oriented error estimation and adaptivity for the finite element method. Comput. Math. Appl.41, 735–756. 10.1016/s0898-1221(00)00317-5
50
PridmoreD.HohmannG.WardS.SillW. (1981). An investigation of finite-element modeling for electrical and electromagnetic data in three dimensions. Geophysics46, 1009–1024. 10.1190/1.1441239
51
PuzyrevV.KoldanJ.de la PuenteJ.HouzeauxG.VazquezM.CelaJ. M. (2013). A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling. Geophys. J. Int.193, 678–693. 10.1093/gji/ggt027
52
RabczukT.BelytschkoT. (2005). Adaptivity for structured meshfree particle methods in 2D and 3D. Int. J. Numer. Methods Eng.63, 1559–1582. 10.1002/nme.1326
53
RenZ.KalscheuerT.GreenhalghS.MaurerH. (2013). A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling. Geophys. J. Int.194, 700–718. 10.1093/gji/ggt154
54
RochlitzR.SkibbeN.GüntherT. (2019). custEM: customizable finite-element simulation of complex controlled-source electromagnetic data. Geophysics84, F17–F33. 10.1190/geo2018-0208.1
55
SchulzK. J.DeYoung JrJ. H.Seal IIR. R.BradleyD. C. (2017). Critical mineral resources of the United States—an introduction. Tech. Rep. U. S. Geol. Surv. 10.3133/pp1802A
56
SchwarzbachC.BörnerR.-U.SpitzerK. (2011). Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example. Geophys. J. Int.187, 63–74. 10.1111/j.1365-246X.2011.05127.x
57
SiH. (2015). TetGen, a delaunay-based quality tetrahedral mesh generator. ACM Trans. Math. Softw.41, 1–36. 10.1145/2629697
58
SlakJ.KosecG. (2019). On generation of node distributions for meshless PDE discretizations. SIAM J. Sci. Comput.41, A3202–A3229. 10.1137/18m1231456
59
SmithR. (2014). Electromagnetic induction methods in mining geophysics from 2008 to 2012. Surv. Geophys.35, 123–156. 10.1007/s10712-013-9227-1
60
SpitzerK. (2024). Electromagnetic modeling using adaptive grids–Error estimation and geometry representation. Surv. Geophys.45, 277–314. 10.1007/s10712-023-09794-9
61
StrangwayD.Swift JrC.HolmerR. (1973). The application of audio-frequency magnetotellurics (AMT) to mineral exploration. Geophysics38, 1159–1175. 10.1190/1.1440402
62
StrattonJ. A. (2007). Electromagnetic theory. USA: John Wiley & Sons.
63
StreichR. (2009). 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: direct solution and optimization for high accuracy. Geophysics74, F95–F105. 10.1190/1.3196241
64
TafloveA.UmashankarK. R. (1990). The finite-difference time-domain method for numerical modeling of electromagnetic wave interactions. Electromagnetics10, 105–126. 10.1080/02726349008908231
65
TakekawaJ.MikadaH.ImamuraN. (2015). A mesh-free method with arbitrary-order accuracy for acoustic wave propagation. Comput. Geosciences78, 15–25. 10.1016/j.cageo.2015.02.006
66
WaitJ. R. (1960). Propagation of electromagnetic pulses in a homogeneous conducting earth. Appl. Sci. Res. Sect. B8, 213–253. 10.1007/bf02920058
67
WangT.HohmannG. W. (1993). A finite-difference, time-domain solution for three-dimensional electromagnetic modeling. Geophysics58, 797–809. 10.1190/1.1443465
68
WardS. H.HohmannG. W. (1988). 4. Electromagnetic theory for geophysical applications. Electromagn. Methods Appl. Geophys.1 (4), 130–311. 10.1190/1.9781560802631.ch4
69
WittkeJ.TezkanB. (2014). Meshfree magnetotelluric modelling. Geophys. J. Int.198, 1255–1268. 10.1093/gji/ggu207
70
YeeK. S. (1966). Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans. Antennas Propag.14, 302–307. 10.1109/TAP.1966.1138693
71
ZengS.HuX.LiJ.FarquharsonC. G.WoodP. C.LuX.et al (2019). Effects of full transmitting-current waveforms on transient electromagnetics: insights from modeling the Albany graphite deposit. Geophysics84, E255–E268. 10.1190/geo2018-0573.1
72
ZhangB.YinC.RenX.LiuY.QiY. (2018). Adaptive finite element for 3d time-domain airborne electromagnetic modeling based on hybrid posterior error estimation. Geophysics83, WB71–WB79. 10.1190/geo2017-0544.1
73
ZhdanovM. S. (2010). Electromagnetic geophysics: notes from the past and the road ahead. Geophysics75, 75A49–75A66. 10.1190/1.3483901
Summary
Keywords
mineral exploration, electromagnetic, resistivity, magnetotelluric, controlled-source, meshfree, numerical modelling
Citation
Long J (2024) Meshfree modelling of magnetotelluric and controlled-source electromagnetic data for conductive earth models with complex geometries. Front. Earth Sci. 12:1432992. doi: 10.3389/feart.2024.1432992
Received
15 May 2024
Accepted
13 June 2024
Published
11 July 2024
Volume
12 - 2024
Edited by
Nannan Zhou, Chinese Academy of Sciences (CAS), China
Reviewed by
Hualiang Zhao, Shandong University, China
Jianhui Li, China University of Geosciences Wuhan, China
Updates
Copyright
© 2024 Long.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jianbo Long, jl7037@mun.ca
† Present Address: Jianbo Long, Department of Earth Sciences, Memorial University of Newfoundland, St. John's, NL, Canada
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.