- Email: [email protected]
Adaptive modeling of waterflooding process in oil reservoirs
Author’s Accepted Manuscript Adaptive Modeling of Waterflooding Process in Oil Reservoirs Farzad Hourfar, Behzad Moshiri, Karim Salahshoor, Mojtaba Zaare-Mehrjerdi, Peyman Pourafshary www.elsevier.com/locate/petrol
PII: DOI: Reference:
S0920-4105(16)30256-X http://dx.doi.org/10.1016/j.petrol.2016.06.038 PETROL3530
To appear in: Journal of Petroleum Science and Engineering Received date: 4 November 2015 Revised date: 28 April 2016 Accepted date: 22 June 2016 Cite this article as: Farzad Hourfar, Behzad Moshiri, Karim Salahshoor, Mojtaba Zaare-Mehrjerdi and Peyman Pourafshary, Adaptive Modeling of Waterflooding Process in Oil Reservoirs, Journal of Petroleum Science and Engineering, http://dx.doi.org/10.1016/j.petrol.2016.06.038 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Adaptive Modeling of Waterflooding Process in Oil Reservoirs Farzad Hourfara*, Behzad Moshiria, Karim Salahshoorb, Mojtaba Zaare-Mehrjerdia, Peyman Pourafsharyc1 a
b
School of Electrical & Computer Engineering, University of Tehran, Tehran, Iran Department of Automation & Instrumentation Engineering, Petroleum University of Technology, Ahwaz, Iran c School of Chemical Engineering, University of Tehran, Tehran, Iran *
Corresponding Author: Farzad Hourfar, Tel.: +98-912-571-5925. [email protected]
Abstract: Waterflooding is a common solution for enhancing oil recovery in oil reservoirs in which water is injected to the reservoir. For this purpose, injection rates should be determined to optimally control the desired oil production from each well to accomplish reservoir management strategies. This goal can be achieved by precisely modeling the waterflooding phenomena. In this paper, an applicable algorithm for modeling the waterflooding process in oil reservoirs has been developed. Using proper system identification techniques will help to consider any change in the inherent reservoir specifications in addition to operating point variations due to injection alternation rate and also severe nonlinear nature of the reservoir and consequently to obtain adaptive model for the reservoir which is able to reflect all the mentioned phenomena. The reservoir has been considered as a Multi-Input/Multi-Output (MIMO) plant with several manipulated variables and outputs to mimic the variables dynamics of the reservoir during its operational life. A Recursive Least Square (RLS) approach has been utilized with proper forgetting factor and sliding window length selection to extract the reservoir dynamic model in the content of a MIMO model structure. The developed approach has been evaluated for SPE10 model as a well-known reservoir case study in this paper. The observed results demonstrate that the proposed method has acceptable capability, compared to the conventional off-line Least Square (LS) and Extended Least Square (ELS) approaches which may lead to AutoRegressive model with eXogenous input (ARX) and AutoRegressiveMoving-Average model with eXogenous inputs (ARMAX) structures, for modeling the dynamical behaviors of the reservoir to be employed as a powerful tool for control and optimization applications. Keywords: Waterflooding Process; Recursive Least Square Method (RLS); Adaptive Reservoir Modeling; Multi-Input-Multi-Output(MIMO) Systems ; System Identification; Control & Optimization 1
Present Address: Department of Chemical & Petroleum Engineering Sultan Qaboos University, Oman
1
1. Introduction In order to fill the gap between energy demand and supply, it is necessary to enhance oil recovery from existing fields. Nowadays, the average for the recovery factor is about 35% (Rossi et al. 2002), which could decrease to about 15% for complex reservoirs including naturally fractured ones (Sarma, 2006). The production duration life of an oil reservoir generally lasts a few decades and maximum three different production phases may be defined during that period. It should be noted that theoretically between 20%-70% of the oil can be recovered by secondary production techniques (Van den Hof et al. 2009). Oil production process from reservoirs can be considered as a large scale dynamical system. Unfortunately, although an oil reservoir is a complex system, most of the time the operators manage the fields based on their experiences and consequently the complex physics of the system would be ignored. Obviously, besides operating the production based on the operator experience, there are chances to optimize the production by model-based control and optimization techniques, on the basis of predictions of mathematical reservoir models and simulators (Van den Hof et al. 2009). Propitiously, with recent developed technology for drilling and also reservoir instrumentation, there is a capability for utilizing model-based control and optimization approaches to develop applicable reservoir management strategies in order to optimize the economic efficiency of the production process (Van den Hof et al., 2009; Jansen et al. 2008). For choosing the most beneficial scenario of water injection in secondary recovery approach, a comprehensive numerical tool or mathematical model for predicting the performance is necessary (Fanchi, 2001). Reservoir simulation prepares an estimation tool of reservoir performance. There are several techniques of simulation and modeling from simple to complex. The selection of each of the available approaches depends on the accessible data, the level of required accuracy and the calculated time-cost (Ertekin, 2001). Methods such as Decline Curve Analysis (DCA)(Baker, 2003) and the Capacitance Resistance Model (CRM) (Delshad et al., 2009; Moghadam et al., 2011), require less data to model the reservoirs performance. Grid-based simulation (Aziz & Settari 1979; Fanchi, 2001) and streamline simulation (Sayyafzadeh et al., 2010) are the most precise approaches however their simulation time is too much and they need large quantity of data. All the required information for reservoir modeling such as porosity and permeability distributions are not always in hand, in addition they contain high degree of 2
uncertainty; so costly and time-consuming history matching task is required for model validation (Oliver & Chen, 2010). It is clear that an accurate reservoir model is not always necessary for all applications and a fast approach with acceptable accuracy is sufficient. Furthermore, all the required data for a grid-based simulation are not always accessible. As a result, finding a quick modeling technique which needs less data, but has reasonable accuracy is an active research area in the field of reservoir modeling and researchers look for faster methods which need less data, but has reasonable accuracy (Sayyafzadeh, 2011). In Van Essen et al. (2012), Lee et al. (2010), Gildin et al. (2008) and Tafti et al. (2013), suitable proxy/surrogate models for control and optimization applications have been introduced. However, the mentioned references have relied based upon a fixed linear (structure) presentation for the whole range of reservoir operation and consequently this approach will limit their applications within the range of the induced injection rates around the operating point. On the other hand, Using Neural Networks for reservoir modeling has been discussed by Mohaghegh & Abdulla (2014), Aifa (2014), and Bruyelle and Guérillot (2014). However, the obtained models based on neural network are usually in the category of nonlinear models. Therefore, controller or optimizer designing based on neural network methods is not as simple as handling the linear models. In this paper, a systematic approach for adaptive and online modeling of oil reservoirs applicable for the secondary production phase is presented via linear techniques. In spite of the method’s linearity, the results demonstrates that the complicated waterflooding process can be modeled accurately. Based on the introduced modeling mechanism, the field experts will be able to plan for reservoir management and apply advanced control and optimization strategies during the operation in order to achieve to the main goal which is improving the economic performance over the production life of an oil reservoir. 2. Problem statement: Availability of a valid and reliable model is required for any controlling or optimization target in the reservoir management approach. Consequently, obtaining a proper and precise model is a prerequisite for any controlling application.
3
Inherently, a reservoir system model consists of partial differential equations which are based on mass and momentum conservation laws. For solving the model equations numerically, it is needed to discretize the equations in space and time and as a result a large scale system including thousands of grids and even millions of states would be constructed. 2.1. Reservoir Model Reservoir model, based on conservation of mass and momentum is generally used to describe the fluid flow in the reservoir and consequently the dynamics of waterflooding process can be obtained based on the involved equations (Aziz & Settari 1979). By considering oil and water in the reservoir and neglecting the presence of gas, the model equations could be simplified (Jansen et al. 2008). For solving the reservoir equations available in (Aziz & Settari 1979)numerically, it is necessary to discretize them in the space which leads to a system built up of finite number of grids, called as grid blocks. After discretization in time, the following state space representation will be obtained:
V ( xk ) xk 1 T ( xk ) xk qk ,
x0 x0 , (1)
where k is the time index and x is the state vector including the oil pressure po and water saturation S w in all grids. Vector x0 contains the initial conditions, which are assumed to be known. For modeling the effect of the wells on the dynamic of the reservoir, a source vector qk is included in (1). The source terms are usually represented well model, which connects the source term to the pressure difference between the well and grid block pressure:
qkj w j ( pbhj , k pkj ) (2) j
in which pbh , k is the well's bottom hole pressure, j the index of the grid block containing the well j
and pk is the grid block pressure in which the well is located. The term w is a well constant that represents the well's geometric factors and also the rock and fluid properties of the reservoir near the well. The geological properties of each grid block are assumed to be constant. The heterogeneity of the reservoir can be specified by assigning different property values to each of the blocks. Usually a very large number of grid-blocks is needed to describe the dynamic of a 4
real oil reservoir. In the commercial numerical simulators the above equations are solved for all grids in each time step which result massive calculation load. 2.2. MIMO Reservoir Model Structure Generally, the typical oil reservoir models have MIMO structure and in practice it may include 10 to 1000 wells. Furthermore, these models are nonlinear due to the relative permeability effects and also possessing high order because of spatial discretization. The conventional equation-based methodology for reservoir modeling which has been mentioned in section 2.1 is often used in commercial and numerical simulators, leading to a thorough model which can simulate most behaviors of the reservoir. However, in many cases a high percentage of the provided information by the mentioned simulator-based approach is not essential during the actual operation. A hydrocarbon field is practically operated just based on the monitored variables such as total injection and production rates or water and oil production rates of all of the wells. As a result, there is no need for modeling all the possible phenomena in an oil reservoir for different field usages and each application requires only limited and related variables. From the control theory point of view, proper input and output ports should be defined for the system. In the other words, the inputs and outputs of a system are essential to be accessible. The inputs of the reservoir which influence the outputs are called well controls which can be manipulated by flow rate control or BHP control. On the other hand, the main target in reservoir management is controlling the oil production rate of each well and also the produced water. So, the mentioned variables can be considered as the system outputs. Therefore, the reservoir model can be considered as depicted in Figure. 1.
5
Fig. 1. Schematic of oil reservoir as a controlled system with input/output determination 2.3. Application of System Identification in Reservoir Modeling “System Identification” (SI) is a well-known approach for practical modeling of an unknown system for controlling and optimization applications. SI modeling approaches normally consider the plant as a black-box with specified inputs and outputs. In this approach, the ultimate goal is to find a suitable mapping between inputs and outputs through ignoring the internal details of the model for providing the ability of appropriate system outputs estimation in the presence of different and challenging conditions. After finishing model training stage, based on the validated model outputs actions such as prediction, controlling or optimization tasks can be perform. In discrete-time system modeling, one alternative is defining the desired output as a function of past outputs and also all the inputs at previous sampling times. The memory of the model and also the time-delay of the system should be obtained as one of the important stage in the modeling approach. Although in this paper all the results are obtained in discrete-time domain, it should be noted that switching from discrete-time to continuous-time modeling is a straight-
6
forward task which could be easily done according to the required applications (Sung et al., 2009). 2.4. Problem definition Since a petroleum reservoir can be considered as a nonlinear system with different operating points based on the value of the inputs, one of the most applicable solutions is to linearly model its dynamic in the vicinity of each operating point. This approach profits from linear system theories which have strong backbone to be implemented in practice. According to “Parsimony Principle” (Verakas & Vlahavas 2008), for identification of the plant’s behavior, it is recommended to start from the simplest structure and if the model outputs are far from the reality, then improving and upgrade the model structure will be done. As explained previously, in case of considering the reservoir as a MIMO system, the outputs may be defined as oil and water production rate of each production well which can be controlled by two independent inputs named rate of water injection of each injection well and also BHP of the production wells. In practice, it is not an invalid assumption to fix the BHP of each production well as a fixed control input in the desired time interval and consequently the main manipulated variable would become the flow rate of injection wells. At this time the problem can be defined as follows: “Extracting the value of the outputs at fixed BHPs and in presence of injection rate variation”. It is necessary to state that fixing BHP does not have negative effect on the generality of the problem. The proposed method in this paper has the capability to generate the reliable outputs at each BHP and consequently by BHP variation during the operation the updated model will be realized. In addition, by applying MIMO identification procedure which is explained in next sections, general modeling of waterlooding process is feasible.
3. Identification methodology 7
Identification is system modeling based on different experiments. If the tests are not selected correctly then no identification algorithm even the sophisticated ones, can generate a valid model of the process. On the other hand, if the input/output data are obtained from a well-designed experiment, least-squares approach which is the simplest method in this domain can often present a good model. Generally, the purpose of experiments is collecting relevant information about the process dynamics. In the next step, the obtained information should be transformed to valid mathematical models of the process by appropriate identification algorithms. Usually, the collected data from the experiments are not suitable for immediate use in identification algorithms. Possible deficiencies such as high frequencies disturbances; occasional spikes and outliers; drift and offset, low-frequency disturbance; and different order of magnitude may exist. To avoid encountering problems at the estimation steps data pre-processing actions like: Peak Shaving, Trend Correction, Scaling and Offset Correction, Delay Correction, Filtering and Sampling Rate Reduction can be performed (Zhu & Backx, 1993). 3.1 ARX structure modeling approach based on LS After finishing pre-processing task, the first practical step in identifying a physical system is performing offline Least Square (LS) algorithm on the batch input-output data by assuming the system as linear time invariant (LTI). In this method model parameters will be adjusted such that a cost function which is sum of square errors (SSE) be minimized. This approach has been frequently used in industrial processes and successful results have been reported in practical applications (Kon et al. 2013, Chetouani 2008) . A discrete dynamic system can be described by the following equations:
yˆ (k ) T (k 1) (3)
y(k ) yˆ (k ) (k ) (4)
T (k 1) [ y(k 1),..., y(k na ), u(k 1),..., u(k nb )] (5)
8
T (1) (k ) (6) T (k )
T [a1 ,..., ana , b1,..., bnb ] (7) in which u(k) and y(k) are the plant’s input and output at time k (sampling-time), respectively. In reservoir modeling literature, inputs are considered as injection flowrates and bottom hole ^
pressure of wells and outputs are assumed as oil and water production of wells. y (k ) is the model’s output and ε(k) is the model error (difference between the output of the reservoir and the model) and its nature is a zero-mean white Gaussian noise term, which can represent the measurement noise or even the model uncertainties. Regular LS algorithm suppose that the error
T (k 1) between plant output and the model output can be considered as white noise.
is the
vector which contains the available information and it is called the Regressor and Ф (k) is the information matrix which contains the previous information of injection/production values in oil reservoir. is the identification parameter vector which will be determined in the identification process. The parameters in can either be constant or variable. It should be noted that in reservoir modeling parameter The main goal is properly estimation of the values for model parameters: s, hence, all available measurements including inputs and outputs up to the current instant j = 1, 2, …, k, should be utilized in order to minimize the following cost function which is summation of the existing error between model and real output at each sample-time: k
J (k ) (i) (8) 2
j 1
which can be written as:
J ( , k ) j 1[ y( j ) yˆ ( j )]2 (9) k
9
and as a result: J ( , k ) j 1[ y( j ) T ( j 1) ]2 (10) k
Differentiating J(θ,K) with respect to θ and setting the obtained equation equal to zero, the minimum cost function value would be available when the estimate parameter vector is as (Lakshminarayanan, 2001): follows ^
[ j 1 ( j 1) T ( j 1)]1[ j 1 ( j 1) y( j )] (11) k
k
^
where is considered as the least squares estimate of (Ljung & Soderstrom, 1983; Jiang & Zhang, 2004) . This approach results to Auto-Regressive with eXogenous input (ARX) model which has the following format in discrete linear system modeling literature (Zhu & Backx, 1993):
y (k )
B(q 1 ) 1 u (k ) n(k ) (12) 1 A(q ) A(q 1 )
which can be written as:
y (k ) a1. y (k 1) b1.u (k 1)
ana . y (k na ) bnb .u (k nb ) (13)
n( k ) 3.2. ARMAX structure modeling approach based on ELS The above modeling methodology implies that the output of the supposed linear system is a function of previous outputs, inputs and also an unknown term which has the characteristics of the white noise. Hence, it is supposed that the nature of the model error or uncertainty is white noise. However in many cases the dynamic of the error is much complicated and it should be modeled differently which one option is assuming that dynamic as colored noise. In this case, 10
Extended Least Square (ELS) modeling approach can be the next step alternative for being tested and it results to Auto Regressive Moving Average with eXogenous input (ARMAX) model which has the following structure:
B(q 1 ) C (q 1 ) y (k ) u (k ) n(k ) (14) A(q 1 ) A(q 1 ) By considering C(q-1), a new degree of freedom would be added to the model in comparison with ARX model . The estimated error in this approach is defined similar to ARX modeling as:
(k ) y(k ) T (k 1)ˆ(k 1) (15) However, the parameter vector includes c:
[a1
ana b1
bnb c1
cnc ] (16)
and Regressor vector is as follows:
T (t 1) [ y (k 1) u (k 1)
y (k na ) u (k nb )
(k 1) (k nc )]
(17)
and consequently the model parameters can recursively obtained based on (18):
ˆ(k ) ˆ(k 1) P(k ) (k 1) (k ) (18) in which
P1 (k ) P1 (k 1) (k 1) T (k 1) (19) where P(k) is called the covariance matrix. However, since this type of modeling has inherently linear structure, there would be cases in which the obtained model by this approach is not acceptable for critical applications such as controlling the desired outputs or optimization of an objective function.
11
As most of the systems/plants in nature, are nonlinear and also time varying, in some occasions, pure ARX and even ARMAX system identification modeling structures are unable to present the main dynamics of the required system. In this case, other approaches such as Recursive Least Square (RLS) method which is a good choice for practical applications are needed. 3.3 RLS modeling approach
The main idea in the RLS algorithm is re-computing of the parameter in each sample time (i=1,2,…k) by considering a correction term based on new measured data at kth sample. This strategy decreases the computational complexity of RLS and make it suitable for on-line regime and adaptive applications (Jiang & Zhang, 2004). A typical presentation of the RLS is as follows: ^
^
(k ) (k 1) P(k ) (k 1) (k ) (20) in which, ^
(k ) y(k ) T (k 1) (k 1) (21) and also,
P(k 1) (k 1) T (k 1) P(k 1) P(k ) P(k 1) (22) 1 T (k 1) P(k 1) (k 1)
where P(k) is covariance matrix and ε(k) is one-step ahead prediction error. It is obvious that the larger prediction error would cause an increase in the correction term (Niedzwiecki, 2000). In time varying systems, for estimating the new values of the parameter, in addition to recursively calculating the system parameters, a mechanism is needed to ignore the past
12
information that represent the former dynamics which may not be valid at the current time. As a result, techniques such as forgetting factor or finite length sliding window may be used in association with the RLS. By adding term
and L to the cost function, the influences of previous data in parameter
updating would be decreased and the estimated parameters would be find in such a way that (24) be minimized.
J ( , k ) j k L1 k j [ y( j 1) T ( j 1) ]2 (23) k
In (24) L shows the length of the sliding window. In addition, λ which its value is in the interval of [0, 1] is called the forgetting factor for exponentially discarding the obsolete data. Now, the LS problem for parameter estimation of θ is minimizing the modified cost function with the most recent data. In this case, updated parameter vector with utilizing sliding exponentially weighted window least squares algorithm is as follows:
1
k k k j T L (k ) ( j 1) ( j 1) k j T ( j 1) z ( j ) (24) j k L1 j k L1 ^
The performance of the above algorithm will depend on the length of the sliding window. For time-invariant systems the better estimation accuracy will be achieved by increasing the window length. However, for time-variant systems, the window length should be adjusted so that the obsolete information from the past measurements can be discarded properly and heavier weight on the latest measurement be assigned for tracking the system parameter changes. Since in the presented study the system parameters are time-varying due to the nature of the reservoir and also operating point changes, the appropriate selection of window length and forgetting factor value is an important stage (Paleologu, 2008).
13
If the responsibility of discarding old information be assigned to forgetting factor λ then L=k in (26) and RLS algorithm for parameters updating can be written that: ^
^
^
(k ) (k 1) K (k )( y(k ) T (k ) (k 1)) (25) in which:
K (k ) P(k ) (k ) P(k 1) (k )( I T (k ) P(k 1) (k )) 1 (26)
and
P(k ) ( I K (k ) T (k )) P(k 1) / (27) Due to the nonlinear nature of relative permeability and some other parameters of an oil reservoir, an oil reservoir is strongly nonlinear system. In addition, the characteristics of a reservoir would change during oil production and so the reservoir is a nonlinear time-variant plant by the aspect of system theory. Fortunately, RLS technique has the ability to consider both nonlinearity and also time varying nature of the plant simultaneously. In addition, since changing the inputs of a nonlinear system, may cause operating point alternation, the RLS method is also able to locally linearize the plant near the operating point and produce a linear model at each sample time. Furthermore, the parameters of the identified models would be updated in each sample time and consequently any variation in the reservoir specification or nonlinear system’s operating point can be reflected in the new updated parameters. 3.4 MIMO System Identification:
Since, the above mentioned methodologies have been developed for SISO (Single-Input-SingleOutput) plants; in this section the extension to MIMO case would be introduced. This is due to 14
most physical systems in nature including oil reservoir can be considered as MIMO system and several inputs and outputs exist.
A suitable presentation of a MIMO process is a set of difference equations for identification purpose as follows:
y(k ) A1. y(k 1)
An . y(k n) B1.u(k 1)
Bn .u(k n) (28)
in which y is the p-dimensional output vector, u is the m-dimensional input vector and A1 (p×p), ... An (p×p), B1 (p×m) ... Bn (p×m) are constant matrices. Using the unit delay operator q-1,
A(q) y(k ) B(q)u(k ) (29) in which A(q) and B(q) are polynomial matrices:
A(q) I A1q 1
B(q) B1q 1
Anq n (30)
Bn q n (31)
generally, in MIMO Systems ,
y G(q)u (32) Now, the transfer matrix of the system can be written as:
G(q) A(q)1 B(q) (33)
Expressing the transfer matrix by the product of two polynomial matrices is called matrix fraction description (MFD). The parameters of matrices A(q) and B(q) can be estimated from the input/output data. In contrast to the SISO case, for a given process, equation (29) is in general not unique even by fixing the leading coefficient matrix of A(q) at I. Hence, some special forms of representations such as canonical forms are required for obtaining unique models of a given process. Model 15
uniqueness is a necessary condition for an identification algorithm to find a proper solution. Using diagonal canonical form MFD is an acceptable option for MIMO processes. The diagonal form MFD of a given G(q) is
A11 (q) 0 0 A22 (q) A(q) 0 0
B11 (q) , B(q) 0 B p1 (q) App (q) 0
B1m (q) (34) B pm (q)
where A 11 (q), …, App(q) are all monic polynomials and the degrees of Bi1(q), …, Bim(q) are equal to or less than that of Aii(q). with this arrangement, m-input, p-output process is decoupled into p m-input single output sub-processes:
A11 (q) y1 (t ) B11 (q)u1 (t )
B1m (q)um (t ) (35)
App (q) y p (t ) B p1 (q)u1 (t )
B pm (q)um (t )
For physical processes, the transfer matrix can be presented in the diagonal form MFD and consequently most of SISO processes identification algorithms are applicable for MIMO cases in this framework. On the other hand, the diagonal form MFD of a given process is not always minimal. So, the order of state space realization can be higher than the McMillan degree of the process. This problem can easily be solved by model reduction techniques (Benner & Faßbender, 2014).
A11o (q) y1 (t ) B11o (q)u1 (t )
B1m o (q)um (t ) A11o (q ) 1 (t ) (36)
App o (q) y p (t ) B p1o (q)u1 (t )
B pm o (q )um (t ) App o (q ) p (t )
Without loss of generality, it has been supposed that the degrees of all polynomials in a submodel are equal. The orders of sub-models are n1, …nm. The Input/Output data sequence is as follows in which n=max{n1, …, nm}. 16
Z N : y1 (1)
y p (1) u1 (1)
um (1)
y1 ( N n)
y p ( N n) u1 ( N n)
u m ( N n)
(37) Assuming that N samples of measurements of outputs and inputs are made at time 1, 2,…, N. Now, the loss function which should be minimized in order to estimate the required parameters is:
V1
1 N n ( A1 (q) y1 (t ) [ B11 (q)u1(t ) N t n1
B1m (q)um (t )]
(38)
and consequently,
ˆ1 [
1 T 1 1 1 ]1 1T y1 (39) N N
where
y1 (n 1) y (n 2) , y1 1 y1 (n N )
y1 (n1 ) y (n 1) 1 1 1 y1 (n1 N )
a1,1 a1,n 1 b11,1 (40) 1 b11,n 1 b 1m ,1 b1m ,n1 y1 (1) | u1 (n1 )
|
| um (n1 )
| u1 (n1 1)
|
| um (n1 1)
|
|
|
y1 ( N 1)| u1 (n1 N )
u1 (1)
u1 ( N 1) |
17
| um (n1 N )
(41) um ( N 1) um (1)
The same procedure can be performed for other sub- models. In addition, it has been supposed that the degrees of all polynomials in a sub-model are equal. By utilizing the above presented approach, it will be possible to extract a MIMO model for any type of oil reservoir with several inputs and outputs. However, although there exists different techniques for directly MIMO modeling of real systems, simplifying the problem to several MISO identification can be a proper alternative for modeling the system (Ding & Chen, 2005a; Ding & Chen, 2005b; El-Anes, 2013).
4. Methodology of Modeling Oil Reservoir: At this section the developed methodology for modeling oil reservoirs is introduced. This technique can be break-down to three sub-stages as follows: 4.1. Finalizing the General Structure In this phase, the source of data would be defined. The data source may either be the real reservoir data which have been gathered from the field or data which have been generated by a valid commercial or academic reservoir. In addition, all manipulated variables which can act as the system inputs and also the desired outputs should be fixed at this step. The main goal of the algorithm is assigning reliable model/models which can demonstrate the behavior of the system appropriately. 4.2. Data Generation As explained in the previous sections, recording proper input-output data is one of the most vital phases in system modeling and identification. The richer collected data would result to the better obtained model. So, it is necessary to stimulate the plant (here; the reservoir) such that all the hidden dynamics become apparent. For this purpose, Persistently Exciting (PE) signals should be applied as system inputs. A good choice of rich inputs is white noise or even Pseudo Random Binary Signal (PRBS). However, since in practice, applying those inputs to the reservoir is 18
somehow impossible and the setting of the control valves cannot be modified with high frequencies, an appropriate signal which has the most similarity to PRBS and also meets the operational constraints should be applied to the system and all inputs/outputs be recorded simultaneously. In the next step before starting the main identification process, according to the above explanations, data pre-processing is necessary to be performed in order to enhance the numerical robustness of identification algorithms and also accelerate the speed of convergence. 4.3. Model Identification In the final step, the proper identification algorithm should be utilized for achieving to the suitable model. According to the “Parsimony Principle” in system modeling and identification, the candidate structures starts from the simplest ones and if there exists a large discrepancy between the model outputs and real plant outputs, a more complicate structure will be examined. So, in the proposed methodology, the start point is evaluating the obtained results from an offline batch ARX structure based on least square technique and if the results are unsatisfying then ARMAX structure and adaptive RLS will be tested respectively. It is important to notice that it is helpful to initialize the parameters of more complicated structures based on the obtained results of the former tested strategy and this approach will significantly increase the convergence speed of the algorithm to find the proper model. The flowchart of this methodology has been illustrated in Figure 2.
19
Fig. 2. Flowchart of the Methodology.
5. Selected Reservoir Case Study The selected model 10th SPE-Model#2 has a simple geometry, with no faults and top structure however, it is one of the well-known case study in the field of reservoir researches. At the fine geological model scale, the model is described on a regular Cartesian grid. The model dimensions are 1,200×2,200×170 ft. The top 70 ft (35 layers) represents the Tarbert formation, and the bottom 100 ft (50 layers) represents Upper Ness. The fine-scale cell size is 20×10×2 ft which results that the model has 60×220×85 grid-cells.. The model consists of part of a Brent sequence. The model was originally generated for use in the PUNQ project. Although in the original model a uniform kV/kH in the whole domain is considered, the used model has been changed such that kV/kH is equal to 0.3 in the channels and 10-3 in the background. The porosity and permeability maps and other model parameters are available in the SPE homepage (Christie & Blunt, 2001; Islam & Sepehrnoori, 2013). 20
In this synthetic model, all wells are vertical and completed throughout the whole formation. In the original model the central injector has an injection rate of 5000 Barrles/ Day. There are four producers in the four corners of the model; In the initial example each well produces at 4,000 psi bottom-hole pressure. The X-Y location of each well is presented in Table 1 in Cartesian coordinate: Table 1. Well Locations in Cartesian coordinate. well name
X location
Y location
(ft)
(ft)
Injection Well I1
600
1100
Production Well P1
0
0
Production Well P2
1200
0
Production Well P3
1200
2200
Production Well P4
0
2200
6. Implementation of the presented algorithm and results: For implementing of the developed algorithm, Matlab reservoir simulation toolbox (MRST) has been used (Lie, 2014). This simulator provides the facility to control each well by adjusting flowrate or BHP. Since, the ultimate goal of system modeling is controlling the desired outputs or even optimization the whole reservoir performance, the manipulated variables would be selected as the injection rate of injection well(s). In addition, it is supposed that the production wells are running under the fixed BHP due to the operational recommendations. Furthermore, the output of the system which are going to be identified can be selected as oil production rate and water production rate of producing wells. In order to maintain the pressure of the reservoir constant, it is supposed that the total injection rates of injection wells are equal to the total production rates of producing wells. 21
m
p
i 1
j 1
qi _ inj q j _ prod (42) Where qi_inj is the flowrate of each injection well (here, flow rate of injected water) and qj_prod is the flowrate of each producing well (here, summation of produced oil and water). When the inputs are set, the simulator can simultaneously calculate and generate total production rate, oil production rate, water production rate and water cut of each well. The water cut is the ratio of water produced compared to the volume of total liquids produced from an oil well.
It is obvious that availability of only two of the mentioned variables, may lead to calculate the others. In practice, each of them can be measured uniquely by applying proper instruments in the line. Now, step by step the proposed algorithm is applied to the well-known 10th SPE-Model#2 for identification of the waterflooding process in an oil reservoir supposing that the BHPs of the production wells are constant which conclude to a Single Input, Multi Output (SIMO) system (one injection well and four production wells). In simulation, 2000 days of reservoir operation with a proper injection rate properties have be analyzed. For stimulating the dynamic of the reservoirs and also for identification the hidden modes, a PE input is required. However, in reservoir management, it is not suitable to change the value of each manipulated variable continuously and in specified times which may be dictated by the experts those values would be changed. So, a proper candidate input sequence which is also applicable in practice and has the most analogy with a pure PE signals can be seen in table 2. Table 2. Sequence of Injection rate Injection # 1 2 3 4 5 6 7 8
Injection value (bbls/day) 5000 4500 6000 4000 5000 6000 3000 5500 22
Starting day
Ending day
1 201 301 501 801 901 1001 1301
200 300 500 800 900 1000 1300 1500
3500 4500 5000
9 10 11
1501 1601 1801
1600 1800 2000
In this case, it has been supposed that the injection rate may only change in each 100 days randomly (10 samples) for having the most similarity to the practical experiments. In addition, it is assumed that the values of injection rates are between 3000bbl/day and 6000bbl/day. Since an oil reservoir is a nonlinear system, these values would help to identify its different operating points and also facilitate to identify the reservoir’s dynamic close to each operating point. The output of the system has been monitored and recorded every 10-days and base on the introduced algorithm, waterflooding modeling and identification procedure has been started. The results show that ARX and ARMAX modeling approaches are unable to precisely model the dynamic of the reservoir, since they try to fit batch models for whole production life. The best obtained result for ARX model with structure in (43), is in (44):
A(q) y(k ) B(q)u(k ) e(k ) (43) and
A(q) 1 q 1 0.0641q 2 (44)
B(q) 0.0298q 2 which tries to model the relationship between total production of well #3 and injected water. Figure 3. Illustrates a sample result for evaluating the performance of ARX structure for estimating the total production of Well #3. It is clear that although the output of ARX model can somehow track the trend of the real value of the production rate, it is far from the desired simulator’s output.
23
Fig. 3. Performance of the ARX model output versus simulator’s output. In the next step, ARMAX structure should be evaluated. The best result of this model structure for the same input/output variable is:
A(q) 1 1.657q 1 0.6816q 2 B(q) 0.01143q 2
(45)
C (q) 1 0.657q 1 0.0001q 2
With the following structure:
A(q) y(k ) B(q)u(k ) C (q)e(k ) (46) Although Figure 4. presents that the performance of this model is slightly better than ARX model, however it cannot be considered as a reliable waterflooding model during the production period.
24
Fig. 4. Performance of the ARMAX model output versus simulator’s output.
The quality of obtained results is predictable since modeling a strongly nonlinear plant with linear approaches in the whole domain is theoretically impossible. However, RLS approach can be an appropriate candidate for extracting and identification of the behavior of the system and predict the required output. By applying this technique, at each sampling time (in this case, every 10 days), a new localized linear model will be generated based on the available data up to that time and consequently any parameter or structural change in the system can be detected and be reflected in the new model. As the following results presents, by adjusting the forgetting factor of RLS to 0.95, this technique can identify all of the required outputs for all production wells precisely and the results obtained of the RLS strategy are very similar to the outputs of the simulator. It is important to notice that in this condition, for each sample time (each 10 days) a new ARX model is valid for indicating the behavior of the plant which its parameters have been also illustrated in the figures. Now this set of models is valid for representing the whole production period behaviors and it has also the ability for detecting any variation in the reservoir’s characteristics. It can be seen in most of the results that any change in the input, which may cause the operating point change in the plant can cause the model’s parameter variation. Furthermore, this approach of waterflooding 25
modeling has the capability to regenerate the new models whenever the BHP of the production wells also alter and its application is not limited to only constant BHP cases. Table 3. shows the value of sum of square error (SSE) between simulator output and each of the 3 modeling strategies’ output related to the total production rate of well#3.
Table 3. Well Locations Modeling approach
SSE
ARX
7.9050e7
ARMAX
6.6519e7
RLS
3.6437e4
Figures 5, 6, 7 & 8 illustrate the simulator outputs and generated models’ outputs in which the outputs are defined as total production rate of the wells and system input is the water injection rate. For another example, Figures 9, 10, compare the estimated and real oil production rates for two selected wells. In figures 11, 12, the estimated produced water of well #2 & well #4 can be observed and finally figure 13 showes the water cuts of well #1 for instance. In addition, the estimated model parameters for all desired outputs have been depicted in each picture in which ai ‘s are the denominator and bi’s are the nominator coefficients of the dedicated model according to (12) and as a result in each time step a unique model is available.
26
Fig. 5. Injection rate, Total Production Rate & Model coefficients related to well #1.
Fig. 6. Injection rate, Total Production Rate & Model Coefficients related to well #2.
27
Fig. 7. Injection rate, Total Production Rate & Model Coefficients related to well #3.
Fig. 8. Injection rate, Total Production Rate & Model Coefficients related to well #4.
28
Fig. 9. Injection rate, Oil Production Rate & Model Coefficients related to well #1.
Fig. 10. Injection rate, Oil Production Rate & Model Coefficients related to well #3.
29
Fig. 11. Injection rate, Water Production Rate & Model Coefficients related to well #2.
Fig. 12. Injection rate, Water Production Rate & Model Coefficients related to well #4.
30
Fig. 13. Injection rate, Water cut & Model Coefficients related to well #1.
7.
Conclusion:
In this paper an adaptive identification algorithm has been developed for modeling the waterflooding process in oil reservoirs. The presented approach has the ability to reflect the timevarying nature of the reservoir with its inherent nonlinearity. The developed approach, following the parsimony principle conducts a step by step model scheme to adaptively identify a simple candidate model structure through a preliminary off-line ARX model. Then, the approach goes through more complex ARMAX model structure in an offline manner if the obtained model accuracy is not enough for application target. An adaptive online model will subsequently adopted to capture the reservoir dynamics in case the offline batch model is not adequate to accomplish the required model accuracy. The observed results on SPE10 benchmark case study have shown that the RLS approach has the ability to estimate each of the required reservoir outputs with acceptable accuracy. Forgetting factor parameter adjustment in RLS method adds the ability of discarding old data for obtaining the most reliable model in each time. Since the presented technique is generally developed for MIMO plants, it can easily include any type of oil reservoir with several injection/production wells and also adjusting flowrate or BHP of each well simultaneously. The ultimate objective of this kind of modeling is to generate proper models for 31
control and optimization goals. As a result, availability of such adaptive models enables to apply useful adaptive control approaches for desired output(s) and also perform optimization tasks on the reservoir during the operation
References: 1. Aifa, T., Neural network applications to reservoirs: Physics-based models and data models, Journal of Petroleum Science and Engineering, 123, pp.1-6, 2014. 2. Aziz, K., Settari, A.N., 1979. Petroleum reservoir simulation. London, Applied Science Publishers. 3. Baker, R.O., Anderson, T., Sandhu, K., Using Decline Curves to Forecast Waterflooded Reservoirs: Fundamentals and Field Cases. Canadian International Petroleum Conference. Petroleum Society of Canada, Calgary, Alberta. 2003. 4. Benner, P., Faßbender, H., Model Order Reduction: Techniques and Tools, Encyclopedia of Systems and Control, Springer-Verlag, 2014. 5. Bruyelle, J., Guérillot, D., Neural networks and their derivatives for history matching and reservoir optimization problems, Computational Geosciences, Volume 18, Issue 3, pp 549-561, 2014. 6. Chetouani, Y., Using ARX approach for modelling and prediction of the dynamics OFA reactor-exchanger, Symposium Series NO. 154, IChemE 2008. 7. Christie, M. A., and M. J. Blunt. Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques. SPE 66599, SPE Reservoir Simulation Symposium, Houston, Feb 11‐14, 2001. 8. Delshad, M., Bastami, A., Pourafshary, P.,The Use of Capacitance-Resistive Model for Estimation of Fracture Distribution in the Hydrocarbon Reservoir. SPE Saudi Arabia Section Technical Symposium. Society of Petroleum Engineers, AlKhobar, Saudi Arabia. 2009. 9. Ding, F., Chen, T., Hierarchical gradient-based identification of multivariable discretetime systems, Automatica 41: 315 – 325, 2005. 10. Ding, F., Chen, T., Hierarchical Least Squares Identification Methods for Multivariable Systems, IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 50, NO. 3. 2005. 11. El-Anes, A., Bouzrara, K., Garna, T., Messaoud, H., Expansion of MIMO ARX model on Laguerre orthonormal bases. International Conference on Electrical Engineering and Software Applications (ICEESA), 2013. 12. Ertekin, T., Abou-Kassem, J.H., King, G.R., Basic applied reservoir simulation. Society of Petroleum Engineers, Richardson, Tex. 2001. 13. Fanchi, J.R., Principles of Applied Reservoir Simulation, Second Edition. Gulf Professional Publishing, Burlington. 2001. 14. Gildin, E., Wheeler, M.F., Control of Subsurface Flow Using Model Predictive Control Techniques, in Proceedings of the International Conference on Engineering Optimization, Brazil, June, 2008. 32
15. Islam, A. W., Sepehrnoori, K. A Review on SPE’s Comparative Solution Projects (CSPs), Journal of Petroleum Science Research (JPSR) Volume 2 Issue 4, October 2013. 16. Jansen, J.D, Bosgra, O.H, and Van den Hof, P.M.J, Model- based control of multiphase flow in subsurface oil reservoirs. Journal of Process Control, 18(9):846-855, 2008. 17. Jiang, J., Zhang, Y., A revisit to block and recursive least squares for parameter estimation. Computers & Electrical Engineering; 30(5): 403-416. 2004. 18. Kon, J., Yamashit, Y., Tanaka, T., Tashiro, A. and Daiguji, M. Practical Application of Model Identification based on ARX Models with Transfer Functions, Proceedings of the 6th International Conference on Process Systems Engineering (PSE ASIA), 2013. 19. Lakshminarayanan, S., Emoto, G., Ebara, S., Tomida, K., and Shah, S.L., Closed loop identification and control loop reconfiguration: An industrial case study. Journal of Process Control 11: 587–599, 2001. 20. Lee, K.H., Ortega, A., Jafroodi, N., & Ershaghi, I., A Multivariate Autoregressive Model for Characterizing Producer-producer Relationships in Waterfloods from Injection/Production Rate Fluctuations. Society of Petroleum Engineers, 2010. 21. Lie, K. A., An Introduction to Reservoir Simulation Using MATLAB: User guide for the Matlab Reservoir Simulation Toolbox (MRST). SINTEF ICT, May 2014. 22. Ljung, L., and Soderstrom, T., Theory and practice of recursive identification. MIT Press, Cambridge. MA. 1983. 23. Moghadam, J.N., Salahshoor, K., Kamp, A. M., Evaluation of Waterflooding Performance in Heavy Oil Reservoirs Applying Capacitance-Resistive Model, Petroleum Science and Technology, 29:17, 1811-1824, 2011. 24. Mohaghegh, Sh., Abdulla, F., Production Management Decision Analysis Using AIBased Proxy Modeling of Reservoir Simulations – A Look-Back Case Study, Society of Petroleum Engineers, 2014. 25. Niedzwiecki, M., Identification of time-varying processes. Wiley 2000, Chichester. 26. Oliver, D., Chen, Y., Recent progress on reservoir history matching: a review. Computational Geosciences 1–37. 2010. 27. Paleologu,C., Benesty, J., Ciochina,S., A Robust Variable Forgetting Factor Recursive Least-Squares Algorithm for System Identification, IEEE SIGNAL PROCESSING LETTERS, VOL. 15: 597-600, 2008. 28. Rossi, D., Carney, M., Kontchou, J.N., Lancaster, D., McIntyre, S., Reservoir and Production Optimization, white paper from www.slb.com, 2002. 29. Sarma, P., Efficient Closed-Loop Optimal Control of Petroleum Reservoirs under Uncertainty, PhD Thesis, Stanford University, Stanford, CA, 2006. 30. Sayyafzadeh, M., Pourafshari, P., Rashidi, F. Increasing ultimate oil recovery by infill drilling and converting weak production wells to injection wells using streamline simulation. International Oil and Gas Conference and Exhibition in China. Society of Petroleum Engineers, Beijing, China., 2010. 31. Sayyafzadeh, M., Pourafshary, P., Haghighi, M., Rashidi, F., Application of transfer functions to model water injection in hydrocarbon reservoir, Journal of Petroleum Science and Engineering 78 ,139–148, 2011. 32. Sung, S.w., Lee, J., Lee, I. B., In-Beum Lee, Process Identification and PID Control, John Wiley & Sons, 2009.
33
33. Tafti, T., Ershaghi, I., Rezapour, A., & Ortega, A. Injection Scheduling Design for Reduced Order Waterflood Modeling. Society of Petroleum Engineers, 2013. 34. Van den Hof, Paul M.J.; Jansen, Jan-Dirk; Essen, Gijs van; Bosgra, Okko H. Modelbased control and optimization of large scale physical systems - Challenges in reservoir engineering, IEEE Chinese Control and Decision Conference (CCDC) – Guilin, 2009. 35. Van Essen G.M., Van den Hof, P.M.J. and Jansen JD., A two-level strategy to realize life-cycle production optimization in an operational setting. SPE Journal 18 (6), 10571066, 2012. 36. Vlahavas, I., Vrakas, D., Artificial Intelligence for Advanced Problem Solving Techniques, Information Science Reference, 2008. 37. Zhu, Y., and Backx, T., Identification of Multivariable Industrial Processes for Simulation, Diagnosis and Control, Springer-Verlag, 1993.
Highlights:
An adaptive modeling approach for describing waterflooding in oil reservoirs is proposed. Time-varying nature of reservoirs & inherent nonlinearities can be considered simultaneously. Using system identification techniques, the reservoir can be presented as a MIMO linear plant. Developed Methodology has been applied & evaluated on SPE10 model#2 case study. The obtained linearized adaptive models are suitable for control & optimization purposes.
34
PII: DOI: Reference:
S0920-4105(16)30256-X http://dx.doi.org/10.1016/j.petrol.2016.06.038 PETROL3530
To appear in: Journal of Petroleum Science and Engineering Received date: 4 November 2015 Revised date: 28 April 2016 Accepted date: 22 June 2016 Cite this article as: Farzad Hourfar, Behzad Moshiri, Karim Salahshoor, Mojtaba Zaare-Mehrjerdi and Peyman Pourafshary, Adaptive Modeling of Waterflooding Process in Oil Reservoirs, Journal of Petroleum Science and Engineering, http://dx.doi.org/10.1016/j.petrol.2016.06.038 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Adaptive Modeling of Waterflooding Process in Oil Reservoirs Farzad Hourfara*, Behzad Moshiria, Karim Salahshoorb, Mojtaba Zaare-Mehrjerdia, Peyman Pourafsharyc1 a
b
School of Electrical & Computer Engineering, University of Tehran, Tehran, Iran Department of Automation & Instrumentation Engineering, Petroleum University of Technology, Ahwaz, Iran c School of Chemical Engineering, University of Tehran, Tehran, Iran *
Corresponding Author: Farzad Hourfar, Tel.: +98-912-571-5925. [email protected]
Abstract: Waterflooding is a common solution for enhancing oil recovery in oil reservoirs in which water is injected to the reservoir. For this purpose, injection rates should be determined to optimally control the desired oil production from each well to accomplish reservoir management strategies. This goal can be achieved by precisely modeling the waterflooding phenomena. In this paper, an applicable algorithm for modeling the waterflooding process in oil reservoirs has been developed. Using proper system identification techniques will help to consider any change in the inherent reservoir specifications in addition to operating point variations due to injection alternation rate and also severe nonlinear nature of the reservoir and consequently to obtain adaptive model for the reservoir which is able to reflect all the mentioned phenomena. The reservoir has been considered as a Multi-Input/Multi-Output (MIMO) plant with several manipulated variables and outputs to mimic the variables dynamics of the reservoir during its operational life. A Recursive Least Square (RLS) approach has been utilized with proper forgetting factor and sliding window length selection to extract the reservoir dynamic model in the content of a MIMO model structure. The developed approach has been evaluated for SPE10 model as a well-known reservoir case study in this paper. The observed results demonstrate that the proposed method has acceptable capability, compared to the conventional off-line Least Square (LS) and Extended Least Square (ELS) approaches which may lead to AutoRegressive model with eXogenous input (ARX) and AutoRegressiveMoving-Average model with eXogenous inputs (ARMAX) structures, for modeling the dynamical behaviors of the reservoir to be employed as a powerful tool for control and optimization applications. Keywords: Waterflooding Process; Recursive Least Square Method (RLS); Adaptive Reservoir Modeling; Multi-Input-Multi-Output(MIMO) Systems ; System Identification; Control & Optimization 1
Present Address: Department of Chemical & Petroleum Engineering Sultan Qaboos University, Oman
1
1. Introduction In order to fill the gap between energy demand and supply, it is necessary to enhance oil recovery from existing fields. Nowadays, the average for the recovery factor is about 35% (Rossi et al. 2002), which could decrease to about 15% for complex reservoirs including naturally fractured ones (Sarma, 2006). The production duration life of an oil reservoir generally lasts a few decades and maximum three different production phases may be defined during that period. It should be noted that theoretically between 20%-70% of the oil can be recovered by secondary production techniques (Van den Hof et al. 2009). Oil production process from reservoirs can be considered as a large scale dynamical system. Unfortunately, although an oil reservoir is a complex system, most of the time the operators manage the fields based on their experiences and consequently the complex physics of the system would be ignored. Obviously, besides operating the production based on the operator experience, there are chances to optimize the production by model-based control and optimization techniques, on the basis of predictions of mathematical reservoir models and simulators (Van den Hof et al. 2009). Propitiously, with recent developed technology for drilling and also reservoir instrumentation, there is a capability for utilizing model-based control and optimization approaches to develop applicable reservoir management strategies in order to optimize the economic efficiency of the production process (Van den Hof et al., 2009; Jansen et al. 2008). For choosing the most beneficial scenario of water injection in secondary recovery approach, a comprehensive numerical tool or mathematical model for predicting the performance is necessary (Fanchi, 2001). Reservoir simulation prepares an estimation tool of reservoir performance. There are several techniques of simulation and modeling from simple to complex. The selection of each of the available approaches depends on the accessible data, the level of required accuracy and the calculated time-cost (Ertekin, 2001). Methods such as Decline Curve Analysis (DCA)(Baker, 2003) and the Capacitance Resistance Model (CRM) (Delshad et al., 2009; Moghadam et al., 2011), require less data to model the reservoirs performance. Grid-based simulation (Aziz & Settari 1979; Fanchi, 2001) and streamline simulation (Sayyafzadeh et al., 2010) are the most precise approaches however their simulation time is too much and they need large quantity of data. All the required information for reservoir modeling such as porosity and permeability distributions are not always in hand, in addition they contain high degree of 2
uncertainty; so costly and time-consuming history matching task is required for model validation (Oliver & Chen, 2010). It is clear that an accurate reservoir model is not always necessary for all applications and a fast approach with acceptable accuracy is sufficient. Furthermore, all the required data for a grid-based simulation are not always accessible. As a result, finding a quick modeling technique which needs less data, but has reasonable accuracy is an active research area in the field of reservoir modeling and researchers look for faster methods which need less data, but has reasonable accuracy (Sayyafzadeh, 2011). In Van Essen et al. (2012), Lee et al. (2010), Gildin et al. (2008) and Tafti et al. (2013), suitable proxy/surrogate models for control and optimization applications have been introduced. However, the mentioned references have relied based upon a fixed linear (structure) presentation for the whole range of reservoir operation and consequently this approach will limit their applications within the range of the induced injection rates around the operating point. On the other hand, Using Neural Networks for reservoir modeling has been discussed by Mohaghegh & Abdulla (2014), Aifa (2014), and Bruyelle and Guérillot (2014). However, the obtained models based on neural network are usually in the category of nonlinear models. Therefore, controller or optimizer designing based on neural network methods is not as simple as handling the linear models. In this paper, a systematic approach for adaptive and online modeling of oil reservoirs applicable for the secondary production phase is presented via linear techniques. In spite of the method’s linearity, the results demonstrates that the complicated waterflooding process can be modeled accurately. Based on the introduced modeling mechanism, the field experts will be able to plan for reservoir management and apply advanced control and optimization strategies during the operation in order to achieve to the main goal which is improving the economic performance over the production life of an oil reservoir. 2. Problem statement: Availability of a valid and reliable model is required for any controlling or optimization target in the reservoir management approach. Consequently, obtaining a proper and precise model is a prerequisite for any controlling application.
3
Inherently, a reservoir system model consists of partial differential equations which are based on mass and momentum conservation laws. For solving the model equations numerically, it is needed to discretize the equations in space and time and as a result a large scale system including thousands of grids and even millions of states would be constructed. 2.1. Reservoir Model Reservoir model, based on conservation of mass and momentum is generally used to describe the fluid flow in the reservoir and consequently the dynamics of waterflooding process can be obtained based on the involved equations (Aziz & Settari 1979). By considering oil and water in the reservoir and neglecting the presence of gas, the model equations could be simplified (Jansen et al. 2008). For solving the reservoir equations available in (Aziz & Settari 1979)numerically, it is necessary to discretize them in the space which leads to a system built up of finite number of grids, called as grid blocks. After discretization in time, the following state space representation will be obtained:
V ( xk ) xk 1 T ( xk ) xk qk ,
x0 x0 , (1)
where k is the time index and x is the state vector including the oil pressure po and water saturation S w in all grids. Vector x0 contains the initial conditions, which are assumed to be known. For modeling the effect of the wells on the dynamic of the reservoir, a source vector qk is included in (1). The source terms are usually represented well model, which connects the source term to the pressure difference between the well and grid block pressure:
qkj w j ( pbhj , k pkj ) (2) j
in which pbh , k is the well's bottom hole pressure, j the index of the grid block containing the well j
and pk is the grid block pressure in which the well is located. The term w is a well constant that represents the well's geometric factors and also the rock and fluid properties of the reservoir near the well. The geological properties of each grid block are assumed to be constant. The heterogeneity of the reservoir can be specified by assigning different property values to each of the blocks. Usually a very large number of grid-blocks is needed to describe the dynamic of a 4
real oil reservoir. In the commercial numerical simulators the above equations are solved for all grids in each time step which result massive calculation load. 2.2. MIMO Reservoir Model Structure Generally, the typical oil reservoir models have MIMO structure and in practice it may include 10 to 1000 wells. Furthermore, these models are nonlinear due to the relative permeability effects and also possessing high order because of spatial discretization. The conventional equation-based methodology for reservoir modeling which has been mentioned in section 2.1 is often used in commercial and numerical simulators, leading to a thorough model which can simulate most behaviors of the reservoir. However, in many cases a high percentage of the provided information by the mentioned simulator-based approach is not essential during the actual operation. A hydrocarbon field is practically operated just based on the monitored variables such as total injection and production rates or water and oil production rates of all of the wells. As a result, there is no need for modeling all the possible phenomena in an oil reservoir for different field usages and each application requires only limited and related variables. From the control theory point of view, proper input and output ports should be defined for the system. In the other words, the inputs and outputs of a system are essential to be accessible. The inputs of the reservoir which influence the outputs are called well controls which can be manipulated by flow rate control or BHP control. On the other hand, the main target in reservoir management is controlling the oil production rate of each well and also the produced water. So, the mentioned variables can be considered as the system outputs. Therefore, the reservoir model can be considered as depicted in Figure. 1.
5
Fig. 1. Schematic of oil reservoir as a controlled system with input/output determination 2.3. Application of System Identification in Reservoir Modeling “System Identification” (SI) is a well-known approach for practical modeling of an unknown system for controlling and optimization applications. SI modeling approaches normally consider the plant as a black-box with specified inputs and outputs. In this approach, the ultimate goal is to find a suitable mapping between inputs and outputs through ignoring the internal details of the model for providing the ability of appropriate system outputs estimation in the presence of different and challenging conditions. After finishing model training stage, based on the validated model outputs actions such as prediction, controlling or optimization tasks can be perform. In discrete-time system modeling, one alternative is defining the desired output as a function of past outputs and also all the inputs at previous sampling times. The memory of the model and also the time-delay of the system should be obtained as one of the important stage in the modeling approach. Although in this paper all the results are obtained in discrete-time domain, it should be noted that switching from discrete-time to continuous-time modeling is a straight-
6
forward task which could be easily done according to the required applications (Sung et al., 2009). 2.4. Problem definition Since a petroleum reservoir can be considered as a nonlinear system with different operating points based on the value of the inputs, one of the most applicable solutions is to linearly model its dynamic in the vicinity of each operating point. This approach profits from linear system theories which have strong backbone to be implemented in practice. According to “Parsimony Principle” (Verakas & Vlahavas 2008), for identification of the plant’s behavior, it is recommended to start from the simplest structure and if the model outputs are far from the reality, then improving and upgrade the model structure will be done. As explained previously, in case of considering the reservoir as a MIMO system, the outputs may be defined as oil and water production rate of each production well which can be controlled by two independent inputs named rate of water injection of each injection well and also BHP of the production wells. In practice, it is not an invalid assumption to fix the BHP of each production well as a fixed control input in the desired time interval and consequently the main manipulated variable would become the flow rate of injection wells. At this time the problem can be defined as follows: “Extracting the value of the outputs at fixed BHPs and in presence of injection rate variation”. It is necessary to state that fixing BHP does not have negative effect on the generality of the problem. The proposed method in this paper has the capability to generate the reliable outputs at each BHP and consequently by BHP variation during the operation the updated model will be realized. In addition, by applying MIMO identification procedure which is explained in next sections, general modeling of waterlooding process is feasible.
3. Identification methodology 7
Identification is system modeling based on different experiments. If the tests are not selected correctly then no identification algorithm even the sophisticated ones, can generate a valid model of the process. On the other hand, if the input/output data are obtained from a well-designed experiment, least-squares approach which is the simplest method in this domain can often present a good model. Generally, the purpose of experiments is collecting relevant information about the process dynamics. In the next step, the obtained information should be transformed to valid mathematical models of the process by appropriate identification algorithms. Usually, the collected data from the experiments are not suitable for immediate use in identification algorithms. Possible deficiencies such as high frequencies disturbances; occasional spikes and outliers; drift and offset, low-frequency disturbance; and different order of magnitude may exist. To avoid encountering problems at the estimation steps data pre-processing actions like: Peak Shaving, Trend Correction, Scaling and Offset Correction, Delay Correction, Filtering and Sampling Rate Reduction can be performed (Zhu & Backx, 1993). 3.1 ARX structure modeling approach based on LS After finishing pre-processing task, the first practical step in identifying a physical system is performing offline Least Square (LS) algorithm on the batch input-output data by assuming the system as linear time invariant (LTI). In this method model parameters will be adjusted such that a cost function which is sum of square errors (SSE) be minimized. This approach has been frequently used in industrial processes and successful results have been reported in practical applications (Kon et al. 2013, Chetouani 2008) . A discrete dynamic system can be described by the following equations:
yˆ (k ) T (k 1) (3)
y(k ) yˆ (k ) (k ) (4)
T (k 1) [ y(k 1),..., y(k na ), u(k 1),..., u(k nb )] (5)
8
T (1) (k ) (6) T (k )
T [a1 ,..., ana , b1,..., bnb ] (7) in which u(k) and y(k) are the plant’s input and output at time k (sampling-time), respectively. In reservoir modeling literature, inputs are considered as injection flowrates and bottom hole ^
pressure of wells and outputs are assumed as oil and water production of wells. y (k ) is the model’s output and ε(k) is the model error (difference between the output of the reservoir and the model) and its nature is a zero-mean white Gaussian noise term, which can represent the measurement noise or even the model uncertainties. Regular LS algorithm suppose that the error
T (k 1) between plant output and the model output can be considered as white noise.
is the
vector which contains the available information and it is called the Regressor and Ф (k) is the information matrix which contains the previous information of injection/production values in oil reservoir. is the identification parameter vector which will be determined in the identification process. The parameters in can either be constant or variable. It should be noted that in reservoir modeling parameter The main goal is properly estimation of the values for model parameters: s, hence, all available measurements including inputs and outputs up to the current instant j = 1, 2, …, k, should be utilized in order to minimize the following cost function which is summation of the existing error between model and real output at each sample-time: k
J (k ) (i) (8) 2
j 1
which can be written as:
J ( , k ) j 1[ y( j ) yˆ ( j )]2 (9) k
9
and as a result: J ( , k ) j 1[ y( j ) T ( j 1) ]2 (10) k
Differentiating J(θ,K) with respect to θ and setting the obtained equation equal to zero, the minimum cost function value would be available when the estimate parameter vector is as (Lakshminarayanan, 2001): follows ^
[ j 1 ( j 1) T ( j 1)]1[ j 1 ( j 1) y( j )] (11) k
k
^
where is considered as the least squares estimate of (Ljung & Soderstrom, 1983; Jiang & Zhang, 2004) . This approach results to Auto-Regressive with eXogenous input (ARX) model which has the following format in discrete linear system modeling literature (Zhu & Backx, 1993):
y (k )
B(q 1 ) 1 u (k ) n(k ) (12) 1 A(q ) A(q 1 )
which can be written as:
y (k ) a1. y (k 1) b1.u (k 1)
ana . y (k na ) bnb .u (k nb ) (13)
n( k ) 3.2. ARMAX structure modeling approach based on ELS The above modeling methodology implies that the output of the supposed linear system is a function of previous outputs, inputs and also an unknown term which has the characteristics of the white noise. Hence, it is supposed that the nature of the model error or uncertainty is white noise. However in many cases the dynamic of the error is much complicated and it should be modeled differently which one option is assuming that dynamic as colored noise. In this case, 10
Extended Least Square (ELS) modeling approach can be the next step alternative for being tested and it results to Auto Regressive Moving Average with eXogenous input (ARMAX) model which has the following structure:
B(q 1 ) C (q 1 ) y (k ) u (k ) n(k ) (14) A(q 1 ) A(q 1 ) By considering C(q-1), a new degree of freedom would be added to the model in comparison with ARX model . The estimated error in this approach is defined similar to ARX modeling as:
(k ) y(k ) T (k 1)ˆ(k 1) (15) However, the parameter vector includes c:
[a1
ana b1
bnb c1
cnc ] (16)
and Regressor vector is as follows:
T (t 1) [ y (k 1) u (k 1)
y (k na ) u (k nb )
(k 1) (k nc )]
(17)
and consequently the model parameters can recursively obtained based on (18):
ˆ(k ) ˆ(k 1) P(k ) (k 1) (k ) (18) in which
P1 (k ) P1 (k 1) (k 1) T (k 1) (19) where P(k) is called the covariance matrix. However, since this type of modeling has inherently linear structure, there would be cases in which the obtained model by this approach is not acceptable for critical applications such as controlling the desired outputs or optimization of an objective function.
11
As most of the systems/plants in nature, are nonlinear and also time varying, in some occasions, pure ARX and even ARMAX system identification modeling structures are unable to present the main dynamics of the required system. In this case, other approaches such as Recursive Least Square (RLS) method which is a good choice for practical applications are needed. 3.3 RLS modeling approach
The main idea in the RLS algorithm is re-computing of the parameter in each sample time (i=1,2,…k) by considering a correction term based on new measured data at kth sample. This strategy decreases the computational complexity of RLS and make it suitable for on-line regime and adaptive applications (Jiang & Zhang, 2004). A typical presentation of the RLS is as follows: ^
^
(k ) (k 1) P(k ) (k 1) (k ) (20) in which, ^
(k ) y(k ) T (k 1) (k 1) (21) and also,
P(k 1) (k 1) T (k 1) P(k 1) P(k ) P(k 1) (22) 1 T (k 1) P(k 1) (k 1)
where P(k) is covariance matrix and ε(k) is one-step ahead prediction error. It is obvious that the larger prediction error would cause an increase in the correction term (Niedzwiecki, 2000). In time varying systems, for estimating the new values of the parameter, in addition to recursively calculating the system parameters, a mechanism is needed to ignore the past
12
information that represent the former dynamics which may not be valid at the current time. As a result, techniques such as forgetting factor or finite length sliding window may be used in association with the RLS. By adding term
and L to the cost function, the influences of previous data in parameter
updating would be decreased and the estimated parameters would be find in such a way that (24) be minimized.
J ( , k ) j k L1 k j [ y( j 1) T ( j 1) ]2 (23) k
In (24) L shows the length of the sliding window. In addition, λ which its value is in the interval of [0, 1] is called the forgetting factor for exponentially discarding the obsolete data. Now, the LS problem for parameter estimation of θ is minimizing the modified cost function with the most recent data. In this case, updated parameter vector with utilizing sliding exponentially weighted window least squares algorithm is as follows:
1
k k k j T L (k ) ( j 1) ( j 1) k j T ( j 1) z ( j ) (24) j k L1 j k L1 ^
The performance of the above algorithm will depend on the length of the sliding window. For time-invariant systems the better estimation accuracy will be achieved by increasing the window length. However, for time-variant systems, the window length should be adjusted so that the obsolete information from the past measurements can be discarded properly and heavier weight on the latest measurement be assigned for tracking the system parameter changes. Since in the presented study the system parameters are time-varying due to the nature of the reservoir and also operating point changes, the appropriate selection of window length and forgetting factor value is an important stage (Paleologu, 2008).
13
If the responsibility of discarding old information be assigned to forgetting factor λ then L=k in (26) and RLS algorithm for parameters updating can be written that: ^
^
^
(k ) (k 1) K (k )( y(k ) T (k ) (k 1)) (25) in which:
K (k ) P(k ) (k ) P(k 1) (k )( I T (k ) P(k 1) (k )) 1 (26)
and
P(k ) ( I K (k ) T (k )) P(k 1) / (27) Due to the nonlinear nature of relative permeability and some other parameters of an oil reservoir, an oil reservoir is strongly nonlinear system. In addition, the characteristics of a reservoir would change during oil production and so the reservoir is a nonlinear time-variant plant by the aspect of system theory. Fortunately, RLS technique has the ability to consider both nonlinearity and also time varying nature of the plant simultaneously. In addition, since changing the inputs of a nonlinear system, may cause operating point alternation, the RLS method is also able to locally linearize the plant near the operating point and produce a linear model at each sample time. Furthermore, the parameters of the identified models would be updated in each sample time and consequently any variation in the reservoir specification or nonlinear system’s operating point can be reflected in the new updated parameters. 3.4 MIMO System Identification:
Since, the above mentioned methodologies have been developed for SISO (Single-Input-SingleOutput) plants; in this section the extension to MIMO case would be introduced. This is due to 14
most physical systems in nature including oil reservoir can be considered as MIMO system and several inputs and outputs exist.
A suitable presentation of a MIMO process is a set of difference equations for identification purpose as follows:
y(k ) A1. y(k 1)
An . y(k n) B1.u(k 1)
Bn .u(k n) (28)
in which y is the p-dimensional output vector, u is the m-dimensional input vector and A1 (p×p), ... An (p×p), B1 (p×m) ... Bn (p×m) are constant matrices. Using the unit delay operator q-1,
A(q) y(k ) B(q)u(k ) (29) in which A(q) and B(q) are polynomial matrices:
A(q) I A1q 1
B(q) B1q 1
Anq n (30)
Bn q n (31)
generally, in MIMO Systems ,
y G(q)u (32) Now, the transfer matrix of the system can be written as:
G(q) A(q)1 B(q) (33)
Expressing the transfer matrix by the product of two polynomial matrices is called matrix fraction description (MFD). The parameters of matrices A(q) and B(q) can be estimated from the input/output data. In contrast to the SISO case, for a given process, equation (29) is in general not unique even by fixing the leading coefficient matrix of A(q) at I. Hence, some special forms of representations such as canonical forms are required for obtaining unique models of a given process. Model 15
uniqueness is a necessary condition for an identification algorithm to find a proper solution. Using diagonal canonical form MFD is an acceptable option for MIMO processes. The diagonal form MFD of a given G(q) is
A11 (q) 0 0 A22 (q) A(q) 0 0
B11 (q) , B(q) 0 B p1 (q) App (q) 0
B1m (q) (34) B pm (q)
where A 11 (q), …, App(q) are all monic polynomials and the degrees of Bi1(q), …, Bim(q) are equal to or less than that of Aii(q). with this arrangement, m-input, p-output process is decoupled into p m-input single output sub-processes:
A11 (q) y1 (t ) B11 (q)u1 (t )
B1m (q)um (t ) (35)
App (q) y p (t ) B p1 (q)u1 (t )
B pm (q)um (t )
For physical processes, the transfer matrix can be presented in the diagonal form MFD and consequently most of SISO processes identification algorithms are applicable for MIMO cases in this framework. On the other hand, the diagonal form MFD of a given process is not always minimal. So, the order of state space realization can be higher than the McMillan degree of the process. This problem can easily be solved by model reduction techniques (Benner & Faßbender, 2014).
A11o (q) y1 (t ) B11o (q)u1 (t )
B1m o (q)um (t ) A11o (q ) 1 (t ) (36)
App o (q) y p (t ) B p1o (q)u1 (t )
B pm o (q )um (t ) App o (q ) p (t )
Without loss of generality, it has been supposed that the degrees of all polynomials in a submodel are equal. The orders of sub-models are n1, …nm. The Input/Output data sequence is as follows in which n=max{n1, …, nm}. 16
Z N : y1 (1)
y p (1) u1 (1)
um (1)
y1 ( N n)
y p ( N n) u1 ( N n)
u m ( N n)
(37) Assuming that N samples of measurements of outputs and inputs are made at time 1, 2,…, N. Now, the loss function which should be minimized in order to estimate the required parameters is:
V1
1 N n ( A1 (q) y1 (t ) [ B11 (q)u1(t ) N t n1
B1m (q)um (t )]
(38)
and consequently,
ˆ1 [
1 T 1 1 1 ]1 1T y1 (39) N N
where
y1 (n 1) y (n 2) , y1 1 y1 (n N )
y1 (n1 ) y (n 1) 1 1 1 y1 (n1 N )
a1,1 a1,n 1 b11,1 (40) 1 b11,n 1 b 1m ,1 b1m ,n1 y1 (1) | u1 (n1 )
|
| um (n1 )
| u1 (n1 1)
|
| um (n1 1)
|
|
|
y1 ( N 1)| u1 (n1 N )
u1 (1)
u1 ( N 1) |
17
| um (n1 N )
(41) um ( N 1) um (1)
The same procedure can be performed for other sub- models. In addition, it has been supposed that the degrees of all polynomials in a sub-model are equal. By utilizing the above presented approach, it will be possible to extract a MIMO model for any type of oil reservoir with several inputs and outputs. However, although there exists different techniques for directly MIMO modeling of real systems, simplifying the problem to several MISO identification can be a proper alternative for modeling the system (Ding & Chen, 2005a; Ding & Chen, 2005b; El-Anes, 2013).
4. Methodology of Modeling Oil Reservoir: At this section the developed methodology for modeling oil reservoirs is introduced. This technique can be break-down to three sub-stages as follows: 4.1. Finalizing the General Structure In this phase, the source of data would be defined. The data source may either be the real reservoir data which have been gathered from the field or data which have been generated by a valid commercial or academic reservoir. In addition, all manipulated variables which can act as the system inputs and also the desired outputs should be fixed at this step. The main goal of the algorithm is assigning reliable model/models which can demonstrate the behavior of the system appropriately. 4.2. Data Generation As explained in the previous sections, recording proper input-output data is one of the most vital phases in system modeling and identification. The richer collected data would result to the better obtained model. So, it is necessary to stimulate the plant (here; the reservoir) such that all the hidden dynamics become apparent. For this purpose, Persistently Exciting (PE) signals should be applied as system inputs. A good choice of rich inputs is white noise or even Pseudo Random Binary Signal (PRBS). However, since in practice, applying those inputs to the reservoir is 18
somehow impossible and the setting of the control valves cannot be modified with high frequencies, an appropriate signal which has the most similarity to PRBS and also meets the operational constraints should be applied to the system and all inputs/outputs be recorded simultaneously. In the next step before starting the main identification process, according to the above explanations, data pre-processing is necessary to be performed in order to enhance the numerical robustness of identification algorithms and also accelerate the speed of convergence. 4.3. Model Identification In the final step, the proper identification algorithm should be utilized for achieving to the suitable model. According to the “Parsimony Principle” in system modeling and identification, the candidate structures starts from the simplest ones and if there exists a large discrepancy between the model outputs and real plant outputs, a more complicate structure will be examined. So, in the proposed methodology, the start point is evaluating the obtained results from an offline batch ARX structure based on least square technique and if the results are unsatisfying then ARMAX structure and adaptive RLS will be tested respectively. It is important to notice that it is helpful to initialize the parameters of more complicated structures based on the obtained results of the former tested strategy and this approach will significantly increase the convergence speed of the algorithm to find the proper model. The flowchart of this methodology has been illustrated in Figure 2.
19
Fig. 2. Flowchart of the Methodology.
5. Selected Reservoir Case Study The selected model 10th SPE-Model#2 has a simple geometry, with no faults and top structure however, it is one of the well-known case study in the field of reservoir researches. At the fine geological model scale, the model is described on a regular Cartesian grid. The model dimensions are 1,200×2,200×170 ft. The top 70 ft (35 layers) represents the Tarbert formation, and the bottom 100 ft (50 layers) represents Upper Ness. The fine-scale cell size is 20×10×2 ft which results that the model has 60×220×85 grid-cells.. The model consists of part of a Brent sequence. The model was originally generated for use in the PUNQ project. Although in the original model a uniform kV/kH in the whole domain is considered, the used model has been changed such that kV/kH is equal to 0.3 in the channels and 10-3 in the background. The porosity and permeability maps and other model parameters are available in the SPE homepage (Christie & Blunt, 2001; Islam & Sepehrnoori, 2013). 20
In this synthetic model, all wells are vertical and completed throughout the whole formation. In the original model the central injector has an injection rate of 5000 Barrles/ Day. There are four producers in the four corners of the model; In the initial example each well produces at 4,000 psi bottom-hole pressure. The X-Y location of each well is presented in Table 1 in Cartesian coordinate: Table 1. Well Locations in Cartesian coordinate. well name
X location
Y location
(ft)
(ft)
Injection Well I1
600
1100
Production Well P1
0
0
Production Well P2
1200
0
Production Well P3
1200
2200
Production Well P4
0
2200
6. Implementation of the presented algorithm and results: For implementing of the developed algorithm, Matlab reservoir simulation toolbox (MRST) has been used (Lie, 2014). This simulator provides the facility to control each well by adjusting flowrate or BHP. Since, the ultimate goal of system modeling is controlling the desired outputs or even optimization the whole reservoir performance, the manipulated variables would be selected as the injection rate of injection well(s). In addition, it is supposed that the production wells are running under the fixed BHP due to the operational recommendations. Furthermore, the output of the system which are going to be identified can be selected as oil production rate and water production rate of producing wells. In order to maintain the pressure of the reservoir constant, it is supposed that the total injection rates of injection wells are equal to the total production rates of producing wells. 21
m
p
i 1
j 1
qi _ inj q j _ prod (42) Where qi_inj is the flowrate of each injection well (here, flow rate of injected water) and qj_prod is the flowrate of each producing well (here, summation of produced oil and water). When the inputs are set, the simulator can simultaneously calculate and generate total production rate, oil production rate, water production rate and water cut of each well. The water cut is the ratio of water produced compared to the volume of total liquids produced from an oil well.
It is obvious that availability of only two of the mentioned variables, may lead to calculate the others. In practice, each of them can be measured uniquely by applying proper instruments in the line. Now, step by step the proposed algorithm is applied to the well-known 10th SPE-Model#2 for identification of the waterflooding process in an oil reservoir supposing that the BHPs of the production wells are constant which conclude to a Single Input, Multi Output (SIMO) system (one injection well and four production wells). In simulation, 2000 days of reservoir operation with a proper injection rate properties have be analyzed. For stimulating the dynamic of the reservoirs and also for identification the hidden modes, a PE input is required. However, in reservoir management, it is not suitable to change the value of each manipulated variable continuously and in specified times which may be dictated by the experts those values would be changed. So, a proper candidate input sequence which is also applicable in practice and has the most analogy with a pure PE signals can be seen in table 2. Table 2. Sequence of Injection rate Injection # 1 2 3 4 5 6 7 8
Injection value (bbls/day) 5000 4500 6000 4000 5000 6000 3000 5500 22
Starting day
Ending day
1 201 301 501 801 901 1001 1301
200 300 500 800 900 1000 1300 1500
3500 4500 5000
9 10 11
1501 1601 1801
1600 1800 2000
In this case, it has been supposed that the injection rate may only change in each 100 days randomly (10 samples) for having the most similarity to the practical experiments. In addition, it is assumed that the values of injection rates are between 3000bbl/day and 6000bbl/day. Since an oil reservoir is a nonlinear system, these values would help to identify its different operating points and also facilitate to identify the reservoir’s dynamic close to each operating point. The output of the system has been monitored and recorded every 10-days and base on the introduced algorithm, waterflooding modeling and identification procedure has been started. The results show that ARX and ARMAX modeling approaches are unable to precisely model the dynamic of the reservoir, since they try to fit batch models for whole production life. The best obtained result for ARX model with structure in (43), is in (44):
A(q) y(k ) B(q)u(k ) e(k ) (43) and
A(q) 1 q 1 0.0641q 2 (44)
B(q) 0.0298q 2 which tries to model the relationship between total production of well #3 and injected water. Figure 3. Illustrates a sample result for evaluating the performance of ARX structure for estimating the total production of Well #3. It is clear that although the output of ARX model can somehow track the trend of the real value of the production rate, it is far from the desired simulator’s output.
23
Fig. 3. Performance of the ARX model output versus simulator’s output. In the next step, ARMAX structure should be evaluated. The best result of this model structure for the same input/output variable is:
A(q) 1 1.657q 1 0.6816q 2 B(q) 0.01143q 2
(45)
C (q) 1 0.657q 1 0.0001q 2
With the following structure:
A(q) y(k ) B(q)u(k ) C (q)e(k ) (46) Although Figure 4. presents that the performance of this model is slightly better than ARX model, however it cannot be considered as a reliable waterflooding model during the production period.
24
Fig. 4. Performance of the ARMAX model output versus simulator’s output.
The quality of obtained results is predictable since modeling a strongly nonlinear plant with linear approaches in the whole domain is theoretically impossible. However, RLS approach can be an appropriate candidate for extracting and identification of the behavior of the system and predict the required output. By applying this technique, at each sampling time (in this case, every 10 days), a new localized linear model will be generated based on the available data up to that time and consequently any parameter or structural change in the system can be detected and be reflected in the new model. As the following results presents, by adjusting the forgetting factor of RLS to 0.95, this technique can identify all of the required outputs for all production wells precisely and the results obtained of the RLS strategy are very similar to the outputs of the simulator. It is important to notice that in this condition, for each sample time (each 10 days) a new ARX model is valid for indicating the behavior of the plant which its parameters have been also illustrated in the figures. Now this set of models is valid for representing the whole production period behaviors and it has also the ability for detecting any variation in the reservoir’s characteristics. It can be seen in most of the results that any change in the input, which may cause the operating point change in the plant can cause the model’s parameter variation. Furthermore, this approach of waterflooding 25
modeling has the capability to regenerate the new models whenever the BHP of the production wells also alter and its application is not limited to only constant BHP cases. Table 3. shows the value of sum of square error (SSE) between simulator output and each of the 3 modeling strategies’ output related to the total production rate of well#3.
Table 3. Well Locations Modeling approach
SSE
ARX
7.9050e7
ARMAX
6.6519e7
RLS
3.6437e4
Figures 5, 6, 7 & 8 illustrate the simulator outputs and generated models’ outputs in which the outputs are defined as total production rate of the wells and system input is the water injection rate. For another example, Figures 9, 10, compare the estimated and real oil production rates for two selected wells. In figures 11, 12, the estimated produced water of well #2 & well #4 can be observed and finally figure 13 showes the water cuts of well #1 for instance. In addition, the estimated model parameters for all desired outputs have been depicted in each picture in which ai ‘s are the denominator and bi’s are the nominator coefficients of the dedicated model according to (12) and as a result in each time step a unique model is available.
26
Fig. 5. Injection rate, Total Production Rate & Model coefficients related to well #1.
Fig. 6. Injection rate, Total Production Rate & Model Coefficients related to well #2.
27
Fig. 7. Injection rate, Total Production Rate & Model Coefficients related to well #3.
Fig. 8. Injection rate, Total Production Rate & Model Coefficients related to well #4.
28
Fig. 9. Injection rate, Oil Production Rate & Model Coefficients related to well #1.
Fig. 10. Injection rate, Oil Production Rate & Model Coefficients related to well #3.
29
Fig. 11. Injection rate, Water Production Rate & Model Coefficients related to well #2.
Fig. 12. Injection rate, Water Production Rate & Model Coefficients related to well #4.
30
Fig. 13. Injection rate, Water cut & Model Coefficients related to well #1.
7.
Conclusion:
In this paper an adaptive identification algorithm has been developed for modeling the waterflooding process in oil reservoirs. The presented approach has the ability to reflect the timevarying nature of the reservoir with its inherent nonlinearity. The developed approach, following the parsimony principle conducts a step by step model scheme to adaptively identify a simple candidate model structure through a preliminary off-line ARX model. Then, the approach goes through more complex ARMAX model structure in an offline manner if the obtained model accuracy is not enough for application target. An adaptive online model will subsequently adopted to capture the reservoir dynamics in case the offline batch model is not adequate to accomplish the required model accuracy. The observed results on SPE10 benchmark case study have shown that the RLS approach has the ability to estimate each of the required reservoir outputs with acceptable accuracy. Forgetting factor parameter adjustment in RLS method adds the ability of discarding old data for obtaining the most reliable model in each time. Since the presented technique is generally developed for MIMO plants, it can easily include any type of oil reservoir with several injection/production wells and also adjusting flowrate or BHP of each well simultaneously. The ultimate objective of this kind of modeling is to generate proper models for 31
control and optimization goals. As a result, availability of such adaptive models enables to apply useful adaptive control approaches for desired output(s) and also perform optimization tasks on the reservoir during the operation
References: 1. Aifa, T., Neural network applications to reservoirs: Physics-based models and data models, Journal of Petroleum Science and Engineering, 123, pp.1-6, 2014. 2. Aziz, K., Settari, A.N., 1979. Petroleum reservoir simulation. London, Applied Science Publishers. 3. Baker, R.O., Anderson, T., Sandhu, K., Using Decline Curves to Forecast Waterflooded Reservoirs: Fundamentals and Field Cases. Canadian International Petroleum Conference. Petroleum Society of Canada, Calgary, Alberta. 2003. 4. Benner, P., Faßbender, H., Model Order Reduction: Techniques and Tools, Encyclopedia of Systems and Control, Springer-Verlag, 2014. 5. Bruyelle, J., Guérillot, D., Neural networks and their derivatives for history matching and reservoir optimization problems, Computational Geosciences, Volume 18, Issue 3, pp 549-561, 2014. 6. Chetouani, Y., Using ARX approach for modelling and prediction of the dynamics OFA reactor-exchanger, Symposium Series NO. 154, IChemE 2008. 7. Christie, M. A., and M. J. Blunt. Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques. SPE 66599, SPE Reservoir Simulation Symposium, Houston, Feb 11‐14, 2001. 8. Delshad, M., Bastami, A., Pourafshary, P.,The Use of Capacitance-Resistive Model for Estimation of Fracture Distribution in the Hydrocarbon Reservoir. SPE Saudi Arabia Section Technical Symposium. Society of Petroleum Engineers, AlKhobar, Saudi Arabia. 2009. 9. Ding, F., Chen, T., Hierarchical gradient-based identification of multivariable discretetime systems, Automatica 41: 315 – 325, 2005. 10. Ding, F., Chen, T., Hierarchical Least Squares Identification Methods for Multivariable Systems, IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 50, NO. 3. 2005. 11. El-Anes, A., Bouzrara, K., Garna, T., Messaoud, H., Expansion of MIMO ARX model on Laguerre orthonormal bases. International Conference on Electrical Engineering and Software Applications (ICEESA), 2013. 12. Ertekin, T., Abou-Kassem, J.H., King, G.R., Basic applied reservoir simulation. Society of Petroleum Engineers, Richardson, Tex. 2001. 13. Fanchi, J.R., Principles of Applied Reservoir Simulation, Second Edition. Gulf Professional Publishing, Burlington. 2001. 14. Gildin, E., Wheeler, M.F., Control of Subsurface Flow Using Model Predictive Control Techniques, in Proceedings of the International Conference on Engineering Optimization, Brazil, June, 2008. 32
15. Islam, A. W., Sepehrnoori, K. A Review on SPE’s Comparative Solution Projects (CSPs), Journal of Petroleum Science Research (JPSR) Volume 2 Issue 4, October 2013. 16. Jansen, J.D, Bosgra, O.H, and Van den Hof, P.M.J, Model- based control of multiphase flow in subsurface oil reservoirs. Journal of Process Control, 18(9):846-855, 2008. 17. Jiang, J., Zhang, Y., A revisit to block and recursive least squares for parameter estimation. Computers & Electrical Engineering; 30(5): 403-416. 2004. 18. Kon, J., Yamashit, Y., Tanaka, T., Tashiro, A. and Daiguji, M. Practical Application of Model Identification based on ARX Models with Transfer Functions, Proceedings of the 6th International Conference on Process Systems Engineering (PSE ASIA), 2013. 19. Lakshminarayanan, S., Emoto, G., Ebara, S., Tomida, K., and Shah, S.L., Closed loop identification and control loop reconfiguration: An industrial case study. Journal of Process Control 11: 587–599, 2001. 20. Lee, K.H., Ortega, A., Jafroodi, N., & Ershaghi, I., A Multivariate Autoregressive Model for Characterizing Producer-producer Relationships in Waterfloods from Injection/Production Rate Fluctuations. Society of Petroleum Engineers, 2010. 21. Lie, K. A., An Introduction to Reservoir Simulation Using MATLAB: User guide for the Matlab Reservoir Simulation Toolbox (MRST). SINTEF ICT, May 2014. 22. Ljung, L., and Soderstrom, T., Theory and practice of recursive identification. MIT Press, Cambridge. MA. 1983. 23. Moghadam, J.N., Salahshoor, K., Kamp, A. M., Evaluation of Waterflooding Performance in Heavy Oil Reservoirs Applying Capacitance-Resistive Model, Petroleum Science and Technology, 29:17, 1811-1824, 2011. 24. Mohaghegh, Sh., Abdulla, F., Production Management Decision Analysis Using AIBased Proxy Modeling of Reservoir Simulations – A Look-Back Case Study, Society of Petroleum Engineers, 2014. 25. Niedzwiecki, M., Identification of time-varying processes. Wiley 2000, Chichester. 26. Oliver, D., Chen, Y., Recent progress on reservoir history matching: a review. Computational Geosciences 1–37. 2010. 27. Paleologu,C., Benesty, J., Ciochina,S., A Robust Variable Forgetting Factor Recursive Least-Squares Algorithm for System Identification, IEEE SIGNAL PROCESSING LETTERS, VOL. 15: 597-600, 2008. 28. Rossi, D., Carney, M., Kontchou, J.N., Lancaster, D., McIntyre, S., Reservoir and Production Optimization, white paper from www.slb.com, 2002. 29. Sarma, P., Efficient Closed-Loop Optimal Control of Petroleum Reservoirs under Uncertainty, PhD Thesis, Stanford University, Stanford, CA, 2006. 30. Sayyafzadeh, M., Pourafshari, P., Rashidi, F. Increasing ultimate oil recovery by infill drilling and converting weak production wells to injection wells using streamline simulation. International Oil and Gas Conference and Exhibition in China. Society of Petroleum Engineers, Beijing, China., 2010. 31. Sayyafzadeh, M., Pourafshary, P., Haghighi, M., Rashidi, F., Application of transfer functions to model water injection in hydrocarbon reservoir, Journal of Petroleum Science and Engineering 78 ,139–148, 2011. 32. Sung, S.w., Lee, J., Lee, I. B., In-Beum Lee, Process Identification and PID Control, John Wiley & Sons, 2009.
33
33. Tafti, T., Ershaghi, I., Rezapour, A., & Ortega, A. Injection Scheduling Design for Reduced Order Waterflood Modeling. Society of Petroleum Engineers, 2013. 34. Van den Hof, Paul M.J.; Jansen, Jan-Dirk; Essen, Gijs van; Bosgra, Okko H. Modelbased control and optimization of large scale physical systems - Challenges in reservoir engineering, IEEE Chinese Control and Decision Conference (CCDC) – Guilin, 2009. 35. Van Essen G.M., Van den Hof, P.M.J. and Jansen JD., A two-level strategy to realize life-cycle production optimization in an operational setting. SPE Journal 18 (6), 10571066, 2012. 36. Vlahavas, I., Vrakas, D., Artificial Intelligence for Advanced Problem Solving Techniques, Information Science Reference, 2008. 37. Zhu, Y., and Backx, T., Identification of Multivariable Industrial Processes for Simulation, Diagnosis and Control, Springer-Verlag, 1993.
Highlights:
An adaptive modeling approach for describing waterflooding in oil reservoirs is proposed. Time-varying nature of reservoirs & inherent nonlinearities can be considered simultaneously. Using system identification techniques, the reservoir can be presented as a MIMO linear plant. Developed Methodology has been applied & evaluated on SPE10 model#2 case study. The obtained linearized adaptive models are suitable for control & optimization purposes.
34