New Horizons: Novel Precise TVD Assurance Workflow

Konstantin Bulychenkov

Konstantin Bulychenkov, Victor Ploskikh

Abstract

High-fidelity geological modeling, primarily vertical bedding, is critical in the economics of development drilling. The precise measurement of the true vertical depth (TVD) is an essential part of building such geological models. However, the TVD is subject to various errors. In deviated wells, particularly in extended reach wells, the main sources of TVD errors are the bottom hole assembly (BHA) sag, misalignment, Stockhausen effect [11], and axial accelerometer error. The oil and gas industry employs several approaches to compensate for these errors:

  • The BHA sag is compensated for using algorithms based on the Lubinski general solution [1] or finite element methods (FEM);
  • The Stockhausen effect is addressed by fusing static surveys and continuous inclination measurements or fusing static surveys with a slide sheet [9];
  • The misalignment remains unaddressed;
  • Accelerometer errors are compensated for using in-lab tool calibration.

However, existing approaches to TVD correction have many problems, and their development lags behind the ever-tightening geosteering requirements. The accuracy analysis by Williamson [2] revealed a considerable residual error of the BHA sag correction, and a study by Codling [3] showed a significant discrepancy in the inclination corrected for the BHA sag from different tools in the same well. Bulychenkov [10] showed that the use of the fusion algorithm of static measurements with the slide sheet can result in a larger TVD error than that of the uncorrected trajectory.

The aim of this study was to develop a comprehensive TVD assurance workflow that addresses the highlighted problems and, as a result, improves vertical borehole positioning.

To achieve this the following work has been done:

  1. The existing approaches were analyzed for potential problems: The limitations of the BHA sag correction are assumptions on the initial constraints [1] and handling the collar-wellbore points of contacts; Stockhausen effect compensation struggles with problems of continuous inclination data density and instability of the algorithm when utilizing slide sheet information [4]; despite the claimed accelerometer error being insignificant, the wide boundaries for the total gravity field acceptance criteria indicate that the actual accelerometer error is considerable and must be considered.
  2. The novel multi-station analysis (MSA) proposed by Bulychenkov [15] was utilized for axial accelerometer correction in the horizontal section. The method was validated using simulated data, and it exhibited a significant improvement in accuracy.
  3. A novel 2D generic deflection calculation method based on an energy minimization approach [5] was developed. The proposed method overcomes the typical problems of conventional approaches: it does not require initial constraints for the column–beam, calculates the contact points automatically, handles continuous contact naturally, and provides better solution stability than commercial software based on a finite element (FE) method.
  4. The proposed method exhibited high accuracy in comparative tests:
    • The algorithm was compared with the Lubinski method in a synthetic environment using a simple BHA;
    • The algorithm was verified for a highly deviated well drilled in a laboratory and measured using a laser probe [6];
    • The algorithm prediction was compared with dual inclination measurements acquired while drilling.
  5. A theoretical error propagation analysis for the proposed sag correction method was performed, and the results were similar to those of [7]. In addition, dog leg severity (DLS) deviation factors were explored, and we observed that the major error source for the BHA sag correction is the wellbore local curvature, which conforms to [3].
  6. To compensate for this, we developed a novel TVD assurance workflow. The workflow encapsulates the novel MSA, proposed BHA sag correction algorithm, high-definition (HD) trajectory calculation, and a novel method for wellbore micro-profile approximation. The true trajectory is restored through sequentially iterating the fusion of the HD trajectory with the slide sheet and the subsequent BHA sag correction of the HD trajectory.
  7. The workflow was tested on two typical drilling cases: a build-up section and a tangent section. In both cases, the workflow exhibited rapid convergence and accurate compensation for the BHA sag error, BHA misalignment and Stockhausen effect. The proposed workflow demonstrated TVD accuracy improvement of 7–15 times compared with the conventional TVD correction workflow.

In addition, because the proposed TVD assurance workflow requires processing data from multiple sources, it was built into a cloud-based automatic survey management system, MWD STD Advanced.

Introduction

Well productivity directly depends on the quality of geological models used for geosteering and reservoir simulations. Consequently, the quality of such models is largely determined by the TVD measurement accuracy. In deviated wells, the major part of the TVD uncertainty is caused by inclination measurement errors and the Stockhausen effect [11].

In MWD tools, inclination is measured using three-axis accelerometers that are subject to various sensor errors, such as bias, scale factor, misalignment, and nonlinearity. These errors affect the measurements and, consequently, contribute to the TVD error, particularly at high inclinations. However, these errors are assumed to be effectively compensated for through the in-lab calibration of the directional sensors and that they are of an order of magnitude less than other inclination error sources.

Owing to gravity, mechanical properties of BHA, and wellbore conditions, both MWD and gyroscopic inclination measurements are subject to the BHA sag and misalignment errors. The error magnitudes of the uncorrected BHA sag and misalignment are 0.2° and 0.1°, respectively, both of which are systematic [12]. In the horizontal section, these errors cause a TVD error (2 sigma) of 7.8 m per 1000 m of drilled depth. The industry employs two approaches to compensate for the BHA sag: the Lubinski solution and FE methods. The Lubinski method is based on an analytical solution of stress equations. FE methods are discrete solutions based on small deformation theory. However, both approaches are used to compensate for the BHA sag only. The misalignment error is de facto ignored.

Williamson’s analysis [2] revealed large residual errors in the sag correction methods (1 sigma of 0.08°). Codling [3] showed a significant discrepancy of up to 0.5° between the sag-corrected inclination of the MWD and gyro surveys. This value cannot be caused by an accelerometer error, which has a magnitude of 0.02°. These facts may be explained by the limitations of the currently employed methods, particularly the necessity of specifying points of contact.

The Stockhausen effect is an interpolation TVD error caused by a low survey frequency. Trajectory interpolation methods assume a constant curvature between consecutive surveys. In practice, a sequence of curved and straight drilling intervals is typically employed to achieve the target DLS on a trajectory segment. Consequently, this can result in a significant TVD error (more than 10 m) and incorrect interpretation of azimuthal logging-while-drilling measurements [11].

The correction of the Stockhausen effect is based on the creation of artificial surveys to increase frequency. Currently, two approaches are used: fusion of static and continuous inclination surveys and utilization of a slide sheet [9]. However, both approaches have limitations. A fusion highly depends on the uniformity of the continuous inclination data, which is unachievable in real-time applications. Interpolation using slide sheet data uses very strict assumptions, which makes the method unstable in some cases [10].

In this study, we analyzed existing approaches to high-fidelity TVD assurance. The main problems with the existing methods have been identified. To solve these problems, we developed novel algorithms for BHA deflection modeling, HD trajectory calculation, and wellbore micro-DLS approximation. The proposed method of calculating the BHA deflection overcomes the typical problems of the conventional approaches, as it does not require initial constraints for the column–beam and does not have limitations on BHA design; it automatically calculates the contact points and naturally handles the continuous contact. The wellbore micro-DLS approximation fuses the HD trajectory with the slide-sheet information to restore the true wellbore micro-profile. The proposed algorithms were integrated into a comprehensive high-fidelity TVD assurance workflow. The workflow was tested using typical drilling cases: build-up and tangent sections. As a result, a significant improvement in TVD accuracy and correction stability was demonstrated: the residual TVD error was reduced by a factor of 7 to 15 compared with the conventional TVD correction approach. The TVD assurance workflow was integrated into an automatic cloud-based survey management system to provide real-time services with minimal labor costs.

Current Approach Analysis

TVD correction algorithms for deviated wells (excluding depth correction) can be divided into four groups: BHA sag correction, misalignment correction, Stockhausen effect compensation, and accelerometer error.

BHA Sag Correction

According to [4], the main challenges of accurate BHA modeling can be defined as follows:

  • System of nonlinear differential equations
  • Unknown boundary conditions at the bit, stabilizers, and top of the BHA
  • Unknown contact points between the BHA and borehole walls
  • Considering the bend in the motor or rotary steerable system (RSS)

FE methods based on the small deformation theory are the most widely used because of their relative simplicity. However, the large deformation caused by the bend in the BHA (at the motor or RSS) results in low accuracy in BHA modeling [8]. In addition, FE methods do not provide sufficient accuracy and stability in determining the contact points between the BHA and borehole walls, which often leads to erratic results [4].

Modern solutions based on the Lubinski approach [1] and the “self-structured” BHA model overcome some of the disadvantages of FE methods but also have some limitations [4]:

  • Assumption of boundary conditions at the top of the BHA
  • Assumption of the bit centering in the wellbore
  • BHA design limitation: at least a stabilizer is present

Thus, we can state that despite significant progress in BHA modeling since the publication of the Lubinski approach, many unsolved problems remain, which hinders the improvement of TVD correction methods. The development of a novel BHA modeling method that solves the aforementioned disadvantages is a logical step in this direction.

Misalignment Correction

The inclination error due to misalignment between the direction and inclination (D&I) probe and collar is considered random and averaged to zero because its effect is toolface-dependent, and survey toolfaces are uniformly scattered during the drilling run [14].

The misalignment between the MWD tool and wellbore is considered a systematic error [12]. The authors did not find any mention of a correction for this effect.

Stockhausen Effect Compensation

Two generic approaches are available: the fusion of static surveys with continuous inclination and the fusion of static surveys with a slide sheet.

The fusion of static surveys with continuous inclination provides satisfactory accuracy in compensating for the Stockhausen effect if the density of the continuous inclination data points is high. Otherwise, if the density is insufficient, it is impossible to effectively compensate for the measurement errors of the continuous inclination. This makes the accurate calculation of the HD trajectory impossible in real time. Additionally, the MWD tool may not have the option of recording continuous inclination data in memory.

The fusion of the static surveys with the slide sheet does not have a data density problem; however, the analysis performed in [10] demonstrated the significant potential instability of the method under certain conditions: owing to BHA misalignment, the fusion algorithm incorrectly evaluates the rotor build-rate and motor yield. Although these errors are small, because they are systematic, they cause TVD error accumulation along the entire trajectory. As a result, the error of the corrected TVD may be several times greater than that of the TVD of the uncorrected trajectory.

Accelerometer Error

The accelerometer errors according to the model [12] are assumed to be 0.004 m/s2 and 500 ppm for the bias and scale factor, respectively. In the horizontal section, the axial accelerometer error results in an inclination error of 0.02° per 1 sigma. In terms of the overall error budget for inclination, this error is insignificant and ignored.

However, an analysis of the field acceptance criteria (FAC) for the total gravity field demonstrated that the actual accelerometer errors are higher. For the worst-case scenario, the gravity FAC calculated using the MWD error model must be approximately 1.25 mG for the 2 sigma error level accepted by the oil and gas industry [16]. The gravity reference error calculated by the latitude gravity model is negligible ( most errors due to gravity anomalies are within +/-0.1 mG) in terms of sensor accuracy. However, oil and gas service companies typically use FAC for gravity between 2.5 and 3.0 mG. This means that, on average, the error of the accelerometer is 2.2 times higher than the budget of the error model. In terms of inclination accuracy, the 1 sigma error is 0.045°, which is comparable to the residual uncertainty of the BHA sag correction (0.08°) and must be addressed. However, multi-station analysis (MSA) cannot calculate the bias and scale factor of the axial accelerometer in the horizontal section because of the fundamental limitations of the algorithm, similar to the MSA no-go zone for azimuth correction [15].

Intermediate Conclusion

Despite the many existing approaches to TVD improvement, the BHA sag and Stockhausen effect corrections have significant limitations. Moreover, accelerometer errors in the horizontal section and misalignment remain unaddressed because of the fundamental limitations of MSA correction and the limitations of BHA deflection models, respectively. Thus, the assurance of precise TVD while drilling requires resolving the aforementioned problems and constructing a reliable workflow orchestrating the TVD correction sub-algorithms.

Accelerometer Correction

During the drilling of a horizontal section, the axial accelerometer bias contributes primarily to the TVD error. However, no generic procedure exists to eliminate this type of error: the MSA algorithm cannot compensate for the error in the horizontal section because of the fundamental limitations of the method. This error cannot be detected because the field acceptance criteria are useless in this case. Nonetheless, the MSA algorithm can be modified to overcome this fundamental limitation by utilizing additional information. Please refer to [17] for further details. Hereinafter, the algorithm is denoted as novel_msa(…).

The novel MSA method can effectively compensate for the axial accelerometer bias in the horizontal section (Figure 1). The simulation tests show that a 2 mG axial accelerometer bias (that is approximately 2 sigma of the revised error budget for the bias) can cause an inclination error up to 0.12° in the horizontal borehole that translates to a TVD error of 2.1 m per 1000 m drilled (red line in Figure 1). The conventional MSA correction fixes the total gravity field reading; however, it is ineffective for the axial bias correction (dashed line). The novel MSA algorithm reduces the residual inclination error to 0.015° (green line).

Figure 1. Axial accelerometer bias, which cannot be detected using FAC and corrected using conventional MSA, causes a TVD error. Using additional data and a modified MSA algorithm can significantly reduce this error

Novel BHA Modeling

The foundation of the high-fidelity TVD assurance workflow should be a BHA deflection algorithm that overcomes the common limitations listed above.

Algorithm

The algorithm is an FE method based on the minimization of the drillstring energy:

where xopt is an optimal BHA deflection corresponding to the minimum BHA energy.

The optimal deflection is calculated using a gradient descent method: the BHA is divided into uniform tubular elements, and then the system of the elements sequentially evolves from the initial position to the optimal one in the direction of the reverse energy gradient. The contact points of the BHA are automatically calculated during the evolution. This eliminates the necessity for initial constraints such as having at least one stabilizer in the BHA or the presumption of fixed contact points (e.g., centered bit position). Consequently, the minimum number of assumptions increases the accuracy and stability of the solution, particularly in a highly deviated borehole. The algorithm can naturally consider BHA bending by specifying an offset curve for the zero energy of BHA deformation. The algorithm is based on only three general assumptions:

  • The model is 2D
  • The model is static
  • All deformations are elastic

Therefore, the proposed algorithm overcomes all the limitations of the Lubinski approach and the common problems of FE methods.

The following are the key features of the algorithm:

  • Calculation of contact points between the BHA and borehole walls in a natural manner
  • No special restrictions on BHA design (e.g., the presence of a stabilizer)
  • No assumptions about boundary conditions (e.g., fixed bit position)
  • Consideration of the bend angle of the steerable tool
  • High solution stability

A full description of the algorithm is given in Appendix A.

Hereinafter, the algorithm is denoted as sag_correction(…).

The algorithm exhibits increased stability of sag estimation (Figure 2).

Figure 2. Comparison between generic and proposed sag correction methods: A – deflection modeling of the BHA; B – trajectory profile; C – calculated sag: generic FE method – red, proposed FE method – green.

Algorithm Validation

Several tests were conducted to ensure algorithm validity.

Algorithm vs Lubinski solution

The first test compared the algorithm with a theoretical solution to the steel pipe deformation problem. We considered a 5 m long steel pipe (outer diameter: 120 mm, inner diameter: 40 mm), which was fixed at one end at a height of 80 mm, and the other end was in contact with the surface. The calculated deflections for both methods were in good agreement, as shown in Figure 3.

Figure 3. Comparison between the closed form solution and proposed FE method

Test on Lab Data

The second test of the sag correction method was conducted on a case described in [6]: A deviated wellbore was drilled through concrete blocks in a laboratory and then measured using an MWD tool for inclination and a laser probe for absolute displacement. An image of the vertical displacements [6, Fig. 16] as the wellbore boundaries was captured. The borehole boundaries are indicated by a black line in Figure 4. The MWD measurements were extracted from the image [6, Fig. 18]. The BHA was composed of an MWD (9.6 m long, 171 mm in diameter) fitted with a bullnose (0.1 m long, 151 mm in diameter) at the bottom and a crossover to the drill pipes at the top; the total BHA length was 10 m with a D&I distance to the bottom of 4.1 m.

Figure 4. Vertical profile of the test trajectory.

The MWD inclination reading was modeled using the proposed correction algorithm (Appendix A), BHA, and actual borehole boundaries. Figure 5 shows a comparison of the modeled inclination (green line) against the inclination measured using the MWD tool (red line): The modeled data reproduced the measured data with good accuracy through the entire drilled interval; The average difference between the calculated and measured inclinations was 0.015°. The blue dashed line represents the smoothed inclination of the borehole centerline.

Figure 5. Measured vs modeled MWD inclinations

Dual Inclination Test

The third test considered the data from two cases. In both cases, a second stream of inclination measurements from an RSS was available in addition to static MWD surveys. The simplified BHA designs are presented in Table 1. Ideally, the observed difference in inclination from different sensors at the same depth should be equal to the difference in the calculated sags. In practice, random errors still occur in the data; therefore, we compared the mean differences.

Table 1. List of the tested BHAs

# Dia. BHA
1 12¼” BIT-ST.RSS.ST-LWD-MWD-SST-XO-NMDC-NMDC-SUB-XO-HWDP-DC
2 6” BIT-RSS.ST-MWD-SST-LWD-XO-SST-LWD-SST-XO-NMDC-XO-NMDC

In the first case, the observed average inclination difference was 0.505° and the average difference between the computed sags was 0.522°. In the second case, the observed average inclination difference was -0.386°, and the average difference of the computed sags was -0.349°. The average difference between the simulated and measured values, normalized to the standard deviation of the accelerometer sensor error (sqrt(2)*0.023°), was 0.94 sigma, which indicated good agreement between the measured and simulated values.

Error Propagation Analysis

A BHA sag error propagation analysis was conducted using the proposed BHA deflection algorithm. The analysis was qualitative. Its objective was to identify the major error sources and not to obtain exact numerical estimates. The main workflow was based on [7]. Four typical BHAs were considered: two with a motor and two with an RSS, with hole sizes of 6, 8, 12, and 16 inches. The BHA compositions are listed in Table 2.

Table 2. List of the BHAs used for analysis

Dia. BHA
6” BIT-RSS.ST-MWD-SST-LWD-XO-SST-LWD-SST-XO-NMDC-XO-NMDC
8⅝“ BIT-ST.DHM-FS-SST-NMDC-MWD-SST-XO-NMDC-DC
12¼” BIT-ST.RSS.ST-LWD-MWD-SST-XO-NMDC-NMDC-SUB-XO-HWDP-DC
16” BIT-ST.DHM-SST-LWD-MWD-LWD-XO-SST-DC-SUB-XO-DC-DC

Several parameters were tested to assess their influence on BHA deflection. The parameters are listed in Table 3. They generally replicated the modeling parameters of [7] with the only difference being that error terms with relatively small effects were omitted, and an additional parameter of local curvature was added. The error modeling mechanisms are listed in Table 3. N(m, s) indicates a normal distribution around m with a standard deviation of s. U(a, b) indicates a uniform distribution in the [a, b] interval. The distributions were selected such that they represented physical conditions: the stabilizer OD can only decrease, the toolface can have any value because of the inability to deduce it from surveys, etc.

Table 3. Error modeling parameters

Name Distribution Prop. Mode Max Error, deg
Borehole overgauge Bit OD * (1 + U(0, 0.05)) Systematic 0.034
D&I position N(Position, 0.15 m) Systematic 0.034
Mud weight N(Mud weight, 100 g/l) Systematic 0.005
Motor bend orientation (toolface) U(-180 deg, 180 deg) Random 0.053
Stabilizer OD Stab OD – U(0, 0.0075) Systematic 0.054
Motor bend angle N(Bend angle, 0.1 deg) Systematic 0.027
Local curvature U(0 deg/10m, 1 deg/10m) Systematic 0.29

The local curvature was modeled as follows: A typical borehole profile with an intensity of 3° per 30 m was considered. A uniform curvature was selected as the baseline. Subsequently, a typical rotor-slide-rotor pattern was simulated by redistributing the intensity between intervals (a decrease for the rotor and an increase for the slide). In the extreme case, the DLS on the rotor intervals became zero, and the entire intensity was on the slide interval.

The results of the error propagation analysis are shown in Figure 6. Because the analysis was qualitative, for illustration, the maximum sag deviations averaged over the BHAs are presented. As shown in the graph, the local curvature error contribution was almost five times greater than that of the runner-up.

Figure 6. Maximum sag error averaged between different BHAs

This large error was in good agreement with [3], which demonstrated a significant difference (up to 0.5°) between MWD and gyro surveys, and the magnitude of this difference correlated with the borehole dogleg. This is because the deflection profile is heavily dependent on the local curvature of the borehole (Figure 7).

Figure 7. Two deflection profiles of the same BHA in the boreholes with the same average DLS: green – uniform curvature, red – rotor-slide-rotor sequence. A) Centered view; B) absolute displacement.

The BHA sag correction using only static surveys cannot address this effect. Therefore, a novel method that considers local curvature is required.

TVD Assurance Workflow

To compensate for the BHA sag and misalignment effectively, we must know the true trajectory; however, to obtain the true trajectory, we must compensate for the BHA sag and misalignment. Here is a paradox: to obtain the true trajectory, we must know the true trajectory. This is because the measurement process in a stiff collar is a type of convolution process. A convolution example is shown in Figure 8, where the true inclination (blue line) is smoothed by the stiff BHA (green line) to the measured inclination (red line), and the true inclination is lost.

Figure 8. Convolution example

To restore the true trajectory profile, we must deconvolve the HD trajectory with respect to the BHA design. However, deconvolution is unstable, sensitive to noise, and requires a high density of continuous measurements, which is not always achievable, particularly in real time. To overcome these problems, we utilize information from the slide sheet: the slide/steering sheet data are fused with the HD trajectory derived from the static surveys and continuous inclination. Information from the slide/steering sheet enables us to completely stabilize the algorithm, minimize the effect of noise in the input data, and solve the problem of a low continuous inclination density.

Sub-Algorithm: Trajectory Approximation

An important element of the deconvolution or true-trajectory restoring algorithm is the piecewise linear approximation of the HD trajectory with the slide sheet information; the algorithm calculates the values and their uncertainties for the BHA DLS performance in the rotor and slide modes and approximates the HD trajectory piecewise linearly to maximize confidence (the optimal fitting of the slide sheet model to the HD trajectory). This algorithm is described in detail in Appendix D.

Ref. [6] demonstrated that the inclination changes nonlinearly: a noticeable step increase occurs in inclination at the beginning of the slide and a similar step decrease at the end. However, our analysis shows that laboratory conditions differ significantly from downhole conditions.

First, in the laboratory test, the borehole diameter was significantly larger than the diameter of the bit in the rotor section (242 mm vs. 216 mm). According to the caliper measurements, such a large overgauge is not observed when drilling a conventional well. This is caused by the significantly higher stiffness of the rock formation compared with that of concrete. This results in a significant decrease in the motor bend angle owing to deformation (see Figure 9). As the slide progresses, the deformation recedes, and the slide intensity drifts to a nominal value (calculated using the three-point model).

Figure 9. In the ingauged borehole, the bit attack angle is almost zero (green), and in the overgauged borehole, this angle can be higher the motor bend angle (red).

Second, when drilling an wellbore with a depth greater than 500 m, it is impossible to begin drilling a slide with the desired toolface. Because the reactive torque on the bit can twist the drill string by 180 °or more, the initial toolface of a slide is frequently very different from the desired one. Moreover, during the drilling process, the toolface changes smoothly until the torque on the bit is balanced by the torque of the drill string (Figure 10).

Figure 10. The motor toolface is changed by approximately 300° until it is set to the desired value.

Thus, in a downhole environment, there is no basis for a sharp change in the inclination of the wellbore at the beginning of the slide. Based on these considerations, we have selected a linear fit as the simplest and most accurate method in the absence of complete information; approximation errors will have a zero mean.

As shown in Figure 11, this type of approximation is rather crude. Moreover, it overestimates the rotor build rate. This is caused by the pick-up effect [6] when the MWD perceives the slide ahead owing to the stiffness of the collar. However, this algorithm still enables us to extract information from the slide sheet regarding the true trajectory profile, albeit with an error. It also minimizes the effect of a low continuous inclination density. Using it in conjunction with the BHA sag correction algorithm in an iterative process minimizes this error.

Figure 11. Trajectory approximation: blue line – based on the low density continuous inclination with the slide sheet; green dotted line – based on the high density continuous inclination

Hereinafter, this algorithm is denoted as aprx_traj(…). The algorithm uses the output of the HD trajectory algorithm as its input.

Sub-Algorithm: HD Trajectory

Parts of the HD trajectory algorithm are presented in [13]. They resolved the most important problems in HD trajectory calculation:

  1. Noise in continuous inclination surveys owing to drillstring vibration, poor demodulation, and telemetry signal delay
  2. DLS overestimation owing to arbitrary data density
  3. Low accuracy of continuous inclination measurement owing to gravity reference and directional sensor calibration errors
  4. Determining the optimal representation of the HD trajectory.

The main part of the algorithm is continuous-static fusion. The algorithm is described in detail in Appendix C and is denoted as hd_traj(…).

Figure 12 shows the operation of the algorithm for the case in [11]. The section was drilled using a motor, and static and continuous inclination surveys were acquired during drilling. Subsequently, the borehole was surveyed using a wireline gyroscopic tool in the continuous mode. Gyroscopic measurements are consistent with static MWD but provide a much higher survey density, and the gyroscope assembly is more flexible and better repeats the micro-profile of the well. Gyro measurements were selected as a reference to evaluate the performance of the HD algorithm. The TVD calculated using only the continuous inclination had an error of approximately 2 m per 300 m section (the average continuous inclination error was approximately 0.4°). The TVD error in the static measurements was about -1 m owing to the Stockhausen effect. The TVD error of the HD trajectory was approximately 0.15 m.

Figure 12. HD trajectory algorithm builds the optimal representation of the true trajectory (blue), which has a minimal TVD error compared with static surveys (green) and continuous inclination readings (red)

TVD Assurance Workflow

Here, we have all the necessary blocks to build the algorithm that reconstructs the true trajectory from continuous inclination, static surveys, BHA design, and slide sheet:

  1. Correct the static surveys for accelerometer errors: Scor = novel_msa(Sraw, BHA, SS, ref)
  2. Calculate the initial HD trajectory: Shd = hd_traj(Scor, Sci, BHA)
  3. Calculate the HD trajectory corrected for sag: IF first time THEN S1no_sag = sag_correction(Shd, Shd, BHA) ELSE Sjno_sag = sag_correction(Shd, Sj-1aprx, BHA, MW)
  4. Restore the approximated trajectory: Sjaprx = aprx_traj(Sjno_sag, SS, BHA)
  5. Check stop criteria as: norm(Sjno_sag – Sj-1no_sag) < ε
  6. Go to #3

where Sraw is a set of raw static surveys, Scor is the corrected static surveys using the Novel MSA, SS is the slide sheet, ref is the geomagnetic reference, Sci is the continuous inclination measurement, Shd is the original HD trajectory, Sjno_sag is the HD trajectory corrected for BHA sag/misalignment for the jth iteration, Sjaprx is the piecewise linear trajectory approximation for the jth iteration, ε is a small number (stop criteria), and MW is the mud weight.

Figure 13. Flowchart of the TVD assurance algorithm

Figure 13 shows a block scheme of the TVD assurance workflow. The output of the workflow is an approximated trajectory.

The workflow was automated and embedded into MWD STD Advanced, an automatic cloud survey management platform.

Verification of Workflow

To validate the TVD assurance workflow, we conducted a series of experiments using synthetic data. Two typical drilling cases were considered: inclination build-up and inclination hold.

The “inclination build-up” case was modeled as a periodical sequence of rotor and slide intervals with a slightly dropping motor BHA (rotor build rate of -0.75°/30 m and motor yield of 6°/30 m). The true trajectory profile is indicated by the blue line in Figure 14 (top left).

The “inclination hold” case was modeled as a periodical sequence of rotor and slide intervals with a dropping motor BHA (rotor build rate of -1.25°/30 m and motor yield of 5°/30 m. The true trajectory profile is indicated by the blue line in Figure 14 (bottom left).

Static surveys were acquired every 30 m. The continuous inclination points were randomly distributed along the trajectory with an average density of 1 point per 3 m. These parameters were selected to simulate real-time continuous inclination measurements, considering the effects of the penetration rate and telemetry signal quality. The true trajectory was used to calculate the inclination error due to sag and misalignment. The simulated static and continuous inclination surveys with added errors are shown as red dots in Figure 14.

Three options were considered: 1) HD trajectory corrected for BHA sag and misalignment using static surveys; 2) HD trajectory corrected for BHA sag and misalignment using HD surveys; 3) the result of the TVD assurance workflow. The first two options were used as benchmarks. The first option is the standard TVD correction procedure utilized by the industry, and the second option is a more advanced version of the first and is not frequently used because of difficulties in software processing a trajectory with a high density of points. The residual TVD errors of the correction options for both cases are shown in Figure 14.

Figure 14. Plots of the inclination, and TVD error for the “inclination build” case (top row) and “inclination hold” case (bottom row).

As shown in Figure 14, the conventional TVD correction workflow (BHA sag correction and HD trajectory calculation) still had a noticeable residual error in TVD: 1.5–2.5 m per 1000 m drilled. Furthermore, this simulation did not consider the axial accelerometer bias effect discussed in the Accelerometer Correction section, which can additionally create a TVD error of 1–2 m per 1000 m drilled. The proposed workflow can significantly improve TVD accuracy (by 7–15 times), particularly at low continuous inclination densities, which are encountered in real-world conditions with poor telemetry signal quality, high rate of penetration, or low-frequency telemetry.

Automation

Maintaining the proposed workflow requires data from different sources and numerous complex calculations. Because manual execution is problematic, the workflow was automated and integrated into the survey data management cloud system, MWD STD Advanced (Figure 15). The system executes the TVD assurance workflow in real time and provides data consistency control and analysis of the quality of calculations in the automatic mode [13].

Figure 15. Architecture of MWD STD Advanced, a survey data management cloud system

Error Model

The present study did not seek to develop an error model; however, the residual error can be roughly estimated from error propagation analysis. If we take the conservative estimate that the error maximum corresponds to 2 sigma and calculate the root mean square of all the analyzed error terms, then the total residual error is 0.0385°. The contribution of the internal parameters (grid step and convergence value) of the method itself is 0.02°, which agrees with the laboratory data test. Therefore, the magnitude of the SAG error term is 0.043°, and the propagation mode is systematic.

Because XYM1, XYM3, and XYM4 encapsulate MX and MY error sources from [2] (misalignment between tool housing and sensor frame) and these errors are not addressed by the workflow, the residual magnitude of misalignments XYM1, XYM3, and XYM4 is 0.06°, and the propagation mode is random because these errors depend on a randomly distributed toolface.

A covariance analysis by [17] showed that systematic ABZ error terms can be reduced from 0.004 to 0.001 m/s2.

The updated uncertainties are summarized in Table 4.

Table 4. Updated error magnitudes after correction

Error Term Description Prop. Mode Magnitude
SAG Residual BHA sag error Systematic 0.043 deg
XYM1, XYM3, and XYM4 BHA misalignment error Random 0.06 deg
ABZ Axial accelerometer bias Systematic 0.001 m/s2

Conclusion

The proposed TVD assurance workflow orchestrates a comprehensive set of sub-algorithms: novel MSA correction, novel BHA deflection modeling, HD trajectory fusion, and novel trajectory approximation methods. An accurate well profile is the product of the fusion of six-axis static surveys, continuous inclination measurements, slide/steering sheets, and BHA design information. The workflow effectively compensates for the major TVD error sources: D&I accelerometer biases, BHA sag and misalignment, and a low survey frequency. As a result, the accuracy of the TVD correction increases by 7–15 times that of the conventional TVD correction workflow (BHA sag correction + HD trajectory), in which the TVD error may exceed 5 m per 1000 m drilled. Furthermore, the workflow works effectively with the low density of continuous inclination, which is a common occurrence. The proposed TVD assurance workflow was integrated into the cloud automatic survey management system MWD STD Advanced, which enables real-time processing without compromising effective well placement.

This paper also presents in detail a novel 2D BHA deflection modeling method based on FE analysis and energy minimization, which overcomes the main limitations of the BHA modeling methods, such as contact point designation, low solution stability, requiring stabilizers, and strict initial constraints. The method was tested on a laboratory-drilled borehole and showed a high accuracy. In the future, the method can be further elaborated for a more comprehensive BHA and drillstring modeling: tendency, side forces, and fully stiff torque and drag analyses.

References

  1. Lubinski, A., and Woods, H. B., “Factors Affecting the Angle of Inclination and Dog-Legging in Rotary Bore Holes”, API Drilling and Production Practice, 1953, pp222-250
  2. Williamson, H. S. 2000. Accuracy Prediction for Directional Measurement While Drilling. SPE Drill & Compl 15 (4): 221–233. SPE- 67616- PA. doi: https://doi.org/10.2118/ 67616- PA
  3. Codling, J., and Zhang, S. 2022. Survey Data Comparison and Data Science. Paper presented at the 56th General ISCWSA Meeting, Houston, USA, 6 October. https://www.iscwsa.net/media/files/events-event/96eaac81/jerrycodling-iscwsa-56th-presentation.pdf
  4. Wu, M., and D. C. -K. Chen. “A Generic Solution to Bottomhole-Assembly Modeling.” Paper presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, September 2006. doi: https://doi.org/10.2118/101186-MS
  5. Brands, S., and R. Lowdon. “Scaled Tortuosity Index: Quantification of Borehole Undulations in terms of Hole Curvature, Clearance and Pipe Stiffness.” Paper presented at the IADC/SPE Drilling Conference and Exhibition, San Diego, California, USA, March 2012. doi: https://doi.org/10.2118/151274-MS
  6. Stockhausen, Ed, Lowdon, Ross, and Bill Lesso. “Directional Drilling Tests in Concrete Blocks Yield Precise Measurements of Borehole Position and Quality.” Paper presented at the IADC/SPE Drilling Conference and Exhibition, San Diego, California, USA, March 2012. doi: https://doi.org/10.2118/151248-MS
  7. Studer, R., and L. Macresy. “Improved BHA Sag Correction and Uncertainty Evaluation Brings Value to Wellbore Placement.” Paper presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, September 2006. doi: https://doi.org/10.2118/102088-MS
  8. Chen, David C-K, and Min Wu. “State-of-the-Art BHA Program Produces Un-Precedented Results.” Paper presented at the International Petroleum Technology Conference, Kuala Lumpur, Malaysia, December 2008. doi: https://doi.org/10.2523/IPTC-11945-MS
  9. Jamieson, Angus. “Introduction to Wellbore Positioning version”, 09.10.17, University of the Highlands and Islands. pp55-58
  10. Bulychenkov, K. “Bad News for High-Def Trajectory”, 2022. https://mwdstd.com/high-definition-trajectory-tvd-issue/
  11. Stockhausen, E. J., and W. G. Lesso. “Continuous Direction and Inclination Measurements Lead to an Improvement in Wellbore Positioning.” Paper presented at the SPE/IADC Drilling Conference, Amsterdam, Netherlands, February 2003. doi: https://doi.org/10.2118/79917-MS
  12. ISCWSA MWD Error Model Rev4, June 2015. https://www.iscwsa.net/files/486/
  13. Ploskikh, Viktor, and Konstantin Bulychenkov. “Building a Bulletproof System for Automatic MWD Survey Processing.” Paper presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, USA, October 2022. doi: https://doi.org/10.2118/210385-MS
  14. Jamieson, Angus: “Introduction to Wellbore Positioning”, version 09.10.17, University of the Highlands and Islands, p132 https://www.uhi.ac.uk/en/t4-media/one-web/university/research/eBook-V9.10.2017.pdf
  15. Bulychenkov, Konstantin. “Overcoming Fundamental Restrictions: Novel Multistation Measurement While Drilling Survey Correction.” SPE J. (2022;): doi: https://doi.org/10.2118/210084-PA
  16. Maus, Stefan, and Croke, Ryan. Field Acceptance Criteria Based on ISCWSA Tool Error Models. Oral presentation given at the 39th ISCWSA Meeting, Long Beach, May 9th, 2014, https://www.iscwsa.net/media/files/files/2965f750/meeting-39-field-acceptance-criteria-based-on-iscwsa-tool-error-models.pdf
  17. Bulychenkov, Konstantin. “Overcoming Fundamental Restrictions: Novel Multistation Measurement While Drilling Survey Correction.” SPE J. (2022;): doi: https://doi.org/10.2118/210084-PA

Appendix A – BHA Sag Correction Algorithm

BHA sag correction algorithm sag_correction(). The x-z coordinate system is selected such that the z-axis is aligned with the average inclination of the considered wellbore section, and the x-axis is perpendicular to the z-axis. The algorithm consists of the following steps:

  1. Divide the drillstring into uniform elements along the z-axis with a fixed step Δz
  2. Calculate the boundaries of the well walls:

where Dibh and xitr are the borehole diameter and well centerline x-coordinate at the ith element, respectively.

where incz is the inclination of z-axis of the selected coordinate system.

  1. Set the initial drillstring deflection vector:

  1. Calculate the gradient of the energy function (see the derivation of the gradient in Appendix B).

where qi is the unit weight of the ith element of the drillstring in the mud, g is gravity, Ii is the inclination at the ith element, dlsi is the DLS at the ith element, EIi is the momentum of inertia multiplied by the Young’s modulus of the ith element, and wob is the weight-on-bit force.

where qi0 is the linear weight of the ith element, ρi is the material density of the ith element, and ρm is the drilling-fluid density.

DLS is calculated as the second central derivative:

where dlsbend is the DLS of bend angle if applicable.

EI for a tubular element:

EI for a tubular element with blades (like string stabilizer):

where Ei is the Young’s modulus for the ith element, ODi is the outer diameter of the ith element, IDi is the inner diameter of the ith element, and ODblade is the blade diameter of the ith element, if applicable.

  1. Update the deflection vector:

where αj-1 specifies the rate of gradient descent and is selected based on the optimal rate of convergence of the algorithm, and j is the number of iterations.

  1. Update the deflection vector with the wellbore wall constraints:

  1. Check the algorithm stop criteria:

  1. Calculate the BHA sag:

where j is the index of the BHA element with the D&I module.

Appendix B – Derivation of Energy Gradient

To calculate the gradient of the energy function, we must consider the full energy of the drillstring. The full energy can be represented as the sum of the energies of all the elements:

where Ei is the energy of the ith element, which is the sum of the gravity, bending, and deflection energies of the weight-on-bit force (WOB). To simplify the algorithm, the contribution of the compression and shear energies are not considered because of their small values (compression energy fraction is of the order of 1e-6; shear energy fraction is of the order of 1e-3 to 1e-2):

The gravity energy can be determined as follows:

The bending energy is defined as follows:

The deflection energy is defined as follows:

DLS for ith element is determined as

The gradient function ∇E:

Because xi only affects the ith element itself and the bending energy of its nearest neighbors (owing to the smallness of xi, we can ignore the change in the deflection energy), the partial derivative of the full energy with respect to xi can be expressed as

Expand Equation 8b:

Appendix C – High-Definition Trajectory Algorithm

The HD algorithm is based on the fusion of static surveys and continuous inclination (CI) measurements. For algorithm stability, we expect the worst CI quality. We assume that

  • The axial accelerometer is subject to bias, scale factor, and misalignment errors;
  • These errors may drift with time;
  • The total gravity field reference for CI can change between static surveys.

The only constructive assumption is that the change in the CI is proportional to the change in the true inclination of the trajectory. As a result, the algorithm becomes immune to all potential errors listed above. The fusion algorithm consists of the following steps:

  1. Filter the CI measurements [13]
  2. Calculate additional points by interpolating at the static survey depths with the CI:

  1. Calculate the weighted sum of the static and interpolated inclinations:

where cijint is the interpolated inclination at the depth of the jth static survey, sij is the inclination of the jth static survey, σci is the uncertainty of the CI measurement, and σint is the interpolation uncertainty, calculated as follows:

where Δmdjmin is the depth distance to the closest CI point, and DLS is the BHA steering capacity.

  1. Add the fused continuous points to the set of the measured ones
  2. Split the CI points into groups between pairs of adjacent static surveys:

where cj is a subset of CI points {(md1, ci1)…(mdmj, cimj)} with mj number of measurements.

  1. Calculate a build rates subset for each subset cj:

  1. Calculate build rate offset broj for the jth interval:

  1. Calculate the HD inclination:

Appendix D – Trajectory Approximation Algorithm

In the vertical plane, the trajectory is approximated by a piecewise linear function I(D) using the least squares method, where Ii is the inclination and Di is the measured depth. Such a function can be defined by using a vector to solve a system of nonlinear fusion equations. Additionally, during the solution process, the BHA performance is calculated as follows, where dls is the steering DLS and br is the rotor build rate:

The observation vector for the system of equations is of the form

where the a priori vector Vapr consists of two parts: the a priori trajectory (Iiapr, Diapr) is interpolated with the HD trajectory for the depth values from the slide/steering sheet; dlsapr and brapr are a priori DLS performances calculated from the BHA design:

The delta steering observation vector ΔVstr represents the difference between the calculated and observer BHA build rates.

Vhd is the HD trajectory points measured by the tool (not interpolated):

The operator f(•)for a nonlinear system of equations is defined as follows:

where Djhd is the depth of the jth survey; αi is the gravity toolface of the ith slide interval.

The covariance matrix is represented as

where Capr is the covariance matrix for the a priori vector Vapr defined by the uncertainty of the BHA performance-prediction model and the interpolation accuracy, CΔstr is the covariance matrix for the uncertainty of the recalculated BHA steering performance, and Chd is the covariance matrix for the HD trajectory.

The solving vector Vslv is determined using the Gauss–Newton algorithm:

  1. Set the initial solution vector Vslv
  2. Calculate the Jacobian matrix J = ∇f(Vslv)
  3. Determine the delta of the observation vector ΔVobs = Vobs – f(Vslv)
  4. Calculate the posteriori covariance matrix Capost = (JTC-1J)-1
  5. Calculate the delta of the solution vector ΔVslv = Capost JTC-1ΔVobs
  6. Update the solution vector Vslv = Vslv + ΔVslv
  7. Check stop criteria as ΔVslvTΔVslv < ε
  8. Go to #2