Assuming harmonic time dependence e − i ω t , Maxwell’s equations for the EM fields in the FD can be written as (Streich, 2009) ∇ × E = − i ω μ 0 H , (1) ∇ × H = J + σ * E , (2) where i² = −1, ω = 2 π f is the angular frequency, f is the frequency, E and H denote the total electric and magnetic fields, μ ₀ within the Earth is assumed to be constant and set to that of free space, J is the source current distribution, σ * = σ + i ω ε is the ohmic conduction term, σ is the electric conductivity, ε represents the air dielectric constant, and σ * E describes the induced currents inside the Earth. According to Eqs 1, 2, the expression of the FD vector Helmholtz equation for the scattered electric field is as follows (Alumbaugh et al., 1996): ∇ × ∇ × E s − k 2 E s = k 2 − k p 2 E p , (3) where k 2 = − i ω μ 0 σ + i ω ε , the superscripts “p” and “s” represent the primary and scattered fields, respectively, and E ^s and E ^p denote the scattered and primary electrical fields, respectively. Eq. 3 forms the basis of the presented modeling scheme and is well suited for simulations in CSEMs. Eq. 3 can overcome the difficulty of source description and the influence of source singularity on the numerical stability. Regardless of the type of source used to create the field, the procedure for calculating the scattered field remains the same. Then, we obtain the following PDE for the component i ^ : ∂ 2 E y s ∂ x ∂ y − ∂ 2 E x s ∂ y 2 − ∂ 2 E x s ∂ z 2 + ∂ 2 E z s ∂ x ∂ z − k 2 E x = k 2 − k p 2 E x p . (4)

In this paper, we adopt the ED or MD as the basic composition unit for any complex geometry source (Figure 1A). For example, the signal emission source of the LOTEM can be decomposed into many EDs (Figure 1B), and vector decomposition is applied to each ED (Figure 1C). The EM fields at any given position can be obtained by superimposing the EM fields generated by the ED component along the x, y, and z directions. Therefore, the primary objective of 3D forward modeling in CSEMs, regardless of the specific source type used, is to solve the electric field of an electric dipole or a magnetic dipole within the background models.

The primary field, as described by Weng et al. (2016), employs a virtual interface technique (Das and De Hoop, 1995) to solve the whole-space EM fields in a 1D layer model for different types of sources (Figure 2). These sources include vertical and horizontal electric dipoles (VED and HED) as well as vertical and horizontal magnetic dipoles (VMD and HMD).

By separating the partial wave solutions of the Helmholtz equations into upward and downward waves within certain boundaries, the potentials for Green’s function are obtained. Starting from the source level, the amplitudes of the potentials in each layer are derived recursively based on the initial amplitudes. For different types of sources, only the initial terms that are associated with the transmitting sources need to be modified, and the kernel connected to the layered media remains the same. Hence, the aforementioned scheme can be easily applied to EM transmitting sources with slight modifications.

To calculate the scattered electric field in the medium, the geoelectric model is discretized by the cuboid cell (including the air layer) (Figure 3A), each denoted by a subscript (i, j, and k). The resistivity of each cell is represented by the symbol ρ (i, j, and k). i, j, and k represent the mesh indices in the x, y, and z directions, ranging from 1 to Nx, Ny, and Nz, respectively. Nx, Ny, and Nz denote the numbers of cells in the x, y, and z directions, respectively. Therefore, Eq. 4 can be discretized using the 3D staggered-grid finite-difference method, and the discretized expression is given in Eq. 5. The electric field components are defined on cell edges, while the magnetic field components naturally correspond to the cell faces (Figure 3B). The discrete finite difference equations of the j ^ and k ^ components are given in Appendix A. [ 1 Δ x i E y i + 1 , j + 1 / 2 , k s − E y i , j + 1 / 2 , k s − 1 Δ y j E x i + 1 / 2 , j + 1 , k s − E x i + 1 / 2 , j , k s − 1 Δ x i E y i + 1 , j − 1 / 2 , k s − E y i , j − 1 / 2 , k s − 1 Δ y j − 1 E x i + 1 / 2 , j , k s − E x i + 1 / 2 , j − 1 , k s ⋅ 1 Δ y ¯ j + [ 1 Δ x i E z i + 1 , j , k + 1 / 2 s − E z i , j , k + 1 / 2 s − 1 Δ z k E x i + 1 / 2 , j , k + 1 s − E x i + 1 / 2 , j , k s − 1 Δ x i E z i + 1 , j , k − 1 / 2 s − E z i , j , k − 1 / 2 s − 1 Δ z k − 1 E x i + 1 / 2 , j , k s − E x i + 1 / 2 , j , k − 1 s ⋅ 1 Δ z ¯ k + i ω μ p y ^ p E x i + 1 / 2 , j , k s = − i ω μ p y ^ i + 1 / 2 , j , k − y ^ p E x i + 1 / 2 , j , k p , (5) where y ^ = σ + i ω ε , w = x, y, z, and l = i, j, k denote the length of the lth cell in the w direction, and ∆ w l ¯ is the distance between the center of the lth and l1th cells in the w direction.

Equation 5 can be assembled into the following system: K E s = b , (6) where K is the coefficient matrix, and b in the right-hand side consists of the dot products between the primary/background field and conductivity abnormalities, as well as the appropriate boundary conditions.

Generally, the coefficient matrix K for Eq. 6 is a large, sparse, and ill-conditioned matrix that is difficult to solve. To reduce the condition number of the coefficient matrix K, Eq. 6 can be written in the preprocessing form as follows: K M − 1 ⋅ y = b . (7) Here, y = ME ^s is the modified unknown vector, where the matrix M is called the (right-hand side) preprocessor, and the matrix KM ^-1 is considered very close to the identity matrix. Once Eq. 7 is solved approximately, we can obtain E ^s using the relationship between y and E ^s. In this study, an incomplete Cholesky decomposition is used as the precondition to accelerate the convergence and improve the accuracy of the iterative solution. Eq. 7 can be solved by using the quasi-minimal residual method (QMR) (Mackie et al., 1994).

Due to the significant difference between the solution region and the anomalous bodies, the secondary field will decay to zero at the boundary far from the anomalous body. As a result, we applied typical homogeneous Dirichlet’s boundary conditions. With n being the normal vector on the domain boundary δΩ, it is defined as: n × E s = 0 . (8)

The grid generation is centered on the area of interest and is divided into finer grids near the anomaly bodies. As depicted in Figure 3A, the grid size gradually increases as the distance from the area of interest expands. The air layer is set to seven layers, with the uppermost layer having a thickness of 30 km. In our paper, the mini software CSEM mesh is used to generate the mesh, which was developed by teachers and students in our research group. It greatly reduces the cost of generating meshes.

The divergence of the secondary electric field E ^s is zero at any point except the source (Shen, 2003; Chen et al., 2011). Due to the accuracy of the numerical calculation, the divergence of E ^s does not disappear during an iteration solution process which can be calculated using φ = ∇ ⋅ E s , considering the additional electrical field generated by the source φ. To calculate the additional electrical field, it is necessary to use the Neumann boundary conditions to solve the potential obtained by Poisson’s equation as follows: ∇ ⋅ ∇ ψ = φ . (9)

When Eq. 6 is solved, the E ^s values required for divergence correction can be expressed as follows: E c o r r e c t s = E o l d s − ∇ ψ , (10) specifically, to satisfy the following equation: ∇ ⋅ E c o r r e c t s = ∇ ⋅ E o l d s − ∇ ⋅ ∇ ψ = φ − φ = 0 . (11)

The preconditioned conjugate gradient (PCG) algorithm is used to solve the divergence in Eq. 11. The iterative solution’s convergence rate is greatly improved, particularly at low frequencies (Smith, 1996a; Smith, 1996b).

Once E ^s on the center of the grid edge is obtained by Eq. 6, the total electric field E can be estimated by adding a primary field E ^p, that is, E = E p + E s , (12) while the total magnetic field normal to the surface confined by the grid edge at the center can be approximated by Faraday’s law: − i ω B = ∇ × E , (13) from the estimated EM fields. Then, the EM fields at the position of interest can be interpolated by E = L e E (14) and H = L h H , (15) where L _e and L _h represent bilinear splines from the 3D grid nodes and edges, respectively, to the data sites.

According to Weng et al. (2017), the TD EM signal h(t) through h t = 1 2 π ∫ 0 ∞ I ω H ω e − i ω t d ω (16) is related to the FD response H(ω) from a source excited by current I(ω) over a conductive model; in the aforementioned equation, i² = −1, ω = 2 π f is the angular frequency, f is the frequency, and t is the time. If a step-off current is assumed, using the Euler formula, Eq. 16 can be further split into a sine transformation as follows: d B / d t = μ 0 d h / d t = − μ 0 2 π ∫ 0 ∞ I m H ω sin ω t d ω , (17) where B is the magnetic flux density and Im [·] denotes the imaginary part of the FD EM fields. We sample the exact six points per decade FD data between 10^–3 Hz and 10⁸ Hz and obtain 67 frequencies (Liu et al., 2016). After obtaining the FD results, the TD EM responses for the step wave are calculated using Eq. 17. Based on the aforementioned theory, we have programmed a 3D forward modeling code using Fortran 90 for FD/TD CSEMs with multitype sources.

In real field surveys, surface obstacles and topography can preclude laying out the source in a theoretical shape such as a straight line, rectangle, or circle. The aforementioned scheme can segment a source with an arbitrarily complex shape into a large number of EDs and MDs. To illustrate this, we conducted a study taking the LOTEM and CSAMT as examples, focusing on the differences between the TD and FD EM fields for both straight and non-straight wire sources. We computed the TD and FD responses for the bent and straight wires. Both wires centered at (0, 0) have the same grounding points; the straight wires of 1 km have a 1 A current on the surface, and the bent wire is the real source consisting of 22 segments (Figure 6A). The 3D model of a 100 Ωm homogeneous half-space is shown in Figure 6B.

In Figure 7, we display the TD magnetic field dB/dt for a bent wire compared to dB/dt for a straight wire at different times: 0.065 ms, 0.9 ms, 4 ms, and 10 ms. As time increases, the difference in dB/dt between the bent and straight wires decreases. dB/dt generated by the straight wire exhibits a symmetrical distribution around the source center and is concentrated near the surface (Figure 7A). At the early stage, there is a significant disparity in the dB/dt between the bent and straight wires (Figure 7C). Therefore, when analyzing early-time TEM data, it is crucial to consider the shape of the source. With a time delay, the center of dB/dt propagates downward and outward, and the strength of the field gradually decays. At the later stage, dB/dt continues to propagate and becomes more uniform, and the relative differences in dB/dt between the bent and straight wires are minimal (Figure 7C). Spatially, the closer the distance to the source, the greater the difference in dB/dt between the bent and straight wires. Conversely, the farther the distance to the source, the smaller the difference in dB/dt between the bent and straight wires. Therefore, when processing data of TEMs observed close to the source, such as the SOTEM, the source’s shape must be taken into account.

Figure 8 shows the FD electric field component Ey for the bent wire relative to Ey for a straight wire at frequencies of 0.1 Hz, 2 Hz, 16 Hz, and 128 Hz. The relative differences in Ey between the bent and straight wires increase with increasing frequency and decrease with increasing distance from the source. At low frequencies, the difference in Ey between the bent and straight wires was very small, and it was mainly concentrated near the source and both sides along the emitting source direction. At high frequencies, the difference in Ey was most pronounced in the survey area. These results highlight the significant influence of the grounding points’ locations and the entire wire layout on the EM fields at these frequencies. Therefore, during the inversion of the high-frequency data from a wire source, such as the CSAMT, it is impossible to ignore the impact of the source’s shape on the measured data.

The comparison between the TD and FD EM methods shows that the FD CSEM appears to be more significantly influenced by the entire wire layout than the TD CSEM. This difference can be attributed to the nature of the measurements. In the TD CSEMs, observations are made after the current is turned off, capturing the induced eddy field or secondary field. In this case, the shape of the emission source only affects the early data. On the other hand, in FD CSEMs, observations are made under the harmonic excitation of the source, representing the total field. The impact of the source’s shape on the field becomes more noticeable, particularly at higher frequencies. At low frequencies, the electric field resembles a potential field that occurs in the static limit. Here, the only factors affecting the electric field are the locations of the grounding points.

According to the aforementioned analysis, we suggest that the shape of the source must be considered for TD and FD CSEMs with a wire source during data processing and inversion. However, in regard to the influence of the source shape, we recommend TD CSEMs when processing data based on the theory of an ideal straight wire or electric dipole source for the same target.