Prediction Focused Analysis With Canonical Functional Component Analysis: An Application To A Libyan Reservoir

Author: Lewis Li. Date: May 6th 2016

Introduction

In this notebook, we will demonstrate the application of Prediction Focused Analysis on an analog to Winstershall's Concession C97-I, which is located in the N-97 field in the Western Hameimat Trough of the Sirte basin of north-central Libya. The reservoir is thus comprised of 5 distinct compartments resulting in 4 uncertain fault transmissibilities. A large aquifer is located in the lowest compartment, but the depth of this oil water contact also remains uncertain. Other uncertain reservoir parameters are relative permeabilites for the oil and water phases and oil viscosity.

We consider the situation where 5 producers and 3 injectors have already been drilled at the locations depicted in Figure 5. The field has been in production for 3500 days, and production data is available for all 5 wells. A decision needs to be made regarding the economic feasibility drilling of a 6th producer in the smallest reservoir compartment (denoted PNEW in Figure 5). Specifically, this decision will be made based on the forecasted performance of this new well over the next 4000 days. Therefore, we will seek to estimate the P10-P50-P90 forecasts of PNEW based on the first 10 years of production data from the existing 5 producers.

Loading Prior Model Response Data

A prior set of models is required for CFCA to establish a statistical relationship between the data and forecast. In this case, a set of 500 prior models are generated by applying Monte Carlo. The prior models were forward modelled using a streamline simulator (3DSL) over all 7500 days, to encompass both the 3500 days of production data denoted by d, as well as the 4000 days of forecast h required to make the decision regarding the new well.

We will next load the data. The historical and forecast data are stored as two structs with the following fields:

Name	Data Type	Dimensions	Description
time	double	$1 \times N_{t}$	A vector containing the time (days, years) at each time step
data	double	$N_{real}\times N_{t}\times N_{responses}$	Data array that contains responses. Each row is one realization, each y is a time step, and each z is a different well
name	string	1	The name of the response, ex: 'Oil Rate'
spline	double	$2$	First element specifies order of spline, second the number of knots
ObjNames	CellArray	$N_{responses} \times 1$	Contains the name of each well ex: {'P1','P2'}
types	string	1	'Historical' or 'Forecast'

For benchmarking, we will set one realization aside and use that as the truth. In this particular case, realization 150 was chosen as the reference. The production data from each of the 5 producers is shown below. The reference is shown in red $d_{obs}$, while the historical responses from the prior model are shown in gray ($d_{prior}$).

Likewise, we can plot the simulated forecast data for the well that is to be drilled for all the prior models $h_{prior}$. The reference (and unknown in reality) forecast is shown in red.

Dimension Reduction With Functional Data Analysis

Each of the response time series are technically infinite dimension, but have been discretized to 200 time steps for purposes of flow simulation. We require a dimension reduction before we can build a statistical relationship between forecast and historical data. One method such method for dimension reduction is Functional Data Analysis (ref). The basic idea is to simply express each response as a linear combination of basis functions. $$d(t) = \sum_{i=1}^K \alpha_i B_{i,n}$$ We use B-splines as basis functions in this example. Using a spline-basis has the advantage of computational ease of evaluation as well as establishing derivatives. The number of knots and order is a modeling choice and will need to be tuned for each case, usually using cross-validation. In this example, we will use a 6th order spline with 20 knots to perform FDA.

Solving for the $\alpha_i$ coefficients effectively reduces the dimension of each response curve to $K$. However, we can achieve further compression by applying Principal Component Analysis on the coefficients. This is termed Functional Princpal Component Analysis (FPCA). The effectiveness of the compression can be verified by inspecting the eigenvalues after PCA. In this case, we will look at the decomposition of Producer #1's historic oil production rate.

This shows that only the first 5 eigencomponents are required to capture over 98% of the variation in the response of all the prior models for Producer 1. The choice of spline is then iteratively adjusted to minimize the average RMS between the original. It is also important to ensure that oscillations caused by Runge’s phenomena do not occur over any of the models, when high order splines are used. We verify the effectiveness of the compression by reconstructing the original curves from the basis functions and coefficients.

Likewise, we can repeat this process for all producers, as well as the well that is to be drilled. We thus have low dimension reduction of all the required responses in the functional space $(d^f_{P1},d^f_{P2},...d^f_{P5})$ in addition to $h^f_{PNEW}$ We inspect the reconstruction of the forecasts 2 randomly selected models for the new well (PNEW).

We must then decide how many eigenvalues we should like to keep for constructing the statistical relationship. These two parameters are used to control a) the minimum number of eigenvalues required and b) the percentage of variation we would like to be kept.

Mixed PCA

We next keep only the eigenvalues that meet the specified conditions and generate a combined matrix of all the producers (d^f*_Producers) and (h^f*_PNEW)

Since the production data is comprised of multiple wells, the relationship between each of the existing wells and the new well must be quantified. However, there exists a level of redundancy between the data from each production well. For instance, a shallower water oil contact in the aquifer would cause both P3 and P4 to experience water breakthrough earlier. Consequently, a second dimension reduction is applied to the production data in the form of a PCA on the matrix

Maximizing the Linear Correlation Between Data and Forecast

So we have now effectively reduced our dimension from 5 producers and 200 element time series into 12 dimensions (8 from the producers (d^f*_{MixedProducers}) and 4 from the forecast (h^f*_PNEW). This enables the possibility of performing regression between data and forecast. However, this is contingent on a sufficiently strong linear correlation between the variables, or else the resulting estimates from regression will be inconsistent. In this example, the correlation between first functional components of data and forecast is rather poor. This is usually attributed to the presence of cross-correlations among the functional data variables.

To fully capitalize on the full multi-variate correlation between all component of h and all components of d a Canonical Correlation Analysis (CCA) is performed to transform the models into a canonical space. This is a form of partial least squares. CCA relies on a linear transformation of both and such that the components in such transformation are maximally correlated, or in terms of notation: $$d^c = d^fA^T$$ and $$h^c = h^fN^T$$ where A and B are obtained as solution of: $$max_{A,B} \frac{A\Sigma_{DH}B^T}{\sqrt{A\Sigma_{DD} A^T}\sqrt{B\Sigma_{HH} B^T}}$$ with $$D = \begin{bmatrix} d_{1}^f \\ d_{2}^f \\ \vdots \\ d_{N}^f \end{bmatrix}^T, H = \begin{bmatrix} h_{1}^f \\ h_{2}^f \\ \vdots \\ h_{N}^f \end{bmatrix}^T$$

The covariances are defined as: $$\Sigma_{DH}=cov(D,H); \Sigma_{DD}=cov(D,D);\Sigma_{HH}=cov(H,H)$$

thereby maximizing the correlations between pairwise components of $h^c_i$ and $d^c_i$ and while constraining all the inter-component correlations between $d^c_i$ and $d^c_{i\neq j}$, between $h^c_i$ and $d^h_{i\neq j}$, $h^c_i$ and $d^h_{i\neq j}$ to be 0.

To apply linear Gaussian regression between $h_c$ and $d_c$, a normal score transform is required for $h_c$.

Linear Gaussian Regression

The linear relationship between data and forecast implies we can express: $$d^c =Gh^c$$ Likewise, a Gaussian likelihood model formulated as follows: $$L(h^c) = exp(\frac{1}{2}(Gh^c-d^c_{obs})^TC^{-1}_{d^c} (Gh^c-d^c_{obs})$$ where $C^{-1}_{d^c}$ is the data covariance matrix of the canonical components. (Satija & Caers, 2015) derive a relationship between the stated covariance of the data variables. Since prior and likelihood are multi-variate Gaussian, the posterior is also Gaussian and the posterior mean and covariance can be readily estimated using classical methods (Tarantola, 2005).

With likelihood and prior multi-variate Gaussian and a linear model, the posterior distribution $f(h^c|d_{obs})$ is also multivariate Gaussian with mean $\tilde{h^c}$ and covariance model $\tilde{C^H}$

$$\tilde{h^c} = (G^TC_{d^c}^{-1}G + C_H^{-1})^{-1}(G^TC_{d^c}^{-1}d^c_{obs} + C_H^{-1}\bar{h^c})$$$$\tilde{C^H} = (G^TC_{d^c}^{-1}G + C_H^{-1})^{-1}$$

Note:

Steps 3 to 7 have been implemented in the ComputeCFCAPosterior.m function. We could just call that function.

Sampling The Posterior

Gaussian regression thus produces a multi-variate normal posterior $f(h^c|d_{obs})$, that is easily sampled from to produce forecast components conditioned to $d^c_{obs}$. We can thus generate posterior by simply sampling from this multivariate Gaussian.

We are not however, interested in posterior forecasts in the canonical space; rather we need to transform them back into the time-domain. We thus need to undo each of the transformations. We start by undoing the normal score transform. We can then observe the reduction in variane in the posterior models (red) versus prior (gray). This represents that a reduction in forecast uncertainty should be achieved.

We next need to undo the canonical transform: $$h^f_{posterior} = h^c_{posterior}B^{-1} + \overline{h^f_{prior}}$$

Finally, we need to reconstruct the time series, by undoing the FPCA and multiplying the coeffcients by the basis functions.

Note:

Step 7 has been coded as a function with the following usage:

We can now plot the posterior forecasts in the time doimain, as well as compute the updated P10,P50,P90. The posterior quantiles exhibit a significant reduction in uncertainty in comparison to the prior.

Compare with Rejection Sampling

To compare the reduction in uncertainty provided by CFCA to an idealized case, where multiple history matched models are available, 5,000 realizations were sampled from the prior distribution, and forward simulated. Rejection sampling is performed on the resulting responses to identify the realizations that match the production data using the same methodology as (Satija & Caers, 2015), resulting in 127 history matched models, from which (only after history matching), a forecast is produced. Note that this is performed in the time domain rather than in functional canonical space.

How Informative Are Our Wells?

An uncertainty statement for a forecast, can be as simple as stating a posterior PDF of that forecast. Since all reservoir modeling forecasts are obtained from limited amount of samples from that PDF, a question of confidence on the stated uncertainty is required. For example, if an uncertainty interval based on quantiles is provided (often in terms of P10-P90), then the confidence on those quantiles needs to be calculated. First, this is relevant for testing whether the posterior quantiles are different from the prior quantiles, as this would indicate predictivity of the data towards that particular forecast. Note that uncertainty and confidence intervals need not be related. One may be very confident on a wide uncertainty, and less confident on a narrow uncertainty. Since, this depends on a number of factors (the prior uncertainty, the particular data and forecast, the flow model, etc), a quantitative method is required to establish this relationship.

In this section we develop a hypothesis test, which tests whether the data informs the forecast (prediction variables h) using the above procedure. Based on the p-value of this hypothesis test, we then plot confidence vs uncertainty for several combinations of field data, to investigate what combination of historical data should be used for that forecast. It is clear that the most optimial situation occurs when we have small posterior uncertainty in combination with high confidence.

In general, the data is informative when there is a significant difference between the prior distribution $f(h)$ and the posterior $f(h|d_{obs})$. Since both these pdfs are multivariate and functional, it would be difficult to develop a simple measure. As a proxy, we propose using the first functional component $h^f_1$ containing the maximum variability of H. We therefore define a theoretical difference between prior and posterior functional component as:

$$\delta=\Delta(f(h^f_1),f(h^f_1|d_{obs})$$

$\Delta$ can be any cdf or pdf based difference such as L1 norm difference (Fenwick et al., 2014), the Kolmogorov-Smirnov difference ( (Justel, 1997); (Smirnov, 1948)), Jensen-Shannon or Kullback-Liebler divergence (Kullback & Leibler, 1951). For the purpose of hypothesis testing, the null hypothesis is that the prior and the posterior distributions of the first functional component are not different. For bootstrapping, B datasets are drawn from the existing N scoping runs without replacement, resulting in a B bootstrap estimates of the difference:

$$\hat{\hat{\delta_b}} = \Delta_{CDF}((\hat{\hat{F_b}}(h^f_1),\hat{\hat{F_b}}(h^f_1|d_{obs})), b = 1,...,N$$

The number of times $\hat{\hat{\delta_b}} < \hat{\delta}$ is a measure of how strong the predictivity is coinciding with the number of times the null-hypothesis is rejected. The more the null-hypothesis is rejected, the more confident we are on the informativeness of the data, hence a good measure is:

$$\omega = \frac{1}{B} \sum_{b=1}^B i(\hat{\hat{\delta_b}} < \hat{\delta})$$

where $i$ is the indicator function.

Bootstrap Confidence When Using All Wells

We can apply the bootstrap test to the case where all 5 wells are used. The number of bootstraps should be set to be sufficiently large. This will take some time to run.

Depending on the number of bootstrap runs, we should obtain an omega between 0.65 and 0.75, which corresponds to a fairly high confidence in our posterior uncertainty reduction.

Bootstrap Confidence Using Only P5

We next consider the case where only Producer 5 is used as the data for forecasting the new well. We first re-load the data such that only P5 is used.

From examining the prior models, we observe that most models have yet to experience water breakthrough, and thus have oil rates equal to the production limits. This means that we should not be able to extract much information about the forecast from this data; all models have the same data, but different forecasts.

We can apply the bootstrap test to the case where only P5 is used. The number of bootstraps should be set to be sufficiently large. This will again take some time to run.

This should yield an omega value of approximately 0.15, suggesting a low confidence

Combinatorial Selection of Wells

We can repeat the process by performing a combinatorial selection of wells, and running the bootstrap. We compute a $\omega$ for each set of wells, as well as the reduction in forecast uncertainty (P90-P10). We load the results and plot them below:

From this we observed that in general lower confidence correlates with larger posterior forecast uncertainty, and it can be seen that using just P4 and P5 results in smaller uncertainty reduction and poor confidence. Conversely, using all 5 wells produces both higher confidence and lower uncertainty. It should also be noted that some combinations of wells such as (P1,P2,P4) do produce reduced uncertainty, but with low confidence, and this must to be taken into account when decisions are made.

References

Ahlbrandt, T. S. (2001). The Sirte Basin Province of Libya: Sirte-Zelten total petroleum system. US Geological Survey. Denver, CO: US Department of the Interior.

Fenwick, D. S. (2014). Quantifying Asymmetric Parameter Interactions in Sensitivity Analysis: Application to Reservoir Modeling. Mathematical Geosciences , 46 (4), 493–511.

Kullback, S., & Leibler, R. (1951). On information and sufficiency. Annals of Mathematical Statistics , 22 (1), 79–86.

Ramsay, J. O. (2002). Applied functional data analysis: methods and case studies. New York: Springer.

Satija, A., & Caers, J. (2015). Direct forecasting of subsurface flow response from non-linear dynamic data by linear least-squares in canonical functional principal component space. Advances in Water Resources , 77, 69-81.

Scheidt, C. R. (2014). Prediction-Focused Subsurface Modeling: Investigating the Need for Accuracy in Flow-Based Inverse Modeling. Mathematical Geosciences , 1-19.

Tarantola, A. (2005). Inverse problem theory and methods for model parameter estimation. Philadelphia: SIAM.