Uncertainty assessment of soil organic carbon content spatial distribution using geostatistical stochastic simulation.
Abstract: Soil organic carbon (SOC) affects many processes in soil. The main objective of this study was the prediction and uncertainty assessment of the spatial patterns of SOC through stochastic simulation using 2 simulation algorithms, sequential Gaussian simulation (sGs) and sequential indicator simulation (sis). The dataset consisted of 158 point measurements of surface SOC taken from an 18-ha field in Lower Austria. Conditional stochastic simulation algorithms were used to generate 100 maps of equiprobable spatial distribution for SOC.

In general the simulated maps represented spatial distribution of SOC more realistically than the kriged map, i.e. overcoming the smoothing effect of kriging. Unlike sGs, sis was able to preserve the connectivity of extreme values in generated maps. The SOC simulated maps generated through sGs reproduced the sample statistics well. The reproduction of class-specific patterns of spatial continuity of SOC for the simulated model produced through sis was also reasonably good. The results highlight that when the class-specific patterns of spatial continuity of the attribute must be preserved, sis is preferred to sGs. For local uncertainty, standard deviations obtained using kriging varied much less across the study area than those obtained using simulations. This shows that the conditional standard deviations achieved through simulations depend on data values in addition to data configuration for greater reliability in reporting the estimation precision. Further, according to accuracy plots and goodness statistic, G, sis performs the modelling uncertainty better than sGs. The simulated models can provide useful information in risk assessment of SOC management in Lower Austria.

Additional keywords: prediction, kriging, sequential Gaussian simulation, sequential indicator simulation.
Article Type: Report
Subject: Soil management (Methods)
Stochastic analysis (Methods)
Geology (Statistical methods)
Geology (Usage)
Soils (Carbon content)
Soils (Observations)
Authors: Delbari, Masoomeh
Loiskandl, Willibald
Afrasiab, Peyman
Pub Date: 02/01/2010
Publication: Name: Australian Journal of Soil Research Publisher: CSIRO Publishing Audience: Academic Format: Magazine/Journal Subject: Agricultural industry; Earth sciences Copyright: COPYRIGHT 2010 CSIRO Publishing ISSN: 0004-9573
Issue: Date: Feb, 2010 Source Volume: 48 Source Issue: 1
Geographic: Geographic Scope: Australia Geographic Code: 8AUST Australia
Accession Number: 222313804
Full Text: Introduction

Soil organic carbon (SOC) is an important soil characteristics, affecting many processes such as soil degradation, surface crusting, runoff, and erosion in soil (Tisdall and Oades 1982; Le Bissonnais and Arrouays 1997). Soil organic matter improves soil structure and workability, enhances aeration and water penetration, increases water-holding capacity, and stores and supplies nutrients for growth of plants and soil microorganisms. SOC increases the cohesion of aggregates, stabilises, holds soil particles together against disruption upon rainfall, and reduces erosion hazard (Le Bissonnais and Arrouays 1997; Chenu et al. 2000). Accurate mapping of SOC content can help us to identify the relative importance of the above-mentioned processes. Accurate and detailed mapping of SOC can also be used in precision farming to optimise the crop production inputs such as fertilizers and chemicals (Chen et al. 2005) or in environmental research related to terrestrial sequestration of atmospheric carbon (Simbahan et al. 2006).

Geostatistics has been used for spatial characterisation and prediction of SOC content (Van Meirvenne et al. 1996; Chen et al. 2000; Mueller and Pierce 2003; Simbahan et al. 2006). There is an increasing awareness of assessing the uncertainty attached to the estimates of SOC at unsampled locations, because whatever prediction method is used, some degree of uncertainty in the estimated values will be introduced. This uncertainty affects subsequent decision-making or will propagate to the response uncertainty through a modelling process (Mowrer 2000). Ordinary kriging is the most popular geostatistical interpolation method. The estimated values at unsampled locations are obtained by linear weighted moving averaging of the values at known locations. A kriged estimate can be combined with the byproduct estimation variance to build a Gaussian-type confidence interval at any unsampled location, assuming a normal distribution for the error (Goovaerts 1997). Ordinary kriging (OK) aims, however, to minimise local error variance (i.e. looking for local accuracy) without regard to the global statistics of the estimates. Moreover, the estimation variance is a function of data arrangement and semivariogram model characteristics and is not dependent on data values. Therefore, in many cases, OK estimation variance cannot be used as a reliable measure of estimation precision (Goovaerts 1997).

An alternative approach is geostatistical stochastic simulation (Goovaerts 1997; Deutsch and Journel 1998). Stochastic simulation generates several alternative, equally probable maps as realisations of the spatial distribution of attributes, in contrast to a map of local best estimates produced by kriging. Conditional stochastic simulation, unlike kriging, reproduces the sample statistics (histogram and model of spatial continuity, e.g. semivariogram model) and honours sample data at their original locations. The simulation is based on both spatial variation of true values and variation of estimated values, whereas kriging takes into account only the spatial variation of observed data at sampled locations, ignoring variation of the estimated value. Therefore, a stochastic simulation map represents the spatial distribution of the attribute more realistically than a smoothed kriged map. Furthermore, a set of alternative realisations can be used to provide a visual and quantitative measure of space of uncertainty, which is then incorporated in related decision-making and risk analysis (Goovaerts 1997).

Simulation algorithms differ in the underlying random function model (multiGaussian or nonparametric), the amount and type of information that can be accounted for, and the computer requirements (Goovaerts 1997). Sequential Gaussian simulation (sGs) and sequential indicator simulation (sis) are the 2 most commonly used algorithms. The sGs is very fast and mathematically simple, but its multi-Gaussianarity assumption of random function model imposes maximum entropy, i.e. maximum randomness beyond the input variogram and data histogram (Deutsch and Journel 1998). That is, generated realisations through sGs show no significant correlation of extreme values. The sGs algorithm is thus not appropriate in applications where it is critical to preserve the connectivity of low or high values. In contrast to Gaussian-related models, sis has the flexibility to reproduce patterns of spatial connectivity specific to extreme values through different indicator semivariogram models (Goovaerts 1997).

Geostatistical stochastic simulation has been applied in many scientific disciplines, such as environmental studies (Englund 1993; Rossi et al. 1993; Juang et al. 2004; Schnabel et al. 2004) and recently in soil science (Pachepsky and Acock 1998; Goovaerts 2001; Bourennane et al. 2007). Stochastic simulation has been also used for mapping and uncertainty assessment of spatial pattern of SOC density (Zhao et al. 2005).


The main objective of this study is the prediction and uncertainty assessment of the spatial patterns of SOC through stochastic simulation. For this purpose, the feasibility of using 2 simulation algorithms, sGs and sis, is examined. The results are compared with those obtained by the more classical approach of ordinary kriging (OK).

Materials and methods

Study area and dataset

The study was conducted on an 18-ha sloped field (16[degrees]34'E, 48[degrees]34'N) in a 40-ha hilly agricultural watershed. The field is in Mistelbach, a town in Lower Austria, one of Austria's 9 Federal provinces. It is ~40km north-east of Vienna and ~25-30km from the Czech and the Slovak boarders (Fig. 1a). Geologically the area is in the northern part of so-called Viennese basin formed by alluvial deposits (Molasse zone; Zotl 1997). Soil texture ranges from silt loam to loam. Mean annual temperature in Mistelbach for 1994-2001 was ~9.6[degrees]C with minimum and maximum mean monthly temperature of -1.1 and 20[degrees]C, respectively. During this time the average annual precipitation was ~665 mm. Sampling was performed on an almost regular grid at 158 locations, ~30m apart (Fig. 1b). During 2001-02, soil samples were taken from different depth intervals, and the SOC content was determined for each sample. For this study, SOC content data of depth 0-0.20 m were considered.


Ordinary kriging

For kriging estimation, the spatial dependence of the SOC content was first quantified using semivariogram (Isaaks and Srivastava 1989). The SOC content at an unsampled location x0 was then estimated using kriging estimator [z.sup.*]([x.sub.0]), defined as:

[z.sup.*]([x.sub.0]) = [n.summation over (i=1)] [[lambda].sub.i] x z([x.sub.i]) (1)

where [z.sup.*]([x.sub.0]) is the estimated SOC content at unsampled location [x.sub.0], [[lambda].sub.i] is the weight assigned to the known SOC content at location [x.sub.i] determined based on a semivariogram model, and n is the number of observations to be used. The weights are allocated to the unsampled locations in such a way that they sum to unity and minimise the kriging estimation variance (Isaaks and Srivastava 1989).

The OK estimation variance at location [x.sub.0], [[sigma].sup.2.sub.OK] ([x.sub.0]), is given as:

[[sigma].sup.2.sub.OK]([x.sub.0]) = [n.summation over (i=1)] [[lambda].sub.i][gamma]([x.sub.0], [x.sub.i]) + [mu] (2)

where [gamma]([X.sub.0], [x.sub.i]) is the average semivariance between the location to be estimated ([x.sub.0]) and the ith sample point and [mu] is the Lagrange parameter for minimisation of the kriging variance. The error variance [[sigma].sup.2.sub.OK]([x.sub.0]) or error standard deviation [[sigma].sub.OK] ([x.sub.0]) can then be combined with the estimate to derive a Gaussian-type confidence interval around the estimated value at [x.sub.0] (Goovaerts 1997).

Stochastic simulation

The main features of the 2 simulation algorithms used to generate simulated models of SOC spatial distribution are reviewed briefly in this section. A complete presentation of stochastic simulation is available in Goovaerts (1997) and Deutsch and Journel (1998).

Consider the simulation of the continuous variable of SOC content, z, at N grid nodes [X.sub.0j], j = 1, ..., N, conditional to n sample data available {z([x.sub.[alpha]]), [alpha] = 1, ..., n}. Sequential simulation (Goovaerts 1997) involves modelling the ccdf F([x.sub.0j]; z[absolute value of (n)) = Prob{Z([x.sub.0j]) [less than or equal to] z](n)} and then sampling it at each node of the grid visited along a random sequence. Each ccdf is estimated conditional to both the original n data and all previously simulated values. The sGs and sis are 2 major classes of the sequential simulation algorithms in which the series of ccdfs are determined using multiGaussian or the indicator formalism, respectively.

Sequential Gaussian simulation

The sGs procedure after defining a regularly spaced grid covering the region of interest involves the following steps:

(1) Transform the sample data values of SOC into standard normal data using a normal score transformation.

(2) Compute and model the semivariogram of the normal score transformed data.

(3) Establish a random path through all of the grid nodes, in a way that each node is visited only once in each sequence.

(4) At each node [x.sub.0]:

(a) Estimate the parameters (mean and variance) of the Gaussian ccdf of SOC by simple kriging estimator with the normal score semivariogram model. The conditional dataset includes a specified number of both the original data and already simulated data within a neighbourhood of the location to be simulated.

(b) Draw a simulated normal score value from the estimated ccdf and then add to the conditional dataset to be used for simulating other nodes.

(c) Proceed to the next grid node along the random path and repeat steps (a) and (b) until the entire grid nodes are simulated.

(5) Back-transform the simulated normal scores into the simulated values of SOC in the original space.

These sequential steps build up only the first realisation, {[z.sup.(1)])([x.sub.0j]), j = 1, ..., N}, which is only one equally probable model of SOC spatial distribution. To generate multiple, say L, realisations, {[z.sup.(1)]([x.sub.0j]), j = 1,..., N, l = 1, ..., L}, steps 3-5 should be repeated with different random paths passing through all nodes. In fact the only difference among several individual simulations is the random-number seed, which initiates the simulation process and affects the random drawing of a value from the ccdf used to generate each simulated value. This is why, in different realisations, the same location will be most likely assigned different values, producing maps that are variable in their specific appearance.

Sequential indicator simulation

The sis, relying on indicator kriging to infer the ccdfs, proceeds as follows:

(1) Discretise the range of variation of SOC data values into (K + 1) classes using K threshold values [z.sub.k] and then transform each SOC observed value z ([x.sub.[alpha]]) into a vector of indicator values defined as:


(2) Compute and model the K indicator semivariograms.

(3) Establish a random path passing through all nodes to be simulated, visiting each point only once.

(4) At each node Xo along the path:

(a) Estimate the indicator random variable I ([x.sub.0]; [z.sub.k]) for each of the K thresholds using indicator kriging (Goovaerts 1997). This provides the probability of SOC value z ([x.sub.0]) not exceeding a given threshold value:

F[[x.sub.0]; [z.sub.k][absolute value of [(n)] = Prob[Z([x.sub.0]) [less than or equal to] [z.sub.k]](n)] (4)

where (n) is the conditioning data including indicator transforms of neighbouring original data and previously simulated values. Given K threshold values, K probability values are estimated at unsampled location [x.sub.0].

(b) Correct for any order relation deviation, and construct a continuous ccdf using interpolation/extrapolation algorithms (Deutsch and Journel 1998).

(c) Draw randomly a simulated value of SOC content from the prior ccdf for that location and add its indicator code to the conditioning indicator variables of Eqn 4 for modelling the prior ccdf at the next location.

(d) Proceed to the next location along the random path and repeat steps (a)-(c).

To obtain several realisations, the sequential processes are repeated with different random paths. One hundred equally probable realisations of SOC spatial distribution are generated using both algorithms. The simulated realisations numerically approximate the ccdf for every single location, and after being post-processed, different measures of optimal estimate and uncertainty are obtained. For example, pixel-by-pixel averaging several generated models can provide a map of the expected value (E-type estimate; Journel 1983) at any location. A corresponding conditional variance map, which represents the spread of estimated ccdfs, can also be produced. The estimation and simulation are performed on a grid of 5150 nodes (5 m by 5m) covering the study area. Kriging and simulations are performed using S-GeMS (Remy 2004).

Validation of uncertainty model

To evaluate the simulation performance, the accuracy and goodness of the reproduction of l-point conditional probability distribution by the set of realisations is examined. Indeed, at any location x, knowledge of the ccdf F (x, z |(n)) allows the computation of a series of symmetric p-probability intervals (PI) bounded by the (1 -p)/2 and (1 + p)/2 quantiles of that ccdf. According to Goovaerts (2001), a probability distribution is considered accurate if the fraction of values of a validation dataset falling in the symmetric p-PI exceeds p for all p [member of] [0, 1]. Considering a validation set, z ([x.sub.j]), j = 1, ..., [N.sub.v] = 51, and the corresponding independently predicted ccdfs F ([x.sub.j], z |(n)), the fraction of true values falling into a given symmetric p-PI interval can be computed as:


The scattergram of the calculated fractions v. the set of probabilities p is called 'accuracy plot' (Goovaerts 2001). The accuracy of a probabilistic model can be visually assessed through an accuracy plot; for an accurate case, most of the points must fall above the 45[degrees] line, i.e. [bar.[xi]](p) > p for most p. Deutsch (1997) also proposed to assess the closeness of the estimated and theoretical fractions using the so called 'goodness' statistic G [member of] [0, 1] as:


where the indicator function a(p) is defined as:


As seen in Eqn 6, twice as much importance is given to deviations when [bar.[xi]](p) < p (inaccurate case) because the weight [absolute value of 3a(p) 2] is equal to 2 instead of 1 for the accurate case. Another characteristic of the probabilistic models is the spread of the local distributions of L (=100) simulated values at each validation site [x.sub.j], a measure of local uncertainty. The average variance of the [N.sub.v](=51) distributions of simulated values, 1 per each validation site, is defined by the statistic U as (Goovaerts 1999):


where [z.sup.*][sub.E]([x.sub.j]) is the E-type estimate, i.e. the arithmetic average of the 100 simulated values at [x.sub.j], that is the mean of the local distribution. Statistic U should be as small as possible while preserving accuracy and precision; the smaller U, the less is local uncertainty. Validation of uncertainty models is performed using GSLIB (Deutsch and Journel 1998).

Results and discussion

Descriptive statistics and spatial continuity

Basic statistics and histogram of SOC raw data are presented in Fig. 2a. The high coefficients of skewness and kurtosis show a positively skewed distribution of surface SOC data. The spatial pattern of SOC measurements in Fig. 1 b reveals that the spatial distribution of SOC content is not homogeneous. There are a few locations in the east of the study area where SOC content is high, because they belong to a small forest area. The experimental semivariogram of SOC is calculated for 4 directions with 45[degrees] angular increments and [+ or -] 22.5[degrees] angular tolerance. They show no signs of anisotropy, and the omni-directional semivariogram of SOC depicted in Fig. 3a is therefore considered for further analysis. Topsoil organic carbon shows moderate correlation in space with an exponential model of spatial structure having a nugget effect of 0.026 ([%.sup.2]). The semivariogram increases until the effective range of ~120 m and flattens out afterwards. The sill of semivariogram is 0.09 ([%.sup.2]).

For the sGs implementation, raw SOC data were first transformed to a set of normal scores. The normal score transformed data are normally distributed around a mean of 0 with a variance of ~1 (Fig. 2b). The experimental semivariogram of SOC normal scores (Fig. 3b) follows an exponential model, which reaches the sill of unity at an effective range of 95 m. For the sis implementation, the experimental semivariograms of the indicators are computed for 7 indicator cutoffs [0.95, 1.02, 1.15, 1.29, 1.41, 1.56, and 1.67(%)] corresponding, respectively, to 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, and 0.95 quantiles of SOC data distribution. The best model in the least-square sense was fitted to the indicator semivariograms and the results are presented in Table 1. Indicator semivariograms have spherical structure of spatial continuity except for 0.1 quantile. Regarding the ratio of nugget effect to the sill, the semivariograms of indicators show a moderate within-class spatial correlation for SOC content.

Mapping of SOC

Kriging yields a single smoothing prediction map of SOC content (Fig. 4). Unlike kriging, simulation algorithms produce several possible models of SOC spatial distribution (Fig. 5); however, there are some differences between the models obtained by the 2 simulation algorithms. Generated realisations through sGs show no significant correlation of extreme values. This is because of the multi-Gaussianarity assumption of random function model in sGs, which imposes maximum entropy, i.e. maximum randomness beyond the input semivariogram and data histogram (Deutsch and Journel 1998). Unlike sGs, sis is able to preserve the connectivity of extreme values in generated maps. The sGs algorithm is thus not appropriate in applications where it is critical to preserve the connectivity of low or high values. The E-type estimate maps (Fig. 6) are, however, similar to the kriged map. This is expected because it is well known that the mean of simulated maps converges to a kriged map when the number of generated realisations is very large (Chiles and Delfiner 1999). The validity of the simulation results was examined by testing the reproduction of the sample histogram and (indicator) semivariogram models by the realisations (Leuangthong et al. 2004). As seen in Figs 7 and 8, the simulated maps reproduce the sample histogram much better than the kriged map. For instance, the distribution of estimated values has a much smaller variance than the sample variance (0.035 against 0.085). The reproduction of the sample semivariogram model by the maps obtained through kriging and sGs are presented in Figs 7 and 9, respectively. The sample semivariogram model reproduction by kriging is not satisfactory; it underestimates the short-range variability represented in a semivariogram with a smaller relative nugget effect. The simulated maps reproduce the sample semivariogram model better than kriging. However, sGs does not guarantee the reproduction of the semivariogram model if the sample histogram is too skewed (Caers 2000). As sis, by theory, reproduces only the indicator semivariogram models (Goovaerts 1997), the first realisation was examined for the reproduction of the 3 indicator semivariogram models corresponding to 0.25, 0.5, and 0.75 quantiles of the cumulative distribution function of SOC observations (see Fig. 9). For comparison, the first realisation produced though sGs was also examined for the reproduction of class-specific patterns of spatial continuity (Fig. 9). The results indicate that sis reproduces better the spatial continuity of the small, median, and large values. Reproduction of indicator semivariograms is important when the connectivity of extreme values existing in data needs to be reproduced in the simulated realisations (Caers 2000). Therefore sis, unlike sGs, allows the patterns of spatial continuity specific to different classes of values to be taken into account, so the simulated models produced look more realistic.



Uncertainty assessment of SOC

Accuracy plots related to sGs and sis algorithms were computed using Eqn 5 and depicted in Fig. 10. For the sGs algorithm, almost all the points fall above the 45[degrees] line. This reflects the accuracy of sGs in modelling estimation uncertainty. For the sis algorithm, some points are below the bisector line. This indicates that the probabilistic model provided by sis is inaccurate in modelling estimation uncertainty for large p-values (P > 0.7). However, departures from the bisector line are smaller for the non-parametric algorithm, resulting in a smaller goodness statistic G (Fig. 10). The statistic U, the average variance of the 51 distributions of simulated values at the validation sites, was computed using Eqn 7 and represented in Fig. 10. Accordingly, sis algorithm achieves a smaller magnitude of U, indicating smaller local uncertainty.








The performance of estimation and simulations is also examined by comparing the spread of the standard deviation maps (Fig. 11) generated through each algorithm. Standard deviation maps provided by both simulation algorithms show greater certainty than that generated by the OK algorithm, mainly in the upper part of the study area. The local standard deviation associated with kriged estimates varies much less across the study area (Fig. 11a). However, the smallest standard deviation reflecting smallest uncertainty is seen at sampling locations and for nearby locations, whereas largest uncertainty belongs to areas of more distant or no samples. This drawback of kriging standard deviation of independence from the actual sample value is overcome using simulation algorithms. Figure 11b and c shows that higher uncertainty (wider ccdfs) is seen mostly in high-valued areas of SOC content, e.g. south-east of the study area; compare with the map of E-type estimates in Fig. 6. These locations are usually considered for the additional sampling. This relation is confirmed by plotting the scattergram of E-type estimate v. local standard deviation (Fig. 12). Indeed, there is a good correlation between the local standard deviations and local mean values obtained through both simulation algorithms (Fig. 12b, c). For example, small to intermediate local mean values of SOC exhibit great certainty across the study area. Unlike simulation, OK results in no correlation between SOC estimates and their associated standard deviations (Fig. 12a). This confirms once again that for estimation algorithm, standard deviation of the errors is independent of the actual data values, and depends mainly on data configurations. Finally, according to the standard deviation maps and the scattergrams of local standard deviation v. local mean values, sis preserving the class-specific pattern of spatial continuity achieves smaller local uncertainty than other approaches mainly for parts of the study area where small to intermediate SOC values were measured. The findings obtained here highlight the need for precise assessment of estimation uncertainty, especially when this information should be incorporated in subsequent decision-making such as identification of low- and high-valued areas of SOC for better management of SOC in Lower Austria.




An empirical comparison of estimation- and simulation-based algorithms for spatial prediction and modelling uncertainty of SOC content is provided. Conditional stochastic simulation algorithms including sequential Gaussian and indicator simulation are used to generate 100 realisations of equally probable spatial distribution of SOC content for an experiment field in Lower Austria. In general, the simulated maps obtained through both simulation algorithms represent spatial distribution of SOC more realistically than a single local best estimation map, i.e. it overcomes the smoothing effect of kriging. Unlike sGs, sis has the ability to preserve the connectivity of extreme values in generated maps. The SOC simulated models generated through sGs reproduce, on average, the sample statistics well. The reproduction of class-specific patterns of spatial continuity of SOC by the simulated map produced through sis is reasonably good.

With regard to estimation uncertainty, the map of standard deviations associated with the kriged estimates indicates that kriging standard deviation depends mainly on data configuration and therefore cannot be used as a reliable measure to report estimation precision. In contrast, the conditional standard deviation provided by simulation algorithms represents the estimation uncertainty more reliably because it depends on data values as well. This shows the risks of decision-making based only on kriging estimates without assessment of uncertainties. Further, the accuracy plots and statistic G are used to evaluate, respectively, the accuracy and goodness of probabilistic models of local uncertainty provided by the simulation algorithms. The sis algorithm, having a larger statistic G, performs the modelling uncertainty slightly better than sGs.

The simulation algorithms, which better reproduce the data statistics, are deemed to be more appropriate for the prediction and uncertainty assessment of SOC than the classical kriging method. Further, when the class-specific patterns of spatial continuity of the attribute must be preserved, sis has more advantages than the mathematically simpler algorithm of sGs. The model of uncertainty produced through the set of alternative realisations generated by the simulation algorithms can be used by decision-makers in risk assessment of many processes such as SOC management in Lower Austria.


Manuscript received 31 January 2009, accepted 7 October 2009


Bourennane H, King D, Couturier A, Nicoullaud B, Mary B, Richard G (2007) Uncertainty assessment of soil water content spatial patterns using geostatistical simulations: an empirical comparison of a simulation accounting for single attribute and a simulation accounting for secondary information. Ecological Modelling 205, 323-335. doi:10.1016/j.ecolmodel.2007.02.034

Caers J (2000) Direct sequential indicator simulation. In 'Proceedings of the 6th International Geostatistics Congress'. Cape Town, South Africa. (Eds W Kleingeld, D Krige) pp. 39-48.

Chen F, Kissel DE, West LT, Adkins W (2000) Field-scale mapping of surface soil organic carbon using remotely sensed imagery. Soil Science Society of America Journal 64, 746-753.

Chen F, Kissel DE, West LT, Rickman D, Luvall JC, Adkins W (2005) Mapping surface soil organic carbon for crop fields with remote sensing. Journal of Soil and Water Conservation 60, 51-57.

Chenu C, Le Bissonnais Y, Arrouays D (2000) Organic matter influence on clay wettability and soil aggregate stability. Soil Science Society of America Journal 64, 1479-1486.

Chiles JP, Delfiner P (1999) 'Geostatistics modeling spatial uncertainty.' Wiley Series in Probability and Statistics. (Wiley: New York)

Deutsch CV (1997) Direct assessment of local accuracy and precision. In 'Geostatistics Wollongong '96'. (Eds EY Baafi, NA Schofield) pp. 115-125. (Kluwer Academic Publishing: Dordrecht, The Netherlands)

Deutsch CV, Journel AG (1998) 'GSLIB: Geostatistical software library and user's guide.' 2nd edn. (Oxford University Press: New York)

Englund EJ (1993) Stochastic simulation: environmental applications. In 'Environmental modeling with GIS'. (Eds MF Goodchild, BO Parks, LT Steyaert) pp. 432-437. (Oxford University Press: New York)

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

Goovaerts P (1999) Impact of the simulation algorithm, magnitude of ergodic fluctuations and number of realizations on the spaces of uncertainty of flow properties. Stochastic Environmental Research and Risk Assessment 13, 161-182. doi:10.1007/s004770050037

Goovaerts P (2001) Geostatistical modeling of uncertainty in soil science. Geoderma 103, 3-26. doi:10.1016/S0016-7061(01)00067-2

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

Journel AG (1983) Non-parametric estimation of spatial distributions. Mathematical Geology 15, 445-468.

Juang KW, Chen YS, Lee DY (2004) Using sequential indicator simulation to assess the uncertainty of delineating heavy-metal contaminated soils. Environmental Pollution 127, 229-238. doi:10.1016/j.envpol.2003.07.001

Le Bissonnais Y, Arrouays D (1997) Aggregate stability and assessment of crustability and erodibility. 2. Application to humic loamy soils with various organic carbon content. European Journal of Soil Science 48, 39-48. doi:10.1111/j.1365-2389.1997.tb00183.x

Leuangthong O, McLennan JA, Deutsch CV (2004) Minimum acceptance criteria for geostatistical realizations. Natural Resources Research 13, 131-141. doi:10.1023/B:NARR.0000046916.91703.bb

Mowrer HT (2000) Uncertainty in natural resource decision support systems: sources, interpretation, and importance. Computers and Electronics in Agriculture 27, 139-154. doi:10.1016/S0168-1699(00)00113-7

Mueller TG, Pierce FJ (2003) Soil carbon maps: enhancing spatial estimates with simple terrain attributes at multiple scales. Soil Science Society of America Journal 67, 258-267.

Pachepsky Y, Acock B (1998) Stochastic imaging of soil parameters to assess variability and uncertainty of crop yield responses. Geoderma 85, 213-229. doi:10.1016/S0016-7061(98)00021-4

Remy N (2004) 'S-Gems: Geostatistical earth modeling software: user's manual.' (Stanford University: Stanford, CA)

Rossi RE, Borth PW, Tollefson JJ (1993) Stochastic simulation for characterizing ecological spatial patterns and appraising risk. Ecological Applications 3, 719-735. doi:10.2307/1942103

Schnabel U, Tietje O, Scholz RW (2004) Uncertainty assessment for management of soil contaminants with sparse data. Environmental Management 33, 911-925. doi:10.1007/s00267-003-2971-0

Simbahan GC, Dobermann A, Goovaerts P, Ping J, Haddix ML (2006) Fine-resolution mapping of soil organic carbon based on multivariate secondary data. Geoderma 132, 471-489. doi:10.1016/j.geoderma. 2005.07.001

Tisdall JM, Oades JM (1982) Organic matter and water-stable aggregates in soils. Journal of Soil Science 33, 141-163. doi:10.1111/j.1365-2389. 1982.tb01755.x

Van Meirvenne M, Pannier J, Hofman G, Louwagie G (1996) Regional characterization of the long-term change in soil organic carbon under intensive agriculture. Soil Use and Management 12, 86-94. doi:10.1111/ j.1475-2743.1996.tb00964.x

Zhao Y, Shi X, Yu D, Wang H, Sun W (2005) Uncertainty assessment of spatial patterns of soil organic carbon density using sequential indicator simulation, a case study of Hebei Province, China. Chemosphere 59, 1527-1535. doi:10.1016/j.chemosphere.2005.01.002

Zotl JG (1997) The spa Deutsch-Altenburg and the hydrogeology of the Vienna basin (Austria). Environmental Geology 29, 176-187. doi:10.1007/s002540050116

Masoomeh Delbari (A,C), Willibald Loiskandl (B), and Peyman Afrasiab (A)

(A) Department of Water Engineering, Faculty of Agriculture, University of Zabol, Zabol, Iran. Email: p_afrasiab@yahoo.com

(B) Institute of Hydraulics and Rural Water Management, Department of Water, Atmosphere and Environment, University of Natural Resources and Applied Life Sciences, Vienna, Austria. Email: willibald.loiskandl@boku.ac.at

(C) Corresponding author. Email: mas_delbari@yahoo.com
Table 1. Characteristics of indicator semivariograms for SOC

            Model         effect     Sill     Range   Effective
Quantile    structure     (% (2))   (% (2))    (m)    range (m)

0.05        Spherical      0.021     0.043      87        87
0.10        Exponential   0.0355     0.095      40       120
0.25        Spherical     0.0811    0.2012      60        60
0.50        Spherical     0.1269    0.2548      80        80
0.75        Spherical     0.0924    0.1858      89        89
0.90        Spherical      0.045    0.0934      90        90
0.95        Spherical     0.0161    0.0537     215       215
Gale Copyright: Copyright 2010 Gale, Cengage Learning. All rights reserved.