On Uncertainty Analysis of the Rate Controlled Production (RCP) Model

Rate controlled production (RCP) model is used to simulate and investigate the performance of the oil wells which are completed by autonomous inflow control devices. In order to quantify the performance of the RCP model, a dimensionless version of the model is considered, and its parameters are estimated. We demonstrate how the model and the measurement uncertainties can be quantified within the Bayesian statistical inference framework. In this relation, Hamilton Monte Carlo (HMC) is used to draw samples from the joint posterior probability distribution. We demonstrate that at the calibration step the modified model is able to capture the variations in the measurements. However, the cross- validation with the new data has revealed that the modified model tends to overpredict the pressure drop. This inadequacy cannot be explained by the measurement noise or the uncertainty in the estimated parameters. These results also imply that the original RCP model needs revision.


Introduction
Increase in oil production and recovery have been always the main objective of the oil industry. Hence, different methods and technologies have been developed to achieve this goal. One of the proven methods is to drill long horizontal wells, which increases the reservoir contact and consequently makes the oil production feasible and more economical.
However, long horizontal wells are likely to experience more pressure differences between the heel and toe section. This leads to non-uniform flow and consequently breakthrough of unwanted fluids in the heel section of the well, as shown in Figure 1.This phenomenon is known as heel to toe effect .
Autonomous Inflow Control Valve (AICV) together with Autonomous Inflow Control Devices (AICD's) like RCP valves are among the newest technologies that have been developed for Increased Oil Recovery (IOR). By balancing reservoir drawdown, these valves delay the onset of water and/or gas breakthrough and in case of breakthrough, it will restrict the production of these unwanted fluids significantly. A mathematical model describing the performance of the RCP valve was originally developed by Mathiesen et. al. in 2011(Mathiesen, et al., 2011.This model is later being used to describe the AICV performance too. In recent years, both lab and production data from various oil wells have been used to check the validity of the model (Mohd Ismail, et al., 2018;Langaas, et al., 2020). This model has been implemented in reservoir simulators such as NETool and Eclipse in order to simulate the performance of the valve under static and dynamic conditions.
In order to be able to employ the model, one needs to estimate the model parameters prior to its deployment. It appears that one of the methods used by many practitioners for parameter estimation prior to utilization of the model in NETool is the trial-and-error method Halvorsen, et al., 2016). Nevertheless, if one assumes that the model is correct, in most practical cases, the classical least square or similar methods are sufficient to produce good estimates for the model parameters (Moradi, et al., 2021). However, there are some evidences that the model does not explain all the variations in the data (Langaas, et al., 2020). There has also been attempts to modify the model (Voll, et al., 2014).
In order to be able to verify model inadequacy, two pre-conditions are needed to be satisfied. The first one is accurate and precis measurements of the valve behaviour and the second one is the quantification of the different sources of the uncertainty. In this short paper, we will demonstrate how results from accurate measurements can be used within the Bayesian statistical inference framework to quantify and model the sources of the uncertainty and check how good the model explains the variations in the measurements.

AICV Principle
AICV utilizes viscosity and density differences between the reservoir fluids in such a way that it will keep the valve open for oil and closed for unwanted fluids like gas and water. Figure 2 illustrates AICV in open and closed position, respectively. This is achieved by taking advantage of the pressure differences in the Laminar Flow Element (LFE) and Turbulent Flow Element (TFE). These two flow restrictors are connected in series, which is illustrated in Figure 3 . AICV consists of two flow paths: the main flow path and the pilot flow path. Pilot flow path consists of two flow restrictors of LFE and TFE. When reservoir fluid enters the main path, a small portion of the flow is guided through the pilot flow, which is located near the main path. If a fluid with high viscosity enters the AICV, its flow through LFE will lead to a higher pressure drop over LFE. This phenomenon can be explained by Darcy-Weisbachs equation: where ΔP is the pressure drop. is the friction factor (64/Re) Re is Reynolds number.
is the fluid density. µ is the fluid viscosity.
is the fluid velocity. and are the length and diameter of the LFE respectively.
After passing through the LFE, which is a pipe segment, fluid enters a chamber. The second flow restrictor TFE, which is a nozzle, is placed in this chamber. The pressure drop across the TFE as described by Bernoulli, is calculated using the equation: in which, is a geometrical constant. Combination of these two flow restrictors results in a pressure drop, which determines how the AICV functions. As it is shown in Figure 3, high P2 will move the piston upwards closing the AICV for unwanted fluids while low P2 will keep the piston at its neutral position that maintains the oil production. The concept and principle of AICV is described in detail in earlier SPE papers (Taghavi, et al., 2019;.

RCP Model
The RCP model for the valve can be described as: where ∆ is the differential pressure across the AICV and µ are the calibration fluid density and viscosity, and and µ are the mixture fluid density and viscosity. The parameter is a valve characteristic given by the ICD strength, is the volumetric mixture flow rate, and and are constants (Mathiesen, et al., 2011).
In order to reduce the complexity in this short article, we will concentrate our efforts on a single-phase oil flow. In addition, the model will be evaluated for three types of oil with different densities and viscosities.
The model described by Eq.
(3), is dimensionally inconsistent. In order to avoid handling this inconsistency and its consequences, we study the flow rate vs. pressure drop with respect to a reference fluid at the same temperature. Therefore, we have chosen water at 20 degrees and a flow rate around 120 l/h. The measured pressure drop for water under these conditions is around 10 bar. Consequently, since is a geometric parameter and hence independent of the fluid type, it will not play a role in the analysis. Then from Eq. (3) follows that the relative pressure drop with respect to water is As it was mentioned earlier, there are some indications that the RCP model does not explain all the variations in the data. For this reason, we propose to use a multiplicative noise term in order to quantify possible model discrepancies. The modified dimensionless RCP model is where denotes the multiplicative noise term. Since the relative pressure drop is positive, we assume that a priori is distributed according to Gamma distribution, with its mode at one. An -value very close to one is an indication that the model can adequately describe the variations in the data. The statistical inference will reveal the probable values of . In the following, Eq.
(5) along with the experimental data are used to estimate the parameters , and . These estimates are used to evaluate the performance of the modified RCP model.

Experimental Setup
The experiments were performed on the AICV prototype test rig at InflowControl's multiphase test facility located in Porsgrunn, Norway. A simplified schematic of the test rig showing the key elements of equipment and key measurement locations is shown in Figure 4. The tests can be carried out with water, pressurized air, and silicone oil as test fluid. The test facility is designed for single-and multiphase oil, water, and gas. A multistage centrifugal pump increases the water/oil pressure from the water/oil supply. Compressed air at room temperature can be regulated to the desired pressure for each case, up to maximum 200 bar. Flow rates, density and temperature are measured close to the inlet of the test vessel by a Coriolis flowmeter. A pressure transmitter measures the inlet pressure, whereas a differential pressure transmitter measures the differential pressure over the test vessel. Multiphase flow tests can be performed by injecting the desired oil flow rate to the test vessel, which is already filled with gas. The desired oil flow rate is injected from a separate test rig, which is connected to the singlephase test rig. The green dashed line in Figure 4 show the multiphase test flow path.

Test Conditions and Data
Single-phase flow tests were performed with silicone oil as test fluid. The system conditions such as temperature and pressure, flow rates, pressure drops over the AICV and fluid properties, such as viscosity and density are controlled and measured in each test. The data obtained during the tests are listed in Table 1 in the Appendix. Temperature, density, and mass flow rate were measured using a Coriolis flow meter and the differential pressure across the AICV were measured by using a high precision pressure transmitter. Viscosity was measured and calculated manually using an Ubbelohde type viscometer. Viscosity measurements were performed several times under stable conditions in order to minimize the uncertainties. The accuracy of the different measuring tools employed in the tests are listed in Table 2 in the Appendix.

Bayesian Inference
The calculus of Bayesian inference is based on the application of two rules, the product, and the sum rules of the probability theory. One of the useful forms of the ( | , , ) is the posterior distribution over the possible values of consistent with the measurements, the model and any other available and relevant background information denoted by ; like any information about the valve construction. On the righthand side of the above equation, ( | , , ) is the likelihood, which is a statement about how likely it is to measure D given the model and specific values for . The term ( | , ) is known as the prior distribution. In the present context, it models the expert opinion about the possible values of the . The last term ( | ) functions as the normalization constant and is independent of and hence not relevant in the present context. We remind the reader that the letters in the parentheses stand for logical propositions and "," denotes the logical "AND" operation. However, in calculations we work with algebraic expressions. The context will determine the use.
Often, as in the present case, the inference on also depends on some other parameters, for which neither their true values are known nor are they of primary interests. Nevertheless, due to dependency of inference on them, they must be part of the estimation process. These parameters are known as the nuisance parameters. Here the sum rule of the probability theory can be useful. Let denote the vector of the nuisance parameters, then by the sum rule we have ( | , , ) = ∫ ( , | , , ) Ω (7) The above operation is called marginalization. Basically, calculating the above integral is the same as averaging the integrant over all possible values of . Marginalization is a very powerful concept and will be used in the next section. The reader is referred to (Kruschke, 2015) for further reading on Bayesian inference.

Statement of the Inference
In the following, let the model parameters, the nuisance parameters, data, and the model be denoted by , , and , respectively. That is, The description of each symbol is listed in Table 3 in the Appendix. The main reason for the choice of the nuisance parameter vector , is that we are uncertain about the true values of these parameters. For example, even though we have taken great care in measuring the viscosity, there is no guaranty that the conditions under which the oil flows through the valve are exactly the same as the viscometer. Therefore, we have chosen to include as one of the nuisance parameters. Similar reasons are behind the choice of other components of . We emphasise that this is an important component in quantification of the sources of the uncertainty. The lack of knowledge about the true values of the parameters under different test conditions, which are not possible to be controlled during the experiments, constitute an important source of the uncertainty.
By the Bayes theorem, the joint posterior distribution is ( , | , , ) ∝ ( | , , , ) × ( , | , ). (8) The choice of the likelihood is determined by the measurements noise, while the choice of the prior distribution is based on the uncertainty in the expert knowledge about the true values of the parameters, before considering the measurements. For example, as was mentioned previously, represents the multiplicative noise. It is positive and we expect its value to be one. However, there are reasons to believe that the model tends to overestimate the pressure drop over the valve. Therefore, we suspect that there is a good chance that can attain values below one. For these reasons, we choose (2,2), the gamma distribution with the parameters (2,2), to represent our prior knowledge about . The expected value of this distribution is one and its mode is at one-half. However, after seeing the data, the posterior distribution of might be different, which as we shall see, is indeed the case. Note that (2,2) has non-zero mass for all > 0. That is, the prior distribution does not exclude any positive values of . It only makes some values less probable. The marginal posterior distribution of will allow the data to modify the belief represented by the prior. In a similar manner, the expert knowledge on the other parameters can be incorporated in the inference process through appropriate choice of the prior distributions for each parameter. We have summarized the choices of the priors and the likelihoods for each parameter in Figure  5. Due to logical independence between the parameters, the joint posterior distribution in Eq. (8) is the product of all the distributions listed in Figure 5. All the parameters are positive and in case of , it is larger than 2. This means that all the normal distributions are truncated at zero. In the case of , we have a truncated gamma distribution with lower limit being 2. The marginal posterior distribution is found by integrating over the domain of .

Markov Chain Monte Carlo Simulation
It is difficult to find an analytical expression for the joint-and the marginal posterior distributions of the parameters. This is generally a challenging task in Bayesian statistics. A common approach is to approximate the joint posterior distribution by large number of samples. The generation of samples are often conducted by a class of dependent sampling methods known as Markov Chain Monte Carlo (MCMC). Roughly explained, the method works by sampling the distribution relative to the height of the distribution function on its domain. The frequency distribution of these samples will on the long run converge to the true distribution. Computationally, one starts with a random sample and generates a chain of samples following certain sets of rules, which will guaranty that the chain will eventually visit all the regions relative to their probability mass. Since in practice one can only generate finite number of samples, it is important to check if the chain has found the regions of highest probability mass. There is a so-called burn-in period, below which all the samples are discarded. The reason for this is to make sure that in a set containing a finite number of samples, the samples from regions with low probability mass are not over-represented.
For the purpose of this study, we run a MCMC method known as Hamiltonian Monte Carlo (HMC), using the statistical software known as Stan, which comes also as a R package known as RStan (Stan development team, 2019). We have run four chains, each with different starting points. Figure 6 shows the output of the chains for each of the model parameters.
As it can be seen, regardless of the initial starting point of the chain, after a burn-in period of roughly 10K, all the chains are stabilized and converged. For more details, we refer the reader to (Kruschke, 2015).  The histogram density of the parameter reveals that the model is hugely over predicting the relative pressure drop over the valve. More specifically, the pressure drops over the valve have to be scaled down to 2.2%-3.2% of their predicted values by the model in order to be consistent with the measurements.

Calibration and Validation
The dataset D used in the MCMC simulation is generated by running experiments on two different oil types with viscosities 6.6 cP and 36.4 cP (see Table 4 in the Appendix). The measurements and the posterior samples from the MCMC with 99% credible intervals are plotted in Figure 8. Except for two points, for the case of 36.4 cP oil, all the pressure drops predicted by the model are within the 99% credible interval. That is, at the calibration step, the model can describe most of the variations in the measurements. At this stage, without further measurements, it is difficult to explain the reason(s) for the two borderline outliers observed in the dataset for 36.4 cP oil.
The validation is conducted on a new dataset, which was not used in the estimation of the model parameters. This second dataset is generated by running the experiment on an oil with viscosity 12.6 cP. For this, we need to find the posterior predictive distribution.
Indeed, let = (∆ , ) denote the unobserved new data. Then the posterior predictive distribution is defined as ( | , , ). In order to be able to use the model M, one needs to know the model parameters. Note that by the product rule, the integrant can be expressed as ( , θ, ω| , , ) = ( |θ, ω, , , ) × (θ, ω| , , ) The observant reader recognizes that the second term on the right-hand side is the joint posterior distribution defined by Eq. (8). The first term on the right-hand side is called the sampling distribution and its functional form is same as the likelihood. The difference is that unlike likelihood, which is a function of the model parameters, the sampling distribution is a function of and is normalized to unity over the domain of . By applying the following algorithm, one can generate samples from the posterior predictive distribution, 1. Generate (θ i , ω i ) from (θ, ω| , , ) 2. Generate from ( |θ , ω , , , ) 3.
The above algorithm is iterated a given number of times. The histogram of the generated samples can then be considered as an estimate for the posterior predictive distribution defined by Eq. (9). Note that we are already in disposition of the samples (θ i , ω i ). They are the samples generated from the joint posterior distribution during the calibration step. Thus, we only need to conduct the step 2 in the above algorithm. From the product rule, and the nature of measurements noise, follows that the sampling distribution can be expressed as product of two normal distributions In the above expressions is the given flow rate for which one seeks to calculate the corresponding pressure drop. The algorithm for generating samples from the sampling distribution can be formulated as follows 2.1. Generate from ( , ) 2.2. Generate ∆ from (∆ |M(θ, ω), , ) The result of the cross-validation with 99% credible error-bars is given in Figure 8. For low flow rates or equivalently low-pressure drops, the model prediction is within the 99% credible interval of the measurements. However, it appears that for high flow rates, the model has some tendency to over-predict the differential pressure over the valve.

Conclusions
In this paper, we demonstrated how the model and the measurement uncertainties can be quantified within the framework of the Bayesian statistics. In order to avoid complications due to dimensional inconsistency of the original model, we proposed a dimensionless version of the model. The result of our analysis revealed discrepancies, which could not be explained by the measurement noise or the uncertainty in the estimated parameters. The model inadequacy can be divided into global and local categories. The most serious problem observed was at the global level. Indeed, the predictions of the dimensionless model given by Eq. (5) had to be scaled down to 2.2%-3.2% of their values in order to be at the same level as the measurements. This has not been observed before or reported in literature. We believe that the main reason for this is that this type of scaling would in general be absorbed into the factor and hence would slip away unnoticed. Further studies are needed to determine the source(s) of this inconsistency. If one accepts the correction factor and hence the modified dimensionless model given by Eq. (5), the deviation at the local level is less significant. The model validation has revealed that there is a tendency for the modified model to over-predict the pressure drop. A closer study of the results has revealed that a slight increase in oil viscosity during its passage through the valve can explain most of the overestimated pressure drop tendencies by the model. Further studies under more stringent conditions will be conducted in order to uncover the causes of these observations.