A finitemdifference algorithm for the transient electromagnetic response of a three-dimensional body

SUMMARY We attempted to develop a direct time-domain, finite-difference solution for the electromagnetic response of a 3-D model. The algorithm is an extension of our 2-D modelling technique, which uses the Du Fort-Frankel finite-difference scheme. However, the vector nature of the field makes the 3-D problem much more complicated than its 2-D counterpart, and a supercomputer is required for computations. Unlike the 2-D case, where we solve for the electric field, the solution is formulated in terms of secondary magnetic field to avoid dealing with a discontinuous normal component of electric field. However, that difficulty is replaced with a problem involving the gradient of conductivity, which is discontinuous at interfaces. We experimented with both smoothing the conductivity variation, so that the gradient is well defined; and integrating the gradient terms, which results in a tangential current density contrast. Our limited experiments indicate that the magnetic field step response is computationally more stable than the impulse response, because it is a smoother function with smaller dynamic range. However, the step response takes longer to compute because its source, which involves the primary fields in the earth, requires very accurate numerical integration. Due to computer time and memory limitations on the supercomputer available to us, we were not able to develop an accurate numerical solution. We were able to carry out only a few tests on small models for which results were not in good agreement with values computed using a 3-D integral equation solution.


INTRODUCTION
In many cases a numerical solution for a 3-D model is required for interpreting transient electromagnetic ( E M ) data in terms of realistic Earth structures. Although a 2-D model can be useful, its application is limited to special cases, and its response can be very different from its 3-D counterpart (Adhidjaja & Hohmann 1988). However, because there are conductivity boundaries in all directions and the unknown field is a vector quantity, the 3-D problem is much more difficult.
Methods that have been used for formulating numerical solutions are the differential equation (finite element and finite difference), integral equation, and hybrid techniques. Differential equations are easier to implement and are more suitable for modelling complex geology, but they require more computer storage. In contrast, integral-equation formulations involve more difficult mathematics, but the unknown fields only need to be found in anomalous regions. Thus, integral-equation solutions are less expensive for * Formerly Department of Geology and Geophysics, University of Utah.
calculating the response of one or a few small bodies (Hohmann 1983). Hybrid techniques are attractive because they take the best features from both the differential-and the integral-equation methods, but they have their own difficulties (Lee, Pridmore & Morrison 1981;Petrick 1984;Best et al. 1985). A complete review of those methods and assessments of some of the solutions are given by Hohmann (1983).
Most EM solutions in the geophysical literature have been formulated in the frequency domain. Integral-equation solutions were developed, for example by Hohmann (1975), Weidelt (1975), Raiche (1974), and Das & Verma (1982). Differential equations were solved by Lines & Jones (1973) with a finite-difference method, and by Reddy, Rankin & Phillips (1977) with a finite-element method applied to geomagnetic-field and magnetotelluric problems, respectively, which typically involve measurements at low frequencies. Zhdanov et al. (1982) presented an algorithm for computing the responses of 2-D and 3-D models with a finite-difference method, but showed no 3-D results. Pridmore et al. (1981) used a finite-element method for direct current (DC) and controlled-source EM methods applied to a 3-D model in a mining problem.
Two approaches have been used to obtain time-domain solutions. The indirect method uses Fourier transformation of frequency-domain results (Lamontagne 1975;Tripp 1982;Newman, Hohmann & Anderson 1986). One advantage of the indirect method is that it can be based on an existing and well-tested frequency-domain computer program. The direct method requires solving an initial value problem directly in the time domain, which can be done via a time-stepping method.
The direct method may provide more insight into the transient (TEM) process within the whole spectrum, from very early to very late times. However, the early-time solution must be calculated accurately to insure correct late-time results. For a differential-equation method the computation can be very time consuming. San Filipo  developed an efficient direct time-domain solution based on an integral equation formulation, but their solution was for a simple prism in an otherwise homogeneous half-space. Yee (1966) presented a finite-difference scheme for solving the 3-D transient electromagnetic problem; it was based on time-stepping Maxwell's equations and solving for both electric and magnetic fields. His scheme has been applied successfully to compute scattering by a perfect conductor. Success seems to be due to an accurate representation of the boundary conditions of fields at interfaces, which is achieved by carefully placing grid points. The scheme utilizes a staggered grid; the electric and magnetic fields at a node are defined at alternate time steps. Later, this scheme was adapted and used successfully for solving various 3-D electromagnetic scattering problems in free-space (Holland 1977;Kunz & Lee 1978a,b). The scheme has also been tested for a geophysical problem by Bergeal (1982). He solved for a 3-D body in a whole-space, but his results, computed for only lops, were not satisfactory. The main difficulty involved propagating the electric field via the displacement current, which we normally neglect in geophysical problems. Including the displacement current requires a very small time step, which is impractical in computing a solution to the later times of interest in geophysics.
Our goal was to develop a direct time-domain solution based on a differential equation for simulating a complex 3-D model with realistic conductivity contrasts. We formulated the solution in terms of the secondary field, based on experience with a successful algorithm for the 2-D problem (Oristaglio & Hohmann 1984;Adhidjaja, Hohmann & Oristaglio 1985). In this paper, we present the theoretical formulation, the finite-difference scheme, the limited numerical results we were able to compute, and discussions of the problems encountered. Following Hohmann (1983), we begin with Maxwell's equations,

THEORETICAL FORMULATION
with magnetic permeability p = po everywhere and with displacement current neglected. Here mP is the dipole moment per volume of an impressed source. In our case we sum vertical magnetic dipoles to create a horizontal loop source.
As in the 2-D case we formulate the problem in terms of primary and secondary fields through the decompositions e = ep + es, where superscript s denotes secondary and p denotes primary.
Primary fields are defined as the fields in the earth if the earth were homogeneous. They satisfy the equations ahP amP VxeP=-po--po-, at at (3) and V X hP = a*ep where (J* is the conductivity of the half-space.
homogeneities in the half-space, satisfy The secondary fields, which are caused by in-
Solving for h" and setting V * h" = 0, assuming no magnetic permeability contrasts, we obtain an inhomogeneous vector diffusion equation for the secondary magnetic field, where the superscript s has been dropped for clarity. The field is a vector quantity whose components, in general, are coupled, which makes the 3-D equation much more complicated than the 2-D equation. Inside a homogeneous medium the vector components become decoupled, because the terms containing the gradient of conductivity vanish.

FINITE-DIFFERENCE APPROXIMATION
For time stepping, we use the Du Fort-Frankel finite-difference method. The method was derived by Du Fort and Frankel (1953) for solving a homogeneous diffusion equation in one dimension. It possesses the desirable properties of being unconditionally stable and explicit. It uses a staggered grid which, when combined with a leap-frog method in time, forms an efficient scheme for time-stepping computations.
The Du Fort-Frankel technique was used successfully for computing responses of 2-D conductors in geophysical environments by Oristaglio & Hohmann (1984), and later by Adhidjaja et al. (1985). Oristaglio & Hohmann made an Downloaded from https://academic.oup.com/gji/article-abstract/98/2/233/638652 by guest on 25 July 2018 Transient electromagnetic response 235 analysis of a stability condition for the method applied to 2-D problems when computed with uniform or non-uniform grid spacings. Here, we extend their 2-D finite-difference formulation to a 3-D formulation.
First, we consider the equation in the host where ( I = u*, so that u, = 0 and V ( l / u ) = 0. Hence, equation (7) reduces to a homogeneous diffusion equation, and in the Cartesian coordinate system each component of the magnetic field satisfies a scalar equation, To derive the finite-difference approximation, we extend the method used by Oristaglio & Hohmann (1984).
Integrating equation (8) over a volume, and applying the divergence theorem to the first term, we get For illustration we derive a finite-difference expression for one component of the vector field at node j , i, k and at time t = n + 1. Fig. 1 shows one of the typical finite-difference molecules, including the distribution of its physical properties, for the general case. In the case we are considering, all of the conductivities are the same, but we retain the general notation so that the results can be applied at boundaries later.
First, evaluate equation (9) over a small volume ABCDEFGH as shown in Fig. 1  showing odd (0) and even ( x ) nodes, conductivity distribution, and volume used in computing the field at the centre node. and using superscript n to denote the nth time we obtain 1 h"(j, i + 1, k )h"(j, i, k ) Applying the Du Fort-Frankel method by substituting Next, we consider the equation inside an inhomogeneity where u is constant, but u # u, and a, # 0. Thus equation (7) where *h"+'(j, i, k) is computed similarly to equation (10) with Finally, we consider the equation in a general region, specifically at interfaces. All three vector components are coupled, and satisfy equation (7). If we solved equation (7) directly, we would face a mathematical problem because we would like to compute results for models with large conductivity contrasts (1:lO to 1:lOOO). For a finitedifference grid with a block-wise conductivity variation the gradient of the conductivity becomes infinite at interfaces.
Brewitt-Taylor & Weaver (1976) suggested smoothing such conductivity variations by averaging, to avoid dealing with a discontinuous function. This approach is very appealing if it works, because it would be easy to compute a model with an arbitrary conductivity distribution.
Following Brewitt-Taylor & Weaver (1976), we solved equation (7) using a straightforward approach. The conductivities are averaged everywhere and assumed to vary smoothly, and terms involving derivatives of conductivities are approximated with central differences. With this procedure, a finite-difference expression for the x component of the magnetic field at node j , i, k, is where *h:+'(j, i, k) is computed similarly to equation (10) with a, = ( q i , k -u,)/p. Here dh$/dt is the x component of the primary magnetic field and e$' is the y component of the primary electric field computed analytically at node j , i, k and time t = n. The other variables are given by Finite-difference expressions for the other components are derived in a similar manner.
In our experiment, this scheme diverged for a model with the transmitting loop not centred over the body, which was 20 times as conductive as the surrounding earth. We believe the problem was caused by the inaccuracy in approximating the gradient of the conductivity. In contrast to frequencydomain computation, time-domain computation demands higher accuracy, because a small error in the solution at early times is propagated and compounded at later times.
Alternatively, we can avoid averaging the conductivity if we consider two media of different conductivities with a bounding surface, and perform the V( l / u ) operation, resulting in a 2-D delta function (6) normal to the interface.
Since the delta function is non-vanishing only at the interface, the corresponding volume integrations collapse into surface integrations, where superscript and subscript 1 and 2 indicate the conductivities of media 1 and 2, respectively, and n is the normal unit vector pointing from medium 1 to medium 2. However, this equation still poses a problem because the conductivity a, and the tangential component of current density, (V x h),, are discontinuous at the interface. To circumvent this problem we perform the surface integrations after taking the limiting values of those terms by approaching the interface from left (medium 1) and right (medium 2). Then, after performing the cross-products, where superscripts indicate the medium in which the terms are evaluated, and subscripts t indicate tangential components on the interface. In effect, we have applied the condition of continuity of h and e, at the interface. Terms on the right-hand side of equation (15) can be thought of as sources for the diffusion equation; they consist of a volume current (the volume integral) and anomalous surface currents (surface integrals). The sources involving primary fields can be computed analytically, while the sources involving secondary fields must be computed numerically.
As an example, consider an interface normal to the z axis. Then the equation for the x component, h,, is where j; and j; are the y components of current densities in medium 2 and 1, respectively, which can be computed either by directly performing the curl operator on h or by applying Stokes' theorem.
For illustration we derive a finite-difference expression for equation (16). Use of Stokes' theorem to compute the current density gives a finite-difference expression for h,, where **h:+'(j, i, k) is computed similarly to equation (12), but with u,,i,k as the volume averaged conductivity, and area, = (Axi + Axi+l)Azk+l, and e; is the y component of primary field computed analytically at node j , i, k and time t = n. Finite-difference expressions for the other components can be derived in a similar manner.
The disadvantage of this analysis is that the solution is not general, i.e. it applies to media involving only two different conductivities, and generalization of it would lead to a very j ,i-l,k '1 , i , k + l l i n e i n t e g r a l p a t h Figure 2. Finite-difference grid showing four conductivities intersecting at a corner, and the line-integration path to compute an average current density at the corner node. complicated numerical scheme. As an approximation, however, it may be generalized rather simply by replacing the anomalous current density ( j s -j i ) with an averaged current density.
For example, at a comer where four different conductivities intersect, the average current density j: can be computed by integrating h along the line integral path shown in Fig. 2. Here, in fact, we have done an averaging, but we average the current densities instead of the conductivities.

Time step and spatial step
The analysis given by Oristaglio & Hohmann (1984), modified for 3-D, suggests that the Du Fort-Frankel method retains its diffusive nature if the following condition is satisfied: where At is the time step, t is time range of interest and As is the spatial step. The At is also the practical time step for the method with a uniform grid and spatial step As.
Discussions on the condition and time step are given in the Appendix.
In our modelling, where a typical grid is non-uniform with variable a, we take both As and u as the smallest values in the model, and increase At gradually during computation.
As examples, for As = 10 m, cr = 0.01 S m-', and t = 0.1 vs, We use a discretization design similar to that of 2-D modelling: the grid is discretized finely and uniformly within the area of interest, and gradually is made coarse and non-uniform toward the boundaries.

BOUNDARY CONDITIONS
Inside the Earth, we use a homogeneous Dirichlet boundary condition, setting both the normal and tangential components of a vector field to zero. This type of boundary condition has been tested for a 2-D scalar field with good success (Oristaglio & Hohmann 1984). It appears that in a conductive host, the amplitude of a diffusing field is attenuated (damped) strongly enough so that with a reasonably distant grid boundary a reflected field will not interfere in the area of interest.
At the air-Earth interface, we use an upward continuation technique to set the boundary condition, assuming that the fields obey Laplace's equation in the air. We perform upward continuation on the known field (at a given time step) on the surface of the Earth one grid step in the air, and use it as a boundary value for the field in the future (field at the next time step).
The upward continuation equation is m m where h(5, v, 0) is the field at the Earth's surface, and h(x, y, z) is the field at one grid interval z above the surface. In 2-D modelling, the corresponding equation is a 1-D integral that can be computed efficiently by interpolating with a cubic spline function to get evenly spaced samples and then Fourier transforming with a 1-D fast Fourier transform (FIT).
In 3-D modelling we compute equation (19) directly with a trapezoidal rule. For a short continuation distance (small z ) the upward continuation operator has a narrow peak, so that it can be approximated with a small operator. The upward continuation is done point by point, which better suits our non-uniform grid. With this procedure, however, we need to interpolate the field to get the correct time level, and in addition to get evenly spaced fields near the point. For the interpolation, first we used a bicubic spline function. Because the bicubic spline function is very time-consuming to compute, we replaced it with a linear function based on the three surrounding points. However the interpolation is still time-consuming, due to the pointwise upward continuation.
Our test of the upward continuation indicates that for a distance of about 10m, a 21 x 21 operator gives about 10 per cent error when applied to a smooth function. A shorter operator, such as 11 x 11, gives about 20 per cent error. Because the longer operator is very time-consuming to compute, we used a 5 x 5 operator and the bicubic spline function for interpolation for the earlier models, and used an 11 x 11 operator and the linear function for the model in this paper.

NUMERICAL RESULTS FOR STEP RESPONSE
To compute the step response we use fields due t o a step current in the transmitter for the ep and dbp/dr terms in equation (7). As a test we computed results for the coaxial model shown in Fig. 3. It is a 100 x 100 X 100 m body with resistivity 1 a m m -' embedded 80m deep in a 100 Om-' half-space. The transmitter is a 100 X 100 m loop carrying 1 A of current and centred over the body. This model was designed to test the case where the response is dominated by vortex current effects, but, as it turns out, it becomes an important model to test the galvanic current effect also. For this model the terms in equation (7) which contribute to the galvanic-current response would be small, but must be calculated accurately otherwise they may change the character of the response, such as causing a sign change or a wrong decay rate.
We discretized the model with the smallest cells being 10 x 10 x 10 m and used a uniform grid within the region of interest and a non-uniform grid in the region far from the body. The time step ranged from 0.1 to 20 p . . As in the 2-D case, the time step was increased gradually during the computation. In the 2-D case the time step typically ranges from 0.1 to 40 p. 50000Hz, and were sampled at five points/decade. The body was discretized with cell size 12.5 X 12.5 X 12.5 m, or 128 cells/quadrant. Figure 3 compares decay curves of the horizontal component, b,, computed with the finite-difference and the integral-equation methods. The responses were computed for a point P 100m from the centre of the loop and lying outside the loop. The finite-difference solution is slightly higher than the integral-equation solution at early times, but decays slightly slower at later times. It also shows a sign change at about 2 0 p , where the solution of the integral-equation method is not available. However, the overall agreement is reasonable. However, at times later than 0.15 there is a large discrepancy in the results. A sudden change in decay rate suggests that the finite-difference solution breaks down at this point.

Comparison with integral-equation solution
The magnetic field impulse response (emf) of a spheroidal body in a conducting whole-space decays as a power law, t-3.5 , at late times (Kaufman 1981). Accordingly, the step response of the body should decay as t-2.5 at late times. Our results do not seem to have the correct decay rate, but the responses plotted may not be late enough for determining the decay rate.   to both the primary field and the current on top of the conductor, the current density distributions are as expected, flowing in a clockwise direction. The early-time current density distribution on the plane may be real, because its amplitude is much smaller than that of the primary field and Because of our limited resources we did not compute further results. The computation of the solution to 1.5ms took 5 hr of cpu time on a Cray X-MPf48 supercomputer.

Magnetic field and current density plots
We examine the qualitative behaviour of the finitedifference solution by plotting the magnetic field and current density distributions in the Earth. Figure 5 shows snapshots of secondary magnetic field patterns in a vertical cross-section through the centre of the body. The patterns show a magnetic dipole field which starts developing on top of the conductor and grows stronger with time. The dipole has a positive polarity as expected (z-axis positive downward), the same direction as the primary field in the transmitter loop shown in Fig. 3. When the current in the transmitter is shut off, the current in the Earth flows in the same direction as the current in the transmitter did, so as to preserve the original magnetic field. At 1.5 ms the curvature of the magnetic fields nearly disappears, showing almost a uniform field inside the body, which is characteristic of the late-time response. However 1.5 ms is an early time for this model. It seems that the solution breaks down at this time as evidenced also from the current-density plots discussed later. Figure 6 shows snapshots of current density distributions associated with the dipole, in a horizontal section through the centre of the body. Except at very early time when the current density on this plane flows in the opposite direction also the distribution at 0.035 ms indicates a transition period when the current slowly switches direction. Opposite flows of current on the two planes suggest that more complicated dipoles are induced in the conductor at early times. Figure 7 shows snapshots of current density distributions associated with the total electric field in a horizontal section through the centre of the body. The current flows in the same direction as the electric field in the earth when the current in the transmitter is shut off. However, at 1.5 ms the ratio of tangential current density outside and inside the body is far from the ideal 1 : 100.
Nevertheless, the general pattern of both magnetic field and current density are reasonable.

DISCUSSION
Initially, we computed the magnetic field impulse response by substituting electric and magnetic fields due to an impulse current for the ep and dhP/dt terms in equation (7). We used finite-difference schemes both with conductivity averaging, and with current-density averaging, and computed results for two models: a coaxial model which is designed to test the solution for a vortex-dominated case, and a model where the body is located with its strike parallel to the long side of the loop, designed to test the galvanic-dominated solution. However, the impulse responses computed with both finite-difference schemes were not accurate. They showed amplitudes that were too large compared with those of the integral-equation solution (San Filipo & Hohmann 1985), and showed an unexpected sign change indicating that the solution oscillates in time.
The problem may be caused by the rapidly varying primary fields in the earth, which are the source in the diffusion equation. Unlike that for the 2-D case, the diffusion equation for the 3-D model has three source terms that interact. For a single source such as in 2-D modelling, the finite-difference method seems to propagate the field correctly, but it has a problem when more and interacting sources are involved. Alford, Kelly & Boore (1974) presented an analysis for accuracy of a finite-difference solution for the acoustic wave equation with regard to the frequency content of a source function. They discussed the necessary number of samples per wavelength in order to avoid dispersion from undersampling, which is known as "grid dispersion".
Another possible reason for the inaccuracy problem, which is more serious, is the mixture of numerical and analytical computation for the source terms. Accuracy of the central differences in approximating the curl operation of the line integral for Stokes' theorem is limited by the grid interval, while the sources with analytical expressions can be computed very accurately. If this is the case, then we face a fundamental problem with our numerical scheme: to achieve accuracy we may need to refine the discretization beyond what is practical. Therefore, we have to seek an alternative way to solve the problem.
One alternative is to solve for the total field, whose equation is simpler and without the analytical sources, but, as we learned from 2-D modelling, the solution requires a larger grid and significantly more computation time. The alternative that we pursued was to compute the step response, i.e. the magnetic field due to a step current. Computation of the step response involves smoother functions; also the analytical source terms have fewer sign reversals, which may alleviate the problem of inaccuracy due to coarse discretization.
Our experiments show that, in general, the shapes of profiles or decay curves of the step response are simpler than those of the impulse response. More importantly, the solution shows no oscillation of sign as in the case of the impulse response, confirming that the step response is better behaved. However, there is still a large discrepancy in amplitude for the vertical field component of the finite-difference solution and the integral-equation solution in the last experiment. Possibly the computer codes may contain errors or may suffer from numerical inaccuracy because the source functions are still too sharp. Smoothness of the source functions seems to be an important factor in our 3-D solution. For the next experiment we should concentrate on analyzing the source function, and try a smoother source waveform, such as a step with ramp turn-off.
In addition, we face an efficiency problem. At present the 3-D code is very time-consuming to run even for a small model. However, symmetry may be applied to some special-case models which could be useful for testing purposes.
The computation for the boundary condition at the air-Earth interface is one of the most time-consuming parts of the solution, because it involves interpolation and integration that must be done at each time step. Use of a uniform grid, although it requires very large memory, would reduce the computation times because it would simplify the interpolations, or it would permit a direct use of a highly optimized 2-D FFT routine to perform the convolution in the transform domain. In addition, it often makes the computer codes amenable for vectorization to take full advantage of the supercomputer features. Utilization of a non-uniform grid makes the computer codes more complicated, which often prohibits vectorization.

Transient electromagnetic response 241
Another time-consuming computation that must be done at each time step is the computation of the analytical sources of the diffusion equation, especially for a source other than an impulse current. The step response takes longer t o compute than the impulse response because of extra numerical integrations required for computing the analytic sources. This problem might be reduced by tabulating the source fields and using a simple interpolation to compute the desired field. Tabulation of 3-D space-time functions, however, would require a large amount of memory, such as that of a Cray-2.
The results presented here were computed on the Cray X-MP/48 of the San Diego Supercomputer Center (SDSC). Unfortunately, the maximum memory available is only three million words; running our program with a large amount of memory is either very costly on a high priority or extremely slow on a low priority. For all models computed the maximum size of the grid used was 50X5Ox50, which required memory of just under one million words.
With the present version of our 3-D code, to compute the solution for the small model shown in Fig. 3 to 10 ms would take 10h (estimated) of cpu time on the Cray X-MP/48 supercomputer.
Faster computation might be achieved by changing the grid several times, as is commonly done in frequencydomain computationsfine discretization for high frequencies and coarse for low frequencies. However, in the time domain, this procedure would require an accurate 3-D interpolation when changing the grid size.
W e finally conclude that the development of the 3-D finite-difference TEM program is not practical given our available resources. However, the result of the stepresponse test, although not completely satisfactory, is somewhat encouraging. In the future, when adequate computational resources are more readily available, the scheme may be worth pursuing.

Stability condition and time step for the Du Fort-Frankel finite-difference method
The Du Fort-Frankel finite-difference method was derived for the solution of a homogeneous diffusion (parabolic) equation in 1-D by Du Fort & Frankel (1953). They also showed that the method is unconditionally stable. Stability proofs of the method for higher dimensions are given by Birtwistle (1968). Analysis of stability criteria for the leapfrog Du Fort-Frankel method applied to an advectivediffusive equation is given by Cushman-Roisin (1984). All of the analyses were done for a homogeneous diffusion equation discretized with a uniform grid.
Oristaglio & Hohmann (1984) gave analysis of stability in 2-D as applied to the geophysical problem, including generalization of it for a non-uniform grid. The Du Fort-Frankel method is unconditionally stable, but it is still limited by the consistency condition as the mesh spacing is refined (Lapidus & Pinder 1982). As shown by Du Fort & Frankel (1953) the method approximates a hyperbolic equation instead of a diffusion equation if A x and At approach zero with constant ratio.
Following Oristaglio & Hohmann (1984), we will show that the condition required in the Du Fort-Frankel method for our purpose can be satisfied. This analysis is an approximation based on scalar wave and diffusion equations in a whole-space, which are approximately applicable to our equation in the larger part of a typical grid, in the Cartesian coordinate.