Adaptive Moving Horizon Estimator for Return Flow Rate Estimation using Fluid Levels of a Venturi Channel

Real-time estimation of the return drilling fluid during oil well drilling is investigated in this study. Online fluid level measurements from a Venturi channel which can be placed on the return flowline is used with a model-based estimator. A reduced order, 1-D, mathematical equation is used for the open flow in the Venturi channel for Newtonian or non-Newtonian fluid types. The volumetric fluid flow rate is estimated using a moving horizon estimator in real-time. The friction factor is also estimated together with the fluid flow rate. The effect of the variation of the channel slope on the flow rate estimation induced by the vibration of the channel during its operation is also studied. The method requires only two level measurements in the Venturi channel together with the channel geometry. The method is validated using a laboratory scale Venturi flow system. The proposed method shows promising potential to be used as a real-time return flow rate measurement in conventional drilling systems.


Introduction
The return flow rate and the mud pit level are the two main kick indicators in conventional drilling systems. The return flow consists of drill cuttings and gases which make flow measurement difficult and inaccurate. Although there exist many flowmeters that can measure the return flow rate, most of the on-shore and off-shore conventional drilling systems use paddle flow sensors. This is just an indicator rather than a real-time flowmeter, thus early kick and loss detection cannot be expected. Therefore, the development of a real-time flow sensor for conventional drilling systems is advantageous for reduction of risk and non-productive time via early kick detection.
The return flow rate can be measured by placing a flow sensor on the return flowline which is an open channel. Use of open channels for fluid flow measurements is a well-known practice in the hydropower industry and agricultural industries. Generally, the flow measurement proposition using open channels is quite complex compared to the pressurized flow in pipes, especially because there are many variables associated with these flow channels, mainly the changes of the free surface. However, for most of the cases, it is possible to approximate and express these variables utilizing the continuity and energy equations of fluid mechanics. There exist a number of established approaches for open channel flow measurement, mostly based on steady flow conditions such as the use of a weir, the volumetric tank method or the area velocity method (Basu, 2019).
Further, empirical equations for prismatic open channels such as the Chezy equation or Manning equation (Chanson, 2004;Chow, 1959;Chaudhry, 2008) for Newtonian fluids and Haldenwang equations (Haldenwang and Slatter, 2006;Burger et al., 2010Burger et al., , 2015aHaldenwang et al., 2010) for non-Newtonian fluids are used in literature. These equations use one level measurement from the channel and are developed for steady flows for specific geometrical channels. For nonprismatic channels, especially for Venturi channels, a steady flow equation based on the Bernoulli principle can be derived for two level measurements with several empirical coefficients that need to be tuned (Pirir et al., 2017;Chhantyal, 2018). Flume equations based on one level measurement and empirical coefficients developed for specific Venturi flumes that are designed according to the ISO 4359:2013 standards are available (International Organization for Standardization [ISO], 2013;Basu, 2019;Chhantyal, 2018;Baker, 2016). However, these perform best for steady or slowly varying flow conditions only. Flow equations can be developed based on the critical level of the channel, yet these are quite impractical because the critical level position can vary along the channel due to the changes in the flow rate and flow conditions (Chhantyal, 2018;Chhantyal et al., 2017). Similar techniques and equations are discussed in detail in Boiten (2002); Basu (2019); Baker (2016); Henderson (1966); Alderman and Haldenwang (2007); Haldenwang et al. (2002). Most of these methods are quite specific and only works for certain types of fluids and flow conditions. Therefore, a versatile, complete and dynamic open channel flow model should be investigated.
Nevertheless, the dynamic modeling and simulation of free surface flows are complex and challenging. Generally, the open channel hydraulics are modeled by the well known and efficient shallow water equations, which are a set of nonlinear, hyperbolic partial differential equations (Chaudhry, 2008). Although these equations are widely used throughout history, they are still difficult to solve. There are many numerical methods of high precision for solving the shallow water equations, but these usually take a considerable amount of computational time, which makes them not suitable for realtime applications. A model reduction method to calculate the fluid flow using these models with an application to oil well drilling process has been studied earlier (Jinasena et al., , 2017. These reduced order models have also been used for real-time estimation using various estimation methods . However, the parameters such as friction factor are considered to be unknown and may vary with time, depending on the operational conditions such as the flow types and different fluid properties (Jinasena et al., 2019). Thus, the friction factor should be estimated continuously (in real-time) together with the state estimation. The motivation for this work is to estimate fluid flow rate and parameters using a moving horizon estimator in real-time, thus making it an adaptive estimator.
The study presents a moving horizon estimation approach for state and parameter estimations for a class of nonlinear hyperbolic system with application to open channel flow. The measurements are two fluid levels in the channel. The estimations converge to their actual values in finite time. Furthermore, the estimator is validated using experimental data. The main application of the method proposed in this paper is the real-time estimation of the return fluid flow rate during an oil well drilling process. The idea is to use a Venturi channel in the return flowline to estimate the flow rate and to aid early kick and loss detection during oil well drilling. For more details of the use of a Venturi channel for oil well drilling, refer to ; . These channels are usually placed on offshore drilling platforms which are subjected to movements due to ocean waves. This could not only generate ripples and waves in the channel which could give false level readings but also could indicate a false channel bottom slope. The reduced order model is sensitive to the noise in the level measurements and to the channel bed slope , and this might affect the estimation significantly. Therefore, the effect of waves via change of the bottom slope of the channel on the estimation is also studied.
The paper is organized as follows. The mathematical model and the estimator for the return flow rate are stated in detail in section 2. Then the experimental set up is summarized in section 3, followed by detailed results and discussion in section 4. Finally, the conclusions drawn from the results and discussion are summarized in section 5.

Mathematical Model for Fluid Flow
Hydraulics models are used with the objective of using a simple model, which removes unnecessary complexities of the use of multiphase flow models including the drill cuttings or gases. This is based on the fact that as far as the differential volumetric flow in and out of the well is concerned, a simple hydraulic model is accurate enough for estimating the net volumetric out flow rate.
from Navier-Stokes equations. There are different ways of expressing these models based on the physical natures those are assumed upon (Chaudhry, 2008). For shallow water flows i.e. the horizontal length scale is much greater than the vertical length scale, the shallow water equations are used.

Full Order Model
The shallow water equations for a 1D, unsteady, open channel system, is expressed as follows (Chow, 1959;Chaudhry, 2008), Here U is a vector of conserved variables A and Q, where A(x, h, t) is the wetted cross-sectional area normal to the flow and Q(x, t) is the volumetric flow rate. Further, h(x, t) is the depth of the flow, x is the distance along the flow direction and t is the time. F , the conservative flux vector consists of the force terms and S is the source term. A Venturi channel is a non-prismatic channel where the cross-sectional area changes with the flow direction (see Fig. 1). For a non-prismatic channel with a trapezoidal cross-section, these U, F and S terms can be expressed as follows , Here, q is the lateral inflow or outflow rate (assumed zero), g is the gravitational acceleration and β is known as the momentum correction coefficient or the Boussinesq coefficient which corresponds to the deviations of the local velocity over the mean velocity of the flow. W (x) is the bottom width of the channel and γ is the angle (in radian) that the channel bed makes with the horizontal axis, where it is considered positive when sloping downwards. I 1 is the first moment of area which represents the hydrostatic pressure term, and is expressed as follows for a channel with an isosceles trapezoidal cross-section, where S s is the side slope of the channel. T f is the frictional stress over the channel solid surface in the channel cross section. Generally, open channel models of water and other Newtonian fluid types have been developed throughout history (Chow, 1959;French, 1985;Chaudhry, 2008). Due to the low viscosity of these fluids, usually, the laminar flow conditions are applicable to turbulent flow conditions as well. The method of calculating the velocity based on the hydraulic grade line developed by Manning in 1890 is still accurate and well used in this regard (Chow, 1959). This method (also known as Gauckler-Manning model) is used in dynamic open channel flow models today, as a friction term. This friction term T f for a Newtonian fluid, when written using the Gauckler-Manning friction model is as follows (Chow, 1959;Chaudhry, 2008), Here n M is the Manning's friction coefficient, P is the wetted perimeter, andū is the average velocity of the flow. However, it is also duly noted that many other fluids exhibit complex flow behaviors which cannot be explained by the Newtonian theory. Nonetheless, significant theoretical advances have been made in this field of non-Newtonian fluid flow behavior in open channels (Kozicki and Tiu, 1967;Haldenwang et al., 2010Haldenwang et al., , 2002Burger et al., 2015a). Mostly for the non-Newtonian fluids, the flow conditions and the frictional relationship are presented in the form of a Moody chart. These charts include a plot of Fanning friction factor f for different Reynolds numbers in different flow regimes such as laminar, transitional, or turbulent. Similarly, the friction term T f in (4) can be written using the Fanning friction factor f for a non-Newtonian fluid as follows (Burger et al., 2015a), However, there are a limited number of experimental and/or theoretical values for f , which can only be applied to specific open channel types (Haldenwang et al., , 2002Burger et al., 2015a). Therefore, the Fanning friction factor can be considered to be unknown and hence it needs to be estimated. More specifically in a drilling operation, the rheological properties of the return fluid can change with time depending on the type of formation that is being drilled. Thus, the value for friction factor cannot be considered as a constant. It needs to be estimated in real-time.
The state A in (2), can be replaced with the fluid level h from the relationship A = hW + h 2 S s , to obtain the state h. It is possible to expand the F term using the product rule for simplicity. However, it will lose the conservative form of the equation. Conservation of the states is important in this study, therefore, no simplification is done here.
This system of equations can be solved by specifying the initial conditions and the boundary conditions. The specification of boundary conditions at either upstream or downstream end of the channel (or at both) depend on the characteristics of the fluid flow. For a super-critical flow, two boundary conditions need to be selected at the upstream of the flow, whereas for a sub-critical flow, two boundaries can be selected by choosing one at downstream and one at upstream. For trans-critical 1D flows, one boundary condition is used at the upstream of the flow (Chow, 1959;Chaudhry, 2008).
The geometry of the channel that is used for the study is shown in Figure 1. The placement of the Venturi section (the narrowest section also known as the throat) in the channel ensures that the flow conditions obtain critical or super-critical conditions at the throat, thus always having trans-critical conditions in the channel. Therefore, the upstream fluid level is taken as the boundary condition for solving the equations.

Reduced Order Model
The shallow water equations are nonlinear hyperbolic system of partial differential equations. These are commonly regarded as complicated systems to solve due to non-smooth solutions and possible discontinuities (Kurganov, 2018). In order to obtain the solutions that are fast enough with less computational resources to aid the online estimation, these partial differential equations are then converted into ordinary differential equations by the orthogonal collocation method. Here, the states h and Q are approximated by Lagrange interpolating polynomial at specific spatial positions along the channel which are known as collocation points. The positions of these collocation points are selected using an orthogonal polynomial, where the points lie at the roots of this polynomial. These points are nonequispaced throughout the domain and tend to cluster towards the edges of the interval when the degree of the interpolating polynomial gets higher. This will necessarily reduce the oscillatory problems given by the Runge's phenomenon. The Legendre functions of the first kind are selected for this study due to the low numerical oscillations given by the Legendre functions. For this study, two collocation points are used for simplicity. The model reduction method is described in detail in (Jinasena et al., , 2017.
The reduced order equations can be written as follows,ḣ anḋ for 2 collocation points (i = 1, 2). Here, subscript i represents the corresponding variable at the i th specific position along the channel. The first collocation point is the place where the level sensor LT 1 is placed for measuring the level h 1 as shown in Figure 1. The second collocation point is at the middle of the throat section where level sensor LT 2 is placed for measuring level h 2 . L ij refers to an element in L matrix at the i th row and j th column, where, Here, l is the length of the channel.

State and Parameter Estimation
The sensitivity of different terms of a similar model has been studied by . The two collocation point model is highly sensitive to the bottom slope, which can slightly change due to vibrations.
It is also sensitive to other factors as well and most of these are geometrical constants which can be measured properly. The level measurements are noisy and the friction factor is considered to be unknown and not measured. The flow rate through the channel is considered to be unknown. The main idea is to use the two level measurements to estimate both the flow rate and the friction factor using a suitable estimator. In reality, the fluid levels should always be positive values and the friction factor should have finite positive limits. This can be achieved in the estimator by applying constraints on the states and parameters, which is possible in a moving horizon estimator (MHE). Therefore, a MHE is selected in this study.
The system needs upstream level measurement as a boundary condition. Therefore, h 1 is taken as an input to the system and h 2 , Q 1 and Q 2 are considered as states while h 2 is measured. The system of equations for the estimation can be written in discrete time form as follows, where Here, x = [h 2 Q 1 Q 2 ] T , u = h 1 and y = h 2 . The function f n represents the nonlinear model given by (8) and (9) forḣ 2 ,Q 1 andQ 2 , respectively. The θ is n M and f for Newtonian and non-Newtonian fluids, respectively. k is the discrete time index. The measurement noise in the measured output at t k is denoted by v k ∈ R 1 , where v k ∼ N (0, R). Similarly, w k ∈ R 3 , where w k ∼ N (0, Q) accounts for process noise. Here, R and Q are the covariance matrices for measurement noise and process noise, respectively. In order to estimate the states using the MHE approach, we start by formulating an optimization statement using the Bayesian maximum a posteriori criterion for states as follows (Haseltine and Rawlings, 2005), Here, the idea is to compute the most likely values of the states {x 0 , . . . , x T }, given the measurements {y 0 , . . . , y T }. Assuming that both noise terms w k and v k are white Gaussian noise, the maximum a posteriori criterion of (14) can be written as, Here R k and Q k are the covariance matrices for the measurement noise and the process noise, respectively. P 0 represents the uncertainty in the initial estimate of x 0 .
Using (11), where x k+1 − f n (x k , t k , u k , θ k ) = w k , the variable w k can be introduced to (15) as follows, In a compact form, (16) can be written as, which is subjected to dynamics and constraints, for positive-definite functions L k and Γ (Rao et al., 2003). All the states {x 0 , . . . , x T } are to be found using all the measurements on the time interval [0, T ], which is known as the full information estimation. However, the main drawback is that the number of decision variables grows linearly with time T , which becomes computationally intractable for continuous processes. To reduce the computational cost with growing T , a finite moving horizon consisting of the most recent N measurements is considered, instead of the full information estimation. With a finite horizon backward in time of N samples and for T > N , (17) can be written as follows, The state at the start of the horizon x T −N and the process noises {w T −N , . . . , w T −1 } are the variables to optimize. The term Γ(x T −N ) is called the arrival cost, and it accounts for all the past data and the process dynamics from k = 0 to k = T − N − 1 (Rawlings et al., 2017). This arrival cost expression is written as follows, The arrival cost is often hard to compute, especially for nonlinear systems. Therefore, approximations such as using a uniform prior, extended Kalman filter covariance formula and MHE smoothing are used to express the arrival cost (Rao et al., 2003). The MHE estimation can also be used to estimate the unknown parameters of the process together with the states. The estimations of the states are adapted to the unknown parameters by continuously estimating both. In this sense, the MHE becomes an adaptive estimator. To estimate the parameters, these are assumed to be constant throughout the horizon. This reduces the number of variables to optimize. However, it is worth mentioning that the parameters themselves can change with time, and are considered constant only within the horizon at any given time. For adaptive MHE, the expression in (18) can be modified as follows, subject to: The third term in the objective function is the estimation error for the parameters, with P θ being the covariance matrix for the parameter noise or uncertainty. The last term is the regularization term for the parameters, with S θ being the weighting matrix for the parameters. Here, superscripts lo and up denote the lower and upper limit of the constraints, respectively. The nonlinear optimization problem for the MHE is solved using fmincon 1 in optimization toolbox of MATLAB2018b to obtain the estimatesx T −N , θ T −N andŵ k for k = {T −N, . . . , T −1}. The optimization algorithm 2 named as 'sqp' 3 in MATLAB is used in the function fmincon. The algorithm is based on sequential quadratic programming and satisfies all constraints at each iteration (Nocedal and Wright, 2006). Then using these estimated values as the initial condition, the past inputs {u T −N , . . . , u T } and the model (11) and (12), the current estimate for the states is calculated. Theta is considered to be constant throughout the horizon.

Experimental Setup
An experimental rig of an open Venturi channel with two level measurements and a pump flow rate measure-   ment is available at the University of South-Eastern Norway. The schematics of the flow loop is shown in Figure 2. The flow rate is measured before it enters the channel using a Coriolis mass flow meter (Promass 63F-uncertainty ± 0.1%). Two level sensors are placed in the Venturi channel (LT 1 and LT 2 in Figure 1). The level sensors are either radar sensors (Krohne Optiwave 7500 and 8300C) with ± 2 mm accuracy or ultrasonic sensors (Rosemount 3108) with ± 2.5 mm accuracy. The channel at operation is shown in Figure 3. Venturi channels are designed to have a critical flow condition at the throat. The criticality is determined by the Froude number which is the ration of inertia to gravity. When Fr < 1, flow is sub-critical, when Fr > 1, flow is super-critical and when Fr = 1, flow is critical. In order to estimate the flow rate using the fluid level, the upstream of the Venturi channel must be sub-critical. This can be achieved by keeping the channel at a lower angle. Therefore, the channel was kept horizontal throughout the experiments. Further, the changes of the bottom slope due to vibrations and

Results and Discussion
The main reason for the selection of MHE and the estimation of the fluid flow rate using MHE is discussed in this section. The effects of bottom slope and the friction factor are also discussed. The details of the simulation and experimental validation are stated below in separate sections.

Comparison with Other Estimators
Jinasena and Sharma (2018) tested a few different estimators for state estimation using the reduced order model shown in (8) and (9) with a Newtonian friction model instead of the non-Newtonian friction model shown in (7). Here, they tested a linear observer, a linear Kalman filter, an extended Kalman filter (EKF) and an unscented Kalman filter (UKF). The simulation results for the estimations of flow rate of this study are shown in Figure 4. The linear estimators were not suitable as evident from Figure 4. Therefore, they tested the other two nonlinear estimators, EKF and UKF, against the experimental data, and concluded that both the EKF and UKF are suitable for estimation of flow rate.
However, these estimators do not allow constraints. Since the full order model fails at zero fluid levels (zero flow rates), being able to have constraints is a necessity to obtain a full operating range of flow rates. Thus the MHE is selected to overcome this problem.

Simulation Results
The results and discussion of the simulations are stated here under different scenarios.

Effect of Friction Factor
In many different applications, the friction factor can change with time. For example, during an oil well drilling process, the rheological properties of the return fluid may change depending on the type of formation being drilled, and hence the friction factor also changes with time. The effect of the changes in the friction factor on the flow rate estimation is studied by comparing the results using a constant and a varied friction fac- tor. Here, the reference model is simulated using a varied Fanning friction factor for a non-Newtonian fluid. Then the estimation was done using the MHE without parameter estimation i.e. with constant friction factor (see Figure 5). Both the reference model and the estimator are excited with the same input level measurement h 1 . The volumetric flow rates obtained from the reference model and the flow rate estimated by the MHE are plotted together in Figure 6. The offset between the reference flow rate and estimator results are clearly evident. The estimation error of the flow rate when using a constant friction factor is clearly seen in Figure 7. The results are only shown for the non-Newtonian fluid model using (7). The results for the Newtonian model using (6) are similar. Therefore, it is essential to estimate the friction factor together with the states in order to obtain an improved flow rate estimation in a practical application because the fluid flow conditions change over time.

Estimation of Flow Rate and Friction Factor
The estimation of fluid flow rate and the friction factor was done using the MHE, where a horizon of 10 samples was considered. The results are shown in Figure 8. Here, the number of samples in a horizon is selected based on the accuracy of the estimation and the computational time. When the number of samples are high, the accuracy is high but the calculation is computationally heavy. Therefore, a trade off between these two is required, hence the number of samples was selected as 10, using trial and error method. Only the non-Newtonian model results are shown here in terms of the Fanning friction factor. Similar results are obtained for the Newtonian model results. The root mean square error (RMSE) value for the estimation error is 0.75 l min . When compared to the estimation error without adaptation to the friction factor in Figure 7, the online estimation of the friction factor has vastly improved the flow rate estimation.
Usually, MHE is known to be computationally heavy, due to the need for solving a nonlinear optimization problem at each time step. However, due to the reduced order model used in this application, the MHE can still be used in real-time. This is clearly seen in the CPU time plot in (d) in Figure 8, where the time taken for estimation is significantly lower compared to the sampling time of 0.5 seconds.

Effect of Channel Bottom Slope
In many different applications, the Venturi channels used for flow rate measurement are subjected to external disturbances such as channel bed vibrations and movements. This could generate false fluid level measurements due to waves in the channel flow and a false channel bed slope due to vibrations. This effect is highly significant in non-stationary process plants such as in ships or in offshore oil drilling platforms where the platforms are subjected to ocean waves. The Venturi channel is subjected to pitch, roll and heave motions in these cases. Since the model is sensitive to the changes in the bottom slope of the channel and the noise in the level measurements, the effect of change of bottom slope on fluid flow estimation is investigated here. The reference model is subjected to random changes of bottom slope (standard deviation .5 • ), and the level measurement under these changed conditions are used to estimate the flow rate and the friction factor. The affected level measurements due to the change of slope, the corresponding estimated flow rate and the friction factor are shown in Figure 9.
The figures illustrate that even though the level mea- surements have become noisy due to the changes in the channel slope, the flow rate through the channel can still be accurately estimated using MHE. The RMSE for the fluid flow rate is 5.86 l min . Therefore, a small change of bottom slope does not diminish the ability to accurately estimate the fluid flow rate. This further suggests that the Venturi channel could potentially be used in a moving platform under limited pitch conditions.

Experimental Results
The measurements which were taken during the experiments at the rig for the bottom slope angle are shown in Figure 10. This shows that at the experimental rig the channel bed is indeed subjected to vibrations. Due to vibrations small waves in the channel appear and contribute to the noise in the level measurements. The upstream fluid level measurement from the experimental data is taken as the input to the system and both the flow rates and the friction parameter are estimated using the MHE for two different fluid types. The level measurements for both the fluid types are shown in Figure 11. The ultrasonic measurements are noisier than the radar measurements.
The estimation for the flow rate of water through the channel is done using the Newtonian friction model given by (6). The estimated flow rate, estimated Manning's friction factor and the estimation error are shown in Figure 12, Figure 13 and Figure 14, respectively. Experimental results show that the flow rate can be estimated reasonably well. The RMSE value for the estimation error is 15.47 l min . Further, the estimated Manning's friction factor is within the range for the channel roughness type (Chow, 1959;French, 1985).
Similarly, the estimations for the non-Newtonian drilling fluid are done using the non-Newtonian friction model given by (7). The estimated flow rate, estimated Fanning friction factor and the estimation error for non-Newtonian drilling fluid are shown in Figure 15, Figure 16 and Figure 17, respectively. As shown in the figures, the fluid flow rate estimation using MHE is reasonably accurate. The RMSE value for the estimation error is 4.55 l min . The reduced order model used in the study is based on the simplest selection of the number of collocation points due to the simplicity. However, it could reduce the accuracy due to the approximation by a polynomial of the lowest possible degree (linear interpolation).  have shown that the least number of collocation points for open loop simulations of this model is three. However, the results of this study with MHE show that the use of two points is also enough when used for closed loop simulations i.e. with an estimator with the feedback information. This suggests that the real-time estimation of fluid flow rate and the friction factor using a simple reduced order model with a MHE can be potentially used in open channel flow applications.

Conclusions
A reduced order model is used for a moving horizon estimation for an open Venturi flow. Using two on-line level measurements in the channel, both the fluid flow rate and the friction factor is estimated in real-time. Adapting the state estimation to the friction factor has improved the fluid flow rate estimation, significantly. Further, two types of friction models with two different friction factors have been used based on the fluid type; i.e. for Newtonian fluids and non-Newtonian fluids.
The effect of change of bottom slope of the channel when the channel is placed on offshore or moving platforms is also investigated. The results show that the estimator is capable of handling the noise and the errors created by the false fluid levels due to variations in the channel slope. The results are validated by laboratory scale experimental results for two types of fluids, one Newtonian and one non-Newtonian fluid. The results suggest that this method can be used for real-time fluid flow rate estimation for open channel flows. In addition, it also has the potential to be used in offshore drilling processes for online estimation of return flow of the drilling fluid.