Spatial correlation and the variogram

Quick refresher

A histogram describes which values occur. A variogram describes how they are arranged. Two reservoirs with the same permeability histogram can behave quite differently in development.
The experimental variogram is half the mean squared difference between the values in pairs at each lag (distance step). Its points are unreliable where there are few pairs. We usually build lags up to half the size of the area.
There are three parameters. The nugget effect is the jump at the origin (measurement errors and variability smaller than the data spacing). The sill is the level that the variogram reaches (usually the variance of the data). The range is the distance at which it reaches the sill.
The models differ in their behaviour near the origin. The spherical and exponential models grow linearly (the property changes sharply), and the Gaussian model grows parabolically (it changes smoothly).
The variogram depends on direction. The ranges along elongated bodies are larger than across them, and in the vertical they are tens to hundreds of times smaller than in the horizontal.
Kriging weights the wells by their closeness to the point and to one another, and it gives an estimate with an error variance. A kriging map is smoothed, so for flow calculations we build realisations of stochastic simulation (several equally likely versions of the reservoir).

Key relationships

\hat\gamma(h) = \dfrac{1}{2N(h)}\sum_{i=1}^{N(h)}\left[z(x_i + h) - z(x_i)\right]^2

experimental variogram

\gamma(h) = C_0 + C\left[\dfrac{3h}{2a} - \dfrac{1}{2}\left(\dfrac{h}{a}\right)^3\right],\ h \le a

spherical model; for $h > a$ it equals the sill $C_0 + C$

C(h) = \sigma^2 - \gamma(h),\qquad \rho(h) = 1 - \gamma(h)/\sigma^2

covariance and correlogram

\sum_j \lambda_j \gamma_{ij} + \mu = \gamma_{i0},\quad \sum_j \lambda_j = 1,\quad \sigma_K^2 = \sum_i \lambda_i\gamma_{i0} + \mu

the ordinary kriging system and its variance

Common pitfalls

Fitting the model to the far lags, where there are few pairs.
Treating the nugget effect as measurement error and removing it.
Carrying the vertical variogram from wells over to the horizontal directions.
Building a variogram from data with a trend without removing the trend.
Using a kriging map as a realistic picture of the reservoir for a displacement calculation.

After this lesson you can

build an experimental variogram from a series of data, and explain what each point shows and how much we can trust it
read a variogram (nugget effect, sill, range, trend and cyclicity) and link these features with the geology and with the measurements
fit a spherical, exponential or Gaussian model, and explain how they differ near the origin (the start of the graph)
move from the variogram to the covariance and the correlogram, and estimate the effective number of independent data
describe anisotropy (change with direction) with directional variograms and an ellipse of ranges
calculate an ordinary kriging estimate from several wells, explain its weights and variance, and say how it differs from sequential Gaussian simulation

Spatial correlation: what the histogram does not show

Descriptive statistics tell us which values occur in a reservoir (Descriptive statistics in geology): the mean, the spread, the shape of the distribution (Distributions and sample size). But they say nothing about where these values are. High-permeability layers may be scattered in a fine patchwork, or they may stretch into continuous strips along which water breaks through from an injection well to a production well. The histograms of these two reservoirs are the same, but their development is different (correlation length).

Some reservoir properties do not change at random through space. Neighbouring values are similar. In geostatistics, such a property is called a regionalised variable. Values at neighbouring points are alike, and values at distant points are not. The variogram describes how fast the similarity disappears. We compare pairs of points a lag $h$ apart. The lag is a distance, and over an area or a volume it also has a direction.

The experimental variogram

Experimental variogram

\hat\gamma(h) = \frac{1}{2N(h)}\sum_{i=1}^{N(h)}\left[z(x_i + h) - z(x_i)\right]^2

where:

$\hat\gamma(h)$ experimental (semi)variogram at lag $h$ , square of the units of the property, for example %²
$h$ lag: the distance (and, over an area, the direction) between the points of a pair, m
$N(h)$ number of pairs of points a distance $h$ apart, —
$z(x_i),\ z(x_i + h)$ values of the property at the points of a pair, units of the property

Units: values in their own units, distances in metres

Holds when

The property is stationary in the studied area: the mean does not change in a regular way, and the differences depend only on the lag.
Irregularly placed data are grouped into classes of lag ( $h$ plus or minus a tolerance) and of direction (plus or minus an angular tolerance).
An estimate from lags with few pairs is unreliable. A practical rule is to use lags no longer than half the size of the area and to keep several tens of pairs per lag.

We put in the factor $1/2$ so that the variogram equals the variance at large distances, where the values are independent. The mean squared difference of two independent values is twice the variance. The variogram does not need the mean, which is its practical advantage over the covariance.

Worked exampleVariogram of porosity along a row of wells

Given (assumed data)

Wells: 14 wells of a row, on a straight line 300 m apart, row azimuth 80° (the compass direction, measured from north)
Mean porosity of the layer, % (in order): 18.5, 17.2, 16.7, 14.8, 15.3, 15.8, 16.6, 16.7, 19.4, 16.1, 16.9, 18.9, 19.1, 19.0

Solution

The mean is $\bar\phi = 17.21$ %, and the sample variance is $s^2 = 2.289$ %².
Lag 300 m: neighbouring wells, 13 pairs. The squared differences are 1.69, 0.25, 3.61, 0.25, 0.25, 0.64, 0.01, 7.29, 10.89, 0.64, 4.00, 0.04, 0.01. Their sum is 29.57 %², and $\hat\gamma(300) = 29.57/(2\cdot13) = 1.137$ %².
Lag 600 m: wells with one well between them, 12 pairs. The sum is 41.60, and $\hat\gamma(600) = 41.60/24 = 1.733$ %². Lag 900 m, 11 pairs: the sum is 50.22, and $\hat\gamma(900) = 2.283$ %².
Lags 1200, 1500 and 1800 m (10, 9 and 8 pairs): $\hat\gamma$ = 2.3075, 2.3261 and 2.2844 %².
As fractions of the variance: 0.50, 0.76, 1.00, 1.01, 1.02, 1.00. From 900 m on, the variogram stays near one.
A spherical model with a nugget effect of 0.3 %², a partial sill of 2.0 %² (sill 2.3 %²) and a range of 1100 m gives 1.10, 1.77 and 2.21 %² at lags of 300, 600 and 900 m, and 2.3 %² beyond the range. This is close to the points.

Answer: The variogram grows from 1.14 %² at 300 m and reaches the sill of about 2.3 %², the variance of the data, by 900–1200 m. The range along the row is about 1.1 km, and the nugget effect is about 0.3 %².

What it means. Wells more than 1–1.2 km apart are practically independent, and neighbouring wells are noticeably alike. The nugget effect here is not an error. With wells 300 m apart, it holds all the variability on a smaller scale, from laminae (very thin layers) to separate sand bodies. These wells cannot show that variability. Its estimate is an extrapolation to zero from the first points, so it is uncertain. The points beyond the range fluctuate around the variance: with 8–10 pairs the variogram is known only roughly, and we should not fit the model to the far lags. The data in this example are assumed for illustration.

CalculatedExperimental variogram of the row in the example (points) and the fitted spherical model: nugget effect 0.3 %², sill 2.3 %² (dashed; the variance of the data is 2.29 %²), range 1100 m (vertical dashed line).

On an irregular network of wells we group the pairs into classes of distance, for example 0–150, 150–450, 450–750 m and so on. For a directional variogram we also group them by angle: the class of azimuth 60° contains the pairs whose direction differs from 60° by no more than a tolerance, say 22.5°. The wider the tolerances, the more pairs we have, and the rougher the resolution in distance and direction. Vertically there are many data (well logs every 0.1–0.2 m), so we can build the vertical variogram reliably. Horizontally there are few wells. So we take horizontal ranges from other sources: geology, outcrop analogues (the same kind of rock exposed at the surface), horizontal wells, and seismic attributes (properties calculated from seismic data).

Nugget, sill and range: how to read a variogram

Nugget effect $C_0$: the jump of the variogram at the origin. The name comes from mining: a gold nugget in one sample makes the neighbouring samples unlike it. The nugget effect is made up of measurement errors and of variability over distances smaller than the data spacing. The larger it is, the less use neighbouring measurements are for interpolation (estimating values between measurements). A variogram that consists wholly of nugget means that, at the scale of the data, the values are not correlated.
Sill $C_0 + C$: the level that the variogram reaches. Under stationarity it is the variance of the property. If the area is small compared with the range, the variogram may not have time to reach the sill.
Range $a$: the distance at which the variogram reaches the sill, and the values stop being alike. It is the same correlation length as in the lesson on the heterogeneity of permeability.

Typical shapes of a variogram and what stands behind them

Shape: At the sill at once (pure nugget)
What it means: at the scale of the data the values are not correlated: the correlation length is smaller than the lag, or the data are noisy
What to do: get denser data; in the model, use random placement, and for averaging use the geometric mean (random patchwork)

Shape: A smooth rise to the sill
What it means: ordinary correlation with range $a$
What to do: fit a model; compare the range with the well spacing and the cell size

Shape: A rise without a sill, accelerating
What it means: trend: the mean changes in a regular way
What to do: remove the trend and build the variogram of the residuals

Shape: Dips below the sill at certain lags (hole effect)
What it means: cyclicity: repeating layers, units or bodies with a period equal to the position of the dip
What to do: include it in the model, or model the bodies as objects

Shape: Different sills in different directions
What it means: zonal anisotropy: for example, in the vertical direction different layers are crossed, and in the horizontal direction they are not
What to do: add a structure that acts in one direction only

Shape	What it means	What to do
At the sill at once (pure nugget)	at the scale of the data the values are not correlated: the correlation length is smaller than the lag, or the data are noisy	get denser data; in the model, use random placement, and for averaging use the geometric mean (random patchwork)
A smooth rise to the sill	ordinary correlation with range $a$	fit a model; compare the range with the well spacing and the cell size
A rise without a sill, accelerating	trend: the mean changes in a regular way	remove the trend and build the variogram of the residuals
Dips below the sill at certain lags (hole effect)	cyclicity: repeating layers, units or bodies with a period equal to the position of the dip	include it in the model, or model the bodies as objects
Different sills in different directions	zonal anisotropy: for example, in the vertical direction different layers are crossed, and in the horizontal direction they are not	add a structure that acts in one direction only

We take a dip as significant if it is clearly below the sill and supported by the geology: with a small number of pairs a variogram fluctuates even without any cyclicity.

Permeability data have a large spread and a log-normal distribution (its logarithm follows the bell-shaped normal curve, see log-normal distribution). So we build the variogram of permeability from its logarithm. Otherwise a few high values would decide the whole picture. The nugget effect also depends on the support, the volume to which a measurement refers. A probe that measures permeability over a centimetre of core sees laminae that a plug does not distinguish, and the plug sees those that disappear at the scale of a well test (scale of measurement).

Variogram models

The experimental variogram is known only at a few lags, but kriging and simulation need its value at any distance. So we draw a model through the points. It is a function chosen from a small set of permitted functions (the advanced section explains why we cannot use just any function). Three models are used most often.

Spherical, exponential and Gaussian models

\begin{aligned} \gamma_{\text{sph}}(h) &= C_0 + C\left[\frac{3h}{2a} - \frac{1}{2}\left(\frac{h}{a}\right)^{3}\right]\ \text{for}\ h \le a;\qquad C_0 + C\ \text{for}\ h > a\\ \gamma_{\text{exp}}(h) &= C_0 + C\left[1 - e^{-3h/a}\right]\\ \gamma_{\text{gauss}}(h) &= C_0 + C\left[1 - e^{-3h^2/a^2}\right] \end{aligned}

where:

$C_0$ nugget effect, units²
$C$ partial sill: the variability that is related to distance, units²
$C_0 + C$ sill, units²
$a$ range: for the spherical model, the distance at which it reaches the sill; for the exponential and Gaussian models, the distance at which it reaches 95% of the partial sill, m
$h$ lag, $h > 0$ ; at $h = 0$ the variogram is zero, m

Units: the variogram in squared units of the property, distances in metres

Holds when

A model is a sum of permitted structures with non-negative coefficients (nested structures).
The nugget effect is a jump at any arbitrarily small $h > 0$ , and the variogram itself is zero at zero.
The exponential and Gaussian models are often written with the parameter $a/3$ or $a/\sqrt3$ . When we move parameters between programs, we check which range is meant.

The models differ near the origin. The spherical and exponential models grow linearly: the property changes sharply from point to point, as the porosity or permeability of rocks does. The Gaussian model grows parabolically: the property changes smoothly, like the depth to the top of a layer or its thickness. Without a nugget effect, the Gaussian model makes the kriging system ill-conditioned (numerically unstable). The spherical model reaches the sill exactly at the range, and the exponential model approaches it slowly.

CalculatedThree models with the same sill of 1 and the same range of 1000 m: spherical (solid blue), exponential (dashed) and Gaussian (orange). The spherical and exponential models grow linearly near the origin, and the Gaussian model grows parabolically. At small lags the exponential model is the highest and the Gaussian model is the lowest.

We fit the model so that it passes close to the first points. There we have the most pairs, and there the fit matters most for kriging. The weights of the points are taken proportional to the number of pairs. We check the sill against the variance of the data, and the range against the geology: the size of the bodies and the distances between them (sizes of sand bodies). If the variogram rises in two steps, for example laminae and then units, the model is built from two structures, a short one and a long one.

Covariance, correlogram and the number of independent data

The same spatial correlation can be described by similarity instead of difference. The covariance $C(h)$ is the average product of the deviations from the mean, for values a lag $h$ apart. The correlogram $\rho(h)$ is their correlation coefficient (a number from −1 to 1 that shows how closely the values go together), as a function of the lag.

Relation between covariance, correlogram and variogram

C(h) = C(0) - \gamma(h),\qquad \rho(h) = \frac{C(h)}{C(0)} = 1 - \frac{\gamma(h)}{C(0)},\qquad C(0) = C_0 + C = \sigma^2

where:

$C(h)$ covariance of values a distance $h$ apart, units²
$C(0) = \sigma^2$ variance of the property, the sill of the variogram, units²
$\rho(h)$ correlogram: the correlation coefficient of values a distance $h$ apart, —
$\gamma(h)$ variogram, units²

Units: in squared units of the property; $\rho$ is dimensionless

Holds when

Second-order stationarity: a constant mean and a finite variance, which the variogram reaches at the sill.
With a trend, or with a variogram that has no sill, the covariance is not defined but the variogram is, and we use the variogram.

The variogram and the covariance are mirror descriptions: where the variogram rises, the covariance falls, and at the range the covariance becomes zero. With a nugget effect, the correlogram drops by a jump from 1 at $h = 0$ to $1 - C_0/\sigma^2$ at any arbitrarily small lag.

Correlated data carry less information than independent data: neighbouring values partly repeat each other. The effective number $N/(1 + (N-1)\bar\rho)$ shows how many independent values a series of $N$ correlated ones is worth. Here $\bar\rho$ is the mean correlation coefficient over all pairs (what similar plugs are worth). The variogram gives $\bar\rho$ without guesses “by eye”.

Worked exampleHow many independent wells there are in a row of fourteen

Given (assumed data)

Row: the 14 wells 300 m apart from the example above; standard deviation of porosity 1.513%
Model: spherical: nugget 0.3 %², sill 2.3 %², range 1100 m

Solution

Correlation coefficients from the model: $\rho(300) = 1 - 1.0979/2.3 = 0.523$ ; $\rho(600) = 1 - 1.7741/2.3 = 0.229$ ; $\rho(900) = 1 - 2.2068/2.3 = 0.041$ ; beyond the range, zero.
There are 91 pairs in the row: 13 neighbouring, 12 with one well between, 11 with two wells between. The mean correlation is $\bar\rho = (13\cdot0.523 + 12\cdot0.229 + 11\cdot0.041)/91 = 0.110$ .
The effective number of independent wells is $14/(1 + 13\cdot0.110) = 5.8$ .
The standard error of the mean porosity of the row is $1.513/\sqrt{5.8} = 0.63$ %, instead of $1.513/\sqrt{14} = 0.40$ % if we treat the wells as independent.

Answer: The fourteen wells of the row are worth about 5.8 independent ones, and the mean porosity from them is known to within about 0.6%, not 0.4%.

What it means. Correlation works both ways. For estimating means, neighbouring wells partly repeat each other and are “worth” less. For predicting between wells, this very similarity is useful, and kriging is based on it. The mean correlation was calculated from the model, not from the data themselves, and it depends on the model: the longer the correlation at the same spacing, the fewer independent values. The data in this example are assumed for illustration.

Anisotropy: the variogram by direction

Sand bodies are elongated along the direction of the palaeocurrent (the direction in which the water flowed when it laid down the sand). Reefs are elongated along the edge of a platform (a shallow marine shelf), and layers along the bedding. So the variogram depends on direction. We build it for several azimuths (directional variograms) and for the vertical. If the sills in different directions are the same but the ranges differ, the anisotropy is called geometric: the ranges in different directions form an ellipse, and in three dimensions an ellipsoid. If the sills differ as well, the anisotropy is called zonal. This happens when the vertical variogram crosses different layers and the horizontal one runs within a single layer.

Reduced distance and range in a direction for geometric anisotropy

h' = \sqrt{\left(\frac{h_{\parallel}}{a_{max}}\right)^2 + \left(\frac{h_{\perp}}{a_{min}}\right)^2},\qquad a(\alpha) = \left[\frac{\cos^2(\alpha - \alpha_0)}{a_{max}^2} + \frac{\sin^2(\alpha - \alpha_0)}{a_{min}^2}\right]^{-1/2}

where:

$h'$ the distance between the points as a fraction of the range: the model is evaluated at $h/a = h'$ , —
$h_{\parallel},\ h_{\perp}$ components of the distance along the major axis of the ellipse and across it, m
$a_{max},\ a_{min}$ correlation ranges along the major and minor axes, m
$\alpha,\ \alpha_0$ azimuth of the direction and azimuth of the major axis, °
$a(\alpha)$ correlation range in the direction $\alpha$ , m

Units: distances in metres, azimuths in degrees

Holds when

Geometric anisotropy: the sill and the nugget effect are the same in all directions.
The axes of the ellipse coincide with the principal directions, usually along the elongation of the bodies and across it. In three dimensions we add a vertical axis with a much smaller range.
For zonal anisotropy, we add to the model a structure that acts in one direction only.

We convert the distance into fractions of the range along each axis, so a pair of wells along the bodies turns out to be “closer” than the same pair across them. The ellipse of ranges in plan view shows how far, in each direction, the values are still alike.

Worked exampleHow alike are two wells 600 m apart

Given (assumed data)

Directional variograms from all the wells of the field: range 1600 m at azimuth 60°, along the sand bodies, and 500 m at azimuth 150°; nugget 0.3 %², sill 2.3 %², spherical model
Pair of wells: 600 m apart at azimuth 90°
Row from the first example: azimuth 80°

Solution

Range along the row (azimuth 80°, which is 20° from the major axis): $a = \left[\cos^2 20°/1600^2 + \sin^2 20°/500^2\right]^{-1/2} = 1109$ m. This is almost the same as the variogram of the row gave (1100 m).
The pair of wells is 30° off the major axis: the components of the distance are $600\cos30° = 519.62$ m and $600\sin30° = 300.0$ m.
The reduced distance is $h' = \sqrt{(519.62/1600)^2 + (300.0/500)^2} = 0.6823$ .
The spherical structure: $1.5\cdot0.6823 - 0.5\cdot0.6823^3 = 0.8646$ ; $\gamma = 0.3 + 2.0\cdot0.8646 = 2.029$ %², and the correlation coefficient is $1 - 2.029/2.3 = 0.118$ .
The isotropic model of the row (range 1100 m in all directions) would give $\gamma(600) = 1.774$ %² and a correlation of 0.229.

Answer: Wells 600 m apart at azimuth 90° are only weakly correlated:

\rho \approx 0.12

, half of what the isotropic model of the row gives (0.23). The range in this direction is only 879 m.

What it means. Anisotropy decides how we interpolate: along the bodies, values are carried for a kilometre and a half, and across them, for half a kilometre. A row of wells that runs almost along the bodies “sees” the long range, and carrying its variogram over to other directions would overstate the continuity. So we build directional variograms from all pairs with an angular tolerance, and we check the major axis against the geology: the direction of the palaeocurrent and the elongation of the bodies (shape of bodies). The data in this example are assumed for illustration.

CalculatedTop: the ellipse of correlation ranges from the example in plan view: 1600 m at azimuth 60°, 500 m at azimuth 150°. The arrows show the direction of the row of wells (80°, range about 1110 m) and of the pair of wells (90°, about 880 m). Bottom: variogram models in three directions: across the bodies (150°, steep), along the row (80°) and along the bodies (60°, gentle). The points are the experimental variogram of the row.

The vertical ranges are usually tens to hundreds of times smaller than the horizontal ones: the layers are thin but extensive (correlation length in the vertical and the horizontal). So we cannot carry the vertical variogram from logs and core over to the horizontal, and we make the model grid thin in the vertical and coarse across the area (Grids in reservoir simulation models).

Kriging: estimating between wells from the variogram

Kriging is a way to estimate a property at a point with no measurement, as a weighted average of the data, choosing the weights from the variogram. The method is named after a mining engineer who applied statistics to the estimation of ore reserves, and it was developed in geostatistics in the 1960s. Unlike inverse distance weighting (a simpler method in which nearer data get more weight, see gridding), the weights of kriging depend on more than the distance to the point. They also depend on how the data are placed relative to one another, and on the anisotropy. Together with the estimate we get the expected error.

The ordinary kriging system

\sum_{j=1}^{n}\lambda_j\,\gamma(x_i - x_j) + \mu = \gamma(x_i - x_0),\quad i = 1, \dots, n;\qquad \sum_{j=1}^{n}\lambda_j = 1

where:

$\lambda_j$ weight of the $j$ th data point, —
$\mu$ Lagrange multiplier, a number that accounts for the condition that the weights add up to one, units²
$\gamma(x_i - x_j)$ variogram between data points, units²
$\gamma(x_i - x_0)$ variogram between a data point and the point of estimation $x_0$ , units²
$n$ number of data points in the neighbourhood, —

Units: variogram in squared units of the property

Holds when

The mean is unknown but constant in the neighbourhood of the estimate (ordinary kriging). When the mean is known, we use simple kriging: without the condition on the sum of the weights, and with the covariance instead of the variogram.
The variogram is known and refers to the same quantity and the same support as the data.
The estimate is a linear combination of the data. Measurement errors are included in the nugget effect.

The condition on the sum of the weights makes the estimate unbiased (not too high or too low on average), and the system minimises the expected squared error of the estimate. The right-hand side reflects how close each data point is to the point of estimation, and the matrix on the left reflects how close the data points are to one another: data that repeat each other share the weight.

The kriging estimate and its variance

z^*(x_0) = \sum_{i=1}^{n}\lambda_i\,z(x_i),\qquad \sigma_K^2 = \sum_{i=1}^{n}\lambda_i\,\gamma(x_i - x_0) + \mu

where:

$z^*(x_0)$ estimate of the property at the point $x_0$ , units of the property
$\sigma_K^2$ kriging variance, the expected squared error of the estimate, units²

Units: as for the property and the variogram

Holds when

The weights and $\mu$ are the solution of the ordinary kriging system.
The kriging variance is correct if the variogram is correct. The distribution of the errors need not be normal.

The kriging variance depends only on how the data are placed and on the variogram, not on the values themselves. We can calculate it before drilling and choose the well location that would reduce the uncertainty most. The downside is that it does not show that the values are more scattered in one part of the reservoir than in another. At a data point, kriging returns the measured value, and the variance is zero.

Worked examplePorosity at a planned well from three neighbours

Given (assumed data)

Well 1: 450 m east of the planned well; porosity 16.2%
Wells 2 and 3: to the west: 420 m west and 150 m north, 18.9%; 500 m west and 320 m north, 19.4%
Variogram: the model of the row: spherical, nugget 0.3 %², sill 2.3 %², range 1100 m; for simplicity we ignore anisotropy

Solution

Distances to the planned well: $d_1 = 450$ m, $d_2 = \sqrt{420^2 + 150^2} = 446.0$ m, $d_3 = \sqrt{500^2 + 320^2} = 593.6$ m. Between the wells: $d_{12} = 882.8$ , $d_{13} = 1002.4$ , $d_{23} = 187.9$ m.
The variogram from the model: to the planned well $\gamma_{10} = 1.459$ , $\gamma_{20} = 1.450$ , $\gamma_{30} = 1.762$ %²; between the wells $\gamma_{12} = 2.191$ , $\gamma_{13} = 2.277$ , $\gamma_{23} = 0.807$ %².
The system: $2.191\lambda_2 + 2.277\lambda_3 + \mu = 1.459$ ; $2.191\lambda_1 + 0.807\lambda_3 + \mu = 1.450$ ; $2.277\lambda_1 + 0.807\lambda_2 + \mu = 1.762$ ; $\lambda_1 + \lambda_2 + \lambda_3 = 1$ . The solution: $\lambda_1 = 0.483$ , $\lambda_2 = 0.426$ , $\lambda_3 = 0.091$ , $\mu = 0.318$ %².
The estimate: $0.483\cdot16.2 + 0.426\cdot18.9 + 0.091\cdot19.4 = 17.64$ %.
The kriging variance: $0.483\cdot1.459 + 0.426\cdot1.450 + 0.091\cdot1.762 + 0.318 = 1.80$ %², so the standard deviation is 1.34%. Without any wells it would be $\sqrt{2.3} = 1.52$ %.
Inverse distance weighting ( $w \propto 1/d^2$ ) would give the weights 0.3857, 0.3927 and 0.2216 and the estimate 17.97%.

Answer: Kriging estimates the porosity at the planned well as 17.6% with a standard deviation of 1.3%. Inverse distance weighting would give 18.0%.

What it means. Well 1 received almost half of the weight, although it is not closer than well 2. Wells 2 and 3 are 188 m apart and largely repeat each other, so their joint weight is 0.52, not 0.61 as in inverse distance weighting. Kriging takes weight away from clusters of data (this is called declustering). Well 3 is also hidden behind well 2 and received only 0.09. This is the screening effect. The standard deviation of 1.3% differs little from the 1.5% that we have without wells: with a range of 1.1 km and a noticeable nugget effect, wells 450–600 m away say little about the planned point. Ignoring anisotropy changes the weights: wells that lie in the direction of elongation of the bodies would in fact weigh more. The data in this example are assumed for illustration.

CalculatedPlan view of the wells of the example and their weights: for kriging (the thickness of the lines to the planned well is proportional to the weight) and for inverse distance weighting (IDW). Kriging gives the pair of western wells less weight, and the farther of the pair is shielded by the nearer one.

Exactness at data points. Kriging passes through the data if the nugget effect is assigned to spatial variability. If instead it describes measurement error, the estimate at a data point is smoothed.
Negative weights. Shielded points can receive a small negative weight. The estimate can then go outside the range of the data, and for porosity or the net-to-gross ratio (the share of reservoir rock) we sometimes have to limit it.
Smoothing. The estimate is an average over the possible values, so a kriging map is smoother than reality: its variance is smaller than the variance of the data, and the extreme values and the connectivity of high-permeability zones are understated. This is not a fault of the method but a property of the best estimate, and it is the reason why we build simulations for flow calculations.
A second variable. If we have dense data on a related quantity, such as a seismic attribute or a thickness from seismic, we use cokriging or kriging with an external trend (methods that use the second variable): the sparse well values are refined with the dense seismic data.

Stochastic simulation and property models

For a reservoir simulation model (a flow model), it is not only the average values that matter but also the variability: continuous high-permeability strips, tight barriers, the connectivity of bodies (reservoir connectivity). A kriging map smooths them out, so water breakthrough calculated from it comes too late and the sweep (the share of the reservoir that the water reaches) comes out too good. Stochastic simulation (simulation that uses random numbers) builds many versions of the distribution of a property, called realisations. Each of them matches the data at the wells, reproduces their histogram and variogram, and is equally likely. The differences between realisations show the uncertainty that the wells do not remove.

Sequential Gaussian simulation

1
Normalise
Transform the data so that they follow the standard normal distribution (the bell curve with mean 0 and standard deviation 1), matching them quantile by quantile (the normal score transform). Then build the variogram of the transformed values.
2
Set a path
Visit the grid nodes in a random order.
3
Estimate
At the next node, find the mean and the variance by simple kriging from the wells and the neighbouring nodes that are already filled.
4
Draw
Choose a random value from the normal distribution with this mean and variance, and write it into the node. From then on the node counts as data.
5
Transform back
When all the nodes have been visited, transform the values back to the units of the property. A different random path and different random numbers give a different realisation.

CalculatedPorosity along a line 3 km long between three wells (dots): the simple kriging estimate (thick line) with a band of ±1 standard deviation (dashed), and three realisations of sequential Gaussian simulation (thin lines) using the variogram of the lesson; nodes every 50 m, and the values are taken to be normally distributed. The realisations pass through the wells and keep the variability. Their jaggedness from node to node is set mainly by the nugget effect of the model. Between the wells, kriging tends towards the mean of 17.2%.

In a geological model we build the properties step by step. First we distribute facies or rock types, by indicator simulation (simulation of categories) using variograms, or by object modelling when the bodies have a known shape (channels, lobes) (Reservoir architecture and performance). Then, inside each facies, we simulate porosity, with its own simulation and its own variogram. We model permeability together with porosity: using the relation between them with its scatter (Reservoir quality: geological controls on porosity and permeability), or by joint simulation with porosity as the second variable. Seismic attributes serve as the second variable for porosity and the net-to-gross ratio. The finished model is upscaled (Permeability upscaling for single-phase flow), and flow is calculated on several realisations. The range of the production forecast depends on this (Reservoir simulation in practice: case studies, Probabilistic volumetrics by Monte Carlo simulation).

We compare the cell size with the range. If the range is smaller than the cell, the properties of neighbouring cells are independent, and the model looks like a random patchwork. The variability inside the cell is then lost, and we take it into account in upscaling. If several cells or more fit into the range, continuous bodies are visible in the model, and we can judge connectivity from it. A variogram found from core or logs (a support of centimetres to decimetres) is brought to the size of the cell before simulation (see the advanced section).

Terms and symbols

Variogram $\gamma(h)$: Half the mean squared difference between the values of a property at pairs of points a lag $h$ apart. It is also called the semivariogram.
Lag $h$: The distance (and direction) between the points of a pair. For irregular data, it is a class of distance with a tolerance.
Nugget effect $C_0$: The jump of the variogram at the origin, made up of errors and variability smaller than the data spacing.
Sill: The level that the variogram reaches. Under stationarity, it is the variance.
Range $a$: The distance at which the variogram reaches the sill. It is the correlation length.
Correlogram $\rho(h)$: The correlation coefficient of values a distance $h$ apart: $1 - \gamma(h)/\sigma^2$ .
Geometric and zonal anisotropy: In geometric anisotropy, the ranges differ between directions while the sill stays the same. In zonal anisotropy, the sills themselves differ.
Ordinary and simple kriging: Estimation by a weighted average, with weights taken from the variogram. Ordinary kriging is used when the constant mean is unknown (the weights add up to one), and simple kriging when the mean is known.
Kriging variance $\sigma_K^2$: The expected squared error of the estimate. It depends on how the data are placed and on the variogram, but not on the values.
Sequential Gaussian simulation: A way to build realisations of a property that have a given histogram and variogram and match the data at the wells.
Support: The volume to which a value refers. It can be a point on core, an interval of a log, or a model cell.

Common mistakes and misconceptions

Wrong

“The variogram at far lags shows the real structure of the reservoir.”

In fact

At far lags there are few pairs, and the variogram fluctuates at random. We fit the model to the first lags, which have many pairs. A practical rule is to use lags no longer than half the size of the area.

Wrong

“The sill is always equal to the variance of the data.”

In fact

Only under stationarity, and if the area is much larger than the range. With a trend the variogram rises without a sill, with zonal anisotropy the sills in different directions are different, and a small area does not let the variogram reach the sill.

Wrong

“The nugget effect is measurement error, and we should remove it.”

In fact

It also contains variability on a scale smaller than the data spacing. If we underestimate it, kriging and simulation will show the reservoir as more continuous than it is.

Wrong

“A kriging map is a realistic picture of the reservoir.”

In fact

It is the best estimate at each point, but it is smoothed: the extreme values and the connectivity are understated. For displacement calculations we use realisations of simulation.

Wrong

“A vertical variogram from logs is good for the horizontal directions too.”

In fact

Horizontal ranges are tens to hundreds of times larger than vertical ones. We take them from geology, analogues, horizontal wells and seismic data.

In depth: advanced detail

Where the kriging system comes from

The error of the estimate is $e = \sum_i \lambda_i Z(x_i) - Z(x_0)$ . If $\sum_i \lambda_i = 1$ , the unknown mean cancels in the error, and the estimate is unbiased. The variance of such an error can be expressed through the variogram: $\mathrm{Var}\,e = 2\sum_i \lambda_i\,\gamma_{i0} - \sum_i\sum_j \lambda_i\lambda_j\,\gamma_{ij}$ .

We look for the smallest variance under the condition $\sum_i \lambda_i = 1$ by the method of Lagrange multipliers (a standard way to find a minimum when a condition must hold). We take the function $\mathrm{Var}\,e - 2\mu\left(\sum_j \lambda_j - 1\right)$ and set its derivative to zero for each $\lambda_i$ . This gives $n$ equations $\sum_j \lambda_j\gamma_{ij} + \mu = \gamma_{i0}$ and the condition on the sum of the weights, which is the system of the lesson. If we multiply each equation by $\lambda_i$ and add them up, we find $\sum_i\sum_j\lambda_i\lambda_j\gamma_{ij} = \sum_i\lambda_i\gamma_{i0} - \mu$ , and the error variance with the weights found is $\sigma_K^2 = \sum_i \lambda_i\gamma_{i0} + \mu$ .

Why we cannot draw the variogram just as we like

The variance of any combination of values $\sum_i w_i Z(x_i)$ with $\sum_i w_i = 0$ is equal to $-\sum_i\sum_j w_i w_j\,\gamma_{ij}$ , and it cannot be negative. An arbitrary curve drawn through the experimental points may not satisfy this condition, and then kriging gives a negative variance or unstable weights. So we use proven functions: the spherical, exponential and Gaussian ones, the power function (without a sill), the nugget effect, and their sums with non-negative coefficients. The spherical model is permitted in one, two and three dimensions. The power model $h^{\omega}$ is permitted for $0 < \omega < 2$ , and it describes a property without a sill, including one with variability at all scales.

The variogram and the cell size: the support effect

The value in a model cell is an average over its volume $v$ , not a point measurement. The variance of the cell averages is smaller than the point variance: $\sigma_v^2 = \sigma^2 - \bar\gamma(v, v)$ , where $\bar\gamma(v, v)$ is the mean value of the variogram over all pairs of points inside the cell. If the cell is much smaller than the range, $\bar\gamma$ is only a little larger than the nugget effect, and the variance falls almost only by that amount. If the cell is comparable to the range or larger, the values in the cells are smoothed, and the variogram of the cells is lower and gentler than the point variogram. The nugget effect practically disappears on averaging over a cell.

So before simulation we recalculate a variogram found from core or logs to the cell size (we regularise it), and when the model is upscaled, the variance of the properties falls once more (Permeability upscaling for single-phase flow). The same reason explains why permeability from a well test is less scattered than permeability from core: the test averages over a huge volume.

Connections to other topics

Reservoir heterogeneity and permeability averaging — the correlation length and the placing of areas of different permeability.
Distributions and sample size — distributions, the standard error and the effective number of independent data.
Descriptive statistics in geology — the variance and the measures of spread, to which the variogram comes at the sill.
Structure maps and contouring — map building: inverse distance weighting and minimum curvature.
Reservoir architecture and performance — the shape and size of bodies, connectivity and object models.
Grids in reservoir simulation models — the model grid and the cell size.
Permeability upscaling for single-phase flow — upscaling of models and averaging of properties over cells.
Special topics in upscaling — when upscaling does not work: correlation length against cell size.
Probabilistic volumetrics by Monte Carlo simulation — probabilistic estimation of hydrocarbons in place, and uncertainty.
Reservoir simulation in practice: case studies — calculations on sets of geostatistical realisations.
Practical structural geology — the structure of a field and data of different scales.

Check yourself

1
The variogram of the mean porosity from wells 400 m apart starts with a jump equal to 40% of the sill. What does this mean?
1. A.Porosity has a trend
2. B.The laboratory made a mistake, and the data must be thrown out
3. C.A noticeable part of the variability is at distances smaller than 400 m, or is measurement error; neighbouring wells are alike only partly
4. D.The range is 400 m
Show the answer
Correct answer: C. A noticeable part of the variability is at distances smaller than 400 m, or is measurement error; neighbouring wells are alike only partly
The nugget effect combines measurement errors and variability smaller than the data spacing. With a spacing of 400 m, it contains all the bodies and layers smaller than this size. Neighbouring values are still correlated, but more weakly: even at a small lag the correlation coefficient is no more than $1 - 0.4 = 0.6$ .
2
Porosity at points 100 m apart: 20, 22, 21, 18, 19%. Find the experimental variogram at lags of 100 and 200 m, and compare it with the sample variance.
Show the answer
Lag 100 m: the differences are 2, −1, −3, 1; the squares are 4, 1, 9, 1, with a sum of 15; $\hat\gamma(100) = 15/(2\cdot4) = 1.875$ %². Lag 200 m: the differences are 1, −4, −2; the sum of squares is 21; $\hat\gamma(200) = 21/6 = 3.5$ %². The sample variance: the mean is 20, the sum of squared deviations is 10, and $s^2 = 10/4 = 2.5$ %². The variogram at 200 m is above the variance. With three pairs this is an ordinary fluctuation, and we must not fit a model to such a point.
3
Which variogram model describes a property that changes most smoothly over short distances?
1. A.Gaussian
2. B.Pure nugget effect
3. C.Spherical
4. D.Exponential
Show the answer
Correct answer: A. Gaussian
The Gaussian model grows parabolically near the origin: at small lags the values hardly differ, as with the depth to the top of a layer or its thickness. The spherical and exponential models grow linearly, and a pure nugget effect means there is no correlation at the scale of the data.
4
The sill of a variogram is 4 %², and at a lag of 300 m the variogram is 3 %². What are the covariance and the correlation coefficient of values 300 m apart?
Show the answer
$C(300) = 4 - 3 = 1$ %², $\rho(300) = 1 - 3/4 = 0.25$ . The values at this distance are only weakly correlated: the covariance is just a quarter of the variance.
5
Two wells stand 50 m apart, both 400 m west of the planned point; a third is 400 m east. How will kriging distribute the weights, compared with inverse distance weighting, which gives the pair 2/3 in total and the single well 1/3?
1. A.The pair gets noticeably less than 2/3 in total, and the single well noticeably more than 1/3
2. B.The pair gets more than 2/3: two wells are more reliable than one
3. C.The single well gets zero
4. D.The same: the weights depend only on the distance to the point
Show the answer
Correct answer: A. The pair gets noticeably less than 2/3 in total, and the single well noticeably more than 1/3
Close wells repeat each other, and kriging divides the weight between them as if they were almost one source of information. With the model of the lesson the pair gets about 0.53 and the single well about 0.47. Inverse distance weighting does not see this and overrates clusters of data.
6
Why do we build several realisations of stochastic simulation for a reservoir simulation model, and not a kriging map?
Show the answer
Kriging gives the best estimate at each point, but a smoothed one: its variance is smaller than the variance of the data, and the extreme values and continuous high-permeability strips disappear, although they decide water breakthrough and sweep. The realisations reproduce the histogram and the variogram and match the wells, and their differences show the uncertainty of the forecast. We calculate flow on several realisations and get a range of production, not a single number.

Updated: 19 October 2026