Real-Time Interpretation Model of Reservoir Characteristics While Underbalanced Drilling Based on UKF

This study presents a novel interpretation model for reservoir characteristics while underbalanced drilling (UBD), by incorporating an unscented Kalman filter (UKF) algorithm in a three-phase variable mass flow model of oil, gas, and liquid. In the model, the measurement parameters are simplified to bottomhole pressure and liquid outlet flow, for decreasing the amount of the computation and time. By taking into account real-time measurements, the permeability and reservoir pressure along the well can be continuously updated. Three cases including single-parameter and doubleparameter estimations have been simulated, and the performance is tested against the extended Kalman filter (EKF). The results show that single-parameter estimation of reservoir permeability or pressure achieves superior performance. The filtered values of bottomhole pressure and outlet flow trace the measured values in real time. When a new section of a reservoir is opened, the estimated reservoir permeability or pressure can always be quickly and accurately returned to its true value. However, it is not possible for the double-parameter estimation to obtain good results; its interpretation accuracy is low. UKF is superior to EKF in both estimation accuracy and convergence speed, which further illustrates the superiority and accuracy of the novel interpretation model based on UKF. Benefits from this model are seen in accurate bottomhole pressure and reservoir characteristic predictions, which are of major importance for safety and economic reasons during UBD and follow-up completion operations.


Introduction
Underbalanced drilling (UBD) is a technical means of drilling under the condition of negative pressure in a wellbore [1,2]. UBD began in the 1930s and developed rapidly in the recent years [3], because of its great advantages [4] in improving rate of penetration (ROP), discovering reservoirs, and reducing reservoir damage. When a reservoir is opened, formation fluid will enter the wellbore under the action of negative pressure difference, which will lead to changes in annulus pressure, flow rate, and other parameters [5,6]. Because of the coupling effect between the wellbore and formation, reservoir characteristics such as reservoir pressure and permeability may be obtained by using measured data acquired during drilling operation, such as bottomhole pressure, flow rate, and injection rate.
Estimation of the near wellbore characteristics of a reservoir gives important information in the drilling and completion process, helping the technical crew make better decisions. It is of great significance for early identifying a reservoir, reducing drilling risk, and improving drilling time. On this basis, Kardolus and Kruijsdijk [7] put forward the idea of formation characteristic interpretation while underbalanced drilling. Differently from the traditional approach, it needed not to go down the test string; the permeability profile near the wellbore was estimated while drilling based on a simplified analytical model. Actually, the joint hydraulics model and estimation algorithm is the basic theory of parameter interpretation. Several researches have addressed the basic theory works, respectively, for a two-phase flow modeling or estimation algorithm of UBD. Rommetveit et al. [8] presented a transient UBD simulator, DynaFloDrill, and the predictions were validated by the full-scale experiments performed with either parasite string or drill string gas injection. Lage et al. [9] proposed a composite discrete scheme for dynamic two-phase flow modelling, combining the firstorder Lax-Friedrichs and the second-order MacCormack format, to describe the transient behavior in UBD. Fjelde et al. [10] established a new multiphase flow model of UBD, which used the MUSCL format to modify the classical upwind format, and could better describe the pressure fluctuation law during pipe connection. Perez-Tellez et al. [11] proposed a mechanical model for predicting the pressure of two-phase flow in an annulus and a drill pipe, in which a drift flux model is coupled. Khezrian et al. [12] developed a two-fluid model in the Eulerian frame of reference for simulation of gas-liquid two-phase flow in the UBD operation. On this basis, Vefring et al. [13] evaluated the performance of the nonlinear least square methodology. The Levenberg-Marquardt optimization algorithm was used to estimate reservoir pressure and permeability in the UBD process. Tang [14] applied the damped least square method to study the interpretation of reservoir characteristics in underbalanced drilling. However, those approaches are both nonrecursive, which need to use all historical data, so they are best suited for postanalysis of data.
Kalman filter (KF) estimation methods have been widely used in oil industry in recent years, because of their better performance in real time. Vefring et al. [15] are the representative researchers in reservoir characteristic interpretation while underbalanced drilling. The measurements were set up as pump pressure, bottomhole pressure, and gas and liquid outlet flow. Based on the Levenberg-Marquardt and extended Kalman filter (EKF), interpretation models of reservoir pressure and permeability were established, in which EKF was a realtime recursive algorithm. Nazari et al. [16] introduced the unscented Kalman filter (UKF) to estimate the gas-liquid mixing velocity in the drill string and annulus. Several pressure sensors were installed in the annulus to improve the accuracy and robustness. Lorentzen et al. [17] designed an EKF method based on a two-phase flow model to regulate the key parameters of a drift flux model in UBD operation. Nikoofard et al. [18,19] introduced a kind of nonlinear horizontal moving observer, which was used to estimate the liquid quality and liquid production coefficient of annular air in an underequilibrium state. Nygaard et al. [20] established a wellbore pressure control method based on the UKF algorithm to address the problem of pressure fluctuation caused by single connection and flow change when a gas reservoir was opened and conducted inversion interpretation for reservoir permeability. Gravdal et al. [21] established a pressure interpretation model by using the UKF algorithm, but only the friction coefficient was used as the estimated parameter.
Based on the previous researches, most of the Kalman filter studies on UBD use pump pressure, bottomhole pressure, and gas and liquid outlet flow parameters as the measurements. However, with more measurement parameters, the amount of the computation and time will be increased accordingly. The measurement parameters in the model need to be further optimized and evaluated due to their internal relations. Furthermore, currently there is still less research of using UKF in the UBD process.
In the present paper, a novel interpretation model for reservoir characteristics while underbalanced drilling is developed, by incorporating the UKF algorithm in a three-phase variable mass flow model of oil, gas, and liquid. The measurement parameters are simplified to bottomhole pressure and liquid outlet flow, to continuously update the permeability and reservoir pressure along the well. Based on this interpretation, the future state of the reservoir system near the wellbore can be predicted, helping the management of the wellbore during UBD and follow-up completion operations.
The outline of the paper is as follows. First, the threephase variable mass flow model of oil, gas, and liquid, including a drift flux model and frictional pressure loss model, is presented. Then, the real-time estimation method of UKF is described. Finally, results from simulations performed with three cases, like reservoir permeability and pressure singleparameter estimation and double-parameter estimation, are conducted. The performance of UKF is also evaluated against EKF in permeability single-parameter estimation.

Three-Phase Variable Mass Flow Model of Oil, Gas, and Liquid
Underbalanced drilling technology is realized by gas injection into a drill string. When underbalanced drilling encounters reservoirs, the oil and gas in the formation will continue to flow into the wellbore. With the elapse of time, reservoir gas and oil rise upward in the wellbore and reservoir opening length prolongs, resulting in a gradual increase in reservoir production. Therefore, the wellbore is actually a variable mass flow system consisting of injected gas, drilling fluid, produced gas, produced oil, and cutting multiphase components.
A mathematical model of three-phase variable mass flow of oil, gas, and liquid is established based on the theory of a wellbore multiphase flow and dynamic reservoir model [22,23]. The basic assumptions are as follows: (1) The wellbore fluid flows in one dimension, ignoring the radial flow change (2) Drilling fluid is water-based mud (WBM), set as the Herschel-Bulkley model, without considering the mass transfer between oil, gas, and liquid phases (3) Ignoring the influence of heat transfer between the wellbore and formation, one can calculate the temperature in the wellbore by a linear geothermal gradient (4) The effect of cuttings on wellbore flow is small, so it is not considered The mass conservation equation of injected gas is The mass conservation equation of drilling fluid is The mass conservation equation of produced gas is The mass conservation equation of produced oil is The three-phase momentum conservation equation of oil, gas, and water is The P-V-T equation is The gas-liquid drift flux model is calculated as The model of frictional pressure loss along the well is calculated as where A is the annular area. ρ g , ρ l , and ρ o are the density of gas, drilling fluid, and oil, respectively. α g , α l , and α o are the volume fraction of gas, drilling fluid, and oil, respectively. v g , v l , and v o are the actual flow rates of gas, drilling fluid, and oil, respectively. x ig and x fg are the mass fraction of injected gas and produced gas, respectively. q fg and q o are the influx rate of gas and oil phases. g is the gravitational acceleration. θ is the angle between the wellbore and horizontal direction. p f is the pressure drop. p ac is the acceleration pressure drop. M g is the molar mass of the gas. p is the wellbore pressure. Z is the deviation factor, solved by the PR-EOS model [24]. R is a general gas constant. c 0 is the gas phase distribution coefficient. v gr is the gas slip velocity. f is the fanning friction coefficient. ρ m is the density of gas-liquid mixture. v m is the velocity of gas-liquid mixture. D is the equivalent diameter.
Distribution coefficients, slip velocity parameters, and pressure drop model will change with the change of a gasliquid two-phase flow pattern. Therefore, accurate flow pattern identification is an important prerequisite for establishing a comprehensive multiphase flow mathematical model. The flow pattern transition of gas-liquid two-phase flow is a complex physical process. With the change of volume fraction, velocity, pressure, and relative position of the twophase medium, the shape of the interface changes, which leads to the change of the flow pattern. At present, there is no mature theoretical support for the physical mechanism of flow pattern transition. For different flow pattern changes, empirical formulas of relevant parameters are often fitted by means of experiments, and then, critical conditions of flow pattern transition are determined. According to previous research results [25][26][27], two-phase flow patterns in the vertical wellbore are divided into five categories: bubble flow, dispersed bubble flow, slug flow, churn flow, and annular flow. For different flow patterns, the distribution coefficients, slip velocity, and pressure drop along the well are determined. The multiphase flow model is solved by a simple front tracking technique and finite difference numerical method. The specific algorithm and formula are not described here. A previous paper by He et al. [28] has introduced the algorithms and formulas in detail.

Dynamic Reservoir Model.
When reservoir oil and gas influx occurs in the drilling process, the influx mode is negative pressure influx, which can be regarded as plane radial flow in isotropic homogeneous elastic porous media. The Dake model [29] with an analytic solution is used to describe the flow process. In the process of underbalanced drilling, the opened zones of the reservoir are all involved in the coupled flow with the wellbore, thus forming the whole variable mass flow process. The opened reservoir is divided into n units along the axis of the wellbore, as shown in Figure 1. The dynamic reservoir model can be expressed as follows: where p a is the annular pressure. p r is the reservoir pressure. q is the influx rate. k is the reservoir permeability. S is the skin factor. Δh is the length of the open reservoir interval per unit time. t is the duration of the open reservoir interval. φ is the reservoir porosity. μ is the viscosity of reservoir fluid. c is the compressibility coefficient of reservoir fluid. r w is the borehole radius. γ is the Euler constant.

Interpretation Model of Reservoir
Characteristics Based on UKF 3.1. Determination of Interpretation Model Parameters. The description of the gas injection in the drill string of underbalanced drilling technology is as follows: gas passes through the gas injection pipeline, mixes with the drilling fluid pumped by a mud pump, and jointly injects into the drill string. With the opening of the reservoir section during drilling, reservoir oil and gas fluids gradually influx into the wellbore and return to the wellhead together with the injected gas and drilling fluids. After passing through the gas-liquid separator on the ground, the gas phase and liquid phase are finally separated. The schematic diagram is shown in Figure 2.
In the previous literature, pump pressure, bottomhole pressure, and gas and liquid outlet flow rate were often regulated as measurement parameters. The bottomhole pressure is measured by pressure while drilling (PWD) instruments. However, according to the principle of the U-tube, the pressure transfer relation can be constructed in the drill string between bottomhole pressure and pump pressure. In addition, there is a mathematical relationship between the gasliquid outlet flow rate and wellbore-reservoir coupled flow system. Because the mass transfer of gas and liquid is not considered, according to the superposition principle, the sum of gas-liquid outlet flow is equal to the sum of inlet flow, reservoir influx rate, and gas expansion rate. The formulas are as follows: Therefore, the measurement parameters can be simplified into two categories: bottomhole pressure and liquid outlet flow. The expressions are as follows: where p b is the bottomhole pressure. p s is the pump pressure. p h is the hydrostatic column pressure in the drill string. Δp d is the circulating pressure loss in the drill string. Δp t is the pressure loss of the drill bit. q out,l and q out,g represent liquid and gas outlet flow rates, respectively. q in,l and q in,g indicate liquid and gas inlet flow, respectively. q inf is the reservoir influx flow rate. q exp is the gas expansion rate. t 1 , t 2 , ⋯, t N is the N points of time corresponding to the measured data. Generally, reservoir lithology is complex, and the internal pore and fracture distribution is random and nonuniform. Therefore, the reservoir to be drilled can be subdivided into several units. Reservoir permeability k and reservoir pressure p r of each unit may be different, as shown in Figure 2. Reservoir permeability and pressure are determined as interpretation parameters of the model, and the expressions are as follows: 3.2. UKF Algorithm. The initial Kalman filtering technique is only suitable for linear systems [30], but most of the real problems are essentially nonlinear. Then, an extended Kalman filter and unscented Kalman filter are developed successively for nonlinear systems [31][32][33]. The main difference between EKF and UKF is the way Gaussian random variables are represented for propagating through system dynamics. UKF adopts a deterministic sampling method with a minimal set of carefully chosen sample points, instead of the local linearization in EKF. These sample points are propagated through the nonlinear system and capture the posterior mean and covariance accurately to the third order for all nonlinearities. In contrast to EKF, UKF appears later and requires no calculation of the complex Jacobian matrix and is a more advanced nonlinear filtering technology. Using UKF, it is possible to combine the information obtained from the measurements with the model to get an improved real-time estimate of the state vector of the system. The state vector is the interpretation parameter. Combined with the discussion in Section 3.1, the state vector x and the measurement z in the model are determined as follows,

Geofluids
where the subscript k denotes the time k.
The state space form of the nonlinear system based on state estimation is as follows: where x k is a state vector. f ð⋅Þ is a nonlinear system state function. z k is a measurement vector. hð⋅Þ is a nonlinear measurement function, which represents the variable mass flow model in this work. w k and v k are, respectively, the process noise and measurement noise of the system, which satisfy the zero mean white noise distribution, w k~N f0, Q k g, and v k~N f0, R k g, and they are not related to each other. The structure of UKF is basically the same as that of EKF, which is mainly divided into two processes: state update and measurement update. The process of state update means that the estimated value at the previous time is known, the current state variables and the estimated error covariance matrix are predicted by using the state model, and a prior estimate is constructed for the next time state. The measurement update process is responsible for feedback, and the posterior estimate of the current state is corrected by combining the prior estimate with a new measurement data. Therefore, this algorithm can also be called a predictor-corrector algorithm.
Firstly, the prior estimation mean of the state vector at the initial time is set tox 0 , and the estimation error covariance matrix P 0 is set.
According to the state vector of timestep k − 1, the estimation of timestep k can be obtained.
2L + 1 sigma points x i are constructed, where L is the dimension of the state vector, and the set of sigma points can be expressed aŝ where ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðL + λÞP k−1 p Þ i is the ith column of the matrix square root, which is calculated by the Cholesky decomposition. α represents the spread of the sigma points around the current state vector and is usually set as a small positive value, ranging from ½0, 1. λ = α 2 ðL + κÞ − L is a scaling parameter. κ is a secondary scaling parameter.
For each sigma point, the transformed result is obtained by nonlinear transformation f ð⋅Þ.
x i,k/k−1 is weighted to predict the prior estimates of mean x k/k−1 and covariance matrix P k/k−1 .
where W ðmÞ i and W ðcÞ i are weighted by mean and covariance, respectively. β is a scaling parameter, used to include the information about the distribution. For the Gaussian distribution, β = 2 is optimal.
Similarly, the nonlinear measurement function hð⋅Þ is used to transfer the sigma point to z i,k/k−1 and to predict the measurement valueẑ k/k−1 , the autocovariance matrix P z k , and the crosscovariance matrix P x k z k .
Finally, calculate the gain matrix K k , and update the state meanx k and covariance matrix P k at timestep k based on a new measured value z k .

Numerical Simulation Analysis
Based on the interpretation model while underbalanced drilling established above, real-time estimation of reservoir characteristics is carried out. Because the interpretation parameters of the model are reservoir pressure and permeability, they are difficult to be obtained under the current laboratory and field test conditions. Thus, this paper uses the synthetic data instead of the experimental data as measurements. Calculation steps are as follows: Firstly, a set of reservoir characteristic parameters (pressure and permeability) are preset as the estimates of interpretation parameters of the open reservoir section at the initial time. Then, the bottomhole pressure and outlet flow are simulated by using the established hydraulics model as measurements. Finally, the interpretation parameters are online estimated by using the UKF algorithm.
Well X is drilled in underbalanced condition by gas injection in the drill string. The injected gas phase is nitrogen. The reservoir is drilled at 3000 m, which is a pure oil reservoir. The reservoir is a typical medium-permeability sandstone reservoir with permeability ranging from 50 mD to 500 mD. Shale is the upper caprock, and low-permeability mudstone is the interlayer. Vertical distribution of sand and mud is interactive. The basic data needed for the well calculation are shown in Table 1.
The drilling process of the 3000-3100 m reservoir section is simulated. By dividing the reservoir into five units, each unit is 20 m in length. The reservoir pressure and permeability of each unit are constant. In order to analyze the influence of interpretation parameter dimension on the results of model estimation, we have carried out three cases, each of which has different interpretation parameters. Case 1 and case 2 are single-parameter estimates. The single parameter of reservoir permeability or formation pressure is estimated as the interpretation parameter, respectively. Case 3 is a double-parameter estimate, and the reservoir permeability and pressure are both estimated as interpretation parameters. Note that the proper specification of the covariance matrix for the modeling error is crucial to get better performance of the filter. The covariance matrix for state parameter error Q and measurement error R are almost diagonal, which depends on the specific case of the research object, while the uncertainty parameters may be related to each other. Although the relationship between those parameters and the effect on the interpretation performance are not the focus of this paper, they can be a topic for further research. We assume that the errors in the measurements and state parameter are statistically independent and are set as a suitable value with personal experience. The standard deviations of reservoir pressure and permeability process noise are 0.15 MPa and 0.5 mD, respectively, while the accuracy of bottomhole pressure and outlet flow measurement is 0.05% and 0.1%, respectively; i.e., the standard deviations of measurement noise are 0.02 MPa and 0.00003 m 3 /s, respectively. Thus, the covariance matrix for parameter and measurement error in three cases is as follows: Cases 1, 2, and 3 show that gas is injected from the drill string and then returned from the annulus wellhead, reaching to a steady state (0-75 min). The gas injection process is the common part of three cases, which are described in detail here. 0-15.5 min is the process of injected gas from the drill string to the drill bit. At this stage, the bottomhole pressure increases slightly and basically keeps unchanged, and the outlet flow increases with the gas injection at the initial time, then gradually decreases. 15.5-51.1 min is the process of gas migration from the bottomhole to the wellhead annulus. In this stage, the bottomhole pressure continues to decrease, while the outlet flow gradually increases due to the gas migration. 51.1-75 min is the process of gas front reaching wellhead annulus to wellbore stable-state flow. Bottomhole pressure continues to decrease at first, then tends to be constant, and the outlet flow quickly falls back to the inlet pump rate of 0.03 m 3 /s. For case 1, the reservoir pressure is known to be 38 MPa and the reservoir permeability is estimated. The calculation results are shown in Figures 3-5. Figures 3 and 4 show the comparison between the estimated data and measurements of bottomhole pressure and outlet flow rate, respectively. It can be seen that 75 min is the initial point of drilling in the reservoir section and real-time interpretation starts synchronously at this time. 450 minutes is the end point, when the reservoir is just opened 100 m in length. The estimated value of bottomhole pressure and outlet flow basically coincides with the measured value. Because of the influence of measurement noise, the measured value is somewhat burr, while the filtered curve (i.e., the blue line) is smoother. However, the calculated values of bottomhole pressure and outlet flow rate without filter processing (i.e., the green line) deviate from the measurements. The setting condition of the nonfilter is that the reservoir permeability is a constant value. Figure 5 shows the comparison between the true and estimated permeability values. The simulation results show that the true values of permeability of each reservoir unit are [250, 500, 350, 150, 50] mD and the initial value is set to 100 mD. 75-150 min is the time period of opening reservoir unit 1. With the elapse of time, the estimated value of reservoir permeability gradually rises from 100 mD to its true value near 150 mD and then basically fluctuates slightly around the true value. When unit 2 is opened (150-225 min), the reservoir permeability gradually approaches its true value of 500 mD from 150 mD, and then, analogies are made until the whole 100 m reservoir section is opened. Therefore, the estimated value of reservoir permeability is in good agreement with the true value. When a new section of the reservoir is opened, the estimated value can always be accurately and quickly returned to its true value.
For case 2, the reservoir permeability is preset as 400 mD and the reservoir pressure is estimated. The calculation results are shown in Figures 6-8. Figures 6 and 7 show the comparison between the estimated data and measurements of bottomhole pressure and outlet flow rate, respectively. The estimated values of bottomhole pressure and outlet flow are basically in agreement with the measurements, while the nonfiltered data deviate greatly from the measurements. Figure 8 shows the comparison between the true and estimated reservoir pressure values. It can be seen that the true values of reservoir pressure in each unit are [38,44,40,36,42] MPa and the initial value is set to 34 MPa. The estimation of reservoir pressure is in good agreement with the real value. When a new section of the reservoir is opened, the estimation of reservoir pressure can quickly approach its true value. Compared with example 1, it takes shorter time and converges faster.
The double parameters of reservoir permeability and pressure are estimated simultaneously in case 3. The simulation results are shown in Figures 9-12. Figures 9 and 10          In summary, the single-parameter estimation of reservoir permeability or pressure in case 1 and case 2 can achieve realtime and accurate interpretation. However, the doubleparameter estimation of reservoir permeability and pressure in case 3 cannot be realized and the interpretation accuracy is low. The reason is that for the two-parameter estimation, the reservoir permeability and pressure are only related to the crude output according to the dynamic reservoir model. Two measurements in this model are too few to obtain enough effective information. However, it is also difficult to obtain effective information for the four measurements in the previous researches. It is not feasible to simply increase the number of measurement parameters, which are both related to pressure or flow rate. Therefore, by increasing measurement dimensions (e.g., dielectric constant or resistivity) and digging useful measurements in depth, synchronous estimation of two or even multiple parameters of reservoir characteristics may be realized.
Taking the single-parameter estimation of reservoir permeability in case 1 as an example, the UKF algorithm and EKF algorithm are compared and analyzed, as shown in Figure 13. The simulation results are shown in Table 2. The average reservoir permeability while drilling through each reservoir unit is used as the estimated value. Compared with EKF, the estimated values of UKF in units 1-5 are both closer to the true value. Moreover, it takes less time for UKF estimation error to converge to 10% than it is for EKF. Only when reservoir unit 5 is used, time of convergence to 10% is higher than that for EKF, which may be related to the randomness of model noise. The numerical results show that UKF is superior to EKF in both estimation accuracy and convergence speed, which further indicates the superiority and accuracy of this model.

Conclusion
(1) A novel interpretation model while underbalanced drilling based on the UKF and three-phase variable mass flow model of oil, gas, and liquid is established. The measurement parameters are simplified to bottomhole pressure and liquid outlet flow, to continuously update the permeability and reservoir pressure along the well, helping the management of the wellbore during the UBD and follow-up completion operations (2) Three cases including reservoir permeability and pressure single-parameter estimation and doubleparameter estimation are simulated. The results show that the interpretation accuracy of single-parameter estimation is high and the filtered values of 12 Geofluids bottomhole pressure and outlet flow trace the measured values in real time. When a new section of the reservoir is opened, the estimated reservoir permeability or pressure can always be quickly and accurately returned to its true value. However, the double-parameter estimation cannot be realized, and the interpretation accuracy is low (3) The comparison between the UKF and the EKF is carried out. The results show that UKF is superior to EKF in both estimation accuracy and convergence speed, which further illustrates the superiority and accuracy of the new model for interpretation while underbalanced drilling based on UKF (4) Although the presented study is promising, for addressing the problem of joint state-andparameter estimation of a multiphase flow during drilling, further development of the methodology is needed. Future works will include how to achieve the double-parameter or even multipleparameter estimation for reservoir characteristics and the application of the methodology to actual oilfield experiments

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.