Study of a Reservoir¶

Import packages¶

Reading data¶

The well information corresponds to the coordinates of their intercepts with the top and bottom surfaces of the target reservoir. It is contained in the ASCII file Wells.dat. In addition, a selection is provided which distinguishes the appraisal wells from the other ones.

The ASCII file is loaded. The file is composed of 189 wells covering a field with an extension of 10 by 15 kilometers. The depth is oriented downwards with a top varying from -983m to -1030m and the bottom from -1036m to -1067m.

The data set is represented where the well intercepts with the top surface are displayed either in black or in yellow for the appraisal wells, whereas the intercepts with the bottom surface are displayed in purple.

Obviously, all wells have penetrated the whole reservoir layer and therefore we can derive the layer thickness at each vertical well location. The thickness varies from 12.21m to 64.25m, as represented in the following histogram.

This step requires the prior definition of a grid used for mapping the surfaces whose extension and cell dimensions must be chosen carefully: it results from a trade-off between the time consumption used for mapping the surfaces and an unnecessary finesse of the grid mesh.

For this study, the optimal choice is a grid with a square mesh of 100m edge, with 101 nodes along X and 151 nodes along Y, covering a surface of 10km along X and 15km along Y. The origin of the grid (lower left corner) is located at 50m along X and 50m along Y.

The grid (only one cell out of 5 for better legibility) and the well headers are displayed.

Estimation¶

In this phase, we are facing the usual problem of estimating linked variables, such as the top, the bottom and the thickness of a layer. The principle consists in processing jointly the set of correlated variables and the uncorrelated one separately. Usually the top and the bottom are non correlated. Then we recommend estimating either the top (respectively the bottom) and the thickness separately and then simply deriving the bottom (respectively the top). Another solution would be to jointly estimate the top and the bottom using a multivariate procedure: this may become difficult if both variables are non-stationary.

Here, and for illustration sake, we will process the three variables (top, bottom and thickness) independently.

Structural analysis of the variables¶

Analysis of the top variable¶

We compute the experimental variogram of the top variable, calculated in 4 conventional directions (N0, N45, N90 and N135), for 15 lags of 500m each.

The variogram is displayed for all calculation directions

We check that the top variable is anisotropic with calculation directions overlapping two by two: 90 and 135 on one hand for the higher variability and 0 and 45 on the other hand for lower variability. As the curves are grouped two by two, this prooves that the main anisotropy axis does not coincide with an already calculated directions, but rather with an intermediate one.

Another possibility is to calculate the variogram map. We obtain the following figure:

This clearly shows the low variability direction around 20 degrees.

The directional experimental variogram is then used for fitting the model.

The model is composed of a single spherical basic anisotropic structure: the anisotropy ellipsoid is rotated (its major axis is oriented 21 degrees) with its ranges are respectively equal to 24km and 5km, as shown in the next figure.

Analysis of the bottom variable¶

We compute the experimental variogram of the thickness variables, calculated in 4 conventional directions for 15 lags of 500m each.

The corresponding variogram map is represented next.

This directional experimental variogram is used for fitting the model.

The nested model is composed of two anisotropic basic structures with (different) rotations of the anisotropy ellipsoids. The fitting is shown in the next figure.

Analysis of the thickness variable¶

We compute the experimental variogram of the thickness variables, calculated in 4 conventional directions for 15 lags of 500m each.

Next comes the calculation of the Variogram map

The directional experimental variogral is used for fitting a Model

The experimental variogram and the model are represented

Traditional Estimation methods¶

The estimation of the different variables will be performed using several interpolation methods to estimate the value of the variables on the node of the grid created beforehand. We will focus on the thickness variable only. Before kriging, we will use some traditional techniques, such as:

• the nearest neighbor

• moving average

• moving median

All these methods are based on linear techniques: for each target site, the estimated value is obtained as the linear combination of the data values at the surrounding information data points.

Note that all these method require the definition of a neighborhood (selection of data used for the estimation of each grid node). For a better understanding of the resulting map, the data information is overlaid using a proportional representation.

Moreover, if a Model is defined, we can also evaluate the variance of the estimation error for each method (in fact, we display the standard deviation instead). This requires the definition of a Model.

The nearest neighbor method¶

The principle is to set the value of each grid node at the value of the closest data point.

As mentioned, we can produce the map of the standard deviation of the estimation error.

The moving average method¶

The method assigns the average of the subset of data selected in the neighborhood centered on the target grid node. The first neighborhood consists in two five closest data selected within a circle of 5km centered on grid grid node.

The same technique is used with a neighborhood increased to the 20 closest data points. The result is an obviously smoother map.

The inverse (squared) distance¶

The method assignes the values obtained as the weighted average of the information. The weights are calculated as the inverse distance beween a datum and the target, raised to a given exponent (generally 2). The neighboring points have been limited to adistance of 2000 from the target in order to lead to reasonable time for computing the variance map.

The map of the standard deviation of estimation error.

The moving median¶

The method is quite similar to the moving average, except that the median is assigned rather than the average. As in the previous paragraph, we first define a small neighborhood composed of the 5 closest data.

The variance map

The same technique is used with a neighborhood increased to the 20 closest data points. The result is an obviously smoother map.

The variance map

Kriging Method¶

The Kriging method is a procedure which uses the model fitted above to produce the optimal estimation at each target grid node, based on a linear combination of the values measured at the samples selected in the neighborhood.

The model used for fitting the thickness variable has been fitted earlier.

Neighborhood search¶

We must now define the parameters of the neighborhood search, where the neighborhood designates the subset of samples used for performing the estimation of one target grid node.

We could imagine to use systematically all the available information when processing each grid node: this would refer to the unique neighborhood. But as we know that the dimension of the kriging system is given by the number of samples in the neighborhood, this unique neighborhood is no more relevant when the number of samples is too large (a decent limitation is around 300 samples).

We must then define a moving neighborhood for which the selection of the subset of data will move together with the target grid node, hence the name of moving neighborhood. The moving neighborhood is characterized by a set of parameters that must be defined.

In the following figures, we will check the selection of the samples in the neighborhood for the target grid node (IX=45 and IY=80), with the corresponding kriging weights (expressed in percentage).

In the first trial, we limit the extension of the neighborhood to a circle (of radius equal to 2000m) centered on the target grid node, and an optimum number of samples equal to 8. In the next figure, we can check that the 8 closest data are selected. Some samples are not selected although they lie within the circle.

Among a large number of other parameters, let us mention the possibility of ensuring that the selected samples are spread more regularly around the target by splitting the neighborhood circle into a number of angular sectors. As an example, we define a moving neighborhood composed of the 2 closest samples per angular sector within a circle of radius of 2000m, centered on the target grid node and split into 6 angular sectors.

We can check that:

• at most two samples are selected in each angular sector; in some angular sectors, the only sample present is selected,

• some samples are still not selected, although they lie within the neighborhood circle,

• the two closest data still receive large weights

• there are some negative weights.

Kriging Weights¶

We recall the Model that has been used for the Estimation of the thickness

We notice that the first lags are calculated based on few pairs. We can imagine a second fit introducing a nugget effect component and a spherical basic structure

We can use the second model for the estimation of the same target node and check for the new kriging weights.

We can check that:

• the weights attached to the closest data have decreased

• there is no more negative weight: this is one of the usual impacts of a model containing a nugget effect component

Cross-Validation¶

The cross-validation is a procedure which allows checking the consistency between the data information, the fitted model and the neighborhood. For the following test, we will use the model with no nugget effect.

The principle is to consider each datum in turn and to estimate it (by kriging) from the other samples in its neighborhood. We can then compare the true value of the target sample to its estimated value. We can also consider the standardized error calculated as the ratio between this experimental error and the standard deviation of the estimation error (or modeled error). In a Gaussian framework, only 3% of the values should lie outside the interval [-2.5;2.5]: they will be called outliers.

The cross-validation procedure results into a set of graphics such as:

• the base map pointing out the outliers,

• the histogram of the standardized errors: it should be close to a Gaussian distribution with few samples in the extreme quantiles,

• the scatter plot of the true values versus the estimated values: it should be close to the first bisector,

• the scatter plot between the standardized error and the estimated values: it should show no correlation (in the simple kriging case).

Estimation map¶

The thickness variable is estimated over the whole grid using a neighborhood slightly enlarged (5 data per angular sector within a circle of 5km radius).

The estimation is performed using the first model (without nugget effect) and presented. We clearly see the anisotropy which is responsible for the main elongated bodies throughout the grid.

Standard deviation map¶

We now display the map presenting the standard deviation of the estimation error.

Obviously, the standard deviation of the error varies from 0 at data locations up to 9.14 in extrapolated areas. The elongated aspect comes from the anisotropy of the model.

Processing the information jointly¶

As mentioned before, in order to get consistent results, it makes sense to estimate two of the three variables (from top, bottom and thickness variables) jointly. In practice, we usually process the top (or the bottom) together with the thickness.

Multivariate structural analysis¶

Both simple variograms and the cross-variogram are modeled in the scope of the linear model of coregionalization.

The model is composed of a spherical basic structure and an exponential component. Both of them are anisotropic.

Joint estimation¶

The Top and Thickness variables are estimated jointly using a Co-Kriging procedure. As for Kriging in a previous paragraph, we use the moving neighborhood composed of 5 closest points per angular sector, 6 sectors and a neighborhood maximum circle of 5km. This leads to co-kriging systems of maximum dimension equal to 60 equations (+ the two universality conditions).

The thickness cokriged map is displayed in the next figure together with the layer thickness information calculated at wells.

The bottom surface is simply derived by adding the resulting thickness map to the top map. Unfortunately, this technique does not provide any information regarding the standard deviation of the estimation of the bottom surface.

Jointly estimated surfaces can now be displayed using the 3-D viewer. In addition, the data information is overlaid to testify that the result of this cokriging procedure produce exact interpolators (top and bottom data are honored).