A sampling scheme for estimating the mean extractable phosphorus concentration of fields for environmental regulation

Geoderma 89 Ž1999. 129–148

A sampling scheme for estimating the mean extractable phosphorus concentration of fields for environmental regulation D.J. Brus a

a,)

b , L.E.E.M. Spatjens , J.J. de Gruijter ¨

a

DLO Winand Staring Centre for Integrated Land, Soil and Water Research, P.O. Box 125, 6700 AC Wageningen, Netherlands b Laboratory for Soil and Crop Testing, P.O. Box 115, 6860 AC Oosterbeek, Netherlands Received 19 January 1998; accepted 31 August 1998

Abstract A soil sampling scheme for estimating the mean extractable P concentration of fields is designed to be used as a tool for environmental regulation of the application rates of manure. The field to be sampled, is split up into geographically compact blocks of equal area that are used as strata. From each stratum one sampling point is selected by Simple Random Sampling. These samples are bulked into one composite for the field. The geographical stratification is performed by restricted least-squares clustering of raster cells using the coordinates of the midpoints as classification variables and the within-group sum of squares as the minimisation criterion. Using a variance model and a cost model, the numbers of sample points and laboratory analyses are optimised simultaneously, given a maximum allowed variance of the total error Žsampling error plus measurement error.. To predict the sampling variance, variograms have been estimated for 16 fields differing in land-use, soil parent material and phosphate level. A pooled relative variogram was used to predict the sampling variance for various sample sizes Ž5 to 50., field-areas Ž1 to 10 ha. and phosphate levels Žfor grassland 20 to 80 mg P2 O5 extracted in ammonimum lactate per 100 g soil, for arable land 20 to 80 mg P2 O5 extracted in water per 1 dm3 soil.. The cost model consists of three components: Ži. fieldwork cost; Žii. field equipment cost and, Žiii. laboratory cost. For the 16 fields, the predicted sampling variance of the Stratified Sampling design is 0.8 to 0.4 times the predicted variance of Simple Random Sampling if 40 points were sampled. To estimate the mean extractable P concentration with a total variance F 9, replicate measurement of the composite only pays if the mean extractable P concentration of the field exceeds 40 to 50. This

)

Corresponding author. Tel.: q31-371-474-237.

0016-7061r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 6 - 7 0 6 1 Ž 9 8 . 0 0 1 2 3 - 2

130

D.J. Brus et al.r Geoderma 89 (1999) 129–148

critical phosphate level increases with the maximum allowed variance of the total error. q 1999 Elsevier Science B.V. All rights reserved. Keywords: soil testing; Stratified Random Sampling; geostatistics; cluster analysis; optimisation; composite samples

1. Introduction In the Netherlands, phosphorus disposal in areas with intensive animal husbandry has become a public concern because of water-quality issues. To prevent further degradation of the environment, the application of manure is regulated by law at present. Although at the moment, the maximum allowed application rates are not related to the phosphate level of the fields, modification of the regulation in this direction is being considered. In that case, it is important to include a sampling protocol in the regulation so that discussions on the quality of the soil test can be avoided. In the current sampling scheme for soil testing in the Netherlands, used to advise farmers on application rates of fertiliser from an agricultural production point of view, the sampling locations are selected on the basis of their assumed representativeness and convenience ŽFerrari and Vermeulen, 1955. . For environmental regulation purposes this procedure is inappropriate because selection bias is not excluded and because points are selected more-or-less ‘subjectively’. The aim of this study is to design an objective, unbiased and efficient soil sampling scheme appropriate to test statistically hypotheses on the mean extractable phosphorus concentration Žhereafter shortened to mean P concentration. of any agricultural field in the Netherlands, and to predict the precision and cost of the new sampling scheme. The cost of the new scheme must be in the same order of magnitude as that of the current scheme Žapproximately US$37 per field. . Special attention is paid to the optimisation of the numbers of sample points and laboratory analyses. We will use the phrase ‘soil sampling scheme’ for a combination of a sampling strategy Žthat is itself a combination of a sampling design and an estimator., sampling method, compositing scheme and measurement method. The paper is organised as follows. The new soil sampling scheme is described in Section 2.1. In this scheme, the field under study is stratified in geographically compact strata. In Section 2.2, we describe how this is done. The error in the estimated mean P concentration of a given field can be split up into the sampling error and the measurement error. To predict the sampling error we need prior information on the spatial variation modelled as a variogram. To obtain this information, 16 fields of various land-use, soil parent material and phosphate level were intensively sampled Ž Section 2.3. . In Section 2.4, we describe the models used to predict the sampling error and the measurement

D.J. Brus et al.r Geoderma 89 (1999) 129–148

131

error. The cost of sampling was predicted for the situation that a Differential Global Positioning System is used to map the boundary of the field and to navigate the fieldworker to randomly selected sample points within it Ž Section 2.5.. Section 2.6 describes how we optimised the sample size and the number of analyses. In Section 3 we present the predicted variance and cost for various sample sizes, field-areas and mean P concentration, as well as the optimal numbers of sample points and analyses.

2. Materials and methods 2.1. The new soil sampling scheme We propose to select samples randomly, i.e., by probability sampling, because this sampling procedure eliminates selection bias and selection subjectivity ŽBrus and De Gruijter, 1997. . In probability sampling points can be selected with various sampling techniques such as Simple Random Sampling, Stratified Sampling, Cluster Sampling, Systematic Sampling and so on Žsee for instance Cochran, 1977; Sarndal et al., 1992. . The question is which of these techniques ¨ is most appropriate. In general, given the sample size, the precision can be increased by spreading the sample points over the field. This will also increase somewhat the time to travel Žwalk or drive by moped. to the sample points, and therefore the cost of fieldwork. However, the contribution of travel time to the total cost is relatively small Žsee Section 2.5.. An obvious sampling design to reach a good coverage is Systematic Sampling with the sample points on a regular grid. However, a problem with this sampling design is that in general the number of sample points that falls within the field is not constant but depends on the randomly selected starting location. It is hard to sell to employers of soil tests that lot has decided that for instance only 32 points will be sampled if on average 40 points are sampled. In the proposed sampling scheme, a good coverage of the randomly selected points is achieved by splitting up the field to be sampled into blocks of equal area, and selecting from each block one point by Simple Random Sampling. A similar approach is mentioned by Black Ž1992.. In sampling terms, this sampling design type is referred to as Stratified Ž Simple. Random Sampling. To ensure that the cost of the new scheme more or less equals the current cost, we propose to bulk the samples at points Ž cores. into one composite sample. In general, the composite sample mean is an unbiased estimator of the field mean only if the blocks have equal area. Note that the sampling variance cannot be estimated from the data obtained by the new soil sampling scheme. To test hypotheses on the mean P concentration, we propose to use the sampling variance as predicted by a validated model.

132

D.J. Brus et al.r Geoderma 89 (1999) 129–148

Fields can be split up into blocks of equal area in many ways. For instance, an elongated rectangle with width w can be split up into n rectangular strips with a width wrn. If these strips were used as strata, sample points still may cluster, which would give a low precision. To avoid this, we require that the strata are geographically as compact as possible, i.e., they approximate circles as far as possible. In Section 2.2 we describe the method we used to form such strata. The sampling method Žsampling device, sampling depth. of the new scheme equals that of the current scheme. The sampling device and sampling depth differs for grassland and arable land, resulting into a different support. For grassland, the support is a cylinder with a diameter of 23 mm and a height Žsampling depth. of 5 cm, for arable land this is a cylinder with a diameter of 13 mm and a height of 25 cm. In the laboratory, the composite sample is subsampled in two steps. In the first step, 300 g is subsampled. The second step differs between grassland and arable land. For grassland, 2.5 g is subsampled, for arable land and maize 1.2 ml is subsampled. To analyse the subsamples, the phosphorus in samples from grassland is extracted in ammonium lactate Ž Egner ´ et al., 1960., that in samples from arable land and maize in water Ž Sissingh, 1971.. In the first method, known as the P-AL method, the P concentration is expressed as mg P2 O5 per 100 g soil, in the latter method, the Pw method, as mg P2 O5 per dm3 soil. 2.2. A method for geographical stratification To avoid bias in the estimates from the composite sample, we require that the strata are of equal size, and to minimise the sampling variance, we want them to be as compact as possible. The rationale for this is suggested by the theory developed by Matern ´ Ž1986. for stationary isotropic random fields. The withinstratum variance of compact strata is generally low. We used the following heuristic method to construct compact strata of equal size. The field is represented by a large number, say N, raster cells. Then these raster cells are clustered by the method of least squares Ž k-means, Hartigan, 1975., using the co-ordinates of the midpoints as classifying variables, under the condition that the clusters have the same number of raster cells. This method starts with a random partition of the set of all N points into L Žnoncontiguous. clusters of equal size, if L is the required number of strata. Of course, in general, N will not be a multiple of L. Therefore, the clusters differ in size by at most one raster cell, which can be neglected because the number of raster cells in the clusters is large. This starting solution is then iteratively improved to minimise the objective function: the sum of the squared distances from the points to the centres of the clusters to which they are allocated. Improvement takes place by re-allocating the points to the cluster with the nearest centre, and recalculating the centres. To ensure that the clusters keep the same size, the

D.J. Brus et al.r Geoderma 89 (1999) 129–148

133

re-allocations are limited to swopping two points between clusters. From this it follows that, apart from small discontinuities due to the rasterization, the clusters which are the sampling strata are least-squares Thiessen polygons of equal area. This method does not guarantee that a global minimum of the objective function is reached. To check this, several random starting points were used, and also the method was applied to problems with a known best solution. In these tests, local minima were rarely found and these were nearly as good as the global ones. Fig. 1 shows the result for Wesepe Ž 2.1 ha, N s 5278 raster cells of 2 = 2 m. and Oude Tonge Ž2.9 ha, N s 7278 raster cells of 2 = 2 m.. In both cases L s 40 strata were constructed. 2.3. Data on spatial Õariation The sampling error of the new scheme for a given field can be estimated by repeated sampling of this field according to this scheme. A drawback of this method is that, in this way, we obtain information about this scheme only. We do not know the precision for other sampling design classes, or even for other sample sizes with the same sampling design class. Therefore, we predicted sampling variances by the method described by Domburg et al. Ž 1994. . This method makes use of a model of spatial variation, more precisely a variogram. Prior information on the spatial variation of the study variable was scarce, so we decided to sample 16 fields intensively. Table 1 shows some properties of the fields. The fields form eight combinations of land-use and parent material. Of each combination, a field with a relatively high mean P concentration and a field with a relatively low mean P concentration was sampled. From each field, samples were taken at 52 points forming a rectangular grid covering the entire field Žsmall fields. or several cross-shaped clusters Ž large fields. . The 52 samples were analysed separately by taking two subsamples and analysing both subsamples as described earlier. Experimental variograms were estimated for the 16 fields, using the mean values of the two subsamples Žduplo-means. . These experimental variograms were then corrected for the measurement error. As all 52 samples of a field were subsampled and analysed twice, we have 52 estimates of the variance of a single measurement per field. This variance divided by two is an estimate of the variance of the duplo-mean. The estimated variance of the duplo-mean was strongly correlated with the duplo-mean itself. We fitted a simple linear model for Pw and for P-AL. These models were used to predict the variance of the measurement error for the 16 fields, using the mean of the 52 duplo-means as a predictor value. This predicted variance was subtracted from the experimental semi-variances. Next, we fitted omnidirectional models to the experimental variograms by iterative regression, using the number of pairs divided by the fitted value in the

134

D.J. Brus et al.r Geoderma 89 (1999) 129–148

Fig. 1. Stratification of Wesepe Ža. and Oude Tonge Žb. into 40 compact blocks of equal area.

D.J. Brus et al.r Geoderma 89 (1999) 129–148

135

Table 1 Some properties of the 16 fields sampled for variogram estimation Field name

Land-use

Parent material

Area Mean extractable P Žha. concentration Žmg P2 O5 per 100 g soil Žgrassland. or per dm3 soil Žarable land..

Wesepe Lemelerveld Oude Tonge Den Bommel Exloermond 1 ¨ Exloermond 2 ¨ Dreumel 12 Dreumel 13 Oldelamer 5 Oldelamer 11 Oostrum Rotteveel Oostrum De 5 Haaren Achter Broek Haaren Bremakker Nuth Slenaken

Maize Maize Arable land Arable land Arable land Arable land Grassland Grassland Grassland Grassland Grassland Grassland Grassland Grassland Grassland Grassland

Aeolian sand Aeolian sand Marine clay Marine clay Humose sandqcut-over peat Humose sandqcut-over peat Fluvial clay Fluvial clay Peat Peat Marine clay Marine clay Aeolian sand Aeolian sand Loess Loess

2.1 3.8 2.9 4.4 2.0 2.0 2.4 1.3 3.2 2.8 1.8 1.7 3.4 1.0 0.9 1.7

23 137 16 38 90 38 40 44 22 55 18 65 24 64 15 46

P extracted in ammonium lactate Žgrassland. or in water Žarable land..

previous fit as a weight ŽCressie, 1985. . Only spherical, exponential and linear models were tried out. The best model was selected on the basis of the residual mean sum of squares. Fig. 2 and Table 2 show the result. Note that for some fields ŽLemelerveld, Slenaken, Dreumel 12, Dreumel 13. the estimated semivariances for the largest distances were not used in the fitting because this would lead to a systematic error in the fitted values for the short distances. For estimating the mean semi-variances of the strata, we only need the first part of the variogram because the strata are rather small. If we compare the two variograms of a combination, the semi-variance at a given distance is generally larger for the field with the high mean P concentration. The semi-variance seems to be related to the mean value. Therefore, we fitted relative variograms Ž Isaaks and Srivastava, 1989. . The experimental variograms were scaled as follows:

gˆ R Ž h . s

gˆ Ž h . m ˆ Ž h.

b

Ž1.

where gˆ RŽ h. is the experimental relative semi-variance at lag h, gˆ Ž h. is the experimental semi-variance at lag h, m ˆ Ž h. is the mean of all data that contribute

136

D.J. Brus et al.r Geoderma 89 (1999) 129–148

Fig. 2. Experimental variograms and fitted models for the 16 fields.

to lag h, and b is the power. The parameter b was fitted simultaneously with the remaining parameters of the model. For six combinations, the fitting algorithm converged. Table 3 and Fig. 3 show the result. The parameter values of the bounded variograms are of the same order of magnitude, except for the grassland on peat combination, which has a much larger value for b Ž2.205. and a shorter range Ž11.95 m. . Also, for the linear relative variograms, the parameter values are more or less similar. The samples from grassland and arable land have slightly different support, therefore, ceteris paribus, a different variogram can be expected for these land-use types. However, the effect of the support on the variogram seems to be overruled by other, more dominant effects. Therefore, we finally fitted a relative variogram for all the fields Ž Table 3. . To this end, we dropped the experimental

D.J. Brus et al.r Geoderma 89 (1999) 129–148

137

Table 2 Function type and estimated parameters of fitted variograms Žcorrected for measurement error. for 16 fields Field

Function type

Nugget

Sill Žslope.

Range parameter Žm.

Wesepe Lemelerveld Oude Tonge Den Bommel Exloermond 1 ¨ Exloermond 2 ¨ Dreumel 12 Dreumel 13 Oldelamer 5 Oldelamer 11 Oostrum Rotteveel Oostrum De 5 Haaren Achter B. Haaren Brem. Nuth Slenaken

Linear Linear Spherical Spherical Exponential Exponential Spherical Spherical Spherical Spherical Spherical Spherical Spherical Spherical Linear Linear

17.7 118.1 31.6 13.3 105.1 10.7 94.1 53.0 36.4 – 71.1 112.9 61.1 82.6 63.1 139.2

0.492 7.98 62.6 96.3 465.5 228.1 284.9 183.3 72.6 499.7 105.4 246.6 252.0 167.2 0.875 6.61

– – 72.7 64.6 39.9 22.5 69.4 35.8 30.7 12.2 73.2 60.4 38.2 176.7 – –

–: Not fitted.

variogram of Lemelerveld because the experimental semi-variances were clearly outliers, even after dividing them by their mean raised to power 1.061. Note that, also, the mean P concentration of Lemelerveld is a clear outlier Ž Table 1. . 2.4. Prediction of precision 2.4.1. Sampling error The variance of the total error V Ž e tot . is equal to the sum of the variance of the sampling error V Ž e fld . and the variance of the measurement error V Ž e lab .: V Ž e tot . s V Ž e f ld . q V Ž e lab .

Ž2.

For the stratified sampling strategy described in Section 2.1, the sampling variance can be predicted by Ž Domburg et al., 1994. : L

L

Vˆ Ž e f ld . s Ýw l2 Vˆl Ž e f ld . s Ýw l2 l

l

gl nl

1

s

L

Ýg l L2

Ž3.

l

where VˆŽ e fld . is the predicted sampling variance, L is the number of strata, w l is the weight Žrelative area. of stratum l, Vˆl Ž e fld . is the predicted sampling variance of stratum l, g l is the mean semi-variance between all pairs of points in stratum l, and n l is the number of points in stratum l. The mean semi-variance between all pairs of points in a stratum, is equal to the spatial variance

138

D.J. Brus et al.r Geoderma 89 (1999) 129–148

Table 3 Function type and estimated parameters of relative variograms for six combinations of land-use and soil parent material and for all fields Land-use

Soil parent material

Function type

Nugget

Sill or slope

Range parameter Žm.

Power

Maize Arable land Arable land

Aeolian sand Marine clay Humose sandq cut-over peat Peat Marine clay Loess

Linear Spherical Exponential

0.709 5.79 3.25

0.0394 18.56 16.23

– 69.0 37.2

1.061 0.449 0.744

Spherical Exponential Linear Spherical

– 7.82 1.954 4.1

0.073 16.63 0.0322 11.4

11.95 26.1 – 68.4

2.205 0.652 1.237 0.830

Grassland Grassland Grassland All fields

within this stratum averaged over all realisations of the model Ž Domburg et al., 1994.. The mean semi-variance of a stratum depends on the variogram and the geometry Žsize and shape. of the stratum. 2.4.2. Measurement error Previous experiments of the Laboratory for Soil and Crop Testing, a private company founded in 1928 by a group of Dutch agricultural organisations, showed that the standard deviation of the measurement error for a single measurement can be predicted by the linear model Ž V Ž e lab .. 1r2 s 0.707 q 0.039 P, where P is the mean P concentration.

Fig. 3. Relative variogram for the arable land–humose sand combination Žleft. and arable land–marine clay combination Žright..

D.J. Brus et al.r Geoderma 89 (1999) 129–148

139

The variance of the mean of several measurements Ž on different subsamples. is simply the variance of a single measurement divided by the number of analyses. The measurement error can be split up into the subsampling error and the error of the chemical analysis. We assumed that the subsampling error is independent of the size of the composite sample for the sizes considered in this study. 2.5. Prediction of cost The total cost Ž expressed in US$. can be split up into the cost of fieldwork Ž Cf ld ., the cost of the equipment for fieldwork Ž Ceqp . and the cost of the measurements in the laboratory Ž C lab .: Ctot s Cf ld q Ceqp q C lab

Ž4.

The cost of the fieldwork and of the field equipment can be related to the time needed for fieldwork: Cf ld s c f ld t Ceqp s ceqp t

Ž5.

where c fld and ceqp are the cost of fieldwork per hour and the cost of field equipment per hour, respectively ŽUS$ per hour., and t is the time needed for fieldwork Ž h. . The cost of fieldwork per hour is largely determined by the salaries of the personnel. In 1997, a fieldworker of the Laboratory for Soil and Crop testing cost US$45.00 per hour. The cost of the equipment is largely determined by the cost of a Differential Global Positioning System ŽDGPS. and a field computer. In 1997, it cost US$21.25 per hour to rent this equipment. The time for fieldwork can be split up into Ži. time for drawing a random sample using the sampling design described above Ž tdrw ., access time, i.e., time to travel Žwalk. to the sample points Ž tacc . , the time for sampling Ž tsam ., and the time for travelling to the field and administration Ž t tra .: t s tdrw q tacc q tsam q t tra

Ž6.

To draw a random sample, first a sampling frame has to be constructed. We assumed that this is done by mapping the field boundary in reality, i.e., by walking around the field with a DGPS. The result is a GIS-file in vector format, which must be converted to raster format. Then a compact geographical stratification must be calculated, and lastly, from each stratum one point Ž raster cell. must be drawn. The time for mapping the field boundary can be calculated by the ratio of the perimeter Žkm. and the speed of walking Žkmrh.. To predict the cost of

140

D.J. Brus et al.r Geoderma 89 (1999) 129–148

fieldwork for fields with unknown perimeter but known area, we calculated the ratio between the perimeter in km and the area in ha for the 16 fields. The mean of this ratio was 0.348. We assumed that the speed of walking is 5 kmrh. The time for calculating the strata depends on the hardware, the software Že.g., optimisation algorithm, language., and the number of raster cells. Given the hardware and the software, computing time increases with number of raster cells. On the other hand, the more raster cells we have for a field, the smaller the rasterising error ŽBurrough, 1986. . It seems reasonable to tune the size of the raster cells to the area to obtain a more or less constant relative rasterising error. We assumed that it takes 3 min to calculate the strata. The time for drawing a sampling point from each stratum is negligible. The access time depends on the distance to walk. The average distance between two neighbouring points for the proposed sampling design with 1 point per stratum can be approximated by the square root of the area of a stratum. As the strata have equal area, this equals Ž ArL. 1r2 , where A is the area of the field in ha. To get the total distance, this average distance is multiplied by the number of points: Ž ArL. 1r2 n s Ž Arn. 1r2 n s Ž An. 1r2 . To obtain the access time this distance is divided by the speed of walking between the points. We assumed that this is 4 kmrh, which is slightly less than when digitising because it costs some time to locate the points. The time to collect the samples at the randomly selected points can be calculated as the product of the number of points and the time per point, tsam . The time per point was set to 10 s. The time for fieldwork should be raised by the time of travelling to the fields, to collect name and address of owners, etc. Ž t tra .. For the current scheme, the average time for travelling and administration, etc., is 10 min. We used this figure to predict the cost of the new scheme. The time for travelling and administration, and the time for drawing a sample Ž tdrw . are independent of the sampling scheme and therefore do not influence the optimisation result. We need estimates of these time components only to predict the total cost. The cost of the measurements in the laboratory can be split up into the cost of the sample preparation plus one analysis Ž Cpre ., and the cost of extra analyses: C lab s Cpre q Ž r y 1 . cana

Ž7.

where r is the number of analyses Žreplicates. and cana are the cost per extra analysis. The cost of the first analysis is included in the cost of sample preparation because this analysis is always done, whereas an extra measurement is incidentally only. This makes that the cost per analysis is different for the first and following measurements. In 1997, it cost US$18.50 to prepare and analyse a sample at the Laboratory for Soil and Crop Testing. The cost of replicates is US$5 for each extra analysis. Finally, from a practical point of view, it is relevant to split up the cost into the once-only and repetitive cost. The field needs to be digitised and stratified

D.J. Brus et al.r Geoderma 89 (1999) 129–148

141

only once, whereas it must be sampled and the sample must be analysed every year we want to estimate the mean. 2.6. Optimisation of sample size and number of analyses We optimised the sampling scheme by searching for the combination of number of sample points and number of analyses that minimises the total cost for a given constraint on the precision. This was done by calculating for a range of values of the mean P concentration and field-area, the minimum number of sample points that gives a predicted total variance equal to or less than the maximum allowed variance if the composite is analysed once, if analysed twice and so on. Next, the cost of the sampling scheme was predicted, given the minimum sample size and the number of analyses, and the cheapest scheme was selected.

3. Results and discussion 3.1. Predicted precision We predicted the precision for the 16 fields for 40 points which is the current sample size, using the field-specific variograms and the pooled relative variogram ŽTable 4.. For most fields, the predicted standard errors are rather close

Table 4 Standard deviation of predicted sampling error for 40 points using the field-specific variogram and the pooled relative variogram Field

Field variogram

Pooled relative variogram

Wesepe Lemelerveld Oude Tonge Den Bommel Exloermond 1 ¨ Exloermond 2 ¨ Dreumel 12 Dreumel 13 Oldelamer 5 Oldelamer 11 Oostrum Rotteveel Oostrum De 5 Haaren Achter Broek Haaren Bremakker Nuth Slenaken

0.8 2.5 1.0 1.1 2.2 1.5 1.9 1.6 1.2 3.3 1.4 1.9 2.0 1.5 1.3 2.3

1.4 – 1.2 1.9 2.5 1.8 1.8 1.8 1.4 2.1 1.3 2.1 1.5 2.1 1.1 1.9

142

D.J. Brus et al.r Geoderma 89 (1999) 129–148

Fig. 4. Predicted spatial variance within strata Žleft. and predicted standard deviation of sampling error Žright. as a function of the number of strata Žsample size., for Dreumel 12 and Oldelamer 5.

which implies that the pooled relative variogram can be used safely to predict the sampling variance for these fields. Exceptions are Oldelamer 11 Ž standard error with pooled relative variogramy standard error with field-specific variogram s y1.2., Den Bommel Ž0.8., Wesepe Ž0.6. and Haaren Bremakker Ž 0.6. . The large difference for Oldelamer 11 can be explained by the low value for the estimated power of the pooled relative variogram Ž0.830, see Table 3. compared to the estimated power for the grassland–peat combination Ž 2.203. . For Den Bommel, Wesepe, and Haaren Bremakker the field-specific variograms show more spatial structure than the pooled relative variogram, which may partly explain the overestimation with the pooled relative variogram. We also predicted the standard errors for sample sizes 20, 25, 30 and 35. Fig. 4 shows the result for Oldelamer 5 and Dreumel 12. For the latter field, the predicted sampling variance decreases much stronger with sample size than for the former field. This can be explained by the shape of the variogram. The variogram of Oldelamer 5 is close to a pure nugget variogram Ž see Fig. 2. , and consequently, the predicted spatial variance within the strata is almost constant with the area of the strata ŽFig. 4., whereas for Dreumel 12, this spatial variance clearly increases with the area due to the spatial structure as shown by the variogram. To predict the sampling variance, the geometry of the field Ž strata. must be known. However, if we assume that the strata are squares with an area equal to the area of the field divided by the number of strata, the error in the predicted standard error is very small, especially if the number of strata is large ŽTable 5. . Table 5 Error in predicted standard deviation of sampling error due to the assumption of square strata Number of strata Mean error

20 y0.0029

25 y0.0028

30 y0.0021

35 y0.0012

40 y0.0004

D.J. Brus et al.r Geoderma 89 (1999) 129–148

143

This assumption leads to a slight underestimation, which means that the squares are a bit more compact than the Thiessen polygons. The relative error Ž mean

Fig. 5. Predicted standard deviation of sampling error as a function of the area and the sample size, for fields with a mean extractable P concentration of 30 Ža. and 60 Žb..

144

D.J. Brus et al.r Geoderma 89 (1999) 129–148

error divided by the standard deviation calculated with the Thiessen polygons as strata. is smaller than 1% in all cases and therefore can be neglected. Consequently, it is possible to predict quite precisely the sampling variance for fields without knowing its geometry. The only thing we must know is its area. We predicted the standard deviation of the sampling error, using the pooled relative variogram, for fields with an area ranging from 1 ha to 10 ha, a sample size ranging from 5 to 50, and a mean P concentration ranging from 20 to 80. Fig. 5 shows the results for a mean P concentration of 30 and 60. Both figures show that, given the sample size, the standard deviation of the sampling error increases with the area, especially when small fields are intensively sampled Župper left parts of the figures. . This result contradicts the conclusion of Ferrari and Vermeulen Ž1955. that the ‘‘the size of the field has, on the average, no effect on the magnitude of the sampling error’’. 3.2. Predicted cost The once-only cost Ž Conc . for digitising and stratifying the field can be expressed by the simple formula Conc s 3.313 q 4.611 A. Fig. 6 shows the repetitive cost of sampling and laboratory analysis for fields as a function of the area Ž1–10 ha. and sample size Ž5–50. if the composite is

Fig. 6. Predicted repetitive cost in US$ as a function of the area and the sample size, for one laboratory analysis on composite.

D.J. Brus et al.r Geoderma 89 (1999) 129–148

145

analysed once. For every extra measurement, US$5 should be added to the values in Fig. 6. 3.3. Stratification effect We compared the predicted sampling variance of the stratified design with the predicted sampling variance from Simple Random Sampling, i.e., without stratification. We expect a gain in precision due to the stratification because the stratification precludes clustering of sample points. The sampling variance for Simple Random Sampling can be predicted by inserting L s 1 into Eq. Ž 3. . We used the relative precision Ž Webster and Oliver, 1990, p. 44. to quantify the stratification effect. The result for five fields is shown in Fig. 7. It shows that the Stratified Sampling design is always more precise than the Simple Random design. The gain in precision differs strongly between the fields, which is caused by the differences in the spatial dependency: the stronger the dependency, the more redundancy we have in case of clustering of sample points, the more gain

Fig. 7. Relative precision expressed as the ratio of the variance for Simple Random Sampling ŽSI. and the variance for the proposed Stratified Simple Random Sampling design ŽSTSI.. For values )1, STSI is more precise than SI.

D.J. Brus et al.r Geoderma 89 (1999) 129–148

146

Table 6 Optimal sample size and number of analyses for a maximum allowed total variance of 9, for field-areas ranging from 1–10 ha and mean extractable P concentration ranging from 20–80; 12r1s12 sample points and 1 analysis Area P-cont.

1

2

3

4

5

6

7

8

9

10

20 30 40 50 60 70 80

12r1 19r1 31r1 27r2 38r2 37r3 49r3

13r1 21r1 34r1 30r2 41r2 40r3 53r3

14r1 22r1 24r2 32r2 43r2 43r3 56r3

15r1 23r1 25r2 33r2 45r2 44r3 58r3

15r1 24r1 26r2 34r2 36r3 46r3 48r4

16r1 25r1 26r2 35r2 37r3 47r3 49r4

16r1 25r1 27r2 36r2 38r3 48r3 51r4

16r1 26r1 28r2 37r2 39r3 49r3 52r4

17r1 26r1 28r2 37r2 40r3 50r3 53r4

17r1 27r1 29r2 38r2 41r3 51r3 54r4

can be expected. For all fields, the gain in precision increases with the sample size. 3.4. Optimal sample size and number of analyses Tables 6 and 7 show the optimal schemes for a maximum variance of 9 and 16. The optimal sample size increases slightly with the field-area Ž given the mean P concentration. and increases strongly with the mean P concentration Žgiven the field-area.. The downward jumps in the optimal sample size coincide with upward jumps in the optimal number of analyses. This optimal number of analyses is hardly influenced by the field-area. Replication of the analysis only pays if the mean P concentration exceeds a critical level. Below this critical level it is more efficient to sample additional points. The critical level increases with the maximum allowed variance. For a maximum allowed variance of 9, it is 40 Žfield-area F 2 ha. or 50 Žfield-area ) 2 ha., for a variance of 16, it is 60 Žfield-area F 3 ha. or 70 Žfield-area ) 3 ha. . For a maximum variance of 25,

Table 7 Optimal sample size and number of analyses for a maximum allowed variance of 16, for field-areas ranging from 1–10 ha and mean extractable P concentration ranging from 20–80; 7r1s 7 sample points and 1 analysis Area P-cont.

1

2

3

4

5

6

7

8

9

10

20 30 40 50 60 70 80

7r1 10r1 13r1 18r1 26r1 20r2 26r2

8r1 11r1 14r1 20r1 28r1 22r2 28r2

8r1 11r1 15r1 21r1 30r1 24r2 30r2

8r1 12r1 16r1 22r1 20r2 25r2 31r2

9r1 12r1 17r1 23r1 21r2 26r2 32r2

9r1 13r1 17r1 23r1 22r2 27r2 33r2

9r1 13r1 18r1 24r1 22r2 27r2 34r2

9r1 13r1 18r1 25r1 23r2 28r2 35r2

9r1 13r1 18r1 25r1 23r2 28r2 35r2

9r1 13r1 19r1 25r1 24r2 29r2 36r2

D.J. Brus et al.r Geoderma 89 (1999) 129–148

147

these were 80 or 90 Žnot shown.. A third or fourth measurement pays only for the maximum variance of 9.

4. Conclusions and concluding remarks Geographical stratification is a simple way to obtain a good coverage of randomly selected sample points within a field. The advantage over Systematic Grid Sampling is that the sample size is fixed. To enable composite sampling and unbiased estimation, the strata must be of equal size. In addition, it seems sensible to make the strata as compact as possible. This will generally result in low within-stratum variances and high sampling efficiencies. We used restricted least squares cluster analysis Ž k-means. of raster cells to divide a field into geographically compact strata of equal size. This method iteratively re-allocates the raster cells to the clusters Žstrata. , minimising the sum of the squared geographic distances of the midpoints of the raster cells to the centroides. The resulting strata are Thiessen polygons in the least squares sense. The increase in precision due to this stratification depends on the variogram. Ceteris paribus, the longer the range and the smaller the nugget to sill ratio, the stronger the stratification effect. Also, given a variogram, the stratification effect generally increases with the sampling density. For the 16 fields, the variance ratio varied from 1.2 to 2.5 if 40 points were sampled. To optimise sampling schemes we need a model to predict the precision and a cost model. The extra cost of a more precise measurement method should be weighed against the gain in precision if more sample points were sampled. To estimate the mean P concentration, replicate measurement of the composite only pays if the mean P concentration of the field to be sampled exceeds a critical level that increases with the maximum allowed variance of the total error. The sampling variance and the laboratory measurement variance are modelled as a function of the mean P concentration. Therefore, if we have no prior information on this mean, we cannot calculate the optimal scheme. However, many farmers test their soils repetitively. We suggest the use of the most recent test to obtain a priori information on the phosphate level. The quality of the optimisation result also depends on the quality of the precision model and cost model. These models still must be validated. Special attention should be paid to the subsampling error and its relation to the size of the composite sample. In this study, we considered the efficiency of one type of sampling design only, viz. Stratified Simple Random Sampling. For large fields sampling design types in which points are clustered, not by chance but purposively, may become more efficient than Stratified Simple Random Sampling if the reduction in the time for fieldwork outweighs the loss in information due to the clustering. Possible sampling design types to bring about clustering are cluster sampling and two-stage sampling. In cluster sampling, to be sure that the sample size is

148

D.J. Brus et al.r Geoderma 89 (1999) 129–148

fixed, the clusters must be of equal size. In two-stage sampling sample points are selected in two steps: firstly, primary units are selected by any sampling design, and secondly, within the selected primary units the secondary units Žpoints. are selected by any design. For instance, in the first stage we select one or more of the compact blocks that served as strata in the proposed sampling design, by Simple Random Sampling with replacement. In the second stage, we select points in the selected blocks, for instance by the Stratified Simple Random Sampling as described in Section 2.1. This is a self-weighting design, so again the composite mean is an unbiased estimator of the field-mean.

Acknowledgements The authors gratefully acknowledge the contributions of P. Barak, A.B. McBratney and A.J. Papritz to this paper.

References Black, C.A., 1992. Soil Fertility Evaluation and Control. Lewis Publishers, Boca Raton, 746 pp. Brus, D.J., De Gruijter, J.J., 1997. Random sampling or geostatistical modelling. Choosing between design-based and model-based sampling strategies for soil. Geoderma 80, 1–59, with Discussion. Burrough, P.A., 1986. Principles of Geographical Information Systems for Land Resources Assessment.Clarendon Press, Oxford, 194 pp. Cochran, W.G., 1977. Sampling Techniques, 3rd edn. Wiley, 428 pp. Cressie, N., 1985. Fitting variogram models by least squares. Mathematical Geology 17, 563–586. Domburg, P., De Gruijter, J.J., Brus, D.J., 1994. A structured approach to designing soil sampling schemes with prediction of sampling error from variograms. Geoderma 62, 151–164. Egner, ´ H., Riehm, H., Domingo, W.R., 1960. Untersuchungen uber ¨ die chemische Bodenanalyse als Grundlage fur der Boden: II. Chemische Extraction¨ de Beurteilung des Nahrstoffzustandes ¨ smethoden zur Phosphor- und Kaliumbestimmung. Kunliga Landbrukshogskolans Annaler 26, ¨ 199–215. Ferrari, Th.J., Vermeulen, F.M.B., 1955. Soil heterogeneity and soil testing. Neth. J. Agric. Sci. 3, 265–275. Hartigan, J.A., 1975. Clustering Algorithms. Wiley, New York, 351 pp. Isaaks, E.H., Srivastava, R.M., 1989. An Introduction to Applied Geostatistics. Oxford Univ. Press, New York, 561 pp. Matern, ´ B., 1986. Spatial Variation, 2nd edn. Springer-Verlag, New York, 151 pp. Sarndal, C.E., Swensson, B., Wretman, J., 1992. Model Assisted Survey Sampling. Springer¨ Verlag, New York, 694 pp. Sissingh, H.A., 1971. Analytical technique of the Pw method, used for the assessment of the phosphate status of arable soils in the Netherlands. Plant and Soil 34, 446–483. Webster, R., Oliver, M.A., 1990. Statistical Methods in Soil and Land Resource Survey. Oxford Univ. Press, 316 pp.