A computer program for interpolation of geohydrologic data using the method of "kriging"



             UNITED STATES DEPARTMENT OF THE INTERIOR
                        GEOLOGICAL SURVEY
                        RESTON, VA. 22092

In Reply Refer To:                            October 7, 1980
EGS-Mail Stop 411


GROUND WATER BRANCH TECHNICAL MEMORANDUM NO. 81.01


Subject:  A computer program for interpolation of geohydrologic
          data using the method of "kriging"

The purpose of this memo is to describe briefly the interpolation
procedure known as kriging, some of its applications and
limitations in hydrology, and the availability of a documented
computer program for kriging.  Although kriging techniques have
potential applicability to a variety of fields, we particularly
ask that this memo be brought to the attention of those persons
involved in ground-water studies.

Frequently there is a need to estimate the areal characteristics
of geohydrologic variables such as aquifer thickness,
transmissivity, storage coefficient, etc., continuously over a
region given the value of the variable at discrete points.  The
characteristics of the kriging method are that:  (1) correlation
information inferred from the data themselves are used in the
analysis; (2) the interpolation estimates are "statistically"
unbiased with minimum variance; and (3) such variance can be
computed and provides an assessment of the estimation error.  The
ability to give an assessment of the interpolating error is the
reason the method has attracted much attention.

Geohydrologic applications:

Kriging is applicable to regionalized variables which have some
spatial correlation structure (see next section).  Applications in
hydrology to date include estimating the areal distribution of:
(1) precipitation, by single event or annual mean; (2)
transmissivities; and (3) potentiometric head.  Examples of other
uses include estimating areal thickness of confining beds, ET
distribution in a basin, or spatial distribution of water-quality
parameters.  Delhomme (1978) gives several applications of kriging
to hydrology and describes the mathematical development.  He also
discusses the applicability of kriging to the analysis of data
networks.  Karlinger and Skrivan (1980) applied kriging to
estimation of the areal distribution of mean annual precipitation.

Development of the kriging system:

A kriging estimate (see original) of a "regionalized variable at
an arbitrary location is determined as a weighted sum of the
observed data.

          (see original for equation)

where the n observations z(subscript i) and their locations
(x(subscript i), y(subscript i)) are given by the hydrologist, and
the gamma(subscript i) terms are unknown weights to be determined.
The term "regionalized" is used to convey that the spatial
variation of the variable possesses a certain structure.  In other
words, that the variation is not completely random.

The kriging procedure determines the n weights, gamma(subscript
i), such that these conditions hold:  (1) the estimate is
unbiased--in the sense that the average error will be zero for a
large number of such estimates, i.e., (see original), where E[] is
the expectation; and (2) the variance of the estimation, i.e. (see
original) is minimal.  These conditions lead to a system of linear
equations involving the weighting terms gamma(subscript i).

An inherent part of the resulting equations is a function called
the semi-variogram which is a measure of the average difference of
data values for points a given distance apart.  The semi-variogram
obtained directly from the basic data is considered empirical.
Partly because of the varying degrees of uncertainty of the basic
data, the empirical semi-variogram often will not be a smooth
function.  An important part of the kriging analysis are the
inferences made from the empirical semi-variogram concerning:  (1)
the distance beyond which there is little or no correlation among
the data; and (2) the regularity or continuity of the data.  If
the semi-variogram indicates "essentially" no spatial correlation,
then kriging offers no advantage over other data fitting
procedures.  However, if a correlation structure is indicated,
several theoretically acceptable smooth semi-variogram models are
tried before selecting a best fit to the empirical variogram.

With a fitted semi-variogram, the linear system of equations is
solved for the unknown gamma(subscript i) values.  As a byproduct
of the solution, the variance of the kriging estimate, or kriging
error, is also calculated.  Interpolated values for the variable
can be computed at a grid of arbitrary locations for contouring of
estimates and kriging errors.  The errors can be used for defining
confidence intervals of the estimates.  For each point where an
interpolated value is needed a linear system of equations must be
solved for a corresponding set of gamma(subscript i) values.

Validation procedure:

A validation procedure is used to test the choice of the fitted
semi-variogram on the resulting interpolations.  This involves:

     (1)  The data for point (x(subscript 1), y(subscript 1)) is
          ignored.  Based on the data at the other (n-1) points a
          kriged estimate z*(subscript 1) is made at (x(subscript
          1), y (subscript 1)) and compared to z(subscript 1).
          The same process is followed for each of the data points
          being ignored one at a time.  In this way a residual
          (z*(subscript i)-z(subscript i)) is obtained for each of
          the n data points.

     (2)  The average residual for the n observations is
          calculated and compared to zero for unbiasedness.

     (3)  The squared residuals are compared to the kriging
          errors, which should approximately equal each other for
          theoretical consistency.

If these comparisons are favorable, then the fitted semi-variogram
is acceptable.  If not, then the semi-variogram model is adjusted.

Some practical considerations:

To define the semi-variogram, probably at least 20-30 data points
are necessary.  Because the inferences from the semi-variogram are
a significant part of the analysis, enough points are needed to
properly select a final form of semi-variogram model.  On the
other hand, because the size of the resulting system of linear
equations for validation is dependent on the number of
observations, it may not be practical to use more than 100-200
points in the validation procedure.  Cases with many more data
might use some subsets of the observational data for validation.

If there is an independent reason to expect a general trend (i.e.,
drift) in the variable, that trend can be made part of the
analysis.

Because a system of equations must be solved for a set of
gamma(subscript i) values for each points where an interpolation
is needed, it is practical to seek ways to reduce the number of
data points used to compute the kriged estimate.  In some cases it
is acceptable to consider only those points within a certain
neighborhood.

Limitations:

Gambolati and Volpi (1979) state that the kriged estimates are not
necessarily more reliable than those gotten from other methods of
estimating.  The primary appeal to the kriging technique is its
ability to provide some assessment of the estimation error.  The
criteria used to obtain kriging estimates are that the estimation
be unbiased and the variance of the estimation be minimized.
Satisfaction of those criteria does not assure that the resulting
estimates reproduce the spatial fluctuations of the true variable.
The spatial variability of the kriged estimate tends to be a
smooth version of the real variability.  The density of data
points probably influence how good a representation of the true
variability is possible with kriging.  The use of the kriged
estimates will influence the extent to which this "smoothing"
characteristic may limit their usefulness.  Hydrologists have only
very recently begun to develop experience with the kriging method.
It has some powerful features--but as all methods, does have
limitations.  Recognizing those limitations and being guided by
basic hydrologic principles, kriging should be a useful addition
to the set of tools used to interpolate between data points.

Computer Program K603:

The FORTRAN program K603 is documented in the enclosed Computer
Contribution (Skrivan and Karlinger, 1980) and is available on
SYS1.LOADLIB on both RE1 and RE3.  The documentation includes
details on data preparation and time estimates.  In addition, the
report has annotated bibliography on kriging development and
applications.  Additional copies of the report are available on
request from the Ground Water Branch.

References:

     Delhomme, J. P., Kriging in the hydrosciences:  Advances in
Water Resources, v. 1, no. 5, p. 251-266.

     Gambolati, G., and Volpi, G., 1979, A conceptual
deterministic analysis of the kriging technique in hydrology:
Water Resources Research, v. 15, no. 3, p. 625-629.

     Karlinger, M. R., and Skrivan, J. A., 1980, Kriging analysis
of mean annual precipitation, Powder River Basin, Montana and
Wyoming:  U.S. Geological Survey, Water-Resources Investigations
(Approved May, 1980).

     Skrivan, J. A., and Karlinger, M. R., 1980, Semi-variogram
estimation and universal kriging program; U. S. Geological Survey
Computer Contribution, 98 p.



                              (s) Charles A. Appel
                              for Eugene P. Patten
                              Acting Chief, Ground Water Branch

Enclosure

WRD Distribution:

A (Memo only)
B (NR, 12 copies w/enclosures)
  (SR,  6 copies w/enclosures)
  (CR, 21 copies w/enclosures)
  (WR, 12 copies w/enclosures)
S (Memo only)
Copy w/enclosure to District Offices
Copy (Memo only) to Subdistrict Offices