A FORTRAN Program to Model Magnetic Gradient Tensor at High Susceptibility Using Contraction Integral Equation Method

The magnetic gradient tensor provides a powerful tool for detecting magnetic bodies because of its ability to emphasize detailed features of the magnetic anomalies. To interpret field measurements obtained by magnetic gradiometry, the forward calculation of magnetic gradient fields is always necessary. In this paper, we present a contraction integral equation method to simulate the gradient fields produced by 3-D magnetic bodies of arbitrary shapes and high susceptibilities. The method employs rectangular prisms to approximate the source region with the assumption that the magnetization in each element is homogeneous. The gradient fields are first solved in the Fourier domain and then transformed into the spatial domain by 2-D Gauss-FFT. This calculation is performed iteratively until the required accuracy is reached. The convergence of the iterative procedure is ensured by a contraction operator. To facilitate application, we introduce a FORTRAN program to implement the algorithm. This program is intended for users who show interests in 3D magnetic modeling at high susceptibility. The performance of the program, including its computational accuracy, efficiency and convergence behavior, is tested by several models. Numerical results show that the code is computationally accurate and efficient, and performs well at a wide range of magnetic susceptibilities from 0 SI to 1000 SI. This work, therefore, provides a significant tool for 3D forward modeling of magnetic gradient fields at high susceptibility.


Introduction
Magnetic surveying is a significant and widely used geophysical exploration technique. Using magnetic measurements, three kinds of data can be obtained: total magnetic intensity (TMI) data, three-components field data and full tensor magnetic gradient data [1]. Because of easy acquisition, the TMI data and the three-components data were most commonly measured in the past decades [2][3][4][5]. However, these measurements are seriously affected by the Earth's background magnetic field and very sensitive to instrument orientation, hence not just functions of the target's magnetic susceptibility [6]. Recently, full tensor magnetic gradient measurements become available with the development of SQUID-based sensors [1,4,7]. In comparison with the traditional magnetic field measurements, gradient measurements offer many potential advantages, which includes less sensitivity to instrument orientation, improved resolution of shallow targets, no necessity for base stations and diurnal corrections, and suppression of regional anomalies [4,[8][9][10]. Mathematically, the gradient tensor is a second rank tensor with 3 × 3 = 9 components, which is defined by derivatives of the magnetic field vector in each of the three directions of three-dimensional space [11]. Each component of the gradient tensor represents a directional filter and hence makes certain structures and characteristics of the magnetic bodies more noticeable [10]. Therefore, the magnetic gradient tensor can work as a useful tool for detecting magnetic bodies by suppressing specific undesirable contributes and emphasizing certain features of the magnetic fields.
Magnetic forward modeling of gradient fields plays a significant role in the interpretation of field measurements obtained from magnetic gradiometry. During the past decades, various analytic and numerical methods have been developed to simulate magnetic fields. Generally, analytic methods are limited to provide closed-form field expressions for magnetic bodies with regular geometric shapes [12][13][14][15]. For magnetic anomalies with arbitrary and complex susceptibility distribution, numerical methods have to be used. The spatialdomain method is one of the most commonly used numerical methods for computing magnetic fields due to bodies of given shapes [16][17][18][19][20][21][22]. These methods often represent the magnetized body using a set of cells (such as cubical cells) and approximate the total field by the sum of the elementary fields [23,24]. One key limitation of the spatial-domain method is that its computational time dramatically increases with the size and complexity of the model. The frequency-domain method is another important approach to model magnetic anomalies [25][26][27]. This approach is performed in the frequency domain and thereby can achieve high computational efficiency by using fast Fourier transform techniques. On the whole, however, most existing numerical methods ignore the effects of self-demagnetization and are only applicable to magnetic bodies with low susceptibilities.
To overcome the problem,Ref. [28] develop an iterative method and use a segmented model consisting of spherical voxels with arbitrary diameter to calculate the magnetic fields at high susceptibilities through an iterative procedure. Ref. [29] conduct a comprehensive study of the self-demagnetization effects on magnetic data, and compare the capability of two existing inversion methods in interpretation of data from highly magnetic areas. More recently,Ref. [30] develop an efficient and accurate frequency-domain iterative method that can be used to simulate magnetic fields from magnetic bodies with arbitrary shapes in a wide range of magnetic susceptibilities (0∼1000 SI). This strategy is based on a contraction integral equation and can achieve fast convergence. Although these methods perform well at high susceptibilities, they mainly focus on the calculation of the magnetic fields, and take no account of the gradient tensor.
In this paper, we present an algorithm to compute the magnetic gradient fields produced by 3-D magnetic bodies of arbitrary shapes and high susceptibilities based on the contraction integral equation method developed by [30]. In order to simulate the magnetic gradient fields for high susceptibility and to facilitate application of the algorithm of [30], we develop a computer program using FORTRAN language. This FORTRAN program generates multi-component fields, which includes six components of the magnetic gradient tensor and three components of the magnetic field vector, and has a good performance at a wide range of magnetic susceptibilities (0 < χ ≤ 1000 SI), particularly applicable for strongly magnetic bodies.
The remainder of the paper is organized as follows. In Section 2, the algorithm for calculating the magnetic gradient fields from strongly magnetic bodies is developed. Subsequently, a detailed introduction is given to the subroutines, inputs and outputs of the FORTRAN program. Then, in Section 4 the performance of the algorithm is tested using several simple models, and the applicability of the code is demonstrated with a synthetic two-dike model. In Section 5, some conclusions are drawn from the work.

Contraction Integral Equation Method
The full magnetic gradient tensor T can be obtained by taking derivatives of the magnetic vector B with respect to the coordinates x, y and z, where ∂ x , ∂ y , and ∂ z denote the partial derivation with respect to x, y and z, respectively. B = µ 0 H and µ 0 is the permeability of free space. Using we can also write T as It is apparent from Equation (3) that the magnetic gradient tensor T is symmetric. This indicates that using six components of T, namely T xx , T yy , T zz , T xy , T xz and T yz , will be sufficient in the forward calculation. Actually, because of ∇ · B = 0, the magnetic gradient tensor has only five independent components.
To obtain the gradient fields at high susceptibility, we will first calculate the anomalous magnetic fields using the iterative contraction integral equation method, and then derive frequency-domain expressions for the anomalous gradient tensor based on the resulting magnetic field vector in the following sections.

The Integral Equation
In a 3D Cartesian coordinate system with downward positive z direction, the integral equation with respect to the total magnetic field H can be written as [30] where r and r are the observation point and the source point in the source region, respectively. χH denotes the magnetization M. G = ∇∇G is the Green's function tensor, ∇ is the gradient operator, and G = 1 (4π|r − r |) is the scalar Green's function, which satisfies the following differential equation, where δ denotes the Dirac delta function [19].
The anomalous fields are defined as Equation (4) indicates that the total magnetic field H results from two parts: the background field H 0 , which is associated with the Earth's magnetic field (EMF), and the anomalous fields H a , which results from magnetic bodies in the source region. In general, the anomalous fields can be neglected at low magnetic susceptibilities ( χ < 0.01 SI ) [31]. Thus, the magnetization reduces to M = χH 0 in this case, and accurate solutions for the integral Equation (4) can be obtained by direct methods [32,33]. For high susceptibility, however, the induced magnetization from neighboring material significantly affects the magnetic field at any point in the medium, which consequently reduces the resultant magnetic field. This phenomenon is the so-called self-demagnetization [34]. In the case of large susceptibility, the demagnetization effect is significant and cannot be neglected. Therefore, the anomalous magnetic fields must be taken into consideration, and the integral Equation (4) has to be solved by an appropriate iterative procedure.

Iterative Scheme
To obtain magnetic fields at high susceptibilities, [30] developed a convergent iterative scheme that efficiently calculates accurate magnetic fields produced by strongly magnetic bodies. Here, we follow this method to determine the Fourier-domain anomalous magnetic fields and then extend the algorithm to compute the magnetic gradient tensor.
In our algorithm, rectangular prisms are used as the building blocks. The whole computational domain is considered as the source region and is evenly divided into N x × N y × N z prisms in the x, y and z directions (see Figure 1). The geometric center of the prismatic element is (x l , y m , z n ), and the length is ∆x, ∆y and ∆z in the three directions, respectively. The magnetization is assumed to be homogeneous in each prism. To simplify mathematical operations, we start with the anomalous magnetic potential U a [30,35], where G U is the negative gradient of the scalar Green's function, i.e., G U = −∇G and M = χH. Using the spatial discretization strategy aforementioned, we can get the discrete We deal with Equation (8) in the frequency domain of the horizontal coordinates so that the triple integral involved in the equation can reduce to a single integral only associated with z.
Taking the 2D Fourier transform of Equation (8) with respect to x and y, we havê with where k x and k y are wavenumbers in the x and y directions.Û a andM are the corresponding spectrum of U a and M, respectively. According to Here, sign(z − z ) is the sign function and k = k 2 x + k 2 y . Substituting Equations (10)-(12) into Equation (9), the anomalous magnetic potential in the Fourier domain is obtained with whereM x ,M y andM z are three components ofM. Note that P, F and H in Equation (14) are functions ofM, whereM = 2DFT[χH] with 2DFT[·] being a 2D Fourier transform operator. This implies that the values of P, F and H will change with the number of iterations, since the magnetic field vector H is renewed iteratively. In contrast, A 1 , A 2 and e ±kz n are all independent ofM and hence can be pre-computed so as to improve the computational efficiency.
According to Equation (2), we can also write the Fourier-domain anomalous magnetic field asĤ As can be seen, the vertical componentĤ a z presents an implicit form due to the spatial derivative ∂ zÛ a . In order to make it explicit, we substitute Equation (13) into Equation (15) and then haveĤ Therefore, onceĤ a is determined, the 2D inverse Fourier transform ofĤ a immediately gives its spatial-domain counterpart H a .
To solve H a , we rewrite the integral Equation (4) as From Equations (14)- (17), we can find that H a is a function of M (or H, since M = χH). Thus, a convergent iterative algorithm is necessary, in order to obtain the magnetic fields caused by strongly magnetic bodies. But the fact is that the iterative calculation of Equation (4) does not always converge, because 2 -norm of the linear integral operator in Equation (4) is bigger than 1 in general cases [30]. It indicates that a contraction operator must be used to ensure convergence stability and convergence rates.
Therefore, instead of Equation (17), we adopt the contraction integral equation developed by [30] By doing this, an iterative format which makes the successive calculation of the total magnetic field always convergent is established using Equations (13)- (16) and (18).
Our final purpose is to calculate the anomalous magnetic gradient tensor produced by magnetic bodies with high susceptibilities. This can be achieved by taking 2D Fourier transform of Equation (1) and making use of Equation (13). It gives the anomalous gradient tensor asT with whereĤ a is given by Equations (15) and (16). It is noteworthy that the calculation of the gradient tensor should be carried out only when the iterative computation of the magnetic field is totally completed or the required convergence accuracy is reached. Certainly, the gradient tensor can also be calculated iteratively, but this only leads to an unnecessary computational cost. Finally, we obtain the spatial-domain gradient fields through 2D inverse Fourier transform.

Workflow for Modeling Gradient Fields
The workflow for modeling magnetic gradient tensor at high susceptibility is summarized in Figure 2. In this study, the Gauss-FFT technique with 4 nodes is applied, and the root mean square (rms) difference is adopted as the convergence criterion, that is, where H (j) denotes the total magnetic field after the j-th iteration.

Figure 2.
Workflow for modeling the magnetic gradient tensor produced by 3D magnetic bodies with high susceptibility (modified after [30]).

Subroutines
The developed FORTRAN program for calculating magnetic gradient fields at high susceptibilities consists of a main program, six major subroutines and four supporting subroutines. The functions of these subroutines are listed in Table 1.
All subroutines are called by the main program (Main). The computer code initially reads the input file para.txt and sets the required parameters in the subroutine ReadIn. Some significant matrices are pre-computed to improve the computational efficiency in the subroutine Storage. Next, the loop over the iteration number starts to operate. In each iteration, the magnetic fields are calculated layer by layer. For low magnetic susceptibility, having the first iterative cycle done is adequate enough to get accurate results. For high magnetic susceptibility, the total magnetic fields are calculated iteratively until the maximum number of iterations is reached or the required convergence accuracy is met. Finally, the magnetic gradient tensor is computed and then recorded in the output file MagneticField_3D.dat. The structure of the main program is shown in Figure 3.

Inputs and Outputs
The FORTRAN program requires a set of input parameters (see Table 2), such as the size of the computation space, the maximum iteration number and the strength and direction of the inducing magnetic field, etc. All these parameters are given in the input file para.txt. The magnetic susceptibility distribution of the studied magnetic bodies can be set freely in the subroutine Model. In the original version of the subroutine Model, we provide four types of magnetic models, including a sphere, a spherical shell, an ellipsoid and a two-dike model. One can modify the code in this subroutine to establish any expected magnetic model with arbitrary susceptibility distributions. It is also noteworthy that the default values of the minimum x and y coordinates of the prism's center are set to be zero in the FORTRAN program, i.e., x0 = y0 = 0. Table 2. The required input parameters in the file para.txt.

Nx, Ny, Nz
Number of prismatic elements in the x, y and z directions. z0 The minimum z coordinate of geometric center of prisms (m).
x1, y1, z1 The maximum x, y and z coordinates of geometric center of prisms (m). The output of the FORTRAN program is saved in a single file MagneticField_3D.dat. This file contains ten columns: x-coordinate, y-coordinate, z-coordinate, three magnetic field components (B a x , B a y , B a z , nT) and six magnetic gradient tensor components (T a xx , T a yy , T a zz , T a xy , T a xz , T a yz , nT/m).

Numerical Tests
In this section, the performance of the algorithm is evaluated and several numerical examples are presented to reveal the validity and applicability of the code.

Performance
To test the performance of the algorithm, we analyze its computational accuracy, convergence behavior and computational efficiency using several models (a spherical shell, a prolate ellipsoid and a sphere). The center of these magnetic bodies is located in the middle of the computational domain which extends from 0 km to 1 km in the x−, y− and z− directions. The strength, inclination and declination of the inducing magnetic field are 50,000 nT, 60 • and 45 • , respectively. Analytical solutions for the spherical shell model are provided in Appendix A. Because no analytical solution is available for the prolate ellipsoid, the reference gradient tensor of this model is computed by numerically taking derivatives of the magnetic field vector obtained from the routines of [36]. All tests are carried out on a personal computer with 2.3 GHz CPU and 12 GB RAM.
In order to verify the accuracy of the algorithm, we compare the six components of the magnetic gradient tensor calculated using the contraction integral method with the theoretical solutions. In the test, the computational domain is divided into 200 × 200 × 200 regular prisms and the magnetic anomalies are a spherical shell with susceptibility of 50 SI and a prolate ellipsoid with susceptibility of 10 SI (see Figure 4).  Figures 5 and 6 show the numerical and analytical solutions of the two models and the differences between them. The statistical properties of the misfits for both models are listed in Table 3. As can be seen from Figure 5, the modeled results for the spherical shell have a perfect agreement with the analytical solutions. The rms error in this case is 0.082 for T a xx and T a yy , 0.094 for T a xz and T a yz , 0.046 for T a xy and 0.142 for T a zz . The relative root mean square (rrms) errors in six magnetic gradient tensor components are all less than 0.5% (Table 3). For the prolate ellipsoid model, there is also a good agreement between the computed and reference solutions, but with a relatively larger rrms error (see Table 3). Part of such error is considered to be caused by the numerical derivation of the magnetic fields in order to obtain the reference magnetic gradient tensor.  . Surface anomalous magnetic gradient fields produced by a magnetic spherical shell with susceptibility of 50 SI. T a * xx , T a * xy , T a * xz , T a * yy , T a * zz , T a * yz denote the analytical solutions, and T a xx , T a xy , T a xz , T a yy , T a zz , T a yz denote the numerical results calculated using the proposed algorithm.  . Surface anomalous magnetic gradient fields produced by a magnetic prolate ellipsoid with susceptibility of 10 SI. T a * xx , T a * xy , T a * xz , T a * yy , T a * zz , T a * yz denote the analytical solutions, and T a xx , T a xy , T a xz , T a yy , T a zz , T a yz denote the numerical results calculated using the proposed algorithm. Next, we investigate the convergence behavior of the algorithm. To this end, we compute the convergence errors in three total magnetic field components based on a magnetic sphere model with different susceptibilities ranging from 1 SI to 1000 SI. Figure 7 shows the trend of the convergence error versus the number of iterations at different magnetic susceptibilities. We can see that the convergence behavior of the algorithm is significantly influenced by the magnitude of the magnetic susceptibility. At lower susceptibilities (χ ≤ 10 SI), the algorithm converges with a high speed. But when the magnetic susceptibility becomes larger (χ > 10 SI), the convergence rate slows down. Although the convergence speed of the algorithm decreases with an increase of magnetic susceptibility, the algorithm exhibits a good and stable convergence behavior at a rather wide range of susceptibilities on the whole (0 < χ ≤ 1000 SI). Now, we test the computation efficiency of the code. According to the theory of the algorithm mentioned in Section 2, the calculation of the magnetic gradient fields mainly consists of two parts. The first part is the computation of the frequency domain magnetic field at each wavenumber point (k x , k y ). The efficiency of this part is mainly affected by the number of prismatic elements (N z ) used in the z direction. The second part is the implementation of the 2D inverse Gauss-FFT applied to obtain the spatial domain wavefields. For this part, the calculation time is determined by the number of prismatic elements (N x × N y ) on the x − y plane and the Gauss nodes (which is set to be 4 in our test) that one uses. Therefore, N z and N x × N y have different influence on the calculation speed.
To assess these factors, two cases are considered in this test. In case 1, we hold N h = 100, where N h = N x = N y , and set N z = 60, 80, 100, 150, 200, 250, 300. In case 2, we hold N z = 100 and set N h = 60, 80, 100, 150, 200, 250, 300. Figure 8 shows the average computer time per iteration for calculating six magnetic gradient components and three magnetic field components in case 1 and case 2. As can be seen, when the number of the prismatic elements increases, the average computation time exhibits an approximately linear growth trend in both cases.

Example
In this subsection, we present a synthetic example with a slightly more complicated geometry to further demonstrate the capability of the algorithm (Figure 9). This synthetic example is a two-dike model with a vertical dike striking east-west, and a stepped dike that strikes north-south and dips to the west at 45 • . Both dikes spans a depth of 50-150 m and are assigned a high magnetic susceptibility value of 3 SI. The inducing field for this study has strength of 50,000 nT, inclination of 45 • and declination of 45 • . The whole computational space, with cells of 5 m cubed, extends from 0 m to 1000 m in the x and y directions and from 0 m to 250 m in the z direction. Figure 10 presents the surface total-field anomaly and the magnetic amplitude for the two-dike model with a magnetic susceptibility of 3 SI subject to a uniform inducing field. To verify the reliability of the computed results, we also display the magnetic data from a similar two-dike model calculated by [29] (see Figure 11). The geometric parameters of the west-dipping dike used in the work of [29] are slightly different from ours, but the values of other parameters are set the same. By comparing Figures 10 and 11, we can see that the magnetic fields obtained from both methods have consistent shape characteristics, which further validates our algorithm.   In addition to the total-field anomaly, we also present the surface anomalous magnetic fields and gradient fields generated by the two-dike model in Figures 12 and 13, respectively. The white lines in these figures outline the projection of the magnetic bodies on the surface. As expected, the gradient tensor is superior to the magnetic field vector in enhancing the directional features of the magnetic anomaly.

Conclusions
We have developed a computer program in FORTRAN programming language and provided a detailed description of the code for modeling magnetic gradient fields at high susceptibilities. The program generates multi-component fields produced by magnetic bodies of arbitrary shapes and high susceptibilities, including six components of the magnetic gradient tensor and three components of the magnetic field vector. A user friendly interface (i.e., input file) is also established to facilitate wide application of the FORTRAN program. Computational performance and capability of the scheme have been illustrated through several numerical examples (a spherical shell, a prolate ellipsoid, a sphere and a synthetic two-dike model). It is shown that the code performs well at a wide range of susceptibilities (0 < χ ≤ 1000 SI) and is particularly applicable for strongly magnetic bodies. This work, therefore, provides a significant open tool for modeling magnetic gradient fields at high susceptibility.