Chapter 1

Geostatistics in Hydrology: Kriging

interpolation

Hydrologic properties, such as rainfall, aquifer characteristics (porosity, hydraulic conductivity,

transmissivity, storage coefficient, etc.), effective recharge and so on, are all functions of space
(and time) and often display a high spatial variability, also called heterogeneity. This variability
is not in general random. It is a general rule that these properties display a so called scale
effect, i.e., if we take measurements at two different points the difference in the measured
values dicreases as the two points come closer to each other.
It is convenient in certain cases to consider these properties as random functions having
a given spatial structure, or in other word having a given spatial correlation, which can be
conveniently described using appropriate statistical quantities. These variables are called re-
gionalized variables [5].
The study of regionalized variables starts from the ability to interpolate a given field starting
from a limited number of observation, but preserving the theoretical spatial correlation. This is
accomplished by means of a technique called kriging developed by [4, 5] and largely applied in
hydrology and other earth science disciplines [3, 6, 7, 1] for the spatial interpolation of various
physical quantities given a number of spatially distributed measurments.
Although theoreticallly kriging cannot be considered superior to other surface fitting tech-
niques [7] and the use of a few arbitrary parameters may lead absurd results, this method is
capable of obtaining objective interpolations evaluating at the same time the quality of the
results.
In the next sections we discuss first the statistical hypothesis that are needed to develop the
theory of kriging and then we proceed at describing the method under different assumptions
and finally report a few applications.

1.1 Statistical assumptions

Let Z(~x, ) be a random function (RF) simulating our hydrological quantity of interest. We
recall that with this notation we imply that ~x denotes a point in space (two or three dimensional
1
2 CAPTER 1. GEOSTATISTICS IN HYDROLOGY

space), while denotes the state variable in the space of the realizations. In other words:

Z(~x0 , ) is a random variable at point ~x0 representing the entire set of realizations of the
RF at point ~x0 ;

Z(~x, 1 ) is a particular realization of the RF Z;

Z(~x0 , 1 ) is a measurement point.

Given a set of sampled values Z(~xi , 1 ), i = 1, 2, . . . we want to reconstruct Z(~x, 1 ), i.e. a

possible realization of Z(~x, ).

1.2 Kriging for Weak or Second Order stationary RF

An RF is said to be second order stationary if:

1. Constant mean:
E[Z(~x, )] = m (1.1)

2. the autocovariance (another name for the covariance) is a function of the distance between
the reference points ~x1 and ~x2 :

cov[~x1 , ~x2 ] = E[(Z(~x1 , ) m)(Z(~x2 , ) m)] = C(~h) (1.2)

In practice second order stationarity implies that the first two statistical moments (expected
value and covariance) be translation invariant. Note that by saying that E[Z(~x, )] = m we
imply that effectively the expected value taken on all the possible realizations does not vary
with space. However, for a given realization Z(~x, 1 ) is a function of ~x.
Because of (1.1), the covariance (1.2) can be written as:

cov[~x1 , ~x2 ] = E[(Z(~x1 , ) m)(Z(~x2 , ) m)]

= E[Z(~x1 , )Z(~x2 , )] mE[Z(~x2 , )] E[Z(~x1 , )]m + m2
= E[Z(~x1 , )Z(~x2 , )] m2 (1.3)

Obviously if ~h = 0 we have the definition of variance, also called the dispersion variance:

C(0) = var[Z] = Z2

For simplicity of notation, from now on, we will denote x = ~x, h = ~h and we will drop the
variable in the RF.
1.2 SECOND ORDER STATIONARY RF 3

Case with m and C(~h) known. If E[Z(x)] = m and C(h) are known, then we can define a
new variable Y (x) with zero mean:

Y (x) = Z(x) m
E[Y (x)] = 0

Given the observed values:

x1 x2 . . . x n
Y1 Y2 . . . Y n
with Yi = Y (xi ) being the observation at point xi , we look for a linear estimator Y (x0 ) of
Y (x0 ) at point x0 using the observed values. The form of the estimator is:
n
X

Y (x0 ) = i Yi (1.4)
i=1

Note that the estimator (1.4) is a realization of the RF

n
X
Y (x0 , ) = i Y (xi , )
i=1

The weights i are calculated by imposing that the statistical error

(x0 ) = Y (x0 ) Y (x0 )

has zero expecetd value and minimum variance:

E[(x0 )] = 0
var[(x0 )] = E[(Y (x0 ) Y (x0 ))2 ] = minimum

Substituting eq. (1.4) in the expression of the variance we have:

n
X
2
E[(Y (x0 ) Y (x0 )) ] = E[( i Yi Y 0 ) 2 ]
i=1
n
X n
X
= E[( i Yi Y0 )( i Yi Y0 )]
i=1 i=1
n
X n
X n
X
= E[( i Yi )( i Yi )] 2E[ i Yi Y0 ] + E[Y02 ]
i=1 i=1 i=1
n
XX n n
X
= i j E[Yi Yj ] 2 i E[Yi Y0 ] + E[Y02 ]
i=1 j=1 i=1

but, since m = 0
E[Yi Yj ] = C(xi xj ) + m2 = C(xi xj )
4 CAPTER 1. GEOSTATISTICS IN HYDROLOGY

and
E[Y02 ] = C(0) = var[Y ]
is the dispersion variance of Y . Then:

n X
X n n
X
2
E[(Y (x0 ) Y (x0 )) ] = i j C(xi xj ) 2 i C(xi x0 )
i=1 j=1 i=1

+ C(0)

The minimum is found by setting to zero the first partial derivatives:

Xn


E[(Y (x0 ) Y (x0 ))2 ] = 2 j C(xi xj ) 2C(xi x0 ) = 0
i j=1

j = 1, . . . , n

This yields a linear system of equations:

C = b (1.5)

where matrix C is given by:

C(0) C(x1 x2 ) . . . C(x1 xn )

C=




C(xn x1 ) C(0)

and the right hand side vector b is given by:

C(x1 x0 )


b=


C(xn x0 )

Matrix C is the spatial covariance matrix and does not depend upon x0 . It can be shown
that if all the xj s are distinct then C is positive definite, and thus the linear system (1.5)
can be solved with either direct or iterative methods. Once the solution vector is obtained,
equation (1.4) yields the estimation of our regionalized variable at point x0 . Thus the calculated
value for is actually function of the estimation point x0 . If we want to change the estimation
point x0 , for example if we need to obtain a spatial distribution of our regionalized variable,
we need to solve the linear system (1.5) for different values of x0 . In this case it is convenient
to factorize matrix C using Cholesky decomposition and then proceed to the solution for the
different right hand side vectors.
1.3. KRIGING WITH THE INTRINSIC HYPOTHESIS 5

Evaluation of the estimation variance The estimation variance is defined as the variance
of the error:
 = Y0 Y0
Hence:
var[Y0 Y0 ] = E[(Y0 Y0 )2 ] E[(Y0 Y0 )]2
but: X
E[(Y0 Y0 )] = E[Y0 ] E[Y0 ] = i E[Yi ] E[Y0 ] = 0 (E[Y ] = m = 0)
i
and since X
j C(xi xj ) = C(xi x0 )
j

we finally obtain:
var[Y0 Y0 ] = E[(Y0 Y0 )2 ]
XX X
= i j C(xi xj ) 2 i C(xi x0 ) + C(0)
i j i
X
= i C(xi x0 ) + C(0)
i
X
= var[Y ] i C(xi x0 )
i
P
which, since i i C(xi x0 ) > 0, shows that the estimation variance of Y0 is smaller than the
dispersion variance of Y (the real variance of the RF). In statistical terms, we can interpret
this result by saying that since we have observed Y at some points xi then the uncertainty on
Y decreases.
It is important to remark the difference between estimation and dispersion variance. The
latter is representative of the variation interval of the RF Y within the interpolation domain,
while the estimation variance represents the residual uncertainty in the estimation of the real-
ization Y0 of Z when n observations are available. The dispersion variance is a constant, while
the estimation variance varies from point to point and is zero at the observation points.
Our original variable was Z = Y + m and its estimate is thus:
X
Z0 = m + i (Zi m)
i
X
var[Z0 Z0 ] = var[Z0 ] i C(xi x0 )
i

1.3 Kriging with the intrinsic hypothesis

The hypothesis of second order stationarity of the RF is not always satisfied, for example C(0)
increases with the distance, violating hypothesis (1.2). In this case the Intrinsic hypothesis
must be used, in which we assume that the first order increments
= Y (x + h) Y (x)
6 CAPTER 1. GEOSTATISTICS IN HYDROLOGY

range

High spatial correlation

Low spatial correlation range

(h)
C(h)

sill

h h

Figure 1.1: Behavior of the covariance as a function of distance (left) and the corresponding
variogram (right).

are second order RF:

E[Y (x + h) Y (x)] = m(h) = 0 (1.6)

var[Y (x + h) Y (x)] = 2(h) (1.7)

where the function (h) is called the variogram. If the mean m(h) is not zero an obvious
change of variable is required. The variogram is defined as the mean quadratic increment of
Y (x) (divided by 2) for any two points xi and xj separated by a distance h:

1 1
(h) = var[Y (x + h) Y (x)] = E[(Y (x + h) Y (x))2 ] (1.8)
2 2

and is related to the covariance function by:

1
(h) = E[(Y (x + h) Y (x))2 ]
2
1 1
= E[Y 2 (x + h)] E[Y (x + h)Y (x)] + E[Y 2 (x)]
2 2
= C(0) C(h)

The intrinsic hypothesis requires a finite value for the mean of Y (x) but not for its variance.
In fact, hypothesis (1.2), as changed into (1.3), implies (1.8), but not viceversa.
The covariance C(h) has a decreasing behavior as shown in Fig. 1.1. When C(h) is known
then the variogram can be directly calculated. When C(0) is finite, the variogram (h) is
bounded asymptotically by this value. The value of h at which the asymptot can be considered
achieved is called the range, while C(0) is called the sill (see Fig. 1.1).
1.3. KRIGING WITH THE INTRINSIC HYPOTHESIS 7

1.3.1 The variogram

The variogram is usually calculated from the experimental observations, and describes the
spatial structure of the RF. It can be shown [5] that given x1 , . . . , xn are n points belonging to
the domain of interpolation, for all the coefficients 1 , . . . , n R satisfing:
n
X
i = 0
i=1

then n X
n
X
i j (xi xj ) 0
i=1 j=0

i.e. (h) is a positive semidefinite function, and

(h)
lim =0
|h| h2

i.e. (h) is tends to infinity slowlier than |h|2 as |h| .

In principle there (h) could assume different behaviors also with the direction of vec-
tor h (anisotropy), but this is in general not easily verifiable due to the limited number
of data points usually available for hydrologic variables. If the experimental variogram dis-
plays anisotropy, then the intrinsic hypothesis is not verified and one has to use the so called
universal kriging [2, 3].
The most commonly used isotropic variograms are shown in Fig. 1.2 and are of the form:

1. polinomial variogram:
(h) = h 0<<2

2. exponential variogram:  
(h) = 1 eh

3. gaussian variogram: h i
2
(h) = 1 e(h)

4. spherical variogram: ( h i
1 3h


h 3
2


h
(h) =
h>

where and are real constants.

The variogram is estimated from the available observations in the following manner. The
data points are subdivided into a prefixed number of classes based on the distances between
the measurement locations. For each pair i and j of points and for each class calculate:

1. the number M of points that fall within the class;

8 CAPTER 1. GEOSTATISTICS IN HYDROLOGY

1
0.95
10
0.8
8
0.6
6 =1.5

0.4
=1.0
4

=0.5 0.2
2

1 2 3 4 5 1 2 3 4 5
3

1 1
0.95

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

1 2 3 4 0.2 0.4 0.6 0.8 1

2/3
3

Figure 1.2: Behavior of the most commmonly used variograms: polinomial (top left); exponen-
tial (top right); gaussian (bottom left); spherical (bottom right)
1.3. KRIGING WITH THE INTRINSIC HYPOTHESIS 9

2. the average distance of the class;

3. the half of the mean quadratic increment

1X
(Yi Yj )2 /M
2

In general the pairs are not uniformly distributed among the different classes as usually there are
more pairs for the smaller distances. Thus the experimental variogram will be less meaningful
as h increases. A best fit procedure together with visual inspection is then used to select the
most appropriate variogram and evaluate its optimal parameters.

Remark 1. An experimental variogram that is not bounded above (for example the polinomial
variogram with 1) implies an infinite variance, and thus the covariance does not exist. Only
the intrinsic hypothesis is acceptable. If the variogram achieves a sill then the phenomenon
has a finite variance and the covariance exists.

Remark 2. For every variogram, (0) = 0. However sometimes the data may display a
jump at the origin. This apparent discontinuity is called the nugget effect and can be due to
measurement errors or to microregionalization effects that are not evidenced at the scale of the
actual data points. If the nugget effect is present the variogram needs to be changed into:

(h) = + 0 (h)

where is the jump at the origin and 0 (h) the variogram without the jump. In Fig. 1.3 we
report a variogram with the nugget effect together with some other special cases that may be
encountered.

Remark 3. When we look for an estimator of the type

n
X

Y (x0 ) = i Yi
i=1

from the intrinsic hypothesis, i.e. eq. (1.6), we see that our random variable Y (x) must have a
constant mean m, so that:
Xn

E[Y (x0 )] = E[ i Yi ] = m
i=1

and, since E[Y (xi )] = E[Yi ] = m, we readily obtain:

n
X
i = 1
i=1
10 CAPTER 1. GEOSTATISTICS IN HYDROLOGY

(h) (h)

random structure
measured values
(no correlation)
nugget effect

h h

(h) (h)

nested structure periodic structure

h h

(h)

hole effect

h
Figure 1.3: Example of typical spatial correletion structures that may be encountered in ana-
lyzing measured data
1.4. REMARKS ABOUT KRIGING 11

The condition of minimum variance becomes now a constrained minimization problem which
can be solve by introducing a lagrange multiplier . The system (1.5) must be rewritten as:

0 (x1 x2 ) . . . (x1 xn ) 1 1 (x1 x0 )
.. .. .. ..
. .
. = .
(xn x1 ) 1 n (xn x0 )
1 1 ... 1 0 1

1.4 Remarks about Kriging

a. Kriging is a BLUE (Best Linear Unbiased Estimator) interpolator. In other word it is a
Linear estimator that matches the correct expected value of the population (Unbiased)
and that minimizes the variance of the observations (Best).

b. Kriging is an exact interpolator if no errors are present. In fact, if we set x0 = xi in (1.5)

we obtain immediately i = 1, j = 0, j = 1, . . . , n, j 6= i.

c. If we assume that the the error  is Gaussian, then we can associate to the estimate Y (x0 )
a confidence interval. For example, the 95% confidence interval is 20 where:
p
0 = var[Y (x0 ) Y (X0 )]

Then the kriging estimator (1.4) becomes:

n
X

Y (x0 ) = i Yi 20
i=1

d. The solution of the linear system does not depend on the observed value but only on xi
and x0 .

e. A map of the estimated regionalized variable, and possibly its confidence intervals, can
be obtained be defining a grid and solving the linear system for each point in the grid.

1.5 Kriging with uncertainties

We now assume that the observations Yi are affected by measurement errors i , and that:
1. the errors i have zero mean:

E[i ] = 0 i = 1, . . . , n

2. the errors are uncorrelated:

cov[i , j ] = 0 i 6= j
12 CAPTER 1. GEOSTATISTICS IN HYDROLOGY

3. the errors are not correlated with the RF:

cov[i , Yi ] = 0

4. the variance i2 of the errors is a known quantity and can vary from point to point.

The new coefficient matrix C of the linear system (1.5) is changed by adding to the main
diagonal the quantity i2 : 2
1 0
C := C 0 0
0 n2
and everything proceeds as in the standard case.

1.6 Validation of the interpolation model

The chosen model (in practice the variogram) can be validated by interpolating observed values.
If n observations Y (xi ), i = 1, . . . , n are available, the validation process proceeds as follows:

For each j, j = 1, . . . , n:
discard point (xj , Y (xj ));

estimate the Y (xj ) by solving the kriging system having set x0 = xj and using the
remaining points xi , i 6= j for the interpolation;

evaluate the estimation error j = Yj Yj ,

The model can be considered theoretically valid if the error distribution is approximately
gaussian with zero mean and unit variance (N (0, 1)), i.e. satisfies the following:

1. there is no bias: n
1X
i 0
n i=1

2. the estimation variance i is coherent with the error standard deviation:

n  2
1 X Yi Yi
1
n i=1 i

One can also look at the behavior of the interpolation error at each point looking at the mean
square error of the vector : v
u n
u1 X
Q=t 2
n i=1 i
1.7. COMPUTATIONAL ASPECTS 13

The uncertainties connected to the choice of the theoretical variogram from the experimental
data can be minimized by anaylizing the validation test. In fact, among all the possible vari-
ograms (h), that close to the origin display a slope compatible with the observations and gives
rise to a theoretically coherent model, one can choose the variogram with the smallest value of
Q.

1.7 Computational aspects

In the validation phase n linear systems of dimension n 1 need to be solved. The system
matrices are obtained by dropping one row and one column of the complete kriging matrix.
This can be efficiently accomplished by means of intersections of n 1-dimensional lines with
appropriate coordinate n-dimensional planes.
Note that the kriging matrix C is symmetric, and thus its eigenvalues i are real. However,
since n n
X X
i = Tr(C) = cii = 0
i=1 i=1

where Tr(C) is the trace of matrix C, it follows that some of the eigenvalues must be negative
and thus C is not positive definite. For this reason, the solution of the linear systems is usually
obtained by means of direct methods, such as Gaussian elimination or Choleski decomposition.
Full Pivoting is often necessary to maintain stability of the algorithm.

1.8 Kriging with moving neighborhoods

Generally the experimental variogram is most accurate for small values of h, with uncertainties
growing rapidly when h is large. The influence of this problem may be decreased by using
moving neighborhoods. With this variant, only the points that lie within a prefixed radius
R from point x0 are considered, provided that we are left with an adequate number of data
(Fig. 1.4). The radius R is selected so that the lag h will remain within the range of maximum
certainty for (h).
This approach leads also to high saving in the computational cost of the procedure be-
cause each linear system is now much smaller (Gaussian elimination has a computational cost
proportional to n3 , where n is the number of equations).

1.9 Detrending
In certain cases the observations display a definite trend that needs to be taken into account.
This is the case, for example, when one needs to interpolate piezometric heads observed from
wells in an aquifer system where a regional gradient is present. To remove the trend from the
data it is possible to work with residuals (= measurements - trend) that have a constant mean.
However this detrending procedure may be dangerous as it may introduce a bias in the results.
For this reason it is important that the trend be recognized not only from the raw data but
14 CAPTER 1. GEOSTATISTICS IN HYDROLOGY

R
x0

Figure 1.4: Example of interpolation with moving neighborhood

also from the physical behavior of the system from which the data come. If the trend cannot
be removed from the observations by simple subtraction, the universal kriging approach [2, 3]
can be used.
 1.10

1.10. EXAMPLE OF APPLICATION OF KRIGING FOR THE RECONSTRUCTION OF AN ELEVATION

Example of application of kriging for the recon-
struction of an elevation map
([DPSOHRI. U L J L Q J

70 52 75 50 75 45
7
60

75 60
6
75 88 102
78
80 70 104 55
5 71
66 75 73
90 80
64
4
60 95 90
60 48
3 80
55 80
105 70 54 70
2
70
50 66 51
1
60

90 45 65 40 55 25
0
0 1 2 3 4 5 6 7 8 9

The range of the observations is 25 - 104

 16
1100
Experimental Experimental
1100
variogram variogram
1000 10 1000
using 25 classes using 10 classes
900 900 20
5
800 800

Variogramma (m2)
Variogramma (m2)
8
700 700
600 600

CAPTER 1. GEOSTATISTICS IN HYDROLOGY

500 48 500
n. coppie 41 149
400 400 111
74 53 33 188
77 14
300 73 73 7959 300
47 180 30
43 13 116 80
200 78 61 200 177
72
39 46 21
100 100
25
0 0
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
Lag Distance (km) Lag Distance (km)

The shape is similar: robustness with respect to

different class subdivisions
 1.10 EXAMPLE OF APPLICATION
First try Estimation
7

Linear variogram = h variance 6

1.2
5 1.1
1
0.9
4
0.8
0.7
3
7 0.6
0.5
2
6 0.4
0.3
1 0.2
5
meaningless 0.1
4
0 0
0 1 2 3 4 5 6 7 8 9

3 0.00 -0.55 0.21 -0.63 0.07 0.00

7
-0.21

2 -0.00 -0.36
6
0.38 0.34 0.68
0.41
0.05 -0.07 0.37 -0.00
1
5 -1.17
-0.25 0.01 -0.23
0.02 -0.08
0 -0.00
4
0 1 2 3 4 5 6 7 8 9 -0.49 1.68 0.73
elevation 3
-0.24
0.40
-0.00
-0.48 0.54
0.17 0.04 -0.00 0.20
2
0.31
Kriging reconstruction 1
-0.27 0.21 -0.00
0.37
is robust evenwhen the 0.03 -0.53 0.25 -0.53 0.00 0.00
variogram is not accurate 0
0 1 2 3 4 5 6 7 8 9

interpolation error

17
 18
Second try

Linear variogram with nugget effect

nugget effect = 41.3

1100
1000
900

CAPTER 1. GEOSTATISTICS IN HYDROLOGY

800

(m2)
700
600

Variogram
500
400
300
200
100
0
0 1 2 3 4 5 6 7 8 9 10
Lag Distance (km)
 1.10 EXAMPLE OF APPLICATION
7

Estimation
variance 6

12
5 11
7
10
4
9
6
8
7
3
5 6
5
2
4 4
3
3
1 2
1
0 0
2 0 1 2 3 4 5 6 7 8 9

1 1.09 -5.87 4.86 -6.61 1.21 -4.20

7
-2.14
0 2.98 -3.71
0 1 2 3 4 5 6 7 8 9 6 3.83
2.70 7.82
1.92
-0.92
Elevation 5
0.94
-7.42
9.17 -2.59
-2.05 0.25 -1.79
1.34 1.22
-6.12

Reconstructed values
4
-4.04 10.57 7.34
-5.60 -3.29
are smoothed because of 3
-4.69 4.50
4.61

nugget effect: large 5.86

2
-0.03 -5.41 3.41
2.50
estimation variance -4.10 2.45 -3.19
1
3.90

1.27 -5.62 5.25 -6.32 3.51 -5.10

0
0 1 2 3 4 5 6 7 8 9

19
 20
Third try

Variogram without nugget effect and with large slope at the origin

polinomial variogram: = 136

= 0.5

900

800

CAPTER 1. GEOSTATISTICS IN HYDROLOGY

700

(m2)
600

500

Variogram
400

300

200

100

0
0 1 2 3 4 5 6 7 8 9 10
Lag Distance (km)
 1.10 EXAMPLE OF APPLICATION
7

Estimation 6

variance 12
5 11
7
10
4
9
6
8
7
3
5 6
5
2
4 4
3
3
1 2
1
0 0
2 0 1 2 3 4 5 6 7 8 9

1 -6E-005 -3 1 -3 0.6 -1E-005

7
-1
0 2E-006 -2
0 1 2 3 4 5 6 7 8 9 6
1 2 4
0.9
-0.4
Elevation 5
0.5
-1 0.2 -0.8
-3
2 -9E-005

0.6 0.5
-7E-007
4
Estimation variance is -2 5 3
-1 -5E-006
large far away from 3
-2 2
2

observation points 2
2 -0.08 3E-007 1
1
-2 0.9 5E-005
1

0.7 -3 1
2
-2 4E-005 -1E-005
0
0 1 2 3 4 5 6 7 8 9

21
 22
Fourth try

Gaussian variogram: small slope at the origin

1100 = 317
1000 = 2.49
900

CAPTER 1. GEOSTATISTICS IN HYDROLOGY

800

(m2)
700

600

Variogram
500

400

300

200

100

0
0 1 2 3 4 5 6 7 8 9 10
Lag Distance (km)
 1.10 EXAMPLE OF APPLICATION
Estimation
7

variance 6

12
7 5 11
10
6
9
4
8
7
5 3
6
5
4 2
4
3
3 1 2
1
2 0 0
0 1 2 3 4 5 6 7 8 9

1
-5E-009 -0.03 0.02 -0.2 -0.009 3E-010
7
-0.03
0
0 1 2 3 4 5 6 7 8 9 5E-005 -0.03
6
0.04 0.05 0.08
Elevation 0.02 0.006
0.04
0.03 7E-005
5 -0.08
0.02 -0.02 -0.07
small and smooth estimation -0.001 -0.02
-2E-005
variance; small errors 4
-0.03 0.1 -0.02
-0.02 -3E-005
3 0.04
-0.03 -0.01
0.004 -0.009 -8E-0050.02
2
0.03
-0.02 0.03 -0.0001
1

-0.02
0.01 -0.07 0.02 -0.042E-005 8E-008
0
0 1 2 3 4 5 6 7 8 9

23
24 CAPTER 1. GEOSTATISTICS IN HYDROLOGY
Bibliography

[1] de Marsily, G. (1984). Spatial variability of properties in porous media: A stochastic

approach. In J. Bear and M. Corapcioglu, editors, Fundemantals of Transport Phenomena
in Porous Media, pages 721769, Dordrecht. Martinus Nijhoff Publ.

[2] Delhomme, J. P. (1976). Applications de la th`eorie des varaibles r`egionalis`ees dans les
sciences de leau. Ph.D. thesis, Ecoles des Mines, Fontainebleaus (France).

[3] Delhomme, J. P. (1978). Kriging in the hydrosciences. Adv. Water Resources, 1, 251266.

[4] Matheron, G. (1969). Le krigeage universel. Technical Report 1, Paris School of Mines.
Cah. Cent. Morphol. Math., Fontainbleau.

[5] Matheron, G. (1971). The theory of regionalized variables and its applications. Technical
Report 5, Paris School of Mines. Cah. Cent. Morphol. Math., Fontainbleau.

[6] Volpi, G. and Gambolati, G. (1978). On the use of a main trend for the kriging technique
in hydrology. Adv. Water Resources, 1, 345349.

[7] Volpi, G., Gambolati, G., Carbognin, L., Gatto, P., and Mozzi, G. (1979). Groundwater
contour mapping in venice by stochastic interpolators. 2. results. Water Resour. Res., 15,
291297.

25