Three-Dimensional Anisotropic Inversions for Time-Domain Airborne Electromagnetic Data

Rocks and ores in nature usually appear macro-anisotropic, especially in sedimentary areas with strong layering. This anisotropy will lead to false interpretation of electromagnetic (EM) data when inverted under the assumption of an isotropic earth. However, the time-domain (TD) airborne EM (AEM) inversion for an anisotropic model has not attracted much attention. To get reasonable inversion results from TD AEM data, we present in this paper the forward modeling and inversion methods based on a triaxial anisotropic model. We apply three-dimensional (3D) finite-difference on the secondary scattered electric field equation to calculate the frequency-domain (FD) EM responses, then we use the inverse Fourier transform and waveform convolution to obtain TD responses. For the regularized inversion, we calculate directly the sensitivities with respect to three diagonal conductivities and then use the Gauss–Newton (GN) optimization scheme to recover model parameters. To speed up the computation and to reduce the memory requirement, we adopt the moving footprint concept and separate the whole model into a series of small sub-models for the inversion. Finally, we compare our anisotropic inversion scheme with the isotropic one using both synthetic and field data. Numerical experiments show that the anisotropic inversion has inherent advantages over the isotropic ones, we can get more reasonable results for the anisotropic


Introduction
Time-domain (TD) airborne electromagnetic (AEM) is an effective and efficient geophysical exploration method that is being widely used in mineral exploration, groundwater, environmental, and engineering investigations [1,2]. Due to its dense spatial sampling and multi-time channels, the large-scale three-dimensional (3D) modeling and inversion for TD AEM is formerly impossible due to slow computations and large memory requirement. In the past two decades, one-dimensional (1D) inversion and imaging methods have mainly been used in the interpretation of TD AEM data [3][4][5][6][7][8][9]. However, for complex 3D earth structures, 1D inversion and imaging methods cannot deliver reasonable results.
To get more accurate inversion results, 3D inversion has become a research focus for AEM data interpretation in recent years. Sasaki [10] achieved the 3D inversion of frequency domain (FD) EM data based on the finite-difference method and the Gauss-Newton (GN) optimization. However, this inversion code is only suitable for simple models with small datasets and will fail when dealing with the large-scale AEM datasets or complex models. The moving footprint method proposed by [11,12] is a breakthrough of 3D AEM inversion. Based on the compactness of AEM system, this method splits a large-scale model into a series of small sub-models, so that the computation can be applied with reasonable computational requirements. When the AEM response and sensitivity matrix of submodels are calculated, they are mapped to a large-scale model for the unified inversion.
Due to the reduced spatial scale of the computational model, the efficiency of 3D AEM inversions can be significantly improved. Liu and Yin [13] further used this technique to realize the 3D inversion of multi-pulse TD AEM data. To overcome the difficulties imposed by using multi-source and multi-channel of TD AEM systems, Oldenburg et al. [14] adopted a direct-solver to quickly solve AEM responses that greatly improves the efficiency of 3D inversions. Based on this research of [14], Yang et al. [15] further proposed the local mesh method that effectively reduces the model dimension and accelerates 3D inversion for large-scale models.
Although the 3D inversions of TD AEM data have become increasingly mature, the inversion models are still based on the assumption of isotropic media. Many studies have shown that anisotropy of subsurface media has serious impacts on the EM data [16][17][18][19][20]. Therefore, anisotropy has become an important research focus for many EM methods, such as in marine controlled-source EM (MCSEM) [21][22][23][24], well-logging [25][26][27][28][29], long-offset transient EM (LOTEM) [30], etc. These anisotropic inversion methods have demonstrated that they can provide more reliable results. However, the effect of electrical anisotropy on AEM methods has been less studied, except for [31][32][33] who discussed the effect of anisotropy on AEM responses. Their results show that using an isotropic model to interpret anisotropic data will deliver erroneous results. Further studies show that AEM signal is affected by earth anisotropy in a quite specific way. If the underground structure is 1D, the source field is purely TE-mode and the underground current mainly flows in the horizontal direction that is insensitive to the vertical conductivity, thus only the lateral resistivity will affect the AEM response. The presence of 3D anomalies will induce vertical current around them but it is usually much weaker than the horizontal one, so that only a small amount of information about the underground electrical structure can be observed.
To improve the accuracy of TD AEM data interpretation, we propose in this paper a 3D inversion algorithm for AEM data inversion in TD using a triaxial anisotropic model. We use the 3D finite-difference algorithm based on the secondary scattered electric field approach [13] and the moving footprint technique [12] to simulate large-scale AEM problems. In the inversion process, we use the adjoint forward modeling of [13] and GN optimization to recover the resistivity of the earth. We test our algorithm on synthetic AEM data and field data acquired from Musgraves, South Australia.

Governing Equations
We first set up a Cartesian coordinate system in which the vertical z-axis is positive downward. A time harmonic dependence e iωt is assumed, where i = √ −1, ω = 2πf, the second order of partial differential equation for the scattered electric field in a 3-axial anisotropic earth can be expressed as [34] ∇ × ∇ × E s + iωµ 0 σE s = −iωµ 0 J p (1) where µ 0 = 4π × 10 −7 H/m is the magnetic permeability of the free space, E s is the scattered electric field, J p is the effective source current densities for the scattered field that is expressed as J p = (σ − σ p )E p . σ is the model conductivity, while σ p is the background conductivity. For the 3-axial anisotropic case, the σ and σ p are both diagonal conductivity tensors composed of three principal conductivities, i.e., Minerals 2021, 11, 218 3 of 17 E p is the primary electric field induced by the current in the background structures. The total electric field is E = E s + E p . The magnetic field can be obtained by the Faraday's law B = (−iω) −1 ∇ × E. Expansion of the curl operation in Equation (1) yields Using the finite-difference formulation to replace the derivatives in Equations (3)-(5), we obtain the final equations system for the secondary electrical field, i.e., where K is a large sparse matrix and s is the source term. The magnetic field is calculated by B = LE, where L is a spatial interpolation operator. The AEM modeling for a fully anisotropic earth based on the finite-difference formulation of Equation (1) can be found in [33]. In our presented anisotropic inversion, we only focus on 3-axial anisotropic cases, so we use a simplified version of the code developed by [33] to solve the forward problem. After obtaining the FD AEM responses, we adopt inverse Fourier transform to calculate the TD AEM step responses, i.e., where F(ω) are the FD responses and B s (t) is the TD step responses. Then, the TD AEM responses for arbitrary transmitting wave can be calculated via the convolution [35].
where I(t) is the transmitting current and dB(t)/dt is the magnetic induction.

Regularized Inversion
The objective functional ϕ(m, d) for our anisotropic inversion is defined as where m is a 3 × M matrix for model parameters containing three principal conductivities, F is the forward modeling operator. F(m) denotes the synthetic data vector calculated from the model m, while d denotes the AEM field data. C d and C m are the data covariance matrix and model roughness operator, respectively. λ is the regularization parameter used to balance the relative contributions of data misfit ϕ d and model smoothness ϕ m . Assuming that m n is the model parameters at the nth iteration, J is the sensitivity matrix evaluated at m n , and r = d − F(m n ) is the data residual. Then, the GN equation derived from objective functional in Equation (9) can be written as which can be solved by preconditioned conjugate gradient (PCG) method for the model update ∆m. The model parameters are updated by m n+1 = m n + α∆m, where α is determined by a linear search [36].

Calculation of Sensitivity Matrix for 3-Axial Anisotropic Media
Since the footprint of AEM systems is much smaller than the survey area, a small number of meshes is required for numerical modeling. To reduce the computational cost we introduce the footprint concept into our AEM modeling and inversions. In the following, we first divide the entire model into N m sub-models based on the AEM footprint and calculate the sensitivity matrix for each sub-model. Then, we rearrange these sensitivity matrices into the entire sensitivity matrix. In the Appendix A, we will show how to divide the model domain into sub-models. For the nth sub-model, following [21] on their transversely anisotropic inversion, we derivate both side of Equation (6) with respect to parameters m x k , m y k , m z k of the kth cell and get Considering that B = LE and ∂E p /∂m x,y,z k = 0, we obtain the elements for the sensitivity matrixĴ Assuming a N edge × 3M matrix where N edge is the number of edges, then the sensitivity matrix can be expressed aŝ Let v = K −1 L T , and multiply both sides by matrix K, we get the final equation for the adjoint forward modeling, i.e., Kv = L T From Equation (16) we can see that no more work needs to be done for calculating the sensitivity matrix than in an isotropic case. The difference between anisotropic and isotropic inversions is the formation of matrix G, which can be done analytically. There are in total N d columns in v and L T , so we select the same column of v and L T each time to set up an equation, and then solve the equation to obtain the corresponding column of v. After getting the matrix v, we use Equation (15) to calculateĴ T n , and resemble eachĴ T n for each sub-model to obtain the entire sensitivity matrix in Equation (10). After finishing the calculation in the FD, we use the inverse Fourier transform and the convolution with the transmiting wave to get the sensitivity matrix for specific waveforms in TD. The workflow for TD AEM inversion is given in Figure 1.

Numerical Experiments
In the synthetic data inversion, we use the transmitting wave of the VTEM (versatile time-domain electromagnetic) system from Geotech company (Aurora, ON, Canada) (c.f. Figure 2). The transmitter and receiver assumes a concentric geometry. The height of the transmitting coil is 30 m and the receiver is 0.5 m above the transmitter. To establish the first model for testing, we assume a tilted ore-body buried in a half-space with a distinct strike as shown in Figure 3a. The ore-body is 60 meter wide in the x-direction, 360 m long in the y-direction, the top surface is located at z = 10 m, and the bottom surface is located at z = 208 m. The entire model is discretized into 48 × 48 × 13 cells (refer to Appendix for specific grid design). Figure 3b is the distribution of the measuring points. The distance between measuring points is 40 m. For easy understanding, we use the resistivity instead of the conductivity in our following discussions. For synthetic and field survey data, the initial and reference model are both selected as isotropic half-space of 100 ohm-m to implement our anisotropic and isotropic inversions (for the anisotropic cases, we assume the same initial resistivities along three principal axes). The data misfit is calculated by the root mean square (RMS), i.e., where di, obs and di,calc are respectively the observed and predicted data, εi, obs is the estimate of data error, N is the number of data. We use 3% of the syntetic data as the data error estimate in our isotropic and anisotropic inversions. In the iterative inversion process, the

Numerical Experiments
In the synthetic data inversion, we use the transmitting wave of the VTEM (versatile time-domain electromagnetic) system from Geotech company (Aurora, ON, Canada) (c.f. Figure 2). The transmitter and receiver assumes a concentric geometry. The height of the transmitting coil is 30 m and the receiver is 0.5 m above the transmitter. To establish the first model for testing, we assume a tilted ore-body buried in a half-space with a distinct strike as shown in Figure 3a. The ore-body is 60 m wide in the x-direction, 360 m long in the y-direction, the top surface is located at z = 10 m, and the bottom surface is located at z = 208 m. The entire model is discretized into 48 × 48 × 13 cells (refer to Appendix A for specific grid design). Figure 3b is the distribution of the measuring points. The distance between measuring points is 40 m. For easy understanding, we use the resistivity instead of the conductivity in our following discussions. For synthetic and field survey data, the initial and reference model are both selected as isotropic half-space of 100 ohm-m to implement our anisotropic and isotropic inversions (for the anisotropic cases, we assume the same initial resistivities along three principal axes). The data misfit is calculated by the root mean square (RMS), i.e., where d i, obs and d i,calc are respectively the observed and predicted data, ε i, obs is the estimate of data error, N is the number of data. We use 3% of the syntetic data as the data error estimate in our isotropic and anisotropic inversions. In the iterative inversion process, the initial regularization parameter is set to be 1000. When the RMS difference between the last two iterations becomes less than 0.02, the regularization parameter is reduced to 1/10 of the previous one. When the RMS falls below 1 and the RMS difference between the last two iterations becomes less than 0.002, or the iterations reach the maximum number, the inversion is terminated.
Minerals 2021, 11, 218 6 of 1 initial regularization parameter is set to be 1000. When the RMS difference between th last two iterations becomes less than 0.02, the regularization parameter is reduced to 1/1 of the previous one. When the RMS falls below 1 and the RMS difference between the las two iterations becomes less than 0.002, or the iterations reach the maximum number, th inversion is terminated.

Synthetic Data Inversion for an Anisotropic Ore-Body
For the first synthetic example, we assume an isotropic background of 100 ohm-m and an anisotropic ore-body (c.f. Figure 3a) with the resistivity tensor of m ohm-20 0 0 Figure 4 shows the inversion results in slices in the xz-section for three principal re sistivities. It is seen that the resistivities in the x-and y-directions are well recovered, bu the resistivity in the z-direction is not recovered. This can be easily explained by the char acteristic EM wave induced by a TD AEM system. If the earth is 1D, the source field i purely TE-mode and the underground current mainly flows in the horizontal direction thus the EM field is only sensitive to the resistivity in the x-and y-directions, while it i insensitive to the resistivity in the z-direction for the vertical resistivity ρz has little effec on AEM responses. This explains that although the model is 3D, the current densities i the x-and y-direction are much larger than that in the z-direction. This results in the in sensitive EM responses and poor resolvability to the resistivity in the z-direction. Figure 4d presents the isotropic inversions for the same dataset. Since the syntheti data in this example is calculated from an anisotropic model, it is difficult to fit the dat by an isotropic model. Thus, artifacts are introduced into the inversion results. This mean that when the data are collected in an area with anisotropy, an anisotropic inversion i initial regularization parameter is set to be 1000. When the RMS difference between the last two iterations becomes less than 0.02, the regularization parameter is reduced to 1/10 of the previous one. When the RMS falls below 1 and the RMS difference between the last two iterations becomes less than 0.002, or the iterations reach the maximum number, the inversion is terminated.

Synthetic Data Inversion for an Anisotropic Ore-Body
For the first synthetic example, we assume an isotropic background of 100 ohm-m and an anisotropic ore-body (c.f. Figure 3a) with the resistivity tensor of m ohm-20 0 0 Figure 4 shows the inversion results in slices in the xz-section for three principal resistivities. It is seen that the resistivities in the x-and y-directions are well recovered, but the resistivity in the z-direction is not recovered. This can be easily explained by the characteristic EM wave induced by a TD AEM system. If the earth is 1D, the source field is purely TE-mode and the underground current mainly flows in the horizontal direction, thus the EM field is only sensitive to the resistivity in the x-and y-directions, while it is insensitive to the resistivity in the z-direction for the vertical resistivity ρz has little effect on AEM responses. This explains that although the model is 3D, the current densities in the x-and y-direction are much larger than that in the z-direction. This results in the insensitive EM responses and poor resolvability to the resistivity in the z-direction. Figure 4d presents the isotropic inversions for the same dataset. Since the synthetic data in this example is calculated from an anisotropic model, it is difficult to fit the data by an isotropic model. Thus, artifacts are introduced into the inversion results. This means that when the data are collected in an area with anisotropy, an anisotropic inversion is

Synthetic Data Inversion for an Anisotropic Ore-Body
For the first synthetic example, we assume an isotropic background of 100 ohm-m and an anisotropic ore-body (c.f. Figure 3a) with the resistivity tensor of Figure 4 shows the inversion results in slices in the xz-section for three principal resistivities. It is seen that the resistivities in the xand y-directions are well recovered, but the resistivity in the z-direction is not recovered. This can be easily explained by the characteristic EM wave induced by a TD AEM system. If the earth is 1D, the source field is purely TE-mode and the underground current mainly flows in the horizontal direction, thus the EM field is only sensitive to the resistivity in the xand y-directions, while it is insensitive to the resistivity in the z-direction for the vertical resistivity ρ z has little effect on AEM responses. This explains that although the model is 3D, the current densities in the xand y-direction are much larger than that in the z-direction. This results in the insensitive EM responses and poor resolvability to the resistivity in the z-direction. Figure 4d presents the isotropic inversions for the same dataset. Since the synthetic data in this example is calculated from an anisotropic model, it is difficult to fit the data by an isotropic model. Thus, artifacts are introduced into the inversion results. This means that when the data are collected in an area with anisotropy, an anisotropic inversion is necessary. To better locate the ore-body and eliminate the false resistivity distribution in surrounding rocks, we calculate an equivalent resistivity in the horizontal direction using the inversion results of ρ x and ρ y via the following formula [37] From Figure 4h, we can see that the equivalent resistivity effectively suppresses the fake anomalies brought by the "over anisotropy" phenomenon.
From Figure 4h, we can see that the equivalent resistivity effectively suppresses th fake anomalies brought by the "over anisotropy" phenomenon.   Figure 3b). The data fittings for isotropic and anisotropic inversions are both very well, the RMSs are at the same level. Thus, for aniso tropic data, using isotropic model can also achieve a good data fitting, but the under ground model is not in line with the actual situation. Thus, it is necessary to implemen anisotropic inversion when the earth is anisotropic. From Figure 6a, we can see that both isotropic and anisotropic inversions achieve the same level of RMS with stable conver gence. The regularization parameter decreases from the initial 1000 to 0.00001 (Figure 6b)   Figure 3b). The data fittings for isotropic and anisotropic inversions are both very well, the RMSs are at the same level. Thus, for anisotropic data, using isotropic model can also achieve a good data fitting, but the underground model is not in line with the actual situation. Thus, it is necessary to implement anisotropic inversion when the earth is anisotropic. From Figure 6a, we can see that both isotropic and anisotropic inversions achieve the same level of RMS with stable convergence. The regularization parameter decreases from the initial 1000 to 0.00001 (Figure 6b).

Synthetic Data Inversion for Anisotropic Host Rock
In the second synthetic data inversion, we assume that the first layer has a thickness of 100 m and a resistivity of the second layer is an isotropic half-space with a resistivity of 100 ohm-m. The isotropic anomalous body has a resistivity of 10 ohm-m. Figure 7 shows the inversion results obtained from the synthetic VTEM data. It is seen from the figures that the anisotropic inversion well recovers the structure and resistivity values of the anomalous body in the xand y-directions. However, like the previous example, the resistivity in the z-direction is not recovered accurately. From Figure 7, we can also see that the anisotropic character of the first layer is not well recovered. The main reason is that we assume a concentric TD AEM system for our synthetic data inversion. If the ore-body is absent, the model can be viewed as a horizontally uniform anisotropic medium and the same response can still be obtained by swapping the resistivity in the x-and y-directions. The presence of a 3D orebody results in the generation of a vertical current, but the resolution of resistivity in the

Synthetic Data Inversion for Anisotropic Host Rock
In the second synthetic data inversion, we assume that the first layer has a thickness of 100 m and a resistivity of the second layer is an isotropic half-space with a resistivity of 100 ohm-m. The isotropic anomalous body has a resistivity of 10 ohm-m. Figure 7 shows the inversion results obtained from the synthetic VTEM data. It is seen from the figures that the anisotropic inversion well recovers the structure and resistivity values of the anomalous body in the xand y-directions. However, like the previous example, the resistivity in the z-direction is not recovered accurately. From Figure 7, we can also see that the anisotropic character of the first layer is not well recovered. The main reason is that we assume a concentric TD AEM system for our synthetic data inversion. If the ore-body is absent, the model can be viewed as a horizontally uniform anisotropic medium and the same response can still be obtained by swapping the resistivity in the x-and y-directions. The presence of a 3D orebody results in the generation of a vertical current, but the resolution of resistivity in the

Synthetic Data Inversion for Anisotropic Host Rock
In the second synthetic data inversion, we assume that the first layer has a thickness of 100 m and a resistivity of the second layer is an isotropic half-space with a resistivity of 100 ohm-m. The isotropic anomalous body has a resistivity of 10 ohm-m. Figure 7 shows the inversion results obtained from the synthetic VTEM data. It is seen from the figures that the anisotropic inversion well recovers the structure and resistivity values of the anomalous body in the xand y-directions. However, like the previous example, the resistivity in the z-direction is not recovered accurately. From Figure 7, we can also see that the anisotropic character of the first layer is not well recovered. The main reason is that we assume a concentric TD AEM system for our synthetic data inversion. If the ore-body is absent, the model can be viewed as a horizontally uniform anisotropic medium and the same response can still be obtained by swapping the resistivity in the xand y-directions. The presence of a 3D ore-body results in the generation of a vertical current, but the resolution of resistivity in the xand y-directions in the inversion does not increase much (see Figure 7e-f). From the inversion results in xand y-direction in Figure 7, we find that in the first layer, the x-direction resistivity has a fake conductive anomaly on the left side, while the y-direction resistivity at the same location is resistive, which means that we find a "over anisotropy" phenomenon in the inversion. Figure 8 shows the horizontal slices of the anisotropic inversion results of ρ x and ρ y and the equivalent resistivity.
x-and y-directions in the inversion does not increase much (see Figure 7e-f). From the inversion results in x-and y-direction in Figure 7, we find that in the first layer, the xdirection resistivity has a fake conductive anomaly on the left side, while the y-direction resistivity at the same location is resistive, which means that we find a "over anisotropy" phenomenon in the inversion. Figure 8 shows the horizontal slices of the anisotropic inversion results of ρx and ρy and the equivalent resistivity. From Figures 7h and 8c, we can see that the equivalent resistivity effectively suppresses the fake anomalies brought by the "over anisotropy". Although the definition of equivalent resistivity has little exact physical meaning, we can get from its distribution a more reliable model structure. Similarly, from the isotropic inversion results given in Figure 7d, we can draw the same conclusions as in the previous example. The resistivity distribution in the isotropic inversion is relatively messy, and the resistivity varies seriously. The reason for these is that to fit the data, many false anomalies are created to fit the anisotropic structure.  From Figures 7h and 8c, we can see that the equivalent resistivity effectively suppresses the fake anomalies brought by the "over anisotropy". Although the definition of equivalent resistivity has little exact physical meaning, we can get from its distribution a more reliable model structure. Similarly, from the isotropic inversion results given in Figure 7d, we can draw the same conclusions as in the previous example. The resistivity distribution in the isotropic inversion is relatively messy, and the resistivity varies seriously. The reason for these is that to fit the data, many false anomalies are created to fit the anisotropic structure.
From Figure 9, we can see that both isotropic and anisotropic inversions have good data fitting. From the RMS curves given in Figure 10a, we can see that both isotropic and anisotropic inversions achieve stable convergence. From Figure 10b, the regularization parameter for the isotropic inversion decreases from the initial value of 1000 to 0.00001, while the regularization parameter for the anisotropic inversion decreases from 1000 to 1. From Figure 9, we can see that both isotropic and anisotropic inversions have good data fitting. From the RMS curves given in Figure 10a, we can see that both isotropic and anisotropic inversions achieve stable convergence. From Figure 10b, the regularization parameter for the isotropic inversion decreases from the initial value of 1000 to 0.00001, while the regularization parameter for the anisotropic inversion decreases from 1000 to 1.    From Figure 9, we can see that both isotropic and anisotropic inversions have g data fitting. From the RMS curves given in Figure 10a, we can see that both isotropic anisotropic inversions achieve stable convergence. From Figure 10b, the regulariza parameter for the isotropic inversion decreases from the initial value of 1000 to 0.00 while the regularization parameter for the anisotropic inversion decreases from 1000 t   From Figure 9, we can see that both isotropic and anisotropic inversions have good data fitting. From the RMS curves given in Figure 10a, we can see that both isotropic and anisotropic inversions achieve stable convergence. From Figure 10b, the regularization parameter for the isotropic inversion decreases from the initial value of 1000 to 0.00001, while the regularization parameter for the anisotropic inversion decreases from 1000 to 1.

Field Data Example
To further demonstrate the effectiveness of our anisotropic inversions, we invert a survey dataset acquired by the VTEM system in Musgrave province, South Australia ( Figure 11). This dataset was acquired for mineral prospecting and regional hydro-geological assessment. We select this dataset because this area presents abrupt lateral variations, product of an outcropping mafic intrusive and more subtle lateral continuous sedimentary sequences [38]. The subsurface structure in this area have obvious anisotropy.
For our inversion, we select 210 surveying points along 7 flight lines, marked as white lines in Figure 11. The inverse model is discretized into 178 × 48 × 19 cells. In the central area, the grid sizes in the x-and y-directions are 20 m and 40 m, respectively, while the dimensions of the 4 external grid sizes increase by a factor of 2, symmetrically to both sides of the x-and y-directions. In the z-direction, the thickness of first layer is 10 m, and the thickness of the underlying layers increase by a factor of 1.2 until the 13th layer, the 14th-19th layers are extended layers with thicknesses of 100, 100, 100, 200, 300, and 400 m, respectively. We add 8 external grids in the air layer. The moving footprint is used to divide the entire model into sub-models. The division of sub-model is the same as that of synthetic examples. In the inversion process, we reduce the regularization parameter to 1/10 of the previous value in each iteration and keep it unchanged when it reaches the threshold of 0.001, which can ensure the stable convergence in the inversion. The inversion is terminated when the RMS difference between the last two iterations is less than a threshold of 0.03. We use 3% of the field data as the data error estimate in our isotropic and anisotropic inversions. Figure 11. Flight lines in Musgraves of South Australia. The TD AEM data are acquired by the VTEM system (modified from [38]). Figure 12 presents the comparison of 3D anisotropic and isotropic inversion results. It can be seen from the figures that the resistivity distributions of the anisotropic and isotropic inversions are similar. They show a layered structure, which is consistent with the Figure 11. Flight lines in Musgraves of South Australia. The TD AEM data are acquired by the VTEM system (modified from [38]).
For our inversion, we select 210 surveying points along 7 flight lines, marked as white lines in Figure 11. The inverse model is discretized into 178 × 48 × 19 cells. In the central area, the grid sizes in the xand y-directions are 20 m and 40 m, respectively, while the dimensions of the 4 external grid sizes increase by a factor of 2, symmetrically to both sides of the xand y-directions. In the z-direction, the thickness of first layer is 10 m, and the thickness of the underlying layers increase by a factor of 1.2 until the 13th layer, the 14th-19th layers are extended layers with thicknesses of 100, 100, 100, 200, 300, and 400 m, respectively. We add 8 external grids in the air layer. The moving footprint is used to divide the entire model into sub-models. The division of sub-model is the same as that of synthetic examples. In the inversion process, we reduce the regularization parameter to 1/10 of the previous value in each iteration and keep it unchanged when it reaches the threshold of 0.001, which can ensure the stable convergence in the inversion. The inversion is terminated when the RMS difference between the last two iterations is less than a threshold of 0.03. We use 3% of the field data as the data error estimate in our isotropic and anisotropic inversions. Figure 12 presents the comparison of 3D anisotropic and isotropic inversion results. It can be seen from the figures that the resistivity distributions of the anisotropic and isotropic inversions are similar. They show a layered structure, which is consistent with the sedimentary lithology in this area. However, a detailed analysis shows that there exist significant differences in details. First, the anisotropic inversions show a smoother and more continuous layered structure than the isotropic ones. In opposite, the isotropic inversions deliver random results among these lines. The continuity is poor and there is a thick lowresistivity layer in the central region of Line 5, which is quite different from the surrounding lines. Second, in the near surface between x = 1500−2500 m, the anisotropic inversion shows a clear high-resistivity layer in the y-direction, but from the isotropic inversion we cannot see this layer clearly. By comparing Figure 12a,b, we can further see that there exists a big difference between the resistivities ρ x and ρ y from x = 1000 to 2500 m. However, this cannot be seen from the isotropic inversion results. The artifacts appearing in the anisotropic and isotropic inversions at the beginning of the profile may result from the poor data fitting. sedimentary lithology in this area. However, a detailed analysis shows that there exist significant differences in details. First, the anisotropic inversions show a smoother and more continuous layered structure than the isotropic ones. In opposite, the isotropic inversions deliver random results among these lines. The continuity is poor and there is a thick low-resistivity layer in the central region of Line 5, which is quite different from the surrounding lines. Second, in the near surface between x = 1500−2500 m, the anisotropic inversion shows a clear high-resistivity layer in the y-direction, but from the isotropic inversion we cannot see this layer clearly. By comparing Figure 12a,b, we can further see that there exists a big difference between the resistivities ρx and ρy from x = 1000 to 2500 m. However, this cannot be seen from the isotropic inversion results. The artifacts appearing in the anisotropic and isotropic inversions at the beginning of the profile may result from the poor data fitting.  We take as example the horizontal slices of isotropic and anisotropic inversions at the depth of 44m and show them in Figure 13. It is seen that the inverted anisotropic models of ρx in Figure 13a and ρy in 13b are quite different from each other between x=1500 and x=2500m, showing the anisotropic characteristics in this area. The equivalent resistivity in Figure 13d is similar to the isotropic inversion results in Figure 13c. This is due to the horizontal current in the earth excited by an AEM system in TD, the EM signal is mainly We take as example the horizontal slices of isotropic and anisotropic inversions at the depth of 44m and show them in Figure 13. It is seen that the inverted anisotropic models of ρ x in Figure 13a and ρ y in 13b are quite different from each other between x = 1500 and x = 2500m, showing the anisotropic characteristics in this area. The equivalent resistivity in Figure 13d is similar to the isotropic inversion results in Figure 13c. This is due to the horizontal current in the earth excited by an AEM system in TD, the EM signal is mainly created by the averaged resistivities in the xand y-directions, the equivalent resistivity takes most effect, so we can see a coincidence between Figure 13c,d.
example. Figure 14 shows that these two inversions have the same level of data fitting a overall RMS. However, the data fitting around 1000 m is not good, which is coincide with the strong resistivity contrast in the inversion results (c.f. Figures 12 and 13). Fig  15 shows the RMSs and regularization parameters for the isotropic and anisotropic inv sions. Both inversions have stable convergence and good data fitting. The anisotropic version converges faster than the isotropic one. From Figure 15b, the regularization rameter of each iteration decreases by 1/10 and remains unchanged when it reaches 0.0   To check the data fitting of anisotropic and isotropic inversions, we take Line 5 as an example. Figure 14 shows that these two inversions have the same level of data fitting and overall RMS. However, the data fitting around 1000 m is not good, which is coincidence with the strong resistivity contrast in the inversion results (c.f. Figures 12 and 13). Figure 15 shows the RMSs and regularization parameters for the isotropic and anisotropic inversions. Both inversions have stable convergence and good data fitting. The anisotropic inversion converges faster than the isotropic one. From Figure 15b, the regularization parameter of each iteration decreases by 1/10 and remains unchanged when it reaches 0.001.
Minerals 2021, 11,218 13 of 17 created by the averaged resistivities in the x-and y-directions, the equivalent resistivity takes most effect, so we can see a coincidence between Figure 13c,d.
To check the data fitting of anisotropic and isotropic inversions, we take Line 5 as an example. Figure 14 shows that these two inversions have the same level of data fitting and overall RMS. However, the data fitting around 1000 m is not good, which is coincidence with the strong resistivity contrast in the inversion results (c.f. Figures 12 and 13). Figure  15 shows the RMSs and regularization parameters for the isotropic and anisotropic inversions. Both inversions have stable convergence and good data fitting. The anisotropic inversion converges faster than the isotropic one. From Figure 15b, the regularization parameter of each iteration decreases by 1/10 and remains unchanged when it reaches 0.001.

Conclusions
In this paper, we have successfully applied the anisotropic model to TD AEM data inversion. The numerical examples for both synthetic and field data proved its effectiveness. From the inversions to the synthetic models we find that it is difficult to invert the resistivity in the vertical z-direction from TD AEM data because for an AEM system the induced current in z-direction is much smaller than those in the horizontal x-and y-directions. Another finding from our synthetic examples is that the inversions based on an ansiotropic model may create "over anisotropy", which, however, can be supressed by the equivalent resistivity calculated from the anisotropic inversion results. Thus, by combining all those inversion results together, we can obtain a reasonable resistivity distribution inside the earth. Further numerical experiment with field survey data demonstrated that compared to the isotropic inversion, the inversion based on an anisotropic model can provide more genuine results than the isotropic one for the TD AEM data collected in an anisotropic area.
Author Contributions: Conceptualization, Y.S. and C.Y.; formal analysis, C.Y. and Y.L.; funding acquisition, C.Y.; investigation, X.R.; methodology, Y.S. and Y.L.; project administration, C.Y.; resources, X.R. and Y.L.; software, Y.L. and Y.S.; supervision, B.Z. and X.R.; validation, B.Z. and B.X.; visualization, Y.L. and Y.S.; writing-original draft, Y.S. and Y.L. All authors have read and agreed to the published version of the manuscript. Acknowledgments: Thanks to Alan Yusen Ley-Cooper for providing the field data. We want to thank three reviewers and editors for their constructive comments and suggestions that improve the clarity of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

A.1. Sub-Model Based on the Moving Footprint
For our forward modeling and calculation of the sensitivity matrix based on the "moving footprint", we extract an area (200 m × 200 m) under each survey site to construct the central area of the sub-model ( Figure A1), and then add 4 external grids in the x-and y-directions with the sizes of 40, 80, 160, and 320m, symmetrically to both sides. We use

Conclusions
In this paper, we have successfully applied the anisotropic model to TD AEM data inversion. The numerical examples for both synthetic and field data proved its effectiveness. From the inversions to the synthetic models we find that it is difficult to invert the resistivity in the vertical z-direction from TD AEM data because for an AEM system the induced current in z-direction is much smaller than those in the horizontal xand y-directions. Another finding from our synthetic examples is that the inversions based on an ansiotropic model may create "over anisotropy", which, however, can be supressed by the equivalent resistivity calculated from the anisotropic inversion results. Thus, by combining all those inversion results together, we can obtain a reasonable resistivity distribution inside the earth. Further numerical experiment with field survey data demonstrated that compared to the isotropic inversion, the inversion based on an anisotropic model can provide more genuine results than the isotropic one for the TD AEM data collected in an anisotropic area. Acknowledgments: Thanks to Alan Yusen Ley-Cooper for providing the field data. We want to thank three reviewers and editors for their constructive comments and suggestions that improve the clarity of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.  Figure A1), and then add 4 external grids in the xand y-directions with the sizes of 40, 80, 160, and 320 m, symmetrically to both sides. We use the volume-averaging technique to calculate the resistivity in the external region. The sub-model area for calculating the sensitivity at each survey site is as large as the footprint, so that there are only few non-zero elements in a row of the overall matrix. The global sensitivity matrix is formed by putting all the sensitivity matrices together. Different from the conventional sensitivity matrix, the sensitivity matrix calculated based on the footprint is sparse and easy to store. The reduction of the model scale greatly improves the model efficiency.
Minerals 2021, 11,218 15 of 17 the volume-averaging technique to calculate the resistivity in the external region. The submodel area for calculating the sensitivity at each survey site is as large as the footprint, so that there are only few non-zero elements in a row of the overall matrix. The global sensitivity matrix is formed by putting all the sensitivity matrices together. Different from the conventional sensitivity matrix, the sensitivity matrix calculated based on the footprint is sparse and easy to store. The reduction of the model scale greatly improves the model efficiency. Figure A1. Sub-model based on the moving footprint.

A.2. Grid Design for Our Synthetic Models
For synthetic modelings, we design the grids as follows. The entire model is discretized into 48 × 48 × 13 cells. The grid sizes in the central area in the x-and y-directions are assumed to be 20 m, and the dimensions of the 4 external grids in the x-and y-directions are 40, 80, 160, and 320 m, symmetrically to both sides. The total extensions in the x-and y-direction are 2000 m. In the z-direction, the thickness of the first layer is 10 m, and the thickness of the underlying layers increases sequentially by a factor of 1.2 until the 9th layer, the 10th-13th layers are extended layers with thicknesses of 100, 200, 300, and 400 m, respectively. The total extension in the z-direction is 1208 m. We add 8 external grids in the air layer. The first two grids in the air have a size of 10 m, and the rest are increased by a factor of 3. The total extension of the air layer in the z-direction is 32,800 m. The final model grids are shown in Figure A2.

Appendix A.2 Grid Design for Our Synthetic Models
For synthetic modelings, we design the grids as follows. The entire model is discretized into 48 × 48 × 13 cells. The grid sizes in the central area in the xand y-directions are assumed to be 20 m, and the dimensions of the 4 external grids in the xand y-directions are 40, 80, 160, and 320 m, symmetrically to both sides. The total extensions in the xand y-direction are 2000 m. In the z-direction, the thickness of the first layer is 10 m, and the thickness of the underlying layers increases sequentially by a factor of 1.2 until the 9th layer, the 10th-13th layers are extended layers with thicknesses of 100, 200, 300, and 400 m, respectively. The total extension in the z-direction is 1208 m. We add 8 external grids in the air layer. The first two grids in the air have a size of 10 m, and the rest are increased by a factor of 3. The total extension of the air layer in the z-direction is 32,800 m. The final model grids are shown in Figure A2.