In EM LWD measurements, fixed frequencies are typically used. In this paper, the time convention of e ^-iωt is assumed, and the variation in rock magnetic permeability is ignored. The Maxwell equations in the frequency domain for inhomogeneous medium can be expressed as ∇ × E = i ω μ 0 H (3) ∇ × H = J s + σ − i ω ε E (4) where, E is the Electric field, V; H is the magnetic fields, A⋅m⁻¹; ω is the angular frequency, rad⋅s⁻¹; J _s is the current source, A⋅m⁻²; σ is the conductivity tensor, S⋅m⁻¹; μ ₀ is the permeability of vacuum, 4π⋅10⁻⁷ H⋅m⁻¹; ε is the dielectric constant, F⋅m^-1. By combining Equations 3, 4, the wave equation of electric field can be obtained as shown in Equation 5. Moreover, the conductivity-permittivity coupling term (σ − iωε) can exhibit anisotropic properties, as characterized in Equation 7. ∇ × ∇ × E = i ω μ 0 J s + i ω μ 0 σ − i ω ε E (5)

The electromagnetic fields generated in the formation is the superposition of the incident field produced by a magnetic dipole source and the scattered field. Specifically, in a homogeneous isotropic formation, only the incident field, denoted as E ⁱ, exists. The scattered field, denoted as E ^s, is due to the heterogeneity of the formation. Therefore, the total field can be expressed as the sum of E ⁱ and E ^s. Substitute the total field E and the incident field E ⁱ into Equation 3 and take the difference to obtain Equation 6: ∇ × ∇ × E s − i ω μ 0 σ ^ E s = i ω μ 0 σ ^ − σ ^ 0 E i (6) where σ − i ω ε = σ x x − i ω ε x x σ x y − i ω ε x y σ x z − i ω ε x z σ y x − i ω ε y x σ y y − i ω ε y y σ y z − i ω ε y z σ z x − i ω ε z x σ z y − i ω ε z y σ z z − i ω ε z z (7)

Since subsurface rocks are often lossy media, the electromagnetic field will decay to zero far from the source. Therefore, the boundary condition of the EM LWD modeling problem becomes Equation 8: E x → ∞ , o r , y → ∞ , o r , z → ∞ = 0 (8)

In fact, considering the exponential decay of the electromagnetic field, using a Perfect Electric Conductor (PEC) boundary at a distance greater than 2.3 times the skin depth from the magnetic dipole source ensures that the boundary’s influence on the simulation results is less than 1%.

Taking into account factors such as the borehole condition, mud invasion and formation structure, a three-dimensional model as shown in Figure 2 is established. For the three-dimensional simulation of EM LWD, we typically need to focus on factors such as mud resistivity, wellbore dimensions, mud invasion, and the environment of deviated or horizontal wells. The dashed lines in Figure 2 represent the finite difference grid, which divides the formation model into a series of small cells. By assigning corresponding electrical parameters to each cell, further difference simulations can be conducted. The figure shows a vertical well model under inclined formation conditions. By rotating the coordinates, it can also be used to simulate inclined well conditions in horizontal formations.

It should be noted that formation conductivity is usually expressed in the formation coordinate system. In EM LWD modeling, it is necessary to rotate the formation coordinate system to the tool coordinate system as expressed in Equation 9: σ t = R − 1 σ f R (9) where, σ _t and σ _f are the conductivity matrices in the tool coordinate system and the formation coordinate system, respectively. R is the rotation matrix, which is shown in Equation 10. R = cos ⁡ θ ⁡ cos ⁡ φ ⁡ cos ⁡ γ − sin ⁡ φ ⁡ sin ⁡ γ − cos ⁡ θ ⁡ cos ⁡ φ ⁡ sin ⁡ γ − sin ⁡ φ ⁡ cos ⁡ γ sin ⁡ θ ⁡ cos ⁡ φ cos ⁡ θ ⁡ sin ⁡ φ ⁡ cos ⁡ γ + cos ⁡ φ ⁡ sin ⁡ γ − cos ⁡ θ ⁡ sin ⁡ φ ⁡ sin ⁡ γ + cos ⁡ φ ⁡ cos ⁡ γ sin ⁡ θ ⁡ sin ⁡ φ − sin ⁡ θ ⁡ cos ⁡ γ sin ⁡ θ ⁡ sin ⁡ γ cos ⁡ θ (10) where θ , φ and γ are the three Euler angles.

To apply the finite difference method, the electric field wave equation is first expanded into a scalar equation: ∂ 2 E y s ∂ x ∂ y + ∂ 2 E z s ∂ x ∂ z − ∂ 2 E x s ∂ y 2 − ∂ 2 E x s ∂ z 2 − i ω μ 0 σ x x E x s + σ x y E y s + σ x z E z s = i ω μ 0 σ x x − σ 0 E x i + σ x y E y i + σ x z E z i ∂ 2 E x s ∂ x ∂ y + ∂ 2 E z s ∂ y ∂ z − ∂ 2 E y s ∂ x 2 − ∂ 2 E y s ∂ z 2 − i ω μ 0 σ y x E x s + σ y y E y s + σ y z E z s = i ω μ 0 σ y x E x i + σ y y − σ 0 E y i + σ y z E z i ∂ 2 E x s ∂ x ∂ z + ∂ 2 E y s ∂ y ∂ z − ∂ 2 E z s ∂ x 2 − ∂ 2 E z s ∂ y 2 − i ω μ 0 σ z x E x s + σ z y E y s + σ z z E z s = i ω μ 0 σ z x E x i + σ y z E y i + σ z z − σ 0 E z i (11)

It can be observed that the incident and scattered fields of the electric field are separated to the left and right sides of the equation, respectively. The scattered field is the unknown to be solved, while the incident field can be obtained through an analytical solution. Therefore, a linear system of equations can be formed by Equation 11. By using the finite difference method, the differential on the left side of the equation is converted into a difference, thereby obtaining the coefficient matrix for the linear system of equations. For the second-order difference, the coordinates at node (i, j, k) are defined as (x _i , y _j , z _k ). The coordinates x _i+1/2, y _j+1/2, and z _k+1/2 correspond to the positions at the midpoints of the cell edges between nodes in the x, y, and z directions, respectively. The distances between two nodes in the x, y, and z directions are Δx _i = x _i+1-x _i ,Δy _j = y _j+1-y _j ,Δz _k = z _k+1-z _k . The distances from the node centers are Δx _i-1/2 = x _i+1/2-x _i-1/2,Δy _j-1/2 = y _j+1/2-y _j-1/2,Δz _k-1/2 = z _k+1/2-z _k-1/2. The derivation of finite difference schemes is a nontrivial process. Without loss of generality, we take the first term on the left side of Equation 11 as an example. Its finite difference scheme is formulated in Equation 12: ∂ 2 E y s ∂ x ∂ y = ∂ ∂ y E y i + 1 , j , k s − E y i , j , k s ∆ x i = 1 ∆ x i ∆ y j − 1 2 × E y i + 1 , j + 1 2 , k s − E y i , j + 1 2 , k s − E y i + 1 , j − 1 2 , k s − E y i , j − 1 2 , k s (12)

The other terms on the left side of the equation can be obtained using the same method. For the terms on the right side, we only need to calculate the incident field at the corresponding nodes using an analytical solution. Another issue in constructing the linear system of equations is the assignment of conductivity values. Since the electric field is located on the edge of a cell and there are four cells surrounding that edge, it is necessary to use the properties of the four cells to equivalently assign the conductivity of the electric field at that location. Here, we use the xx component as an example to derive the equivalent conductivity. Assume that the electric field within any cell is uniform and equal to the electric field at (i+1/2, j, k). For the electric field in the x-direction, the current generated in the x-direction is given by J x x i , j − 1 , k − 1 S i , j − 1 , k − 1 + J x x i , j , k − 1 S i , j , k − 1 + J x x i , j − 1 , k S i , j − 1 , k + J x x i , j , k S i , j , k = J x x i + 1 2 , j , k S (13) where S represents the cross-sectional area. By combining the differential form of Ohm’s law with Equation 13, we obtain the equivalent conductivity as expressed in Equation 14. σ x x i + 1 2 , j , k = ∆ z k − 1 2 ∆ y j − 1 σ x x i , j − 1 , k − 1 + ∆ y j σ x x i , j , k − 1 4 ∆ y j − 1 2 ∆ z k − 1 2 + ∆ z k ∆ y j − 1 σ x x i , j − 1 , k + ∆ y j σ x x i , j , k 4 ∆ y j − 1 2 ∆ z k − 1 2 (14)

The equivalent conductivity for other components can be derived in a similar manner and will not be repeated here.

Typically, due to the complex formation structures and measurement environments, the electromagnetic field distribution generated by the EM LWD tools is quite intricate. However, when the complexity of the formation structure is reduced, and only horizontal layered formations are considered, the electromagnetic field distribution will exhibit certain symmetric or antisymmetric properties. Consider a deviated borehole model within a horizontally layered formation. The top and bottom shoulder beds are modeled as semi-infinite, isotropic formations, with resistivities of 3 Ω·m and 4 Ω·m, respectively. The middle layer is characterized by anisotropic properties, with a horizontal resistivity of 10 Ω·m and a vertical resistivity of 40 Ω·m. The wellbore deviates at an angle of 60° relative to the formation normal, with a radius of 4.25 inches (10.8 cm), and is filled with drilling mud of resistivity 0.1 Ω·m. Mud invasion occurs within the middle layer, extending to a depth of 0.5 m. The invaded zone exhibits a horizontal resistivity of 1 Ω·m and a vertical resistivity of 4 Ω·m. Figure 3 illustrates the cross-sectional view of the electric field in the YOZ plane generated by a magnetic dipole source oriented in the z-direction. It is evident that all components of the electric field are symmetric about the y = 0 plane. Much more numerical simulations indicate that similar results can be observed for magnetic dipole sources oriented in other directions. This symmetry arises from the fact that, although the current model is three-dimensional, it is symmetric about the y = 0 plane. Therefore, the symmetry of the field is an inherent property of the model.

To further validate the symmetry of the electromagnetic field, a three-dimensional model is established. Specifically, a vertical well is embedded in a homogeneous formation with a resistivity of 10 Ω·m. The borehole has a diameter of 8.5 inches and is filled with mud of resistivity 0.1 Ω·m. The transmitter coil is positioned at the center of the well axis at a relative depth of 0 m, with coordinates (0, 0, 0). A low-resistivity anomaly, characterized by a resistivity of 1 Ω·m and a side length of 2 m, is placed adjacent to the well, with its center located at coordinates (2 m, 0 m, 0 m). Figure 4 presents the cross-sectional distribution of the Ez component on the XOY plane (z = 0.5) generated by different directional dipole sources. The results indicate that the electric field maintains symmetry with respect to the y = 0 plane.

Since the electromagnetic field exhibits symmetrical properties, we can simplify the simulations by computing the field in half of the space. Figure 5a presents a schematic of the electric field nodes on the YOZ plane within a full-space computation domain, while Figure 5b illustrates the corresponding electric field node distribution on the YOZ plane in a half-space computation domain. In this case, we assume the midpoint of the y-axis coordinates y N y + 1 2 and y N y + 1 2 + 1 of the symmetry plane. In half-space modeling, a crucial consideration is that by eliminating half of the spatial domain, the boundary conditions of the computational model are consequently altered. Specifically, in the full-space 3D finite difference simulation of EM LWD, PEC boundary conditions are used. However, in the case of half-space modeling, one of the boundaries is transformed into a symmetry plane (i.e., the YOZ plane). On this symmetry plane, the electric field cannot be treated with PEC boundary conditions. In contrast, we can leverage the symmetry of the field to establish new boundary conditions. Assuming the number of grid points in the x, y, and z directions are denoted as Nx, Ny, and Nz, respectively, with each being an odd number. The boundary conditions on the symmetry plane can be set according to Equation 15: E x i + 1 2 , N y + 1 2 , k = s g n M · E x i + 1 2 , N y + 1 2 + 1 , k E z i , N y + 1 2 , k + 1 2 = s g n M · E z i , N y + 1 2 + 1 , k + 1 2 (15) where sgn (·) is a function of the magnetic dipole moment. For magnetic dipole sources oriented in the x and z directions, sgn (·) equals −1, while for those oriented in the y direction, sgn (·) equals +1. The relative position relationship of the electric field in the equation can be seen in Figure 4B. It should be noted that on the YOZ plane, only electric field components in the x and z directions are present. Therefore, boundary conditions need only be considered for E _x and E _z . Since the computational domain has been reduced to half of its original size, the number of electric field nodes to be solved is similarly halved. Consequently, the computational efficiency can be significantly improved.

During the drilling process, the wellbore is typically filled with drilling mud. As a result, the electromagnetic waves generated by EM LWD tools are inevitably influenced by the borehole environment. In this section, we take the Hzz (coaxial component) and Hzx (cross-coupling component) components as examples to analyze the effects of the borehole on the tool response. Consider a three-layer formation model, where the middle resistive layer exhibits a resistivity of 10 Ω·m and a thickness of 3 m. The top and bottom layers are infinitely thick with a resistivity of 1 Ω·m. The formation is penetrated by a vertical borehole. In the first scenario, we ignore the influence of the borehole and consider only the variation in formation resistivity. In the second scenario, we take the borehole effects into account, with a borehole radius of 6 inches and a mud resistivity of 0.2 Ω·m. In this case, we assume a one-transmitter, one-receiver configuration, with an operating frequency of 20 kHz and a transmitter-receiver spacing of 40 inches. Figure 7 illustrates the variation of the imaginary part of the Hzz component with depth in the formation. In Figure 7a, the comparison between the FDFD results and the analytical solutions demonstrates a good agreement, confirming the validity of the 3D finite difference algorithm. Figure 7b compares the response results with and without the presence of a borehole, revealing that the existence of the borehole leads to deviations in the Hzz curve. Since the Hzz component is typically used for measuring apparent resistivity in formations, it is important to appropriately consider the borehole environment during apparent resistivity analysis. Figure 8 illustrates the variation of the Hxz component with depth (in this case, the formation dip angle is 45°, and the other parameters are the same as in Figure 7). Separation of the curves can be observed as the tool approaches the formation boundary. The Hxz component is commonly used for detecting bed boundaries, and significant borehole effects may lead to errors in determining the distance to the boundary. It should be noted that, in most cases, a borehole radius of 4.25 inches is typically used rather than 6 inches. However, in this example, our primary goal is to demonstrate the validity of the method and the overall impact of the borehole on the curves. Therefore, slight variations in borehole size will not affect the conclusions.

When discussing borehole effects, it is essential to consider the influence of mud resistivity. The resistivity of drilling mud can vary significantly; for instance, saline-based muds exhibit very low resistivity, while freshwater and oil-based muds demonstrate much higher resistivity. Consider a homogeneous formation with a resistivity of 10 Ω·m, penetrated by a wellbore with a radius of 4.25 inches. The resistivity of the mud within the wellbore varies from 0.002 to 100 Ω·m. Here, we focus on examining the impact of wellbore mud on the apparent resistivity curve. Figure 9 illustrates the variation in amplitude ratio and phase shift with mud resistivity for electromagnetic logging-while-drilling (LWD) configurations of 2 MHz–28 in and 400 kHz-40 in. As shown in Figure 9a, the amplitude ratio for the 2 MHz-28 in configuration decreases rapidly with increasing mud resistivity, eventually leveling off, indicating the response of a homogeneous medium. In particular, when the mud resistivity is as low as 0.002 Ω·m, the amplitude ratio of the 400 kHz-40 in configuration remains approximately equal to that of a homogeneous formation. In contrast, when the mud resistivity falls below 0.01 Ω·m, the amplitude ratio of the 2 MHz-28 in configuration is significantly affected by the borehole. In Figure 9b, the two curves exhibit singularities in the low-resistivity (high-conductivity) mud region, which may be attributed to the nonlinear phase shift induced by the presence of highly conductive mud. For the 400 kHz-40 in phase difference, the measurement results remain largely unaffected by borehole conditions when the mud resistivity exceeds 0.05 Ω·m. In contrast, the 2 MHz-28 in phase difference remains susceptible to borehole effects unless the mud resistivity exceeds 0.2 Ω·m. Overall, high-frequency phase difference signals exhibit greater sensitivity to low-resistivity mud conditions, whereas low-frequency amplitude ratio signals demonstrate higher robustness against such influences.

In the process of oil and gas well drilling, overbalanced drilling is commonly employed, where the wellbore pressure exceeds the formation pressure, leading to drilling mud invasion into the formation. Considering that logging-while-drilling (LWD) incorporates real-time measurement during the drilling process, the effects of mud invasion are typically neglected. In this paper, we utilize numerical simulation to rigorously investigate the impact of mud invasion on the response characteristics of various logging tools.

First, we consider the impact of mud invasion on apparent resistivity measurements obtained from a conventional EM LWD tool. Two mud invasion models are established as follows: (1) Conductive invasion model: The wellbore radius is 4.25 in, with a mud resistivity of 0.5 Ω·m, a formation resistivity of 10 Ω·m, and an invasion zone resistivity of 1 Ω·m. (2) Resistive invasion model: The wellbore radius is 4.25 in, with a mud resistivity of 1,000 Ω·m, a formation resistivity of 1 Ω·m, and an invasion zone resistivity of 10 Ω·m. Figure 10 illustrates the variation of A28H, P28H, A40L, and P40L with respect to the invasion depth. In this context, A28H and P28H represent the amplitude ratio and phase shift signals from the 2 MHz-28 in configuration, while A40L and P40L represent the amplitude ratio and phase shift signals from the 400 kHz-40 in configuration, respectively. Here, ‘A’ and ‘P’ denote the amplitude ratio and phase shift, respectively, to distinguish the types of signals, while ‘H’ and ‘L’ stand for high frequency and low frequency, respectively, to differentiate between the 2 MHz and 400 kHz signals. It can be observed that conductive invasion increases both the amplitude ratio and phase shift, while resistive invasion results in a decrease in both parameters. The amplitude ratio shows the most significant variation at invasion depths between 0.5 and 1.0 m, while the phase shift exhibits the most pronounced variation between 0 and 0.5 m. When the invasion depth exceeds 1.0 m, neither the amplitude ratio nor phase shift undergoes significant changes with further increases in invasion depth, primarily because the invasion zone extends beyond the detection range of the EM LWD tools.

Using the same mud invasion models as previously described, we further analyze the variation in the responses of an extra-deep EM LWD tool with different measurement modes as a function of invasion depth, taking the 24 kHz-13 m coil configuration as an example. The numerical simulation results are shown in Figure 11.

The figure demonstrates that mud invasion does not significantly affect the responses of USDA, USDP, UADA, UADP, UHAA, and UHAP. Although mud invasion causes a slight variation in the values of UHRA and UHRP, the effect is minimal, especially when the invasion depth is less than 1 m, where the influence can be considered negligible.

Numerical simulations indicate that the effect of mud invasion on tool responses is strongly dependent on the tool parameters and the type of invasion. For high-frequency, short-spacing signals, mud invasion leads to a deviation of the apparent resistivity from the true formation resistivity. In contrast, for low-frequency, long-spacing signals, the impact of mud invasion is minimal and can generally be disregarded in well logging analysis.

To further illustrate the applicability of the proposed method, we developed the model presented in Figure 12 to investigate the effect of a resistivity anomalous body adjacent to the borehole on EM LWD responses. In this model, a homogeneous sandstone formation with a resistivity of 10 Ω·m is assumed, incorporating an infinitely long, low-resistivity shale inclusion with a cross-sectional dimension of 2 m × 2 m. A horizontal well traverses beneath the shale, with the borehole positioned at a distance D2A of 1 m from the lower boundary of the shale.

Figure 13 illustrates the tool responses obtained from EM LWD, where GA96L and GA96M denote the amplitude ratio signals at 96 in-100 kHz and 96 in-400 kHz, respectively, while GP96L and GP96M represent the phase shift signals at 96 in-100 kHz and 96 in-400 kHz, respectively. It is well established that when the tool operates in a homogeneous formation, geosignals are equal to zero. The detection capability of the tool is typically defined by threshold values of 0.25 dB for the amplitude ratio and 1.5° for the phase shift. As shown in Figure 13a, when adopting 0.25 dB as the threshold, GA96L can detect the presence of the anomalous body at a distance of 2.5 m from its left boundary, whereas GA96M can only detect it within 1.5 m of the left boundary. This observation indicates that lower-frequency signals exhibit greater detection depth. In Figure 13b, due to the relatively small magnitude of the phase shift signal, GP96L fails to detect the anomalous body when using 1.5° as the threshold, while GP96M is capable of detecting it only when the distance to the left boundary is reduced to 0.5 m. These results demonstrate that the amplitude ratio signal provides superior detection capability compared to the phase shift signal.