A Comparison Study of Explicit and Implicit 3-D Transient Electromagnetic Forward Modeling Schemes on Multi-Resolution Grid

: This study compares the efﬁciency of 3-D transient electromagnetic forward modeling schemes on the multi-resolution grid for various modeling scenarios. We developed time-domain ﬁnite-difference modeling based on the explicit scheme earlier. In this work, we additionally implement 3-D transient electromagnetic forward modeling using the backward Euler implicit scheme. The iterative solver is used for solving the system of equations and requires a proper initial guess that has signiﬁcant effect on the convergence. The standard approach usually employs the solution of a previous time step as an initial guess, which might be too conservative. Instead, we test various initial guesses based on the linear extrapolation or linear combination of the solutions from several previous steps. We build up the implicit scheme forward modeling on the multi-resolution grid, which allows for the adjustment of the horizontal resolution with depth, hence improving the performance of the forward operator. Synthetic examples show the implicit scheme forward modeling using the linearly combined initial guess estimate on the multi-resolution grid additionally reduces the run time compared to the standard initial guess approach. The result of comparison between the implicit scheme developed here with the previously developed explicit scheme shows that the explicit scheme modeling is more efﬁcient for more conductive background models often found in environmental studies. However, the implicit scheme modeling is more suitable for the simulation with highly resistive background models, usually occurring in mineral exploration scenarios. Thus, the inverse problem can be solved using more efﬁcient forward solution depending on the modeling setup and background resistivity.


Introduction
The transient electromagnetic (TEM) method is widely used in many geophysical applications, including exploration for minerals [1][2][3]. In order to fully exploit the observed data, the full 3-D quantitative data inversion is warranted. The inversion techniques rely heavily on accurate and efficient forward modeling algorithms.
The TEM forward modeling can be implemented in frequency-domain and thereafter transformed to time-domain using the Fourier transform. However, the frequency-domain approach generally requires a sufficiently dense and wide frequency range to ensure the accuracy of the solution [4]. Alternatively, the simulation can be implemented directly in the time-domain so that the electromagnetic (EM) fields are recursively advanced using time-stepping [5]. The time-domain TEM modeling can be further classified into explicit and implicit schemes, depending on how the time-stepping is approached.
The classical explicit scheme discretizes Maxwell equations directly in a leapfrog fashion. It involves lightweight computations in each time step, alternating between electric and magnetic fields from Ampere's and Faraday's laws, which are merely matrixvector operations. However, the time-stepping is conditionally stable in that the length of the time step is constrained by the Courant-Friedrichs-Lewy (CFL) condition [6], i.e., the time that the EM wave spends to pass the minimum length of the model grid at light speed. This usually small time step needs to be fixed during the whole modeling process, which is impractical to implement for the time range of the TEM problem. Oristaglio and Hohmann [7] and Wang and Hohmann [8] applied a modified version of the Du Fort-Frankel method, which replaces the displacement currents term with a variable artificial term that can manually slow down the EM wave speed in the late time. Thus, the time step is allowed to increase gradually with time, whilst the diffusion property of the EM field propagation is maintained. The stable conditions are determined by the minimum grid size and conductivity of the model.
The implicit scheme is usually implemented using the backward Euler method to discretize the time-stepping. This approach is unconditionally stable, allowing arbitrary time steps [9,10]. However, each time step requires the solution of a system of linear equations that is often solved using the direct solver (e.g., Cholesky decomposition method) or the iterative solver (e.g., Krylov methods). The direct solver is expensive at the matrix factorization stage. However, it is only done once if the model and time step are fixed. The advantage is that the subsequent computations are simple and fast, and the solution is accurate [10,11]. The saved factorization can also be utilized when solving the adjoint problem in the sensitivities computations.
The iterative solver requires relatively lower computational resources. Iterations start from an initial guess to approximate the solution. A standard initial guess approach uses the solution from the previous time step [12]. Domnikov [13] proposed starting from a linear combination of two previous time steps. The coefficients of the linear combination are determined by minimizing a residual function. A better initial guess together with a proper preconditioner significantly improves the solving process. The time-stepping generally begins from a short time step to ensure the solution's accuracy at the early time. As time progresses, the variations of the EM fields are smoother, and the larger time step can be used to reduce the number of steps. However, when the coefficient matrix is altered using a different time step, additional matrix preconditioning or factorization is required. Um et al. [12] proposed an Adaptive Time Step Doubling (ATSD) scheme, which is a rather efficient strategy for enlarging time steps with fewer times of changing the coefficient matrix.
In order to model EM fields accurately around the source location, fine grid discretization is usually required. However, EM fields propagate into the subsurface in a diffusive manner, and the field variations become smoother with depth. Therefore, the horizontal grid resolution can be coarser with depth. Several approaches with variable model discretization were proposed, e.g., Octree grid [14,15], unstructured grid discretized by the tetrahedral cells [4,10] for the finite-element method modeling, etc. Our previous study developed the Multi-Resolution (MR) gird approach, which has been initially implemented in the 3-D magnetotelluric forward modeling [16] and extended to the direct current resistivity method [17] and TEM method using the explicit scheme [18]. The MR grid can be thought of as a subclass of Octree grids adequately discretizing the media.
In this study, we focus on the implicit scheme TEM forward modeling implemented on the MR grid. The system matrix is symmetric and positive definite. Therefore, a preconditioned conjugate gradient method together with incomplete Cholesky factorization as the preconditioner are employed to solve the governing system of linear equations. A proper initial guess affects the convergence performance of iterative solvers. Hence, we test various initial guess estimates. One can take two previous time steps to estimate the initial guess at current time step in a linear extrapolation manner. Based on the linearly combined approach, we also evaluate the effects of involving more previous steps. Those initial guess estimates are compared to demonstrate the most efficient one. Finally, a comparison between the explicit and implicit schemes modeling on various resistivity scenarios is present to show the preferable condition for them.
The paper is organized as follows. We first briefly introduce the explicit scheme modeling. Then we describe the implicit scheme solution in more detail and explore the effects of the initial guess. The MR grid approach is also shortly explained. Three synthetic examples are presented to verify the accuracy and efficiency of 3-D TEM forward modeling using the implicit scheme on the MR grid. Finally, we compare the explicit and backward Euler implicit schemes for various modeling scenarios.

Methods
The electromagnetic field variations in the time-domain are governed by Maxwell equations: where e and h denote the electric and magnetic fields, respectively; σ represents the spatially distributed electric conductivity; ε 0 and µ 0 are the constant dielectric permittivity and magnetic permeability of the vacuum, respectively; j src indicates the source electric current density. The displacement currents term ε 0 ∂e ∂t for the time ranges considered in TEM modeling can often be neglected. By ignoring ε 0 ∂e ∂t and substituting Ampere's law into Faraday's law, we also obtain the second-order diffusion equation for electric field e: Equation (1) or (2) should be discretized both in spatial and time domains. In this study, we use the finite-difference method to discretize the spatial coordinates. Thus, e ∈ R N E and h ∈ R N F are respectively defined on edges and faces of the staggered (SG) grid [19], where N E and N F denote the amounts of edges and faces in the grid. In order to satisfy Ohm's law, the diagonal matrix σ ∈ R N E ×N E of the averaged conductivity on grid edges is also calculated from the cell conductivity σ through the volume-weighted averaging approach. The curl operators are discretized as C ∈ R N F ×N E , which maps from the grid edges to faces, and its adjoint C † ∈ R N E ×N F works in reverse. C and C † are defined as: where C ∈ R N F ×N E and its transpose are expressed respectively as the topology matrices of C and C † . The diagonal matrices L E (edge-lengths),Ã F (dual-face-areas) ∈ R N E ×N E , and A F (face-areas),L E (dual-edge-lengths) ∈ R N F ×N F are the additional Cartesian metric elements defined on the grid. TEM forward modeling starts time-stepping from the initial statues of fields. Here, we study the ground survey ( Figure 1) using square loop source with a step-off excitation, but the loop source can also be any arbitrary geometry and waveform excitation. The −db z /dt response is computed for each site, where b = µ 0 h denotes the magnetic flux density. The initial h, generated by the source current (j src ), is static before the switch-off time and can be calculated via Biot-Savart law [18]. Correspondingly, the initial e field is zero over the modeling domain Ω.

Explicit Scheme in the Modified Version of the Du Fort-Frankel Method
We implemented the explicit scheme using the modified version of the Du Fort-Frankel method, which can be expressed by discretizing and reorganizing Equation (1) in an iterative form: where e n and h n are defined at the time instants t n andt n , respectively. The iterations between e n and h n are progressing in a leapfrog fashion, ∆t n and ∆t n denote the lengths of time steps, respectively. The time axis and iterations are depicted in Figure 2. The displacement currents term often has negligible effects. Therefore, the original ε 0 is replaced with a variable artificial factor γ n to enable increasing ∆t n with time and maintain the stability of the forward modeling [8]. The stability criterion is expressed as: Time step ∆t n is allowed to increase gradually. However, it is determined by the minimum conductivity and minimum edge-length of the grid. Thus, if a model is relatively conductive, ∆t n will be relatively large, and the explicit scheme's performance will be superior (we show comparison thereafter).
The Dirichlet boundary conditions are used for the surrounding and bottom distant boundaries Γ ∞ , where the tangential e n are assumed to be zero, i.e., (n × e) Γ ∞ = 0 [5]. Since involving a low air conductivity (σ air ≈ 0) would significantly suppress ∆t n in Equation (5), we use the upward continuation boundary condition [8,20] for the top boundary. Thus, the Earth's surface is assumed to be flat, and the model requires only one air layer. The horizontal components h x and h y in the air layer (z < 0) are derived from the vertical com-ponent h z on the ground (z = 0) through the Fast Fourier transform (FFT) [18]. Therefore, ∆t n is determined only by the subsurface conductivity.
Because of the rich null space of curl operators, the numerical errors can accumulate during the time-stepping and cause the late time response to be inaccurate. The divergence correction is a standard way to cope with this problem. Wang and Hohmann [8] proposed a scheme to explicitly involve the divergence-free condition of h in the forward modeling (i.e., ∇ · h = 0). At each time instantt n , the vertical component h z is computed from the horizontal h x and h y in the subsurface layers [18], ensuring the divergence-free of the magnetic field.
The explicit scheme is initiated from the static fields e 0 and h 0 and runs through the recursion of Equation (4). Each time step only consists of lightweight matrix computations. When the time-stepping is passing a time gate, the −db z /dt data response is computed

Implicit Scheme Using the Backward Euler Method
The implicit scheme is based on solving Equation (2). In order to improve the accuracy, the second-order backward Euler method is used to discretize time derivatives [4], e.g., ∂e ∂t is approximated as: Consequently, Equation (2) is discretized to an iterative form: where A ∈ R N E ×N E represents the coefficient matrix of the system of linear equations determined by ∆t. The right-hand-side term b n ∈ R N E consists of source term and the field solutions from previous time steps. Equations are only solved for e n during the timestepping ( Figure 2). The nonsymmetric curl-curl A can be symmetrized by multiplying with the diagonal matrix of edge-volumes V E = L EÃF , ∈ R N E ×N E : whereb n denotes the new weighted right-hand-side term, accordingly. The above system of equations is usually large for the 3-D problem and can be solved using Krylov subspace iterative methods. In our implementation, we use the Preconditioned Conjugate Gradient (PCG) method, which is an efficient solver for the equations with a symmetric coefficient matrix. The Incomplete Cholesky factorization (ICHOL) method is used as the preconditioner. An initial guess x 0 of e n solution affects the performance of the iterative solvers. Usually, the standard initial guess approach adopts the solution from the previous time step to start PCG iterations [12], i.e., x 0 = e n−1 . However, we have found that this approach might be too conservative. Instead, we first propose a linearly extrapolated approach that calculates a linear extrapolation of x 0 under the simple assumption that there is a gradual change between adjacent time steps, i.e., e n − e n−1 ≈ e n−1 − e n−2 . Therefore, x 0 = 2e n−1 − e n−2 is taken as the new initial guess vector of e n . Domnikov [13] proposed a linearly combined approach that approximates e n as x 0 = α 1 e n−1 + α 2 e n−2 . Unlike the linearly extrapolated approach, the coefficients α = ( α 1 α 2 ) are determined from solving min Ã e prev α −b n 2 at each time step where e prev = (e n−1 , e n−2 ), and α is obtained using the least-squares solution as α = e T prevÃ TÃ e prev −1 e T prevÃ Tb n . We also consider a linear combination of more previous steps (N prev = 2, 3, 4) as The time-stepping usually starts with a sufficiently small ∆t to ensure accuracy at early times. Using a constant ∆t with fixed model parameters does not require modifying A and recomputing matrix preconditioning (or full matrix factorization when using a direct solver). We only need to updateb n in order to advance the time-stepping. However, the small ∆t is not needed at late times; therefore, it is more efficient to increase ∆t gradually. However, this should be compromised with the number of updates ofÃ. We follow the Adaptive Time Step Doubling (ATSD) scheme proposed by Um et al. [12]. The forward modeling starts with a suitable ∆t at the beginning, e.g., 1 × 10 −7 s, and progresses a certain number of time steps N ∆t . Afterwards, we increase ∆t and verify the solution's accuracy. If it satisfies the predefined tolerance, we accept the new ∆t. This procedure is repeated until the end, as shown in Figure 2.
In contrast to the explicit scheme, the Dirichlet boundary conditions, (n × e) Γ ∞ = 0, are used for all boundaries, including the top face. Therefore, the model contains several air layers to move the top boundary away from the study region. Even though air layers increase the computations, they also allow for the simple implementation of topography into the forward modeling, which is an advantage compared to the upward continuation boundary condition.
The implicit scheme suffers from the same numerical errors during the time-stepping as the explicit scheme. Accumulated errors may cause slow convergence at late time, and this also leads to nondivergent solutions [21]. Therefore, we use an additional step of divergence correction [22] to maintain the divergence-free condition of the current density. This correction significantly accelerates the convergence rate of the iterative solution. When solving Equation (8) in each time step, we compute the correction δe by solving a Poisson type equation: where K ∈ R N N ×N N indicates the coefficient matrix of the above Poisson equations set; N N represents the number of grid nodes; R ∈ R N N denotes the corresponding right-hand-side term; the unknown set φ ∈ R N N represents the potential defined on the grid nodes, and δe = ∇φ; G † ∈ R N N ×N E and G ∈ R N E ×N N are respectively the discrete divergence and gradient operators defined as: where G ∈ R N E ×N N denotes the topology matrix of G, and its transpose also works as the topology matrix of G † ; the diagonal matrixṼ C ∈ R N N ×N N represents dual-cell-volumes (volumes around nodes). Note thatṼ C is multiplied on both sides of Equation (9), and it also makes K symmetric. In this case, the combination of ICHOL and PCG is also used to solve the above equations set. Solving Equation (9) is similar to the direct current resistivity forward modeling problem. More details of implementations can be found in Cherevatova et al. [16] and Gao et al. [17]. Since K is independent of ∆t, the preconditioning for Equation (9) only needs to be done once. Similar to the explicit scheme, we start from the initial condition of e 0 = 0 and run the recursion of solving Equation (8) alternating with solving Equation (9). Each time step involves significantly more computations than the explicit scheme, but it may require fewer time steps. Because only e is estimated during the time-stepping, Faraday's law ∇ × e = −µ 0 ∂h/∂t in Equation (1) is used to calculate the −db z /dt response.

Multi-Resolution Grid
In the previous sections, we described TEM forward modeling on the SG grid. However, the multi-resolution (MR) grid (MR3DMod framework) has been shown to be a better-suited choice for such type of problem [16,18].
MR grid consists of several vertically assembled SG sub-grids ( Figure 1). Each subgrid is a regular SG grid but has a different horizontal resolution N i x × N i y , where i denotes the sub-grid index. MR grid can be easily created from a fundamental SG grid with the fine grid resolution N fine The horizontal resolution of a sub-grid is controlled by a coarseness parameter Cs i and determined as: The most important part of the MR implementation is the definitions of differential operators C, C † , G, and G † on the interfaces between sub-grids. In the interiors of the sub-grids, the structures of those operators are the same as the inside of the SG grid. Following the development of Cherevatova et al. [16], we adopt the Coarse Active scheme and reconstruct the topology matrices C and G on the common interfaces. Therefore, those operators can also be used in the MR grid. More details of the implementations can be found in Cherevatova et al. [16], Gao et al. [17,18].

Results
In this section, we present three synthetic examples to verify the accuracy and efficiency of TEM forward modeling using the implicit scheme on the MR grid. The first and second examples are used to validate the modeling in 1-D and 3-D, respectively. The comparison of SG and MR grids and comparison of using different initial guess (x 0 ) approaches are also shown. In the last example, we compare explicit and implicit schemes to show which one is more suitable under certain conditions. In the examples, the air resistivity is 10 6 Ωm. The ∆t selection strategy in the implicit scheme is as follows: ∆t starts with 10 −7 s, and it is enlarged after every N ∆t = 50 steps with a factor of 5. The tolerance of relative residual in PCG is specified as 10 −6 .

Example 1
We start with a simple 1-D model for which the semi-analytical solution is available [23]. The model consists of four layers with resistivity and thickness from top to bottom as (ρ, ∆z) = (100 Ωm, 80 m; 1000 Ωm, 60 m; 5 Ωm, 60 m; 100 Ωm, ∼). We place a 200 m × 200 m loop source on the ground with a step-off current waveform. A site located at (x, y) = (5 m, 5 m) measures −db z /dt response at 30 time gates covering 10 −5 s, 10 −2 s in a logarithmic interval. The model and survey setup are shown in Figure 3. The SG grid is designed with a resolution of 80 × 80 × 55, where the center of the model is discretized with 10 m uniform cells. The grid has 15 air layers start-ing with a 10 m thickness at the surface and increasing upwards. Furthermore, we create an MR grid from the SG grid with five sub-grids having coarseness parameters as Cs i , N i z ; | i=1,5 = [2, 5; 1, 5; 0, 22; 1, 14 ; 2, 9]. The finest horizontal resolution (80 × 80) is only around the air-earth interface (z ∈ [−310 m, 170 m]) and gradually roughening in the Earth with depth. Horizontal resolution in the air layers is also gradually decreased. The common interfaces between the sub-grids are located at z = −10, 230 m, −310 m, 170 m, and 558 m (Figure 3). The forward modeling responses are shown in Figure 4a, and the anomalies caused by the resistivity changing are clear to see. We compare solutions based on SG and MR grids against the 1-D solution. Our 3-D solutions agree well with the 1-D solution, and relative differences are shown in Figure 4b. The maximum relative differences are around 3.4% at the t = 2.21 × 10 −4 s time gate. The differences might be caused by the high resistivity contrast between layers in the 3-D modeling. The SG and MR grid approaches have similar accuracy in the early time. In the late time, the relative difference of the MR solution is slightly higher than the SG solution. This result is expected, as the MR grid has coarser horizontal resolutions at depth, but the accuracy reduction is negligible.  Next, we compare the efficiency of iterative solver by showing the computational time at each time step t step in Figure 5a, the total run time t total in Figure 5b, and the number of PCG iterations of each time step N iter in Figure 5c. The SG and MR grid cases are presented separately on the left and right columns in Figure 5. The time-stepping involves 270 steps in total, and ∆t has been enlarged five times. The general trend is that, once ∆t is enlarged, t step is also increased, which means using a larger ∆t requires more iterations to solve the system of equations. Subsequently, t step is gradually decreasing in the following time steps. This phenomenon can be interpreted as the drastic variations of e n being gradually attenuated as time passes and the equations becoming easier to solve. The performances of using different initial guesses (x 0 ) to initialize the iterative solver are also shown. The standard approach (x 0 = e n−1 ), the linearly extrapolated approach (x 0 = 2e n−1 − e n−2 ), and the linearly combined approach (x 0 = ∑ N prev i=1 α i e n−i ) arrive in the same modeling accuracy (the maximum relative differences between them are less than 0.005%). However, the latter two approaches require much less time for solving the equations. The linearly combined approach with N prev = 3 works more efficiently that requires much less N iter , and thereby significantly reduces t step and t total compared to the standard approach, e.g., t total = 613.7 s → 201.8 s for the SG grid and t total = 328.6 s → 84.8 s for the MR gird. For the linearly combined approach, the N prev = 4 case was expected to be superior to the N prev = 3 case because the former case involves one more previous time step and can be more flexible. However, the result is on the contrary. One can see that for t step and N iter of the N prev = 4 case, there exists many more perturbations in Figure 5. We may interpret this phenomenon as too much early time information dragging down the convergence on the new time step.
In addition, we compare the computational time of SG and MR grids using the linearly combined x 0 approach (N prev = 3 case). The numbers of Degrees of Freedom (DoF) to solve the main system of equations DoF e and the divergence correction equations DoF φ for the SG gird are 1,025,815 and 337,014 in contrast to 506,615 and 165,014 on the MR grid (i.e., 49.39% and 48.96% of the SG grid).
As shown in Figure 5, the MR grid is undoubtedly more computationally efficient than the SG grid approach. The ratios of the averaged t step (t step ) and t total between the MR and SG modelings are 0.31 s/0.75 s and 84.8 s/201.8 s, respectively. As a result, the MR grid requires only 42% computational time of the SG grid.

Example 2
In the second numerical test, we consider a 3-D example. The model has a 1-D background consisting of 3-layers: The model grid includes 15 air layers with the same discretization as in Example 1. MR grid is designed from the SG grid and has 5 sub-grids with the coarseness parameters: Cs i , N i z ; | i=1,5 = [2, 5; 1, 5; 0, 19; 1, 20 ; 2, 11]. Note that the third common interface is located at the place between two background layers, and it also crosses the deeper conductive block. The MR grid discretization is shown in Figure 6. Since the differences between the SG and MR grid modeling responses are hardly visible, only the MR grid solution is shown in Figure 7a. The anomalous responses caused by the conductive blocks are clearly visible. The 3-D solution using the Lanczos spectral decomposition method (SLDM) [24] was chosen as the reference to verify the accuracy of our solutions. The relative differences between SG and MR solutions against the SLDM solution are shown in Figures 7b,c, respectively. The maximum relative difference of the SG and MR grid solutions are 3.44% at 5.29 × 10 −5 s in Figure 7(b1) and 3.59% at 2.55 × 10 −4 s in Figure 7(c2), respectively.
The same efficiency parameters t step , t total , and N iter in Example 1 were also evaluated here (Figure 8). Similarly, the iterative solver using either the linearly extrapolated x 0 approach or the linearly combined x 0 approach results in much shorter run time than the standard x 0 approach. The linearly combined approach with N prev = 3 works best again, e.g., for the SG grid case, t total : 213.17 s versus 759.32 s, andt step : 0.79 s versus 2.81 s in contrast to the standard x 0 approach.
Benefiting from using the MR grid that involves much less DoF e = 480, 060 and DoF φ = 156, 139 in contrast to the SG grid (DoF e = 1, 120, 220 and DoF φ = 368, 219), we can solve a much smaller forward modeling problem and consequently further reduce the modeling time:t step = 0.79 s → 0.4 s, t total = 213.17 s → 107.67 s.

Example 3
In the last example, we compare explicit and implicit schemes for different background model resistivity. Three models are tested that have a conductive block located in the same position:  Figure 9.
In previous examples, we have found that the implicit scheme modeling using the linearly combined x 0 approach with N prev = 3 has much better performance than the other initial guess approaches. Therefore, only the former approach is presented in this example. The modeling responses are shown in Figure 10a. As expected, the amplitude of the response of the conductive model (e.g., 100 Ωm) is generally higher than the resistive one (e.g., 10,000 Ωm). The anomaly response in the conductive model is observed at a later time, whilst anomaly responses are visible at an earlier time for the resistive models. This phenomenon is due to the different skin depths of EM fields in those three models. The modeling responses of explicit and implicit schemes are compared, and their relative differences are shown in Figure 10b. For the models with 100 Ωm and 1000 Ωm background resistivity, the relative differences between explicit and implicit schemes are small, which are 2.34% and 3.07% respectively. For the 10,000 Ωm model, the relative difference is slightly higher (5.05%). The efficiencies of the forward modeling using explicit and implicit schemes are compared and shown in Table 1. The explicit scheme usually requires more time steps (N step ) than the implicit scheme, but each time step is faster.
It is clearly visible for the conductive (100 Ωm) model that the explicit scheme took only 97.61 s to finish the job, which is more efficient than the implicit scheme (265.72 s). However, the explicit scheme is strongly affected by the background resistivity. As the resistivity increases, N step is increased accordingly (3089 → 9724 → 30, 296), and the modeling speed is degrading (97.61 s → 1205.11 s → 14, 244.9 s). The reason is that the EM fields propagate faster in the resistive media, and ∆t n must be suppressed to stabilize the modeling. In addition, the upward continuation boundary condition used in the explicit scheme requires the inclusion of a large area for 2-D FFT. Since FFT requires an equidistant grid that is usually discretized by the minimum horizontal edge-length to match the short wave variations of the fields, the larger area, especially at the padding cells, significantly increases the computations [5]. The same influence ont step can also be found (0.032 s → 0.12 s → 0.47 s), which is also mainly caused by the heavier FFT computations in the resistive models. In contrast, the implicit scheme is less affected by the model's resistivity. Even though t total is also gradually increasing for the resistive models, which means the equations took a longer time (t step ) to be solved, the run time fluctuations are relatively stable.
Thus, the explicit scheme modeling is more efficient for the model with a conductive background, whilst the implicit scheme modeling is preferable for the resistive one.

Discussion
We have implemented the 3-D TEM implicit scheme forward modeling based on the second-order backward Euler method on the multi-resolution (MR) grid. The implicit scheme forward modeling requires the solution of the system of linear equations, which is often done using iterative solvers. We have found that initializing the iterative solver with a better guess has a profound influence on iterations. The standard initial guess approach usually uses the previous time step's solution, which is too conservative. We have tested the linearly extrapolated and linearly combined approaches based on the solutions from several previous time steps to derive a better initial guess. We show that it improves the convergence rate significantly and therefore reduces the modeling time. Our result also shows the linearly combined initial guess using three previous time steps works the best. Involving more previous time steps to predict the current solution does not improve the performance further.
The implicit scheme modeling is tested using three synthetic examples against other algorithms to verify the accuracy and efficiency. The MR grid substantially decreases the number of degrees of freedom and results in better computational efficiency without compromising the accuracy.
Finally, we compare 3-D TEM forward modeling using the explicit and implicit schemes on different models and derive strategies to choose the more efficient solution. The explicit scheme generally requires more time steps than the implicit scheme, but the computational costs at each time step are lower. The result shows that the explicit scheme is more suitable for the simulation in a conductive environment, e.g., environmental studies or sedimentary basins.
In contrast, the modeling speed of the implicit scheme is less affected by the background resistivity. Hence, the implicit scheme is preferable in a resistive environment, e.g., shield areas. Future studies will implement the forward modeling into the inversion and use the real field data to examine the algorithms. We would use an initial run to verify which modeling scheme would be more efficient for certain prior model conditions. This testing can allow us to choose the modeling scheme and optimize the inversion algorithm with respect to speed and accuracy. In the previous work, we also used time-domain electromagnetic modeling to simulate electric fields at the surface for the geomagnetically induced currents simulation. Based on the results presented in this study, we expect that the implicit scheme modeling would be more suitable for this purpose. Funding: This research was supported by the ARN project, and the project is funded by the European Union, through the Swedish Agency for Economic and Regional Growth and the Norrbotten County Council, as research project no. 20200552.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.