Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Investigating statistical modeling approaches for reservoir characterization in waterfloods from rates fluctuations
(USC Thesis Other)
Investigating statistical modeling approaches for reservoir characterization in waterfloods from rates fluctuations
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
INVESTIGATING STATISTICAL MODELING APPROACHES
FOR RESERVOIR CHARACTERIZATION IN WATERFLOODS
FROM RATES FLUCTUATIONS
by
Kun-Han Lee
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(ELECTRICAL ENGINEERING)
December 2010
Copyright 2010 Kun-Han Lee
Dedication
Po-Hsien Lee
Yun-Ru Chen
Yun-Fang Lu
ii
Acknowledgements
I would like to express my sincere appreciation and gratitude to my supervisor, Dr.
Antonio Ortega, for his support throughout my graduate studies at the University
of Southern California. I have been fortunate to have an advisor who gave me the
freedom to explore on my own, and at the same time the guidance to recover when
my steps faltered. I am grateful to his guidance, patience, and encouragement.
I owe my gratitude to Dr. Iraj Ershaghi, for his support of many useful discus-
sions and suggestions. He has made many insightful comments from the petroleum
engineering aspects at different stages of my research. Without him my thesis would
not have been completed. Thanks are also due to Dr. Mendel for serving as mem-
bers of my dissertation committee and Dr. Kuo and Dr. Jenkins for serving as
committee members of my qualify exam. Their valuable comments and suggestions
were really appreciated.
I would like to thank to several PhD students in Petroleum Engineering: Nelia
Jafroodi, Amir Mohammad Nejad, Mohammad Javaheri and Reza Rastegar Moghadam.
They cooperate with me to finish all the numerical simulations and make many valu-
able discussions. I am also indebted to the staff of the CiSoft project for their help
and continuous support.
Many friends have helped me stay sane through these difficult years. Their
support and care helped me overcome setbacks and stay focused on my graduate
study. I greatly value their friendship.
iii
Most importantly, none of this would have been possible without the love and
patience of my family. They to whom this dissertation is dedicated to, has been a
constant source of love, concern, support and strength all these years. I would like
to express my heart-felt gratitude to my family.
iv
Table of Contents
Dedication ii
Acknowledgements iii
List Of Tables viii
List Of Figures ix
Abstract xi
Chapter 1: Introduction 1
1.1 Background . ...... ..... ..... ..... ...... ... 1
1.2 ANewResearchTrend . ..... ..... ..... ...... ... 3
1.3 ResearchContributions. ..... ..... ..... ...... ... 5
Chapter 2: New Predictive Models 8
2.1 LinearTime-invariantSystem . . ..... ..... ...... ... 9
2.2 RevisitingtheCapacitanceModel ..... ..... ...... ... 10
2.3 FiniteImpulseResponseModel . ..... ..... ...... ... 13
2.3.1 DiscreteModel . ..... ..... ..... ...... ... 14
2.4 DistributedCapacitanceModel . ..... ..... ...... ... 15
2.4.1 InterpretationofModel.. ..... ..... ...... ... 20
2.5 MultivariateARXModel..... ..... ..... ...... ... 22
2.5.1 DiscreteModel . ..... ..... ..... ...... ... 26
2.6 FindingtheModelParameters.. ..... ..... ...... ... 26
2.7 Summary .. ...... ..... ..... ..... ...... ... 29
Chapter 3: Model Validation and Comparison 30
3.1 ValidationofModels . . ..... ..... ..... ...... ... 31
3.2 SimulationSetting ... ..... ..... ..... ...... ... 33
3.2.1 Permeability Map..... ..... ..... ...... ... 35
3.2.2 InjectionRateSetting .. ..... ..... ...... ... 38
3.2.3 OtherSetting .. ..... ..... ..... ...... ... 39
v
3.3 ValidationResultsforInterwellConnectivity.... ...... ... 40
3.3.1 QuantitativeVerification . ..... ..... ...... ... 43
3.4 PhysicalConstraintsonParameters .... ..... ...... ... 44
3.5 ModelComparison ... ..... ..... ..... ...... ... 50
3.5.1 NumberofParameters . . ..... ..... ...... ... 50
3.5.2 PredictionPerformanceComparisons .... ...... ... 51
3.6 ComparisonResults... ..... ..... ..... ...... ... 52
3.6.1 NoiseSensitivityAnalysis ..... ..... ...... ... 57
3.7 Discussions . ...... ..... ..... ..... ...... ... 59
3.7.1 DrawbacksofPurelyStatisticalMethods .. ...... ... 59
3.7.2 LargeScaleProblems... ..... ..... ...... ... 60
3.7.3 SuggestionsforUsingPredictiveModels .. ...... ... 62
3.8 Summary .. ...... ..... ..... ..... ...... ... 63
Chapter 4: Linear Modeling Framework 64
4.1 LinearModelsforWellInteractions .... ..... ...... ... 64
4.2 AFrameworkforPredictiveModels .... ..... ...... ... 66
4.2.1 CapacitanceModel .... ..... ..... ...... ... 67
4.2.2 DistributedCapacitanceModel .. ..... ...... ... 68
4.2.3 DoublePoleModel[39].. ..... ..... ...... ... 69
4.2.4 Finite-Impulse-ResponseModel .. ..... ...... ... 70
4.2.5 MultivariateARXModel. ..... ..... ...... ... 72
4.2.6 HigherOrderModels ... ..... ..... ...... ... 73
4.3 Conclusion. . ...... ..... ..... ..... ...... ... 74
Chapter 5: Prediction Under Controlled Producers 75
5.1 InterpretingM-ARXModel ... ..... ..... ...... ... 76
5.2 PredictionforShut-InProducers ..... ..... ...... ... 78
5.2.1 SimulationResults .... ..... ..... ...... ... 80
5.3 PredictionforConstrainedProducers ... ..... ...... ... 82
5.3.1 SimulationResults .... ..... ..... ...... ... 84
Chapter 6: Injection Rates Design 90
6.1 LiteratureReview.... ..... ..... ..... ...... ... 91
6.1.1 Optimal Input Design . . ..... ..... ...... ... 91
6.1.2 ChannelEstimationinCommunication ... ...... ... 93
6.2 A Novel Deterministic Approach for Input Sequence Design . . . . . 95
6.2.1 Inverse-repeatsignals ... ..... ..... ...... ... 95
6.2.2 PropertyofaSetofInverse-RepeatSignals. ...... ... 96
6.2.3 DesignExample. ..... ..... ..... ...... ... 98
6.2.4 Discussion .... ..... ..... ..... ...... ... 99
6.3 StochasticApproach . . ..... ..... ..... ...... ... 101
6.3.1 CriterionforInjectionRateDesign ..... ...... ... 102
vi
6.3.1.1 MISOcase ... ..... ..... ...... ... 104
6.3.1.2 MIMOcase ... ..... ..... ...... ... 105
6.3.1.3 TheCrestFactor ..... ..... ...... ... 106
6.3.2 ADesignProcedure ... ..... ..... ...... ... 107
6.3.3 DesignExample: ApplicationtoCM .... ...... ... 109
6.3.4 SimulationResults .... ..... ..... ...... ... 112
6.3.5 ComparisontoDeterministicApproach ... ...... ... 116
6.4 Conclusion. . ...... ..... ..... ..... ...... ... 118
Chapter 7: Conclusions and Future Work 119
7.1 Conclusions . ...... ..... ..... ..... ...... ... 119
7.2 FutureWork. ...... ..... ..... ..... ...... ... 121
References 123
vii
List Of Tables
3.1 SimulationSettingsforCMGsimulator. . ..... ...... ... 34
3.2 Estimated values of interwell connectivity for Scenario A. . . . . . . 40
3.3 Average absolute differences of interwell connectivities between CM
andothermodelsforScenarioA. ..... ..... ...... ... 41
3.4 Estimated values of interwell connectivity for Scenario B. . . . . . . 42
3.5 Average absolute differences of interwell connectivities between CM
andothermodelsforScenarioB. ..... ..... ...... ... 43
3.6 Estimated interwell connectivities for peripheral injector pattern (Sce-
narioE,F,andG). ... ..... ..... ..... ...... ... 45
3.7 Numberofparametersusedfordifferentmodels. . . ...... ... 51
4.1 Different predictive models and the characteristics of their transfer
function. ... ...... ..... ..... ..... ...... ... 74
5.1 The number of parameters needed to be retrained when a producer
isshut-in. .. ...... ..... ..... ..... ...... ... 81
6.1 Non-zeroDFTindexesfroinverse-repeatsignals... ...... ... 97
viii
List Of Figures
2.1 TheimpulseresponseofCapacitanceModel. .... ...... ... 11
2.2 Analogous RC circuit for capacitance model. .... ...... ... 12
2.3 Representative element for distributed capacitance model. There is
a high permeability channel between the injector and producer. . . 17
2.4 Impulse response curves of distributed capacitance model for differ-
ent τ
1
and τ
2
. The injected water of injector is also shown. . . . . . 19
2.5 Interpretation of distributed capacitance model in terms of CMs. . . 21
2.6 Scenario with multiple injection and production wells. The red shadow
areas represent the pore volume related to that producer. . . . . . . 23
2.7 DiagramofthemodelingprocessofM-ARXmodel. ...... ... 27
3.1 Permeability map for channel case (Scenario A): channel case. . . . 34
3.2 Permeability map for streak case (Scenario B). . . . ...... ... 35
3.3 Permeability map for homogeneous case (Scenario C)...... ... 35
3.4 Permeability map for multiple fractures case (Scenario D). . . . . . 36
3.5 Permeability map for Scenarios E, F, and G..... ...... ... 36
3.6 InjectionratesforScenariosA-D. ..... ..... ...... ... 37
3.7 InjectionratessetforScenarioF. ..... ..... ...... ... 38
3.8 InjectionratesassignedforScenarioG. .. ..... ...... ... 39
3.9 EstimatedinterwellconnectivitiesforScenarioA. . ...... ... 41
3.10 EstimatedinterwellconnectivitiesforScenarioB. . ...... ... 43
ix
3.11 Estimated interwell connectivities for peripheral injector pattern (Sce-
narioE,F,andG). ... ..... ..... ..... ...... ... 46
3.12 Performance evaluation of different models on the Scenario C . . . . 53
3.13 Performance evaluation of different models on the Scenario D. . . . 54
3.14 Performance evaluation of different models for Scenario A. . . . . . 55
3.15 Performance evaluation of different models for Scenario G. . . . . . 56
3.16 Noise sensitivity analysis for different models for Scenario A. . . . . 58
5.1 The CRMP and M-ARX model in 2-injectors/2-producers scenario. 76
5.2 Performance of prediction error for production rates based on M-
ARX model in Scenario A with P4shut-in...... ...... ... 82
5.3 Prediction of production rates by M-ARX model in Scenario A with
P4shut-in.. . ...... ..... ..... ..... ...... ... 83
5.4 PredictionresultsforConstrainedP1.... ..... ...... ... 85
5.5 Prediction results for Constrained P2. .. ..... ...... ... 86
5.6 Prediction results for Constrained P3. .. ..... ...... ... 87
5.7 Prediction results for Constrained P4with600bbl/day..... ... 88
5.8 Prediction results for Constrained P4 with 2500 bbl/day. . . . . . . 89
6.1 Illustration of non-zero frequency indexes for inputs with different
periods. ... ...... ..... ..... ..... ...... ... 97
6.2 Anexamplefordesignedsetofinverse-repeatsignals. ..... ... 99
6.3 Block diagram for the procedure which generates the set of inputs
forinjectionratedesign. ..... ..... ..... ...... ... 108
6.4 A realization of designed injection rates used for performance evalu-
ation...... ...... ..... ..... ..... ...... ... 114
6.5 Performance plot for injection rates design with different sequences. 115
6.6 Performance plot for injection rates design with both deterministic
andstochasticapproaches. .... ..... ..... ...... ... 117
x
Abstract
Reservoir characterization is important for reservoir management and performance
optimization. For waterflood optimization, traditionally several techniques have
been suggested, most of which are either too time-consuming or the data needed
are often unavailable. There is a new research trend to overcome these limitations
by applying advanced statistical techniques on only injection and production data,
which are often readily available for any waterflood operations.
In this work, we follow this trend to formulate the reservoir characterization and
forecasting problems in waterflood projects using a system identification framework:
the injection rates are seen as inputs; the production rates are seen as the outputs;
and the reservoir is considered as a dynamic system. By addressing the properties
of general linear dynamic systems, we discuss the limitations of previous models and
build three new predictive models to characterize some reservoir behavior, such as
producer-to-producer interactions, which was not considered in previous literature.
Then we discuss a general parameter estimation approach under the prediction-
error framework.
For model evaluation, we propose two techniques: one is based on evaluating
their prediction ability on a fresh data set, while the other is based on comparing the
interpretations they provide about certain reservoir characteristics with the ground
truth. All proposed models are verified by these two approaches. To perform
a comparative analysis, we provide a practical metric to compare the prediction
xi
performance of different proposed models under various scenarios. From the results,
we make several observations and suggestions for reservoir engineers to use the
models.
To clarify the relationship between different models, we develop a general linear
modeling framework and demonstrate that all proposed models can be considered as
special cases within this framework. Moreover, the transfer function of the general
linear model can be interpreted to provide insight on reservoir characteristics. Also,
the relationship between different models can easily be built from this work.
We propose a multivariate autoregressive model to characterize situations in
which a producer is shutting-in or a new producer is being brought online. As
a totally new application, we introduce a novel “constrained producer” approach
which that only requires minimal changes in production rates (e.g., limiting them
to some level below their normal production capacity) to predict the performance
after a producer is shut-in. This allows us to handle various “what if” scenarios in
waterflood management.
Finally, to achieve a better model estimation, the patterns of injection rates
play an important role. We addressed the problem of designing a set of injection
rates to achieve a better estimation of target parameters in the reservoir. Two
different approaches, deterministic and stochastic approaches, are discussed. For
the deterministic approach, we propose a new procedure using a set of inverse-
repeat signals to design a set of signals with zero cross-correlation property. For
stochastic approach, we applied a common approach in system identification and
evaluate all design procedures on some predictive model.
xii
Chapter 1
Introduction
1.1 Background
Waterflooding is a common operation in many petroleum reservoirs. It refers to a
process where some wells, denoted as injectors, inject water into the reservoir in
order to increase the reservoir pressure and displace the oil to some surrounding
wells, which denoted as producers. The pressure gradient between injectors and
producers helps in increasing oil production.
To maximize oil recovery by waterflooding, it is useful to “understand” the un-
derground structure of the reservoirs. Traditionally, there are several approaches to
estimate the reservoir characteristics among wells. Among these methods, multiple-
well tests (interference and pulse tests [27]), a type of pressure-transient test ( [53]
[24]), have been developed to establish communication between wells and deter-
mine the interwell reservoir properties. Basically these approaches rely on matching
pressure response from wells to a theoretical model of reservoir. Using a nonlinear
parameter estimation, some crucial reservoir characteristics can be approximated.
Refer to Kamal [29] for a review on this topic. For these tests, the costs are signifi-
cant because the amplitudes of responses are small (sometimes less than 1 psi) and
1
precise pressure measurements are needed to get meaningful data. On the other
hand, downhole pressure recorders data are not available for most of the field op-
erations. This becomes an impediment for routine application of these techniques.
As a alternative, tracer tests have widely been used for mapping flow commu-
nication among wells. Tracer studies can also be useful for estimation of reservoir
parameters, such as reservoir swept volume, fluid velocities, and flow geometry in a
reservoir. Basically this method requires injecting some chemicals (tracers) in some
wells and observing the tracer concentration curve on the produced fluid from the
surrounding wells. Such information is used to estimate the underground flow in the
reservoir. For instance, Abbaszadeh-Dehghani and Brigham [1] [2], and Oliver [13]
showed some successful flow-simulation and field examples. However, conducting
frequent tracer tests is often uneconomical and time consuming. Moreover, tracer
testing cannot provide a dynamic view of the system: an estimate of the flow
characteristics is obtained after the testing, but the model parameters cannot be
updated unless a new test (e.g., using a different tracer) is conducted.
Besides characterization, engineers sometimes need to have a more complete
understanding of the reservoir as well as the ability to forecast the future flows of
the fields in order to improve the waterflood management. To achieve this, com-
prehensive reservoir simulations are widely used to simulate the whole reservoir (or
region-of-interests in the reservoir) by integrating all information available for the
reservoir. With the reservoir simulation approach, it is possible to predict both the
reservoir and individual well performance and decide on a better strategy for man-
agement and decision-making. However, in this kind of approach, the integration
of numerous data types throughout the life of the reservoir is required, and this
is time-consuming and costly. Some techniques for managing reservoir uncertainty
have been proposed including improved sampling techniques (such as using Monte
2
Carlo simulation [43] [40]) or speeding up the model simulation (such as using the
coarse grid simulation [26] or a predictive model [20]). But these kinds of reservoir
simulations are still not suitable for many field applications because of the high
complexity. Furthermore, some of the inputs for building this model may not avail-
able or determinable. As a result, it is rather impractical for reservoir engineers to
use complex simulation models for daily monitoring of operations.
1.2 A New Research Trend
Injection and production rates are the most abundant data available in any water-
flood operation and they often correlate to each other in some complex manner in
the reservoir. Recently, a variety of methods have been used to express the rate
performance of a production well as a function of injection rates in the surrounding
injectors. In all these works, the reservoir is viewed as a system that converts an
input signal (injection rate at injector) into an output signal (production rate at
producer) and the goal is to analyze input and output signals to estimate some
characteristics of the reservoir. If this can be achieved, the estimated parameters
can be used to facilitate waterflood management and optimization.
As an example, Heffer et al. [22] used Spearman rank correlations to relate
injector-producer pairs and associated these relations with geomechanics. Panda
and Chopra [50] used artificial neural networks to determine the interactions be-
tween injection and production rates. Albertoni and Lake [3] estimated the con-
nectivity between wells based on a linear model using the multiple linear regression
(MLR) method. Gentil [15] explained the physical meaning of MLR constants by
presenting the connectivity weight as a function of transmissibility. This trend
turned to a predictive modeling approach in 2006, in which Yousef et al. [64] [65]
3
improved previous work by building a simple reservoir model, called the capacitance
model, to describe the relationship between injection and production wells. They
used a parameter to describe the effects of compressibility, in addition to trans-
missibility, between the injection-production interwell channels. Liang et al. [38]
used this model, accompanied with a fractional-flow model, to perform the oil rates
optimization procedure in waterfloods. Sayarpour et al. [55] used the model and
presented analytic solutions to three different reservoir-control volume scenarios.
This facilitates the capacitance model’s application for rapid assessment at dif-
ferent levels of a field study, from a single well, to a group of wells, and to an
entire field. Weber et al. [62] extended this work to large scale reservoirs and
Izgec and Kabir [25] made a comparison between the model and transient-pressure
approaches, and proved that it is valid for both before and after the water break-
through. Although this work has been successful, the possible use of alternative
predictive models, together with theoretical discussions from a system identification
point of view, have never been investigated and are proposed for the first time in
this thesis.
The above approaches all assumed the system (reservoir) is time-invariant during
some period of time. In 2007, Liu et al. [39] proposed an extended Kalman filter
approach to characterize the reservoir by assuming that the parameters in the model
will change over time (time-varying system), so that the Kalman filter is used to
track and estimate the model parameters. Independent of these works, Thiele
and Batycky [59] [60] estimated the connectivity between wells using streamline-
based workflow. Their approach requires building a complete stream flow reservoir
model for the region of interest, which is hard to achieve for most fields with daily
operations.
4
1.3 Research Contributions
In our work, we follow the predictive modeling approach by addressing the problem
of estimating reservoir properties in waterfloods from injection/production rates
and note that it can be seen as a system identification problem, where injection
rates are inputs and production rates are outputs. System identification refers to
the process of building mathematical models of systems based on measured data.
Before the parameters are identified, it is necessary to build a relatively simple
but reliable (depending on applications of need) model. We denote these models
as predictive models (because the main objective of these models is to predict the
reservoir behavior) in this thesis. The model can be built by first principles (e.g.,
conservation of momentum, mass, and energy), or empirically, or as a combination
of both approaches.
The state-of-the-art capacitance model (CM) [64] and its applications can be
seen as a pioneering predictive modeling approach. But, in spite of its success, CM
still has some limitations and many other possible predictive models still have not
been investigated. Our first contribution is to build three new predictive models:
the first one is the finite-impulse-response model, which leads to impulse responses
with more general shapes than those possible within the capacitance model. The
second one is the distributed capacitance model, which extends the concepts of
CM to more heterogeneous reservoir scenarios. The third one is the multivariate
autoregressive model. This model takes into consideration the producer-to-producer
relationships, which are not captured in all previous modeling approaches.
As a second contribution, we first propose two approaches for model evaluation:
(1) evaluating their prediction ability on a fresh data set; (2) comparing their
interpretation in terms of reservoir characteristics to the ground truth knowledge
5
for that reservoir. We verify all proposed models using both evaluation methods.
The so-called“grey-box”approach for model estimation procedure is also introduced
and fully investigated in this dissertation. Moreover, we define a practical metric
and compare the prediction performance of all models under various scenarios, and
make some suggestions to reservoir engineers about how to select among these
predictive modeling approaches.
The third contribution is that we provide a unified framework for predictive
modeling. We first show that all predictive models proposed up to now all belong
to a special case of a general linear dynamic model. Moreover, from this framework,
we show that the transfer function of this general linear model can be interpreted
in terms of some reservoir characteristics (by given its poles) in the interwell region
between each well pair. Finally, the relationship between different models can be
easily seen from this framework.
We demonstrate a totally new application of the multivariate autoregressive
model: forecasting the performance when some producer is shutting-in or some new
producer is brought into operation. Using this model, one only needs to constraint
producers to operate at a certain rate in order to estimate the impact of a well
shut-in, while all previous modeling approaches actually require to actually shut-in
the well for forecasting, thus potentially reducing overall production significantly.
This new application makes it much more practical to control producers in order to
predict performance under several possible “What if” scenarios, which is our fourth
contribution.
Our last contribution is to consider the injection rate designs. It can be easily
shown that the shapes of the inputs (injection rates) affect the model estimation
results noticeably. To the best of our knowledge, the improvements for reservoir
parameter estimation that can be achieved by carefully designing injection rates
6
have not been discussed in the literature. There are two different approaches:
deterministic and stochastic. For deterministic approach, we propose a new pro-
cedure based on the property of inverse-repeat signals and that can generate a set
of periodic signals with vanishing cross-correlation to each other. For stochastic
approach, we apply a procedure in system identification to reservoir applications
and demonstrate its performance on some predictive model.
The outline of this dissertation is as follows: In chapter 2, we develop three
new predictive models by that address the limitations of the previous proposed
models. In chapter 3, we discuss the Verification these new models with different
approaches, and making a comparative analysis of these models so that reservoir
engineers have guidelines to use them. We also investigate the grey-box approach
for modeling and summarizing its influence. In chapter 4, we provide a general
and unified framework for all predictive models, and clarifying the relationship
between different models. Moreover, we try to interpret the transfer function of
linear model as some reservoir characteristics in the field. In chapter 5, we develop
a novel application for multivariate autoregressive models. This new application
makes it much more practical to control producers in order to predict performance
under several possible “What if” scenarios, without actually shutting-in the wells.
In chapter 6, we investigate different ways of injection rate design to improve the
parameter estimation for the predictive models. We develop a new deterministic
procedure to generate a set of injection rates with vanishing cross-correlation to
each other, and apply a common stochastic procedure to field applications. Both
methods are evaluated by some predictive model.
7
Chapter 2
New Predictive Models
In mature waterfloods, fluids are produced by two main factors: either they corre-
spond to the primary production rates from the reservoir itself, or they are due to
the pressure gradients caused by fluid injections. As discussed in the introduction,
new techniques aim to capture the well interactions using statistical techniques
based on the fluctuations of injection and production data. The core procedure
relies on modeling for the reservoir, especially between injection and production
wells. These predictive models can be estimated using a common methodology:
(1) building the models based on physical principles and statistical techniques; (2)
estimating the parameters in the models mainly by rate information; (3) using the
trained models to characterize the reservoir or predict the future performance.
In this chapter, we propose and develop several new predictive models. We
first discuss the general linear time-invariant system, and revisit the state-of-the-art
modeling approach: capacitance model (CM). By addressing the limitations of CM,
we propose several new models. The first one is the finite impulse response (FIR)
model, which releases the constraints on the impulse response shape of CM. This
leads to improvement of the prediction performance for some scenarios. Second, we
8
generalize the concepts of CM to deal with more heterogeneous reservoir cases. We
call this approach the distributed capacitance model (DCM).
The above models are all based on using injection rates to predict the future
production rates. All these models do not consider changes in production in certain
wells may affect other wells. To address the producer-to-producer interactions, we
propose a new model: multivariate autoregressive with extra inputs (M-ARX).
Finally, we discuss a general approach for finding the model parameters. This
approach is based on minimizing the residues between actual and estimated produc-
tion rates, which involves some non-linear optimization procedures. The so-called
“linear-in-the-parameter” property is discussed which leads to a much less complex
optimization.
2.1 Linear Time-invariant System
From the system point of view, if we have M inputs and N outputs, and assuming
the system is linear and time-invariant (LTI), the outputs y
j
(t)with j =1, ..., N
can be expressed as
y
j
(t)= y
j
(t
0
)+
M
i=1
u
ij
(t) (2.1)
where u
ij
(t) is a filtered version of the input u
i
:
u
ij
(t)= u
i
(t)∗ h
ij
(t)=
∞
ξ=0
u
i
(t− ξ)h
ij
(ξ)dξ (2.2)
The h
ij
(t) is denoted as the impulse response between input i and output j.Basi-
cally (2.1) and (2.2) mean that the output can be expressed as a linear combination
of some filtered version of the inputs. If we assume that the reservoir in the region-
of-interest (ROI) is a LTI system, the injectors are inputs, and the producers are
9
the outputs, as we will see, the main difference between various linear models is
depending on how they describe the impulse response h
ij
(t), which accounts for
any shape of attenuation and delays between the injector i and producer j.Ifone
injects a water pulse into the reservoir, which leads to changes in production over
time, the impulse response h
ij
(t) describes these changes. As the waterflood project
begins, it denotes secondary oil recovery period. During this period, the primary
production curve stays relatively constant compared to the influence caused from
injection rate fluctuations, so the constant value y
j
(t
0
) can be seen as the primary
depletion effects.
2.2 Revisiting the Capacitance Model
In this section, we revisit the capacitance model, as a predictive modeling approach
from rate fluctuations [63] [64]. The CM is based on deriving a total mass balance
equation with compressibility. Take one injector-producer well pair in a drainage
volume as a representative element, the governed material balance differential equa-
tion is given by:
c
t
V
p
dp
dt
= u(t)− y(t), (2.3)
where c
t
is the total compressibility, V
p
is the drainage pore volume, P is the average
pressure in V
p
, u(t) is the injection rate (input) and y(t) is the total gross production
rate (output). Also, a linear productivity model can be used to relate pressure to
rates:
y = J(P − P
wf
) (2.4)
where P
wf
and J are the flowing bottom hole pressure (BHP) and productivity
index of the producer, respectively. Note that (2.4) holds only for stabilized flow,
10
Figure 2.1: The response of capacitance model by step injection. Different colors
show the responses of different time constant τ.
so its appropriateness can only be established by numerical simulation [64]. By
eliminating the average pressure on (2.3) and (2.4) using these two equations, and
solving the differential equation, the resulting production rates can be expressed as:
y(t)= y(t
0
)e
−(t−t
0
)
τ
+
e
−t/τ
τ
ξ=t
ξ=t
0
e
ξ/τ
i(ξ)dξ
+ J
P
wf
(t
0
)e
−(t−t
0
)
τ
− P
wf
(t)+
e
−t/τ
τ
ξ=t
ξ=t
0
e
ξ/τ
P
wf
(ξ)dξ
(2.5)
where the τ is defined as the “time constant” of the drainage volume:
τ =
c
t
V
p
J
(2.6)
and t
0
is the initial time.
(2.5) summarizes the capacitance model, and is easy to interpret: it shows that
the production rates of the producer can be divided into three components. The
11
() yt
wf
P
() ut
1
R
J
=
tp
CcV =
Z = ∞
Current
Source
Voltage
Source
Resistor
Capacitor
Figure 2.2: Analogous resistor-capacitor circuit for capacitance model.
first term accounts for primary production depletion. The second component is the
contribution from the injection rates. The impulse response shape is determined
by the variable τ, and Fig. 2.1 shows the response curve of different τs. The last
component is the influence of production rates caused by the change on bottom-hole
pressure of the producer.
We can make an analogy between the CM and the resistor-capacitor (RC) cir-
cuit, shown in Fig. 2.2. This analogy highlights the resistive and capacitance effects
of this model: the fluid flowing into the reservoir is like the electric current in the
circuits, and the reservoir behaves like a capacitance with C = c
t
V
p
and a re-
sistance with R =
1
J
, and thetimeconstant τ is equal to R× C. So generally
speaking, the CM tries to model the reservoir with its capacitance effects, together
with its resistive effects. For this reason, several researchers have called this model
“capacitance-resistive model” or CRM for short [54] [62].
This representative element can easily be extended to multiple injectors and
producers using superposition principle. For each injector-producer pair, weights
12
or coefficients, λ
ij
, are defined to incorporate the fact that one injector is shared by
more than one producer, and τ
ij
describes the impulse response shape between the
injector i and producer j. For the formulas and further reference, see [64] for more
discussions.
2.3 Finite Impulse Response Model
For CM with multiple wells, the λ
ij
accounts for the weight of this impulse response,
and τ
ij
accounts for the shape of the response. If we try to model the reservoir
by a general LTI system, the drawbacks of CM is it only uses one parameter τ
to describe the response shape. From a petroleum engineering point of view, it
only uses one parameter to characterize the attenuation and delays between each
injector-producer well pair. Even in numerical simulations, we can easily see that
one parameter is not enough for describing the shape of impulse response [35]. Thus,
we develop a new model by relexing the constraints on the impulse response and only
assume that the response is finite (which is reasonable for physical phenomenon).
This leads us to the finite-impulse-response (FIR) model.
This model can be built directly from the general LTI system description: sup-
pose we have M injectors and N producers in the ROI, the FIR model can be
written the same as (2.7) and (2.8), where now the length of impulse response
is limited to L (which means the response is finite). That is, the model can be
expressed as
y
j
(t)= y
j
(t
0
)+
M
i=1
u
ij
(t) (2.7)
with
u
ij
(t)= u
i
(t)∗ h
ij
(t)=
L
ξ=0
u
i
(t− ξ)h
ij
(ξ)dξ. (2.8)
13
The (2.7) and (2.8) represents the basic FIR model, where the production rate
of each well depends on its surrounding injectors and is captured by the finite
impulse response. The FIR model can be approximated to other predictive models,
provided that the response length of the FIR model L is large enough. But the
penalty of doing this is also obvious: it may requires an unnecessarily large number
of parameters. This will be discussed in more detail in Chapter 3.
2.3.1 Discrete Model
For any predictive model, history matching on performance data is a necessary
step for estimating parameters in the models. Injection and production rates are
measured by sampling, so discrete models are preferable in general. For FIR model,
we discretize the integrals in (2.8) leading to
y
j
(n)= y
j
(n
0
)+
M
i=1
u
ij
(n) (2.9)
and
u
ij
(n)= u
i
(n)∗ h
ij
(n)=
L−1
m=0
u
i
(n− m)h
ij
(m) (2.10)
where L is the length of h
ij
(n), and n is the discrete index, n ∈ Z.Asinthe
continuous case, 2.9 and 2.10 state that the total production rate at time n is a
function of the primary production component and the injection history between n
and n− L+1.
After constructing the FIR model, the next step is to estimate the unknown
parameters in the model. We will discuss a general approach for training the model
parameters in Section 2.6.
14
2.4 Distributed Capacitance Model
In (2.3) and (2.4), the pressure support that helps maintain production rates at pro-
ducers is denoted as the average pressure in the reservoir. The “average” pressure
accounts for the averaging effects which are caused by the whole interwell region
of the reservoir, including any heterogeneity and discontinuity in-between. If the
reservoir between injector-producer pair includes some heterogeneity, such as the
presence of some high permeability channels, fractures or faults, the analytic solu-
tion of CM [64] should be modified to take the heterogeneity into consideration. To
address this, we consider a new representative element, which is constructed by one
injector and one producer with a fracture or high permeability channel in-between.
This is illustrated in Fig. 2.3. The high permeability channel divides the interwell
region into two parts, and for each part, we can derive the material balance dif-
ferential equation with a linear injectivity (productivity) equation, similar to CM.
Mathematically, the equations in first region can be expressed as:
u(t)− u
f
(t)= c
t1
V
p1
dP
1
dt
(2.11)
where c
t1
is the compressibility of the first interwell region; V
p1
is the drainage
pore volume for this region; P
1
is the average pressure in V
p1
; u(t) is the injection
rate and u
f
(t) is the total liquid flowing from this region to the high permeability
channel. The linear injectivity equation for this region can be expressed as:
u
f
= J
1
P
1
(2.12)
15
where J
1
is the injectivity index of the injector. Similarly, we can derive equations
for the second region:
u
f
(t)− y(t)= c
t2
V
p2
dP
2
dt
(2.13)
and
y = J
2
P
2
− P
wf
(2.14)
where J
2
is the productivity index and P
wf
is the flowing bottom hole pressure
(BHP) of the producer. (2.11), together with (2.13), state that the net rates of
mass depletion from the drainage volume are proportional to the change of aver-
age pressure in the first and second regions, respectively. Combining these four
equations and eliminating the average pressure P
1
and P
2
,we get:
τ
1
τ
2
d
2
y
dt
2
+(τ
1
+ τ
2
)
dy
dt
+ τ
1
τ
2
J
2
d
2
P
wf
dt
2
+ τ
2
J
2
dP
wf
dt
= u(t)− y(t) (2.15)
where τ
1
and τ
2
are the “time constant” of the drainage volume similar to the
definition of CM. They are defined by:
τ
1
=
c
t1
V
p1
J
1
and τ
2
=
c
t2
V
p2
J
2
(2.16)
Although the
dP
wf
dt
term appeared explicitly in capacitance model, pressure data are
often not available on a daily basis in real fields. For most applications, researchers
assume that the changes in BHP are slow compared to the fluctuations in injection
and production rates. Thus we assume that the terms relating to
d
2
P
wf
dt
2
and
dP
wf
dt
are
close to zero and can be removed from the equation. Note that our new proposed
model still has the ability to incorporate the BHP data when they are available. It
16
Figure 2.3: Representative element for distributed capacitance model. There is a
high permeability channel between the injector and producer.
only needs to re-derive the equations to include these terms, but this is out of the
scope of most daily applications so we just skip the derivations here.
After removing the BHP terms, (2.15) is a linear second order differential equa-
tion with constant coefficients, and the solution is well-known:
y(t)=
⎧ ⎪⎪⎨ ⎪⎪⎩
c
1
e
−t
τ
1
+ c
2
e
−t
τ
2
+
1
τ
1
−τ
2
t
0
e
ξ−t
τ
1
u(ξ)dξ−
t
0
e
ξ−t
τ
2
u(ξ)dξ
if τ
1
= τ
2
c
1
e
−t
τ
1
+ c
2
te
−t
τ
1
+
1
τ
2
1
t
0
(t− ξ)e
ξ−t
τ
1
u(ξ)dξ if τ
1
= τ
2
(2.17)
The case τ
1
= τ
2
rarely happens in real geological situations. Thus here we discuss
the latter case without the loss of generality. In this case, the first term accounts for
the primary production associated with the total gross production; the second term
is the contribution from the injected water and describes the interactions between
injection and production rates. Note that (2.17) is the basic form of the model.
17
The impact from injected water can be shown more clearly in terms of convolu-
tion between injection rate and a particular shaped function; that is, we can rewrite
(2.17) as:
y(t)=
c
1
e
−t
τ
1
+ c
2
e
−t
τ
2
+ u(t)∗
1
τ
1
− τ
2
e
−t
τ
1
− e
−t
τ
2
(2.18)
and from a system point of view, the impulse response h(t) can be expressed as:
h(t)=
1
τ
1
− τ
2
e
−t
τ
1
− e
−t
τ
2
(2.19)
The influence of water injection in the injection history is described by the time
constants τ
1
and τ
2
. The time constants, similar to CM, can be seen a direct
measurement of dissipation of the pressure waves between injector-to-fracture and
fracture-to-producer. As an illustration, Fig. 2.4 shows h(t) with several different
τ
1
and τ
2
. If one of the values of τ is small, the model behaves just like the original
CM; otherwise, the new model leads to a “smoother” response than CM, as shown
in Fig. 2.4.
As in Section 2.3.1, in order to apply this model to real field applications, it
should be discretized with some selected discretization interval. The discrete version
of the model is given by:
y(n)=
c
1
e
−n
τ
1
+ c
2
e
−n
τ
2
+
1
τ
1
− τ
2
n
k=0
e
k−n
τ
1
u(k)−
n
k=0
e
k−n
τ
2
u(k)
(2.20)
(2.20) can also be written in terms of convolution between injection rate and impulse
response function:
y(n)=
c
1
e
−n
τ
1
+ c
2
e
−n
τ
2
+ u(n)∗ h(n) (2.21)
18
0 50 100 150 200
0
100
200
300
400
500
600
Time
Inj. and Pro. Rates
τ
1
= 1
Impulse Water Injection
τ
2
= 1
τ
2
= 5
τ
2
= 10
τ
2
= 20
0 50 100 150 200
0
100
200
300
400
500
600
Time
Inj. and Pro. Rates
τ
1
= 5
Impulse Water Injection
τ
2
= 1
τ
2
= 5
τ
2
= 10
τ
2
= 20
0 50 100 150 200
0
100
200
300
400
500
600
Time
Inj. and Pro. Rates
τ
1
= 10
Impulse Water Injection
τ
2
= 1
τ
2
= 5
τ
2
= 10
τ
2
= 20
0 50 100 150 200
0
100
200
300
400
500
600
Time
Inj. and Pro. Rates
τ
1
= 20
Impulse Water Injection
τ
2
= 1
τ
2
= 5
τ
2
= 10
τ
2
= 20
Figure 2.4: Impulse response curves of distributed capacitance model for different
τ
1
and τ
2
. The injected water of injector is also shown.
where the impulse response function h(n) is defined by
h(n)=
1
τ
1
− τ
2
e
−n
τ
1
− e
−n
τ
2
(2.22)
As in the continuous case, (2.21) states the total production rate at time step n can
be divided into two components. The first term accounts for the primary depletion
of the production rates, and the second term is the contribution of injected water.
In general, there are many injectors and producers in a reservoir. The production
rates at one producer are often supported by several surrounding injectors. So we
need to generalize our proposed model to the multiple producers and injectors case.
Given M injection wells and N production wells in the ROI, we use coefficients
or weights λ
ij
with i =1, 2, ..., M and j =1, 2, ..., N in order to capture the fact
19
that one producer is supported by more than one injector (similar to CM). Putting
these weights into the expression, the generalized model between the producer j
and the M surrounding injectors can be expressed as
y
j
(n)= c
1
e
−n
τp
+
M
i=1
λ
ij
i
i
(n)∗ h
ij
(n) (2.23)
The first term accounts for the primary production of this producer, and represented
by an exponential decay function of “total effect” time constant τ
p
. The impulse
response function for this well pair, h
ij
(n), is defined as
h
ij
(n)=
1
τ
ij1
− τ
ij2
e
−n
τ
ij1
− e
−n
τ
ij2
(2.24)
τ
ij1
and τ
ij2
are two parameters which represent the dissipation of pressure wave
between injector i and producer j when this well pair is the only active well pair
in the reservoir. (2.23), together with (2.24), represent the distributed capacitance
model for the multiple wells case.
For the estimation of model parameters, see Section 2.6.
2.4.1 Interpretation of Model
This model can be seen as a cascading of two CMs, that is, where the output of first
CM becomes the input of the second CM, and these two CMs are combined together
to form a bigger system, as shown in Fig. 2.5. This is because we can see the each
interwell region in the considered scenario (see Fig. 2.3) as a relative homogeneous
region, and separated by some discontinuous region (e.g., high permeability channel
in Fig. 2.3). It is this interpretation that led us to call this model as Distributed
Capacitance Model (DCM).
20
Figure 2.5: Interpretation of distributed capacitance model in terms of CMs.
Thus, DCM can also be interpreted in terms of some particular geological fea-
tures. For example, if the reservoir has many layers and the injection wells and
production wells are in different layers, they must communicate with each other
by passing some interface between different layers. These layers may be relatively
homogeneous, so the total communication path can be seen as the pressure pass-
ing two different regions, which means the path is heterogeneous even though each
region is homogeneous. In these kinds of situations, the proposed DCM is more
suitable than CM because it accounts for some heterogeneity between the interwell
regions. Moreover, DCM will behave similarly to CM when one of the parameters
τ is small, so DCM can be seen as a generalized version of CM.
21
2.5 Multivariate ARX Model
For all proposed models up to now, the production rates are described as a linear
combination of some filtered version of injection rates from surrounding injection
wells. Generally speaking, in this kind of approaches, the influences from other pro-
ducers are approximated by the injection rates of injection wells. This means that
the influence a producer on neighboring producers is ignored and the production
rates are modeled only as a function of injectors in the reservoir. To achieve a better
estimation, we develop a multivariate autoregressive with extra inputs (M-ARX)
model in [33] to take the producer-to-producer interactions into account.
Most predictive models first consider a 1-injector/1-producer scenario, and then
extend the derived relationship to multiple injectors and producers using the super-
position principle. In this kind of approach, every producer is treated independently.
To account for this producer-to-producer influence, we will follow a derivation sim-
ilar to that leading to the CM in [64], but considering the whole reservoir including
multiple injectors/producers.
We consider a ROI in the reservoir with M injectors and N producers, as shown
in Fig. 2.6. For the region close to a given producer j, the material balance equations
are given by a set of differential equations:
⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ c
t
V
p1
d
¯
P
1
dt
=
M
i=1
α
1i
u
i
(t)−
N
j=1
β
1j
y
j
(t)
c
t
V
p2
d
¯
P
2
dt
=
M
i=1
α
2i
u
i
(t)−
N
j=1
β
2j
y
j
(t)
.
.
.
c
t
V
pN
d
¯
P
N
dt
=
M
i=1
α
Ni
u
i
(t)−
N
j=1
β
Nj
y
j
(t)
(2.25)
22
1
u
3
u
3
y
2
u
2
y
1
y
Figure 2.6: Scenario with multiple injection and production wells. The red shadow
areas represent the pore volume related to that producer.
where V
pj
with j =1, ..., N is the drainage pore volume related to producer j;
¯
P
j
is the average pressure in V
pj
. α
ki
and β
kj
with k =1, ..., N are the weight
factors for injection rates u
i
(t) and production rates y
j
(t), respectively. This set
of equations is similar to the material balance equation of CM as shown in (2.3).
Note that for ROIs with closed boundaries, we will also have that
M
k=1
α
ki
=1and
N
k=1
β
kj
= 1, but for most practical use scenarios for the model, these conditions
often no longer hold because of the open boundary around the ROI. We also use a
linear productivity model to change the pressure variables in the set of equations
into rates. This set of linear productivity model equations can be expressed as
⎧ ⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎩ y
1
(t)= J
1
(
¯
P
1
− P
1,wf
)
y
2
(t)= J
2
(
¯
P
2
− P
2,wf
)
.
.
.
y
N
(t)= J
N
(
¯
P
N
− P
N,wf
)
(2.26)
23
where P
j,wf
with j =1, ..., N is the flowing bottom-hole-pressure (BHP) of producer
j and J
j
is the productivity index of that producer. Substituting for the average
pressure in (2.25) using (2.26), we get
⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎩ τ
1
dy
1
(t)
dt
+ τ
1
J
1
dP
1,wf
(t)
dt
=
M
i=1
α
1i
u
i
(t)−
N
j=1
β
1j
y
j
(t)
.
.
.
τ
N
dy
N
(t)
dt
+ τ
N
J
N
dP
N,wf
(t)
dt
=
M
i=1
α
Ni
u
i
(t)−
N
j=1
β
Nj
y
j
(t)
(2.27)
where τ
j
with j =1, ..., N is the“time constant”of the drainage volume of producer
j andisdefinedby
τ
j
=
c
t
V
pj
J
j
(2.28)
In (2.27), the BHP terms
dP
j,wf
(t)
dt
are often set to zero because we assume
that the BHPs are changing slowly, as compared to the changes in injection and
production rates, as discussed in Section 2.4. This assumption is important in
practical uses of the model because for many fields the BHP information is often
not available. Therefore most injector-producer modeling techniques make this
assumption. For techniques that can be used in situations where large variations
of BHP occur, refer to [30].
Here we assume that changes in BHPs are negligible, and divide each of the
equalities in (2.27) by its corresponding τ
j
, so that (2.27) becomes
⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎩ dy
1
(t)
dt
=
M
i=1
b
1i
u
i
(t)−
N
j=1
a
1j
y
j
(t)
.
.
.
dy
N
(t)
dt
=
M
i=1
b
Ni
u
i
(t)−
N
j=1
a
Nj
y
j
(t)
(2.29)
24
We define a set of variables b
ki
= α
ki
/τ
j
and a
kj
= β
kj
/τ
j
. The equalities in (2.29)
are the 1
st
-order multivariate autoregressive model with some extra inputs (1
st
-
order M-ARX model). The order of the model refers to the maximum order of
the differential operator (a first order differential operator is used in (2.29)) This
modeling approach is well-known in time series analysis and system identification
research [28] [52]. For convenience, the M-ARX model is often expressed in ma-
trix form. Let y and u represent the production rate and injection rate vectors,
respectively; that is,
y(t)=[y
1
(t) y
2
(t) ... y
N
(t)]
T
(2.30)
and
u(t)=[u
1
(t) u
2
(t) ... u
M
(t)]
T
(2.31)
We also define two coefficient matrices A and B as
A
c
=
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ a
11
··· a
1N
.
.
.
.
.
.
.
.
.
a
N1
··· a
NN
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (2.32)
B
c
=
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ b
11
··· b
1M
.
.
.
.
.
.
.
.
.
b
N1
··· b
NM
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (2.33)
where A
c
is N by N and B
c
is N by M and the subscript c denotes the continuous
model. Thus, (2.29) can be expressed in matrix form as
dy(t)
dt
+ A
c
y(t)= B
c
u(t) (2.34)
25
where A
c
represents the “characteristic” of the autoregressive behavior of the sys-
tem, while B
c
represents the weighting coefficients from the system inputs u.
2.5.1 Discrete Model
Discretizing (2.29), we can get the following expression
y(k+1)=−Ay(k)+ Bu(k) (2.35)
The matrices A and B are both coefficient matrices and correspond to the autore-
gressive behavior and the input weights of the system, respectively. (2.35) shows
that the production rate of any producer at time instant t = k+1, can be computed
in terms of two factors: one that depends on ALL production rates at time instant
t = k; the other on the injection rates at t = k. This means that the system out-
puts (production rates) are described by both the new inputs (injection rates) and
a feedback loop from the outputs of the previous time instant, as shown in Fig. 2.7.
After the M-ARX model is built, the next step is to estimate the unknown
coefficients in the model. For the model training process, please refer to Section
2.6.
2.6 Finding the Model Parameters
To estimate the unknown model parameters, the most common approach is to try to
“fit”the observed data using a quadratic function of the fitting-errors [41]. Suppose
26
Reservoir
System
Delay Operator
M-ARX
Model
Injection Rates
Production Rates
() k u
() k y
(1) k= u
(1) k= y
ˆ() k y
Estimated
Production Rates
Figure 2.7: Diagram of the modeling process of M-ARX model.
we have a model M with some parameter vector θ, then we try to minimize the
following criterion:
V
K
(θ)=
1
2K
K
k=1
[y(k)− ˆ y(k|θ)]
T
[y(k)− ˆ y(k|θ)] (2.36)
where ˆ y(k|θ) denotes the vector of predicted production rates at time k under
model M with parameter value θ. For example, for CM, θ is
θ =[λ
11
... λ
MN
τ
11
... τ
MN
]
T
(2.37)
For each particular value ofθ
∗
, we have a corresponding V
K
(θ
∗
), so the optimization
procedure is basically searching the minimum value over a space with dimension
M×N×2. This can be achieved by some iterative search methods, such as gradient
or steepest-descent method. For further reference, see [41].
27
The above procedure discusses the general approach for estimating model pa-
rameters. In our work, this method still involves many computations and we only
follow some simpler procedure from [63]. That is, we first start from τ
ij
=0.1to
30 with increasing 0.1 each time. For each fixed τ
ij
, the remaining parameters in
criterion V
K
(θ)are λ
ij
s, which can be solved by multiple linear regression. After
iterating among all possible τ
ij
, we choose the parameter set which has the mini-
mum value of V
K
as the optimization results. For DCM, We also follow the same
procedure, with iterating among both τ
ij,1
and τ
ij,2
now.
The optimization procedure of (2.36) is a nonlinear optimization problem. To
solve these kinds of problems is often computationally expensive and no global
minimum solutions can be guaranteed. On the other hand, both the FIR and M-
ARX models have the property of being “linear-in-the-parameter”, which means
that the parameters are linear on the criterion function (for a discussion of the
linear-in-the-parameters property of models, see [23]). This allows the use of linear
regression techniques for parameter estimation. To see this, we take the M-ARX
model as an example (the procedure for FIR model is similar). We first define a
parameter matrix as follows:
θ =[AB]
T
(2.38)
and the vector of regressors (regression vector) as
ϕ(k)=
⎡ ⎢ ⎣ −y(k− 1)
u(k− 1)
⎤ ⎥ ⎦ (2.39)
So that the predicted production rate values at time k become
ˆ y(k|θ)=θ
T
ϕ(k) (2.40)
28
and the minimization of the criterion
V
K
(θ)=
1
2K
K
k=1
y(k)−θ
T
ϕ(k)
T
y(k)−θ
T
ϕ(k)
(2.41)
can be seen as N different linear regressions, corresponding to each row of θ
T
,all
with the same regression vector ϕ(k). Thus minimizing (2.41) will be equivalent to
solving the N linear regressions.
2.7 Summary
In this chapter, we have developed three different predictive models. The first one
was the FIR model, which allows very flexible modeling of the shape of impulse
response using a set of unknown parameters. The second model, DCM, was devel-
oped as a generalization of CM, which can deal with more heterogeneous scenarios.
Then we also proposed the third model, M-ARX model, to address the interactions
between producers to producers to facilitate the model building, which is deficient
in consideration for previous models. Finally, we discuss a general approach for
estimating the model parameters, and the linear-in-the-parameters property was
also pointed out.
29
Chapter 3
Model Validation and Comparison
In this chapter, we first discuss the model validation issue. The central question is:
after building the model, how can we determine whether the model is suitable for
field applications or not? Do we have any validation step to double-check the model
we just built? To answer these questions, we will discuss two different approaches
to validate a predictive model in the reservoir: the first is trying to validate models
based on their prediction ability; another is trying to validate models by interpret-
ing the model parameters and comparing them to known geological features (or
characteristics estimated from other approaches).
The information we have used for building models up to now is often purely sta-
tistical, with only limited physical knowledge involved. In order to build a stronger
model, we also investigate the use of “grey-box” modeling approach on predictive
models. The grey-box approach refers to integrating some physical knowledge into
optimization procedures. Finally, we provide a comparative analysis by proposing
a prediction-error metric under various scenarios. Some noise sensitivity analysis is
also performed to deal with situations encountered in real fields. From the simu-
lation results, together with the theoretical discussions, we make some suggestions
for reservoir engineers to use the predictive models on field applications.
30
3.1 Validation of Models
To test whether a predictive model is suitable or not, there are two common ap-
proaches: the first approach is to evaluate its prediction ability. This can be
achieved by validation on a given data set. The second approach is to validate
by the interpretation of parameters. This can be achieved by comparing the es-
timated parameters against the theoretical values of the synthetic data or known
geological features of the real fields. We will discuss both validation procedures in
the section.
To examine whether a built model is suitable for a physical system, the most
natural approach is to evaluate model behavior on a fresh data set that has not
been used for training [57]. Thus, given a historic data set, we can separate it into
two data sets: the estimation (or training) data set and the validation (or testing)
data set. We first estimate the model parameters on the estimation data set, then
compute and measure the error on the validation data set. In addition, we can
use this validation procedure to evaluate different predictive models and favor the
model that shows the better prediction performance on a fresh data set. This leads
to the quantitative comparisons between models, which will be discussed later in
this chapter.
A second way to evaluate a model is to first estimate the model parameters,
which are often related to some characteristics of inter-well region, such as per-
meability or porosity, and then compare the estimated values against the ground
truth or some known geological features (of real field data). Although sometimes
only qualitative comparisons are available, this kind of evaluation highlights one
very important practical usage of model: it can be used to indicate and estimate
some important reservoir characteristics of the reservoir. In our problem, the most
31
important reservoir characteristics are the static gains between each input-output
pair in the system. The static gains between well pairs are interpreted as the ef-
fective contribution from a particular injector to a specific producer, and are given
different names by different authors, e.g., interwell connectivity in the context of
CM [64] [54]; injector-producer relationships for Kalman filter approaches [39] [67];
effective flow units for active method and FIR model [34] [35]. The reader should
keep in mind that although these are different names, they all refer to the static
gain of the system, and we will use these terms interchangeably in our work.
The static gains are easy to calculate for different predictive models. For CM
and DCM, the static gain between the i-th injection well (input i)and the j-th
production well (output j) is denoted as λ
ij
and appears in the models explicitly.
For FIR model, λ
ij
can be estimated as:
λ
ij
=
L−1
n=0
h
ij
(n). (3.1)
For the 1
st
order M-ARX model, the static gain λ
ij
can be estimated easily from
z-domain. To see this, we apply the z-transform in (2.29) and get:
zY(z)=−AY(z)+ BU(z) (3.2)
where Y(z)and U(z) are the z-transforms of y(k)and u(k), respectively. Then we
can write the relationship between Y(z)and U(z)as:
Y(z)=(zI + A)
−1
BU(z)
=
I + Az
−1
−1
Bz
−1
U(z), (3.3)
32
where I denotes the N by N identity matrix. (3.3) uses the matrix fraction descrip-
tion (MFD) for the transfer function matrix between Y(z)and U(z). For detailed
discussions about MFD, see [28]. The frequency domain representation can be ob-
tained by setting z = e
jω
where ω denotes the angular frequency. By definition, the
static gain (interwell connectivity) represents the cumulative effect of total injection
on total production and thus can be computed as the transfer function at ω =0,
i.e., at the zero frequency. In z-domain, the zero frequency corresponds to z =1
and if we substitute this into (3.3) we obtain
Y(z)
U(z)
z=1
=(I + A)
−1
B, (3.4)
which can be interpreted as follows: once the M-ARX model parameter matrices
A and B have been computed, the interwell connectivities are represented by the
N by M matrix (I + A)
−1
B,inwhich the element on j-th row and i-th column
denotes the interwell connectivity between injector i and producer j. This shows
the M-ARX model can be used to estimate the interwell connectivity, as can be
done for other predictive models.
3.2 Simulation Setting
In this section, we will validate and compare different models using using a numer-
ical flow-simulator, named CMG [42], under various scenarios. To make it more
clear, we summarize in Table 3.1.
For the pattern of injection and production wells, we use three different patterns:
(1) five-spot with five injectors and four producers (Scenarios A and B); (2) line-
drive with six injectors and three producers (Scenarios C and D); (3) peripheral
33
Table 3.1: Simulation Settings for CMG simulator.
Scenario Inj. & Pro. Pattern Permeability Map Inj. Rate
A five-spot channel (Fig. 3.1)
From [63] (Fig. 3.6)
B (5-inj./4-pro.) streak (Fig. 3.2)
C line-drive homogeneous (Fig. 3.3)
D (5-inj./4-pro.) multiple fractures (Fig. 3.4)
E peripheral injector hydraulic fractures one injector
F (8-inj./6-pro.) (Fig. 3.5) step (Fig. 3.7)
G variable (Fig. 3.8)
Figure 3.1: Permeability map for channel case (Scenario A). Note that Producer 4
is included in the high permeability area.
injectors with 8 injectors and 6 producers (Scenarios E, F, and G). It constructs
with all producers close to the center and surrounded by injectors.
34
Figure 3.2: Permeability map for streak case (Scenario B). Note that well pairs
I1-P1and I3-P4 are connected by high permeability channels.
100 200 300 400 500 600 700
50
100
150
200
250
300
350
400
450
500
550
Figure 3.3: Permeability map for homogeneous case (Scenario C).
3.2.1 Permeability Map
For the five-spot pattern, we test two different permeability maps: (1) channel case,
where there is a large area passing through Producer 4 with permeability of 500 md
35
100 200 300 400 500 600 700
50
100
150
200
250
300
350
400
450
500
550
Figure 3.4: Permeability map for multiple fractures case (Scenario D).
100 200 300 400 500 600 700 800
100
200
300
400
500
600
Figure 3.5: Permeability map for Scenarios E, F, and G.
and the background permeability is of 10 md; (2) streak case: there are two high
permeability streaks passing through I1to P1and I3to P4, and they are with
permeability of 1000 md and 500 md, respectively (Fig. 3.2).
36
0 500 1000 1500 2000 2500 3000 3500
0
1000
2000
3000
4000
Time(day)
Injection Rate
Injection Rates on I1
0 500 1000 1500 2000 2500 3000 3500
0
500
1000
1500
2000
Time(day)
Injection Rate
Injection Rates on I2
0 500 1000 1500 2000 2500 3000 3500
0
500
1000
1500
Time(day)
Injection Rate
Injection Rates on I3
0 500 1000 1500 2000 2500 3000 3500
0
500
1000
1500
Time(day)
Injection Rate
Injection Rates on I4
0 500 1000 1500 2000 2500 3000 3500
0
1000
2000
3000
Time(day)
Injection Rate
Injection Rates on I5
Figure 3.6: Injection rates for Scenarios A-D. For Scenarios C and D, we randomly
generate the injection rates for Producer 6.
For the line-drive pattern, it also has two different permeability settings: (1)
homogeneous case: there is an isotropic permeability of 100 md for all directions
(Fig. 3.3); (2) multiple fractures case: there are three fractures (high permeabil-
ity channels) of permeability 1000 md passing through the production wells with
different lengths and all lay in about 45 degree direction (Fig. 3.4).
For the peripheral injectors pattern (Scenarios E, F, and G), the reservoir has
a background permeability of 1 md, and there are 14 fractures of 1000 md passing
through all injection and production wells, with a 135 degree orientation (3.5). The
fractures passing through wells are analogous to the hydraulic fractures.
37
Figure 3.7: Injection rates set for Scenario F. At one time, only one injector changes
the injection rate from 500 to 1000 bbl/day.
3.2.2 Injection Rate Setting
For five-spot and line-drive patterns (Scenarios A-D), the injection rates we used
are the same as those used in [63], which were obtained from a real oilfield. These
rates are shown in Fig. 3.6.
Injection rates is the only difference between Scenario E, F, and G. For Scenario
E, we only inject water into one injector at a time and observe the response from
other wells, so we denote this scenario as one injector case. For Scenario F, we
have a step function in each injector, as shown in Fig. 3.7. We denote this scenario
as step case. For Scenario G, each injection rate is set to constant for each 200
days, and all injection rates are made to vary randomly between a predetermined
set of values, as shown in Fig. 3.8. This scenario is denoted as variable case. We
38
500 1000 1500 2000
0
500
1000
Time(day)
Injection Rate
I − 1
500 1000 1500 2000
0
500
1000
Time(day)
Injection Rate
I − 2
500 1000 1500 2000
0
500
1000
Time(day)
Injection Rate
I − 3
500 1000 1500 2000
0
500
1000
Time(day)
Injection Rate
I − 4
500 1000 1500 2000
0
500
1000
Time(day)
Injection Rate
I − 5
500 1000 1500 2000
0
500
1000
Time(day)
Injection Rate
I − 6
500 1000 1500 2000
0
500
1000
Time(day)
Injection Rate
I − 7
500 1000 1500 2000
0
500
1000
Time(day)
Injection Rate
I − 8
Figure 3.8: Injection rates set for Scenario G.
will discuss why we use three different kinds of injection rates for the peripheral
injector pattern in Section 3.3.
3.2.3 Other Setting
In all cases we simulate two component water and oil fluid systems, and have only
vertical wells, with the sampling interval set to 1 day, that is, assuming that daily
data are available. Five-spot and line-drive patterns (Scenarios A-D) simulate a
five-layered reservoir, while for the peripheral injectors pattern (Scenarios E, F,
and G), we consider one-layered reservoir simulation.
39
Table 3.2: Estimated values of interwell connectivity for Scenario A.
CM FIR Model
Pro 1 Pro 2 Pro 3 Pro 4 Pro 1 Pro 2 Pro 3 Pro 4
Inj 1 0.035 0.060 0.020 0.892 0.035 0.060 0.020 0.890
Inj 2 0.182 0.029 0.222 0.567 0.185 0.029 0.226 0.567
Inj 3 0.052 0.045 0.053 0.847 0.048 0.045 0.049 0.847
Inj 4 0.037 0.031 0.018 0.901 0.039 0.031 0.020 0.906
Inj 5 0.042 0.017 0.130 0.809 0.042 0.017 0.129 0.809
DCM M-ARX Model
Pro 1 Pro 2 Pro 3 Pro 4 Pro 1 Pro 2 Pro 3 Pro 4
Inj 1 0.035 0.060 0.020 0.893 0.023 0.063 0.015 0.905
Inj 2 0.182 0.029 0.222 0.566 0.181 0.030 0.222 0.566
Inj 3 0.052 0.045 0.053 0.847 0.083 0.038 0.066 0.812
Inj 4 0.038 0.031 0.018 0.898 0.069 0.025 0.033 0.870
Inj 5 0.042 0.017 0.130 0.811 0.024 0.020 0.121 0.830
3.3 Validation Results for Interwell Connectivity
In order to validate proposed models, we use the models we built to estimate the
static gain (interwell connectivity) of each well pair for Scenarios A and B. The
results are validated via comparisons with the state-of-the-art approach, CM, and
also with the ground truth of the simulators.
For the Scenario A, the total data available are approximately three thousand
days. This three thousand days period includes both before and after water break-
through time. We used 75th to 2074th day, totally 2000 days as the training period,
to estimate the static gains (interwell connectivity) of the system. The results are
shown in Table 3.2 and Fig. 3.9.
From the results, the estimated interwell connectivities from FIR model, DCM
and 1
st
order M-ARX model are all highly consistent with the values estimated
from CM. To be more specific, the average absolute differences between CM and
40
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
I−1 I−2
I−3
I−4 I−5
Estimated by CM
P−1
P−2 P−3
P−4
X
Y
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
I−1 I−2
I−3
I−4 I−5
Estimated by FIR Model (L=30)
P−1
P−2 P−3
P−4
X
Y
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
I−1 I−2
I−3
I−4 I−5
WAF Estimated by DCM
P−1
P−2 P−3
P−4
X
Y
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
I−1 I−2
I−3
I−4 I−5
Estimated by 1st Order M−ARX
P−1
P−2 P−3
P−4
X
Y
Figure 3.9: Estimated interwell connectivities with different predictive models for
Scenario A. Note that the longer arrows, represent the larger values.
Table 3.3: Average absolute differences of interwell connectivities between CM and
other models for Scenario A.
Model Avg. abs. diff. to CM
FIR model 0.0015
DCM 0.0006
M-ARX model 0.0131
other models are shown in Table 3.3. This shows that all predictive models can
generate highly consistent results when they are used to estimate the interwell
connectivities. If we compare the estimated results against the ground truth of
Scenario A (Fig. 3.1), the results are also consistent with qualitative descriptions
that can be inferred from the ground truth. For example, because P4is inside
the high permeability area, we expect that for all injectors the largest interwell
connectivities values will correspond to P4, which is the case for the estimated
values from all models. This further validates the results we got from all proposed
41
predictive models. Of course, the estimated interwell connectivities could be wrong
even if all the results are consistent, so we will provide quantitative verification in
the next section.
For the Scenario B, we used the 75th to 2074th days as the training period to
estimate the interwell connectivity. The results are shown in Table 3.4 and Fig. 3.10.
Table 3.4: Estimated values of interwell connectivity for Scenario B.
CM FIR Model
Pro 1 Pro 2 Pro 3 Pro 4 Pro 1 Pro 2 Pro 3 Pro 4
Inj 1 0.962 0.012 0.009 0.032 0.961 0.012 0.008 0.029
Inj 2 0.519 0.009 0.179 0.287 0.532 0.009 0.182 0.289
Inj 3 0.101 0.005 0.008 0.899 0.084 0.006 0.003 0.890
Inj 4 0.153 0.167 0.030 0.640 0.161 0.169 0.029 0.640
Inj 5 0.140 0.018 0.170 0.655 0.135 0.018 0.173 0.671
DCM M-ARX Model
Pro 1 Pro 2 Pro 3 Pro 4 Pro 1 Pro 2 Pro 3 Pro 4
Inj 1 0.962 0.012 0.009 0.032 0.966 0.012 0.009 0.030
Inj 2 0.519 0.009 0.179 0.287 0.521 0.009 0.180 0.287
Inj 3 0.101 0.005 0.008 0.898 0.091 0.006 0.006 0.899
Inj 4 0.153 0.166 0.029 0.640 0.149 0.169 0.030 0.646
Inj 5 0.140 0.019 0.171 0.655 0.141 0.017 0.170 0.655
Similar to the Scenario A (channel case), the estimated interwell connectivities
from all models are highly consistent with the values estimated from CM. The av-
erage absolute differences between CM and other models are shown in Table 3.5.
If we compare the results with the ground truth of Scenario B, the estimated val-
ues are also consistent with the geological features of the permeability distribution
(Fig. 3.2). To see this, we examine the estimated interwell connectivities between
I1-P1and I3-P4, which are very large values (close to 1), which indicates high
permeability channels exist between these well pairs. In short, both Scenario A
42
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
I−1 I−2
I−3
I−4 I−5
Estimated by CM
P−1
P−2 P−3
P−4
X
Y
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
I−1 I−2
I−3
I−4 I−5
Estimated by FIR Model (L=30)
P−1
P−2 P−3
P−4
X
Y
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
I−1 I−2
I−3
I−4 I−5
WAF Estimated by DCM
P−1
P−2 P−3
P−4
X
Y
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
I−1 I−2
I−3
I−4 I−5
Estimated by 1st Order M−ARX
P−1
P−2 P−3
P−4
X
Y
Figure 3.10: Estimated interwell connectivities with different predictive models for
Scenario B. Note that the longer arrows, represent the larger values.
Table 3.5: Average absolute differences of interwell connectivities between CM and
other models for Scenario B.
Model Avg. abs. diff. to CM
FIR model 0.0050
DCM 0.0001
M-ARX model 0.0021
and B validate that all proposed predictive models can generate estimated inter-
well connectivities similar to CM, and that these estimates are consistent with the
geological features of the ground truth in the simulations.
3.3.1 Quantitative Verification
To quantitatively verify the estimated interwell connectivities, we consider scenarios
E, F, and G, which have different kinds of injection rates. In Scenario E, we shut-in
all injectors except for one target injector each time, injecting some constant water
43
into this target injector, then observing the total fluid obtained from each producers.
The ratio between the produced fluid and the injected water can be calculated and is
a good approximation to the interwell connectivity for the corresponding well pair.
This is because the definition of interwell connectivity is the effective contribution
to a producer from a particular injector.
For Scenario F, we keep constant on all injection rates, and only change one
injection rate to another constant at each time. The injection rates are shown in
Fig. 3.7. As this target injection rate changes to another constant value, we would
expect all production rates to increase by different amounts. The ratio between the
increased fluid of this producer and the injected amount of water in the injector is
also a good approximation to the interwell connectivity for this well pair.
Scenario G is designed for validation of proposed predictive models. The simu-
lation time is 2160 days and we use the data obtained from Scenario G to estimate
the interwell connectivities. The estimated values are compared with the value ob-
tained (directly from the simulation data and divisions on the data) from Scenario E
and F. Here we only show the estimated interwell connectivities for FIR model, and
all other models have similar estimations. The results are shown in Table 3.6 and
also plotted in Fig. 3.11. From the results, we found the estimated interwell con-
nectivities are highly consistent with the simulation results, which strongly justifies
the use of predictive models for the estimation of interwell connectivity.
3.4 Physical Constraints on Parameters
When the statistical identification techniques are applied assuming a “black box”
model structure (i.e., no insights about the physical system are used), sometimes
44
Table 3.6: Estimated interwell connectivities from both simulations and FIR model
for peripheral injector pattern (Scenario E, F, and G).
Scenario Pro 1 Pro 2 Pro 3 Pro 4 Pro 5 Pro 6
Inj 1
E 0.73 0.14 0 0.025 0 0
F 0.73 0.15 0.05 0.05 0.01 0.00
G 0.74 0.15 0.05 0.05 0.01 0.00
Inj 2
E 0.18 0.48 0.27 0. 0 0
F 0.19 0.48 0.28 0.01 0.01 0.00
G 0.19 0.48 0.28 0.00 0.01 0.00
Inj 3
E 0 0.13 0.75 0 0 0
F 0.04 0.14 0.75 0.00 0.01 0.05
G 0.04 0.14 0.75 0.00 0.01 0.05
Inj 4
E 0.37 0.01 0 0.48 0.04 0
F 0.38 0.05 0.01 0.49 0.06 0.01
G 0.38 0.04 0.01 0.50 0.06 0.01
Inj 5
E 0 0.04 0.47 0 0.01 0.38
F 0.01 0.06 0.48 0.01 0.05 0.39
G 0.01 0.06 0.49 0.01 0.05 0.39
Inj 6
E 0.01 0 0 0.75 0.13 0
F 0.05 0.01 0.00 0.74 0.14 0.04
G 0.05 0.01 0.00 0.75 0.15 0.04
Inj 7
E 0 0 0 0.26 0.48 0.18
F 0.01 0.01 0.01 0.28 0.48 0.20
G 0.02 0.01 0.01 0.29 0.48 0.20
Inj 8
E 0 0 0.02 0.01 0.14 0.73
F 0.00 0.01 0.04 0.03 0.13 0.70
G 0.01 0.01 0.05 0.04 0.16 0.74
they produce unrealistic results which are not compatible with the underlying phys-
ical reality. For example, the use of CM sometimes results in negative interwell
connectivity values [64], which contradicts what we expect to be the reservoir be-
havior. While some models can relate some or all their parameters to the reservoir
characteristics, these parameters can also be constrained by some general physical
knowledge. In general, this identification process is often denoted as “grey-box”
modeling [61]. In our problem, an important constraint is the static gain between
45
(a) Calculated by Scenario E.
(b) Calculated by Scenario F.
(c) Estimated by FIR model for Scenario G.
Figure 3.11: Estimated interwell connectivities for peripheral injector pattern (Sce-
narios E, F, and G). Note that the longer arrows, represent the larger values.
46
each input-output pair in the system, whose value should be between zero and one.
As mentioned before, the static gains between well pairs can be interpreted as the
effective contribution from a particular injector to a producer, and should have a
non-negative impact that can be no greater than one. Theoretically speaking, the
constraints on static gains can be applied to all interwell models: the static gains
are the zero frequency behavior (z =1in z-domain) of the system. In other words,
the constraints can be written as
0≤ H(z)|
z=1
≤ 1 (element-wise) (3.5)
where H(z)is the z-transform of the matrix of impulse response h(n). In some
models, such as CM, DCM and FIR, the constraints on the system are equivalent
to constraints on linear combination of parameters. For example, in FIR model,
the constraints become:
0≤
⎡ ⎢ ⎢ ⎢ ⎢ ⎣
L
k=0
h
11
(k) ···
L
k=0
h
M1
(k)
.
.
.
.
.
.
.
.
.
L
k=0
h
1N
(k) ···
L
k=0
h
MN
(k)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ≤ 1 (element-wise) (3.6)
Furthermore, since fluid production cannot be negative, the impulse response can-
not be negative at any time, and their sum should also smaller or equal to 1. So
that in addition to (3.6), we also have constraints for each FIR coefficient:
0≤ h
ij
(k)≤ 1 for all i, j, k. (3.7)
47
For CM and DCM, the static gains correspond to the parameters (λ
ij
), so the
constraints apply to them directly; that is,
0≤ λ
ij
≤ 1 for all i, j,. (3.8)
For these models, the linear constraints on parameters can be integrated into the
optimization procedures to have a more reasonable and reliable estimation of model
parameters. Generally speaking, the linear inequality constraints will form a convex
set, and if the estimation criterion (2.36) is convex on the parameter vector θ,we
can solve the optimization by some general methods from convex optimization [7].
But the optimization work is much too involved. In our work, we extend a simple
methods provided by [63]. Taking CM for example, we first assume the τ is given,
then the criterion function (2.36) becomes linear on the remaining parameters λ
ij
.
So given the linear constraints on λ
ij
, we can perform quadratic programming
to minimize (2.36). Then we iterate on different possible τ to find the smallest
parameter set which has the minimum value of (2.36). For DCM, we can use
a similar procedure to estimate the model parameters. For FIR model, because
of its “linear-in-the-parameter” property, it is particularly easy to use the linear
constraints. Linear objective function with linear constraints can be solved by
quadratic programming method [7].
The “grey-box” modeling approaches are particularly preferable when the pro-
vided system information is very limited, so that not many training data are avail-
able. In these more challenging situations, the black-box approaches often result in
models that violate the constraints imposed by some common physical knowledge.
48
On the contrary, the static gains of the M-ARX model cannot be expressed as
a linear combination of its parameters, which makes these kinds of physical con-
straints much harder to integrate into the M-ARX modeling procedure. For these
models, we can only examine the appropriateness of the models by checking the
system behaviors in order to determine whether it satisfies the physical constraints.
49
3.5 Model Comparison
In this section, we perform a comparative analysis of all proposed predictive mod-
els (including the state-of-the-art CM work). For all models, the “Grey-box” .ap-
proaches are also included if possible. The simulation results, together with some
theoretical analysis, leads to some general suggestions for using the predictive mod-
els for field applications in the next section.
3.5.1 Number of Parameters
The number of model parameters plays a crucial role for the quality of the model.
Basically we would like to include more unknown parameters into the model in
order to describe all possible responses for the reservoir, but increasing the number
of unknowns also increases the uncertainty for parameter estimation. In particular,
the quality of the model may be more easily affected by noise in the observed
data. In general, the best model structure for identifying a system is a trade-off
between [41]:
• Flexibility: The model structure should provide good capabilities of describing
different possible systems. This can be achieved by adding more parameters
or putting the parameters in “strategic positions”.
• Parsimony: Do not use unnecessarily parameters; the model parameterization
should be kept as parsimonious as possible.
Table 3.7 shows the number of parameters used for different interwell models
with an example of 5 injectors and 4 producers. From the table, we can see that
the FIR model has the largest degree of flexibility with a large number of unknown
parameters, while other models have a comparable number of parameters. This
50
Model # of Parameters 5-inj/4-pro
CM M × N × 2 40
DCM M × N × 3 60
z-domain M × N × 2 40
FIR M × N × L 600 (L = 30)
1
st
order M-ARX M × N + N
2
36
Table 3.7: Number of parameters used for different models.
means that the FIR model should have worse noise performance than other models,
as will be shown by a noise sensitivity analysis in the Section 3.6.1.
3.5.2 Prediction Performance Comparisons
For both theoretical and practical reasons, different predictive models should be
evaluated and compared by their prediction ability on observed data that has not
been used for training. This is a general methodology known as model validation
on fresh data set [41,57] or out-of-sample forecasts in time series analysis [9,10].
The quadratic criterion is commonly used for mathematical reasons. Assume
we have trained an interwell model with some estimated parameters, denoted as M.
Let the predicted production rates on the data set not used for training at time k
be denoted as ˆ y(k|M), then the prediction error can be expressed as
J
squ
(M)=
1
T
T
k=1
y(k)− ˆ y(k|M)
2
(3.9)
51
where T is the number of data points in the validation data set. Although this
is mathematically convenient, for practical field applications it is often more ap-
propriate to define an error metric that captures the relative error. Defining the
absolute difference between the predicted and actual rates:
J
abs
(M)=
1
TN
T
k=1
N
j=1
|y
j
(k)− ˆ y
j
(k|M)| , (3.10)
a normalized performance measurement R
abs
can be defined as
R
abs
=
J
abs
(M)
1
TN
T
k=1
N
j=1
|y
j
(k)|
. (3.11)
This measurement can be interpreted as the “average ratio of errors in prediction”
and makes it easy for reservoir engineers to judge if the trained model is “good
enough” to use for prediction. We will use (3.11) as the criterion function for our
comparisons.
3.6 Comparison Results
We evaluated and compared different models using the simulation data under dif-
ferent scenarios, as indicated in Section 3.2. The grey-box extensions of CM, DCM
and FIR model, which impose physical constraints on static gains, are included in
our comparisons. In the simulation, we denote the original CM and DCM as the
“Unconstrained CM” and “Unconstrained DCM”, respectively. The grey-box CM
and DCM, we denote as “Constrained CM” and “Constrained DCM”, respectively.
For the FIR model, we only simulate the grey-box case, which incorporates the
52
400 600 800 1000 1200 1400 1600
0
1
2
3
4
5
6
x 10
−3
Training Period (Day)
R
abs
P − 1
400 600 800 1000 1200 1400 1600
0
1
2
3
4
5
6
x 10
−3
Training Period (Day)
R
abs
P − 2
400 600 800 1000 1200 1400 1600
0
1
2
3
4
5
6
x 10
−3
Training Period (Day)
R
abs
P − 3
Constrained CM
Constrained DCM
FIR
1st−order M−ARX
Unconstrained CM
Unconstrained. DCM
Figure 3.12: Performance evaluation of different models on the Scenario C. All
models have reasonably good predictions in this case. Note that P − 1, P −2and
P − 3 represent the producers 1, 2 and 3, respectively.
physical constraints (3.6) and (3.7) into the model. This is because in some cases,
the FIR model without constraints has much worst results than other approaches.
We first show the estimation results using the data from Scenario C (homoge-
neous case). The simulation time is 2835 days. Various training periods are chosen,
and the days after the 1800th are selected as the validation period. The predic-
tion results are evaluated via R
abs
for all producers, as shown in Fig. 3.12. The
results show that all models have reasonably good predictions in this case, with the
maximum R
abs
≈ 0.006, i.e., average prediction errors of about 0.6%.
Then we turn to estimation results of the Scenario D (multiple fractures case).
The simulation time are the same as we have done in Scenario C, and the results
are shown in Fig. 3.13. In this scenario, the 1
st
-order M-ARX model performs
53
400 600 800 1000 1200 1400 1600
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
Training Period (Day)
R
abs
P − 1
400 600 800 1000 1200 1400 1600
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
Training Period (Day)
R
abs
P − 2
400 600 800 1000 1200 1400 1600
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
Training Period (Day)
R
abs
P − 3
Constrained CM
Constrained DCM
FIR
1st−order M−ARX
Unconstrained CM
Unconstrained. DCM
Figure 3.13: Performance evaluation of different models on the Scenario D. From
the results, the M-ARX model performs the best for all training periods of all
producers, and (constrained) FIR model perform almost in the second place. Note
that P − 1, P −2and P − 3 represent the producers 1, 2 and 3, respectively.
the best for all training periods of all producers, and (constrained) FIR model
perform almost in the second place. The constrained CM and DCM perform slightly
better than unconstrained ones, but not too much gain is achieved when taking into
consideration physical constraints. If we take producer 2 for example, the M-ARX
model achieve R
abs
≈ 0.01 to 0.02, which is less than half the R
abs
achieved by CM
or DCM.
The next simulation is using the setting of Scenario A. Its simulation time is
3040 days, and various training periods of 500, 1000, 1500 days are chosen and
the period after 2000 days is used as the validation period. The R
abs
of prediction
results for all producers are shown in Fig. 3.14. From the results, all models have
good predictions and the maximum R
abs
value is about 0.03, which means the
54
500 1000 1500
0
0.005
0.01
0.015
0.02
0.025
0.03
Training Period (Day)
R
abs
P − 1
500 1000 1500
0
0.005
0.01
0.015
0.02
0.025
0.03
Training Period (Day)
R
abs
P − 2
500 1000 1500
0
0.005
0.01
0.015
0.02
0.025
0.03
Training Period (Day)
R
abs
P − 3
500 1000 1500
0
0.005
0.01
0.015
0.02
0.025
0.03
Training Period (Day)
R
abs
P − 4
Constrained CM
Constrained DCM
FIR
1st−order M−ARX
Unconstrained CM
Unconstrained. DCM
Figure 3.14: Performance evaluation of different models for Scenario A. There is no
model that can always outperform the others. Note that P − 1, P − 2, P −3and
P − 4 represent the producers 1, 2, 3 and 4, respectively.
maximum average prediction-errors are about 3%. If we compare the performance
of different models, there is no model that can always outperform the others, and
all are good enough to satisfy most applications.
The last simulation is for the Scenario G (variable injection case). The simula-
tion time is about 2160 days. The training periods are chosen as 600, 1000, 1400
and the period after 1600 days is used as the validation period. For the results, the
unconstrained CM and DCM are much worse than other models, sometime leading
to totally useless predictions. For example, the R
abs
for producer 2 using the un-
constrained CM is about 33% average prediction errors. This is also the case for
modeling via unconstrained DCM. If we consider the static gains (interwell connec-
tivities), we will find that there are many unreasonable values, such as values much
55
600 800 1000 1200 1400
0
0.02
0.04
0.06
Training Period (Day)
R
abs
P − 1
600 800 1000 1200 1400
0
0.02
0.04
0.06
Training Period (Day)
R
abs
P − 2
600 800 1000 1200 1400
0
0.02
0.04
0.06
Training Period (Day)
R
abs
P − 3
600 800 1000 1200 1400
0
0.02
0.04
0.06
Training Period (Day)
R
abs
P − 4
600 800 1000 1200 1400
0
0.02
0.04
0.06
Training Period (Day)
R
abs
P − 5
600 800 1000 1200 1400
0
0.02
0.04
0.06
Training Period (Day)
R
abs
P − 6
Constrained CM
Constrained DCM
FIR
1st−order M−ARX
Figure 3.15: Performance evaluation of different models for Scenario G. From the
results, the 1
st
M-ARX model performs worse for most of the situations. Note that
P − 1, P − 2, P − 3, P − 4, P −5and P − 6 represent the producers 1, 2, 3, 4, 5
and 6, respectively.
higher than 1 or negative values. For the performance, we only show constrained
CM and DCM, together with (constrained) FIR model and 1
st
order M-ARX model,
which is shown in Fig. 3.15. As the results shows, the 1
st
M-ARX model performs
worse for most of the situations. This is because the 1
st
M-ARX model is the only
one that does not integrate the physical constraints into the model, while other
“constrained” models are built taking this information into consideration. We will
discuss why this affects the prediction results for this case in Section 3.7.
56
3.6.1 Noise Sensitivity Analysis
In real applications, the measured rates often are subject to measurement errors,
especially on the producers. To address this kind of problem, we perform a noise
sensitivity analysis by adding different levels of white noise (uncorrelated, zero mean
and normally distributed) into the production rates and evaluating the performance
of different models from this artificial noisy data. Part of this analysis has been done
in previous literature, see [64] [33] for details. In this work, we extend this prior
work to provide a general evaluation of noise sensitivity for all models. While the
white noise does not correspond to what the observed in the real field data, using
this analysis can still provides us about the robustness of various models. Because
different models have different parameters, we cannot compare the models only
by measuring the difference of estimated parameters between noiseless and noisy
data. To solve this, we choose the static gains (interwell connectivities or injector-
producer relationships) as the comparison basis; that is, we calculate the absolute
differences of the static gains estimated by various models between noiseless and
noisy data and plot them together.
For the analysis, we choose the data from Scenario A. We add white noise of
different energy levels to the production rates. The noise level is measured by the
signal-to-noise ratio (SNR) in dB, which is defined as
SNR(dB)=10· log
P
signal
P
noise
(3.12)
where P
signal
and P
noise
are the signal and noise powers, respectively. Because of
the randomness of the noise, we perform 50 realizations (with respect to the white
noise) and obtain the average measurement errors. We also include the FIR model
57
1.5 2 2.5 3 3.5 4 4.5 5
0.04
0.05
0.06
0.07
0.08
0.09
0.1
SNR (dB)
Absolute Difference of Static Gain
P − 4
Constrained CM
Constrained DCM
Constrained FIR
1st−order M−ARX
Unconstrained CM
Unconstrained DCM
Unconstrained FIR
Figure 3.16: Noise sensitivity analysis for different models for Scenario A. The
comparisons are made by calculating the average absolute differences of the static
gains estimated under noiseless and noisy production rates data.
without using any physical constraints, which is denoted as “unconstrained FIR”
model. The results are shown in Fig. 3.16.
From the results, the constrained approaches of a model all outperform the
unconstrained version of same model. If we compare unconstrained FIR, CM,
DCM and 1
st
order M-ARX model, or the constrained CM, DCM and FIR model,
in both cases the FIR model has the worst performance. This is because the FIR
model uses many more parameters than the other models, which increases the
estimation uncertainty under noisy environment. The constrained FIR model has
similar performance to that of the 1
st
order M-ARX model, which again shows the
power of grey-box modeling approaches. Of course, both constrained CM and DCM
have the best performance.
58
3.7 Discussions
3.7.1 Drawbacks of Purely Statistical Methods
All the interwell models presented in this thesis can be considered as statistical
modeling approaches: we try to understand the whole system from only the statis-
tical analysis of the data sent to the system and the data generated by the system.
For these kinds of approaches, if the measured input and output data are not“infor-
mative” enough to provide sufficient information in order to build reliable models,
the models may not able to track the “true” behavior of the system. In this work,
the injection rates of Scenario G are different from those of Scenario A-D because
they are much more collinear (see [63], [34] and Section 6 for more discussions of
collinearity). If we check the static gains estimated by M-ARX model in this case,
we will find that there are some negative values between some well pairs. But the
static gains estimated by 1
st
M-ARX model in other cases are all between 0 and 1.
In the Scenario G, the designed injection rates are all constant for most of the time,
and only change at some time instants. The statistical modeling approach can only
gain information when the rates are fluctuating, which means injection data with
constant rates is not helpful for building the models. That is the reason why while
the 1
st
M-ARX model has good performance in other cases, but it is not so suitable
for Scenario G.
In real field applications, it is not uncommon to face a situation where the given
injection and production rates are not informative enough, e.g., the injection rates
are often kept relatively constant during normal operations or the training data
available are very limited. For all cases, we suggest that it is always better to use
the “grey-box” modeling approach by combining physical constraints. For some
models it may be difficult to combine the physical constraints into the optimization
59
procedure, e.g., as in the M-ARX model, it is always a good practice to double-
check whether the estimated parameters lead to reasonable reservoir characteristics
or not; if the answer is no, we should discard this model and look for other models.
3.7.2 Large Scale Problems
One of the main advantages of using predictive modeling approaches is their ability
to describe large scale systems, with computation time significantly lower than that
of traditional reservoir simulations. However, applying these predictive models to
large reservoirs still presents several challenges. In [62], a good example of applying
CM to large scale reservoirs is presented. Most of the procedures discussed in [62]
also apply to other predictive interwell models. In this section, we particularly focus
on issues associated with dealing with large scale systems.
For any models, the number of unknown parameters increases with the number
of injection and production wells in the ROI, which means we will face a large num-
ber of system unknowns when applying these techniques to large scale reservoirs.
Statistically speaking, given enough data for parameter estimation, including more
injectors into the optimization procedure of modeling will often generate a better
fit to the historic matching, but sometimes results in an unrealistic and unreliable
model. This is because the built models are statistical (even with some physical
constraints) without involving any geological considerations. For example, if by
coincidence a particular pair of injection and production wells has some common
distinct changes during a period of time, they will have relatively large connectivity
(static gain) even if they are far away from each other in the reservoir. Common
knowledge of the reservoir often assumes that two wells should have very limited
interactions if they are far away from each other, except if there exist some special
60
fractures connecting them. To reduce these kinds of errors, we need to decide a
suitable region of interest (ROI) for each producer (or injector). Only injectors and
producers in this ROI are considered for building models, and the connectivity of
outside wells is assumed to be zero. So now the problem becomes how to decide a
suitable ROI. There are several common approaches for this:
• Set a threshold distance for the ROI. So given the target producer (or injec-
tor), the ROI is a circle with the radius chose as the threshold. [62]
• Use a similar strategy as the first one but with the shape of the ROI chosen
as an ellipse instead of a circle. This situation accounts for some prior field
knowledge available which indicates that there are some parallel fractures
in the field. The ROI along with the directions of fractures should be the
semi-major axis.
• Use a prediction error approach. begin with the target producer, and increase
the radius (or semi-major/semi-minor axis) of ROI, which in order to include
more and more and increasing number of wells, and use these models for
production rate prediction and to calculate the prediction error. When the
prediction errors are below some threshold, we use this circle (or ellipse) as
ROI [67]
The threshold can be decided based on field experience. It often needs some ex-
perience from the reservoir engineers to decide a suitable threshold, and it depends
on the field conditions. In [62], 4000 ft was chosen for their target reservoir. The
prediction error approach is used because it can provide reasonable prediction for
the target producer. This is the approach used by [67]. In short, the prediction
error approach may get better prediction results and no field experience is needed.
61
But as a penalty, much more computation is involved and may still result in large
connectivity for well pairs that are far away from each other.
3.7.3 Suggestions for Using Predictive Models
From the above comparative analysis of different models, we can make some sug-
gestions for the use of linear predictive models:
1. Always use historic data for validation of prediction ability before applying
any predictive models. This will give us some ideas about how good the
models are.
2. Always use the “grey-box” extensions of the models when the physical con-
straints can be imported into the optimization procedures. For models with-
out grey-box extension, the results should be verified by checking some phys-
ical characteristics, e.g., the static gain. If these characteristics are not com-
patible with the physical reality, we need to discard the model and look for a
new one.
3. In terms of computational complexity, we can use the 1
st
M-ARX model for
fast evaluation. Because of the “linear-in-the-parameter” property, the pa-
rameters in M-ARX model can be easily estimated by linear regression tech-
nique. Then we check the physical characteristics and the prediction errors
via cross-validation. If the results are compatible with the physical reality
but the prediction errors remain high, it is likely that the linear assumption
is not suitable for the field. We should change to other modeling approaches
or try the LTV approaches.
62
4. If the measurement errors of the rate information are high or the injection
rates are highly collinear, the “constrained” CM/DCM approaches are prefer-
able as they require fewer parameters and are more robust.
3.8 Summary
In this chapter, we discussed the verification and comparisons of predictive mod-
els. For verification, there are two schemes used in this work. The first one is
verification via prediction ability on fresh data set. The second one is via the inter-
pretation of interwell connectivity against ground truth (synthetic data) or some
geological features (field data). We also discussed the so-called“grey-box”approach
of predictive models, the advantages of using it are investigated. Finally, we de-
fined a prediction-error based metric and compare all proposed predictive models.
Some suggestions for reservoir engineers are made from the results of comparative
analysis.
63
Chapter 4
Linear Modeling Framework
Having evaluated and compared different models, some theoretical questions about
statistical modeling remain open: what is the relation between the different models?
Is there a unified framework that can include all models? Are there alternative
models with similar or better performance?
In this chapter, we develop a linear modeling framework that allows us to in-
tegrate all predictive models proposed to date. Different predictive models are ex-
pressed as special cases of this framework. Also, the relationships between different
models become much clearer when we use this framework.
4.1 Linear Models for Well Interactions
Suppose we have M injection wells with the sequence of injection rates at the i-th
injector denoted as u
i
(k), and N production wells with the sequence of production
rates at the j-th producer denoted as y
j
(k), where k denotes the time index for
the measurement data of rates. We also assume the system behavior can be ap-
proximated as being linear. If the system response does not change for a period of
64
time, i.e., it can be seen as time-invariant, the linear model can be expressed using
a rational transfer function via z-transform:
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Y
1
(z)
.
.
.
Y
N
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ =
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ C
11
(z)
D
11
(z)
···
C
M1
(z)
D
M1
(z)
.
.
.
.
.
.
.
.
.
C
1N
(z)
D
1N
(z)
···
C
MN
(z)
D
MN
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ U
1
(z)
.
.
.
U
M
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , (4.1)
which can be written in matrix form as:
Y(z)= G(z)U(z) (4.2)
where U
i
(z)and Y
j
(z), i =1, ..., M and j =1, ..., N, are the z-transforms of
u
i
(k)and y
j
(k), respectively. C
ij
(z)and D
ij
(z) are polynomial functions of z.
Another commonly used description is the matrix fraction description (MFD) for
multivariate systems [28], which is analogous to the rational transfer function of
univariate (single-input and single-output) systems. In the MFD, (4.2) can be
expressed as:
Y(z)=
D
L
(z)
−1
N
L
(z)
U(z) (4.3)
= U(z)
N
R
(z)D
R
(z)
−1
(4.4)
where all D
L
(z), N
L
(z), N
R
(z)and D
R
(z) are all matrices where each element is
a polynomial function of z. (4.3) is called left MFD, and (4.4) is called right MFD.
D
L
, D
R
are analogous to the denominator in the univariate case, and similarly N
L
,
N
R
are analogous to the numerator. We will use MFD for model description when
it is more convenient.
65
Although various predictive modeling approaches have been proposed since 2006
[64] [39] [35] [33], we propose for the first time that all these predictive modeling
for describing the interactions between wells can be seen as special cases of the
linear time-invariant (LTI) models of (4.2) or (4.3). Although in this thesis only
LTI systems are discussed, all the models can be extended to linear time-varying
(LTV) systems, where now the parameter estimation methods become recursive in
nature [41] and continuously update the current estimate. The transfer function
representation of LTV systems are much more involved and out of the scope for this
paper. For a reference about more details on the transfer function representation
of LTV systems, see [4].
Also we need to note that although the behavior of many physical systems,
including reservoirs, is nonlinear, it is still useful to approximate the system first
using linear models. As we will show in the simulation results, the numerical flow
simulator we use [42] makes no linearity assumptions, but its results can be rea-
sonably matched and predicted using some LTI modeling approaches. Moreover,
the LTV system is related to linearization of a nonlinear system around a certain
trajectory [41], although we do not exploit the possibility of tracking nonlinear be-
havior in this work. In this work, we mainly focus on the linear modeling for LTI
systems.
4.2 A Framework for Predictive Models
In this section, we focus on putting the models introduced in Chapter 3 into the
transfer function matrix representation to show that they are actually a special
case of linear modeling approach. Using this linear modeling framework, we can
66
also extend the existing models to higher order and interpret the model parameters
easily.
4.2.1 Capacitance Model
As introduced in Section 2.2, in discrete form, the CM can be expressed as [64]
y(k)= y(0)e
−k
τ
+
1
τ
k
ξ=0
e
(ξ−k)/τ
u(ξ) (4.5)
where τ accounts for the attenuation and delays between wells. Applying the z-
transform and assuming that the initial state is negligible at the time when a
waterflood project begins, i.e., that production is essentially almost zero at that
time, we can compute the transfer function in z-domain as:
Y (z)=
1
1− e
−1
τ
z
−1
U(z), (4.6)
which shows that in CM, there is a pole located at z = e
−1/τ
.
The CM can be easily extended to the multiple wells case using the superposi-
tion principle [64]. In this case, the z-domain expression for the transfer function
becomes:
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Y
1
(z)
.
.
.
Y
N
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ =
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ λ
11
1−e
−1
τ
11 z
−1
···
λ
M1
1−e
−1
τ
M1 z
−1
.
.
.
.
.
.
.
.
.
λ
1N
1−e
−1
τ
1N z
−1
···
λ
MN
1−e
−1
τ
MN z
−1
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ U
1
(z)
.
.
.
U
M
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (4.7)
where λ
ij
with i=1, ..., M and j =1, ..., N are the weight factors (interwell connec-
tivity). The time constants τ
ij
, extended from single injector/producer case, play
arole as a pole (at z = e
−1
τ
ij
) between the i-th input and j-th output. This shows
67
that CM is a special case of the general LTI model (4.1), which possesses one pole
for each injector-producer well pair.
Besides, because the interpretation of CM in terms of reservoir characteristics
has been fully discussed in [63] [64] (also see Section 2.2 for a brief discussion), this
gives us a way to connect the transfer function representation to the characteristics
in the reservoir. A pole in transfer function representation can be interpreted as the
resistive and capacitance effects of the interwell region, since this is an interpretation
that has been proposed for the CM. Moreover, the location of the pole is e
−1
R
ij
C
ij
,
where R
ij
and C
ij
represent the equivalent resistive and capacitance effects between
injector i and producer j, respectively.
4.2.2 Distributed Capacitance Model
In the distributed capacitance model (DCM) introduced in Section 2.4, for the
single injector and single producer case, the production rates can be expressed as
q(k)=
c
1
e
−k
τ
1
+ c
2
e
−k
τ
2
+
1
τ
1
− τ
2
k
ξ=0
e
ξ−k
τ
1
u(ξ)−
k
ξ=0
e
ξ−k
τ
2
u(ξ)
(4.8)
where τ
1
and τ
2
account for the attenuation and delays of the regions between
wells, and are analogous to the time constant τ in CM. We can easily derive the
corresponding transfer function in the z-domain:
Y (z)=
1
(1− e
−1
τ
1
z
−1
)(1− e
−1
τ
2
z
−1
)
U(z) (4.9)
This shows that DCM leads to a transfer function with two poles at z = e
−1
τ
ij,1
and
z = e
−1
τ
ij,2
in order to describe the injector-to-producer relationships. Comparing
(4.9) to (4.6), it is clear that DCM is equivalent to cascading two CMs, with time
68
constants τ
1
and τ
2
, respectively [32]. Similar to CM, the DCM can be extended
to the multiple wells case using the superposition principle, so that the transfer
function becomes
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Y
1
(z)
.
.
.
Y
N
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ =
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ λ
11
(1−e
−1
τ
11,1
z
−1
)(1−e
−1
τ
11,2
z
−1
)
···
λ
M1
(1−e
−1
τ
M1,1
z
−1
)(1−e
−1
τ
M1,2
z
−1
)
.
.
.
.
.
.
.
.
.
λ
1N
(1−e
−1
τ
1N,1
z
−1
)(1−e
−1
τ
1N,2
z
−1
)
···
λ
MN
(1−e
−1
τ
MN,1
z
−1
)(1−e
−1
τ
MN,2
z
−1
)
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ U
1
(z)
.
.
.
U
M
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (4.10)
where λ
ij
with i=1, ..., M and j =1, ..., N have the same meaning as in CM and
can also be interpreted as the interwell connectivities in the reservoir. The time
constants τ
ij,1
and τ
ij,2
are extended from single well pair case and model the two
poles at z = e
−1
τ
ij,1
and z = e
−1
τ
ij,2
between the i-th input and j-th output. In general,
the DCM behaves as a second-order (two poles) LTI system for each input-output
pair, where the numerators are described by a constant (zero-order) now.
4.2.3 Double Pole Model [39]
The double pole model (initially introduced in [39] with the name of “z-domain
model”) is a two parameter auto-regressive model for the single well pair case. In
this case, the transfer function of the model can be expressed in z-domain as
H(z)=
γz
−1
(1− αz
−1
)
2
(4.11)
where γ and α are the unknown parameters. Obviously, this model describes the
system behavior with a double pole at α and a weight coefficients γ. Similar to CM
69
and DCM, this model can be extended to the multiple wells case by the superposi-
tion principle and expressed as
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Y
1
(z)
.
.
.
Y
N
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ =
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ γ
11
(1−α
11
z
−1
)
2
···
γ
M1
(1−α
M1
z
−1
)
2
.
.
.
.
.
.
.
.
.
γ
1N
(1−α
1N
z
−1
)
2
···
γ
MN
(1−α
MN
z
−1
)
2
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ U
1
(z)
.
.
.
U
M
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (4.12)
If we compare (4.12) with (4.10), we can see that double pole model is a special
case of DCM, where now the two poles in DCM overlap and become a double pole
(e
−1/τ
ij,1
= e
−1/τ
ij,2
= α
ij
). The weight factors γ
ij
in double pole model play similar
roles as λ
ij
in CM and DCM, all indicating the static gains of the system (different
names are used by different authors: interwell connectivity in CM [64] and DCM;
injector-producer relationship in double pole model [39] [67]). This also means that
the z-domain model can be interpreted and related to some reservoir characteristics
in a similar way as the DCM.
4.2.4 Finite-Impulse-Response Model
As the finite-impulse-response (FIR) model introduced in Section 2.3, the produc-
tion rates of j-th producer can be expressed as
y
j
(k)=
M
i=1
u
i
(k)∗ h
ij
(k) (4.13)
where h
ij
(k)with k =0, 1, ..., L− 1 denotes the impulse response of producer j
from the i-th injector, which describe the response curve between this input-output
pair. The shape of the impulse response can be easily interpreted in terms of some
characteristics of the reservoir, see [35] and Section 2.3 for more details.
70
For the multiple wells, the transfer function of FIR model is expressed as
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Y
1
(z)
.
.
.
Y
N
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ =
⎡ ⎢ ⎢ ⎢ ⎢ ⎣
L−1
k=0
h
11
(k)z
−k
···
L−1
k=0
h
M1
(k)z
−k
.
.
.
.
.
.
.
.
.
L−1
k=0
h
1N
(k)z
−k
···
L−1
k=0
h
MN
(k)z
−k
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ U
1
(z)
.
.
.
U
M
(z)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , (4.14)
which is a special case of (4.1). This is easily seen by setting C
ij
(z)equals to
h
ij
(0) + ... + h
ij
(L− 1)z
−L+1
and D
ij
(z) = 1 for all i=1, ..., M and j =1, ..., N.
So the FIR model belongs to general LTI system with no poles in each well pair.
Moreover, the FIR model can approximate any LTI system provided that the
impulse length L is long enough. To see this, we know that each pole can be
expressed as an infinite series. For example, in the discrete model, suppose we have
apoleat z = α with |α| < 1, this pole can be expressed as
1
1− αz
−1
=1+ αz
−1
+ α
2
z
−2
+ α
3
z
−3
+ ... (4.15)
So each fractional polynomial of z for the i-th input and j-th output can be ex-
pressed as
C
ij
(z)
D
ij
(z)
=
C
ij
(z)
(1− α
1
z
−1
)...(1− α
p
z
−1
)
= C
ij
(z)(1 + α
1
z
−1
+ α
2
1
z
−2
+ ...)(...)(1 + α
p
z
−1
+ α
2
p
z
−2
+ ...) (4.16)
= h
ij
(0) + h
ij
(1)z
−1
+ h
ij
(2)z
−2
+ h
ij
(3)z
−3
+ ...
here we assume D
ij
(z) was expressed by p poles α
1
, ..., α
p
and the last equality holds
because C
ij
(z) is also a polynomial of z. For a stable system, the coefficients h
k
converge to zero asymptotically, with the rate depending on the dominant pole (the
absolute value of the pole which is closest to the unit circle). Suppose we truncate
71
the infinite series at k = L− 1, this new set of coefficients {h
i
j(k)|k=0, ..., L− 1}
are the unknown parameters in (4.13), which will be determined from the training
data. So this justifies that the FIR model can approximate any possible response
within some estimation errors provided the FIR length L is long enough. But the
penalty is also obvious: many more unknown parameters are needed than for other
models.
4.2.5 Multivariate ARX Model
As for the Multivariate Auto-Regressive with eXogenous (M-ARX) model intro-
duced in Section 2.5 (here we all refer to first-order M-ARX model) for describing
injection and production relationships, it can be expressed as
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ y
1
(k)
.
.
.
y
N
(k)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = A
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ y
1
(k− 1)
.
.
.
y
N
(k− 1)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + B
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ u
1
(k− 1)
.
.
.
u
M
(k− 1)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (4.17)
where A and B are N by N and N by M coefficient matrices, respectively, which
are independent of time index k for a LTI system. (4.17) shows that the production
rates at time k can be expressed as a linear combination of both injection and
production rates at time k− 1. In z-domain, (4.17) becomes
Y(z)= Az
−1
Y(z)+ Bz
−1
U(z) (4.18)
Then the transfer function matrix representation becomes
Y(z)=
I− Az
−1
−1
Bz
−1
U(z) (4.19)
72
This expression is obviously the left MFD (4.3) with D
L
(z)= I−Az
−1
and N
L
(z)=
Bz
−1
. In other words, the 1
st
order M-ARX model is a special case of (4.3), with
now all elements in denominator matrix D
L
(z) with maximum order of z
−1
.
As we will show in the Section 5.1, when the coefficient matrix A is diagonal,
the 1
st
-order M-ARX model will reduce to the capacitance-resistive model with one
time constant for each producer (CRMP), which is a special case of CM [54]. In
general, the off-diagonal terms of matrix A are not zeros, which means the 1
st
-order
M-ARX model uses the production rates from other producers, together with the
injection rates, in order to estimate the future value of the production rates.
4.2.6 Higher Order Models
By putting all interwell models into the transfer function matrix representation, we
can easily extend the interwell models to higher orders. When we try to character-
ize the reservoir behavior with higher order LTI models (4.2) or (4.3), we basically
introduce more poles to describe the system. From the CM [64] and (4.6), the loca-
tion of each pole can be interpreted as some capacitance and resistive effects in the
reservoir. So when higher order models are involved, this can be seen as modeling
the reservoir by cascading a series of small fluid-capacitors and fluid-resistors, which
implies high heterogeneity along the path of a particular well pair. Although high
order models have the potential to describe a broader range of possible systems,
there is also a penalty for using them: the number of unknown parameters in these
models increases. For more discussions on issues related to number of parameters
issue, please refer to Section 3.5.1.
73
Table 4.1: Different predictive models and the characteristics of their transfer func-
tion.
Model Transfer func. Pole # and location
CM
λ
ij
1−e
−1
τ
ij
z
−1
(each pair)
One at z = e
−1
τ
ij
DCM
λ11
(1−e
−1
τ
ij,1
z
−1
)(1−e
−1
τ
ij,2
z
−1
)
(each pair)
Two at z = e
−1
τ
ij,1
and z = e
−1
τ
ij,2
Double Pole Model
γ11
(1−α11z
−1
)
2
(each pair) Two at z = α
ij
(double-pole)
FIR Model
L−1
k=0
h
ij
(k)z
−k
(each pair) No Pole
M-ARX Model (I− Az
−1
)
−1
Bz
−1
Depends on A
4.3 Conclusion
In this chapter, we use the general transfer function for LTI system to show that each
model belongs to a special case of this model, which means all proposed models are
linear models. Moreover, the relationship between different models is investigated,
so the interpretation of CM for reservoirs can be easily extended to all models. In
summary, Table 4.1 summarizes all proposed models discussed in this chapter.
74
Chapter 5
Prediction Under Controlled Producers
In waterflood projects, the number of producers sometimes changes, e.g., a producer
is shut-in or a new producer is added. To analyze these data sets using conventional
predictive models, which assume a fixed number of producers, we would need to
divide the whole data period into several smaller intervals, so that the number of
producers remains constant in each interval. This is because when the number of
producers changes, it leads to a totally different input-output relationship, so that
the model parameters need to be re-estimated. Thus, every time a producer is
shutting in, we need a new data set to retrain the model parameters.
In this chapter, we investigate the use of M-ARX models for this kind of situ-
ation. The M-ARX model can reduce the number of parameters that needs to be
retrained as compared to other predictive models. More importantly, the M-ARX
model can predict the performance in the shut-in case, only requiring that the pro-
ducer be set to a constant rate for a while. This makes it much more practical to
control producers in order to predict performance under several possible “What if”
scenarios.
75
I - 1
P - 1
I - 2
P - 2
1
τ
2
τ
I - 1
P - 1
I - 2
P - 2
1
τ
2
τ
Figure 5.1: The CRMP and M-ARX model in 2-injectors/2-producers scenario. The
black arrows denote the contributions from other producers, which are considered
in 1
st
order M-ARX model but not in CRMP.
5.1 Interpreting M-ARX Model
In order to describe the producer-to-producer relationships included in 1
st
order
M-ARX model, we need to understand the relationship between M-ARX model
and CM. We first consider a simplified version of CM, named capacitance-resistive
model with one time constant for each producer (CRMP), which was proposed
in [54]. This model uses one time constant τ for each producer by indicating that
for each producer, the pore volume shared with any injectors should be similar,
since it corresponds to the region surrounding the well, and therefore they should
have similar time constants. The CRMP was proposed with the goal of reducing
the number of unknown parameters and simplifying the optimization procedure,
while still achieving reasonable prediction results (see [54] for more discussions).
Here we will show that the 1
st
order M-ARX model can be interpreted as a
CRMP that takes into consideration producer-to-producer interactions. To simplify
the derivation, suppose we have two injectors and two producers in the ROI, as
76
shown in Fig. 5.1. Under the CRMP model, the production rates of each producer
can be expressed as:
y
j
(k)= y
j
(0)e
−k/τ
j
+
2
i=1
λ
i
u
i
(k)·
1
τ
j
e
−k/τ
j
(5.1)
Applying the z-transform in (5.1), we get:
⎡ ⎢ ⎣ Y
1
(z)
Y
2
(z)
⎤ ⎥ ⎦ =
⎡ ⎢ ⎣ λ
11
z
−1
1−e
−1
τ
1 z
−1
λ
21
z
−1
1−e
−1
τ
1 z
−1
λ
12
z
−1
1−e
−1
τ
2 z
−1
λ
22
z
−1
1−e
−1
τ
2 z
−1
⎤ ⎥ ⎦ ⎡ ⎢ ⎣ U
1
(z)
U
2
(z)
⎤ ⎥ ⎦ =
⎡ ⎢ ⎣ 1− e
−
1
τ
1
z
−1
0
01− e
−
1
τ
2
z
−1
⎤ ⎥ ⎦ −1
⎡ ⎢ ⎣ λ
11
z
−1
λ
21
z
−1
λ
12
z
−1
λ
22
z
−1
⎤ ⎥ ⎦ ⎡ ⎢ ⎣ U
1
(z)
U
2
(z)
⎤ ⎥ ⎦ (5.2)
Recalling the M-ARX model expression (3.3) and comparing it with (5.2), this
means that when the coefficient matrices A and B are as follows:
A =
⎡ ⎢ ⎣ e
−
1
τ
1
0
0 e
−
1
τ
1
⎤ ⎥ ⎦ , (5.3)
B =
⎡ ⎢ ⎣ λ
11
λ
21
λ
12
λ
22
⎤ ⎥ ⎦ . (5.4)
the 1
st
order M-ARX model will reduce to the CRMP model, and the coefficients
in matrices A and B can be interpreted the same way as in the CRMP model.
A being a diagonal matrix can be interpreted by stating that no improvements in
modeling accuracy can be achieved by using production rates from other producers
in order to estimate production for a given producer, which may happen when the
producers are far away from each other (thus with little influence on each other).
77
Although only shown in the 2-injectors and 2-producers case, the result can be
easily extended to cases involving more injectors and producers.
In general, the off-diagonal terms in A are not zero, and they can be interpreted
as representing the contribution from the production rates of other producers. This
implies that in the 1
st
order M-ARX model, the producer-to-producer interactions
can be analogous to the injector-to-producer interactions in the CRMP model. To
see this, assume the matrix A in (5.3) has non-zero off-diagonal terms a
11
and a
21
,
then (5.2) in the time-domain can be expressed as
⎡ ⎢ ⎣ y
1
(k+1)
y
2
(k+1)
⎤ ⎥ ⎦ =
⎡ ⎢ ⎣ e
−
1
τ
1
z
−1
a
12
a
21
e
−
1
τ
2
z
−1
⎤ ⎥ ⎦ ⎡ ⎢ ⎣ y
1
(k)
y
2
(k)
⎤ ⎥ ⎦ +
⎡ ⎢ ⎣ λ
11
λ
21
λ
12
λ
22
⎤ ⎥ ⎦ ⎡ ⎢ ⎣ u
1
(k)
u
2
(k)
⎤ ⎥ ⎦ ,
(5.5)
y
1
(k), y
2
(k) have a similar role as the inputs of u
1
(k)and u
2
(k). In short, the 1
st
order M-ARX model can be seen as a generalization of the CRMP model, where
now the producer-to-producer interactions are considered in a similar way to the
injector-to-producer relationships.
5.2 Prediction for Shut-In Producers
When the production of a target producer is forced to zero (shut-in) during some
time period, the question we want to solve is: can we predict the performance of
all other producers in the reservoir? Forecasting is important here because on the
one hand sometimes producers on daily operations for such as due to maintenance
or other reasons. On the other hand, shutting-in producer is sometimes a good
solution for waterflood management and optimization, and we need some forecasting
ability in order to determine whether a specific producer should be shut-in. When
78
some producer is shutting-in, the whole system is considered changed and most
of the predictive modeling approaches, such as CM or FIR model, can only deal
with this kind of situation by considering that the whole fluid distribution has
changed, which requires re-training all parameters in the reservoir model. In [30],
a compensated capacitance model (CCM) was proposed to reduce the number of
parameters that need to be re-trained by introducing a“pseudo-injector”at the well
that is being shut-in. In this work, the relations between the pseudo-injector and
the producers are also interpreted as producer-to-producer interactions between the
shut-in producer and the other producers. Now we will show how to use the 1
st
order M-ARX model to further reduce the data required for retraining, and also
extend the concept of pseudo-injector to deal with a more general scenario.
As shown before, the producer-to-producer interactions in the 1
st
order M-ARX
model can be interpreted as being similar to injector-to-producer interactions in
CRMP, so when a producer, e.g., the α-th producer, is artificially shutting in after
some time instant, the new production rates y
j
(t)onall otherproducers can be
expressed as:
˜ y
j
(t+1) = y
j
(t+1)− y
α
(t)∗
λ
jα
τ
j
e
−t/τ
j
(5.6)
where j =1, ..., N and j = α; y
j
(t+1) and y
α
(t) are the original production rates
(the predicted production rates without any constrained producer at producer j
and producer α, respectively.) (5.6) shows that there is a transient period when
constraints on a producer are first imposed, and each producer will have an in-
creasing rate. (Because the interaction λ
jα
between producers j and α is always
79
negative, the production rate ˜ y
j
(t + 1) will increase.) To highlight the transient
effects, ˜ y
j
(t + 1) can be rewritten as
˜ y
j
(t+1)= y
j
(t+1)+ y
α
(t)· k
j
1− e
−(t−t
sh
)/τ
j
(5.7)
where t
sh
denotes the time at which the shut-in begins, and the coefficients k
j
=
−λ
jα
can be interpreted as the influence of producer α on producer j.Astime t
increases, e
−(t−t
sh
)/τ
j
becomes close to zero, so the ˜ y
j
can be approximated by
˜ y
j
(t+1)= y
j
(t+1)+ y
α
(t)· k
j
(5.8)
We can use (5.7) when the data are in transient period (from the simulation data,
usually t− t
sh
<
20
T
where T is measured in days) or (5.8) with the data after the
transient period. Note that (5.7) and (5.8) can also be derived based on pseudo-
injector concepts [30], but a 1
st
order M-ARX model is being used instead of CM.
For (5.7), only two parameters need to be retrained for each producer (k
j
and τ
jα
),
and if (5.8) is used, only k
j
needs to be estimated. Compared to CCM, which needs
to estimate (M +1)× N unknown parameters (a smaller number of parameters
than required by other models), the 1
st
order M-ARX only needs to re-estimate
(N−1)×2 parameters. Table 5.2 summarizes the number of parameters that need
to be retrained for different predictive models.
5.2.1 Simulation Results
In the simulations, we show the application of M-ARX model to the producer shut-
in case, where it can be used for performance prediction with a very short re-training
period after the shut-in has begun.
80
Table 5.1: The number of parameters needed to be retrained when a producer is
shut-in.
# of Parameters 5-injectors/4-producers Case
CM M × (N − 1)× 2 30
FIR Model M × (N− 1)× L 450 (L = 30)
DCM M × (N − 1)× 3 45
Compensated CM (M +1)× (N − 1) 18
M-ARX Model (N − 1)× 2 6
For the synthetic data, we use the Scenario A introduced in Section 3.2, which is
showninFig.3.1. Inthiscase, P4 is shut-in from the 730-th day to the 1824-th day.
Because P4 has the maximum production rates among all producers before shut-in,
as the simulation data show, the production rates at all other producers increase
significantly during the P4 shutting-in period. We use the prediction procedure
built on M-ARX model (5.7) to predict the gross production rates during the shut-
in period with different training periods after the shut-in begins. The prediction
results are also evaluated via the R
2
measurement, which is defined as
R
2
=
J
2
(M)
1
T
T
k=1
y(k)
2
(5.9)
where
J
2
(M)=
1
T
T
k=1
y(k)− ˆ y(k|M)
2
(5.10)
The results for P1- P3 are shown in Fig. 5.2. Also, Fig. 5.3 illustrates the prediction
results when re-training period is equal to 20 days. As the results show, we can
generate good prediction results when re-training periods are chosen larger than 15
days. This validates this procedure and demonstrates one of the applications on
M-ARX model.
81
0 5 10 15 20 25 30
0
0.05
0.1
0.15
0.2
0.25
Retraining Period
R
2
P − 1
0 5 10 15 20 25 30
0
0.05
0.1
0.15
0.2
0.25
Retraining Period
R
2
P − 2
0 5 10 15 20 25 30
0
0.05
0.1
0.15
0.2
Retraining Period
R
2
P − 3
Figure 5.2: Performance of prediction error for production rates based on M-ARX
model in Scenario A with P4 shut-in from the 705-th day to the 1800-th day. The
performance is measured via R
2
and different retraining periods are evaluated.
5.3 Prediction for Constrained Producers
The main drawback of the prediction with a shut-in producer is obvious: we need to
close this producer at least for a while in order to predict the long term performance.
This heavily limits this application for economic reasons, because any producer
shut-in can potentially lead to a decrease of oil production. To successfully handle
several“what if”scenarios, we would like to predict the performance after a producer
shut-in but requiring minimal changes in production rates so the economic impact
of adjusting production for modeling is small. This can be achieved by the use
of M-ARX model with producer set to constant rate. For example, we first limit
the producer to operate at a certain rate C (we denote the target producer with
82
200 400 600 800 1000 1200
0
500
1000
1500
Before Shut−In←→ Shut−In Period
Time(day)
Production Rate
P − 3
200 400 600 800 1000 1200
0
500
1000
1500
2000
2500
Before Shut−In←→ Shut−In Period
Time(day)
Production Rate
P − 3
200 400 600 800 1000 1200
0
500
1000
1500
2000
Before Shut−In←→ Shut−In Period
Time(day)
Production Rate
P − 3
Simulation Data
Predicted by 1
st
−order M−ARX Model
Figure 5.3: Prediction of production rates by M-ARX model in Scenario A with P4
shut-in. Note that for retraining data between two black lines (solid and dash) are
used.
artificial control, shut-in or set to some rate, as constrained producer in the following
discussions). Based on the M-ARX model, in this case (5.7) can be rewritten as
˜ y
j
(t+1)= y
j
(t+1)+(y
α
(t)− C)· k
j
1− e
−(t−t
sh
)/τ
j
(5.11)
This expression can be seen as a general case of (5.7), in which the constant rate
C was set to zero. By estimating the parameters using (5.11), we can then predict
performance under the scenario where the constrained producer can be set to any
rates, including the shut in condition. This of course gives us much more flexibility
to forecast performance given that some producers are under artificial control.
Besides, in this situation the estimated parameters k
j
can be interpreted as the
weight factors of the influence from producer α to producer j, which highlights the
83
interactions between these two producers. So we can define k
j
as the producer-
producer relationship (PPR), which denotes the weight factors that characterizes
the effective contribution of production decrease of a target production well to the
total gross production of surrounding production wells.
Note that the main novelty of this application is that one does not actually need
to shut-in a well (and thus potentially reduce overall production significantly) in or-
der to estimate the impact of a well shut-in, while all previous modeling approaches
(including CCM) need to actually do the shut-in for forecasting. We summarize
this new constrained producer method as follows:
1. Training phase: the field has been under normal operation (without any con-
strained producer) for some time period. Then we constrain the target pro-
ducer for a while (e.g., 30 days in the simulation data). Only these data are
needed for the training of M-ARX model.
2. Predicting phase: we use the trained M-ARX model to predict the perfor-
mance when the target producer is under ANY constrained rates, including
the shut-in.
5.3.1 Simulation Results
In the simulations, we show the prediction results for the constrained producer. As
the previous section, Scenario A (introduced in Section 3.2) is used, and all settings
are the same. The only difference is now we set the target constrained producer
to some constant rate after day 730. Using production rate data up to day 760
(which means 30 days are used for re-training the parameters), we can estimate the
parameters in (5.11) and use them to predict the performance when the constrained
producer is shutting in. Figures 5.4, 5.5 and 5.6 show the prediction results when
84
200 400 600 800 1000 1200 1400 1600 1800
0
100
200
300
400
500
600
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 1
200 400 600 800 1000 1200 1400 1600 1800
0
50
100
150
200
250
300
350
400
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 2
200 400 600 800 1000 1200 1400 1600 1800
0
100
200
300
400
500
600
700
800
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 3
200 400 600 800 1000 1200 1400 1600 1800
0
1000
2000
3000
4000
5000
6000
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 4
Simulation Data (Shut−In)
Data from Constrained Producer Case
Predicted by M−ARX Model
Figure 5.4: Prediction of production rates for P1 shut-in by M-ARX model in
Scenario A when P1 is set to 200 bbl/day.
P1, P2and P3 are set to constant rate 200, 100, and 200 bbl/day, respectively. For
P4, we simulate two different constant rates, 600 and 2500 bbl/day, and the results
are shown in Figures 5.7 and 5.8. All results show that the constrained producer
application is very promising.
85
200 400 600 800 1000 1200 1400 1600 1800
0
100
200
300
400
500
600
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 1
200 400 600 800 1000 1200 1400 1600 1800
0
50
100
150
200
250
300
350
400
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 2
200 400 600 800 1000 1200 1400 1600 1800
0
100
200
300
400
500
600
700
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 3
200 400 600 800 1000 1200 1400 1600 1800
0
1000
2000
3000
4000
5000
6000
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 4
Simulation Data (Shut−In)
Data from Constrained Producer Case
Predicted by M−ARX Model
Figure 5.5: Prediction of production rates for P2 shut-in by M-ARX model in
Scenario A when P2 is set to 100 bbl/day.
86
200 400 600 800 1000 1200 1400 1600 1800
0
100
200
300
400
500
600
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 1
200 400 600 800 1000 1200 1400 1600 1800
0
50
100
150
200
250
300
350
400
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 2
200 400 600 800 1000 1200 1400 1600 1800
0
100
200
300
400
500
600
700
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 3
200 400 600 800 1000 1200 1400 1600 1800
0
1000
2000
3000
4000
5000
6000
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 4
Simulation Data (Shut−In)
Data from Constrained Producer Case
Predicted by M−ARX Model
Figure 5.6: Prediction of production rates for P3 shut-in by M-ARX model in
Scenario A when P3 is set to 200 bbl/day.
87
200 400 600 800 1000 1200 1400 1600 1800
0
500
1000
1500
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 1
200 400 600 800 1000 1200 1400 1600 1800
0
500
1000
1500
2000
2500
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 2
200 400 600 800 1000 1200 1400 1600 1800
0
500
1000
1500
2000
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 3
200 400 600 800 1000 1200 1400 1600 1800
0
1000
2000
3000
4000
5000
6000
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 4
Simulation Data (Shut−In)
Data from Constrained Producer Case
Predicted by M−ARX Model
Figure 5.7: Prediction of production rates for P4 shut-in by M-ARX model in
Scenario A when P4 is set to 600 bbl/day.
88
200 400 600 800 1000 1200 1400 1600 1800
0
500
1000
1500
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 1
200 400 600 800 1000 1200 1400 1600 1800
0
500
1000
1500
2000
2500
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 2
200 400 600 800 1000 1200 1400 1600 1800
0
500
1000
1500
2000
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 3
200 400 600 800 1000 1200 1400 1600 1800
0
1000
2000
3000
4000
5000
6000
Before Shut−In←→ Shut−In Period
Time(day)
Prod. Rate
P − 4
Simulation Data (Shut−In)
Data from Constrained Producer Case
Predicted by M−ARX Model
Figure 5.8: Prediction of production rates for P4 shut-in by M-ARX model in
Scenario A when P4 is set to 2500 bbl/day.
89
Chapter 6
Injection Rates Design
For parameter estimation in dynamic systems, it is desirable to control the con-
ditions under which the data are collected. The objective of this chapter is to
investigate the problem of designing inputs so that the collected outputs are as
informative as possible with respect to the models to be built using the data. The
ultimate goal is to design a set of injection rates based on some criterion in order
to facilitate the estimation of model parameters and reservoir characteristics.
In this chapter, we first provide a brief literature review focusing on both system
identification and channel estimation in communication systems. Deterministic and
stochastic approaches are discussed and a brief comparison is made. For determin-
istic approach, the novelty of this work is that we propose a novel procedure for
designing a set of signals with zero circular cross-correlation for arbitrary shifts.
This procedure is based on a property of set of inverse-repeat signals. For the
stochastic approach, we first survey input design in system identification and ex-
tend the work in the literature to the multiple well case. Another novelty is that we
apply a well-known procedure to field applications and evaluate the performance
based on some predictive models.
90
6.1 Literature Review
In this section, we review prior work on input design within two different fields:
system identification and channel estimation in communication systems. In both
fields the goal is to design the inputs used for probing the system in order to estimate
some system characteristics. While the goals are similar, the approaches are totally
different. We will focus on the frameworks proposed in each of these fields, namely
as deterministic and stochastic approaches, respectively. Finally, some comparisons
between these two approaches are provided.
6.1.1 Optimal Input Design
Optimal input design for linear dynamic systems was first considered around 1960
( [36] is one of the earliest contributions) and it became an active area of research
in the 1970’s. Different qualify measures for the identified model were used. For
more comprehensive discussions, see Mehra [46] [47], Zarrop [66], Goodwin and
Payne [18] and Goodwin [17]. Until the 1990’s, almost all research focused on the
minimization of some measure of the variance error of the estimated quantity. The
goal was to minimize some objective functions, usually various measures of the
covariance matrix P
θ
,where θ is defined as the parameter vector of the model
that is being estimated. An open-loop structure was used in this case. Even
though some of the optimal input design work of the 1970’s considered closed-loop
experiments [48] [49], the objective functions considered at that time were limited
to functions of the covariance of the open-loop model parameters.
In the system identification literature, the prediction error method (PEM) [41]
is widely used with a full order model structure. The estimated parameter vector,
91
denoted as
ˆ
θ
N
, was proven to converge to a Gaussian distribution under some mild
assumptions:
ˆ
θ
N
−θ
0
N→∞
→ N(0,P
θ
) (6.1)
where θ
0
denotes the “true” value of the system parameters and P
θ
denotes the
asymptotic covariance matrix of parameter estimation, which can be estimated from
the data. The matrix has been used to judge the “goodness” of different designs
because it provides a measure of the average difference between the estimated and
the true value of the parameters. The classical approach has been to minimize some
scalar function of the asymptotic covariance matrix P
θ
with constraints on input
and/or output power. Examples of criteria that are commonly used are [51]:
A-optimality : minTr{P
θ
} (6.2)
D-optimality : min det{P
θ
} (6.3)
E-optimality : min λ
max
{P
θ
} (6.4)
L-optimality : minTr{WP
θ
} (6.5)
where Tr is the trace function of the matrix, det is the determinant of the matrix,
λ
max
represents the maximum eigenvalue of the matrix and W is a nonnegative
weighting matrix.
Both time and frequency domain approaches are considered in the previous
literature. For designing the signals in the time domain, the problem typically
reduces to a nonlinear optimal control problem with N free variables, in which
N is the data length [11]. The resulting complexity was one of the reasons that
motivated researchers to solve the input design problem in the frequency domain
instead. Making some assumptions, it is possible to derive nice expressions for
92
the asymptotic covariance matrix. Moreover, it is easier to interpret the results in
general for design in the frequency domain.
In short, the optimal input design for system identification is typically based
on a stochastic approach where every actual input for the system can be seen as
a realization of a process. The criterion used for input design was the spectrum
of the input signals, with performance guaranteed achievable when N is large (by
asymptotic parameter variance P
θ
). We will discuss how to apply this approach
with field applications where the goal is injection rate design later in this chapter.
6.1.2 Channel Estimation in Communication
In communication systems, the multiple-input multiple-output (MIMO) technique
has been shown to greatly increase the capacity of wireless systems, and it can fit the
growing demand for high data rates in the wireless environment. However, to use the
advantages of MIMO systems, accurate channel state information (CSI) is required.
If space-time coding is used, an accurate CSI is crucial for the performance of
decoders. Therefore, channel estimation plays a key role in MIMO communication
systems [19] [8].
Using training sequences is one of the most widely applied approaches for
MIMO channel estimation. In this approach, the channel is estimated using the
received signal resulting from a predetermined sequence, denoted as training se-
quence, being sent from the transmitter. There has been a growing interest in the
training-based MIMO channel estimation. For example, Hassibi and Hochwald [21]
linked the training sequence problem with channel capacity; Marzetta [44] consid-
ered the BLAST training using maximum likelihood (ML) method; Li [37] devel-
oped a least-square (LS) training-based channel estimation technique for orthogonal
93
frequency-division multiplexing (OFDM) systems with multiple transmit antennas;
Scaglione and Vosoughi [56] improve the LS approach via minimum mean-square-
error (MMSE) symbol estimate. As a comprehensive study, Larsson and Stoica [31]
gave a general discussion on optimal MIMO training schemes based on LS crite-
rion. Biguesh and Gershman [6] discussed the tradeoffs between LS, scaled LS, and
MMSE methods.
The training sequence design was treated as a deterministic approach, that is,
the aim is to select the value of training sequences exactly before transmission.
Some authors, e.g., [12] and [14], constraint the training sequences to be derived
from a finite alphabet (such as BPSK constellation). This work concludes that,
except for some special cases, the deterministic optimal training sequence for any
length could only be found via exhaustive search, but only possible for short se-
quences with small alphabet size because of complexity. One of the main conclu-
sions derived from the training-based channel estimation was the optimal sequences
should have an impulse-like auto-correlation and zero cross-correlation properties.
Actually, the sequence design for good correlation properties was very crucial for
some communication fields, such as code division multiple access (CDMA) system,
where Gold pair and Kasami sequences are commonly used because of their good
correlation properties. For details, see [16]. Because these designs were focused on
minimizing the maximum absolute auto-correlation and cross-correlation values,
they are not optimal for parameter estimation in dynamical systems.
94
6.2 A Novel Deterministic Approach for Input
Sequence Design
In this section, we approach the input design problem using a deterministic ap-
proach. We developed a procedure to generate a set of signals, which have vanish-
ing cross-correlation with each other with arbitrary shifts. To design these inputs,
inverse-repeat signals are first introduced and one property of these signals is dis-
covered. Then a new procedure is proposed to generate a set of signals with all zero
cross-correlation with each other with arbitrary time shifts.
Note that the procedure described in this section is also valid for continuous
signals but here we only discuss the discrete case because for our problem, datasets
(inputs and outputs) are all discrete in some time scale (day, week, or month).
6.2.1 Inverse-repeat signals
Inverse-repeat signals are signals x(n) with period N anddefinedas
x(n +
N
2
)=−x(n)for n∈
0,
N
2
(6.6)
Inverse-repeat signals are such that their even harmonic frequency components are
all equal to zero. That is, they only have non-zero values in the odd-order harmonic
frequencies ±f
0
,±3f
0
,±5f
0
, ... with f
0
=
2π
N
. To see this, we just need to calculate
95
the corresponding N-point discrete Fourier transform (DFT). The DFT of signal
x,denoteas X, can be calculated as
X(k)=
N−1
n=0
x(n)e
j2πnk
N
(6.7)
=
N
2
−1
n=0
x(n)e
j2πnk
N
+
N
2
−1
n=0
x(n +
N
2
)e
j2πk
N
(n+
N
2
)
(6.8)
=
N
2
−1
n=0
x(n)e
j2πnk
N
−
N
2
−1
n=0
x(n)(−1)
n
e
j2πnk
N
(6.9)
=
⎧ ⎪⎪⎪⎨ ⎪⎪⎪⎩ 2
N
2
−1
n=0
x(n)e
j2πnk
N
for k odd
0for k even
(6.10)
6.2.2 Property of a Set of Inverse-Repeat Signals
Now, we investigate the non-overlapping frequency components property of a set of
inverse-repeat signals. To the best of our knowledge, this property has never been
used for input sequence design in the literature.
Considering a signal with period N, when we apply 2N-point DFT, the fre-
quency domain representation can be obtained by inserting a zero in between each
original N-point DFT. This can be written as follows:
⎧ ⎪⎨ ⎪⎩ X
2N
(2k)= X
N
(k)for k=0, ..., N− 1
X
2N
(2k+1)=0 otherwise
(6.11)
where X
N
and X
2N
denote the N-point and 2N-point DFT, respectively. So
for an inverse-repeat signal with period N, when we apply the 2N-point DFT, it
will have non-zero values only in frequency components (±f
0
,±3f
0
,±5f
0
, ...)×2=
±2f
0
,±6f
0
,±10f
0
, ... Following this rule, for an inverse-repeat signal with period
96
Table 6.1: Non-zero DFT indexes fro inverse-repeat signals.
2
α
N-point DFT for Inverse-Repeat Signal
Period Non-zero DFT index
2
α
N ±1,±3,±5,±7, ...
2
α−1
N ±2,±6,±10,±14, ...
2
α−2
N ±4,±12,±20,±28, ...
.
.
.
.
.
.
N 2
α
×{±1,±3,±5,±7, ...}
{} ... 5 , 4 , 3 , 2 , 1
{}{ } ... 7 , 5 , 3 , 1 1 2 = + l
{} l 2
{}{ } ... 14 , 10 , 6 , 2 2 4 = + l
{} l 4
{}{ } ... 20 , 12 , 4 4 8 = + l
{} l 8
Non-zero index for first signal
Non-zero index for second signal
Non-zero index for third signal
{} ... 5 , 4 , 3 , 2 , 1
{}{ } ... 7 , 5 , 3 , 1 1 2 = + l
{} l 2
{}{ } ... 14 , 10 , 6 , 2 2 4 = + l
{} l 4
{}{ } ... 20 , 12 , 4 4 8 = + l
{} l 8
Non-zero index for first signal
Non-zero index for second signal
Non-zero index for third signal
......
Figure 6.1: Illustration of non-zero frequency indexes for inputs with different pe-
riods.
N, if we apply 2
α
N-point DFT, it will have non-zero values only at frequencies
±2
α
f
0
,±2
α
3f
0
,±2
α
5f
0
, ... Now suppose we have an inverse-repeat signal set: x
1
,
x
2
, ..., x
α+1
with different periods N
1
= N, N
2
=2N, N
3
=4N, ..., N
α+1
=2
α
N,
where the number of the set is α + 1. When we apply DFT over 2
α
N to all of the
signals, they will possess non-zero value for some frequency components as shown
in Table 6.1:
Fig. 6.1 shows the non-zero DFT indexes, which obviously form a disjoint set.
So we find that signals with different periods occupy non-overlapping frequency
components in the frequency domain. This means that if we select signals from the
set as inputs to a LTI system, any of two signals will have circular cross-correlation
97
equal to zero with period 2
α
N with any shifts. This design procedure will be
illustrated with an example in the following section.
6.2.3 Design Example
Using the procedure, we can construct a set of discrete α+1 inverse repeat signals,
whose periods are N,2N,4N,8N,..., 2
α
N. To apply this set of signals in practice,
a natural way is to choose the periods as small as possible. The smallest period
is N = 2 (for discrete inverse-repeat signals, the smallest period is 2) and for this
signal within one period, it has the a positive value and a negative value, with the
same amplitude (normally chosen as
1
√
2
to make the energy normalized to 1). For
the signal with period 2N = 4, it is natural to choose two positive values followed
by two negative values with the same amplitude. We can use the same procedure
and extend it to the signal with period 2
α
N =2
α+1
, which has 2
α
positive values
followed by 2
α
negative values, all with the same amplitude. Fig. 6.2 shows an
example for one period of this set of signals with in the α = 3 case. As mentioned
before, if we choose any two signals from this set, they will have zero circular cross-
correlation with period 2
α+1
for arbitrary shifts. This is a very nice property because
that means if we set these signals into a multiple-inputs LTI system, each signal
will not cause interference to other signals (because of vanishing cross-correlation
with period 2
α+1
). As a result, this set of signals can be used for some applications
if the goal is to separate the influence from different inputs. We will apply this
inverse-repeat signal design to some reservoir models and compare the results with
those achieved with some well-known deterministic sequences in Section 6.3.5.
In summary, this new procedure is formulated as follows:
98
Figure 6.2: This figure shows an example for designed set of inverse-repeat signals.
α = 3 (period is equal to 2
α+1
= 16) and one period is displayed.
1. Deciding the whole period of the input signal set, which is 2
α
N where N is
an integer with N ≥2and α + 1 is the number of signals in this set.
2. Choose an arbitrary signal with length N/2 that meets the constraints on
input signals (e.g., finite-alphabet), and transform it into an inverse-repeat
signal by adding the inverse version of itself. Then repeat the signal until its
length becomes 2
α
N
3. Repeat step 2 with signal lengths N,2N, ..., to 2
α−1
N. This will generate
α + 1 inverse-repeat signals with vanishing cross-correlation with each other
at arbitrary shifts.
6.2.4 Discussion
This proposed procedure can be used to design a set of signals with perfect cross-
correlation property. As a well-known result, chosen as inputs to a dynamic system
should be designed by the following rules:
1. Low out-of-phase auto-correlation value (ideal value is zero)
2. Low cross-correlation value (ideal value is zero)
99
As shown in the CDMA literature, in order to achieve the first goal, the maximum-
length-sequence (MLS) or m-sequence was introduced and developed. In order to
solve the second goal, Welch’s bound [45] was derived to describe the upper bound
on performance. Different design techniques, such as Gold pair sequences, Kasami
sequence, were proposed to approach this bound [45]. Our proposed procedure try
to addresses the problem by optimizing the second design rule, that is, it aims to
make the cross-correlation values all equal to zero. In this sense, we reach the ideal
valueforthesecond design rule.
The problem with our proposed method is that is has relatively high out-of-phase
auto-correlation values. This can be shown in either time or frequency domain. In
time domain, the signals in this set have different periods, from N to 2
α
N.When
we calculate the circular auto-correlation over 2
α
N points, the signal of period N
will perform the worst among this set of signals. Its out-of-phase auto-correlation
value will reach the maximum value for every shift by N. Ofcoursethisdoesnot
satisfy the first goal. In frequency domain, the frequency components with non-zero
values for signal with period N are only on 2
α
×{±1,±3,±5,±7, ...}.Asweknow,
a “good” probing signal should occupy (have non-zero value) as many frequency
components as possible in the frequency domain, so the signal with period N is far
away from being “good”.
The input sequence design problem based on a deterministic approach often has
some constraints on the input signals, such that the signals can only be taken in
some ranges (usually positive values with some upper and lower bounds). For some
industrial applications, the input signals may be even more constrained, e.g., they
can only be taken from some discrete values. For all of these cases, it is impossible
to design a set of input signals that reaches both the ideal value of both design
rules. What makes the design even worse is that the only way to find the optimal
100
input signals is via exhaustive search, except for some special cases. Sequences with
“good”properties according to both rules, such as Gold pair or Kasami sequence, are
only optimal for some specific applications, and not always optimal in a more general
situation such as parameter estimation of a dynamic system. Our work can be seen
as an example of designing a set of“good”input signals with perfect performance in
the second design rule, but with suboptimal performance with respect to the first
rule. Also, there may be some applications where the cross-correlation property is
much more important than the auto-correlation property, and our design can be
used to get a good performance in that kinds of applications.
Back to the goal of this chapter, to solve the optimal injection rates problem,
we will turn to the stochastic approach in the next section.
6.3 Stochastic Approach
Because of the drawbacks of the deterministic approach, we now consider the
stochastic approach. As discussed, this framework is used for most optimal in-
put design work in the system identification literature.
From the system identification literature [41], there are two design criteria for
input design. The first consideration is the input signal spectrum. [41] proposed a
criterion for the single-input single-output (SISO) case, and we extend it to multiple-
inputs single-output (MISO) and multiple-inputs multiple-outputs (MIMO) case.
Another design criterion, crest factor, quantifies the amount of input energy into
the system. Based on these two criteria, there is a design procedure to come to
achieve a compromise between them. The main novelty of this Section is we apply
this procedure to field applications with evaluation. An example for estimating
parameters in capacitance model is shown and some general discussions are made.
101
6.3.1 Criterion for Injection Rate Design
For a LTI system, in the system identification literature a complete model is used
to describe the LTI system, which we are going to use for the rest of the discussion.
Consider a discrete-time LTI system, denoted as G(z)inthe z-domain, with inputs
u(t), outputs y(t), and additive disturbance with spectrum σ
2
|H(e
jω
)|
2
, the model
in the z-domain is expressed as
y(t)= G(z)u(t)+ H(z)e(t). (6.12)
A particular model thus corresponds to specification of the functions G, H and the
probability density function (PDF) of e(t). In practice, it is very common to assume
that e(t) is Gaussian, in which case the PDF is entirely specified by its first and
second moments. The specification of (6.12) in terms of a finite number of numerical
coefficients is the most important for the purposes of system identification. That
is, the coefficients in question in model (6.12) will the same as parameters to be
determined. So we can denote the target parameters by the vector θ and have a
model description
y(t)= G(z,θ)u(t)+ H(z,θ)e(t), (6.13)
The θ ranges over a subset of R
d
,where d is the dimension of θ:
θ∈ D
M
⊂ R
d
(6.14)
This is the parameterized LTI model for any modeling approaches.
102
In [41], for prediction-error methods (PEMs), the covariance matrix of the pa-
rameter estimation P
θ
can be expressed as
P
θ
∝ σ
2
E
dˆ y(t|θ)
dθ
dˆ y(t|θ)
dθ
T
−1
, (6.15)
where ˆ y(t|θ) is the predicted value of the outputs based on the parameters. The
expression (6.15) gives a suggestive hint for the choice of input signals: choose the
outputs y(t) and corresponding inputs u(t) so that the predicted output becomes
sensitive with respect to target parameters of our interests.
In frequency domain, the asymptotic covariance matrix P
θ
is given by the inverse
of average information matrix per sample,
¯
M,as
¯
M ∝ σ
2
π
−π
G
θ
(e
jω
,θ
0
)
G
θ
(e
−jω
,θ
0
)
T
Φ
u
(ω)
Φ
v
(ω)
dω
+ σ
2
π
−π
H
θ
(e
jω
,θ
0
)
H
θ
(e
−jω
,θ
0
)
T
σ
2
Φ
v
(ω)
dω (6.16)
provided inputs u and e
0
are independent. Here G
θ
and H
θ
are the d× 1 gradients
of G and H. Introducing
˜
M(ω)=
σ
2
G
θ
(e
jω
,θ
0
)[G
θ
(e
−jω
,θ
0
)]
T
Φ
v
(ω)
(6.17)
M
e
= σ
4
π
−π
H
θ
(e
jω
,θ
0
)[H
θ
(e
−jω
,θ
0
)]
T
Φ
v
(ω)
dω (6.18)
We have
¯
M(Φ
u
)=
π
−π
˜
M(ω)Φ
u
(ω)dω + M
e
(6.19)
This expression helps us understand the influence of input spectrum to the infor-
mation matrix: to achieve a large information matrix (small covariance matrix),
103
the input power should be spent at frequencies where the weight
˜
M(ω) is large,
that is, where the Bode plot is sensitive to parameter variations. This expression,
together with the intuitions described, help us in designing the desired input signal
spectrum. Now we extend the criterion in (6.16) to MISO and MIMO cases.
6.3.1.1 MISO case
The goal now is to extend the design criterion (6.15) and (6.16) to MISO case.
Suppose now we have M inputs, for the LTI system model, the output signal can
be expressed as
y(t)=
M
i=1
G
i
(z,θ)u
i
(t)+ H(z,θ)e(t) (6.20)
Because the output number is the same as in the SISO case (only one output), we
can use the same formula in time domain as in the SISO case:
P
θ
∝ σ
2
E
dˆ y(t|θ)
dθ
dˆ y(t|θ)
dθ
T
−1
(6.21)
In the frequency domain, the expression (6.16) needs to be slightly changed:
¯
M ∝ σ
2
π
−π
G
θ
(e
jω
,θ
0
)Φ
u
(ω)[G
θ
(e
−jω
,θ
0
)]
H
Φ
v
(ω)
dω
+ σ
2
π
−π
H
θ
(e
jω
,θ
0
)
H
θ
(e
−jω
,θ
0
)
H σ
2
Φ
v
(ω)
dω (6.22)
where now G
θ
became a d× M matrix and Φ
u
became a M × M matrix, which
are expressed as
G
θ
(e
jω
)=
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ dG
1
(e
jω
,θ)
dθ
1
...
dG
M
(e
jω
,θ)
dθ
1
.
.
.
.
.
.
.
.
.
dG
1
(e
jω
,θ)
dθ
d
...
dG
M
(e
jω
,θ)
dθ
d
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (6.23)
104
and
Φ
u
(ω)=
⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Φ
11
(ω) ... Φ
1M
(ω)
.
.
.
.
.
.
.
.
.
Φ
M1
(ω) ... Φ
MM
(ω)
⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , (6.24)
where Φ
ij
represents the cross-spectral density (CSD) function between input i and
j. Therefore the matrix Φ
u
is the CSD matrix over inputs 1, 2, ..., M.
The expression in the frequency domain gives us a key intuition: if we want
some measurement (such as Tr
¯
M
) of the information matrix to be large, we
should make the matrix Φ
u
as diagonal as possible, that is, all cross-spectral density
function Φ
ij
should be as small as possible (ideally zero) for i = j,which means
all inputs should be uncorrelated to each other. This perfectly matches the second
design rule described in the deterministic approach section.
6.3.1.2 MIMO case
In MIMO case, the number of outputs no longer one, so we need to derive the
expression for both time and frequency domain. For LTI model, suppose we have
N outputs, the output j can be expressed as
y
j
(t)=
M
i=1
G
ij
(z,θ)u
ij
(t)+ H
j
(z,θ)e
j
(t) (6.25)
for j =1, 2, ..., N. In this expression, we make an assumption that each output is
statistically independent of other outputs. For real physical systems, the outputs
often interfere with each other, and this is also true for the system we are trying
to characterize. The producers in the reservoir will affect each other because the
production rates of each producer will affect the bottom-hole flowing pressure on
105
its well. In our work, however, we still make the independence assumption, which
will be verified via simulation data.
In this model and under the PEMs approach, the criterion (6.15) in the time
domain becomes
P
θ
∝ σ
2
E
N
j=1
dˆ y
j
(t|θ)
dθ
dˆ y
j
(t|θ)
dθ
T
!
−1
. (6.26)
This is based on the independence of outputs. In the frequency domain, the criterion
(6.16) becomes
¯
M ∝
N
j=1
σ
2
π
−π
G
j
(e
jω
,θ
0
)Φ
u
(ω)
G
j
(e
−jω
,θ
0
)
H
Φ
v
(ω)
dω
+ σ
2
N
j=1
π
−π
H
j
(e
jω
,θ
0
)
H
j
(e
−jω
,θ
0
)
H
σ
2
Φ
v
(ω)
dω (6.27)
where now G
j
(e
jω
,θ
0
)is the same as G
θ
(e
jω
,θ
0
) but we add the subscript j to
denote for output j, and so is the gradient of disturbance vector H
j
(e
−jω
,θ
0
).
The equation (6.35) has a very intuitive interpretation: the energy of desired
input signals should be concentrated on the frequencies that are very sensitive to the
target parameters (large G
j
) and small noise energy (small Φ
v
)so thatthe target
parameters θ could have a better estimation. We will use this formula, together
with other considerations, and take the capacitance model as an example to design
a set of signals that fit our need.
6.3.1.3 The Crest Factor
Another consideration for input design is the input power. This is because the
covariance matrix is typically multiplied by a term that is inversely proportional
106
to the input power. For practical situations, the inputs are often limited by some
upper and lower values, that is, the input signals u satisfy u≤ u≤ u where u and
u are defined upper and lower bounds on the signal instantaneous values. Thus a
desired property of the waveform can be defined in terms of the crest factor C
r
.
For a zero-mean signal, it is
C
2
r
=
max
t
u
2
(t)
lim
N→∞
1
N
N
t=1
u
2
(t)
(6.28)
A good signal waveform is one that has a small factor, with a theoretic lower bound
of C
r
≥ 1, which is achieved for binary, symmetric signals.
6.3.2 A Design Procedure
The basic design rule for the input signal design is obvious now: we should design a
signal to achieve the desired input spectrum and as small a crest factor as possible
at the same time. But these properties are somewhat in conflict. A common and
easy choice [41] to achieve the desired spectrum is to pass a white Gaussian noise
through a linear filter. By choosing the filter, we can virtually design any signal
with the desired spectrum. But the problem is that: this filtered Gaussian white
noise may have a large crest factor, so input energy sent to the system remains
small. To overcome this, another common approach [41] is to simply take the sign
of the filtered signal to make it into a binary signal. This can be adjusted to any
desired binary levels. In this design, the crest factor takes the ideal value of 1,
but by taking the sign operation, the spectrum of the signal may change. In this
situation, we can however check the spectrum of the signal before using it as input
signal, to make sure it is acceptable for our need. The main novelty here is that we
107
Figure 6.3: Block diagram of a system to generate the set of inputs obtained from
the rate design procedure. Note this design comes from the stochastic approach
and the signal will be different for each realization.
apply this procedure into field applications and take the state-of-the-art modeling
approach CM as a design example to evaluate this method.
This procedure can be summarized as follows:
1. Calculate the input spectrum criterion (6.27) from the data available
2. Design a filter using the criterion and pass a Gaussian white noise through
the filter
3. Pass the filtered signal to a sign operator
This diagram of this procedure is shown in Fig. 6.3.
108
6.3.3 Design Example: Application to CM
We now apply the design procedure to the reservoir problem. The ranges of the
injection rates vary for different injectors, and we roughly use the average value
of historic data in some time period. The injection rates should also be upper
bounded, this is because if the injected water volume is too high, the pressure will
be very high and it may break some rock layers and create some fractures. This
should be avoided in any case. The injection rates can be as low as zero (shut-in).
Thus the injection rate u
i
(k) for the injector i has the following constraint:
0≤ u
i
(k)≤ u
max
for k =1, 2, ..., K, (6.29)
where K is the time period for our design and this constraint is applied to all
injectors i=1, 2, ..., M. Besides, the injection rates should not be selected far away
from the normal operating points in order to preserve the linearity assumption.
In order to design a set of injection rates according to the CM, we need to
calculate the criterion (6.27) in advance. In CM, the system transfer function g
ij
between injector i and producer j is expressed as
g
ij
(t)=
1
τ
ij
e
−t
τ
ij
(6.30)
In the frequency domain, this becomes
G
ij
(e
jω
)=
1
√
2π(1 + jτ
ij
ω)
(6.31)
109
Suppose now we are interested in the interwell connectivities λ
ij
for all injector-
producer well pairs. Thus, the target parameter vector θ is:
θ =
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Λ
1
Λ
2
.
.
.
Λ
N
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (6.32)
where Λ
i
represents the interwell connectivity related to producer i:
Λ
i
=
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ λ
1i
λ
2i
.
.
.
λ
Mi
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (6.33)
110
Now we can calculate the criterion (6.27) for the design of input spectrum. The
G
j
(e
jω
,θ
0
) in (6.27) becomes
G
j
(e
jω
,θ
0
)=
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0
.
.
.
0
1
1+jτ
1j
ω
0 ... 0
0
1
1+jτ
2j
ω
...
.
.
.
.
.
.
.
.
. 0
0 ... 0
1
1+jτ
Mj
ω
0
.
.
.
0
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (6.34)
If we design a set of sequences with zero CSD function, which is easily achieved by
generating each signal independently, the criterion (6.27) becomes
¯
M∝
π
−π
diag
1
1+τ
2
11
ω
2
Φ
1
(e
jω
), ...,
1
1+τ
2
MN
ω
2
Φ
M
(e
jω
)
Φ
v
(ω)
dω + Noise Term (6.35)
Here we do not write the noise term because we focus on the injection rate design,
which affects only at the first term of this formula. According to (6.35), we need
to design the input spectrum of injector i based on the weights
1
1+τ
2
i1
ω
2
,...,
1
1+τ
2
iN
ω
2
,
which means we need to have a rough estimation of the parameters τ
i1
,...,τ
iN
in
advance. This makes sense because unless we know some rough characteristics on
the system, we can not design a set of optimized inputs for parameters estimation
because the optimal inputs depend on the system behavior. In field application,
111
we can always have some rough estimation of the system from historical data, so it
is not unrealistic to design a set of injection rates optimized for better parameter
estimation in a reservoir.
From this interpretation, it seems that we should put the energy into as low
frequency as possible. In practice, this is not true because the performance of
stochastic approach is only guaranteed by its asymptotic property. To put it into
practice, the stochastic approach is only justified by a large data set K,and this
is usually not the case for our applications. If the daily data are available for both
injection rates and production rates, a typical practical duration for these experi-
ments is on the order of a few months. In this situation, if we put the energy of
injection rates on very low frequency components, the parameter estimation per-
formance will degrade because the injection rates will appear to be almost constant
during the time-frame chosen to observe the system.
6.3.4 Simulation Results
In our simulation, we assume the underground model is capacitance model, and
we consider the five-spot scenario with 5 injectors and four producers. The target
parameters are only the interwell connectivity weights λ
ij
between all injector-
producer well pairs, as we have described above. We use (1) a binary white sequence;
(2) sequences generated from our procedure with different cutoff frequency. The
first-order low-pass filter (LPF) is used in order to design the inputs energy at low
frequency components. The first-order LPF in the frequency domain is expressed
as
H
LPF
(ω)=
K
1+ jωT
(6.36)
112
where K is the passband gain and T is the time constant that controls the cut-off
frequency for this filter. A larger T means a lower cut-off frequency and the signal
passing it will have more energy on the low frequency.
Fig. 6.4 shows one realization of the injection rates for two different sequences:
one is a binary white sequence and another is the filtered binary sequence with the
time constant T set to 20. The operating points of the injection rates are set to 500
bbl/day and the variations are set to plus or minus 100 bbl/day. The experiment
period is set from 60 days to 300 days with a gap of 30 days; that is, from about 2
months to 10 months.
The performance is evaluated with different experiment periods. The noise are
all set to Gaussian white noise with standard deviation 30, to simulate the noisy
environment for the data gathering and other error factors. The results are shown
in Fig. 6.5. Two different sets of system parameters, one with the τ
ij
range from 1
to 3 and another range from 10 to 15, are presented to simulate two totally different
systems. For systems with small τ
ij
, the variations on injection rates cause a nearly
instantaneous and equal change at the producer; on the other hand, the large τ
ij
result in large attenuation and more time delay. So this simulation presents two
extreme cases which are also typical in the CM.
The results show that for these two cases, the filtered binary sequence obtained
using the discussed procedure always outperforms the binary white sequence. For
the first case (τ
ij
range from 1 to 3), the performance is similar, and this is because
when the time constant τ is small, the reservoir behaves like an all-pass filter, so
there is not much improvement for putting inputs energy into the low frequency
components. For the second case (τ
ij
range from 15 to 20), the filtered binary
sequence with T ≥ 5 outperforms white binary sequence noticeably. This is because
now the reservoir behaves like a low-pass filter, if we put the input energy into
113
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 1
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 2
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 3
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 4
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 5
(a) Binary white sequence
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 1
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 2
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 3
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 4
20 40 60 80 100 120 140
0
200
400
600
Time (Day)
Injection Rates (bbl/day)
Injector 5
(b) Filtered binary sequence by our procedure
Figure 6.4: A realization of designed injection rates used for performance evaluation.
Here we only show two cases: binary white sequence and sequence obtained by
proposed procedure with the filter parameter T = 20. The noise level is set to SNR
=10.46 dB.
low-frequency components, a significant gain can be achieved. Another interesting
phenomenon is for the filtered binary sequence with T = 20, which performs much
114
50 100 150 200 250 300
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Period
Relative Error of Parameter λ
White Binary Seq.
Filtered Binary Seq. (T = 1)
T = 2
T = 5
T = 10
T = 20
(a) Reservoir parameter τ
ij
s range from 1 to 3
50 100 150 200 250 300
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Period
Relative Error of Parameter λ
White Binary Seq.
Filtered Binary Seq. (T = 1)
T = 2
T = 5
T = 10
T = 20
(b) Reservoir parameter τ
ij
s range from 10 to 15
Figure 6.5: Performance plot for injection rates design with different sequences.
The candidate sequences are binary white sequence and filtered binary sequence by
our procedure with different cut-off frequencies (controlled by T). The noise level
are all set to SNR = 10.46 dB. The results are averaged on 500 realizations.
worse when the experiment period length is small. As we have described, this is
115
because the performance is only guaranteed asymptotically whereas the experiment
period is too short in this case.
Here we only show the application of injection rate design scheme to CM. Follow-
ing the same procedure, we can apply this scheme on all predictive models, accord-
ing to the models we decide to use. There are also some potential applications for
the injection rate design, such as waterflooding surveillance and monitoring [58] [5].
We will discuss these possible applications in the last chapter.
6.3.5 Comparison to Deterministic Approach
To complete this chapter, we compare all discussed procedures, including both
deterministic and stochastic approaches. For the deterministic approach, we choose
the proposed inverse-repeat signal design, together with well-known Gold pair and
Kasami sequence [16]. For stochastic approach, we show the white binary sequence
and the filtered binary sequence with T = 5. All simulation settings are the same
as in the previous section. The results are shown in Fig. 6.6.
If we only compare deterministic approaches (Inverse-repeat design, Gold pair
and Kasimi sequence), the proposed inverse-repeat design almost always outper-
forms the rest otherdeterministic sequences, especially when the data length is
short or the reservoir setting τ
ij
s are small. For both approaches, the filtered bi-
nary sequence has the best performance for all cases except when data length is
very short. This means that with careful design, the stochastic approaches can
achieve better “average” performance than the deterministic approaches.
We need to note here that for stochastic approaches, the results are the“average”
performance on 500 realizations, but for deterministic approaches, the sequences
are fixed so the performance is identical for each run. In real field applications, the
116
50 100 150 200 250 300
0
0.1
0.2
0.3
0.4
0.5
0.6
Data Length
Relative Error of Parameter λ
White Binary Seq.
Inverse−Repeat Design
Gold Pair Seq.
Kasami Seq.
Filtered Binary Seq. (T = 5)
(a) Reservoir parameter τ
ij
s range from 1 to 3
50 100 150 200 250 300
0
0.5
1
1.5
Data Length
Relative Error of Parameter λ
White Binary Seq.
Inverse−Repeat Design
Gold Pair Seq.
Kasami Seq.
Filtered Binary Seq. (T = 5)
(b) Reservoir parameter τ
ij
s range from 10 to 15
Figure 6.6: Performance plot for injection rates design with both deterministic and
stochastic approaches. The noise level are all set to SNR = 10.46 dB. For stochastic
approaches, the results are averaged on 500 realizations.
reservoir engineers need to be careful about this: the performance is only guaranteed
on the average sense, not for each realization, so when one realization of the filtered
117
binary sequence is selected, one needs to examine if any value of its out-of-phase
auto-correlation or cross-correlation is large. If it is the case, one should discard
this realization and look for another realization to use.
6.4 Conclusion
In this chapter, we investigate the problem of injection rate design for parameter es-
timation. After the literature review, we first propose a new deterministic approach.
We introduce a property for a set of inverse-repeat signals, and use this property
to design a set of injection rates having zero cross-correlation with arbitrary shifts.
For the stochastic approach, we discuss two criterions: spectrum criterion and in-
put power criterion. We extend the spectrum criterion on SISO case [41] to MIMO
case, and apply a design procedure, binary filtered sequence, to field applications.
Finally, both deterministic and stochastic approaches are evaluated and compared
on capacitance model to show the superiority of the filtered binary sequence design.
118
Chapter 7
Conclusions and Future Work
7.1 Conclusions
In Chapter 2, we proposed three new predictive models. We first show a general
LTI system, and address the limitation of CM on its impulse response shape. The
FIR model was developed to release the shape of response curve between injectors
and producers. The second model, DCM, was derived by extending the concepts
of CM to dealing with more heterogeneous scenarios. When building the model,
in order to take the producer-to-producer interaction into consideration, we pro-
posed the M-ARX model to characterize this effect between producers. Finally,
we suggested using the prediction-error method to estimate the model parameters.
This often involves some non-linear optimization procedure, but the FIR and M-
ARX model possess the linear-in-the-parameter property so the optimization can
reduce to some kind of linear regression procedure, which makes the estimation
much less-computationally intensive.
In Chapter 3, we first discussed two approaches for model validation: (1) val-
idation based on its prediction ability on fresh data set; (2) validation based on
119
interpretation of model parameters according to some reservoir characteristics. Fol-
lowing this, we verified all proposed models and compared them with the numer-
ical simulation data. We investigated the grey-box approach for model building
procedures, and the results were also evaluated completely. Finally, we defined a
practical metric based on prediction-errors and provided a comparative analysis of
all predictive models using this metric.
In Chapter 4, a unified linear modeling framework was proposed to integrate
all predictive models. It was shown that they all can be seen as some special cases
of a general linear model, and the transfer function of this linear model can also
be interpreted as some reservoir characteristics. Also, the relationships between
different models were easily shown under this framework.
In Chapter 5, we demonstrated a totally novel application for M-ARX model.
Comparing to other models, hen we use the model to predict the behaviors for
a producer shut-in, the number of parameters needed to be retrained by M-ARX
model is much fewer. More importantly, we showed that the shut-in performance
prediction can be achieved by using the“constrained producer”, which means setting
this producer to some constant rates (instead of shutting it in) during a period of
time. As compared to actual shut-in of a producer, this procedure keeps most of
the produced rates, which makes it much more practical when we want to evaluate
various“what if”scenarios, in particular those involving a potential shut-in of some
wells.
In Chapter 6, we first did a literature survey of the input design problem on
two different fields: system identification and channel estimation in communication
systems. Two frameworks, deterministic and stochastic approaches, are discussed.
For the deterministic approach, we proposed a new procedure for generating a set
of input signals with vanishing cross-correlation property. For stochastic approach,
120
we extended the results in system identification field to MIMO case, and applied a
common procedure to field applications to generate a set of “good” injection rates.
All proposed procedures were evaluated on the capacitance model.
7.2 Future Work
For future research, the emphasis will lie on combining some automatic control
knowledge to extend the current work to waterflood management and optimization.
The following recommendations are suggested for future work:
• Most of the predictive models were only verified with synthetic simulation
data, but not on real field data. The field data, especially the production
rates, are very noisy and with many measurement errors. The linearity as-
sumption of the system and the assumption of constant bottom-hole pressure
also need to be calibrated via the field data. Besides, more field trials should
be performed to decide the parameters, such as the amplitude of variations,
the cut-off frequency and suitable time resolution of the injection patterns,
for the injection rate design experiments.
• The discussion of the predictive models focused on gross fluid production
rate (oil + water) and assumed that gas is negligible, but actually only the
oil production rates that has the important economic values. There are a
number of reference discussing the water-to-oil ratio (WOR), i.e., [38] and [54]
proposed an oil fractional-flow model to separate the oil rates from the total
liquid rates for CM. We should point out here that the model is general to all
predictive models, so it can be applied to other models directly without any
changes. More simulations and evaluations are needed for this.
121
• Up to now, the injection rate design problem only focused on the parameter
estimation part, without considering waterflood management and optimiza-
tion. Actually, one of the ultimate goals for understanding the reservoir is to
help us in decision-making. Once we have a better understanding of the reser-
voir, the next step is to make some “smart” decision or changes to increase
the oil recovery. This work provides the basis for some approaches to real-
time waterflood management and optimization, and it needs to be extended
to cover these topics.
• The injection rate design concepts could be further extended to some ap-
plications of waterflood surveillance. When the injection rates are designed
to have some variations, such as piece-wise constant curves, the correlations
between injection/production rates will reveal some information about the
reservoir. The injection rates often can be scheduled and controlled auto-
matically and precisely, which means using a designed injection rate on daily
field production without disturbing the operation is feasible. Some kinds of
real-time monitoring on the waterflood conditions can be achieved by the use
of injection rate design.
• The concepts of constrained producers for performance prediction is a totally
new idea for the predictive modeling approach, and we demonstrate how to
use it to predict the future reservoir behavior when a particular producer is
shutting in, but without needing to actually shut in this producer. Up to
now, the validation was mainly based on the synthetic simulation data, and
only single producer shut-in case was discussed. Conceptually this can be
extended to multiple-producers shut-in case. This topic looks very promising
and needs more investigation and verification for further studies.
122
References
[1] M. Abbaszadeh-Dehghani and W.E. Brigham. Analysis of well-to-well tracer
flow to determine reservoir heterogeneity. In Annual California Regional Meet-
ing, San Francisco, California, USA, March 1982.
[2] M. Abbaszadeh-Dehghani and W.E. Brigham. Analysis of well-to-well tracer
flow to determine reservoir layering. Journal of Petroleum Technology, pages
2257–2270, 1984.
[3] A. Albertino and Larry W. Lake. Inferring connectivity only from well-rate
fluctuations in waterfloods. SPE Reservoir Evaluation and Engineering Jour-
nal, 6:6–16, 2003.
[4] JA Ball, I. Gohberg, and MA Kaashoek. A frequency response function for
linear, time-varying systems. Mathematics of Control, Signals, and Systems
(MCSS), 8(4):334–351, 1995.
[5] R.P. Batycky, M.R. Thiele, R.O. Baker, and S.H. Chugh. Revisiting reservoir
flood-surveillance methods using streamlines. SPE Reservoir Evaluation and
Engineering Journal, 11:387–394, 2008.
[6] M. Biguesh and A.B. Gershman. Training-based MIMO channel estimation: A
study of estimator tradeoffs and optimal training signals. IEEE Transactions
on Signal Processing, 54:884–893, 2006.
[7] S.P. Boyd and L. Vandenberghe. Convex optimization. Cambridge Univ Pr,
2004.
[8] C. Budianu and L. Tong. On maximal accuracy estimation with output power
constraints. IEEE Transactions on Signal Processing, 50:2515–2528, 2002.
[9] C. Chatfield. Time-Series Forecasting. Chapman & Hall/CRC, 2001.
[10] C. Chatfield. The analysis of time series: an introduction.CRCPrILlc,
2004.
[11] B.L. Cooley and J.H. Lee. Control-relevant experiment design for multivariable
systems described by expansions in orthonormal bases. Automatica, 37:273–
281, 2001.
123
[12] S.N. Crozier and D.D. Falconer. Least sum of squared errors (LSSE) channel
estimation. In Proceedings of Inst. Elect. Eng., volume 138, pages 371–378,
August 1991.
[13] Oliver D.S. The sensitivity of tracer concentration to nonuniform permeability
and porosity. Transport in Porous Media, 30, 1998.
[14] C. Fragouli, N. Al-Dhahir, and W. Turin. Finite-alphabet constant-amplitude
training sequence for multiple-antenna broadband transmissions. In Proceed-
ings of IEEE International Conference of Communication, volume 1, pages
6–10, New York, NY, USA, May 2002.
[15] P.H. Gentil. The use of multilinear regression models in patterned waterfloods:
Physical meaning of the regression coefficients. Master’s thesis, The University
of Texas at Austin, Austin, Texas, USA, 2005.
[16] S.W. Golomb and G. Gong. Signal Design for Good Correlation. Cambridge,
U.K.: Cambridge Univ. Press, 2005.
[17] G.C. Goodwin. Experiment design. In 6th IFAC Symposium on System Iden-
tification, Washington D.C., USA, 1982.
[18] G.C. Goodwin and R.L. Payne. Dynamic System Identification: Experiment
Design and Data Analysis, volume 136 of Mathematics in Science and Engi-
neering. Academic Press, 1977.
[19] A. Grant. Joint decoding and channel estimation for linear MIMO channels.
In IEEE Wireless Communications Networking Conference, Chicago, Illinois,
USA, September 2000.
[20] F.J. Guevara. Improved Methodologies for Stochastically Forecasting Oil Re-
covery Processes. PhD thesis, The University of Texas at Austin, Austin,
Texas, USA, 1997.
[21] B. Hassibi and B.M. Hochwald. How much training is needed in multiple-
antenna wireless links? IEEE Transactions on Information Theory, 49:951–
963, 2003.
[22] K.J. Heffer, R.J. Fox, C.A. McGill, and N.C. Koutsabeloulis. Novel techniques
show links between reservoir flow directionality, earth stress, fault structure
and geomechanical changes in mature waterfloods. In SPE Annual Technical
Conference and Exhibition, Dallas, Texas, USA, October 1995.
[23] P.S.C. Heuberger, PMJ Van den Hof, and B. Wahlberg. Modelling and iden-
tification with rational orthogonal basis functions. Springer-Verlag New York
Inc, 2005.
124
[24] R.N. Horne. Modern Well Test Analysis: A Computer-Aided Approach.
Petroway, Inc., Palo Alto, CA, 2nd edition, 1995.
[25] O. Izgec and C.S. Kabir. Establishing injector/producer connectivity before
breakthrough during fluid injection. In SPE Western Regional Meeting,Son
Jose, California, USA, March 2009.
[26] P.R. JBallin, A.G. Journel, and Aziz. K. Prediction of uncertainty in reservoir
performance forecast. Journal of Canadian Petroleum Technology, 31:52–62,
1992.
[27] C.R. Johnson, R.A. Greenkron, and E.G. Woods. Pulse-testing: A new method
for describing reservoir flow properties between wells. Journal of Petroleum
Technology, pages 1599–1604, 1966.
[28] T. Kailath. Linear systems. Prentice-Hall Englewood Cliffs, NJ, 1980.
[29] M.M. Kamal. Interference and pulse testing - a review. Journal of Petroleum
Technology, pages 2257–2270, 1983.
[30] D. Kaviani, J.L. Jensen, L.W. Lake, and M. Fahes. Estimation of interwell
connectivity in the case of fluctuating bottomhole pressures. In Abu Dhabi In-
ternational Petroleum Exhibition and Conference,AbuDhabi,UAE,November
2008.
[31] E.G. Larsson and P. Stoica. Space-Time Block Coding for Wireless Commu-
nications. Cambridge, U.K.: Cambridge Univ. Press, 2003.
[32] K.H. Lee, N. Jafroodi, A. Ortega, and I. Ershaghi. A distributed capacitance
model for performance evaluation and prediction in waterfloods from injection
and production rates fluctuations. Plan to submitt to Journal of Petroleum
Science and Engineering, 2010.
[33] K.H. Lee, A. Ortega, N. Jafroodi, and I. Ershaghi. A multivariate autoregres-
sive model for characterizing producer-producer relationships in waterfloods
from injection/production rate fluctuations. In SPE Western Regional Meet-
ing, Anaheim, California, USA, May 2010.
[34] K.H. Lee, A. Ortega, A.M. Nejad, and I. Ershaghi. A method for charac-
terization of flow units between injection-production wells using performance
data. In SPE Western Regional and Pacific Section Joint Meeting, Bakersfield,
California, USA, March 2008.
[35] K.H. Lee, A. Ortega, A.M. Nejad, N. Jafroodi, and I. Ershaghi. A novel
method for mapping fractures and high permeability channels in waterfloods
using injection and production rates. In SPE Western Regional Meeting,Son
Jose, California, USA, March 2009.
125
[36] M.J. Levin. Optimal estimation of impulse response in the presence of noise.
IRE Transactions on Circuit Theory, 7:50–56, 1960.
[37] Y. Li. Optimal training sequences for ofdm systems with multiple transmit
antennas. In IEEE GLOBECOM, volume 3, pages 1478–1482, 2000.
[38] X. Liang, B. Weber, T.F. Edgar, L.W. Lake, M. Sayarpour, and A.A. Yousef.
Optimization of oil production in a reservoir based on capacitance model of
production and injection rates. In SPE Hydrocarbon Economics and Evaluation
Symposium, Dallas, Texas, USA, April 2007.
[39] F. Liu, J.M. Mendel, and A.M. Nejad. Forecasting injector-producer relation-
ships from production and injection rates using an extended kalman filter. In
SPE Annual Technical Conference and Exhibition, Anaheim, California, USA,
September 2007.
[40] N. Liu and D.S. Oliver. Critical evaluation of the ensemble Kalman filter on
history-matching of geologic facies. SPE Reservoir Evaluation and Engineering
Journal, 8:470–477, 2005.
[41] L. Ljung. System identification: Theory for the User. Upper Saddle River,
NJ: Prentice Hall, 2nd edition, 1999.
[42] Computer Modeling Group Ltd. CMG numerical simulators.
”http://www.cmgroup.com/software/completesuite.htm”, October 2010.
[43] James A. M. Monte Carlo simulation: Its status and future. Journal of
Petroleum Technology, 49:361–373, 1997.
[44] T.L. Marzetta. BLAST training: Estimating channel characteristics for high
capacity space-time wireless. In 37th Annual Allerton Conference on Com-
munications, Control, and Computing, Monticello, Illinois, USA, September
1999.
[45] J.L. Massey and T. Mittelholzer. Welch
a
es bound and sequence sets for code-
division multiple-access systems. Sequences II: Methods in Communication,
Security and Computer Science, 1991.
[46] R.K. Mehra. Optimal input signals for parameter estimation in dynamic sys-
tems - survey and new results. IEEE Transactions on Automatic Control,
19:753–768, 1974.
[47] R.K. Mehra. Choice of input signals. In Trends and Progress in System Iden-
tification, Pergamon Press, Oxford, 1981.
126
[48] T.S. Ng, G.C. Goodwin, and Payne R.L. On maximal accuracy estimation
with output power constraints. IEEE Transactions on Automatic Control,
22:133–134, 1977.
[49] T.S. Ng, G.C. Goodwin, and Soderstrom T. Optimal experiment design for
linear systems with input-output constraints. Automatica, 13:571–577, 1977.
[50] M.N. Panda and A.K. Chopra. An integrated approach to estimate well in-
teractions. In SPE India Oil and Gas Conference and Exhibition, New Delhi,
India, February 1998.
[51] F. Pukelsheim. Optimal design of experiments. Society for Industrial Mathe-
matics, 2006.
[52] G.C. Reinsel. Elements of multivariate time series analysis. Springer Verlag,
2003.
[53] M.A. Sabet. Well Test Analysis. Gulf Publishing Company, 1991.
[54] M. Sayarpour, E. Zuluaga, C.S. Kabir, and L.W. Lake. The use of capacitance-
resistive models for rapid estimation of waterflood performance and optimiza-
tion. In SPE Annual Technical Conference and Exhibition, Anaheim, Califor-
nia, USA, September 2007.
[55] M. Sayarpour, E. Zuluaga, C.S. Kabir, and L.W. Lake. The use of capacitance-
resistive models for rapid estimation of waterflood performance and optimiza-
tion. In SPE Annual Technical Conference and Exhibition, Anaheim, Califor-
nia, USA, September 2008.
[56] A. Scaglione and A. Vosoughi. Turbo estimation of channel and symbol in
precoded MIMO systems. In ICASSP, volume 4, pages 413–416, Montreal,
PQ, Canada, May 2004.
[57] M. Stone. Cross-validity choice and assessment of statistical predictors. Jour-
nal of Royal Statistical Society, 36(B):111–147, 1974.
[58] M. Terrado, S. Yudono, and G. Thakur. Waterflooding surveillance and mon-
itoring: Putting principles into practice. SPE Reservoir Evaluation and Engi-
neering Journal, 10:552–562, 2007.
[59] Marco R. Thiele. Streamline simulation. In 6th International Forum on Reser-
voir Simulation, Schloss Fuschl, Austria, September 2001.
[60] Marco R. Thiele and Rod. P Batycky. Water injection optimization using a
streamline-based workflow. In SPE Annual Technical Conference and Exhibi-
tion, Denver, Colorado, USA, October 2003.
127
[61] H.J.A.F. Tulleken. Grey-box modelling and identification using physical
knowledge and bayesian techniques. Automatica, 29(2):285–308, 1993.
[62] D. Weber, T.F. Edgar, L.W. Lake, L. Lasdon, S. Kawas, and M. Sayarpour.
Improvements in capacitance-resistive modeling and optimization of large scale
reservoirs. In SPE Western Regional Meeting, Son Jose, California, USA,
March 2009.
[63] A.A. Yousef. Investigating Statistical Techniques to Infer Interwell Connectiv-
ity from Production and Injection Rate Fluctuations. PhD thesis, The Univer-
sity of Texas at Austin, Austin, Texas, USA, 2005.
[64] A.A. Yousef, P.H. Gentil, J.L. Jensen, and L.W. Lake. A capacitance model
to infer interwell connectivity from production and injection rate fluctuations.
SPE Reservoir Evaluation and Engineering Journal, 9:630–646, 2006.
[65] A.A. Yousef and L.W. Lake. Analysis and interpretation of interwell connectiv-
ity from production and injection rate fluctuations using a capacitance model.
In SPE Symposium on Improved Oil Recovery, Tulsa, Oklahoma, USA, April
2006.
[66] M. Zarrop. Design for Dynamic System Identification. Lecture Notes in Con-
trol and Information Sciences. Sci. 21. Springer Verlag, Berlin, 1979.
[67] D. Zhai, J. Mendel, and F. Liu. A new method for continual forecasting of
interwell connectivity in waterfloods using an extended kalman filter. In SPE
Western Regional Meeting, 2009.
128
Abstract (if available)
Abstract
Reservoir characterization is important for reservoir management and performance optimization. For waterflood optimization, traditionally several techniques have been suggested, most of which are either too time-consuming or the data needed are often unavailable. There is a new research trend to overcome these limitations by applying advanced statistical techniques on only injection and production data, which are often readily available for any waterflood operations.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Intelligent signal processing for oilfield waterflood management
PDF
Data based modeling and analysis in reservoir waterflooding
PDF
Real-time reservoir characterization and optimization during immiscible displacement processes
PDF
Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data
PDF
Aggregation and modeling using computational intelligence techniques
PDF
Integrated reservoir characterization for unconventional reservoirs using seismic, microseismic and well log data
PDF
Rate control techniques for H.264/AVC video with enhanced rate-distortion modeling
PDF
Waterflood induced formation particle transport and evolution of thief zones in unconsolidated geologic layers
PDF
A method for characterizing reservoir heterogeneity using continuous measurement of tracer data
PDF
Dependent R-D modeling for H.264/SVC bit allocation
PDF
Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter
PDF
Characterization of visual cells using generic models and natural stimuli
PDF
Robust speaker clustering under variation in data characteristics
PDF
Optimization of coupled CO₂ sequestration and enhanced oil recovery
PDF
Sparse representation models and applications to bioinformatics
PDF
Compression algorithms for distributed classification with applications to distributed speech recognition
PDF
Parameter estimation to infer injector-producer relationships in oil fields: from hybrid constrained nonlinear optimization to inference in probabilistic graphical model
PDF
Geological modeling in GIS for petroleum reservoir characterization and engineering: a 3D GIS-assisted geostatistics approach
PDF
Deep learning for subsurface characterization and forecasting
PDF
History matching for fractured reservoirs with Hough-transform-based parameterization
Asset Metadata
Creator
Lee, Kun-Han
(author)
Core Title
Investigating statistical modeling approaches for reservoir characterization in waterfloods from rates fluctuations
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Electrical Engineering
Publication Date
11/22/2010
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
interwell connectivity,OAI-PMH Harvest,preditive model,statistical modeling,waterflood
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Ortega, Antonio (
committee chair
), Ershaghi, Iraj (
committee member
), Mendel, Jerry (
committee member
)
Creator Email
copperfield543@yahoo.com.tw,kunhanle@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m3554
Unique identifier
UC1303165
Identifier
etd-Lee-4179 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-417046 (legacy record id),usctheses-m3554 (legacy record id)
Legacy Identifier
etd-Lee-4179.pdf
Dmrecord
417046
Document Type
Dissertation
Rights
Lee, Kun-Han
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
interwell connectivity
preditive model
statistical modeling
waterflood