Using logging while drilling resistivity imaging data to quantitatively evaluate fracture aperture based on numerical simulation


 Fractured formations are strongly heterogenous, and thus exhibit a complex logging response mechanism. By using the logging while drilling (LWD) resistivity imaging tool, fractures can be visually identified and their aperture quantitatively calculated. Because physical fracture model simulation is time consuming and costly, we propose using a 3D finite element method (FEM) numerical simulation to interpret the LWD resistivity imaging tool logging responses in conjunction with a new aperture calculation model based on the forward model. First, we used the single fracture model to investigate the effect of fracture aperture and formation resistivity contrast on the maximum current contrast at the fracture. The results showed that the aperture is linearly related to the maximum current contrast, while the formation resistivity contrast exhibits a pronounced exponential relationship with the maximum current contrast. Both of these relationships are affected by the fracture's dip angle, so segmented fitting is required when the fracture dip angles differ. Next, using the forward model, we developed the fracture aperture calculation model based on the maximum current contrast. The aperture calculation model was established in three segments in accordance with the different fracture dips, and the influence factors affecting the fracture inverting inclination were analyzed using multi-fracture simulation images. Finally, the accuracy of the new model was verified with the simulated fracture images. The novel model for calculating fracture aperture is of great significance for processing and interpreting LWD resistivity imaging logging data.


Introduction
Resistivity logging plays an important role in formation evaluation during oil and gas exploration (Chen 2017;. Before the 1980s, resistivity logging measurement results were mainly two-dimensional log-ging curves, which were used to evaluate formation fluid saturation. Since then, microresistivity image logging has been introduced, which involves obtaining three-dimensional, high resolution imagery of microresistivity variations around the borehole wall. The imagery resolution is on the order of a few millimeters in both the vertical and azimuthal directions, which is very useful for evaluating the natural fracture parameters, such as aperture, density and porosity. Thus, this logging method opened a new chapter in logging methodology and interpretation (Ekstrom et al. 1986;Luthi 2001;Keeton et al. 2015;Sun et al. 2016;Brekke et al. 2017;Lai et al. 2018).
With the development of logging while drilling (LWD) technology, the microresistivity image logging tool was superseded by the LWD resistivity image tool. The latter generally uses the toroidal coils for excitation, and its current focus principle is similar to wireline logging of borehole resistivity images (Arps 1967). The tool generally has 1-4 electrode pads, which is substantially fewer than the microresistivity imaging logging tool. Measurements of borehole wall resistivity images are obtained through drill pipe rotation, with borehole coverage that can reach 100%. So far, the established tools include azimuthally focused resistivity or AFR (Halliburton), Geovision or GVR (Schlumberger) and Startrak (Baker Huges) (Akimov et al. 2007;Prammer et al. 2009;Ortenzi et al. 2012). The new tool benefits from advances in microsignal detection technology and uses smaller buttons to obtain higher resolution images matching the Formation Mi-croImager (FMI) (Allouche et al. 2010). Moreover, it offers borehole imaging data from underground to overground in real time, enabling real-time geological guidance and reducing of drilling risk.
Based on these advantages, the tool is commonly used in complex formations, as it gets good application results, especially in carbonate and volcanic formations with high variability where it is difficult to predict logging response characteristics (Prammer et al. 2007;Tan et al. 2011;Fan et al. 2016;Wang & Thagafy 2017). In these reservoirs, fractures are storage spaces and connection channels for fluid and gas (Nelson 2001). This tool is very sensitive to conductive fractures, so it is suitable for use in water-based mud (Maeso et al. 2014). A logging interpretation expert can easily identify the fractures and calculate apertures from the borehole images. This is highly relevant, as the natural fracture's aperture is the basis for calculating fractured reservoir reserves and fracturing (Rajabi et al. 2010;Ponziani et al. 2013;Maeso et al. 2014Maeso et al. , 2015. The current fracture aperture calculation model is mostly based on the 3D forward model. Luthi & Souhaité (1990) used 3D FEM to calculate the formation microresistivity scanner (FMS) logging response in a fractured formation, and the fracture aperture model was established by the abnormal current concepts at a single fracture. Recently, new established physical models have validated the model's accuracy (Ponziani et al. 2015). The existing calculation model of fracture aperture is obtained when the fracture's dip angle is 0°. In fact, fractures with 0°and 40°are uncommon in many reservoirs and steeper angles are particular abundant in extensional and strike-slip situations (Luthi & Souhaité 1990). Therefore, the calculation error increases when the fracture dip angle is higher than 0°. At the same time, the actual calculation needs to integrate the area of abnormal current at the fracture, which complicates the calculation. In this work, we simulated the responses of LWD resistivity imaging tool in a fractured formation with different dip angles on the basis of 3D FEM. Moreover, we also proposed a new model to calculate the fracture aperture via the maximum current contrast at the fracture and background formation. In the new model, the fractures' dip angles were divided into three parts and the fracture aperture was very easy to obtain. In addition, the new model's accuracy was verified by two single fracture models with different dip angles and apertures.

The principle of measurement
There are two types of LWD resistivity imaging tool excitation source: (i) directly loading current to the electrode and (ii) generating an equal amount of potential on the drill collar by the toroidal coils to achieve the auto-focusing results. Because the latter method is easy to achieve with LWD, it is widely used throughout industry (Gianzero et al. 1985;Kang et al. 2018). Figure 1a shows a schematic of the toroid. If the coil's size is ignored, it can serve as a perfect magnetic ring (figure 1b). Moreover, the measurement frequency cannot be considered because it is usually <1 kHz, so it can serve as an extended voltage dipole (figure 1c). Figure 1 parts b and c depict two different mechanisms; figure 1b is based on electromagnetic field theory while figure 1c is based on the Laplace equation. The two mechanisms are called the AC solution and DC solution, respectively (Bittar & Hu 2004). Figure 1c shows the current emitted in the formation by toroid coils when the measurement frequency effect is ignored, which is similar to that of the common lateral logging method. The volt field is satisfied with the Laplace equation. By solving the Laplace equation, the potential = b from the center of the borehole is mathematically expressed as: (1) where V and are constants, I 0 and K 0 are the zero-order first-and second-class Bessel functions, respectively, and z presents the longitudinal coordinates.
In equation (1), the Γ dc ( ) expression is calculated as follows: where R t is the formation resistivity, R m is the mud resistivity and I 1 and K 1 are first-order first-and second-class Bessel functions, respectively. The current collected between two receiving electrodes at different positions, z 1 and z 2 , is mathematically expressed as: Finally, the apparent formation resistivity R a is obtained using Ohm's law, where

Model
LWD resistivity imaging tool responses were simulated in single and network models based on the 3D FEM. Figure 2 is the single fracture model, which consists of a borehole, LWD tool, fracture and formation. The fracture is a natural and infinite extension filled with conductive mud. By moving the tool vertically and rotating it circumferentially, we obtained the tool responses in single point and full borehole images.

Solution method
The structure of toroidal coils is complex and requires simplification to incorporate it in the numerical simulation. In our previous studies, it was simplified by considering it an extended voltage dipole, as shown in figure 1c. The same method was also adopted in this work. Based on 3D FEM, the LWD resistivity imaging response simulation solution method is similar to that of the lateral logging, namely in that it solves the Laplace equation of current density. A fixed positive voltage and the equivalent negative voltage were applied to the drill collar both under and above the transmitting toroid. In the cylindrical coordinate system (r, , z), the current density satisfies the following partial differential equation: The current density J at the button electrodes can be collected by transforming the definite solution of partial differential equation (5) into the functional value. Then, the current is also calculated by integrating the area of current density, as shown in equation (6): By collecting the button electrode's surface current and bringing it into equation (4), the current can be converted into the apparent resistivity of the formation.

Mesh
There is a large contrast between the button, fracture, formation and tool size in the model. If the traditional structured grid generation method is used, the grid freedom will be too large or calculation accuracy will be insufficient. Therefore, we used the local mesh refinement technology of finite element. The main area of grid refinement includes button electrode and fracture, which could greatly improve the current acquisition accuracy. Also, the grid numbers of the whole model gradually decrease from inside the model to its boundary.

Aperture of the fracture
The fracture aperture is defined as the width of the open part of a fracture, measured at a right angle to the fracture walls (Ameen et al. 2012). The resistivity contrast between the drilling muds and matrix is used to assess fracture aperture (Zeng et al. 2013). Due to the fracture's presence, the measurement current is abnormal at the fracture. Luthi & Souhaité (1990) introduced the concept of abnormal current A's area integral at the fracture, and established the aperture interpretation model using A (equation (7)): where W is the fracture aperture, A is the integrated additional current, R m is the mud resistivity, R xo is the formation resistivity, and b and C are the tool dependent coefficients (Luthi & Souhaité 1990). The fracture apertures can be computed from the electrical image logs through equation (7). The model is obtained when the fracture dip angle is 0°. The calculation error increases for fractures with high angles. Simultaneously, the abnormal area integral also needs to be incorporated into the calculation. Besides concept A, the maximum current contrast can be used to analyze the change in the measured current due to the existence of a fracture. The maximum and minimum current on the fracture track and surrounding formation are defined as I max and I min , respectively. The ratio of I max to I min is called the maximum current contrast. A single fracture model with different dip angles was established, and the relationship between the fracture aperture and the maximum current contrast was further quantitatively studied using different apertures. In the model, the fracture's dip angle and aperture ranged from 0°to 85°a nd from 0.1 to 1 mm, respectively, the fracture and mud resistivity were both 1 Ω ⋅ m −1 , and the formation resistivity was 100 Ω ⋅ m −1 . Figure 3 shows the simulation result. Note that for a single fracture with different dip angles, as the fracture aperture increased, the different fracture dip angles showed a strong linear relationship, which can be expressed as I max ∕I min = a + b ⋅ W where a and b are constants. When the aperture was <0.2 mm, the change in the maximum current contrast corresponding to different dip angles was small; but as the aperture progressively increased to >0.2 mm, the effect of the fracture dip angles on the maximum current contrast increased in parallel. The fitting error associated with using one linear formula to fit the aperture and maximum current contrast degree under different dip angles is shown in Table 1. Clearly, the larger the aperture, the larger the fitting error.
According to the values of the dip angles, the fracture can be divided as follows: (i) low-angle or horizontal fractures with dip angles ≤30°; (ii) medium-angle fractures with dip angles ranging from 30°to 60°and (iii) high-angle or vertical fractures with dip angles ≥60° (Lai et al. 2018). To reduce the fitting error, the fracture angle was divided into three sections: 0-30°, 30-60°and 60-85°. Figure 4 shows that the fitting results of different segments is suitable, while Table 2 depicts the relative error between the simulated value and the fitted value of the maximum current contrast under different fracture dip angles. The fitting error can be effectively reduced by adopting the subsection fitting method.

The resistivity ratio of fracture and formation
Although the LWD resistivity imaging tool button electrode was very close to the borehole wall, there was still mud in the wellbore. At the same time, mud invaded the fracture, filling it within the measurement range. Thus, we also established a single fracture model to quantitatively study the influence of mud resistivity and formation resistivity and subsequently compare the results. In the model, mud and fracture resistivity ranged from 0.1 to 20 Ω ⋅ m −1 , formation resistivity was 100 Ω ⋅ m −1 and the fracture aperture was 1 mm. Figure 5 shows that for a single fracture with different dip angles, as the fracture to formation resistivity ratio increases, the I max ∕I min gradually decreases, showing a linear relationship (actually a power exponential relationship) in the double logarithmic coordinate system. The relationship can be 5 Figure 5. The relationship between maximum current ratio and fracture and formation resistivity ratio.
) d , where c and d are constants. At the same time, as the resistivity ratio between the fracture and formation decreases, the maximum current contrast difference under different dip angles increases. The same linear formula is used to fit the resistivity ratio and maximum current contrast between the fracture and the formation. The fitting error is shown in Table 3. Note that as the fracture dip angle increases, the fitting error increases in parallel and the maximum relative error is >100%.
The relationship between the maximum current contrast and the fracture to formation resistivity ratio is again fitted in sections: as seen figure 6, when compared with figure 5, which demonstrates that the fitting effect is optimized using the sections method. Table 4 lists the relative errors between simulated values of different fracture dip angles and fitted values in the different sections. Except for individual data points, the overall error is far less than that of the data listed in Table 3. Furthermore, error analysis shows that the fitting effect is better when the fracture dip angles are smaller and R m ∕R b is larger.

Calculation model
The existing fracture aperture interpretation model is based on wireline logging, such as FMS. The application results show that a 0.1 mm fracture aperture on the image may be much higher than the actual aperture, so it is impossible to detect the fracture aperture directly with the electrical images' measuring electrode (Zhang et al. 2017). Thus, it is necessary to use the forward modeling method, with the concept of maximum current contrast, to establish the aperture calculation model. Based on these simulation results, the single fracture model with certain dip angles satisfies: where a, b, c and d are constants. When y = I max ∕I min , equation (8) becomes: where p = − a b , q = 1 b . The I max and I min were obtained by the resistivity images, so there is no need to calculate the parameters c and d. The quantitative fracture aperture evaluation formula, under any R m ∕R b , can be obtained using the following equations: The constants of p and q in equation (10) were obtained by 3D FEM simulation on the basis of the specific tool structure. The detailed process is to establish the relationship between fracture I max /I min and aperture under different fracture dip angles and R m /R b to obtain p and q, then establish the relationship between I max /I min and R m /R b under different fracture dip angles and fracture apertures to obtain c and d, and   finally establish a database for convenient access. When the data point is not in the database, the interpolation method is used to obtain the relevant data. Establishing a database may take time, so this work should be done before calculating the fracture aperture. In this way, it is useful for fast and accurate interpretation of the fracture aperture. According to the results, these constants vary under different fracture dip angles. In the actual calculation, the fracture dip angle changes significantly, so it is not practical to use the forward simulation to determine the parameters under all the dip angles. Based on the conclusion drawn in section 4, the fracture dip angle can be divided into three sections 0°-30°, 30°-60°and 60°-85°, which simplifies the calculation and increases the operability in the actual treatment. However, it is imperative that the dip angle is accurately calculated.

The effect factors of calculating the dip angles of fracture
The occurrence information, such as fracture tendency and dip angle, can be quantitatively calculated using the LWD resistivity imaging tool. In the borehole wall image, the fracture unfoldment image is a complete sine curve, and the twodimensional image is in different directions from left to right. The fracture tendency can be determined by identifying the fracture image's wave peak and wave trough. The fracture's dip angle can be calculated using the curve's amplitude (the difference between the wave crest and wave trough longitudinal depth) and the borehole diameter. The sine curve's amplitude is defined as h, the borehole diameter is d and the fracture dip angle is Theoretically, as long as there are enough sampling points, the button electrode resolution is high enough, the fracture image quality is high, and the fracture occurrence calculated by the image should be consistent with the real fracture occurrence. However, in the actual measurement process, there are errors between the calculated fracture occurrence and the real occurrence due to the constraints of data transmission technology, tool measurement accuracy and drilling cost. Thus, detailed error analysis needs to be performed by numerical simulation. In this work, two sets of buttons (B1, B2) with different imaging resolution and different sampling rates (the acquisition bin was reduced from 128 bins to 16 bins) were investigated. Five network fracture images were simulated at fracture dip angles of 15°, 30°, 45°, 60°and 75°, respectively. Among them, three fractures tend to be in the west and two in the east. The fracture resistivity was 1 Ω ⋅ m −1 , the formation resistivity was 100 Ω ⋅ m −1 and the fracture aperture was 1 mm. The B1 and B2 button electrode diameters were 10 and 25 mm, respectively. The directions of 0°, 90°, 180°, 270°and 360°in the image correspond with north, east, south, west and north, respectively. A linear scale was used on the image. The fracture resistivity in the image was low and is marked as dark, while the formation is marked as bright.
The network fracture images of bins 128, 64, 32 and 16, respectively, were drawn using the simulation data. Figure 7 shows that with the bin reduction, the fractures' sinusoidal shapes become more and more fuzzy. The resolution gradually decreases, the entire curve gradually zigzags and a stepped image finally develops. At the same time, the bin reduction leads to the identification of the fractures' sinusoidal shapes when processing fracture images. More difficulty is incurred with curve tracing, which greatly influences fracture dip inversion. Because resolution of the B1 is higher than that of the B2, B2 is more susceptible to being influenced by the sampling rate than B1.
The fracture curve's amplitude h was obtained by inverting the resistivity image, and the fracture dip angle was quantitatively calculated by bringing it into equation (11). But this method exhibits a large error. In practical processing, researchers often select the fracture image characteristic points, use the least square method to fit the sinusoidal fracture equation and then get the highest and lowest fracture image peaks, which is the method used in this work (Chen & Chui 1991;Qiao 2005). Table 5 shows the fracture dip angle obtained from fracture image inversion of bins 128, 64, 32 and 16. Note that with gradual reduction of the bins, the calculated fracture curve amplitude accuracy gradually decreases. For low-angle fractures (the fracture dip angle is <45°), the error between the inverted dip angle value and the true dip angle becomes progressively larger. For high-angle fractures (fracture dip angle >45°), the bin reduction has less impact. At the same time, it is noted that when imaging with the same bin, image quality of B1 has higher resolution and is therefore better than that of B2, so the dip angle inversion accuracy is higher in the former.

Model verification in theory
To verify the calculation accuracy of fracture aperture using the new model, two fractures with different dip angles and apertures were simulated using the single fracture model. Most of the fracture dip angles in the actual formation were medium or high-angle fractures, so the model parameters of dip angle were set to 45°and 60°. As shown in figure 8, the  dip angle in figure 8a is 45°and the aperture is 1 mm, while the dip angle in figure 8b is 60°and the aperture is 0.5 mm. The mud fracture resistivity is 1 Ω ⋅ m −1 and that of the formation is 100 Ω ⋅ m −1 . B1 was used as the button electrode and had a 10 mm diameter. In the image, the directions of 0°, 90°, 180°, 270°and 360°correspond to north, east, south, west and north, respectively. A linear scale was also used in this image. First, the two fractures' dip angles were inversed by the least square method. The inversion results were 45.7°and 59.8°, respectively, which are slightly different from the fractures' true dip angles. The inverted fracture dip angle was brought into equation (10) to calculate the fracture aperture. In figure 8 parts a and b, the calculation model adopted medium and high-angle fractures, respectively, which resulted in calculated apertures of 0.95 and 0.51 mm, respectively. The relative errors in the model with the real aperture were 5% and 2%, respectively. Thus, the solution process is rapid and the model accuracy was theoretically verified.

Conclusions and discussion
In this work, we investigated the logging response characteristics of the LWD resistivity imaging tool, based on the 3D FEM, in a fractured formation with different dip angles. This investigation simulated both the fracture center point logging response characteristics and the imaging results. Processing the simulated fracture image helped to determine the fracture dip angle influencing factors and verify the model accuracy. Compared with physical simulation, the interference from other uncertain factors is well controlled.
The results show that there is a linear relationship between the maximum current contrast and the fracture aperture, and that the fitting coefficients differ with different fracture dip angles. The maximum current contrast has a superior power exponential relationship with the fracture and formation resistivity contrast, but the fitting error is large just using one formula when the fracture dip angles vary. Through error analysis, we determined that the fitting relative error can be effectively reduced by fitting the above relationship in sections. Based on the simulation results, a fracture aperture calculation model was established. Compared with the commonly used model, it is suitable for fracture aperture calculation under different dip angles. Meanwhile, only the formation maximum and minimum fracture current need to be detected. Thus, the area integral of abnormal current around the fracture is not needed, making the model parameters easy to obtain. In addition, the model first needs to calculate the fracture dip angles. Through analyzing the factors affecting the dip angles, we determined that increasing the tool's number of acquisition bins corresponds to a higher imaging resolution and a more accurate fracture angle inversion. Therefore, it is very important for improving accuracy of the measurement results. Finally, the accuracy of the new model was verified with a theoretical model. The fracture aperture calculation model obtained in this work provides a theoretical basis for logging interpretation of the LWD resistivity imaging tool. Moreover, it should be noted that the LWD tool cannot always be in the borehole center. In the case of eccentricity, it will be useful to further verify the new model proposed in this paper. Because the cable type resistivity imaging tool measurement principle is similar to that of the LWD imaging tool, the model is also suitable for cable imaging instruments, such as FMI.
However, there are often many intersecting network fractures in the actual formation, and the applicability of a single fracture model needs to be evaluated as well. Although in this work we tried to use a more reasonable method for finite element mesh generation, fracture image simulation is still time consuming. Thus, it is necessary to propose a new mesh generation method to improve the calculation efficiency of the forward model in the future.