# Spatial prediction of topsoil salinity in the Chelif Valley, Algeria, using local ordinary kriging with local variograms versus whole-area variogram.

Introduction

The widespread and growing phenomenon of salinisation in Mediterranean areas is related largely to climatic factors and irrigation development (Isbell et al. 1983; Puigdefabregas and Mendizabal 1998). The problem is acute in Algeria even in the highest productive areas like the Chelif Valley (Boulaine 1957; Daoud 1983; Saidi et al. 1999), but no detailed spatial analysis of its extension had been undertaken.

Geostatistics has been widely applied on electrical conductivity (EC) measurements for evaluating the salinity hazard (Hajrasuliha et al. 1980; Samra et al. 1988; Oliver and Webster 1990; Hosseini et al. 1994; Sylla et al. 1995; Odeh et al. 1998; Utset et al. 1998). An overview of modern geostatistics applied to soil survey has even been exemplified on a salinity data set (Bourgault et al. 1997). Nevertheless, few studies have addressed the regional variability of soil salinity. At that areal extent, spatial trends and non-stationarity of the process may affect the accuracy of classical geostatistical approaches. In addition, variations in the spatial covariance structure should be detected as they may reveal changing factors influencing salinity spatial variability.

The present work explores the application of a particular form of kriging, namely ordinary kriging with local variogram (OKLV), to assess the topsoil EC over the Chelif Valley in Algeria. In this approach, the kriging procedure is fully dependent on the local spatial structure of the salinity process and should be able to adapt to trends in the data. The aims are therefore: (1) to map the topsoil salinity in the lower Chelif Valley, an extensively irrigated area in northern Algeria; (2) to evaluate the use of local ordinary kriging with local versus whole-area variograms, where the data set consists of intensively sampled observations.

Methods

The study area

The Chelif Valley forms a 5000 [km.sup.2] plain (0 [degrees] 60 E, 36 [degrees] 00 N) extending eastward in the northern part of Algeria along the Chief River (Fig. 1). Tertiary marl and limestone foothills bound the valley in the northern and southern parts. The study area is located in the lower part of the valley and covers 38000 ha with mainly irrigated fields used for growing cereals, fruits, and vegetables.

[Figure 1 ILLUSTRATION OMITTED]

The soils of the Chelif Valley are developed from Quaternary colluvial or alluvial material. Boulaine (1956, 1957) identified mainly soils with weak pedological organisation: the main differentiation criteria were the texture ranging from loam to light clay, the calcareous content, and the occurrence of vertic features. He was the first to notice a much larger proportion of saline soils in the lower part of the valley, later confirmed by Daoud (1983) and Saidi et al. (1999).

The lower Chelif Valley falls within the semi-arid climatic zone with an average rainfall of 370 mm in the H'madena area; this rainfall is winter dominant and occurs mainly in the form of heavy showers which may cause flooding in the lower parts of the valley. From May to September, the monthly rainfall is usually [is less than] 20 mm with average temperatures above 25 [degrees] C. During 1997, the year when sampling was conducted, the annual rainfall was 392 mm over 78 raindays, but of that 205 mm fell in only 5 days.

The sampling process

The data set consists of 5141 sampling sites, arranged on a systematic square grid with 250-m spacing between sites (Fig. 2): thirty 1: 5000 topographic maps were considered and the soil was sampled on their grid node locations where agricultural fields were present. Supplementary transects were also supplied for more detailed sampling in the western part of the area. The sampling took place during the dry season of May to September 1997.

[Figure 2 ILLUSTRATION OMITTED]

At each site the soil was sampled at the depth of 0.0-0.20 m, air-dried, and ground to pass through a 2mm sieve. A saturated paste of the sample was then prepared and the EC of the solution suction extract was measured and adjusted to 25 [degrees] C (US Salinity Laboratory Staff 1954; Rhoades 1993).

Spatial prediction: kriging with local v. whole-area variograms

Two ordinary kriging methods (Cressie 1991) were used to predict the topsoil electrical conductivity, EC*, at any unknown location, [x.sub.0], using the measured values EC([x.sub.i]) from n sampled locations [x.sub.i]. Both methods assume that:

(1) MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

These 2 methods differ in the variogram estimation and modelling procedures, i.e. the way the weights [[Lambda].sub.i] in Eqn 1 are computed.

Local ordinary kriging with a whole-area variogram (OKWV)

In this conventional approach, a unique semivariogram is inferred and modelled for the whole data set, providing a global spatial covariance structure over the whole study area. Assuming local stationarity of the mean, the kriging prediction is limited to observations in the neighbourhood of any unknown point and their weights are obtained by minimising the prediction error under a non-bias condition (Cressie 1991). The prediction variance [[Sigma].sup.2]([x.sub.0]), i.e. the minimised mean-square error, is also estimated at each unknown location through:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Where [Psi] is a Lagrange multiplier and [Gamma]() is the semivariogram model. OKWV was performed to predict EC* at points or over 100 m by 100 m blocks by fitting a double exponential function to the mean experimental variogram (Fig. 3) and adopting a neighbourhood of 100 observations around each prediction site.

[Figure 3 ILLUSTRATION OMITTED]

Ordinary kriging with a local variogram (OKLV)

Unlike the previous approach, OKLV involves the estimation and modelling of the variogram in the neighbourhood of each unknown site [x.sub.0]. A local spatial covariance model is therefore established around each point and used in Eqns 1 and 2 to estimate the property and its estimation variance.

This method was initially developed by Haas (1990a, 1990b) as a tool to map properties affected by spatial trend or changing covariance structure and applied to estimate rainfall sulfate deposition at a continental scale with a relatively sparse sample (160 monitoring sites). Further development considered the combined spatial and temporal variability of similar property (Haas 1995). On a quite finer resolution, namely an individual farm field, prediction with local variogram was also found to be well-suited to mapping crop yields using intensively measured data (with thousands of sites) obtained by continuous yield monitoring (McBratney et al. 1996; Whelan 1998). The ability of the method to adapt to situations of marked differences in local structure from the overall spatial variability was quoted as its main interest and explains its growing application in precision agriculture.

OKLV was applied to predict the EC at points or over 100 m by 100 m blocks. For each unknown site, the experimental variogram was computed considering a neighbourhood of 100 observations within a mean radius of 1650 m varying with the sampling density. This procedure was chosen to estimate the variogram with at least 100 sampling points as suggested by Webster and Oliver (1990). Haas (1995) proposed to use cross-validation to find the minimum size of the local data set, but this method appeared too computationally expensive.

An exponential model was fitted automatically using non-linear least-squares optimisation weighted by the number of pairs in the variogram. The exponential model was chosen because of its flexibility in fitting the variogram. The computation was done by using VESPER (Minasny et al. 1999) developed at The University of Sydney(1).

The validation procedures

Two validation procedures were adopted to compare the kriging techniques, testing either their prediction quality at points or their ability to evaluate the proportion of high saline soils within non-overlapping blocks.

Prediction quality at points

400 points have been selected at random and considered as a validation set. The prediction quality was determined by comparing the actual EC values with EC* kriging estimates using the remaining points. Three prediction quality criteria were considered: (i) the mean error (ME), (ii) the root mean-square error (RMSE), and (iii) the mean rank of each method (MR). If m is the number of validation points:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The reliability of the prediction variances estimated by each method was also tested on the validation set by computing the mean standardised squared error (MSSE):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

MSSE should be close to 1 if the kriging variance is an accurate estimate of the actual error variance (Wackernagel 1998).

Another way to assess the prediction variance consisted of combining the prediction estimate and prediction variance to derive a Gaussian-type confidence interval centred on the prediction estimate (Cressie 1991): a 95% confidence interval of EC* was estimated at each validation site [x.sub.1], by [EC*([x.sub.i])-2 [Sigma]([x.sub.i]); EC*([x.sub.i])+2 [Sigma]([x.sub.i])]. Computing the proportion of validation points where the 95% confidence interval actually includes the true value tested this approach. This proportion should be near 95% for an accurate modelling of uncertainty.

Prediction quality of high saline soils proportion within blocks.

We tested also the ability of the kriging methods to evaluate the proportion of high saline soils within large blocks, with the notion of, for example, selecting for remediation the most affected ones. The study area was therefore divided into 4 km by 4 km blocks (Fig. 2) and the data set was randomly equally split into prediction and validation sites. 17 blocks mostly included in the study area were considered; in each block, the actual proportion of soils with EC values exceeding a threshold of 16 dS/m was derived from 80-120 available validation sites. OKWV and OKLV were then performed at the validation sites using the prediction set, and their respective predicted proportion of values exceeding 16 dS/m was calculated. The prediction quality of this proportion was also determined as before through ME, RMSE, and MR criteria.

Results

Exploratory data analysis

The histogram of the 5141 EC data was highly skewed with values ranging from 0.2 to 58.7 dS/m and a median value of 5.4 dS/m (Fig. 4). Compared with the thresholds based on crop salt tolerance (US Salinity Laboratory Staff 1954), 36% of the samples were considered non- or slightly-saline ([is less than] 4 dS/m), 48% saline ([4-16] dS/m) and 16% highly saline ([is greater than] 16 dS/m).

[Figure 4 ILLUSTRATION OMITTED]

The location map of the data (Fig. 2) revealed sharp regional trends with varying mean EC values and spatial variability; for instance, low and almost uniform values were observed in the southern part, whereas the central part exhibits higher values and short distance variability. The semivariogram of the whole area had a large nugget effect, about 50% of the total variance (Fig. 3, circles). Two spatial variation structures could be distinguished from the inflexions of the experimental curve: a local scale (range ~ 500 m) and a regional scale (range ~ 4000 m). This variogram was best fitted up to a distance of 4000 m with a double-exponential model with a nugget value of 28.8 [(dS/m).sup.2], sill parameters of 29 and 153 [(dS/m).sup.2], and distance parameters of 150 and 17000 m:

[Gamma](N) = 28.8 + 29 [1-exp(-h/150)] +153 [1-exp (-h/17000)] if h [is less than or not equal to] 4000

Beyond 4000 m, the semi-variance appeared to be constant. When limiting the variography analysis to sub-areas, the previously observed sharp differences between sectors became apparent. The experimental variogram for the southern area indicated weak variability and lack of spatial structure (Fig. 3, squares); by contrast, the semivariogram of the central area (Fig. 3, discs) showed a huge spatial variability marked by a high nugget effect and a exponential shape increase until 1500 m. These differences between sub-areas may be related to changing topographic and lithologic conditions, but other explanation factors have also to be considered (land use, quality of irrigation water) and the spatial trends are therefore not fully understood.

Topsoil electrical conductivity prediction

The prediction maps of the electrical conductivity in the lower Chelif Valley were obtained on 37906 sites by OKWV (Fig. 5a) and OKLV (Fig. 5b) using 1-ha blocks. Both maps revealed the same trends of salinity, marked by gradual changes from low salinity areas to hot spots of very high salinity which vary from 1 ha to 8000 ha. The variance of the electrical conductivity predictions (Table 1) for both methods was much lower than the initial data, indicating a strong smoothing effect of the kriging processes.

[Figure 5 ILLUSTRATION OMITTED]

Table 1. Summary statistics of the electrical conductivity for the data set and for the kriging predictions over 1-ha blocks

The similarity of the prediction maps was confirmed by the results of the validation procedures. In the first test, i.e. the comparison to 400 validation points, the prediction errors of OKLV and OKWV (Table 2) were identical and there is no preferred method. In the second test comparing their ability to estimate the proportion of high saline soils in blocks, both methods underestimated this proportion (Table 2), but as stated by its lower mean rank, the OKLV method gave better results.

Table 2. Validation statistics of the prediction quality at points and the prediction of the proportion of highly saline soils within 16 [km.sup.2] blocks

Mean error (ME), root mean-square error (RMSE), and mean rank (MR) of the methods

Uncertainty of the predictions

By contrast with the similarity of the prediction results, OKLV and OKWV differed sharply in describing the uncertainty associated with the predictions.

For OKWV, the statistics of the prediction variance (Table 1) were characterised by a high mean value and a low variance. The corresponding map (Fig. 6a) showed slight differences in the prediction variance over the area, which were dependent on the sampling intensity. Such a pattern was expected as OKWV assumes the same covariance model over the area and only the distances between the predicted and the observed sites may affect Eqn 2.

[Figure 6 ILLUSTRATION OMITTED]

For OKLV, the prediction variance had a lower mean and higher variance (Table 1). The resulting map had a radically different pattern (Fig. 6b); large areas with very low prediction variance [([is less than] 4 dS/m).sup.2] are represented in the southern and northern parts, denoting more precise predictions of the electrical conductivity in those areas. On the other hand, areas with high prediction variance were noted both in the cases of low sampling intensity, for instance in the western part, and of high spatial variability, i.e. in central regions.

The consistency of the prediction variance for OKLV was first tested graphically by comparing the prediction errors computed on the 400 validation points against the prediction variance. Box plot statistics of the prediction errors against prediction variance classes (Fig. 7) confirmed that the prediction errors increased significantly according to the prediction variance.

[Figure 7 ILLUSTRATION OMITTED]

Validation statistics of the prediction variance gave similar results for OKLV and OKVW (Table 3); the MSSE is about 1.1 indicating a slight underestimation of the actual prediction error by both methods, which can be attributed to the presence of a few outliers with very high values in the validation set. When considering the proportion of validation points included in the 95% confidence interval, both methods gave also similar results with 92% and 93% of the points included in that interval.

Table 3. Validation statistics of the prediction variance [[Sigma].sup.2] at 400 points

Mean standardised squared error (MSSE) and proportion of true values included in the 95% confidence interval (P95)

The prediction errors on the validation set therefore appeared in agreement with the prediction variances estimated by OKLV or OKWV.

Discussion

Assessment of the topsoil salinity in the Chelif Valley

The topsoil EC in the Chelif Valley showed marked spatial structure differences between sub-areas, ranging from homogeneous low value domains to sub-areas with very high variability and superposed local and regional scales of variation. As OKLV estimates parameters of the covariance model at each prediction site, these changes can be visualised by mapping, for instance, the nugget values to assess short-range variability. Such maps may help to identify the main factors affecting salinity in each area and will be compared in an ulterior work to physiographic and irrigation practice data.

Estimation of areas affected by salinity in the Chelif Valley may be derived in first approximation from the ordinary kriging maps (Table 4); 43% of the valley would be considered as saline (8-16 dS/m) and 9% as highly saline ([is greater than] 16 dS/m). In a more rigorous approach, such an assessment must also consider the uncertainty, which may affect the salinity prediction. This may be done considering the 95% confidence interval of the electrical conductivity estimated at each unknown site [x.sub.0], by [EC*([x.sub.0])-2 [Sigma]([x.sub.0]), EC*([x.sub.0])+2 [Sigma]([x.sub.0])]. An exhaustive inventory of the areas possibly affected by salinity may therefore consider the upper limit of that confidence interval, as done in Table 4 for the OKLV estimates; this leads to a sharp increase of the areas considered as highly saline, 40% of the valley being affected by salinity predictions whose confidence interval includes values above 16 dS/m.

Table 4. Areas of topsoil EC classes according to OKLV estimates and OKLV 95% upper confidence limits

Thus, the differences between OKLV and OKWV prediction variance estimates receive here practical signification, as shown in Fig. 8; according to OKWV, 60% of the valley has a upper confidence limit above a threshold of 16 dS/m, therefore much more than the 40% estimated by OKLV. This is to be related to the high and almost uniform OKWV prediction variance, influenced only by the sampling intensity, whereas OKLV also takes into account the local variability structure.

[Figure 8 ILLUSTRATION OMITTED]

Kriging with local variograms

Until recently, the development of geostatistical approaches using local variograms, which could also embrace non-linear estimators such as multiple indicator kriging, was limited by computational time and scarcity of sampling data. The first limiting factor has been reduced by enhanced computing performance, and to a lesser extent, the development of comprehensive soil databases and increasing data acquisition potential reduced the sampling limitation.

Geostatistical approaches based on local variograms should therefore be reconsidered as they convey 2 types of improvements.

Detection and adaptation to spatial trends

A common situation in soil spatial analysis is where major scales of spatial variation are related to changes in land use, soil type, or lithology inducing differing spatial structures from one entity to another. Current approaches to deal with such situations include within-stratum kriging (Voltz and Webster 1990), kriging with local means (Wackernagel 1998), or kriging with an external drift (Bourennane et al. 1996). Where sufficient sampling intensity is available, local variograms approaches are easier to implement and less dependent on the validity of pre-existing information.

Local uncertainty assessment

Goovaerts (1997, 1999) criticised the usual confidence interval approach for modelling local uncertainty with strong emphasis given to the inability of this approach to adapt to the local variability context. For a more straightforward assessment of the local uncertainty, he promoted therefore more complex approaches such as disjunctive kriging or non-parametric approaches. Kriging with local variograms includes the local spatial structure in the uncertainty assessment and its prediction variance appeared relevant when compared with validation points. This method improves therefore the more accessible conventional approach of uncertainty assessment in case studies where broad indicators of this uncertainty are sufficient.

An alternative approach suggested by Isaaks and Srivastava (1989) and Goovaerts (1999) could also be considered, consisting of rescaling locally the sill of the global variogram to equal the variance of the data within the kriging search neighbourhood and keeping the range constant. Indeed, OKWV and OKLV yielded similar predictions and changing the local shape of variograms had therefore little impact on prediction; hence, the local rescaling of the global semivariogram could be as good as OKLV. Nevertheless, calculating the variance within the local neighbourhood is about as computationally expensive as calculating the local variogram and therefore only limited computing speed gains should be expected.

Conclusion

The OKLV method was found to be better than OKWV to describe the topsoil salinity in the Chelif Valley through the more relevant prediction variance, which was shown to be useful in delineating areas potentially affected by high salinity. Despite enhanced computational time and the need of large data sets, geostatistical approaches of soil spatial variability based on local variograms should be developed as they convey a more informative uncertainty model of the predictions.

Finally, this study confirmed the importance of the salinity hazard in the region, with 10-40% of the areas affected. The kriging maps can serve as guide to remediation and future monitoring. In particular, in selected highly saline and variable areas, more continuous and extensive approaches of soil salinity measurement and monitoring will be developed involving remote sensing and geophysical approaches.

(1) Compiled program available from Australian Centre for Precision Agriculture, University of Sydney: www.usyd.edu.au/su/agric/acpa

References

Boulaine J (1956) Carte des sols des plaines du Cheliff au 1/50.000e, feuilles 1 a 5. Inspection generale de l'Agriculture du Gouvernement General de l'Algerie.

Boulaine J (1957) Etude des sols des plaines du Chelif. These d'Etat de l'Universite al'Alger.

Bourennane H, King D, Chery P, Bruand A (1996) Improving the kriging of a soil variable using slope gradient as external drift. European Journal of Soil Science 47, 473-483.

Bourgault G, Journel AG, Rhoades LI, Corwin DL, Lesh SM (1997) Geostatistical analysis of a soil salinity data set. Advances in Agronomy 58, 241-292.

Cressie N (1991) Part 1. Geostatistical data. In `Statistics of spatial data.' pp. 29-382 (John Wiley & Sons: New York)

Daoud Y (1983) Contribution a l'etude de la dynamique des sols dans un sol irrigue du perimetre du Haut-Chelif (Algerie). PhD thesis, ENSA Rennes, France.

Goovaerts P (1997) `Geostatistics for natural resources evaluation.' (Oxford University Press: New York)

Goovaerts P (1999) Geostatistics in soil science: state-of-the-art and perspectives. Geoderma 89, 1-46.

Haas TC (1990a) Kriging and automated variogram modeling within a moving window. Atmospheric Environment 24A, 1759-1769

Haas TC (1990b) Lognormal and moving window methods of estimating acid deposition. Journal of the American Statistical Association 85, 950-963.

Haas TC (1995) Local prediction of a spatio-temporal process with an application to wet sulfate deposition. Journal of the American Statistical Association 90, 1189-1199.

Hajrasuliha S, Baniabbassi N, Metthey J, Nielsen DR (1980) Spatial variability of soil sampling for salinity studies in Southwest Iran. Irrigation Science 1, 197-208.

Hosseini E, Gallichand J, Marcotte D (1994) Theoretical and experimental performance of spatial interpolation methods for soil salinity analysis. Transactions of the American Society of Agricultural Engineers 37, 1799-1807.

Isaaks EH, Srivastava RM (1989) `An introduction to applied geostatistics.' (Oxford University Press: Oxford)

Isbell RF, Reeve R, Hutton JT (1983) Salt and sodicity. In `Soils, an Australian viewpoint', pp. 107-117. (CSIRO Publishing: Melbourne)

McBratney AB, Whelan BM, Viscarra Rossel RA (1996) Spatial prediction for precision agriculture. In `Proceedings of the 3rd International Conference on Precision Agriculture'. Minneapolis, pp. 331-342. (ASA, CSSA, SSSA: Madison, WI)

Minasny B, McBratney AB, Whelan BM (1999) VESPER version 1.0. Australian Centre for Precision Agriculture, McMillan Building A05, The University of Sydney, NSW 2006, http://www.usyd.edu.au/ su/agric/acpa.

Odeh IOA, Todd AJ, Triantafilis J, McBratney AB (1998) Status and trends of soil salinity at different scales: the case for the irrigated cotton growing region of eastern Australia. Nutrient Cycling in Agroecosystems 50, 99-107.

Oliver MA, Webster R (1990) Kriging: a method of interpolation for geographical information systems. International Journal of Geographical Systems 4, 313-332.

Puigdefabregas L, Mendizabal T (1998) Perspectives on desertification: Western Mediterranean. Journal of Arid Environments 39, 209-224.

Rhoades JD (1993) Electrical conductivity methods for measuring and mapping soil salinity. Advances in Agronomy 49, 201-251.

Saidi D, Douaoui A, Le Bissonnais Y, Walter C (1999) Sensitivity of soil surface layers from Chelif plain (Algeria) to structural degradation. Etude et Gestion des Sols 6, 15-25

Samra JS, Singh VP, Sharma KMS (1988) Analysis of spatial variability in sodic soils. 2: Point and block-kriging. Soil Science 145, 250-256.

Sylla M, Stein A, Van Breemen N, Fresco LO (1995) Spatial variability of soil salinity at different scales in the mangrove rice agro-ecosystem. West Africa Agriculture, Ecosystems and Environment 54, 1-15.

US Salinity Laboratory Staff (1954) `Diagnosis and improvement of saline and alkali soils'. (USDA Handbook No. 60)

Utset A, Ruiz ME, Herrera J, De Leon DP (1998) A geostatistical method for soil salinity sample site spacing. Geoderma 86, 143-151.

Voltz M, Webster R (1990) A comparison of kriging, cubic splines and classification for predicting soil properties from sample information. Journal of Soil Science 41,473-490.

Wackemagel H (1998) `Multivariate geostatistics: An introduction with applications.' (Springer Verlag: Berlin)

Webster R, Oliver MA (1990) `Statistical methods in soil and land resource survey.' (Oxford University Press: Oxford)

Whelan BM (1998) Reconciling continuous soil variation and crop yield. PhD thesis, The University of Sydney, Australia.

Manuscript received 15 September 1999, accepted 29 June 2000

Christian Walter(ABD), Alex B. McBratney(B), Abdelkader Douaoui(C), and Budiman Minasny(B)

(A) Soil Science Laboratory, ENSA-INRA Rennes, 35042 Rennes, France.

(B) Department of Agricultural Chemistry and Soil Science, The University of Sydney, NSW 2006, Australia.

(C) Institute of Agronomy, The University of Chlef, BP 151, Chlef, Algeria.

(D) Corresponding author; email: cwalter@roazhon.inra.fr

The widespread and growing phenomenon of salinisation in Mediterranean areas is related largely to climatic factors and irrigation development (Isbell et al. 1983; Puigdefabregas and Mendizabal 1998). The problem is acute in Algeria even in the highest productive areas like the Chelif Valley (Boulaine 1957; Daoud 1983; Saidi et al. 1999), but no detailed spatial analysis of its extension had been undertaken.

Geostatistics has been widely applied on electrical conductivity (EC) measurements for evaluating the salinity hazard (Hajrasuliha et al. 1980; Samra et al. 1988; Oliver and Webster 1990; Hosseini et al. 1994; Sylla et al. 1995; Odeh et al. 1998; Utset et al. 1998). An overview of modern geostatistics applied to soil survey has even been exemplified on a salinity data set (Bourgault et al. 1997). Nevertheless, few studies have addressed the regional variability of soil salinity. At that areal extent, spatial trends and non-stationarity of the process may affect the accuracy of classical geostatistical approaches. In addition, variations in the spatial covariance structure should be detected as they may reveal changing factors influencing salinity spatial variability.

The present work explores the application of a particular form of kriging, namely ordinary kriging with local variogram (OKLV), to assess the topsoil EC over the Chelif Valley in Algeria. In this approach, the kriging procedure is fully dependent on the local spatial structure of the salinity process and should be able to adapt to trends in the data. The aims are therefore: (1) to map the topsoil salinity in the lower Chelif Valley, an extensively irrigated area in northern Algeria; (2) to evaluate the use of local ordinary kriging with local versus whole-area variograms, where the data set consists of intensively sampled observations.

Methods

The study area

The Chelif Valley forms a 5000 [km.sup.2] plain (0 [degrees] 60 E, 36 [degrees] 00 N) extending eastward in the northern part of Algeria along the Chief River (Fig. 1). Tertiary marl and limestone foothills bound the valley in the northern and southern parts. The study area is located in the lower part of the valley and covers 38000 ha with mainly irrigated fields used for growing cereals, fruits, and vegetables.

[Figure 1 ILLUSTRATION OMITTED]

The soils of the Chelif Valley are developed from Quaternary colluvial or alluvial material. Boulaine (1956, 1957) identified mainly soils with weak pedological organisation: the main differentiation criteria were the texture ranging from loam to light clay, the calcareous content, and the occurrence of vertic features. He was the first to notice a much larger proportion of saline soils in the lower part of the valley, later confirmed by Daoud (1983) and Saidi et al. (1999).

The lower Chelif Valley falls within the semi-arid climatic zone with an average rainfall of 370 mm in the H'madena area; this rainfall is winter dominant and occurs mainly in the form of heavy showers which may cause flooding in the lower parts of the valley. From May to September, the monthly rainfall is usually [is less than] 20 mm with average temperatures above 25 [degrees] C. During 1997, the year when sampling was conducted, the annual rainfall was 392 mm over 78 raindays, but of that 205 mm fell in only 5 days.

The sampling process

The data set consists of 5141 sampling sites, arranged on a systematic square grid with 250-m spacing between sites (Fig. 2): thirty 1: 5000 topographic maps were considered and the soil was sampled on their grid node locations where agricultural fields were present. Supplementary transects were also supplied for more detailed sampling in the western part of the area. The sampling took place during the dry season of May to September 1997.

[Figure 2 ILLUSTRATION OMITTED]

At each site the soil was sampled at the depth of 0.0-0.20 m, air-dried, and ground to pass through a 2mm sieve. A saturated paste of the sample was then prepared and the EC of the solution suction extract was measured and adjusted to 25 [degrees] C (US Salinity Laboratory Staff 1954; Rhoades 1993).

Spatial prediction: kriging with local v. whole-area variograms

Two ordinary kriging methods (Cressie 1991) were used to predict the topsoil electrical conductivity, EC*, at any unknown location, [x.sub.0], using the measured values EC([x.sub.i]) from n sampled locations [x.sub.i]. Both methods assume that:

(1) MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

These 2 methods differ in the variogram estimation and modelling procedures, i.e. the way the weights [[Lambda].sub.i] in Eqn 1 are computed.

Local ordinary kriging with a whole-area variogram (OKWV)

In this conventional approach, a unique semivariogram is inferred and modelled for the whole data set, providing a global spatial covariance structure over the whole study area. Assuming local stationarity of the mean, the kriging prediction is limited to observations in the neighbourhood of any unknown point and their weights are obtained by minimising the prediction error under a non-bias condition (Cressie 1991). The prediction variance [[Sigma].sup.2]([x.sub.0]), i.e. the minimised mean-square error, is also estimated at each unknown location through:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Where [Psi] is a Lagrange multiplier and [Gamma]() is the semivariogram model. OKWV was performed to predict EC* at points or over 100 m by 100 m blocks by fitting a double exponential function to the mean experimental variogram (Fig. 3) and adopting a neighbourhood of 100 observations around each prediction site.

[Figure 3 ILLUSTRATION OMITTED]

Ordinary kriging with a local variogram (OKLV)

Unlike the previous approach, OKLV involves the estimation and modelling of the variogram in the neighbourhood of each unknown site [x.sub.0]. A local spatial covariance model is therefore established around each point and used in Eqns 1 and 2 to estimate the property and its estimation variance.

This method was initially developed by Haas (1990a, 1990b) as a tool to map properties affected by spatial trend or changing covariance structure and applied to estimate rainfall sulfate deposition at a continental scale with a relatively sparse sample (160 monitoring sites). Further development considered the combined spatial and temporal variability of similar property (Haas 1995). On a quite finer resolution, namely an individual farm field, prediction with local variogram was also found to be well-suited to mapping crop yields using intensively measured data (with thousands of sites) obtained by continuous yield monitoring (McBratney et al. 1996; Whelan 1998). The ability of the method to adapt to situations of marked differences in local structure from the overall spatial variability was quoted as its main interest and explains its growing application in precision agriculture.

OKLV was applied to predict the EC at points or over 100 m by 100 m blocks. For each unknown site, the experimental variogram was computed considering a neighbourhood of 100 observations within a mean radius of 1650 m varying with the sampling density. This procedure was chosen to estimate the variogram with at least 100 sampling points as suggested by Webster and Oliver (1990). Haas (1995) proposed to use cross-validation to find the minimum size of the local data set, but this method appeared too computationally expensive.

An exponential model was fitted automatically using non-linear least-squares optimisation weighted by the number of pairs in the variogram. The exponential model was chosen because of its flexibility in fitting the variogram. The computation was done by using VESPER (Minasny et al. 1999) developed at The University of Sydney(1).

The validation procedures

Two validation procedures were adopted to compare the kriging techniques, testing either their prediction quality at points or their ability to evaluate the proportion of high saline soils within non-overlapping blocks.

Prediction quality at points

400 points have been selected at random and considered as a validation set. The prediction quality was determined by comparing the actual EC values with EC* kriging estimates using the remaining points. Three prediction quality criteria were considered: (i) the mean error (ME), (ii) the root mean-square error (RMSE), and (iii) the mean rank of each method (MR). If m is the number of validation points:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The reliability of the prediction variances estimated by each method was also tested on the validation set by computing the mean standardised squared error (MSSE):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

MSSE should be close to 1 if the kriging variance is an accurate estimate of the actual error variance (Wackernagel 1998).

Another way to assess the prediction variance consisted of combining the prediction estimate and prediction variance to derive a Gaussian-type confidence interval centred on the prediction estimate (Cressie 1991): a 95% confidence interval of EC* was estimated at each validation site [x.sub.1], by [EC*([x.sub.i])-2 [Sigma]([x.sub.i]); EC*([x.sub.i])+2 [Sigma]([x.sub.i])]. Computing the proportion of validation points where the 95% confidence interval actually includes the true value tested this approach. This proportion should be near 95% for an accurate modelling of uncertainty.

Prediction quality of high saline soils proportion within blocks.

We tested also the ability of the kriging methods to evaluate the proportion of high saline soils within large blocks, with the notion of, for example, selecting for remediation the most affected ones. The study area was therefore divided into 4 km by 4 km blocks (Fig. 2) and the data set was randomly equally split into prediction and validation sites. 17 blocks mostly included in the study area were considered; in each block, the actual proportion of soils with EC values exceeding a threshold of 16 dS/m was derived from 80-120 available validation sites. OKWV and OKLV were then performed at the validation sites using the prediction set, and their respective predicted proportion of values exceeding 16 dS/m was calculated. The prediction quality of this proportion was also determined as before through ME, RMSE, and MR criteria.

Results

Exploratory data analysis

The histogram of the 5141 EC data was highly skewed with values ranging from 0.2 to 58.7 dS/m and a median value of 5.4 dS/m (Fig. 4). Compared with the thresholds based on crop salt tolerance (US Salinity Laboratory Staff 1954), 36% of the samples were considered non- or slightly-saline ([is less than] 4 dS/m), 48% saline ([4-16] dS/m) and 16% highly saline ([is greater than] 16 dS/m).

[Figure 4 ILLUSTRATION OMITTED]

The location map of the data (Fig. 2) revealed sharp regional trends with varying mean EC values and spatial variability; for instance, low and almost uniform values were observed in the southern part, whereas the central part exhibits higher values and short distance variability. The semivariogram of the whole area had a large nugget effect, about 50% of the total variance (Fig. 3, circles). Two spatial variation structures could be distinguished from the inflexions of the experimental curve: a local scale (range ~ 500 m) and a regional scale (range ~ 4000 m). This variogram was best fitted up to a distance of 4000 m with a double-exponential model with a nugget value of 28.8 [(dS/m).sup.2], sill parameters of 29 and 153 [(dS/m).sup.2], and distance parameters of 150 and 17000 m:

[Gamma](N) = 28.8 + 29 [1-exp(-h/150)] +153 [1-exp (-h/17000)] if h [is less than or not equal to] 4000

Beyond 4000 m, the semi-variance appeared to be constant. When limiting the variography analysis to sub-areas, the previously observed sharp differences between sectors became apparent. The experimental variogram for the southern area indicated weak variability and lack of spatial structure (Fig. 3, squares); by contrast, the semivariogram of the central area (Fig. 3, discs) showed a huge spatial variability marked by a high nugget effect and a exponential shape increase until 1500 m. These differences between sub-areas may be related to changing topographic and lithologic conditions, but other explanation factors have also to be considered (land use, quality of irrigation water) and the spatial trends are therefore not fully understood.

Topsoil electrical conductivity prediction

The prediction maps of the electrical conductivity in the lower Chelif Valley were obtained on 37906 sites by OKWV (Fig. 5a) and OKLV (Fig. 5b) using 1-ha blocks. Both maps revealed the same trends of salinity, marked by gradual changes from low salinity areas to hot spots of very high salinity which vary from 1 ha to 8000 ha. The variance of the electrical conductivity predictions (Table 1) for both methods was much lower than the initial data, indicating a strong smoothing effect of the kriging processes.

[Figure 5 ILLUSTRATION OMITTED]

Table 1. Summary statistics of the electrical conductivity for the data set and for the kriging predictions over 1-ha blocks

Counts EC data or predictions Prediction variance Mean Variance Mean Variance (dS/m) [(dS/m). [(dS/m). [(dS/m). sup.2] sup.2] sup.4] Data 5141 9.0 84.5 -- -- OKWV 37906 9.2 26.6 20.7 32.3 OKLV 37906 9.2 27.1 13.8 269.4

The similarity of the prediction maps was confirmed by the results of the validation procedures. In the first test, i.e. the comparison to 400 validation points, the prediction errors of OKLV and OKWV (Table 2) were identical and there is no preferred method. In the second test comparing their ability to estimate the proportion of high saline soils in blocks, both methods underestimated this proportion (Table 2), but as stated by its lower mean rank, the OKLV method gave better results.

Table 2. Validation statistics of the prediction quality at points and the prediction of the proportion of highly saline soils within 16 [km.sup.2] blocks

Mean error (ME), root mean-square error (RMSE), and mean rank (MR) of the methods

EC prediction at 400 % (EC > 16 dS/m) over blocks validation points ME RMSE MR ME RMSE MR (dS/m) (dS/m) (%) (%) OKWV 0.3 8.3 1.51 -8 13.5 1.88 OKLV 0.3 8.3 1.49 -6 11.4 1.12

Uncertainty of the predictions

By contrast with the similarity of the prediction results, OKLV and OKWV differed sharply in describing the uncertainty associated with the predictions.

For OKWV, the statistics of the prediction variance (Table 1) were characterised by a high mean value and a low variance. The corresponding map (Fig. 6a) showed slight differences in the prediction variance over the area, which were dependent on the sampling intensity. Such a pattern was expected as OKWV assumes the same covariance model over the area and only the distances between the predicted and the observed sites may affect Eqn 2.

[Figure 6 ILLUSTRATION OMITTED]

For OKLV, the prediction variance had a lower mean and higher variance (Table 1). The resulting map had a radically different pattern (Fig. 6b); large areas with very low prediction variance [([is less than] 4 dS/m).sup.2] are represented in the southern and northern parts, denoting more precise predictions of the electrical conductivity in those areas. On the other hand, areas with high prediction variance were noted both in the cases of low sampling intensity, for instance in the western part, and of high spatial variability, i.e. in central regions.

The consistency of the prediction variance for OKLV was first tested graphically by comparing the prediction errors computed on the 400 validation points against the prediction variance. Box plot statistics of the prediction errors against prediction variance classes (Fig. 7) confirmed that the prediction errors increased significantly according to the prediction variance.

[Figure 7 ILLUSTRATION OMITTED]

Validation statistics of the prediction variance gave similar results for OKLV and OKVW (Table 3); the MSSE is about 1.1 indicating a slight underestimation of the actual prediction error by both methods, which can be attributed to the presence of a few outliers with very high values in the validation set. When considering the proportion of validation points included in the 95% confidence interval, both methods gave also similar results with 92% and 93% of the points included in that interval.

Table 3. Validation statistics of the prediction variance [[Sigma].sup.2] at 400 points

Mean standardised squared error (MSSE) and proportion of true values included in the 95% confidence interval (P95)

MSSE P95 (%) OKWV 1.17 93 OKLV 1.13 92

The prediction errors on the validation set therefore appeared in agreement with the prediction variances estimated by OKLV or OKWV.

Discussion

Assessment of the topsoil salinity in the Chelif Valley

The topsoil EC in the Chelif Valley showed marked spatial structure differences between sub-areas, ranging from homogeneous low value domains to sub-areas with very high variability and superposed local and regional scales of variation. As OKLV estimates parameters of the covariance model at each prediction site, these changes can be visualised by mapping, for instance, the nugget values to assess short-range variability. Such maps may help to identify the main factors affecting salinity in each area and will be compared in an ulterior work to physiographic and irrigation practice data.

Estimation of areas affected by salinity in the Chelif Valley may be derived in first approximation from the ordinary kriging maps (Table 4); 43% of the valley would be considered as saline (8-16 dS/m) and 9% as highly saline ([is greater than] 16 dS/m). In a more rigorous approach, such an assessment must also consider the uncertainty, which may affect the salinity prediction. This may be done considering the 95% confidence interval of the electrical conductivity estimated at each unknown site [x.sub.0], by [EC*([x.sub.0])-2 [Sigma]([x.sub.0]), EC*([x.sub.0])+2 [Sigma]([x.sub.0])]. An exhaustive inventory of the areas possibly affected by salinity may therefore consider the upper limit of that confidence interval, as done in Table 4 for the OKLV estimates; this leads to a sharp increase of the areas considered as highly saline, 40% of the valley being affected by salinity predictions whose confidence interval includes values above 16 dS/m.

Table 4. Areas of topsoil EC classes according to OKLV estimates and OKLV 95% upper confidence limits

EC class (dS/m) Total [0 - 4] [4 - 8] [8 - 16] >16 EC prediction [km.sup.2] 39.99 142.82 163.36 32.89 379.06 % 10 38 43 9 100 EC upper 95% confidence limit [km.sup.2] 18.14 36.33 172.89 151.70 379.06 % 5 10 45 40 100

Thus, the differences between OKLV and OKWV prediction variance estimates receive here practical signification, as shown in Fig. 8; according to OKWV, 60% of the valley has a upper confidence limit above a threshold of 16 dS/m, therefore much more than the 40% estimated by OKLV. This is to be related to the high and almost uniform OKWV prediction variance, influenced only by the sampling intensity, whereas OKLV also takes into account the local variability structure.

[Figure 8 ILLUSTRATION OMITTED]

Kriging with local variograms

Until recently, the development of geostatistical approaches using local variograms, which could also embrace non-linear estimators such as multiple indicator kriging, was limited by computational time and scarcity of sampling data. The first limiting factor has been reduced by enhanced computing performance, and to a lesser extent, the development of comprehensive soil databases and increasing data acquisition potential reduced the sampling limitation.

Geostatistical approaches based on local variograms should therefore be reconsidered as they convey 2 types of improvements.

Detection and adaptation to spatial trends

A common situation in soil spatial analysis is where major scales of spatial variation are related to changes in land use, soil type, or lithology inducing differing spatial structures from one entity to another. Current approaches to deal with such situations include within-stratum kriging (Voltz and Webster 1990), kriging with local means (Wackernagel 1998), or kriging with an external drift (Bourennane et al. 1996). Where sufficient sampling intensity is available, local variograms approaches are easier to implement and less dependent on the validity of pre-existing information.

Local uncertainty assessment

Goovaerts (1997, 1999) criticised the usual confidence interval approach for modelling local uncertainty with strong emphasis given to the inability of this approach to adapt to the local variability context. For a more straightforward assessment of the local uncertainty, he promoted therefore more complex approaches such as disjunctive kriging or non-parametric approaches. Kriging with local variograms includes the local spatial structure in the uncertainty assessment and its prediction variance appeared relevant when compared with validation points. This method improves therefore the more accessible conventional approach of uncertainty assessment in case studies where broad indicators of this uncertainty are sufficient.

An alternative approach suggested by Isaaks and Srivastava (1989) and Goovaerts (1999) could also be considered, consisting of rescaling locally the sill of the global variogram to equal the variance of the data within the kriging search neighbourhood and keeping the range constant. Indeed, OKWV and OKLV yielded similar predictions and changing the local shape of variograms had therefore little impact on prediction; hence, the local rescaling of the global semivariogram could be as good as OKLV. Nevertheless, calculating the variance within the local neighbourhood is about as computationally expensive as calculating the local variogram and therefore only limited computing speed gains should be expected.

Conclusion

The OKLV method was found to be better than OKWV to describe the topsoil salinity in the Chelif Valley through the more relevant prediction variance, which was shown to be useful in delineating areas potentially affected by high salinity. Despite enhanced computational time and the need of large data sets, geostatistical approaches of soil spatial variability based on local variograms should be developed as they convey a more informative uncertainty model of the predictions.

Finally, this study confirmed the importance of the salinity hazard in the region, with 10-40% of the areas affected. The kriging maps can serve as guide to remediation and future monitoring. In particular, in selected highly saline and variable areas, more continuous and extensive approaches of soil salinity measurement and monitoring will be developed involving remote sensing and geophysical approaches.

(1) Compiled program available from Australian Centre for Precision Agriculture, University of Sydney: www.usyd.edu.au/su/agric/acpa

References

Boulaine J (1956) Carte des sols des plaines du Cheliff au 1/50.000e, feuilles 1 a 5. Inspection generale de l'Agriculture du Gouvernement General de l'Algerie.

Boulaine J (1957) Etude des sols des plaines du Chelif. These d'Etat de l'Universite al'Alger.

Bourennane H, King D, Chery P, Bruand A (1996) Improving the kriging of a soil variable using slope gradient as external drift. European Journal of Soil Science 47, 473-483.

Bourgault G, Journel AG, Rhoades LI, Corwin DL, Lesh SM (1997) Geostatistical analysis of a soil salinity data set. Advances in Agronomy 58, 241-292.

Cressie N (1991) Part 1. Geostatistical data. In `Statistics of spatial data.' pp. 29-382 (John Wiley & Sons: New York)

Daoud Y (1983) Contribution a l'etude de la dynamique des sols dans un sol irrigue du perimetre du Haut-Chelif (Algerie). PhD thesis, ENSA Rennes, France.

Goovaerts P (1997) `Geostatistics for natural resources evaluation.' (Oxford University Press: New York)

Goovaerts P (1999) Geostatistics in soil science: state-of-the-art and perspectives. Geoderma 89, 1-46.

Haas TC (1990a) Kriging and automated variogram modeling within a moving window. Atmospheric Environment 24A, 1759-1769

Haas TC (1990b) Lognormal and moving window methods of estimating acid deposition. Journal of the American Statistical Association 85, 950-963.

Haas TC (1995) Local prediction of a spatio-temporal process with an application to wet sulfate deposition. Journal of the American Statistical Association 90, 1189-1199.

Hajrasuliha S, Baniabbassi N, Metthey J, Nielsen DR (1980) Spatial variability of soil sampling for salinity studies in Southwest Iran. Irrigation Science 1, 197-208.

Hosseini E, Gallichand J, Marcotte D (1994) Theoretical and experimental performance of spatial interpolation methods for soil salinity analysis. Transactions of the American Society of Agricultural Engineers 37, 1799-1807.

Isaaks EH, Srivastava RM (1989) `An introduction to applied geostatistics.' (Oxford University Press: Oxford)

Isbell RF, Reeve R, Hutton JT (1983) Salt and sodicity. In `Soils, an Australian viewpoint', pp. 107-117. (CSIRO Publishing: Melbourne)

McBratney AB, Whelan BM, Viscarra Rossel RA (1996) Spatial prediction for precision agriculture. In `Proceedings of the 3rd International Conference on Precision Agriculture'. Minneapolis, pp. 331-342. (ASA, CSSA, SSSA: Madison, WI)

Minasny B, McBratney AB, Whelan BM (1999) VESPER version 1.0. Australian Centre for Precision Agriculture, McMillan Building A05, The University of Sydney, NSW 2006, http://www.usyd.edu.au/ su/agric/acpa.

Odeh IOA, Todd AJ, Triantafilis J, McBratney AB (1998) Status and trends of soil salinity at different scales: the case for the irrigated cotton growing region of eastern Australia. Nutrient Cycling in Agroecosystems 50, 99-107.

Oliver MA, Webster R (1990) Kriging: a method of interpolation for geographical information systems. International Journal of Geographical Systems 4, 313-332.

Puigdefabregas L, Mendizabal T (1998) Perspectives on desertification: Western Mediterranean. Journal of Arid Environments 39, 209-224.

Rhoades JD (1993) Electrical conductivity methods for measuring and mapping soil salinity. Advances in Agronomy 49, 201-251.

Saidi D, Douaoui A, Le Bissonnais Y, Walter C (1999) Sensitivity of soil surface layers from Chelif plain (Algeria) to structural degradation. Etude et Gestion des Sols 6, 15-25

Samra JS, Singh VP, Sharma KMS (1988) Analysis of spatial variability in sodic soils. 2: Point and block-kriging. Soil Science 145, 250-256.

Sylla M, Stein A, Van Breemen N, Fresco LO (1995) Spatial variability of soil salinity at different scales in the mangrove rice agro-ecosystem. West Africa Agriculture, Ecosystems and Environment 54, 1-15.

US Salinity Laboratory Staff (1954) `Diagnosis and improvement of saline and alkali soils'. (USDA Handbook No. 60)

Utset A, Ruiz ME, Herrera J, De Leon DP (1998) A geostatistical method for soil salinity sample site spacing. Geoderma 86, 143-151.

Voltz M, Webster R (1990) A comparison of kriging, cubic splines and classification for predicting soil properties from sample information. Journal of Soil Science 41,473-490.

Wackemagel H (1998) `Multivariate geostatistics: An introduction with applications.' (Springer Verlag: Berlin)

Webster R, Oliver MA (1990) `Statistical methods in soil and land resource survey.' (Oxford University Press: Oxford)

Whelan BM (1998) Reconciling continuous soil variation and crop yield. PhD thesis, The University of Sydney, Australia.

Manuscript received 15 September 1999, accepted 29 June 2000

Christian Walter(ABD), Alex B. McBratney(B), Abdelkader Douaoui(C), and Budiman Minasny(B)

(A) Soil Science Laboratory, ENSA-INRA Rennes, 35042 Rennes, France.

(B) Department of Agricultural Chemistry and Soil Science, The University of Sydney, NSW 2006, Australia.

(C) Institute of Agronomy, The University of Chlef, BP 151, Chlef, Algeria.

(D) Corresponding author; email: cwalter@roazhon.inra.fr

Printer friendly Cite/link Email Feedback | |

Author: | Walter, Christian; McBratney, Alex B.; Douaoui, Abdelkader; Minasny, Budiman |
---|---|

Publication: | Australian Journal of Soil Research |

Geographic Code: | 6ALGE |

Date: | Mar 1, 2001 |

Words: | 4400 |

Previous Article: | Tillage and traffic effects on runoff. |

Next Article: | Land suitability assessment in the Namoi Valley of Australia, using a continuous model. |

Topics: |