Numerical simulation and dimension reduction analysis of electromagnetic logging while drilling of horizontal wells in complex structures

Electromagnetic logging while drilling (LWD) is one of the key technologies of the geosteering and formation evaluation for high-angle and horizontal wells. In this paper, we solve the dipole source-generated magnetic/electric fields in 2D formations efficiently by the 2.5D finite difference method. Particularly, by leveraging the field’s rapid attenuation in spectral domain, we propose truncated Gauss–Hermite quadrature, which is several tens of times faster than traditional inverse fast Fourier transform. By applying the algorithm to the LWD modeling under complex formations, e.g., folds, fault and sandstone pinch-outs, we analyze the feasibility of the dimension reduction from 2D to 1D. For the formations with smooth lateral changes, like folds, the simplified 1D model’s results agree well with the true responses, which indicate that the 1D simplification with sliding window is feasible. However, for the formation structures with drastic rock properties changes and sharp boundaries, for instance, faults and sandstone pinch-outs, the simplified 1D model will lead to large errors and, therefore, 2.5D algorithms should be applied to ensure the accuracy.


Introduction
Electromagnetic (EM) logging while drilling (LWD) has been widely used in high-angle and horizontal (HA/HZ) wells for real-time geosteering and formation evaluation (Netto et al. 2013;Hu et al. 2017;Yuan et al. 2018;Lai et al. 2018). A fast and accurate EM LWD modeling algorithm is required for the pre-drilling forward simulation, the real-time inversion and the post-drilling evaluation in directional wells (Tan et al. 2012;Shao et al. 2013;Wang et al. 2019). The simulation of EM LWD in HA/HZ wells is essentially a 3D modeling problem (Deng et al. 2010(Deng et al. , 2015Zhang et al. 2019a, b). A variety of accurate 3D EM LWD modeling algorithms, e.g., finite difference method, finite element method and finite volume method (Gao et al. 2010;Rumpf et al. 2014;Zhang et al. 2019a, b), are therefore developed to address the 3D problem associated with tool structures, borehole and mud invasion. However, the model constructions and solutions of the 3D algorithms are quite complicated and time-consuming (Bensdorp et al. 2016), which restrict their applications in real-time data processing (Li and Wang 2016;Wu et al. 2017). To solve this problem, the formation is simplified to a 1D layered model through the sliding-window strategy in practice Hu and Fan 2018), and the coils are modeled by the magnetic dipoles (Hong et al. 2014). The real-time geosteering Edited by Jie Hao requirement can then be satisfied with the use of EM fields' closed form. This technique, however, lacks for feasibility and accuracy analysis toward complex structures, e.g., folds, faults, which can be formed under the in situ stress in realistic subsurface formation (Tan et al. 2019).
For some geological structures, the rock properties are continuous in strike but discontinuous in lateral direction within the well logging scope. Therefore, the structures like folds and faults can normally be simplified to 2D models in LWD modeling and simulated effectively and efficiently by the 2.5D algorithm (Dyatlov et al. 2015;Rodríguez-Rozas et al. 2018;Thiel et al. 2018). The basic principle of the 2.5D algorithm is to apply the Fourier transform to convert a 3D problem in spatial domain into a series of 2D problems in spectral domain, solve the fields in spectral domain and convert it back to the spatial domain through the inverse Fourier transform (Xu and Li 2018;Liu et al. 2012;Wu et al. 2019). To satisfy the real-time data processing and geosteering requirement, the key issues for the 2.5D modeling to answer are (1) how to rapidly solve the electric and magnetic fields in 2D spectral domain with arbitrary oriented dipoles; (2) how to improve the IFT efficiency based on the characteristics of the fields in spectral domain.
In this paper, a 2.5D finite difference method is proposed to model the EM LWD and applied to analyze the tool responses in complex geological structures. First, the dipoles in spatial domain are converted into incident fields in spectral domain as the effective sources. Second, the scattered field wave equations in spectral domain are derived from Maxwell's equations and then solved by the direct solvers. Finally, we propose a truncated Gauss-Hermite quadrature to accelerate the IFT, which transforms the solution in spectral domain to spatial domain. The 2.5D algorithm is then applied to compare the tool responses in complex reservoirs (folds, faults and sandstone pinch-outs) with those from the 1D model and analyze the 1D simplification's feasibility. Note that in this paper, we focus on modeling the complex geological structures while ignoring the borehole environment, tool structure and mud resistivity; the neglected effects are either negligible or can be corrected by calibration (Kang et al. 2018;Rosa et al. 2018).

Physics of EM LWD and formation model
The EM LWD provides measurements of multiple investigation depths by using multi-frequency and multi-spacing configurations. Figure 1a-I shows the basic structures of one transmitter-two receivers; Fig. 1a Fig. 1 Tool configuration of the EM LWD and the relationship between the phase shift (attenuation) and the formation resistivity: a basic structure with one transmitter and two receivers (I) and the tool structure of ARC675 developed by Schlumberger (II), b phase shift and c attenuation includes five basic measurement units (five axial transmitters and two coaxial receivers) and works at 2 MHz and 400 kHz. For an arbitrary basic measurement unit, the phase shift resistivity (PSR) (attenuation resistivity (ATR)) can then be obtained according to its single-value relationship with the phase shift (PS) (attenuation (AT)) of the voltages at the receivers, as shown in Fig. 1b, c. The PS and AT are defined as (Li and Wang 2016)

2.5D modeling of electric logging
Generally, electric logging modeling requires 3D model since the complex borehole and geological structures should be taken into account. In LWD modeling, the borehole and invasion effects are usually neglected; therefore, the model selection is highly dependent on the depth of investigation (DOI) of the tool and the shape of formation structure. (1) If the formation structure is complex within the tool detection ranges, a 3D model should be used (Fig. 2a).
(2) If the formation properties keep invariant along one direction within the tool detection ranges, a 2D model can be chosen (Fig. 2b). (3) If the bed boundaries are parallel to each other within the tool detection ranges, a 1D model can be a good choice (Fig. 2c). The 1D models have been widely used in real-time geosteering and fast inversion since the tool responses can be computed analytically with low CPU time and memory requirement. However, the application of the 1D model in complex formation reconstruction is restricted due to the over simplification. Therefore, the forward modeling and inversion of EM LWD in the 2D models gain more and more attentions. In this paper, we focus on the fast modeling of EM LWD in 2D formations, e.g., folds and faults, to provide technology supports for EM LWD evaluation in complex formations. The EM LWD modeling can be converted to the EM field computation problem with magnetic dipoles since the EM wavelengths are much larger than the coil sizes. For a lossy medium, the time harmonic Maxwell's equations can be expressed as where ε = (ε′ + iσ/ω) is the complex permittivity, M (J) is the magnetic (electric) current density. The EM fields in spatial domain can be split into background fields and scattered fields. Assuming a homogeneous isotropic background with permeability μ b and permittivity ε b , the background E b and H b field will follow Eliminating the source term by subtracting Eq. 3 from Eq. 2, we can obtain where the equivalent sources M sct and J sct are defined as

Fig. 2 Formation models
The EM wave equations can then be obtained as If the formation properties keep invariant along y direction, we can obtain the wave equations in spectral domain by applying the FT to Eq. 6 (Wu et al. 2019): where the right-hand side term is the equivalent source in spectral domain. At the working frequency of the well logging problems, electric fields and magnetic fields are coupled together. Therefore, either the electric wave equation or the magnetic wave equation can be selected to solve the EM fields accordingly. The selected wave equation is then discretized using a staggered finite difference grid and generates a linear system: where ̃ is the coefficient matrix obtained from the finite difference scheme, ̃ is the electric or magnetic field vector to be solved, and ̃ is the source term. Note that the 2.5D finite difference scheme and matrix assembly process are similar to the 3D method except that ∇ = ( x, y, z) should be replaced by ∇ = x, ik y , z . Finally, the field counterparts in spatial domain can be obtained by applying the IFT to the solved EM fields in spectral domain (Dyatlov et al. 2015):

Inverse Fourier transform
The 2.5D algorithm convert a 3D problem into a series of 2D problems of different wave numbers, and the number of 2D problem determines the simulation efficiency. The straightforward way to address the IFT is the IFFT, which can be accelerated by reducing the number of sampling points through the linear interpolation in spectral domain. Another class of approaches is to use numerical integration methods including Gauss-Legendre (GL) and Gauss-Hermite (GH), by viewing the IFT as an integral problem. GL quadrature is appropriate for definite integral, and therefore, the integration intervals should be determined first. GH demands the integrand decaying as exp(− x 2 ), but do not have to determine the integration ranges. In this paper, we propose a truncated GH quadrature method according to the decaying property of the fields in spectral domain.

Interpolation-fast Fourier transform
The sample density can be improved with cheap time cost by applying the linear interpolation in spectral domain.
Assuming that the interval between the ith and the (i + 1)th sampling point is split into (m + 1) segments, the field at the interpolation point can then be expressed as 2. Gauss-Legendre quadrature Equation 9 shows that the integral interval of IFT ranges from negative infinity to positive infinity, but it can be truncated into [k min , k max ] since the magnetic fields in spectral domain will approach to 0 as the wave number increases (see Sect. 3). The magnetic fields in spatial domain can be approximated as where k min and k max are the cutoff wave numbers.

Gauss-Hermite quadrature
Different from GL, the integration interval of GH quadrature is from negative infinity to positive infinity. Since the decaying properties of the magnetic fields in spectral domain are similar to the requirement of GH integrand, GH can be applied to the IFT of the magnetic fields with a slight modification as where n i is the ith Gauss node and k m is the truncation value of k y . The spatial magnetic fields can then be approximated as: where k yc = k y ·C is the Gaussian quadrature after the truncation and ω i is the weight of the ith GH node.

Analysis of fields in spectral domain
In this section, we analyze the characteristics of fields in spectral domain by taking magnetic dipoles-induced magnetic fields as example. We assume a z-oriented magnetic dipole located at (0,0) m inside a 10 Ω m homogeneous model radiates 400 kHz EM fields. Figure 3 shows the magnetic fields at z = 2.0 m (x is from − 10 to 10 m) vary with wave numbers and lateral distances in spectral domain, from which we can conclude that magnitude of the magnetic fields decay along both spatial direction (x) and wave number direction (k y ). Note that only one 2D problem on x-z plane needs to solve per wave numbers. Figure 4 shows the magnitude of H zz component on y = 0 plane. From the spatial perspective, magnetic fields attenuate along both x and z directions rapidly. From the wave number's perspective, the magnetic fields show similar distributions with different wave numbers, but the magnitude decays as the wave number increases. Aforementioned examples indicate that magnetic fields have good decaying properties in spectral domain. Note that this property is also valid for formations with complex structures. In sum, the wave number can be cut off when applying the IFT since the magnetic fields will reduce to 0 if the absolute of the wave number is large enough. The procedure to truncate the wave numbers is (1) determine the maximum magnitude value of the magnetic fields H m ; (2) find a minimum wave number k 0 to ensure the magnetic fields less than one thousandth of H m when the absolute of wave number is larger than |k 0 |; (3) choose the cutoff wave number k m as k 0 .

Fast solver for IFT
The EM fields in x-z plane can be obtained by solving the linear equations in spectral domain; the spatial domain fields can be further obtained by applying the IFT. To evaluate the efficiency of the IFT methods, we calculate the relative errors with different numbers of sampling points as where H 0 is the analytical solution of magnetic field tensor in spatial domain, H IFT is the magnetic field tensor obtained by IFT and (x, z) are the coordinate values. The resistivity of the formation is fixed at 5 Ω m, and the frequency of the tool is set to 24 kHz. Without loss of generality, an arbitrary point x = 0 and z = 9.56 m is chosen as an example to analyze the efficiency of the IFT methods. Figure 5a shows the relative errors of interpolation IFFT (IFFTI) decrease as the number of samples and number of interpolation points increase, indicating that the linear interpolation is a feasible way to improve the IFFT. Figure 5b shows the relative error comparison among IFFT, IFFTI, GL and GHT. The relative errors decrease as the number of samples increases for all four methods, where the slope of IFFT (GL) closes to IFFTI (GHT). Table 1 presents the fewest samples for all methods with the relative error lower than 1.0%, which shows that the efficiency of the different IFT methods follows GHT > GL > IFFTI > IFFT. The efficiency of GHT proposed in this paper is 45 times faster than the traditional IFFT because the number of samples decreases to 1/45. To further verify the conclusion, we test the IFT methods in another scenario. In this case, the formation resistivity and tool frequency are 10 Ω m and 400 kHz, and the observation point locates at (0, 2.0 m). The IFT results are shown in Fig. 5c, and the same conclusion as aforementioned can be obtained.

Responses in fold model
Fold is one of the most common geological structures that is difficult to be modeled by 1D methods due to its lateral In practice, the formation is split into a series of 1D model by using the sliding-window method: (1) assume a rectangular window with a tiny width along the direction vertical to the tool (Fig. 6b); (2) the bed boundaries within the window can be viewed as a bunch of straight line segments; (3) extend the segments to obtain the simplified 1D model (red dash lines in Fig. 6b). To investigate the feasibility of the simplification of fold, we build a three-layered model as shown in Fig. 6b. The middle bed is resistive of 1-m thickness, and the boundaries vary with cosine function of 20 m period (L). The EM LWD responses with 2 MHz-28 in (0.71 m) configuration are presented in Fig. 6a as the tool penetrates the fold horizontally. The red solid (dashed) curve and blue solid (dashed) curve represent, respectively, the apparent resistivity of attenuation, and the apparent resistivity of phase shift, while the black solid curve is the true model resistivity. We can see that the curves of PSR and model resistivity completely overlap within the horizontal section (0-2 m, 22-24 m). By contrast, ATR is larger than model resistivity, indicating that the ATR is significantly affected by the shoulder bed. The different behavior of PSR and ATR is because the investigation depth of the attenuation is much larger than the phase shift. Besides, both PSR and ATR from simplified 1D model agree well with the original 2D fold model, indicating that the simplification by using the sliding-window method is feasible. Another comparison of tool responses and simplified results for L = 10 m is shown in Fig. 7, from which we can also observe that 1D results agree well with 2D results.

Responses in fault model
Fault plays an important role in oil and gas reservoir; for one thing, hydrocarbons migrate along faults; for another, the fault sealing by mud daubing protects the reservoir. Three kinds of dip-slip fault models, vertical fault (Fig. 8b), reverse fault (Fig. 9b) and normal fault (Fig. 10b)

Responses in sandstone pinch-out
Sandstone pinch-out can be observed as hydrocarbon accommodation space. The main characteristic of the trap is the sedimentary bed thinning along the basin margin direction. A sandstone pinch-out in Fig. 11 includes a 20-Ω m oil/gas layer and a 3-Ω m water layer. A horizontal well penetrates the reservoir from the left to the right, 0.5 m from both the    Fig. 11a, b. Obvious curve separation can be observed when the tool penetrates the sandstone pinch-out, which demonstrates that the simplification introduces large errors.

Conclusions
In this paper, a 2.5D algorithm for arbitrary 2D models is applied to model the EM LWD efficiently in complex geological structures. Based on the decaying properties of the EM fields in spectral domain, we propose the GHT to truncate the integration of wave numbers self-adaptively; the efficiency can be improved several tens of times compared to the traditional IFT.
For the geological fold with smooth lateral changes in formation properties and bed boundaries, the 1D simplification with sliding windows can ensure both the modeling speed and accuracy for EM LWD modeling. By contrast, for the formation structures with drastic property changes or sharp boundaries like faults and sandstone pinch-outs, the simplified 1D model can cause large errors, and therefore, 2.5D algorithms should be applied to secure the accuracy. Besides, when the fault dip is identical with the well trajectory, the fault effects on EM LWD are larger than those when the fault dip is opposite to the well trajectory.
We also observe that when the tool is far away from the faults or non-parallel boundaries, the results from the simplified model can agree well with the original model. Therefore, it is a feasible way to simulate the EM LWD using the cross-dimensional method by combining the 1D and 2.5D algorithm to ensure both efficiency and accuracy.