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

Authors: Miao He; Yihang Zhang; Mingbiao Xu; Jun Li; Jianjian Song
Journal: Geofluids
 (2020)

---

## Abstract

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 double-parameter 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.


---

## Body

## 1. 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 first-order 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 real-time 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 three-phase 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 single-parameter estimation and double-parameter estimation, are conducted. The performance of UKF is also evaluated against EKF in permeability single-parameter estimation.


## 2. 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


### 2.1. Multiphase Flow Equation

The mass conservation equation of injected gas is(1)∂∂txigρgαgA+∂∂zxigρgαgvgA=0.The mass conservation equation of drilling fluid is(2)∂∂tρlαlA+∂∂zρlαlvlA=0.The mass conservation equation of produced gas is(3)∂∂txfgρgαgA+∂∂zxfgρgαgvgA=qfg.The mass conservation equation of produced oil is(4)∂∂tρoαoA+∂∂zρoαovoA=qo.The three-phase momentum conservation equation of oil, gas, and water is(5)∂∂tρgαgvgA+ρlαlvlA+ρoαovoA+∂∂zρgαgvg2A+ρlαlvl2A+ρoαovo2A+ρgαg+ρlαl+ρoαogsinθA+∂pA∂z+A∂pf∂z=0.TheP-V-T equation is
(6)ρg=MgpZRT.The gas-liquid drift flux model is calculated as(7)vg=c0vm+vgr.The model of frictional pressure loss along the well is calculated as(8)∂pf∂z=2fρmvmDvm+∂pac∂z,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. vg, vl, and vo are the actual flow rates of gas, drilling fluid, and oil, respectively. xig and xfg are the mass fraction of injected gas and produced gas, respectively. qfg and qo are the influx rate of gas and oil phases. g is the gravitational acceleration. θ is the angle between the wellbore and horizontal direction. pf is the pressure drop. pac is the acceleration pressure drop. Mg 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. c0 is the gas phase distribution coefficient. vgr is the gas slip velocity. f is the fanning friction coefficient. ρm is the density of gas-liquid mixture. vm 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 gas-liquid 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 two-phase 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–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.


### 2.2. 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:
(9)qiti=4πkiΔhpr,i−pa,iμ2S+ln4kiti/eγϕμcrw2,where pa is the annular pressure. pr 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. rw is the borehole radius. γ is the Euler constant.Figure 1
Diagram of wellbore-reservoir coupling variable mass flow.


## 2.1. Multiphase Flow Equation

The mass conservation equation of injected gas is(1)∂∂txigρgαgA+∂∂zxigρgαgvgA=0.The mass conservation equation of drilling fluid is(2)∂∂tρlαlA+∂∂zρlαlvlA=0.The mass conservation equation of produced gas is(3)∂∂txfgρgαgA+∂∂zxfgρgαgvgA=qfg.The mass conservation equation of produced oil is(4)∂∂tρoαoA+∂∂zρoαovoA=qo.The three-phase momentum conservation equation of oil, gas, and water is(5)∂∂tρgαgvgA+ρlαlvlA+ρoαovoA+∂∂zρgαgvg2A+ρlαlvl2A+ρoαovo2A+ρgαg+ρlαl+ρoαogsinθA+∂pA∂z+A∂pf∂z=0.TheP-V-T equation is
(6)ρg=MgpZRT.The gas-liquid drift flux model is calculated as(7)vg=c0vm+vgr.The model of frictional pressure loss along the well is calculated as(8)∂pf∂z=2fρmvmDvm+∂pac∂z,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. vg, vl, and vo are the actual flow rates of gas, drilling fluid, and oil, respectively. xig and xfg are the mass fraction of injected gas and produced gas, respectively. qfg and qo are the influx rate of gas and oil phases. g is the gravitational acceleration. θ is the angle between the wellbore and horizontal direction. pf is the pressure drop. pac is the acceleration pressure drop. Mg 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. c0 is the gas phase distribution coefficient. vgr is the gas slip velocity. f is the fanning friction coefficient. ρm is the density of gas-liquid mixture. vm 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 gas-liquid 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 two-phase 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–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.


## 2.2. 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:
(9)qiti=4πkiΔhpr,i−pa,iμ2S+ln4kiti/eγϕμcrw2,where pa is the annular pressure. pr 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. rw is the borehole radius. γ is the Euler constant.Figure 1
Diagram of wellbore-reservoir coupling variable mass flow.


## 3. 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 Figure2.Figure 2
Diagram of underbalanced drilling of gas injection in the drill string.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 gas-liquid 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:(10)pb=ps+ph−Δpd−Δpt,qout,l+qout,g=qin,l+qin,g+qinf+qexp.Therefore, the measurement parameters can be simplified into two categories: bottomhole pressure and liquid outlet flow. The expressions are as follows:(11)pb=pbt1,pbt2,⋯,pbtN,qout,l=qout,lt1,qout,lt2,⋯,qout,ltN,where pb is the bottomhole pressure. ps is the pump pressure. ph is the hydrostatic column pressure in the drill string. Δpd is the circulating pressure loss in the drill string. Δpt is the pressure loss of the drill bit. qout,l and qout,g represent liquid and gas outlet flow rates, respectively. qin,l and qin,g indicate liquid and gas inlet flow, respectively. qinf is the reservoir influx flow rate. qexp is the gas expansion rate. t1,t2,⋯,tN 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 permeabilityk and reservoir pressure pr 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:
(12)k=kt1,kt2,⋯,ktN,pr=prt1,prt2,⋯,prtN.


### 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–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 Section3.1, the state vector x and the measurement z in the model are determined as follows, where the subscript k denotes the time k.
(13)x=ktk,prtkT,z=pbtk,pout,ltkT.The state space form of the nonlinear system based on state estimation is as follows:(14)xk=fxk−1+wk−1,zk=hxk+vk,where xk is a state vector. f⋅ is a nonlinear system state function. zk is a measurement vector. h⋅ is a nonlinear measurement function, which represents the variable mass flow model in this work. wk and vk are, respectively, the process noise and measurement noise of the system, which satisfy the zero mean white noise distribution, wk~N0,Qk, and vk~N0,Rk, 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 P0 is set.According to the state vector of timestepk−1, the estimation of timestep k can be obtained.
(15)x^k=x^k−1,Pk=Pk−1+Qk−1.2L+1 sigma points xi are constructed, where L is the dimension of the state vector, and the set of sigma points can be expressed as
(16)x^i,k−1=x^k−1,i=0,x^k−1+L+λPk−1i,i=1,⋯,L,x^k−1−L+λPk−1i,i=n+1,⋯,2L,where L+λPk−1i 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. λ=α2L+κ−L is a scaling parameter. κ is a secondary scaling parameter. If L≤3, κ=3‐L; if L>3, κ=0.For each sigma point, the transformed result is obtained by nonlinear transformationf⋅.
(17)xi,k/k−1=fxi,k−1,i=0,1,⋯,2L.xi,k/k−1 is weighted to predict the prior estimates of mean x^k/k−1 and covariance matrix Pk/k−1.
(18)x^k/k−1=∑i=02LWimxi,k/k−1,Pk/k−1=∑i=02LWicxi,k/k−1−x^k/k−1xi,k/k−1−x∧k/k−1T+Qk−1,Wim=λL+λ,i=0,12L+λ,i≠0,Wic=λL+λ+1+β−α2,i=0,12L+λ,i≠0,where Wim and Wic 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 functionh⋅ is used to transfer the sigma point to zi,k/k−1 and to predict the measurement value z^k/k−1, the autocovariance matrix Pzk, and the crosscovariance matrix Pxkzk.
(19)zi,k/k−1=hxi,k−1,i=0,1,⋯,2L,z^k/k−1=∑i=02LWimzi,k/k−1,Pzk=∑i=02LWiczi,k/k−1−z^k/k−1zi,k/k−1−z∧k/k−1T+Rk,Pxkzk=∑i=02LWicxi,k/k−1−x^k/k−1zi,k/k−1−z∧k/k−1T.Finally, calculate the gain matrixKk, and update the state mean x^k and covariance matrix Pk at timestep k based on a new measured value zk.
(20)Kk=PxkzkPzk−1,x^k=x^k/k−1+Kkzk−z^k/k−1,Pk=Pk/k−1−KkPzkKkT.


## 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 Figure2.Figure 2
Diagram of underbalanced drilling of gas injection in the drill string.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 gas-liquid 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:(10)pb=ps+ph−Δpd−Δpt,qout,l+qout,g=qin,l+qin,g+qinf+qexp.Therefore, the measurement parameters can be simplified into two categories: bottomhole pressure and liquid outlet flow. The expressions are as follows:(11)pb=pbt1,pbt2,⋯,pbtN,qout,l=qout,lt1,qout,lt2,⋯,qout,ltN,where pb is the bottomhole pressure. ps is the pump pressure. ph is the hydrostatic column pressure in the drill string. Δpd is the circulating pressure loss in the drill string. Δpt is the pressure loss of the drill bit. qout,l and qout,g represent liquid and gas outlet flow rates, respectively. qin,l and qin,g indicate liquid and gas inlet flow, respectively. qinf is the reservoir influx flow rate. qexp is the gas expansion rate. t1,t2,⋯,tN 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 permeabilityk and reservoir pressure pr 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:
(12)k=kt1,kt2,⋯,ktN,pr=prt1,prt2,⋯,prtN.


## 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–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 Section3.1, the state vector x and the measurement z in the model are determined as follows, where the subscript k denotes the time k.
(13)x=ktk,prtkT,z=pbtk,pout,ltkT.The state space form of the nonlinear system based on state estimation is as follows:(14)xk=fxk−1+wk−1,zk=hxk+vk,where xk is a state vector. f⋅ is a nonlinear system state function. zk is a measurement vector. h⋅ is a nonlinear measurement function, which represents the variable mass flow model in this work. wk and vk are, respectively, the process noise and measurement noise of the system, which satisfy the zero mean white noise distribution, wk~N0,Qk, and vk~N0,Rk, 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 P0 is set.According to the state vector of timestepk−1, the estimation of timestep k can be obtained.
(15)x^k=x^k−1,Pk=Pk−1+Qk−1.2L+1 sigma points xi are constructed, where L is the dimension of the state vector, and the set of sigma points can be expressed as
(16)x^i,k−1=x^k−1,i=0,x^k−1+L+λPk−1i,i=1,⋯,L,x^k−1−L+λPk−1i,i=n+1,⋯,2L,where L+λPk−1i 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. λ=α2L+κ−L is a scaling parameter. κ is a secondary scaling parameter. If L≤3, κ=3‐L; if L>3, κ=0.For each sigma point, the transformed result is obtained by nonlinear transformationf⋅.
(17)xi,k/k−1=fxi,k−1,i=0,1,⋯,2L.xi,k/k−1 is weighted to predict the prior estimates of mean x^k/k−1 and covariance matrix Pk/k−1.
(18)x^k/k−1=∑i=02LWimxi,k/k−1,Pk/k−1=∑i=02LWicxi,k/k−1−x^k/k−1xi,k/k−1−x∧k/k−1T+Qk−1,Wim=λL+λ,i=0,12L+λ,i≠0,Wic=λL+λ+1+β−α2,i=0,12L+λ,i≠0,where Wim and Wic 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 functionh⋅ is used to transfer the sigma point to zi,k/k−1 and to predict the measurement value z^k/k−1, the autocovariance matrix Pzk, and the crosscovariance matrix Pxkzk.
(19)zi,k/k−1=hxi,k−1,i=0,1,⋯,2L,z^k/k−1=∑i=02LWimzi,k/k−1,Pzk=∑i=02LWiczi,k/k−1−z^k/k−1zi,k/k−1−z∧k/k−1T+Rk,Pxkzk=∑i=02LWicxi,k/k−1−x^k/k−1zi,k/k−1−z∧k/k−1T.Finally, calculate the gain matrixKk, and update the state mean x^k and covariance matrix Pk at timestep k based on a new measured value zk.
(20)Kk=PxkzkPzk−1,x^k=x^k/k−1+Kkzk−z^k/k−1,Pk=Pk/k−1−KkPzkKkT.


## 4. 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 Table1.Table 1
Basic parameters for well X.
ParametersValueParametersValueWell depth (m)3000Yield point (Pa)3.84Casing shoe depth (m)2200Surface temperature (°C)20Casing ID (mm)228.47Geothermal gradient (°C/m)0.0266Bit diameter (mm)215.9Casing pressure (MPa)1Drill pipe OD (mm)127Gas injection rate (Nm3/s)0.5Drill pipe ID (mm)108.6ROP (m/h)16Bit nozzle area (mm2)660Oil density (kg/m3)850Pump rate (L/s)30Oil viscosity (mPa·s)30Mud density (kg/m3)1180Compressibility of liquid (Pa-1)1.0×10−9Consistency index (Pa·sn)0.37Reservoir porosity0.15Flow behavior index0.68Skin factor0.013The 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. Case1 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. The true value of reservoir permeability and pressure of five units is given in each case: case 1: reservoir permeability is [250, 500, 350, 150, 50] mD and reservoir pressure has a fixed value (38 MPa); case 2: reservoir pressure is [38, 44, 40, 36, 42] MPa and reservoir permeability has a fixed value (400 mD); and case 3: reservoir permeability is [250, 500, 350, 150, 50] mD and reservoir pressure is [38, 44, 40, 36, 42] MPa.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 errorQ 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 m3/s, respectively. Thus, the covariance matrix for parameter and measurement error in three cases is as follows:Case 1.
Q=diag0.52, R=diag0.022,0.000032Case 2.
Q=diag0.152, R=diag0.022,0.000032Case 3.
Q=diag0.52,0.152, R=diag0.022,0.000032Cases1, 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 m3/s.For case1, 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 3
Bottomhole pressure variation (case1).Figure 4
Liquid outlet rate variation (case1).Figure 5
True and estimated reservoir permeability (case1).Figure5 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 case2, 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.Figure 6
Bottomhole pressure variation (case2).Figure 7
Liquid outlet rate variation (case2).Figure 8
True and estimated reservoir pressure (case2).The double parameters of reservoir permeability and pressure are estimated simultaneously in case3. The simulation results are shown in Figures 9–12. Figures 9 and 10, respectively, show the comparison between the estimated data and measurements of bottomhole pressure and outlet flow. Similarly, the estimated values of bottomhole pressure and outlet flow rate are highly consistent with the measured values. Bottomhole pressure gradually decreases with the increase in crude oil production during 150-225 min. In 225-300 min, the effects of the increase in crude oil production and the prolongation of the reservoir section on the bottomhole pressure offset each other, thus basically remaining unchanged. In 300-450 min, the prolongation of the reservoir section becomes the main controlling factor, which results in the gradual increase in bottomhole pressure. It is obvious that there are five different stages of slope change in the outlet flow rate, which matches with five different reservoir units drilled.Figure 9
Bottomhole pressure variation (case3).Figure 10
Liquid outlet rate variation (case3).Figure 11
True and estimated reservoir permeability (case3).Figure 12
True and estimated reservoir pressure (case3).Figures11 and 12 show the contrasting results of true and estimated reservoir permeability and pressure, respectively. It can be seen that the true values of reservoir permeability and pressure of each unit are [250, 500, 350, 150, 50] mD and [38, 44, 40, 36, 42] MPa, respectively, and the initial values are set to 100 mD and 34 MPa, respectively. The estimated values of reservoir permeability and pressure are quite different from the true values, which basically do not coincide with each other. When the permeability estimates are lower than the true values, the reservoir pressure estimates are higher and vice versa.In summary, the single-parameter estimation of reservoir permeability or pressure in case1 and case 2 can achieve real-time and accurate interpretation. However, the double-parameter 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 case1 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.Figure 13
Comparison between EKF and UKF of the true and estimated value (case1).Table 2
Comparative analysis of EKF and UKF estimated results.
Reservoir unitReservoir permeability (mD)Time of convergence to 10% (min)True valueUKFEKFUKFEKF1250218213212625004504392126335037938618254150189198252955069712725


## 5. 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 double-parameter estimation are simulated. The results show that the interpretation accuracy of single-parameter estimation is high and the filtered values of 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-and-parameter 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 multiple-parameter estimation for reservoir characteristics and the application of the methodology to actual oilfield experiments



---

*Source: 8967961-2020-05-19.xml*