Uncertainty Analysis of a Simpliﬁed 2D Control-relevant Oil Reservoir Model

In this paper, a simpliﬁed 2D control relevant model for a slightly slanting wedge-shaped black oil reservoir is made more realistic by incorporating model uncertainty. The uncertainty in the model is computed via Monte Carlo simulation. Furthermore, based on this model with uncertainty, a Proportional + Integral (PI) controller is implemented to increase oil production while minimizing water production. A PI controller is used to control the valve opening of the Inlet Control Valves (ICVs) in the production well. Implementation of a PI controller enhanced the oil recovery in 1000days by 1 . 79%, while the total water production is reduced by 2 . 59%.


Background
Norway is one of the leading suppliers of oil and gas to the global market. Revenues from sales of oil and gas have played a vital role in creating modern Norwegian society. Oil and gas are trapped in the subsurface formation of relatively thin slabs of porous rock. Oil wells are drilled into the subsurface with an oil rig to extract the oil and gas from the reservoir. The production of oil can be increased by predicting and managing the future performance of the oil reservoir. However, because of the subsurface complexity and limited data, numerous uncertainties are present in oil reservoir characterization. These uncertainties should be considered for better future prediction of oil reservoir performance. In project no. 308817, "DigiWell", of the Research Council of Norway, it is of interest to combine reservoir models and well transport models under uncertainty; these models operate under very different time scales, which poses a numeric problem. It is of interest to formulate simplified models for conceptual studies. A potential reservoir model in such a conceptual study could be a 2D, wedge-shaped black-oil model.

Previous Work
A number of simulation tools exist for prediction of oil reservoir performance, e.g., ECLIPSE, MRST, INTER-SECT, MEERA, OLGA ROCKX 1 . Most of these tools are commercial; a few of them support simulation under model uncertainty. In project "DigiWell", tool MRST will be a working tool. Zolotukhin and Ursin (2000) give an introduction to petroleum production, how to find experimental data/model parameters from laboratory analysis, and indicate basic model formulation. Chen et al. (2006a) focus more on general model formulation. Lie (2019) presents the modeling framework used in simulation tool MRST. Zhang (2013) developed a simplified 2D, control relevant model of a slightly slanting, 2D wedge-shaped black oil reservoir. The model was implemented in MATLAB, using fixed step-length Explicit Euler discretization. A PI controller was used to control the valve opening of the inflow control valves (ICVs) in the production well. Minimum water saturation S w over the reservoir was taken as a set-point for the PI controller. Because this minimum saturation is not available from measurements, a different approach is needed for a realistic solution. Bhattarai (2021) re-formulated the model from (Zhang, 2013) and implemented the model in computer language Julia, using the DifferentialEquations.jl package with variable step-length solver Tsit5(). Different saturation vs. relative permeability correlations were used, and simple parameter uncertainty was introduced, allowing for Monte Carlos simulations posed as EnsembleProblem in Julia. A more realistic PI controller that reduces the water cut (WC) was implemented. This work is presented here. In Section 2, the simulation model is developed. Section 3 presents model uncertainty and the PI controller. Section 4 provides simulation results. Finally, some conclusions are drawn in Section 5.

Two-phase Flow in a Porous Media
A black oil reservoir model has a water component, as well as hydrocarbon components divided into a gas component and an oil component with no mass transfer between the water phase and the other two phases (oil and gas) (Chen et al., 2006a,b). We further simplify the black oil model by considering a relatively new heavy oil reservoir without gas.

Reservoir Overview
The model in this work is developed for a slightly slanting wedge-shaped horizontal black oil reservoir with homogeneous dispersion of water and homogeneous in its geological features permeabilities, porosities, etc., Figure 1.
A natural aquifer with constant pressure P a and constant relative permeability k wa is located at the bottom of the reservoir, while a horizontal well is located at the top of the wedge-shaped reservoir. We use coordinates length ∈ [0, L] along the well, and radius r ∈ [0, R] from the well towards the aquifer. The boundary conditions are zero flux at ( = 0, r), ( = L, r), and at ( , r) for the slanting angles θ = α ± β 2 . A water flux exists at ( , r = R) and a flux of oil and water mixture at ( , r = 0) into the well. The total production rate of the mixture of oil and water from the reservoir q mix,tot is specified to be constant. The spatialtemporal variables are represented by (r, ), and t, respectively.
The wedge-shaped reservoir model is combined with the well model. From the combined reservoir and well model, we are interested in finding how the water saturation S w , reservoir pressure P, water and oil volumetric production rates (q w,s , and q o,s , respectively), well pressure P rw , well bottom hole pressure P bh , etc., vary with time.

Reservoir Model
With bulk volume V b and pore volume V p of the reservoir, porosity is given as φ ∂V p /∂V b . Fluid saturation S η for fluid η ∈ {o, w} (oil, water) is defined as S η ∂V η /∂V p where V η is the volume taken up by fluid η. It follows that ∑ η S η = 1. For two-phase flow (Zhang, 2013) (1) The velocity terms in Equation 1 are given by Darcy's law: where, λ η = k rel,η µ η K, K is the absolute permeability, k rel,η is a relative permeability for phase η, µ η is the viscosity of phase η, P η is the fluid pressure, and g is the acceleration of gravity.
The z-term has to be projected to the r-coordinate according to Figure 1 Here, angle α is the slope of the wedge, while angle β is the angular width of the wedge, Figure 1.

Well Model
According to Chen and Zhang, the pressure close to the well declines much faster than near the aquifer. Therefore, a small step size ∆r near the well is required for accurate pressure calculation in the reservoir cell at the neighborhood of the well. This can be handled by using local grids refinement in the neighborhood of the well. However, this can lead to restrictions on time steps in the numerical simulation (Chen et al., 2006a). The alternative solution is to derive an analytical solution for the steady-state flow model that yields the Peaceman equation (Peaceman, 1993). As-suming only radial flow in grids near the well, (4) Total specified flow rate q η,s,tot is the sum of the flow rates from all perforated zones where, N v is the total number of perforated zones of the well, and (r, ) Specified total oil and water mixture production rate q mix,tot can be written as where, q w,s,tot and q o,s,tot are the total production rate of water, and oil, respectively.  Figure 2. The block-centered grid system and a five-point stencil scheme (Zhang, 2013).
Total oil, and water mixture production rate for each cell near the well q mix,r, can be written as q mix,r, = q w,s,r, + q o,s,r, .

Simplifying Assumptions
To simplify the implementation of the model, it is convenient to number the grids by i r and i , where we number the grid for (i r , i ) ∈ {1, . . . , n r } × {1, . . . , n } as shown in Figure 2.
The following simplifying assumptions are made: 1. Both rock and fluid are incompressible, leading to constant density and formation volume factor.
3. Capillary pressure is assumed to be equal to zero, i.e., P cow = P o − P w = 0. Here, P o and P w are pressure exerted by oil and water, respectively.
4. Effect of temperature is neglected.
6. Isotropic medium K r = K With these assumptions, and introducing the definition Equation 1 can be simplified to The well model is,

Valve and Pipe
The valves are represented by a homogeneous flow model of sub-critical flow through a pipe containing restriction as where and here, C u is a unit conversion constant, C v is a dimensionless flow coefficient of the valve, A valve is the constriction effective area, ρ mix is the density of the fluid mixture, q mix is the volumetric flow rate of the mixture, P rw is the well pressure, and P bh is the bottom hole pressure. The pipes are modeled as a hydraulic network using the following equation: with where (i r , i ) ∈ n r × {1, . . . , n well }, f is the fanning friction factor, ∆L pipe is the pipe step length which assumed to be equal to well step length ∆ well and reservoir step length ∆ , r p is the radius of a production pipe.

Water Saturation Versus Relative Permeability
The data for water saturation S w and corresponding relative permeabilities are typical values provided by the industry. This is achieved through least-square fit using a standard package in Julia called LsqFit 2 . Comparison of water saturation and permeability relation for actual data and least-square function is shown in Figure 3

Mobility Determination
Most of the works in the literature use an upstream scheme to evaluate the mobilities (Cordazzo et al.).The mobility λ at the integration point is evaluated upstream of the flow.  . Relation between water saturation S w and relative permeabilities.

Numerical Solution
In the numerical simulation, the finite difference method uses finite differences to approximate derivatives of ordinary differential equations. The forward difference method is used in this work.
φ dS η dt i r ,i = a η,1,i r ,i P i r +1,i + a η,2,i r ,i P i r ,i −1 + a η,3,i r ,i P i r ,i + a η,4,i r ,i P i r ,i +1 + a η,5,i r ,i P i r −1,i − a η,6,i r ,i , where,

Pressure Equation
To derive the pressure equation, we sum Equation 16 for the oil and water phases, and using the fact that where we set a w,1 + a o,1 = a 1 a w,2 + a o,2 = a 2 a w,3 + a o,3 = a 3 a w,4 + a o,4 = a 4 (19) a w,5 + a o,5 = a 5 a w,6 + a o,6 = a 6 .
This leads to a 6,i r ,i = a 1,i r ,i P i r +1,i + a 2,i r ,i P i r ,i −1 + a 3,i r ,i P i r ,i + a 4,i r ,i P i r ,i +1 + a 5,i r ,i P i r −1,i The pressure Equation 20 can be written in matrix form where A is a five-diagonal sparse matrix, and P is a vector of unknown pressures, The pressure is solved from this implicit, linear equation. 3 Model Uncertainty and PI Controller

Uncertainty Analysis
Some of the methods developed for uncertainty in the field of the petroleum industry are experimental design, response surface, multiple realization tree, and Monte Carlo simulations. In this work, we use the Monte Carlo method, which implies simulating an ensemble of cases where uncertain parameters are drawn from some statistical distribution (Zhang, 2003). Uncertain reservoir parameters include porosity, saturation, and permeability.
We assume uniform distributions for porosity, permeability, and initial water saturation where these parameters are varied by ±20%, Table 1.
Monte Carlo simulation is performed using Ensem-bleProblem in Julia which interfaces well with the standard differential equation solving package DifferentialEquations.jl.

PI Controller
There are a number of inflow control devices (ICD; passive) with the flow that is contained in an inner pocket in the production pipe. Inflow control valves (ICV; active) give a controlled flow out of this inner pocket and into the production pipe. It is assumed that the ICVs are installed at every 60 m length along the oil production pipe, i.e., at each segment of the pipe length. The group of ICVs receives the same control signals. In other words, there are 20 ICVs installed at the production pipe in which the group of valves at n well = {1, . . . , 5}, n well = {6, . . . , 10}, n well = {11, . . . , 15}, and n well = {16, . . . , 20} receive the same control signals, respectively.
The standard methods for SISO PID controller tuning are Skogestad's method, Ziegler-Nichols method, the Good Gain method, etc. Because the reservoir system is a MIMO system with strong interactions, these standard methods fail, and we instead tune the PI controller manually through trial and error. The values of K p and T i are taken as 7.506 · 10 −6 , and 2.9376 · 10 9 , respectively, for all loops.

Simulation Results
The value of parameters used in the simulation of the model is shown in Table 2. Figure 4 shows 3D water saturation profiles after 50, 100, 300, 500, 800, and 1000 days of production, respectively. Water from the aquifer slowly advances towards the production well. When the water flooding front reaches the wellbore at r = 0, water breakthrough occurs. In Figure 5, it can be observed that the water saturation slowly starts to increase along the reservoir radius with increase in production time. In these water saturation profiles, the water coning effect is not visible because the volumetric flow rate of fluids in the well is assumed to be evenly distributed.
In the water saturation profile after 300 days of production, the water saturation at grids n r = 10 to n r = 20 are the same as the initial water saturation of the reservoir i.e., 0.15. However, after 500 days of production, we can observe that the water saturation at grids n r = 20 is higher than 0.15. This is because the water breakthrough already occurred after 390 days of production. After 1000 days of water production, the average water saturation at n r = 20 is observed to be 0.249. Figures 6 and 7 demonstrate the effect of PI controller in total production of oil and water in 1000 days, respectively. In these figures, we can see that the total oil production, after implementing the PI controller, is increased by 1.79% while the total water production is decreased by 2.59%. Figure 8 demonstrates the total volumetric flow rate of oil and water after implementing the PI controller. The average oil volumetric flow rate per day is increased by  Step length along reservoir radius L n r m ∆ Step length along reservoir length L n m ∆ well Step length along well length L well n well m Figure 5. Saturation profile for model with PI controller.  1.738%, while the average water volumetric flow rate per day is decreased by 2.7715%. Similarly, Figure 9 shows the effect of PI controller on water cut at = 0.5L or n = 10. The average water cut is decreased by 2.755% after implementing the PI controller. Finally, in Figure 10, we can observe that the bottom hole pressure tends to decrease after the use of the PI controller. The PI controller is implemented in a model with uncertainties. The comparison for the total oil and water production volume after implementing PI controller is shown in Figures 11 and 12, respectively. In Figure 11, the solid lines with maroon color represent uncertainties in total oil production volume without PI controller, while the green dot lines represent uncertainties in total oil production with PI controller. In this figure, we can observe that the minimum total oil production volume in 1000 days after implementing the PI controller is 359.0391 · 10 3 m 3 which is higher than that of the model without PI controller i.e., 326.4863 · 10 3 m 3 . However, the maximum total oil production volume in 1000 days after implementing the PI controller is 614.6868 · 10 3 m 3 which is lower than that of the model with PI controller i.e., 657.2351 · 10 3 m 3 .
Similarly, In Figure 12, blue solid lines represent uncertainties in total water production volume without PI con-   troller, while the red dot lines represent uncertainties in total water production with PI controller. In this figure, the maximum total water production volume in 1000 days after implementing the PI controller is 440.9470 · 10 3 m 3 which is lower than that of the model without PI controller i.e., 473.5101 · 10 3 m 3 . However, the minimum total water production volume in 1000 days after implementing the PI controller is 185.3096 · 10 3 m 3 which is higher than that of the model without PI controller i.e., 142.7613 · 10 3 m 3 .

Conclusions
An overview of an oil well holding black oil, with a reservoir, and pipes is given. A simplified 2 D control-relevant model is developed and implemented in Julia programming language, and solved using an efficient, variable time step method. Results such as total production volume of oil and water, water cut, volumetric flow rates of oil and water, bottom hole pressure etc., are compared to those in Zhang (2013). Even when the :parameter values used in this work differ from those in (Zhang, 2013), the results are qualitatively similar to those of (Zhang, 2013). Further-more, uncertainties in the model are discussed and studied by Monte Carlo simulation. Finally, a PI controller is im-plemented in the model based on uncertainties to enhance the oil production, while minimizing the water production. The PI controller helps to increase the oil production by manipulating the ICVs at the production well. Implemen-tation of PI controller improved the total oil production in 1000 days by 1.79%. However, the effect is not very sig-nificant due to the limited capability of a PI controller. In this case, a more effective controller is required, such as MPC.