Comparison of spatial prediction methods for mapping floodplain soil pollution

I

i CATENA

vol. 17, p. 535 550

Cremlingen 1990

C O M P A R I S O N OF SPATIAL P R E D I C T I O N M E T H O D S FOR MAPPING F L O O D P L A I N SOIL P O L L U T I O N H. Leenaers, J.R Okx, RA. Burrough Utrecht Summary

I

The floodplain soils of the River Geul (The Netherlands) are polluted by heavy metals, i.e. lead, zinc and cadmium, Several spatial prediction methods (local trend analysis, mathematical splines, inverse distance weighting, block kriging, point kriging and point co-kriging) were used to map zinc levels in a small and intensively sampled study area. Three subsets of the set of 145 data points, containing 12, 23 and 44 data respectively, were used to estimate zinc levels at 99 test locations. Correlation coefficients of linear relations between observed and estimated zinc levels and contour maps of prediction errors indicate that weighted local averaging methods perform better than the other methods. In the study area, point co-kriging with elevation as a co-variable outperforms all other methods when only few data on zinc are available.

ISSN 034l 8162

@1990 by CATENA VERLAG,

D-3302 Cremlingen-Destedt, W. Germany 0341-8162/90/5011851/US$ 2.00 + 0.25

('A[ENA

A n lnterdiscipfinar~, J o u r n a l o f S O I L S C I E N C E

HYDROLOGY

Introduction

Due to past mining activities, the floodplain soils of the Geul river are polluted with heavy metals, particularly zinc, lead and cadmium ( L E E N A E R S et al. 1988, R A N G et al. 1986). The general pollution pattern consists of a logarithmic decay with distance to the source of contaminants ( L E E N A E R S et al. 1988), with local deviations that arc', caused by variations of flood frequency and of sedimentary conditions during flood events. The pollutants constrain the land use in these areas, therefore detailed maps are required that delineate zones with high concentration levels. Successful attempts have been made by W O L F E N D E N & L E W I N (1977) and R A N G et al. (1987) to relate the pollution level of floodplain soils to floodplain characteristics such as geomorphology, inundation frequency and soil type. These studies led to the production of choropleth maps that delineated broad pollution zones, but did not provide detailed information about the continuous spatial variation of pollution levels within the zones. There are many point interpolation techniques that can provide this type o f information, see for example LAM

GEOMORPHOI.O(J5

536

Leenaers, O k x & Burrough

(1983) and B U R R O U G H (1986). However, to date only few attempts have been made to compare spatial prediction methods quantitatively (BREGT et al. 1987, DAVIS 1976, LASLETT et al. 1987, K U I L E N B U R G et al. 1982). Obviously, one seeks an interpolation technique that gives the best results for a given investment in observations. The aim of this study is to evaluate some of these techniques (local trend analysis, mathematical splines, inverse distance weighting, point kriging, block kriging and point co-kriging) for the purpose of mapping floodplain soil pollution. The quality of the interpolation techniques used is based on the correlation between predicted and observed values at 99 test locations and on the spatial distribution of prediction errors.

2

The study area

Significant occurrences of metal ore are found near Plombi6res and Kelmis, both situated in the Belgian part of the Geul basin (fig. 1). The exploitation of zinc and lead ores in the Geul basin probably began in the thirteenth century. Mining activities were highest from 1820 to 1880. During that period, maximum yearly ore production in the open-cast mines of Kelmis was about 135,000 tonnes, and 1,300 men were employed in the zinc industries. The last mine closed in 1938, but the processing of metal ores continued until the 1950s. Separation techniques exploiting the differences in specific density were used when mining the ores. These techniques were inefficient and resulted in high concentrations of ore particles and metalrich spoil in the effluent, which was discharged directly into the river. The reject material and tailings were dumped in (ATINA

large heaps. Some of these heaps still exist. The presence of these point sources has led to the dispersal of heavy metals through the catchment during and after the active mining period. During high flow stages flooding occurs and the metal-rich suspended sediments are deposited on the floodplains (LEENAERS et al. 1988). Due to erosion processes during periods of rainfall, the spoil heaps are still active sources of heavy metals. Erosion of locally contaminated stream banks is an additional source of heavy metals. 3 3.1

prediction methods

Spatial

Fitting of local trend surfaces

The simplest way to describe gradual variations of a property is to model them by polynomial regression. The idea is to fit a polynomial surface by least squares through the data points, where it is assumed that the spatial coordinates X and Y are independent variables, and that Z, the property of interest, i,; the dependent variable. In two dimensions the polynomials are surfaces of the form: Z(X, Y) =

2

(br~ " X ~ " Y~)

r+s<=p

where r and s are the polynomial coefficients, p is the order of the trend surface and b~ is a constant. A thorough explanation of this technique is given by DAVIS (1986), RIPLEY (1981), CHORLEY & H A G E T T (1965) and WATSON (1971). Two important disadvantages of this general trend analysis, i.e. the susceptibility to outliers and the inability to fit a low order polynomial through complex data ( B U R R O U G H 1986), can be reduced with a moving window approach: only the data that are within a

Al~lnterdisciplinaryJournalc~[SOlI S ( I E N ( ' [

HYDROIfI(~Y

(}FOMORPtI()IO(IY

Spatial Prediction Methods £or Mapping Soil Poflution

,•[•

t" " - "L.

..I J

537

-,,

.,.

1£~

_ Valkenburg( , _ -..

I._ ~--.

I

Schin ~{o Geul

Gulpen •

/

f'e~len~ )L.~ FRG

t

i

I

"~x+ # ÷'t 'r ~a"~ , ~ u g g " e l ~ ÷+ "1"%~, ' - % \ A I '

Plombi~re~sox-~

f--/

~+

] ~ r" LUlUlW_~-/"~r/

Kelmis

% 2 I

4 I

6 I

8 I

.~-I-+++ +./. + ++ ~. ~--~.~'~

\-

\

0 ;

K

1Okra I

, %,...

~..%

~, ~,/

/

I

Fig. 1" The River Geul catchment. Fitting of mathematical splines

specified search radius are used to estimate the property of interest by fitting

3.2

a local trend surface. For our study we fitted first order local trend surfaces by using the nearest points within a search radius, which was derived from an analysis of the spatial structure (see section on variogram models),

Given a set of datapoints, one of the simplest mathematical expressions for a continuous surface that intersects these points is an interpolating polynomial of the lowest order that passes through all datapoints (LAM 1983, R I P L E Y 1981, D U B R U L E 1983, T I P P E R 1979). The major deficiency of this polynomial fit is that since the polynomial is entirely

(J,\l E N A

An Interdisciplinar~ Journal of S O I L S £ ' I E N C E

HYDROLOGk

-GEOMORPHO[

OGY

Leenaers, Okx & Burrough

538

unconstrained, except at the datapoints, the values attained between the points may be drastically different from those at nearby data points (LAM 1983). This problem may be alleviated to a certain extent by employing piecewise polynomial surfaces to cover the area. For our study we used a local algorithm that exploits only the four nearest data points to estimate Z at an unvisited site. 3.3

Inverse distance weighted interpolation

The principle of distance weighting methods (WEAVER 1964, S H E P A R D 1968, M c L A I N 1976) is to assign more weight to nearby points than to distant points, which is based on the idea that observations located close together tend to be more alike than observations spaced further apart. The most common form of the weighting function is the reciprocal function 1/d 2, where d is the distance between a nearby point and the point to be estimated:

~'/_lZ(xi) .do2 n

Z(Xj) =

Zi=I du 2

where the x; is the point at which the surface is to be interpolated. For d = 0 the exact value of the original sampling point has to be preserved, otherwise the function goes to infinity. By changing the value of the inverse power parameter, the distribution of weights can be varied from being highly biased in favour of near data points to giving nearly equal weight of all data points, There are several disadvantages to weighting methods. First, the choice of a weighting function may introduce ambiguity, especially when the characteristics of the underlying surface are not known. Second, the weighting method is easily

CAIENA

affected by uneven distribution of data points since unequal weight will be assigned to each of the points. Finally, because this method is by definition a smoothing technique, maxima and minima in the interpolated surface can occur only at data points. The location and magnitude of extreme values therefore, cannot be detected when they are not included as original sample points (LAM 1983). However, the simplicity of the principle, the speed in calculation, the ease of programming and reasonable results for many types of applications have led to wide adoption of this method, as well as various improvements (LAM 1983).

3.4

Point kriging, block kriging and point co-kriging

Recent publications in soil science have demonstrated the application of regionalized variable theory to mapping soil variables (BURGESS & WEBSTER 1980a, 1980b, B U R R O U G H 1986, DAVIS 1986, WEBSTER 1985). The theory provides a convenient means of summarizing soil spatial variability in the form of a semi-variogram which can be used to estimate weights for interpolating the value of a given soil property at an unsampled location. Kriging is a form of weighted local averaging that is optimal in the sense that it provides estimates of values at unsampled sites without bias and with minimum variance which is estimable. The weighting procedure is similar to that used in inverse distance weighted interpolation except that the weights are derived from a geostatistical spatial analysis based on the sample semivariogram. Point kriging or punctual kriging is an exact interpolator in the sense that

An Interdisciplinary

Journal o[ SOIL S ( I E N ( E

IIYDROI f)GY

(H OMt')RPHOI O(]~

Spatial Prediction Methods for Mapping Soil Pollution

Ee l vao tin

539

~~

?O'o~

Zn ic

M

%

lOOK ~

~ ~~

Fig. 2 a) Digital Elevation Model of the study area (x- and y-coordinates and elevation in m) b) Block diagram of Zn levels (mg/kg) in the study area.

I'¢.1 L N A

A n I n t e r d i s d p l i n a r ? J o u r n a l o f S O I L SCIENC.'E

I~IYDROLO(3Y

GEOMORPHOLO(~

Leenaers, Okx & Burrough

540

at sampled locations, interpolated and measured values coincide. Point kriging implies that all interpolated values relate to an area or volume that is equivalent to the area or volume of an original sample, Given the often large, short-range nature of soil variation, point kriging sometimes results in maps that have many sharp spikes or pits at the datapoints. If average values of soil properties over areas are considered more interesting than values at points, block kriging is the alternative (BURGESS et al. 1980b). The estimation variances obtained for block krigingare usually substantially less than those for point kriging. The results are much smoother surfaces free from the pits and spikes resulting from point kriging. Co-kriging is a logical extension of kriging to situations where two or more variables are spatially interdependent and the one of immediate interest is undersampled ( J O U R N A L & H U I J B R E G T S 1978, M c B R A T N E Y & WEBSTER 1983, M Y E R S 1982, 1984, L E E N A E R S et al. 1989). Co-kriging may be a useful method for interpolating a property that is expensive to measure by making use of a spatial correlation between the property of interest and an attribute that is less expensive to measure.

4

Experimental and analytical methods

4.1

Sampling scheme and laboratory analyses

A regular grid of 5x29 (grid size 5x5 m) was laid out in the floodplain area of the River Geul, near the village of Valkenburg (fig. 1); 100 g-samples of top soil material (depth: 0~10 cm) were collected

('A'IENA

at 145 locations. Elevation relative to mean sea level was determined for each sample location by using a level. All samples were dried for 24 hours at 60°C, and crushed with a mortar. Two gram of this material was gently boiled for 2 hours with 20 ml 30'}/0 HNO3. The extract was then separated from the sediment by centrifugation and brought to 40 ml with distilled water. The concentrations of zinc was determined by direct flame absorption spectometry. Block diagrams of elevation and zinc levels are shown in fig. 2a and b. It is clear that zinc levels are inversely related to elevation, which can be explained by the fact that low-lying areas have a higher risk of inundation. It will be shown that benefit may be drawn from the relation between zinc levels and elevation by means of a co-kriging procedure.

4.2

The test procedure

The set of 145 observations on topsoil zinc levels was divided into two subsets. The larger set contained 99 data points (and one missing value) and was set aside for future validation. The smaller set of 44 data (and one missing value) was used for the estimation of experimental semi-variograms of zinc levels and elevation and a cross semi-variogram of both. The same set was used to design three configurations of sample locations with different densities, containing 12, 23 and 44 observations respectively (fig. 3). For each configuration, zinc levels were estimated using the interpolation methods discussed above for locations in the original sampling grid. For cokriging, the data on zinc levels were supplemented by all elevation data (n=145). Correlation coefficients of linear reid-

An Interdisciplinary Journal nlSOl[

S( [ E N ( |

HYI)ROI,()G'{

(}EOMllRPHOI

O(iY

Spatial Prediction Methods for Mapping Soil Pollution 10

20

,

~o

20

-~

,

• I

,

*

-,-'L

10

,,~

30

|

~

20

~

40 ,

® ,

~L

30

50 ,

,

(~) I

~

4-0

60 ,

-,-J~

50

80

90

100

1 I0

120

130

140

"~"

® [

70

541

* ,

~

® I

60

~L

70

<~) I

~

80

• I

-,-'L

90

* ,

~

1O0

® [

T"L

110

~) I

~

• I

120

-l-'L

130

0 |

140

set of 12 samples O

set of 23 samples set of 44 samples

Fig. 3: Sampling grids in the study area.

attribute(s) Zn elevation Zn-elevation

Co 4588 417 -697

C

C0+C

106694 111282 2657 3074 -15992 -16689

a

a'

17.1 23.8 20.7

29.6 40.8 35.8

Co: nugget; C0+C: sill, a': practical range. Tab. 1 : Parameters o.1 the theoretical Gaussian variogram Jhnctions.

tions between e s t i m a t e d and observed zinc levels were used to evaluate the precision o f the estimates and to c o m p a r e the p r e d i c t i o n methods. In addition, contour m a p s o f the a b s o l u t e p r e d i c t i o n errots were p r o d u c e d in o r d e r to j u d g e their spatial distribution,

4.3

Estimation of the semi-variogram models

In general, a few simple features contribute to the form o f a semi-variogram, The semivariance at zero lag m u s t be zero, b u t in practice the e x t r a p o l a t e d s e m i v a r i o g r a m usually intercepts the ord i n a t e at a positive value k n o w n as the 'nugget variance'. The nugget variance can arise from m e a s u r e m e n t error, discrete r a n d o m v a r i a t i o n a n d spatially d e p e n d e n t v a r i a t i o n occurring over dis-

('ATIINA An Interdisciplinary Journal of SOIL S('IEN(!E

tances much less than the s a m p l i n g interval ( O L I V E R & W E B S T E R 1986). In most instances, the semivariance increases from the nugget variance to a m a x i m u m level b e y o n d which it remains level. Such s e m i - v a r i o g r a m s are transitive, the m a x i m u m value is k n o w n as the sill, which is the a priori variance o f the variable. The lag distance at which the sill is reached is k n o w n as the range. As an o b s e r v a t i o n o f a variable can d e p e n d in the statistical sense on observations o f the same variable at other locations nearby, they can also be spatially related to o b s e r v a t i o n s o f o t h e r variables. Such variables are co-regionalized. A n a l o g o u s to the single variable, the d e p e n d e n c e between two variables can be expressed by a cross s e m i - v a r i o g r a m ( M c B R A T N E Y & W E B S T E R 1983). In this study, elevation was chosen as a co-

HYDROI O(}Y GEOMORPHOI ()Gk

Leenaers, Okx & Burrough

542

7(h) 150000

8 o

D

100000

~

D

C

50000 /

0

i 0

1 10

~

I

i

20

I

i

30

I

I

40

50

tag,h (m) 7(h) 0

b

-5000

-10000

-150OO

-20000

0

I

I 10

I

I I i 20 30 IGg,h (m)

I

I 40

I

50

Fig. 4 a) semi-variogram for Zn (semi-variance in mg/kg 2) b) cross-variogram for Zn and elevation (cross-variance in m.mg/kg).

CATENA

An Interdisciplina D .Iournal of S O l [ ~;(][gNCE

HYDROLC)(]Y

(][iOMORPIIOLO(;'r

Spatial Prediction Methods For Mapping Soil Pollution

543

spatial prediction method

number of input data n - 12 n=23 n=44

first order local trend surfaces mathematical splines inverse distance weighting block kriging point kriging point co-kriging

0.41 0.10 0.77 0.73 0.76 0.85

0.57 0.61 0.75 0.79 0.80 0.89

0.63 0.77 0.83 0.84 0.90 0.89

Tab. 2: Correlation coefficients o['linear relations between estimated and observed Zn

levels (n=99) for sets of input data. variable for estimating zinc levels. Most samples were located in the natural levee area where only small variations in elevation existed; 19 observations were in a a b a n d o n e d channel loop. Multidirectional semi-variograms for zinc (fig. 4a) and elevation were cornputed for the set o f 44 data. Due to the shape of the 5x29 grid, the results are dominated by effects of the spatial variation parallel to the longer axis. The parameters required for kriging were obtained by fitting a gaussian equation:

7(h'C°'C'a)=C°+C(1-exp(-h2/a2)' where 7 h

=

semi-variance; lag (m); nugget variance; sill;

Co co + c a = distance parameter, through the estimates o f experimental semivariance. The practical range (a') was determined to be the lag at which 95% o f the sill is reached. The cross semi-variogram for topsoil zinc and telative elevation (fig. 4b) could also be modelled adequately using a gaussian model. The cross semi-variances are negative because of the inverse relationship between the variables. In the process of cross semi-variogram fitting the CauchySchwarz inequality:

(',",lENA

A n Interdisciphnar'~ J o u r n a l o f S O I L S C I E N C E

HYDROLOGY

712(h) < = V/(Tl(h) * 72(h)) where ]"12 --

cross

for all h > = 0

semi-variance;

~,~ - auto semi-variance zinc; )'2 - auto semi-variance elevation, was checked to guarantee a positive cokriging variance under all circumstances ( J O U R N E L & H U I J B R E G T S 1978, M Y E R S 1982, M Y E R S , 1984). The parameter o f the fitted semi-variogram models are listed in tab. 1. During the block kriging procedure a block size of 20×20 m was used. For the other prediction methods, the parameters were set in such a manner that for each estimate only data fall within the range (a') o f the semivariogram o f zinc are used: for inverse distance weighting and the local trend m e t h o d the search radius was set equal to the range (a') and for the fitting o f mathematical splines the number' o f nearest points was set to 4.

5 5.1

Results and discussion Goodness-of-fit of the modelled zinc content surfaces

Fig. 5 shows the contours of predicted zinc levels based on interpolation from 44 data points. As might be expected from its approximate nature, the fitting

GEOM(JRPHOIOGY

Leenaers, Okx & Burrough

544

of first order local trend surfaces yields a rather smooth zinc surface, whose generalized nature is most striking in the vicinity of the abandoned channel where steep gradients should prevail (fig. 2). Visual examination of the contours suggest that by fitting mathematical splines and by local averaging methods (inverse distance weighting and kriging)more accurate results are obtained, Correlation coefficients of the linear relations between observed and predicted zinc levels were used as a measure of goodness-of-fit (tab. 2). For all methods the correlation coefficients increases with an increasing number of input data. The rate of increase however, is notably different. The increase for the local trend method occurs gradually (from r=0.41 for n=12 to r=0.63 for n=44), while the performance of the splines method changes from a very poor (r=0.10 for n=12) to a relatively high level (r=0.77 for n=44). This pattern illustrates the different qualities of smoothing, approximate methods and flexible, exact prediction methods, When only 12 data points are available, these two methods are outperformed by the weighted local averaging methods (inverse distance weighting and kriging). The correlation coefficient is remarkably high for co-kriging, which illustrates the advantage of using additional, spatially correlated, information when only few data points are available on the variable of interest. The relatire benefit of co-kriging decreases as the number of data points with known zinc levels increases: for n=44 point kriging and point co-kriging perform equally well (r=0.90 and r=0.89 resp.). With these data points, the correlation coefficients for inverse distance weighting and block kriging are 0.83 and 0.84 respec-

CATENA

tively. For block kriging, this can be explained by the fact that this method smoothesextreme values over a 2 0 x 2 0 m area. The relatively small improvements by inverse distance weighting when the number of data points was increased however, are probably due to the choice of the weighting function. Since this function does not relate to the spatial correlation structure of the data, its applicability shows an undesirable dependence on the input data density and/or data geometry. The fact that the correlation coefficient decreased from 0.77 to 0.75 when the number of input data points was increased from 12 to 23 illustrates this behaviour. 5.2

The spatial distribution of estimation errors

For environmental assessment studies not only the goodness-of-fit of the spatial model but also some idea about the magnitude of the prediction errors is of importance. Instead of using measures of central tendency and dispersion, we decided to produce contour maps of prediction errors. For each interpolation technique, the grids containing the estimates were subtracted from the grid containing the 145 data at the data points, and the absolute value of the difference at each grid intersection was computed. The resulting grids have the same dimensions as the sampling grid and contour maps could be produced. As for the approximate methods, errors that occur at the input data locations are included in the maps. By this means, not only the magnitude but also the spatial distribution of prediction errors can be studied. The occurrence of larger errors should preferably be restricted to small, well-defined areas.

A n I n t e r d i s c i p l i n a r y J o u r n a l o1" S O I l N ( ' I E N ( ' E

HYDROIO(~Y

(HOMORPHO[O(~'I

Spatial Prediction Methods for Mapping Soil Pollution o

20

~o

20

. . . .

~o 0

~ I

I

~o

i

,

~o

t

,

I

I

~o

70

. . . . .

~o

90

~t

,,~,/

,oo

,

~

~io

,

,~o

,~

~o

t\,

,

~o

\~,i2°

I

I

~o

i

i

i

I

/~

I

i

i

i

o

I

"L.._.~ -.-y

.o 7 . . . . . . . . . .

I

I

,,~

I

i

i

i

i \

i%~i

i

l/

i

i'rx

i

7/

J,o

i

I

~/

.

I

.

I

.

.._,...W__t_ #

I

.

I

.

!

o

.

.

'

'

yt't'77t,.~

\ \ \

t._.../

\

"

\ \

i

e>

,~D/"- ~

I

'

"~ "'°<'

fT'',.,

K.___~Y

0

°

,o

~°

,

oo~,~o

,o \/%F~-__/~

0

~o

,/

545

o

.I'° 0

,o

~

o

I I'll. . . .

~o

10

10

o

,

~..'Y

1o 0

i .-I---I.

,

,

,

,

,,~--Y-

, /T",,\\~I

'~

\'

'I

~' \ t'

'

'

'

o

'

20

~---"~J m

~J

,

P

o 0

,

10

I

I

,

, 20

1o

30

i

i

,

, ~

4.0

.¢'1'~'~ i

50

,

~

? 60

I

I

i

i

°

, , , , , 70

\

._

,\",,,\\'~

,I/,I~7. 80

90

1 O0

110

,\,

120

'

I

I

0

, 130

Fig. 5: Contours of Zn levels (mg/kg) obtained from 44 data by a).first order local trend analysis; b) mathematical splines; c) inverse distance weighting; d) block kriging; e) point kriging; f ) point co-kriging; and g) contours of Zn levels based on all 143 data..

140

Leenaers, Okx & Burrough

546 20

0

10

20

50

40

50

60

7(3

80

90

100

110

120

1,30

140 20

10

10

0

0

20

20

100

100

ioo

10

0

20

~

o

20

::::::::::::::::::::::::::::::::

..............

~o

0

0

'°

........

!ii i,i~i ~i:ii1o'0

20 ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ..:~:::::::::::

o

K • ::::::::::::::::::::::

o

0

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ....::

~0

...........~::::::'::::: 20

30

20

:

::::::::::::::::::::::

::::::::

::::::::::::::::::::::: 40 ~0

60

o

:::::::::::::: ::::::::::::::

70 o-loo

80

~00

~0

~20

~0

o ~,,0

r~o/kg

~

100--300

i~

300- 500 rag/kg > 500 mg/kg

m

90

mg/kg

Fig. 6: Contours of prediction errors (mg/kg Zn) by a)first order local trend analysis," b) mathematical shines; c) inverse distance weighting; d) block kriging ; e) point kriging; f ) point co-kriging..

( A I ENA An [nterd sciplinary Journ

ol SOIl. S('IEN('E

HYDROI O(;Y

(;E(IMORI'IIOLO(IY

Spatial Prediction Methods for Mapping Soil Pollution

However, one should be aware of the fact that approximate methods (block kriging and fitting of local trend surfaces), due to their nature of generalization, are likely to produce a pattern of regularly distributed errors, Fig. 6 shows contour maps of absolute errors (>100 mg/kg) that were obtained by interpolation from 44 input data points. With an exception for block kriging, all methods yield large errors only in those areas where high zinc levels prevail. The errors produced by fitting local trend surfaces are highest and cover a relatively large area. The error maps obtained by inverse distance weighting and fitting of mathematical splines are much alike: the area that is covered by the occurrence of errors >100 mg/kg are about equally large and in the abandoned channel a number (3M) of errors >500 mg/kg occur. For block kriging the occurrence of errors >500 mg/kg is restricted to one grid cell only, but smaller errors of 100-300 mg/kg can be observed in the entire area. Point kriging and point co-kriging combine both advantages: a relatively small area is covered by errors in the ranges 100 300 and 300500 mg/kg and in only one grid cell an error of >500 mg/kg occurred, 5.3

Diseussion

The performance of spatial prediction methods will likely depend on the type of data that are used and on the spatial correlation structure of the data. It is therefore impossible to draw general conclusions from the results that were obtained in this case study. Nevertheless, this study illustrates some of the qualities of the methods, as will be discussed below, The surface ofzinclevels that was used

( NI tNA

An Interdisciplinary Journal of SOIL SCIENCE

HYDRO[.OGY

547

as the basis for this study shows large short range variations (fig. 2). It is clear that approximate methods, that do not aim at preserving the original data nor at exact prediction at unvisited sites, will yield a smooth zinc surface that ignores extreme values. Therefore, fitting of local trend surfaces yields poor results and is not recommended for this type of data. The main purpose of using this method was to illustrate the difference with more appropriate methods. Fitting of mathematical splines, contrary to trend surfaces, has the advantage that very steep gradients can easily be modelled. However, it lacks the possibility to attach more weight to nearby points. By inverse distance weighting, the problem of attaching equal weights to all data within the search radius is to a certain extent alleviated. However, since the weights are derived from some convenient function and do not relate to the spatial correlation structure of the data, the results depend on the often unknown applicability of the chosen function. In this study the function lid 2 provided excellent results with 12 data points, where the search radius was set equal to the range of the semi-variogram of zinc. However, compared to point kriging relatively small improvements occurred when the number of data points was enlarged, which illustrates the inability of this function to optimally exploit additional information. By choosing a few other weighting functions we found that with the set of 44 data points weighting by 1/d 4 was sligthly more appropriate: the correlation between predicted and observed values increased from 0.83 to 0.86. It is clear that in practice one will not always have the possibility to check the validity of the weighting function. The semi-variogram model of zinc

GEOMORPHOLOGY

Leenaers, Okx & Burrough

548

provide weights for interpolation by point kriging and block kriging. The advantage of this sample-based weighting function is best illustrated by the results of point kriging. Both in terms of goodness-of-fit and the spatial distribution of prediction errors this method performs better than any of the other methods that do not use additional elevation data. The smoothing nature of block kriging is nicely illustrated by the fact that the area covered by small errors (10(~300 mg/kg)is relatively large, The only method that clearly shows major advantages over all other methods is point co-kriging. Its capability to benefit from inexpensive additional information by making use of a strong spatial correlation between zinc levels and elew~tion, accounts for the fact that with 12 input data only an adequate spatial model of zinc levels could be obtained, Taking advantage of this spatial correlation during the design of a sampling strategy, a substantial reduction of costs can be realized. In this study for example, the cost for sampling and laboratory analysis is c. ECU 25,- per sample, and the collection of elevation data costs c. ECU 1,30 per location. Thus, if one neglects the difference in required computer time, for point kriging with 44 data on zinc an amount of c. ECU 1150,is needed, while point co-kriging with 145 elevation data and 23 data on zinc costs ECU 790,-. Without loss of precision --- the two methods yield the same goodness-of-fit (tab. 2) - - a cost reduction of more than 30% could be brought about.

timates of zinc levels when many data points are available in areas with large short range variations. Because of their capability to account for spatial dependence, local weighted averaging methods (inverse distance weighting and kriging) require less data to provide an acceptable estimate of the spatial distribution of zinc levels. However, in this case study the geostatistically based sample semi-variogram provides better weights for interpolating zinc levels than inverse distance weighting functions. Major improvement in the spatial prediction of" zinc levels in the Geul floodplains is provided by point co-kriging with elevation as a co-variable. Inexpensive elevation data exhibit a strong spatial correlation with zinc and contain information on the dispersal mechanism of zinc. The results presented demonstrate that the cokriging method can be greatly improved when there is a physical understanding of the contamination process. Moreover, without loss of precision a cost reduction of 30% could be realized because co-kriging requires less data points. Acknowledgements This study was financed by the Netherlands Organization for Scientific Research (N.E.O.). We wish to thank R Nienhuis (Free University of Amsterdam) and J. van Keulen for geostatistical advice and software support. References

BREGT, A.K., BOUMA, J. & JELLINEK, M.

6

(1987): Comparison of thematic maps derived from soil map and from kriging of point data. Geoderma 39, 281-291.

Conclusion

In the Geul floodplains, the use of mathematical splines only provides reliable es('A'IENA

BURGESS, T.M. & WEBSTER, R. (1980a): Optimal interpolation and isarithmic mapping of

An Interdisciplinary Journal o l S O I l

SCIEN{'L:

HYDR{}L(,(;Y

(}hOMORI'H(IL()GY

Spatial Prediction Methods for Mapping Soil Pollution

soil properties. 1. The semi-variogram and punctual kriging. Journal of Soil Science 31, 315 3 3 1 . BURGESS, T.M. & WEBSTER, R. (1980b): Opritual interpolation and isarithmic mapping of soil propertes, lI. block kriging. Journal of Soil Sience 31,333 341. BURROUGH, P.A. (1986): Principles of Geographical Information Systems for Land Resources Assessment. Monographs on soils and resources survey no. 12. Clarendon Press, Oxt'ord, 193 p.

British Geographers, Publ.

no.

DAVIS, J.C. (1976): Contouring algorithms. In: AUTOCARTO II, Proc. Int. Symp. on Computer-assisted Cartography, Sept. 1 9 7 5 , 352-359. DAVIS, J.C. (1986): Statistics and Data Analysis in Geology. John Wiley and Sons, New York. 646 p. DUBRULE, O. (1983): Two Methods with Different Objectives: Splines and Kriging. Mathematical Geology 15, 2, 245-257. JOURNEL, A.G. & HUIBREGTS (1978): Mining Geostatistics. Academic Press, New York, 600 p. KUILENBURG, J. VAN, DE GRUIJTER, J.J., MARSMAN, B.A. & BOUMA, J. (1982): Accuracy of spatial interpolation between point data on soil moisture supply capacity, compared with estimates from mapping units. Geoderma 27, 311 325.

MeBRATNEY, A.B. & WEBSTER, R. (1983): Optimal interpolation and isarithmic mapping of soil properties. V, Co-regionalization and multi-sampling strategy. Journal of Soil Science 34, 137 162. MeLAIN, D.H. (1976): Two dimensional interpolation from random data. Computer Journal 19, 178 181. MYERS, D.E. (1982): Matrix formulation of cokriging. Mathematical Geology 14, 3. 249 257. MYERS, D.E. (1984): Cokriging: new developments. In: Verly, G., David, M., Journel, A.G. & Marechel, A. (eds.), Geostatistics for Natural Resonrces Characterization (part l). NAT()

CHORLEY, R.J. & HAGETT, P. (1965): TrendsurI'ace mapping in geographical research. Trans. Inst. 37, 47 67.

549

ASI Series, Series C: Mathematical and Physical Sciences 122, 295 305. D. Reidel Publishing Compny, Dordrecht. OLIVER, M.A. & WEBSTER, R. (1986): Semivariograms for modelling the spatial pattern of landform and soil properties. Earth Surlace Processes and Landforms 11,491 504. RANG, M.C., KLEIJN, C.E. & SCHOUTEN, C.J. (1986): Historical changes in the enrichment of fluvial deposits with heavy metals. IAHS Publication 157, 47-59. RANG, M.C., KLEIJN, C.E. & SCHOUTEN, C.J. (1987): Mapping of soil pollution by application of classical geomorphological and pedological field techniques. In: Gardiner, V. (ed.), International geomorphology, part 1, 1029 1044. Wiley & Sons, New York. RIPLEY, B. (1981): Spatial Statistics. Wiley, New York, 252 p.

LAM, N.S. (1983): Spatial interpolation methods: a review. The American Cartographer 10,

SHEPARD, D. (1968): A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 23rd National Conference ACM, 517 524.

2, 129 149. LASLETT, G.M., McBRATNEY, A,B., PAHL, P.J. & HUTCHINSON, M.F. (1987): Cornparison of several spatial prediction methods for soil pH. Journal of Soil Science 38, 325 341. LEENAERS, H., RANG, M.C. & SCHOUTEN, C.J. (1988): Variability of the metal content of flood deposits. Environmental Geology and Water Science I1, 1, 95-106. LEENAERS, H., BURROUGH, P.A. & OKX, J.P. (1989): Efficient mapping of heavy metal pollution on floodplains by co-kriging from elevation data. Chapter 4 in: Raper, J., Three dimensional applications in Geographical Informarion Systems. Taylor & Francis, London.

TIPPER, J.C. (1979): Surface modelling techniques. Kansas Geological Survey Series on Spatial Analysis 4. 108 p. WATSON, G.S. (1971): Trend-surface analysis. Journal of the International Association of Mathematical Geology 3, 3, 215 221. WEAVER, R.C. (1964): Relative merits of interpolation and approximate functions in the grade prediction problem. In: Parks, G.E. (ed.), Computers in the Mineral Industries. Stanford University Publications in the Geological Sience 9, 171-185. WEBSTER, R. (1985): Quantitative Spatial Analysis of Soil in the Field. Springer-Verlag, New York, 69 p.

I ~1 L N A

a, Ea Intcrdisciphnar~ Journal of SOIL S ( I [ N( I

HYDROt()(}Y

GEO[vl()Rt~t[()l.()(i~

Leenaers, Okx & Burrough

550

W O L F E N D E N , P.J. & LEWIN, J. (1977): Distributions of metal pollutants in floodplain sediments. C A T E N A 4, 309 317.

Addresses of authors: H. Leenaers J.P. Okx* P.A. Burrough Department of Physical Geography University of Utrecht p.o. Boxc 80.115 3508 TC Utrecht The Netherlands * present address: T A U W Infra Consult B.V. P.O. Box 479 7400 AL Deventer The Netherlands

CA T [N A

An Interdisciplinary Journal of SOil S C I [ N C E

HYDROI.O(IY

(]EOMORPIIOIO(;Y