Bayes linear uncertainty analysis for oil reservoirs based on multiscale computer experiments
Abstract and Keywords
This article discusses the results of a Bayes linear uncertainty analysis for oil reservoirs based on multiscale computer experiments. Using the Gullfaks oil and gas reservoir located in the North Sea as a case study, the article demonstrates the applicability of Bayes linear methods to address highly complex problems for which the full Bayesian analysis may be computationally intractable. A reservoir simulation model, run at two different levels of complexity, is used, and a simulator of a hydrocarbon reservoir represents properties of the reservoir on a threedimensional grid. The article also describes a general formulation for the approach to uncertainty analysis for complex physical systems given a computer model for that system. Finally, it presents the results of simulations and forecasting for the Gullfaks reservoir.
Keywords: uncertainty analysis, oil reservoirs, linear methods, simulations, multiscale computer experiments, forecasting
10.1 Introduction
Reservoir simulators are important and widely used tools for oil reservoir management. These simulators are computer implementations of highdimensional mathematical models for reservoirs, where the model inputs are physical parameters, such as the permeability and porosity of various regions of the reservoir, the extent of potential faults, aquifer strengths and so forth. The outputs of the model, for a given choice of inputs, are observable characteristics such as pressure readings, oil and gas production levels, for the various wells in the reservoir.
Usually, we are largely uncertain as to the physical state of the reservoir, and thus we are unsure about appropriate choices of the input parameters for a reservoir model. Therefore, an uncertainty analysis for the model often proceeds by first calibrating the simulator against observed production history at the wells and then using the calibrated model to forecast future well production, and act as an information tool for the efficient management of the reservoir.
In a Bayesian analysis, all of our uncertainties are incorporated into the system forecasts. In addition to the uncertainty about the input values, there are three other basic sources of uncertainty. First, although the simulator is deterministic, an evaluation of the simulator for a single choice of parameter values can take hours or days, so that the function is unknown to us, except at the subset of values which we have chosen to evaluate. Secondly, the reservoir simulator, even at the best choice of input values, is only a model for the reservoir, and we must therefore take into account the discrepancy between the model and the reservoir. Finally, the historical data which we are calibrating against is observed with error.
This problem is typical of a very wide and important class of problems each of which may be broadly described as an uncertainty analysis for a complex physical system based on a model for the system (Sacks, Welch, Mitchell, and (p. 242) Wynn 1989; Craig, Goldstein, Seheult, and Smith 1996; O’Hagan 2006). Such problems arise in almost all areas of scientific enquiry; for example climate models to study climate change, or models to explore the origins and generating principles of the universe. In all such applications, we must deal with the same four basic types of uncertainty: input uncertainty, function uncertainty, model discrepancy and observational error. A general methodology has been developed to deal with this class of problems. Our aim, in this chapter, is to provide an introduction to this methodology and to show how it may be applied for a reservoir model of realistic size and complexity. We shall therefore analyse a particular problem in reservoir description, based upon a description of our general approach to uncertainty analysis for complex models. In particular, we will highlight the value of fast approximate versions of the computer simulator for making informed prior judgements relating to the form of the full simulator. Our account is based on the use of Bayes linear methodology for simplifying the specification and analysis for complex highdimensional problems, and so this chapter also serves as an introduction to the general principles of this approach.
10.2 Preliminaries
10.2.1 Model description
The focus of our application is a simulation of a hydrocarbon reservoir provided to us by Energy SciTech Ltd. The model is a representation of the Gullfaks oil and gas reservoir located in the North Sea, and the model is based around a threedimensional grid of size 38 × 87 × 25 where each grid cell represents a cuboid region of subterranean rock within the reservoir. Each grid cell has different specified geological properties and contains varying proportions of oil, water and gas. The reservoir also features a number of wells which, during the course of the simulation, either extract fluids from or inject fluids into the reservoir. The overall purpose of the simulation is to model changes in pressure, and the flows and changes in distribution of the different fluids throughout the reservoir, thereby giving information on the pressures and production levels at each of the wells. A simple map of the reservoir is shown in Figure 10.1.
The inputs to the computer model are a collection of scalar multipliers which adjust the magnitudes of the geological properties of each grid cell uniformly across the entire reservoir. This results in four field multipliers – one each for porosity (ϕ), xpermeability (k_{x}), zpermeability (k_{z}), and critical saturation (crw). There is no multiplier for ypermeability as the (x, y) permeabilities are treated as isotropic. In addition to these four inputs, we have multipliers for aquifer permeability (A_{p}) and aquifer height (A_{h}) giving a total of six input parameters. The input parameters and their ranges are summarised in Table 10.1. (p. 243)
The outputs of the model are collections of time series of monthly values of various production quantities obtained for each well in the reservoir. The output quantities comprise monthly values of oil, water and gas production rates, oil, water and gas cumulative production totals, watercut, gasoil ratio, bottomhole pressure and tubinghead pressure. For the purposes of our analysis, we shall focus exclusively on oil production rate since this is the quantity of greatest practical interest and has corresponding historical observations. In terms of the time series aspect of the output, we shall focus on a threeyear window in the operation of the reservoir beginning at the start of the third year of production. We smooth these 36 monthly observations by taking fourmonth averages. By making these restrictions, we focus our attention on the 10 production wells which were operational throughout this period and so our outputs now consist of a collection of 10 time series with each 12 time points.
Table 10.1 The six input parameters to the hydrocarbon reservoir model.
Description 
Symbol 
Initial range 

Porosity 
ϕ 
[0.5, 1.5] 
xpermeability 
k_{x} 
[0.25, 6] 
zpermeability 
k_{z} 
[0.1, 0.75] 
Critical saturation 
crw 
[0.4, 1.6] 
Aquifer height 
A_{h} 
[50, 500] 
Aquifer permeability 
A_{p} 
[300, 3000] 
(p. 244) For the purposes of our analysis, it will be necessary to have access to a fast approximate version of the simulator. To obtain such an approximation, we coarsened the vertical gridding of the model by a factor of 10. The evaluation time for this coarse model is between 1 and 2 minutes, compared to 1.5–3 hours for the full reservoir model.
10.2.2 Uncertainty analysis for complex physical systems
We now describe a general formulation for our approach to uncertainty analysis for complex physical systems given a computer model for that system, which is appropriate for the analysis of the reservoir model. There is a collection, x^{+}, of system properties. These properties influence system behaviour, as represented by a vector of system attributes, y = (y_{h}, y_{p}), where y_{h} is a collection of historical values and y_{p} is a collection of values that we may wish to predict. We have an observation vector, z_{h}, on y_{h}. We write
(10.1)
and suppose that the observational error, e, is independent of y with E[e] = 0.
Ideally, we would like to construct a deterministic computer model, F(x) = (F_{h}(x), F_{p}(x)), embodying the laws of nature, which satisfies y = F(x^{+}). In practice, however, our actual model F usually simplifies the physics and approximates the solution of the resulting equations. Therefore, our uncertainty description must allow both for the possible differences between the physical value of x^{+} and the best choice of inputs to the simulator, and also for the discrepancy between the model outputs, evaluated at this best choice, and the true values of the system attributes, y.
Therefore, we must make explicit our assumptions relating the computer model F(x) and the physical system, y. In general, this will be problem dependent. The simplest and most common way to relate the simulator and the system is the socalled ‘Best Input Approach’. We proceed as though there exists a value x^{∗}, independent of the function F, such that the value of F ^{∗} = F(x^{∗}) summarizes all of the information that the simulator conveys about the system, in the following sense. If we define the model discrepancy as the difference between y and F^{∗}, so that
(10.2)
then our assumption is that ϵ is independent of both F and x^{∗}. (Here, and onwards, all probabilistic statements relate to the uncertainty judgements of the analyst.) For some models, this assumption will be justified as we can identify x^{∗} with the true system values x^{+}. In other cases, this should be viewed more (p. 245) as a convenient simplifying assumption which we consider to be approximately true because of an approximate identification of this type. For many problems, whether this formulation is appropriate is, itself, the question of interest; for a careful discussion of the status of the best input approach, and a more general formulation of the nature of the relationship between simulators and physical systems, see Goldstein and Rougier (2008).
Given this general framework, our overall aim is to tackle previously intractable problems arising from the uncertainties inherent in imperfect computer models of highly complex physical systems using a Bayesian formulation. This involves a specification for (i) the prior probability distribution for best input x^{∗}, (ii) a probability distribution for the computer function F, (iii) a probabilistic discrepancy measure relating F(x^{∗}) to the system y, (iv) a likelihood function relating historical data z to y. This full probabilistic description provides a formal framework to synthesise expert elicitation, historical data and a careful choice of simulator runs. From this synthesis, we aim to learn about appropriate choices for the simulator inputs and to assess, and possibly to control, the future behaviour of the system. For problems of moderate size, this approach is appropriate, practical and highly effective (Kennedy and O’Hagan, 2001; Santner, Williams, and Notz, 2003). As the scale of the problem increases, however, the full Bayes analysis becomes increasingly difficult because (i) it is difficult to give a meaningful full prior probability specification over highdimensional spaces; (ii) the computations, for learning from data (observations and computer runs), particularly in the context of choosing informative sets of input values at which to evaluate the simulator, become technically difficult and extremely computer intensive; (iii) the likelihood surface tends to be very complicated, so that full Bayes calculations may become highly nonrobust.
However, the idea of the Bayesian approach, namely capturing our expert prior judgements in stochastic form and modifying them by appropriate rules given observations, is conceptually appropriate, and indeed there is no obvious alternative. In this chapter, we therefore describe the Bayes linear approach to uncertainty analysis for complex models. The Bayes linear approach is (relatively) simple in terms of belief specification and analysis, as it is based only on mean, variance and covariance specifications. These specifications are made directly as, following de Finetti (1974, 1975), we take expectation as our primitive quantification of uncertainty. The adjusted expectation and variance for a random vector y, given random vector z, are as follows.
(10.3)
(10.4)
(If Var[z] is not invertible, then an appropriate generalized inverse is used in the above forms.)
(p. 246) For the purpose of this account, we may either view Bayes linear analysis as a simple approximation to a full Bayes analysis, or as the appropriate analysis given a partial specification based on expectation. We give more details of the rationale and practicalities of the Bayes linear approach in the appendix, and for a detailed treatment, see Goldstein and Wooff (2007).
10.2.3 Overview of the analysis
The evaluation of complex computer models, such as the hydrocarbon reservoir simulation, at a given choice of input parameters is often a highly expensive undertaking both in terms of the time and the computation required. This expense typically precludes a largescale investigation of the behaviour of the simulation with respect to its input parameters on the basis of model evaluations alone. Therefore, since the number of available model evaluations is limited by available resources there remains a substantial amount of uncertainty about the function and its behaviour, which we represent by means of an emulator (see Section 10.3.2.1).
For some problems, an approximate version of the original simulation may also be available. This coarse simulator, denoted F^{c}, can be evaluated in substantially less time and for substantially less cost, albeit with a consequent lower degree of accuracy. Since, both this coarse simulator and the original accurate simulator, F^{a}, are models of the same physical system, it is reasonable to expect that there will be strong qualitative and quantitative similarities between the two models. Therefore, with an appropriate belief framework to link the two simulators, we can use a large batch of evaluations of F ^{c} to construct a detailed emulator of the coarse simulator, which we can then use to inform our beliefs about F ^{a} and supplement the sparse collection of available full model evaluations. This is the essence of the multiscale emulation approach.
Our multiscale analysis of the hydrocarbon reservoir model proceeds in the following stages:
1. Initial model runs and screening – we perform a large batch of evaluations of F^{c}(x) and then identify which wells are most informative and therefore most important to emulate.
2. Emulation of the coarse simulator – given the large batch of evaluations of F^{c}(x), we now emulate each of the remaining outputs after the screening process.
3. Linking the coarse and accurate emulators – we use our emulators for F^{c}(x) to construct an informed prior specification for the emulators of F^{a}(x), which we then update by a small number of evaluations of F^{a}(x).
(p. 247) 4. History matching – using our updated emulators of F^{a}(x), we apply the history matching techniques of Section 10.3.4.1 to identify a set X^{∗} of possible values for the best model input x^{∗}.
5. Refocusing – we now focus on the reduced space, X^{∗}, identified by our history matching process. In addition to the previous outputs we now consider an additional time point 12 months after the end of our original time series, which is to be the object of our forecast. We then build our emulators in the reduced space for F^{c}(x) and F^{a}(x) over the original time series and the additional forecast point.
6. Forecasting – using our emulators of F^{a}(x) within the reduced region, we forecast the ‘future’ time point using the methods from Section 10.3.5.1.
10.3 Uncertainty analysis for the Gullfaks reservoir
10.3.1 Initial model runs and screening
We begin by evaluating the coarse model F^{c}(x) over a 1000point Latin hypercube design (McKay, Beckman, and Conover, 1979; Santner et al., 2003) in the input parameters. Since emulation, history matching and forecasting are computationally demanding processes, we choose to screen the collection of 120 outputs and determine an appropriate subset which will serve as the focus of our analysis. In order to identify this reduced collection, we will apply the principal variable selection methods of Cumming and Wooff (2007) to the 120 × 120 correlation matrix of the output vectors $\left\{{F}^{c}\left({x}_{i}\right)\right\},i=1,\dots ,1000$.
10.3.1.1 Methodology – Principal variables
Given a collection of outputs, y_{1:q}, with correlation matrix R, the principal variable selection procedure operates by assigning a score ${h}_{i}={{\displaystyle {\sum}_{j=1}^{q}\text{Corr}\left[{y}_{i},{y}_{j}\right]}}^{2}$ to each output y_{i}. The first principal variable is then identified as the output which maximises this score. Subsequent outputs are then selected using the partial correlation given the set of identified principal variables. This allows for the choice of additional principal variables to be made having removed any effects of those variables already selected. To calculate this partial correlation, we first partition the correlation matrix into block form
where R_{11} corresponds to the correlation matrix of the identified principal variables, R_{22} is the correlation matrix of the remaining variables, and R_{12} and (p. 248) R_{21} are the matrices of correlations between the two groups. We then determine the partial correlation matrix R_{22·1} as
The process continues until sufficient variables are chosen that the partial variance of each remaining output is small, or a sufficient proportion of the overall variability of the collection has been explained. In general, outputs with large values of h_{i} have, on average, large loadings on important principal components of the correlation matrix and thus correspond to structurally important variables.
10.3.1.2 Application and results
The outputs from the hydrocarbon reservoir model have a group structure, with groups formed by the different wells, and different time points. We intend to retain all time points at a given well to allow for a multivariate temporal treatment of the emulation. Therefore we make our reduction in the number of wells in the model output by applying a modified version of the above procedure where, rather than selecting a single output, y_{i}, at each stage, we instead select all 12 outputs corresponding to the well with highest total h_{i} score. We then continue as before, though selecting a block of 12 outputs at each stage. The results from applying this procedure to the 10 wells in the model are given in Table 10.2.
We can see from the results that there is a substantial amount of correlation among the outputs at each of the wells, as the first identified principal well accounts for 77.7% of the variation of the collection. Introducing additional wells into the collection of principal outputs only increases the amount of variation of all the outputs explained by the principal variables by a small
Table 10.2 Table of summary statistics for the selection of principal wells.
Well name 
Well h_{i} 
Cumulative % of variation 

B2 
4526.0 
77.7 
A3H 
26.5 
81.8 
B1 
19.5 
84.6 
B5 
14.1 
87.1 
B10A 
10.5 
94.7 
B7 
7.0 
95.2 
A2AH 
7.3 
99.2 
A17 
1.1 
99.7 
B4 
1.0 
99.7 
A1H 
1.1 
100.0 
(p. 249) amount. On the basis of this information one could choose to retain the first four or five principal wells and capture between 87% and 95% of the variation in the collection. For simplicity, we choose to retain four of the ten wells, namely B2, A3H, B1 and B5.
10.3.2 Representing beliefs about F using emulators
10.3.2.1 Methodology – Coarse model emulation
We express our beliefs about the uncertainty in the simulator output by constructing a stochastic belief specification for the deterministic simulator, which is often referred to as an emulator. Our emulator for component i of the coarse simulator, F^{c}(x), takes the following form:
(10.5)
In this formulation, ${\beta}_{i}^{c}=\left({\beta}_{i1}^{c},\dots ,{\beta}_{i{p}_{i}}^{c}\right)$ are unknown scalars, ${g}_{i}\left(x\right)=\left({g}_{i1}\left(x\right),\dots ,{g}_{i{p}_{i}}\left(x\right)\right)$ are known deterministic functions of x (typically monomials), and ${u}_{i}^{c}\left(x\right)$ is a stochastic residual process. The component ${g}_{i}{\left(x\right)}^{T}{\beta}_{i}^{c}$ is a regression term which expresses the global variation in ${F}_{i}^{c}$ namely that portion of the variation in ${F}_{i}^{c}\left(x\right)$ which we can resolve without having to make evaluations for ${F}_{i}^{c}$ at input choices which are near to x. The residual u^{c}(x) expresses local variation, which we take to be a weakly stationary stochastic process with constant variance.
Often, we discover that most of the global variation for output component ${F}_{i}^{c}$ is accounted for by a relatively small subset, x_{[i]} say, of the input quantities called the active variables. In such cases, we may further simplify our emulator, as
(10.6)
where ${u}_{i}^{c}\left({x}_{\left[i\right]}\right)$ is now a stationary process in the x_{[i]} only, and ${\upsilon}_{i}^{c}\left(x\right)$ is an uncorrelated ‘nugget’ term expressing all of the residual variation which is attributable to the inactive inputs. When variation in these residual terms is small, and the number of inactive inputs is large, this simplification enormously reduces the dimension of the computations that we must make, while usually having only a small impact on the accuracy of our results.
The emulator expresses prior uncertainty judgements about the function. In order to fit the emulator, we must choose the functions g_{ij}(x), specify prior uncertainties for the coefficients β^{c}_{i} and update these by carefully chosen evaluations of the simulator, and choose an appropriate form for the local variation u^{c}_{i}(x). For a full Bayesian analysis, we must make a full prior specification for each of the key uncertain quantities, $\left\{{\beta}_{i}^{c},{u}_{i}^{c}\left({x}_{\left[i\right]}\right),{\upsilon}_{i}^{c}\left(x\right)\right\}$ often choosing a (p. 250) Gaussian form. Within the Bayes linear formulation, we need only specify the mean, variance and covariance across each of the elements and at each input value x. From the prior form (10.6), we obtain the prior mean and variance of the coarse emulator as
(10.7)
(10.8)
where a priori we consider $\left\{{\beta}_{i}^{c},{u}_{i}^{c}\left({x}_{\left[i\right]}\right),{\upsilon}_{i}^{c}\left(x\right)\right\}$ as independent. There is an extensive literature on functional emulation (Sacks et al., 1989; Currin, Mitchell, Morris, and Ylvisaker, 1991; Craig, Goldstein, Rougier, and Seheult, 2001; Santner et al., 2003; O’Hagan, 2006).
As the coarse simulator is quick to evaluate, emulator choice may be made solely on the basis of a very large collection of simulator evaluations. If coarse simulator evaluations had been more costly, then we would need to rely on prior information to direct the choice of evaluations and the form of the collection ${\mathcal{G}}_{i}={\cup}_{i,j}\left\{{g}_{ij}(\cdot )\right\}$ (Craig, Goldstein, Seheult, and Smith, 1998). We may make many runs of the fast simulator, which allows us to develop a preliminary view of the form of the function, and therefore to make a preliminary choice of the function collection G_{i} and therefore to suggest an informed prior specification for the random quantities that determine the emulator for F^{a}. We treat the coarse simulator as our only source of prior information about F^{a}(x). This prior specification will be updated by careful choice of evaluations of the full simulator, supported by a diagnostic analysis, for example based on looking for systematic structure in the emulator residuals.
With such a large number of evaluations of the coarse model, the emulator (10.6) can be identified and wellestimated from the data alone. For a Bayesian treatment at this stage, our prior judgements would be dominated by the large number of model evaluations. In contrast, our prior judgements will play a central role in our emulation of F^{a}(x), as in that case the data are far more scarce.
In general, our prior beliefs about the emulator components are structured as follows. First, we must identify, for each ${F}_{i}^{c}\left(x\right)$ the collection of active inputs which describe the majority of global variation, their associated basis functions G and the coefficients β^{c}. Having such an ample data set allows for model selection and parameter estimation to be carried out independently for each component of F ^{c} and to be driven solely by the information from the model runs. The residual process ${u}_{i}^{c}\left({x}_{\left[i\right]}\right)$ is a weakly stationary process in x_{[i]} which represents the residual variation in the emulator that is not captured by our trend in the active variables. As such, residual values will be strongly correlated for neighbouring values of x_{[i]}. We therefore specify a prior covariance structure (p. 251) over values of u^{c}_{i}(x_{[i]}) which is a function of the separation of the active variables. The prior form we use is the Gaussian covariance function
(10.9)
where ${\sigma}_{{u}_{i}}^{2}$ is the point variance at any given x, θ^{c}_{i} is a correlation length parameter which controls the strength of correlation between two separated points in the input space, and  ·  is the Euclidean norm.
The nugget process υ^{c}_{i}(x) expresses all the remaining variation in the emulator attributable to the inactive inputs. The magnitude of the nugget process is often small and so is treated as uncorrelated random noise with $\text{Var}\left[{\upsilon}_{i}^{c}\left(x\right)\right]={\sigma}_{{u}_{i}}^{2}$. We consider the point variances of these two processes to be proportions of the overall residual variance of the computer model given the emulator trend, ${\sigma}_{i}^{2}$ so that ${\sigma}_{{u}_{i}}^{2}=\left(1{\delta}_{i}\right){\sigma}_{i}^{2}$ and ${\sigma}_{{u}_{i}}^{2}={\delta}_{i}{\sigma}_{i}^{2}$ for some ${\sigma}_{i}^{2}$ and some typically small value of δ_{i}.
10.3.2.2 Application and results
We have accumulated 1000 simulator runs and identified which production wells in the reservoir are of particular interest. Prior to emulation, the design was scaled so that all inputs took the range [−1, 1], and all outputs from F ^{c} were scaled by the model runs to have mean 0 and variance 1. We now describe the emulator of component i = (w, t) of the coarse simulator F^{c}, where w denotes the well, and t denotes the time associated with the ith output component of the computer model.
The first step in constructing the emulator is to identify, for each output component ${F}_{i}^{c}$ the subset of active inputs x_{[i]} which drive the majority of global variation in ${F}_{i}^{c}$. Using the large batch of coarse model runs, we make this determination via a stepwise model search using simple linear regression. We begin by fitting each F^{c}_{i} on all linear terms in x using ordinary least squares. We then perform a stepwise delete on each regression, progressively pruning away inactive inputs until we are left with a reduced collection x_{[i]} of between 3 and 5 of the original inputs. The chosen active variables for a subset of the wells of ${F}_{i}^{c}$ are presented in the third column of Table 10.3. We can see from these results that the inputs ϕ and crw are active on almost all emulators for those two wells, a pattern that continues on the remaining two wells. Clearly ϕ and crw are important in explaining the global variation of F ^{c} across the input space. Conversely, the input variable A_{h} appears to have no notable effect on model output.
The next stage in emulator construction is to choose the functions g_{ij}(x_{[i]}) for each ${F}_{i}^{c}\left(x\right)$ which form the basis of the emulator trend. Again, since we have an ample supply of computer evaluations we determine this collection by stepwise fitting. For each ${F}_{i}^{c}\left(x\right)$ we define the maximal set of basis functions to
Table 10.3 Emulation summary for wells B2 and A3H.
Well 
Time 
x_{[i]} 
No. model terms 
Coarse simulator R^{2} 
Accurate simulator ${\tilde{R}}^{2}$ 

B2 
4 
ϕ, cr w, A_{p} 
9 
0.886 
0.951 
B2 
8 
ϕ, cr w, A_{p} 
7 
0.959 
0.958 
B2 
12 
ϕ, cr w, A_{p} 
10 
0.978 
0.995 
B2 
16 
ϕ, cr w, k_{z} 
7 
0.970 
0.995 
B2 
20 
ϕ, cr w, k_{x} 
11 
0.967 
0.986 
B2 
24 
ϕ, cr w, k_{x} 
10 
0.970 
0.970 
B2 
28 
ϕ, cr w, k_{x} 
10 
0.975 
0.981 
B2 
32 
ϕ, cr w, k_{x} 
11 
0.980 
0.951 
B2 
36 
ϕ, cr w, k_{x} 
11 
0.983 
0.967 
A3H 
4 
ϕ, cr w, A_{p} 
9 
0.962 
0.824 
A3H 
8 
ϕ, cr w, k_{x} 
7 
0.937 
0.960 
A3H 
12 
ϕ, cr w, k_{z} 
10 
0.952 
0.939 
A3H 
16 
ϕ, cr w, k_{z} 
7 
0.976 
0.828 
A3H 
20 
ϕ, cr w, k_{x} 
11 
0.971 
0.993 
A3H 
24 
ϕ, cr w, k_{x} 
10 
0.964 
0.899 
A3H 
28 
ϕ, k_{z}, A_{p} 
10 
0.961 
0.450 
A3H 
32 
ϕ, cr w, k_{z} 
11 
0.968 
0.217 
A3H 
36 
ϕ, cr w, k_{x} 
11 
0.979 
0.278 
include an intercept with linear, quadratic, cubic and pairwise interaction terms in x_{[i]}. The saturated linear regression over these terms is then fitted using the coarse model runs and we again prune away any unnecessary terms via stepwise selection. For illustration, the trend and coefficients of the coarse emulator for well B1 oil production rate at time t = 28 are given in the first row of Table 10.4.
For each component ${F}_{i}^{c}$, we have now identified a subset of active inputs x_{[i]} and a collection of p_{i} basis functions g_{i} (x_{[i]}) which adequately capture the majority of the global behaviour of ${F}_{i}^{c}$. The next stage is to quantify beliefs
Table 10.4 Table of mean coefficients from the emulator of oil production rate at well B1 and time 28.
Model 
Int 
ϕ 
ϕ^{2} 
ϕ^{3} 
crw 
crw^{2} 
crw^{3} 

Coarse 
0.663 
−0.326 
−1.858 
2.313 
−0.219 
0.064 
0.079 
Accurate 
0.612 
−0.349 
−0.599 
0.811 
0.313 
−0.331 
−0.822 
Refocused coarse 
−0.149 
2.204 
−0.614 
−0.858 
−0.586 
0.386 
0.119 
Refocused accurate 
0.678 
0.402 
−0.456 
−0.098 
−0.057 
−0.055 
0.053 
ϕ × crw 
A_{p} 
kz 
w^{c}(x) 
R^{2} 
${\tilde{R}}^{2}$ 


Coarse 
0.407 
−0.008 
0.904 

Accurate 
0.072 
0.111 
0.112 
0.952 

Refocused coarse 
0.206 
0.044 
−0.037 
0.905 

Refocused accurate 
−0.045 
−0.034 
−0.045 
−0.109 
0.945 
(p. 253) about the emulator coefficients ${\beta}_{i}^{c}$. We fit our linear description in the selected active variables using ordinary least squares assuming uncorrelated errors to obtain appropriate estimates for these coefficients. The value of $E\left[{\beta}_{i}^{c}\right]$ is then taken to be the estimate ${\widehat{\beta}}_{ij}^{c}$ from the linear regression and $\text{Var}\left[{\beta}_{i}^{c}\right]$ is taken to be the estimated variance of the corresponding estimates. As we have 1000 evaluations in an approximately orthogonal design, the estimation error is negligible.
The results of the stepwise selection and model fitting are given in the first five columns of Table 10.3. We can see from the R^{2} values that the emulator trends are accounting for a very high proportion of the variation in the model output. We observe similar performance on the emulators of the remaining wells, with the exceptions of the emulators of well B5 at times t = 4 and t = 8, which could not be wellrepresented using any number or combination of basis functions. These two model outputs were therefore omitted from our subsequent analysis leaving us with a total of 34 model outputs.
The final stage is to make assessments for the values of the hyperparameters in our covariance specifications for u_{i}(x_{[i])} and υ_{i}(x). We estimate ${\sigma}_{i}^{2}$ by the estimate of variance of the emulator trend residuals, and then obtain estimates for θ_{i} and δ_{i} by applying the robust variogram methods of Cressie (1991). We then use these estimates as plugin values for ${\theta}_{i},{\delta}_{i}^{2}$ and ${\sigma}_{i}^{2}$.
For diagnostic purposes, we then performed a further 100 evaluations of F^{c}(x) at points reasonably well separated from the original design. For each of these 100 runs, we compared the actual model output, ${F}_{i}^{c}\left(x\right)$ with the predictions obtained from our coarse emulators. For all emulators, the variation of the prediction errors of the 100 new points was comparable to the residual variation of the original emulator trend, indicating that the emulators are interpolating well and are not overfitted to the original coarse model runs. Investigation of residual plots also corroborated this result.
10.3.3 Linking the coarse and accurate emulators
10.3.3.1 Methodology – Multiscale emulation
We now develop an emulator for the accurate version of the computer model F^{a}(x). We consider that F^{c}(x) is sufficiently informative for F^{a}(x) that it serves as the basis for an appropriate prior specification for this emulator. We initially restrict our emulator for component i of F^{a}(x) to share the same set of active variables and the same basis functions as its coarse counterpart ${F}_{i}^{c}\left(x\right)$. Since the coarse model F^{c}(x) is wellunderstood due to the considerable number of model evaluations, we consider the coarse emulator structure as known and fixed and use this as a structural basis for building the emulator of F^{a}(x). Thus we specify a prior accurate emulator of the form
(10.10)
(p. 254) where ${\omega}_{i}^{c}\left(x\right)={u}_{i}^{c}\left({x}_{[i]}\right)+{\upsilon}_{i}^{c}\left({x}_{i}\right),{\omega}_{i}^{a}\left(x\right)={u}_{i}^{a}\left({x}_{[i]}\right)+{\upsilon}_{i}^{a}\left({x}_{i}\right)$ and we have an identical global trend structure over the inputs albeit with different coefficients. On this accurate emulator, we also introduce some unknown multiple of the coarse emulator residuals ${\beta}_{{\omega}_{i}}^{a}{\omega}_{i}^{c}\left(x\right)$ and include a new residual process ${\omega}_{i}^{a}\left(x\right)$ which will absorb any structure of the accurate computer model that cannot be explained by our existing set of active variables and basis functions. Alternative methods for constructing such a multiscale emulator can be found in Kennedy and O’Hagan (2000) and Qian and Wu (2008).
As we have performed a large number of evaluations of F^{c}(x), over a roughly orthogonal design, our estimation error from the model fitting is negligible and so we consider the ${\beta}_{i}^{c}$ as known for each component i, and further for any x at which we have evaluated F^{c}(x), the residuals ${\omega}_{i}^{c}\left(x\right)$ are also known. Thus we incorporate the ${\omega}_{i}^{c}\left(x\right)$ into our collection of basis functions with associated coefficient ${\beta}_{{\omega}_{i}}^{a}$. Absorbing w^{c}(x) into the basis functions and ${\beta}_{{\omega}_{i}}^{a}$ into the coefficient vector ${\beta}_{i}^{a}$ we write the prior expectation and variance for the accurate simulator as
(10.11)
(10.12)
where now ${g}_{i}\left({x}_{\left[i\right]}\right)=\left({g}_{i1}\left({x}_{\left[i\right]}\right),\dots ,{g}_{i{p}_{i}}\left({x}_{\left[i\right]}\right),{\omega}_{i}^{c}\left(x\right)\right)$, and ${\beta}_{i}^{a}=\left({\beta}_{i1}^{a},\dots ,{\beta}_{i{p}_{i}}^{a},{\beta}_{{\omega}_{i}}^{a}\right)$. We also specify the expectation and variance of the residual process w^{a}(x) to be
(10.13)
(10.14)
where the covariance function between any pair of residuals on the accurate emulator has the same prior form and hyperparameter values as that used for u^{c}(x_{[i]}) in (10.9).
We now consider the prior form of ${\beta}_{i}^{a}$ in more detail. If we believe that each of the terms in the emulator trend corresponds to a particular qualitative physical effect, then we may expect that the magnitude of these effects will change differentially as we move from the coarse to the accurate simulator. This would suggest allowing the contribution of each g_{ij}(x_{[i]}) to the trend of F^{a}(x) to be changed individually. One prior form which allows for such changes is
(10.15)
where ρ_{ij} is an unknown multiplier which scales the contribution of ${\beta}_{ij}^{c}$ to ${\beta}_{ij}^{a}$ and γ_{ij} is a shift that can accommodate potential changes in location. We consider ρ_{ij} to be independent of γ_{ij}. In order to construct our prior form for ${F}_{i}^{a}\left(x\right)$, we must specify prior means, variances and covariances for ρ_{ij} and (p. 255) γ_{ij}. We develop choices appropriate to the hydrocarbon reservoir model in Section 10.3.3.2.
As our prior statements about F^{a}(x) describe our beliefs about the uncertain value of the simulator output, we can use observational data, namely the matrix F^{a}(X^{a}) of evaluations of the accurate simulator over the elements of the chosen design X^{a}, to compare our prior expectations to what actually occurs. A simple such comparison is achieved by the discrepancy ratio for ${F}_{i}^{a}\left({X}^{a}\right)$ the vector containing accurate simulator evaluations over X^{a} for the ith output component, defined as follows
(10.16)
which has prior expectation 1, and where $\text{rk}\left\{\text{Var}\left[{F}_{i}^{a}\left({X}^{a}\right)\right]\right\}$ corresponds to the rank of the matrix $\text{Var}\left[{F}_{i}^{a}\left({X}^{a}\right)\right]$. Very large values of $\text{Dr}\left({F}_{i}^{a}\left({X}^{a}\right)\right)$ may suggest a misspecification of the prior expectation or a substantial underestimation of the prior variance. Conversely, very small values of $\text{Dr}\left({F}_{i}^{a}\left({X}^{a}\right)\right)$ may suggest an overestimation of the variability of ${F}_{i}^{a}\left(x\right)$.
Given the prior emulator for ${F}_{i}^{a}\left(x\right)$ and the simulator evaluations ${F}_{i}^{a}\left({X}^{a}\right)$ we now update our prior beliefs about ${F}_{i}^{a}\left(x\right)$ by the model runs via the Bayes linear adjustment formulae (10.3) and (104). Thus we obtain an adjusted expectation and variance for ${F}_{i}^{a}\left(x\right)$ given ${F}_{i}^{a}\left({X}^{a}\right)$.
(10.17)
(10.18)
The constituent elements of this update can be derived from our prior specifications for F^{a}(x) from (10.11) and (10.12), and our belief statements made above.
10.3.3.2 Application and results
We first consider that the prior judgement that the expected values of the fine emulator coefficients are the same as those of the coarse emulator is appropriate, and so we specify expectations E[ρ_{ij}] and E[γ_{ij}] We now describe the covariance structure for the ρ_{ij} and γ_{ij} parameters. Every ρ_{ij} (and similarly γ_{ij}) is associated with a unique well w and time point t via the simulator output component ${F}_{i}^{a}\left(x\right)$. Additionally, every ρ_{ij} is associated with a unique regression basis function g_{ij}(·). Given these associations, we consider there to be two sources of correlation between the ρ_{ij} at a given well. First, for a given well w we consider there to be temporal effects correlating all (ρ_{ij}, ρ_{i}_{′ j′}) pairs to a degree (p. 256) governed by their separation in time. Secondly, we consider that there are model term effects which introduce additional correlation when both ρ_{ij} and ρ_{i}_{′ j′} are multipliers for coefficients of the same basis functions, i.e. ${g}_{ij}(\cdot )\equiv {g}_{{i}^{\prime}{j}^{\prime}}(\cdot )$.
To express this covariance structure concisely, we extend the previous notation and write ρ_{ij} as ρ_{(w,t,k)} where w and t correspond to the well and time associated with ${F}_{i}^{a}\left(x\right)$ and where k indexes the unique regression basis function associated with ρ_{ij}, namely the single element of the set of all basis functions $\mathcal{G}={\cup}_{i,j}\left\{{g}_{ij}(\cdot )\right\}$. Under this notation, for a pair of multipliers (ρ_{(w,t,k)}, ρ_{(w,t’,k’)}) then k = k if and only if both are multipliers for coefficients of the same basis function, say ϕ^{2}, albeit on different emulators. On this basis, we write the covariance function for ρ_{(w,t,k)} as
(10.19)
where ${\sigma}_{\rho 1}^{2}$ governs the contribution of the overall temporal effect to the covariance, and ${\sigma}_{\rho 1}^{2}$ controls the magnitude of the additional model term effect, ${R}_{T}\left(t,{t}^{\prime}\right)=\mathrm{exp}\left\{{\theta}_{T}{\left(t{t}^{\prime}\right)}^{2}\right\}$ is a Gaussian correlation function over time, and ${\mathbb{I}}_{p}$ is the indicator function of the proposition p. Our covariance specification for the γ_{(ω,t,k)} takes the same form as (10.19) albeit with variances ${\sigma}_{\gamma 1}^{2}$ and ${\sigma}_{\gamma 2}^{2}$.
To complete our prior specification over F^{a}(x) we assign ${\sigma}_{\rho 1}^{2}={\sigma}_{\gamma 1}^{2}=0.1$ and ${\sigma}_{\rho 2}^{2}={\sigma}_{\gamma 2}^{2}=0.1$ for all output components, which correspond to the belief that coefficients are weakly correlated with other coefficients on the same emulator, and that the model term effect has a similar contribution to the covariance as the temporal effect. We also assigned θ_{T} = 1/12^{2} to allow for a moderate amount of correlation across time.
We now evaluate a small batch of 20 runs of the accurate simulator. The runs were chosen by generating a large number of Latin hypercube designs and selecting that which would be most effective at reducing our uncertainty about β^{a} by minimising $\text{tr}\left\{\text{Var}\left[{\widehat{\beta}}^{a}\right]\right\}$ by least squares. Considering the simulator output for each well individually, since information on ${F}_{\omega}^{a}\left(x\right)=\left({F}_{\left(\omega ,4\right)}^{a}\left(x\right),{F}_{\left(\omega ,8\right)}^{a}\left(x\right),\dots ,{F}_{\left(\omega ,36\right)}^{a}\left(x\right)\right)$ for each well w is now available in the form of the model runs ${F}_{\omega}^{a}\left({X}^{a}\right)$ over the design X^{a}, we can make a diagnostic assessment of the choices made in specifying prior beliefs about ${F}_{\omega}^{a}$. In the case of our prior specifications for the multivariate emulators for each well, we obtain discrepancy ratio values of 0.86, 1.14, 0.67, and 1.07 suggesting our prior beliefs are broadly consistent with the behaviour observed in the data. For more detailed diagnostic methods for evaluating our prior and adjusted beliefs see Goldstein and Wooff (2007).
Using the prior emulator for ${F}_{\omega}^{a}\left(x\right)$ and the simulator evaluations ${F}_{\omega}^{a}\left({X}^{a}\right)$, we update our beliefs about ${F}_{\omega}^{a}\left(x\right)$ by the model runs via the Bayes linear adjustment formulae (10.17), (10.18). To assess the adequacy of fit for the (p. 257) updated accurate emulators of the outputs updated using the 20 runs of F^{a}(x), we calculate a version of the R^{2} statistic by using the residuals obtained from the adjusted emulator trend ${g}_{i}{\left({x}_{\left[i\right]}\right)}^{T}{E}_{{F}_{i}^{a}\left({X}^{a}\right)}\left[{\beta}_{i}^{a}\right]$, denoted ${\tilde{R}}^{2}$. These are given in the final column of Table 10.3. It is clear that the majority of the accurate emulators perform well and accurately represent the fine simulator, except for the emulators at well A3H and times t = 28, 32, 36, which display poor fits to the fine model due to the behaviour of F^{c}(x) at those locations being uninformative for the corresponding accurate model. For additional illustration, the coefficients for the coarse emulator and the adjusted expected coefficients of the accurate emulator of well B1 oil production rate at time t = 28 are given in the first two rows of Table 10.4. We can see from these values in this case that both emulators have good fits to the simulator despite the different coefficients.
10.3.4 History matching and calibration
10.3.4.1 Methodology – History matching via implausibility
Most models go through a series of iterations before they are judged to give an adequate representation of the physical system. This is the case for reservoir simulators, where a key stage in assessing simulator quality is termed history matching, namely identifying the set X^{∗} of possible choices of x^{∗} for the reservoir model (i.e. those choices of input geology which give a sufficiently good fit to historical observations, relative to model discrepancy and observational error). If our search reveals no possible choices for x^{∗}, this is usually taken to indicate structural problems with the underlying model, provided that we can be reasonably confident that the set X^{∗} is indeed empty. This can be difficult to determine, as the input space over which we must search may be very high dimensional, the collection of outputs over which we may need to match may be very large, and each single function evaluation may take a very long time. We now describe the method followed in Craig, Goldstein, Seheult, and Smith (1997).
History matching is based on the comparison of simulator output with historical observations. If we evaluate the simulator at a value, x, then we can judge whether x is a member of X^{∗} by comparing F(x) with data z. We do not expect an exact match, due to observational error and model discrepancy, and so we only require a match at some specified tolerance, often expressed in terms of the number of standard deviations between the function evaluation and the data. In practice, we cannot usually make a sufficient number of function evaluations to determine X^{∗} in this way. Therefore, using the emulator, we obtain, for each x, the values E[F(x)] and Var [F(x)]. We seek to rule out regions of x space for which we expect that the evaluation F(x) is likely to be a very poor match to observed z.
(p. 258) For a particular choice of x, we may assess the potential match quality, for a single output F_{i}, by evaluating
(10.20)
which we term the implausibility that F_{i}(x) would give an acceptable match to z_{i}. For given x, implausibility may be evaluated over the vector of outputs, or over selected subvectors or over a collection of individual components. In the latter case, the individual component implausibilities must be combined, for example by using
(10.21)
We may identify regions of x with large I_{M}(x) as implausible, i.e. unlikely to be good choice for x^{∗}. These values are eliminated from our set of potential history matches, X^{∗}.
If we wish to assess the potential match of a collection of q outputs F, we use a multivariate implausibility measure analogous to (10.20) given by
(10.22)
where I(x) is scaled to have expectation 1 if we set x = x^{∗}. Unlike (10.20), the calculation of I(x) from (10.22) requires the specification of the full covariance structure between all components of z and F, for any pair of x values.
For comparison, a direct Bayesian approach to model calibration is described in Kennedy and O’Hagan (2001). The Bayesian calibration approach involves placing a posterior probability distribution on the ‘true value’ of x^{∗}. This is meaningful to the extent that the notion of a true value for x^{∗} is meaningful. In such cases, we may make a direct Bayesian evaluation over the reduced region X^{∗}, based on careful sampling and emulation within this region. If our history matching has been successful, the space over which we must calibrate will have been sufficiently reduced that calibration should be tractable and effective, provided our prior specification is sufficiently careful. As a simple approximation to this calculation, we may reweight the values in this region by some function of our implausibility measure. The Bayes linear approach to prediction that we will describe in Section 10.3.5.1 does not need a calibration stage and so may be used directly following a successful history matching stage.
10.3.4.2 Application and results
We now use our updated emulator of F^{a}(x) to history match the reservoir simulator. At a given well, we consider outputs corresponding to different (p. 259) times to be temporally correlated. Thus we apply the multivariate implausibility measure (10.22) to obtain an assessment of the potential match quality of a given input x at each well. Incorporating the definitions of z and y from (10.1) and (10.2) into the implausibility formulation, we can write the implausibility function as
(10.23)
which is a function of the adjusted expectations and variances of our emulator for F^{a}(x) given the model evaluations, combined with the corresponding observational data, z, and the covariances for the observational error, e, and the model discrepancy, ϵ.
We now specify our prior expectation and variance for the observational error e and the model discrepancy ϵ. We do not have any prior knowledge of biases of the simulator or the data therefore we assign E[e] = 0 and E[ϵ] = 0. It is believed that our available well production history has an associated error of approximately ±10%, therefore we assign 2 × sd (e_{i}) = 0.1 × z_{i} for each emulator component ${F}_{i}^{a}\left(x\right)$ and we assume that there is no prior correlation between observational errors. Assessing model discrepancy is a more conceptually challenging task requiring assessment of the difference between the model evaluated at the best, but unknown, input, x^{∗}, and the true, also unknown, value of the system. For simplicity, we assign the variance of the discrepancy to be twice that of the observational error to reflect a belief that the discrepancy has a potentially important and proportional effect. In contrast to observational errors, we introduce a relatively strong temporal correlation over the ϵ_{i} = ϵ_{(w,t)} such that Corr[ϵ_{(w,t)}, ϵ_{(w,t′)})] = exp{−θ_{T}(t − t′)^{2}}, where we assign θ_{T} = 1/36^{2} to allow the correlation to persist across all 12 time points spanning the 36month period. We specify such a correlation over the model discrepancy since we believe that if the simulator is, for example, substantially underestimating the system at time t, then it will be highly likely that it will also underpredict at time t + 1.
To assess how the implausibility of input parameter choices changes with x, we construct a grid over the collection of active inputs spanning their feasible ranges and we evaluate (10.23) for each of the four selected wells, at every x point in that grid. We then have a vector of four implausibilities for every input parameter combination in the grid. To collapse these vectors into a scalar for each x, we use the maximum projection (10.21) where we maximise over the different wells to obtain a single measure I_{M}(x). This gives a conservative measure for the implausibility of a parameter choice, x, since if x is judged implausible on any one of the wells then it is deemed implausible for the collection. Thus the implausibility scores are combined in such a way that (p. 260) a particular input point x is only ever considered a potential match to the simulator if it is an acceptable match across all wells.
Thus we obtain a quantification for the match quality for a representative number of points throughout the possible input space. The domain of the implausibility measure is a fivedimensional cube and, as such, it is hard to visualize the implausibility structure within that space. To address this problem, we project this hypercube of implausibility values down to 2D spaces in every pair of active inputs using the method described in Craig et al. (1997). If we partition x such that x = (x′, x″), then we obtain a projection of $\widehat{\mathcal{I}}\left(x\right)$ onto the subspace x′ of x by calculating
which is a function only of x′.
ℐ_{M}(x) is a Mahalanobis distance over the four time points for each well. We produce the projections of the implausibility surface in Figure 10.2, colouring by appropriate ${\mathcal{X}}_{4}^{2}$ quantiles for comparison. The first plot in Figure 10.2(a) shows the approximate proportion of the implausibility space that would be excluded if we were to eliminate all points x with I_{M}(x) greater than a number of the standard deviations from the rescaled X^{2} distribution. For example, thresholding at three standard deviations, corresponding to ${\mathcal{I}}_{M}\left(x\right)\simeq 4$, would excludes approximately 90% of of input space. The subsequent plots in Figure 10.2(b) to Figure 10.2(f) are a subset of the 2D projections of the implausibility surface onto pairs of active variables. It is clear from these plots that there are regions of low implausibility corresponding to values of ϕ less than approximately 0.8 which indicates a clear region of potential matches to our reservoir history. Higher values of ϕ are much more implausible and so would be unlikely history matches. Aside from ϕ, there appears to be little obvious structure on the remaining active variables. This is reinforced by Figure 10.2(f), which is representative of all the implausibility projections in the remaining active inputs. This plot clearly shows that there is no particular set of choices for k_{x} or k_{z} that could reasonably be excluded from consideration without making very severe restrictions of the input space. Therefore, we decide to define our region of potential matches, X^{∗}, by the set {x : I_{M}(x) ≤ 4}. Closer investigation revealed that this set can be wellapproximated by the restriction that ϕ should be constrained to the subinterval [0.5, 0.79].
10.3.4.3 Reemulation of the model
Given the reduced space of input parameters, we now refocus our analysis on this subregion with a view towards our next major task – forecasting. Our intention for the forecasting stage is, for each well, to use the (p. 261)
(p. 262) last four time points in our existing series to forecast an additional time point located 12months beyond the end of our original time series at t = 48 months. Therefore, we no longer continue investigating the behaviour of well B2 since it ceases production shortly after our original threeyear emulation period.
To forecast, we require emulators for each the four historical time points as well as the additional predictive point, which we now construct over the reduced space X^{∗}. To build the emulators, we follow the same process as described in Section 10.3.2.2 and Section 10.3.3.1, requiring a batch of additional model runs. Of our previous batch of 1000 evaluations of F^{c}(x) 262 were evaluated within X^{∗} and so can be used at this stage. Similarly, of the 20 runs of F^{a}(x) a total of six remain valid. These runs will be supplemented by an additional 100 evaluations of F^{c}(x) and then by an additional 20 evaluations of F^{a}(x).
Adopting the same strategy as Section 10.3.2.2, we construct our coarse emulators from the information contained within the large pool of model evaluations, albeit with two changes to the process. Since we have already emulated these output components in the original input space (with the exception of our predictive output), we already have some structural information in the form of the x_{[i]} and the g_{i}(x_{[i]}) for each ${F}_{i}^{c}\left(x\right)$ that we obtained from the original emulation. Rather than completely reexecuting the search for active variables and basis functions, we shall begin our searches using the original x_{[i]} and g_{i}(x_{[i]}) as the starting point. We allow for the emulators to pick up any additional active variables, but not to exclude previously active inputs; and we allow for basis functions in the new x_{[i]} to be both added and deleted to refine the structure of the emulator.
An emulation summary for F^{c}(x) within the reduced region X^{∗} is given in Table 10.5 for wells A3H and B1. We can see that the emulator performance is still good with high R^{2} indicating that the emulators are still explaining a large proportion of the variation in model output. Observe that many of the emulators have picked up an additional active input variable when we refocus in the reduced input space.
Considering the emulator for F^{a}(x), we make a similar belief specification as before to link our emulator of F^{c}(x) to that of F^{a}(x). We choose to make the same choices of parameters (σ_{ρ1}, σ_{ρ2}, σ_{γ}1, σ_{γ2}, θ_{T}) to reflect a prior belief that the relationship between the two emulators in the reduced space is the similar to that over the original space. Comparing this prior with the data via the discrepancy ratio (10.16) showed that it was again reasonably consistent with Dr(F^{a}(X^{a})) taking values of 2.14, 0.98, and 2.02, although perhaps we may be slightly understating our prior variance. The prior emulator was then updated by the runs of the accurate model. Looking at the final column at Table 10.5 we see that the emulator trend fits the data well, although the ${\tilde{R}}^{2}$ values appear to be decreasing over time. Example coarse and accurate emulator coefficients for (p. 263)
Table 10.5 Refocused emulation summary for wells A3H and B1.
Well 
Time 
x_{[i]} 
No. model terms 
Coarse trend R^{2} 
Accurate trend ${\tilde{R}}^{2}$ 

A3H 
24 
ϕ, cr w, k_{x}, k_{z} 
8 
0.981 
0.974 
A3H 
28 
ϕ, cr w, k_{x}, k_{z}, A_{p} 
11 
0.971 
0.989 
A3H 
32 
ϕ, cr w, k_{x}, k_{z} 
11 
0.973 
0.958 
A3H 
36 
ϕ, cr w, k_{x} 
10 
0.958 
0.917 
A3H 
48 
ϕ, cr w, k_{x} 
10 
0.981 
0.888 
B1 
24 
ϕ, cr w, k_{x}, A_{p} 
11 
0.894 
0.982 
B1 
28 
ϕ, cr w, k_{z}, A_{p} 
9 
0.905 
0.945 
B1 
32 
ϕ, cr w, k_{x} 
11 
0.946 
0.953 
B1 
36 
ϕ, cr w, k_{x} 
11 
0.953 
0.927 
B1 
48 
ϕ, cr w, k_{x}, A_{p} 
11 
0.941 
0.880 
the refocused emulator are also given in the bottom two rows of Table 10.4, which shows more variation in the coefficients as we move from coarse to fine in the reduced space and also shows the presence of an additional active input.
10.3.5 Forecasting
10.3.5.1 Methodology – Bayes linear prediction
We wish to predict the collection y _{p} of future well production using the observed well history z_{h}. This is achieved by making an appropriate specification for the joint mean and variance of the collection (y_{p}, z_{h}), and so our prediction for y_{p} using the history z_{h} is the adjusted expectation and variance of y_{p} given z_{h}. This Bayes linear approach to forecasting is discussed extensively in Craig et al. (2001).
The Bayes linear forecast equations for y_{p} given z_{h} are given by
(10.24)
(10.25)
Given the relations (10.1) and (10.2), we can express this forecast in terms of the ‘best’ simulator run ${F}^{*}={F}^{a}\left({x}^{*}\right)$, the model discrepancy ϵ, the observed history z_{h}, and the observational error e. From (10.2) we write the expectation and variance of y as $E\left[y\right]={E}_{{F}^{a}\left({X}^{a}\right)}\left[{F}^{*}\right]+E[\in ]$ and $\text{Var}\left[y\right]={\text{Var}}_{{F}^{a}\left({X}^{a}\right)}\left[{F}^{*}\right]+\text{Var}\left[\epsilon \right]$, namely the adjusted expectation and variance of the best accurate simulator run F^{∗}, given the collection of available simulator evaluations F^{a}(X^{a}) plus the model discrepancy. For simplicity of presentation, we introduce the shorthand notation ${\mu}^{*}={E}_{{F}^{a}\left({X}^{a}\right)}\left[{F}^{*}\right],{\sum}^{*}={\text{Var}}_{{F}^{a}\left({X}^{a}\right)}\left[{F}^{*}\right],{\sum}^{\epsilon}=\text{Var}\left[\epsilon \right]$ and ${\sum}^{e}=\text{Var}\left[e\right]$ and we again use the subscripts h, p to indicate the relevant subvectors and submatrices of these quantities corresponding to the historical and predictive components. We also assume $E\left[\epsilon \right]=0$ to reflect the belief that there are no systematic (p. 264) biases in the model known a priori. The Bayes linear forecast equations are now fully expressed as follows
(10.26)
(10.27)
Given a specification for ${\sum}^{\epsilon}$ and ${\sum}^{e}$ we can assess the first and second order specifications $E\left[{F}_{i}^{a}\left(x\right)\right],\text{Cov}\left[{F}_{i}^{a}\left(x\right),{F}_{i}^{a}\left({x}^{\prime}\right)\right]$ from our emulator of F ^{a} for every $x,{x}^{\prime}\in {\mathcal{X}}^{*}$. We may therefore obtain the mean and variance of ${F}^{*}={F}^{a}\left({x}^{*}\right)$ by first conditioning on x^{∗} and then integrating out with respect to an appropriate prior specification over X^{∗} for x^{∗}. Hence ${E}_{{F}^{a}\left({X}^{a}\right)}\left[{F}^{*}\right]$ and ${\text{Var}}_{{F}^{a}\left({X}^{a}\right)}\left[{F}^{*}\right]$ are calculated to be the expectation and variance (with respect to our prior belief specification about x^{∗}) of our adjusted beliefs about ${F}^{a}\left(x\right)\text{at}\text{\hspace{0.17em}}x={x}^{*}$ given the model evaluations ${F}^{a}\left({X}^{a}\right)$. Specifically, this calculation requires the computation of the expectations, variances and covariances of all ${g}_{ij}\left({x}_{\left[i\right]}{}^{*}\right)$ and ${\omega}_{i}^{a}\left(x\right)$, which, in general, may require substantial numerical integration.
This analysis makes predictions without a preliminary calibration stage. Therefore, the approach is tractable even for large systems, as are search strategies to identify collections of simulator evaluations chosen to minimize adjusted forecast variance. The approach is likely to be effective when global variation outweighs local variation and the overall collection of global functional forms g(x) for ${F}_{h}^{a}$ and ${F}_{p}^{a}$ are similar. It does not exploit the local information relevant to the predictive quantities, as represented by the residual terms ${\omega}_{i}^{a}\left(x\right)$ in the component emulators. If some quantities that we wish to predict have substantial local variation, then we may introduce a Bayes linear calibration stage before forecasting, whilst retaining tractability (Goldstein and Rougier, 2006).
10.3.5.2 Application and results
We now apply the forecasting methodology, as described in Section 10.3.5.1, to the three wells under consideration from the hydrocarbon reservoir model. The goal of this forecasting stage is to forecast the collection of future system output, y_{p}, using the available historical observations z_{h}. For a given well in the hydrocarbon reservoir, we will consider the vector of four average oil production rates at times t = 24, 28, 32, and 36 as historical values, and the quantity to be predicted is the corresponding production rate observed 12 months later at time t = 48. As we actually have observations z_{p} on y_{p}, these may act as a check on the quality of our assessments.
By history matching the hydrocarbon reservoir, we have determined a region X^{∗} in which we believe it is feasible that an acceptable match x^{∗} should lie. For our forecast, we shall consider that x^{∗} is equally likely to be any input point contained within the region X^{∗}. This means that we take our prior for x^{∗} to (p. 265) be uniform over X^{∗}. As we take the g_{ij}(x) to be polynomials in x, then the expectations, variances and covariances of the g_{ij}(x^{∗}) can be found analytically from the moments of a multivariate uniform random variable which greatly simplifies the calculations of μ^{∗} and ∑^{∗}, the adjusted mean and variance for ${F}^{*}={F}^{a}\left({x}^{*}\right)$. We now refine our previously generous specification for Var [e]. Since each output quantity is an average of four monthly values, we now take Var [e] to be 1/4 its previous value to reflect the reduced uncertainty associated with the mean value.
Before we can obtain a prediction for y_{p} given z_{h} we require an appropriate belief specification for the model discrepancy ϵ, both at the past time points, and also at the future time point to be predicted. The role of the model discrepancy, ϵ, is important in forecasting as it quantifies the amount by which we allow the prediction to differ from our mean simulator prediction, ${\mu}^{*}=E\left[{F}^{*}\right]$, in order to move closer to the true system value y_{p}. If the specified discrepancy variance is too small, then we will obtain overconfident and potentially inaccurate forecasts located in the neighbourhood of μ^{∗}. If the discrepancy variance is too large then the forecast variances could be unfavourably large or we could overcompensate by the discrepancy and move too far from y_{p}. We now briefly consider the specification of Var [ϵ].
Consider the plots in Figure 10.3 depicting the outputs from the reservoir model over the time period of our predictions. Observe that for time points 24 to 36 at wells B1 and B5, the mean values μ^{∗} (indicated by the solid black circles) underestimate the observational data (indicated by the thick solid line). However at well A3H, the simulator ‘overestimates’ observed production. Furthermore, we observe that for well A3H the size of μ^{∗} − z is a decreasing function of time, for well B1 this distance increases over time, and for well B5 μ^{∗} − z appears to be roughly constant.
Given the specification for the observational error used in Section 10.3.4.2 and the observed history, we can compare Var [e_{h}] to the observed values of ${\left({\mu}^{*}{z}_{h}\right)}^{2}$ at each well and historical time point to obtain a simple order of magnitude data assessment for the discrepancy variance at the historical time points. To obtain our specification for Var [ϵ], we took a weighted combination of prior information in the form of our belief specification for Var [ϵ] from Section 10.3.4.2 and these order of magnitude numerical assessments. As the value of z_{p} is unknown at the prediction time t = 48, we must make a specification for $\text{Var}\left[{\epsilon}_{p}\right]$ in the absence of any sample information. To make this assessment, we performed simple curve fitting to extrapolate the historical discrepancy variances to the forecast point t = 48. The resulting specification for Var [ϵ] is given in Table 10.6; the correlation structure over the discrepancies is the same as in Section 10.3.4.2
Given all these specifications, we evaluate the forecast equations (10.26) and (10.27). The results of the forecasts are presented in Table 10.7 alongside the corresponding prior values, and the actual observed production z_{p} at t = 48. (p. 266)
Table 10.6 Table of specified values for 2 × sd ϵ_{(w,t)}.
2 × sd ϵ_{(w,t)} 
t = 24 
t = 28 
t = 32 
t = 36 
t = 48 

A3H 
504.9 
390.4 
124.5 
71.5 
14.7 
B1 
130.5 
142.0 
260.4 
245.1 
408.0 
B5 
284.3 
239.3 
305.8 
214.7 
260.8 
The forecasts and their errors are also shown on Figure 10.3 by a hollow circle with attached error bars.
We can interpret the prediction E_{zh} [y_{p}] as the forecast from the simulator for y_{p} which is modified in the light of the discrepancy between E_{zh} [y_{h}] and the observed z_{h}. In the case of the wells tabulated, the simulator for well A3H overestimates the value of z_{h} during the period t = 24, … , 36 resulting in a negative discrepancy and a consequent downward correction to our forecast. However, interestingly in the intervening period the simulator changes from underestimating to overestimating observed well production, and so on the basis of our observed history alone we underpredict the oil production rate for this well. The wells B1 and B5 behaved in a more consistent manner with a constant underprediction of the observed data compared to the observed history being reflected by an increase to our forecast. Note that whilst some of the best runs of our computer model differ substantially from z_{p}, the corrections made by the model discrepancy result in all of our forecast intervals being within the measurement error of the data.
In practice, the uncertainty analysis for a reservoir model is an iterative process based on monitoring the sizes of our final forecast intervals and investigating the sensitivity of our predictions to the magnitude of the discrepancy variance and correlation, for example by repeating the forecasts using the discrepancy variance values from Table 10.6 scaled by some constant α, for various values of α and varying the degree of temporal correlation across discrepancy terms. If we obtain forecast intervals which are sufficiently narrow to be useful to reservoir engineers then the we can end our analysis at this
Table 10.7 Forecasts of y_{p} at t = 48 using z_{h} for three wells in the hydrocarbon reservoir.
A3H 
B1 
B5 


Observation 
z_{p} 
190.0 
1598.6 
227.0 
2 × sd (e) 
19.0 
159.9 
22.7 

Prior 
E [y_{p}] 
170.2 
1027.1 
69.4 
2 × sd (y_{p}) 
62.9 
595.1 
270.3 

Forecast 
E_{zh} [y_{p}] 
170.3 
1207.1 
299.8 
2 × sd_{zh} (y_{p}) 
39.1 
348.9 
167.4 
(p. 268) stage. If, however, our prediction intervals are unhelpfully large then we return to and repeat earlier stages of the analysis. For example, introducing additional wells and other aspects of historical data into our analysis could be helpful in further reducing the size of X^{∗}, and allow us to refocus again and therefore reduce the uncertainties attached to our emulator, and so narrow the forecast interval.
Appendix
A. Broader context and background
The Bayes linear approach is similar in spirit to conventional Bayes analysis, but derives from a simpler system for prior specification and analysis, and so offers a practical methodology for analysing partially specified beliefs for large problems. The approach uses expectation rather than probability as the primitive for quantifying uncertainty; see De Finetti (1974, 1975). In the Bayes linear approach, we make direct prior specifications for that collection of means, variances and covariances which we are both willing and able to assess. Given two random vectors, B, D, the adjusted expectation for element B_{i}, given D, is the linear combination a_{0} + a^{T} D minimising E[(B_{i} − a_{0} − a^{T} D)^{2}] over choices of a_{0}, a. The adjusted expectation vector, E _{D} [B] is evaluated as
(If Var [D] is not invertible, then we use an appropriate generalised inverse). The adjusted variance matrix for B given D, is
Stone (1963), and Hartigan (1969) were among the first to discuss the role of such assessments in partial Bayes analysis. A detailed account of Bayes linear methodology is given in Goldstein and Wooff (2007), emphasizing the interpretive and diagnostic cycle of subjectivist belief analysis.
The basic approach to statistical modelling within this formalism is through second order exchangeability. An infinite sequence of vectors is secondorder exchangeable if the mean, variance and covariance specification is invariant under permutation. Such sequences satisfy the second order representation theorem which states that each element of such a sequence may be decomposed as the uncorrelated sum of an underlying ‘population mean’ quantity and an individual residual, where the residual quantities are themselves uncorrelated with zero mean and equal variances. This is similar in spirit to de Finetti’s representation theorem for fully exchangeable sequences but is sufficiently weak, in the requirements for prior specification, that it allows us to construct (p. 269) statistical models directly from simple collections of judgements over observable quantities.
Within the usual Bayesian view, adjusted expectation offers a simple, tractable approximation to conditional expectation, and adjusted variance is a strict upper bound for expected posterior variance, over all prior specifications consistent with the given moment structure. The approximations are exact in certain special cases, and in particular if the joint probability distribution of B, D is multivariate normal. Adjusted expectation is numerically equivalent to conditional expectation when D comprises the indicator functions for the elements of a partition, i.e. each D_{i} takes value one or zero and precisely one element D_{i} will equal one. We may therefore view adjusted expectation as a generalization of de Finetti’s approach to conditional expectation based on ‘calledoff’ quadratic penalties, where we remove the restriction that we may only condition on the indicator functions for a partition. Geometrically, we may view each individual random quantity as a vector, and construct the natural inner product space based on covariance. In this construction, the adjusted expectation of a random quantity Y, by a further collection of random quantities D, is the orthogonal projection of Y into the linear subspace spanned by the elements of D and the adjusted variance is the squared distance between Y and that subspace. This formalism extends naturally to handle infinite collections of expectation statements, for example those associated with a standard Bayesian analysis.
A more fundamental interpretation of the Bayes linear approach derives from the temporal sure preference principle, which says, informally, that if it is necessary that you will prefer a certain small random penalty A to C at some given future time, then you should not now have a strict preference for penalty C over A. A consequence of this principle is that you must judge now that your actual posterior expectation, E_{T} [B], at time T when you have observed D, satisfies the relation E_{T} [B] = E _{D} [B] + R, where R has, a priori, zero expectation and is uncorrelated with D. If D represents a partition, then E_{D} [B] is equal to the conditional expectation given D, and R has conditional expectation zero for each member of the partition. In this view, the correspondence between actual belief revisions and formal analysis based on partial prior specifications is entirely derived through stochastic relationships of this type.
Acknowledgements
This study was produced with the support of the Basic Technology initiative as part of the Managing Uncertainty for Complex Models project. We are grateful to Energy SciTech Limited for the use of the reservoir simulator software and providing the Gullfaks reservoir model.
References
Craig, P. S., Goldstein, M., Rougier, J. C., and Seheult, A. H. (2001), Bayesian forecasting for complex systems using computer simulators, Journal of the American Statistical Association, 96, 717–729.Find this resource:
Craig, P. S., Goldstein, M., Seheult, A. H., and Smith, J. A. (1996). Bayes linear strategies for history matching of hydrocarbon reservoirs. In Bayesian Statistics 5 (eds. J. M. Bernardo, J. O., Berger, A. P., Dawid and A. F. M., Smith), Clarendon Press, Oxford, UK, pp. 69–95.Find this resource:
Craig, P. S., Goldstein, M., Seheult, A. H., and Smith, J. A. (1997). Pressure matching for hydrocarbon reservoirs: A case study in the use of Bayes linear strategies for large computer experiments. In Case Studies in Bayesian Statistics (eds. C., Gatsonis, J. S., Hodges, R. E., Kass, R., McCulloch, P., Rossi, and N. D., Singpurwalla), SpringerVerlag, New York, vol. 3, pp. 36–93.Find this resource:
Craig, P. S., Goldstein, M., Seheult, A. H., and Smith, J. A. (1998). Constructing partial prior specifications for models of complex physical systems. Applied Statistics, 47, 37–53.Find this resource:
Cressie, N. (1991). Statistics for Spatial Data. John Wiley, New York.Find this resource:
Cumming, J. A. and Wooff, D. A. (2007). Dimension reduction via principal variables, Computational Statistics & Data Analsysis, 52, 550–565.Find this resource:
Currin, C., Mitchell, T., Morris, M., and Ylvisaker, D. (1991). Bayesian prediction of deterministic functions with applications to the design and analysis of computer experiments. Journal of the American Statistical Association, 86, 953–963.Find this resource:
De Finetti, B. (1974), Theory of Probability, Vol. 1, New York: John Wiley.Find this resource:
De Finetti, B. (1975). Theory of Probability. Vol. 2. John Wiley, New York.Find this resource:
Goldstein, M. and Rougier, J. C. (2006), Bayes linear calibrated prediction for complex systems, Journal of the American Statistical Association, 101, 1132–1143.Find this resource:
Goldstein, M. and Rougier, J. C. (2008), Reified Bayesian modelling and inference for physical Systems. Journal of Statistical Planning and Inference, 139, 1221–1239.Find this resource:
Goldstein, M. and Wooff, D. A. (2007). Bayes Linear Statistics: Theory and Methods. John Wiley, New York.Find this resource:
Hartigan, J. A. (1969). Linear Bayes methods. Journal of the Royal Statistical Society, Series B, 31, 446–454.Find this resource:
Kennedy, M. C. and O’Hagan, A. (2000). Predicting the output from a complex computer code when fast approximations are available, Biometrika, 87, 1–13.Find this resource:
Kennedy, M. C. and O’Hagan, A. (2001). Bayesian calibration of computer models. Journal of the Royal Statistical Society, Series B, 63, 425–464.Find this resource:
McKay, M. D., Beckman, R. J., and Conover, W. J. (1979). A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21, 239–245.Find this resource:
O’Hagan, A. (2006). Bayesian analysis of computer code outputs: A tutorial. Reliability Engineering and System Safety, 91, 1290–1300.Find this resource:
Qian, Z. and Wu, C. F. J. (2008). Bayesian hierarchical modeling for integrating lowaccuracy and highaccuracy experiments. Technometrics, 50, 192–204.Find this resource:
Sacks, J.,Welch,W. J., Mitchell, T. J., and Wynn, H. P. (1989). Design and analysis of computer experiments. Statistical Science, 4, 409–435.Find this resource:
Santner, T. J., Williams, B. J., and Notz, W. I. (2003). The Design and Analysis of Computer Experiments, SpringerVerlag, New York.Find this resource:
Stone, M. (1963). Robustness of nonideal decision procedures. Journal of the American Statistical Association, 58, 480–486.Find this resource: