Printer Friendly

Uncertainty analysis for large-scale prediction of the van Genuchten soil-water retention parameters with pedotransfer functions.


Soil-water retention curves (SWRCs) are the key functions required in hydrological, environmental and ecological modelling. However, direct measurement of SWRCs is costly, time-consuming and labour-intensive, especially considering that a large number of SWRCs is needed for the simulation of the large-scale water flow and solute transport processes. Therefore, pedotransfer functions (PTFs) have been developed as an alternative, indirect method to predict SWRCs from more easily measurable, basic soil properties such as particle-size distribution, bulk density and organic matter content (e.g. Vcreecken et al. 1989; Tomasella et al. 2003; Liao et al. 2011). Although PTFs have obvious advantages, the uncertainty associated with the parameter estimates of PTFs cannot be ignored. Evaluation of the uncertainty can be useful for improving the parameter estimation in hydrological models.

According to Sivapalan (2003), Chirico et al. (2007) and Deng et al. (2009), the overall uncertainty in parameter prediction of PTFs can be identified by analysing the propagation of three major sources: (i) the model structure, (ii) the model parameters, and (iii) the input variables. The uncertainty in the model structure is explained by the different estimations yielded by different PTFs that implement different conceptualisations of physical or statistical realities. Because PTFs are developed by statistical methods such as multiple regression and artificial neural network (ANN) models, the accuracy of PTFs outside of their development dataset is generally unknown. Guber et al. (2006, 2009) suggested that the model structure uncertainty could be addressed by multi-model ensemble prediction methods such as simple averaging, super-ensemble with singular-value decomposition, Bayesian model averaging, and the corrected Akaike information criterion-based multi-model forecasting.

The uncertainty in model parameters, also called the intrinsic uncertainty, is another source of uncertainty in the parameter prediction of PTFs. The nonparametric bootstrap method has been used to simulate this uncertainty in previous studies. For example, Schaap et al. (1998) and Ye et al. (2007) used bootstrap method to investigate uncertainty in soil hydraulic properties predicted by ANN-based PTFs. Deng et al. (2009) further investigated the effect of bootstrap parameter estimation uncertainty on unsaturated flow modelling based on the study by Ye et al. (2007). Santra et al. (2009) coupled the multiple linear regression with bootstrap method to develop PTFs for reducing the uncertainty in predicted hydraulic parameters.

The third source of parameter prediction uncertainty of PTFs is the uncertainty in the input variables (e.g. basic soil properties), which is due to measurement errors at the local scales and spatial interpolation errors at the field and larger scales. This uncertainty has been assessed by stochastic methods such as Monte Carlo modelling, Latin hypercube sampling (LHS), and sequential Gaussian simulation. For example, Minasny and McBratney (2002) used the LFIS method to investigate the effect of measurement errors of input variables on the uncertainty of PTF predictions and soil-water balance simulation. In addition, Minasny et al. (1999) used Monte Carlo modelling to analyse the uncertainty in the input variables, and showed that such uncertainty dominates the intrinsic uncertainty of PTFs. However, Chirico et al. (2007, 2010) used a sequential Gaussian simulation approach to address the input uncertainty of PTFs and found that such uncertainty is smaller than the uncertainty in the model structure and parameters.

A few studies have conducted uncertainty analysis for the parameter estimates of PTFs at scales of point (e.g. Minasny and McBratney 2002), field (e.g. Ye et al. 2007; Deng et al. 2009) and hillslope (e.g. Chirico et al. 2007, 2010); however, to our knowledge little research has been conducted to investigate prediction uncertainty of soil hydraulic properties at the larger scale. Such analysis could help to improve prediction of hydraulic parameters and simulation of large-scale hydrological processes. The intrinsic uncertainty of PTFs can be evaluated by using bootstrap technique as discussed earlier. The input variables of PTFs, such as basic soil properties, should be spatially interpolated (e.g. using kriging and co-kriging) based on a limited number of observations. Thus, the input uncertainty of PTFs, measured by the (co)kriging variance of the predicted basic soil properties (Isaaks and Srivastava 1989), can be assessed with the LHS method.

The objectives of this research were to: (i) quantify the intrinsic uncertainty of PTFs for predicting the van Genuchten soil-water retention parameters at large scale; and (if) evaluate the relative contributions of the intrinsic and input uncertainties to the overall parameter prediction uncertainty of PTFs.

Materials and methods

Study area

This study was conducted in the agricultural region of Pingdu City. The study site (36[degrees]28'-37[degrees]02'N, 119[degrees]31'-120[degrees]19'E), with an area of 3166.54 [km.sup.2], is in the eastern part of Shandong Province, China (Fig. 1). The annual mean temperature is 12.3[degrees]C, with an average annual precipitation of 707 mm, >70% of which occurs from July to September. The landscape of the area is generally sloping. Elevation ranges from >150m above mean sea level (a.m.s.l.) in the north-east part to <50 m a.m.s.l. in the south and central parts of the study area. The typical land-use practices for this region arc winter wheat and summer maize cropping. According to the World Reference Base for Soil Resources (WRB) (FAO/ISRIC/ISSS 1998), soils in the study area can be classified as Luvisols, Fluvisols, Vertisols, Lixisols and Solonchaks. The first three of these soil types cover >96% of the area.

Soil sampling and analyses

Based on the number of sampling points being proportional to each soil type, 58 points were randomly selected--27, 16 and 15 points in the areas of Luvisols, Fluvisols and Vertisols, respectively (Fig. 1). At each point, undisturbed soil cores (100 [cm.sup.3] in volume) were collected from topsoils (0-15 cm) in November 2007. These soil cores were brought to the laboratory for analysis of soil-water retention characteristics and bulk density (BD). Five replicates per site were taken, and the measured soil-water retention data and BD were expressed as averages. At the same time, five surface soil columns were collected by a hand auger at each sampling site. These five soil columns were fully mixed, and ~ 1 kg of the mixed soil was packed in a bag and brought to the laboratory for determining soil organic matter (OM) content and particle-size distribution (PSD). The coordinates of all sampling locations were recorded using a global positioning system (GPS) receiver. For each sampling point, the soil-water retention data were measured using a pressure membrane apparatus at suctions of 2.5, 6, 10, 30, 100, 300 and 1500kPa following the method by Klute (1986). The three-parameter van Genuchten model was fitted to the seven measured retention data (van Genuchten 1980):

[theta](h) = [[theta].sub.s]/[(1 + [[absolute value of [alpha]h].sup.n]).sup.(1-1/n)] (1)

where [theta](h) is the volumetric water content ([cm.sup.3]/[cm.sup.3]) at a specified matric potential head, h (cm); 0S is the soil saturated volumetric water content ([cm.sup.3]/[cm.sup.3]); and a (1/cm) and n are shape parameters of the SWRCs. The unknown van Genuchten parameters ([[theta].sub.s], [alpha] and n) were estimated with the nonlinear least-squares optimisation procedure in the RETC program (van Genuchten et al. 1991) from measured retention data. After measuring soil-water retention data, the same soil cores were then oven-dried at 105[degrees]C for 24 h to determine the BD. The soil OM content was estimated from the organic carbon content obtained by the Walkley-Black wet oxidation method (Nelson et al. 1996), using a constant 1.724 for transformation. The contents of sand (0.05-2 mm), silt (0.002-0.05 mm) and clay (0-0.002 mm) particles were measured with the pipette method. The soil texture was determined by the USDA classification. Table 1 summarises the basic soil properties in this study.

Remote-sensing data

Remote-sensing data can be used as auxiliary variables for spatial prediction of basic soil properties by using co-kriging method. For example, Sullivan et al. (2005) used co-kriging with IK.ONOS multispectral data to predict field-scale surface total soil carbon and clay content. Wu et al. (2009) applied co-kriging with Landsat Enhanced Thematic Mapper (ETM) data to estimate large-scale surface soil OM content. In the present study, a Landsat ETM image acquired on 7 November 2007 was used to predict basic soil properties. Nearly 80% of the pixels with a normalised difference vegetation index value <0.1 in the study area are classified as bare soil. ETM has a swath width of 185 km and offers seven spectral bands: visual (bands 1-3) (450-690 nm), near-infrared (band 4) (760-900 nm) and mid-infrared (bands 5 and 7) (1550-2350 nm) with a spatial resolution of 30 m, and a thermal infrared band (band 6) with a spectral range of 10400-12 500 nm and a spatial resolution of 60 m. The revisit time of Landsat ETM is 16 days. Since band 6 is used mainly for surface temperature inversion, the digital numbers (DNs) of the other six bands (bands 1-5 and 7) were evaluated to obtain the best one for spatial estimation of basic soil properties (Wu et al. 2009). To make good use of remote-sensing data, 402 pixels were sampled in this study. Among these 402 pixels, 344 pixels were sampled using a grid-based sampling scheme with 3-km grid spacing, and the other 58 pixels were sampled based on the corresponding locations of the 58 soil sampling points (Fig. 1).

Multiple stepwise regression

Multiple stepwise regression (MSR) was conducted to predict the van Genuchten parameters using BD, and contents of OM, sand, silt and clay of the 58 samples. Values of a were natural-logarithm-transformed to follow the normal distribution. In MSR analysis, important predictor variables were identified at P=0.05 (Norusis 1994). The MSR model was fitted using the regress function of MATLAB software (The MathWorks Inc., Natick, MA, USA).

Geostatistical analysis

Geostatistical techniques were used to predict the spatial distribution of basic soil properties with the assistance of remote sensing data. Kriging is a basic geostatistical technique that provides the best linear unbiased estimate for a spatially dependent variable. This method, however, may produce unsatisfying results if the sample size is limited (Zhu and Lin 2010). Co-kriging is another geostatistical method that extends kriging of a primary variable to secondary variables based on their correlations with the primary variable. It has been demonstrated that co-kriging is superior to kriging for minimising the estimation variance (Istok et al. 1993; Wu et al. 2009).

(Co)kriging uses the (cross-)semivariograms to quantify the spatial patterns of spatially dependent soil properties. The estimator for the (cross-)semivariogram is:

[[gamma].sub.ij](h) = 1/2N{h) [N(h).summation over (k=1)] {[[z.sub.i]([x.sub.k] + h) - [z.sub.i](x)][[z.sub.j]([x.sub.k] + h)-[z.sub.j] (x)]} (2)

where [[gamma].sub.ij] is the semivariance (when i = j) with respect to random variable [z.sub.i]; h is the separation distance; and N(h) is the number of pairs of [z.sub.i]([x.sub.k]) and [z.sub.j]([x.sub.k]) in a given lagged distance interval of (h + dh). When i [not equal to] j, [[gamma].sub.ij] is the cross-semivariogram, which is a function of h (Yates and Warrick 1987). Four models (spherical, exponential, linear and Gaussian) were used to describe the (cross-)semivariogram, and the best-fitted models were selected with the largest coefficient of determination ([R.sup.2]) between model predicted and the measured soil property variances.

All geostatistical computations were performed using GS + version 9.0 (Gamma Design Software LLC, Plainwell, MI, USA). ArcGIS 9.3 (Esri, Redlands, CA, USA) was used for mapping and spatial analysis of study field data.

Bootstrap method

The bootstrap method, first introduced by Efron and Tibshirani (1993), is a nonparametric technique for estimating the probability distribution of any dataset when the distributional form is unknown, which is the most common situation. In this study, the spatial distribution of the van Genuchten parameters was predicted by MSR-bascd PTFs developed from the 58 samples, which is only one realisation of the original population of soil properties. Hence, developing PTFs on a different realisation of the population would result in slightly different predictions of [[theta].sub.s], [alpha] and n. To predict the probability distribution of the van Genuchten parameters, the bootstrap datasets containing 58 samples can be generated by repeated random resampling with replacement of the original dataset of 58 samples. Since resampling is conducted with replacement, each bootstrap realisation contains ~63% of the original dataset, while the remaining 37% of the dataset is left out. The multiple realisations of the dataset by the bootstrap technique represent the intrinsic uncertainty of PTFs. The number of bootstrap realisations is critical for evaluating the intrinsic uncertainty. Efron and Tibshirani (1993) suggested that the number of bootstrap realisations be 50-200. In the study by Schaap et al. (1998), 100 realisations were used. Therefore, 100 realisations of random resampling were also used in this study.

Latin hypercube sampling

The LHS is a stratified sampling method that allows efficient prediction of any statistical distribution (McKay et al. 1979; Iman et al. 1980). The standard LHS consists of the following steps (Janssen et al. 1992): (7) equi-probable subdivision of the probability distribution of each parameter; (2) stratified sampling within each interval according to the distribution; and (3) random pairing in case of non-correlated parameters or in case of correlation combination of the sampled values in a correlated fashion. In this study, the LHS approach was used to evaluate the input uncertainty of PTFs due to the spatial variability of the input variables (basic soil properties). Before using LHS, the distribution of the (co)kriged variables should be assumed to be Gaussian, which is deemed reasonable (e.g. Cressie 1991; Deutsch and Joumel 1998). On this basis, the (co) kriged value and (co)kriging variance of the distribution are used for random sampling. When using LHS, the number of LHS realisations was decided by the number of parameters to be sampled. Previous studies suggested that the number of LHS realisations should be more than 4/3 of the number of input parameters to be sampled (Janssen et al. 1992; Klier 2012). Therefore, 20 LHS realisations were representatively sampled from the probability distributions of the input variables in this study.


Development of PTFs

Equation 1 was fitted to the 58 measured soil-water retention datasets. The means and standard deviations of the fitted van Genuchten parameters are listed in Table 2. Mean [R.sup.2] values for all of the textural groups were >0.94, indicating that the Eqn 1 can well describe the SWRCs for the 58 samples of the study area. No clear pattern wasobserved for [[theta].sub.s], [alpha] and [R.sup.2] among different soil textures, whereas a decreasing trend was observed for n as the texture varied from moderately coarse to moderately fine (T able 2). The parameter n controls the slope of the SWRCs. A high value of n results in a sand-like SWRC, where a large change in soil moisture content occurs with a small change in suction above the air-entry value.

In order to derive the PTFs, Pearson correlation between the basic soil properties and soil-water retention parameters was conducted (Table 3). The [[theta].sub.s] was significantly (P< 0.05) correlated with the sand content and BD of soils. The ln([alpha]) was significantly (P<0.05) correlated with soil BD. The n was significantly (P<0.05) correlated with sand, clay and OM contents of soil samples. The PTFs for these van Genuchten parameters were derived from the 58 soil samples using MSR (Table 4). The performance of all PTFs was evaluated by their [R.sup.2] values. The BD was the only predicting variable for [[theta].sub.s] and In(cc). Their PTFs explained ~43% and 10% of the variance in [[theta].sub.s] and ln([alpha]), respectively. Clay and OM contents were the predicting variables for n. The regression equation explained ~34% of the variance in n.

Spatial estimation of the van Genuchten function parameters

Spatial distributions of the basic soil properties are needed to predict spatial variations of van Genuchten parameters. The remote-sensing data were treated as auxiliary variables for predicting the spatial distributions of basic soil properties. Pearson correlation was conducted between the basic soil properties and Landsat ETM DNs of six bands (bands 1-5 and band 7) (Table 5). Positive correlations (P<0.05) were observed between sand content and DNs of bands 2, 3,4, 5, and 7, whereas silt had negative correlations (P<0.05) with DNs of bands 2, 3, 4, 5, and 7. Clay and OM contents had negative correlations (P< 0.05) with DNs of all bands. The DNs of band 7 and band 2 had the largest absolute correlation coefficients with texture and OM, respectively, indicating that DNs of band 7 and band 2 could be used as auxiliary variables to predict soil texture and OM content by co-kriging, respectively. However, none of the DNs was significantly correlated with BD. Semivariogram analysis and kriging method can be used to predict the spatial variability of soil BD. Although kriging would produce too smooth a map for BD, it should be acceptable to use predicted BD to estimate van Genuchten parameters because of the very low coefficient of variation (CV) (0.05) for BD (Table 1).

The semivariogram and cross-semivariogram analyses were used to analyse the spatial variability of basic soil properties. The cross-semivariogram of sand content (or fraction) was fitted well by the Gaussian model, with a nugget: sill ratio of 0.1%, and cross-semivariograms of silt, clay and OM contents were well fitted by the spherical model, with nugget: sill ratios of 0.2%, 0.2% and 9.2%, respectively. The semivariogram of BD was fitted well by the exponential model, with a nugget: sill ratio of 16.4%. Spatial patterns of soil texture and OM content were predicted by co-kriging using auxiliary variables of DNs of band 7 and band 2, respectively; the spatial pattern of BD was predicted by kriging. Estimated spatial patterns of sand content and BD are shown in Fig. 2 as examples. The soil sand content was higher in the north-eastern parts of the study area (Fig. 2a). This could be explained by the real distribution pattern of soil types in the study area (Fig. 1). Sand content in Luvisols, formed from coarse-grained residual and slope deposits, is higher than that in Vertisols, formed from fine-textured lacustrine deposits. Therefore, soil sand content in the north-eastern parts of the study area, with most of the Luvisols distribution, was much higher than in the area with Vertisols. Areas with large standard deviations of predicted sand content and BD corresponded to the fringe regions of Pingdu City. These areas also had few soil sampling sites (Fig. 1). The MSR-based PTFs were used to predict the spatial distribution of van Genuchten parameters from spatial distributions of basic soil properties. Estimated spatial patterns of [[theta].sub.s] and n are shown in Fig. 3 as examples.

Intrinsic uncertainty of PTFs

Effects of intrinsic uncertainty of PTFs on large-scale estimation of the van Genuchten soil-water retention parameters were further investigated using bootstrap method. One hundred realisations of random resampling were used in this study. However, no significant correlation was observed between the randomly sampled basic soil properties and ln([alpha]) in 25 of 100 realisations. Thus, only 75 realisations were used for assessing the intrinsic uncertainty of PTF-predicted ln([alpha]). One hundred PTFs have been developed for estimating [[theta].sub.s] and n based on 100 realisations, and 75 PTFs have been derived for predicting ln([alpha]) based on 75 realisations. These PTFs were used to transfer the spatial distributions of basic soil properties to the spatial distributions of van Genuchten parameters. One hundred predicted [[theta].sub.s] and n values and 75 estimated ln([alpha]) values for the 58 samples passed the Kolmogorov-Smimov test at P=0.05 for normal distribution. Figure 4a, b presents the histograms for [[theta].sub.s] and n for a loamy soil sample, indicating that they fit well to the normal distribution. Therefore, the intrinsic uncertainty could be described by the mean and standard deviation of a normal distribution, predicted from the multiple bootstrap realisations.

The CVs of the van Genuchten parameters were calculated from the bootstrap realisations for assessing the intrinsic uncertainty of PTFs. Means and standard deviations of CVs of [[theta].sub.s], ln([alpha]) and n for these 58 samples are shown in Fig. 5a. Because the values of the calculated CVs for ln([alpha]) were negative, the absolute CVs of ln([alpha]) are shown. The CV of ln([alpha]) was greater than those of [[theta].sub.s] and n, indicating that the intrinsic uncertainty of PTFs for predicting ln([alpha]) was greater than those for estimating [[theta].sub.s] and n. Comparisons between fitted and predicted van Genuchten parameters for the 58 samples were made based on bootstrap realisations (Fig. 6). The variability of PTF-estimated parameters was smaller than that of observations. However, the variability of the fitted van Genuchten parameters was partly captured by the MSR-based PTFs when the intrinsic uncertainty of PTFs was considered. For [[theta].sub.s], ln([alpha]) and n, 36%, 29% and 47% of measurements were encompassed within the corresponding error bars (95% confidence intervals of the means), respectively (Fig. 6a-c).

Input uncertainty of PTFs

The input variables of PTFs are only the kriging and co-kriging predictions of the basic soil properties. Therefore, the smoothing effect of kriging and co-kriging should result in smaller variability of the PTF predictions than that of the observations. Standard deviations of sand and BD were almost 1/6 and 1/13 of the corresponding mean values, respectively (Fig. 2), indicating that the uncertainty in the input variables cannot be ignored. The LHS method was used to assess the input uncertainty of PTFs. Twenty LHS realisations were representatively sampled from the probability distributions of the input variables in this study. A small number of random sand values were >100%, and clay, OM, and BD values negative (data not shown). However, these values are not physically meaningful. The minimum and maximum values of input variables for the 58 soil samples, therefore, were used as the bounds for sampling. At the same time, linear correlation, rather than the spatial one, was considered between the input variables during the sampling process.

In order to evaluate both intrinsic and input uncertainties of PTFs in predicting van Genuchten soil-water retention parameters, 2000 realisations of [[theta].sub.s] and n were derived based on 20 LHS samples of the input variables and 100 bootstrap realisations. Similarly, 1500 realisations of ln([alpha]) were generated based on 20 LHS samples and 75 bootstrap realisations. Two thousand predicted [[theta].sub.s] and n values and 1500 estimated ln([alpha]) values for the 58 samples also passed the Kolmogorov-Smimov test at P=0.05 for normal distribution. Figure 4c, d shows the histograms for [[theta].sub.s] and n for a loamy soil, indicating that they fit well to the normal distribution. Therefore, both the intrinsic and input uncertainties could be described by the mean and standard deviation of a normal distribution predicted from the multiple bootstrap and LHS realisations. The mean values and standard deviations of CVs of the van Genuchten parameters for the 58 samples arc shown in Fig. 5b. The absolute CV of [[theta].sub.s] was similar to that of ln([alpha]), but greater than that of n. This indicates that intrinsic and input uncertainties of PTFs for predicting [[theta].sub.s] and ln([alpha]) were greater than those for estimating n. Comparisons between fitted and predicted van Genuchten parameters for the 58 samples were made based on both bootstrap and LHS realisations (Fig. 7). The variability of the observed van Genuchten parameters was well captured by the MSR-based PTFs when both intrinsic and input uncertainties of PTFs were considered. For the [[theta].sub.s], ln([alpha]) and n, 86%, 66% and 88% of measurements were encompassed within the corresponding error bars, respectively (Fig. 7a-c).


PTFs for soil-water retention parameters

The regression equations for predicting ln([alpha]) and n have lower [R.sup.2] values than those for estimating [[theta].sub.s], indicating that the predictive performances for ln([alpha]) and n are poorer than that for [[theta].sub.s]. This is consistent with other studies (e.g. Wosten et al. 1999; Pachepsky and Rawls 2004; Santra et al. 2009). For example, Wosten et al. (1999) used the 5521 hydraulic properties from 1777 soil profiles of the European database HYPRES for deriving PTFs, and showed poor performances of ln([alpha]) ([R.sup.2] = 0.20) and In (n - 1) ([R.sup.2] = 0.54) compared with [[theta].sub.s] ([R.sup.2] = 0.76). A possible explanation is the fact that the fitting parameters a and n are not real soil properties and their values are very sensitive to the fitting method and criteria (Wosten and van Genuchten 1988). In addition, performances of our PTFs are worse than reported by Wosten et al. (1999). This is probably because, in this study, the basic soil properties used for predicting the soil-water retention characteristics are very limited. Another reason could be that PSD is not detailed enough to describe the individual contribution of soil particles to the soil-water retention characteristics. This suggests that caution is needed when applying PTFs outside of their development region.

The PTFs can be improved by introducing additional predictors, such as porosity, particle density, soil chemical properties, soil structure, mineralogy and pedality (Rawls et al. 1991). Additional particle size fractions could also be applied to improve the performance of the PTFs. However, determining additional predictors would take too much time and energy, thereby going against the original intention of the PTF development. In particular, introducing additional predictors will also bring new uncertainty to the PTFs.

Correlations between basic soil properties and remote sensing data

As expected, soil texture has stronger absolute correlations with near-infrared bands (band 4-5 and band 7) than with visible bands (band 1-3) of ETM imagery. This is related to the fact that spectral signatures of the soil texture are typically dominated by the near-infrared spectra (Stenberg et al. 2010). Sand and clay contents have better correlation with the DN values than does silt content. This finding is consistent with that of Wetterlind and Stenberg (2010), who indicated that silt was worse predicted because the silt fraction can be assumed to be a mixture of sand and clay minerals and is therefore difficult to distinguish in the visible and the near-infrared regions. In this study, the DN of band 7 has a significant correlation with soil texture in a large area with different soil types. This is attributed to the fact that the band 7 is referred to as short-wave infrared (SWIR), which is sensitive to soil-available water content. The available water content depends greatly on the soil texture. In addition, soil texture has direct impact on SWIR reflectance, as incoming radiation is scattered differently by coarse particles compared with fine particles (Lagacherie et al. 2008).

In our study, the OM has a stronger correlation with visible bands (especially band 2) in ETM imagery. The results are similar to those of Wu et al. (2009). This may be attributed to the fact that soil OM has broad absorption bands in the visible region that are dominated by chromophores and the darkness of organic matter (Stenberg et al. 2010). Mortimore et al. (2004) also found that adsorptions in the visible region occur mostly from electronic excitations in the chromophores that are often present in various organic compounds. In addition, van der Meer et al. (2001) proposed that if OM in soils falls below 2%, it has only a minimal effect on the reflectance property. However, in this study, the mean OM is 1.36% with a standard deviation of 0.38%. Other studies have also found significant correlations between soil spectral properties and OM when most OM values are <2% (Scharf et al. 2002; Wu et al. 2009). One possible reason is that influences of other soil properties such as mineral composition and water content on soil spectral properties are relatively weak. Thus, soil OM is still the main factor influencing soil spectral properties.

Bulk density was not found to have significant correlations with the DN values. The reason for this issue is still unclear but may be related to fact that the spatial resolution of Landsat ETM is relatively coarse. Moreira et al. (2009) found that the near-infrared reflectance spectroscopy (NIRS) method, which offers the advantages of rapidity, low-cost and reproducibility, produces good estimates of soil BD. Therefore, the NIRS method may be more useful for soil BD assessment. This will be investigated in a future work.

Integrated analysis of PTF intrinsic and input uncertainties

Comparison of the fitted and predicted van Genuchten parameters by the MSR-based PTFs shows root mean-squared errors (RMSEs) of 0.046 [cm.sup.3]/[cm.sup.3], 0.650 [cm.sup.-1] and 0.091 for [[theta].sub.s], ln([alpha]) and n of the 58 soil samples, respectively. This suggests that the performances of the PTFs are similar to those reported by Merdun et al. (2006) and Santra et al. (2009) even though we only have 58 soil samples. However, uncertainties associated the spatial prediction of van Genuchten parameters cannot be neglected. As depicted in Fig. 6, although the intrinsic uncertainty of PTFs is considered, a large portion of measurements is still outside the intervals. This indicates that the relative contribution of the PTF intrinsic uncertainty to the total prediction uncertainty is relatively limited.

Comparison of Figs 6 and 7 shows that the parameter prediction uncertainty increases significantly when both the intrinsic and input uncertainties of PTFs are taken into account for uncertainty assessment. In this case, most of these observations, except some extreme values, are within ranges of the predicted values. This implies that the PTF input uncertainty is more important than the PTF intrinsic uncertainty. The main reason is obviously that the kriged soil BD, co-kriged PSD and OM cannot Hilly describe the spatial heterogeneity reflected in the basic soil properties, because of the smoothing effect of the kriging and co-kriging methods. Zhu and Lin (2010) proposed that sample size can significantly affect the interpolation accuracy and showed that soil properties with 60-100 points can be reliably interpolated using geostatistical method. Webster and Oliver (1992) also found that >100 data points are needed to predict a regular variogram. As a result, the PTF input uncertainty can be reduced efficiently with a relatively large number of sampling points (e.g. 100 data points). In addition, applying more powerful spatial interpolation approaches, such as the multifractal inverse distance weighted interpolation and high-accuracy surface modelling, is also useful for improving the spatial prediction of the basic soil properties. These two interpolation methods have been demonstrated to be much more accurate than the geostatistical method (Lima et al. 2008; Shi et al. 2009).

Some studies have found that PTF input uncertainty dominates over PTF intrinsic uncertainty (Minasny et al. 1999; Matula and Spongrova 2007; Deng et al. 2009). However, other studies have shown that PTF parameter estimation uncertainty is primarily caused by PTF intrinsic uncertainty other than the PTF input uncertainty (To and Kay 2005; Chirico et al. 2007, 2010). The various studies suggest that the answer to the question of 'which type of uncertainty is more important' depends on the amount of data available for the development of PTFs. Wbsten et al. (2001) and van der Keur and Iversen (2006) also concluded that if PTFs are calibrated adequately, then PTF intrinsic uncertainty is usually smaller than PTF input uncertainty.

Although the LHS is efficient for assessing the input uncertainty in this study, a drawback of this method is that its highly structured form makes it difficult to increase the size of an already generated sample (Sallaberry et al. 2008). In previous studies, the sequential Gaussian simulation (SGS) has also been applied to evaluate the input uncertainty of PTFs (Chirico et al. 2007, 2010). Simuta-Champo and Herrera-Zamarron (2010) compared the LHS and SGS methods and found that it is not possible to accumulate the realisations in a sequential manner for the LHS, whereas it is possible to accumulate them for the SGS. In future, comparison of these two methods to assess the PTF input uncertainty should be considered for improving estimation of the van Genuchtcn parameters.


This study was conducted to assess the uncertainty in spatially estimated van Genuchten soil-water retention parameters at the large scale. Results suggest that the MSR-based PTFs cannot capture the variability of the fitted van Genuchten parameters without considering the parameter estimation uncertainty. The investigation of uncertainty contributions reveals that the errors in PTF predictions arc more affected by the input uncertainty than by the intrinsic uncertainty of PTFs. Therefore, reducing the PTF input uncertainty is particularly important in the prediction of van Genuchten parameters. Complete spatial characterisation of basic soil properties and cautious positioning of the sampling points would bring about an optimal scheme to spatially characterise soil-water retention parameters. Future work will be conducted to test the propagation of the uncertainties associated with the spatial prediction of soil-water retention parameters in a large-scale distributed hydrological model, for improving the predictions of spatial soil-water distribution over the study area.

Received 3 January 2013, accepted 23 March 2014, published online 16 June 2014


This study was financially supported by the National Natural Science Foundation of China (Grant Nos. 41301234 and 41171183) and Key 135' Project of Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences (NIGLAS2012135005).


Chirico GB, Medina H, Romano N (2007) Uncertainty in predicting soil hydraulic properties at the hillslope scale with indirect methods. Journal of Hydrology 334, 405-422. doi:10.1016/j.jhydrol.2006.10.024

Chirico GB, Medina H, Romano N (2010) Functional evaluation of PTF prediction uncertainty: An application at hillslope scale. Geoderma 155, 193-202. doi: 10.1016/j.geoderma.2009.06.008

Cressie N (1991) 'Statistics of spatial data.' (John Wiley: New York) Deng HL, Ye M, Schaap MG, Khaleel R (2009) Quantification of uncertainty in pedotransfer function-based parameter estimation for unsaturated flow modeling. Water Resources Research 45, W04409. doi: 10.1029/2008WR007477

Deutsch CV, Joumel AG (1998) 'GSLIB-geostatistical software library and user's guide.' (Oxford University Press: Oxford, UK)

Efron B, Tibshirani RJ (1993) 'An introduction to the bootstrap.' Monographs on Statistics and Applied Probability, Vol. 57. (Chapman and Hall: London)

FAO/ISRIC/ISSS (1998) 'World Reference Base for Soil Resources.' World Soil Resources Reports. (FAO: Rome)

Guber AK, Pachepsky YaA, van Genuchten MTh, Rawls WJ, Simunek J, Jacques D, Nicholson TJ, Cady RE (2006) Field-scale water flow simulations using ensembles of pedotransfer functions for soil water retention. Vadose Zone Journal 5, 234-247. doi:10.2136/vzj2005.0111

Guber AK, Pachepsky YaA, van Genuchten MTh, Simunek J, Jacques D, Nemes A, Nicholson TJ, Cady RE (2009) Multimodel simulation of water flow in a field soil using pedotransfer functions. Vadose Zone Journals, 1-10. doi:10.2136/vzj2007.0144

Iman RL, Davenport JM, Zeigler DK (1980) 'Latin hypercube sampling. Program user's guide.' SAND79-1473. (Sandia National Laboratories: Albuquerque, NM, USA)

Isaaks EH, Srivastava RM (1989) 'Applied geostatistics.' (Oxford University Press: New York)

Istok JD, Smyth JD, Flint AL (1993) Multivariate geostatistical analysis of ground-water contaminant: A case history. Ground Water 31, 63-74. doi: 10.1111/j.1745-6584.1993.tb00829.x

Janssen PHM, Heuberger PSC, Sanders R (1992) UNCSAM 1.1: a software package for sensitivity and uncertainty analysis. Report No. 959101004. National Institute of Public Health and Environmental Protection, Bilthoven, The Netherlands.

Klier C (2012) Use of an uncertainty analysis for genome-scale models as a prediction tool for microbial growth processes in subsurface environments. Environmental Science & Technology 46, 2790-2798. doi: 10.1021/es203461 u

Klute A (1986) Water retention: laboratory methods. In 'Agronomy monograph. No. 9'. (Ed. VT Holliday) pp. 635-662. (American Society of Agronomy and Soil Science Society of America: Madison, WI, USA)

Lagacherie P, Baret F, Feret JB, MadeiraNetto J, Robbez-Masson JM (2008) Estimation of soil clay and calcium carbonate using laboratory, field, and airborne hyperspectral measurements. Remote Sensing of Environment 112, 825-835. doi:10.1016/j.rse.2007.06.014

Liao K, Xu S, Wu J, Ji S, Lin Q (2011) Assessing soil water retention characteristics and their spatial variability using pedotransfer functions. Pedosphere 21, 413-422. doi: 10.1016/S1002-0160(11)60143-4

Lima A, Plant JA, De Vivo B, Tarvainen T, Albanese S, Cicchella D (2008) Interpolation methods for geochemical maps: a comparative study using arsenic data from European stream waters. Geochemistry: Exploration, Environment, Analysis 8, 41-48. doi: 10.1144/1467-7873/07-146

Matula S, Spongrova K (2007) Pedotransfer function application for estimation of soil hydrophysical properties using parametric methods. Plant, Soil and Environment 53, 149-157.

McKay MD, Conover WJ, Beckman RJ (1979) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21, 239-245.

Merdun H, Qinar O, Meral R, Apan M (2006) Comparison of artificial neural network and regression pedotransfer functions for prediction of soil water retention and saturated hydraulic conductivity. Soil & Tillage Research 90. 108-116. doi: 10.1016/j.still.2005.08.011

Minasny B, McBratney AB (2002) Uncertainty analysis for pedotransfer functions. European Journal of Soil Science 53, 417-429. doi: 10.1046/j.1365-2389.2002.00452.x

Minasny B, McBratney AB, Bristow KL (1999) Comparison of different approaches to the development of pedotransfer functions for water retention curves. Geoderma 93, 225-253. doi: 10.1016/S0016-7061(99)00061-0

Moreira CS, Brunet D, Vemeyre L, Sa SMO, Galdos MV, Cerri CC, Bemoux M (2009) Near infrared spectroscopy for soil bulk density assessment. European Journal of Soil Science 60, 785-791. doi: 10.1111/j. 1365-2389.2009.01170.x

Mortimore JL, Marshall LJR, Almond MJ, Hollins P, Matthews W (2004) Analysis of red and yellow ochre samples from Clearwell Caves and Catalhoyuk by vibrational spectroscopy and other techniques. Spectrochimica Acta. Part A: Molecular and Biomolecular Spectroscopy 60, 1179-1188. doi: 10.1016/j.saa.2003.08.002

Nelson DW, Sommers LE, Sparks D, Page A, Helmke P (1996) Total carbon, organic carbon, and organic matter. In 'Methods of soil analysis. Part 3. Chemical methods', pp. 961-1010. (Soil Science Society of America: Madison, WI, USA)

Norusis JM (1994) 'SPSS professional statistics 6.1.' (SPSS Inc.: Chicago, IL, USA)

Pachepsky YA, Rawls WJ (2004) 'Development of pedotransfer functions in soil hydrology.' (Elsevier: Amsterdam)

Rawls WJ, Gish TJ, Brakensiek DL (1991) Estimating soil water retention from soil physical properties and characteristics. Advances in Soil Science 16, 213-234. doi: 10.1007/978-1-4612-3144-8_5

Sallaberry CJ, Helton JC, Hora SC (2008) Extension of Latin hypercube samples with correlated variables. Reliability Engineering & System Safety 93, 1047-1059. doi:10.1016/j.ress.2007.04.005

Santra P, Sahoo RN, Das BS (2009) Estimation of soil hydraulic properties using proximal spectral reflectance in visible, near-infrared, and shortwave-infrared (VIS-NIR-SWIR) region. Geoderma 152, 338-349. doi:10.1016/j.geoderma.2009.07.001

Schaap MG, Leij FJ, van Genuchten MTh (1998) Neural network analysis for hierarchical prediction of soil water retention and saturated hydraulic conductivity. Soil Science Society of America Journal 62, 847-855. doi:10.2136/sssaj1998.03615995006200040001x

Scharf PC, Schmidt JP, Kitchen NR, Sudduth KA, Hong SY, Lory JA, Davis JG (2002) Remote sensing for nitrogen management. Journal of Soil and Water Conservation 57, 518-524.

Shi W, Liu J, Du Z, Song Y, Chen C, Yue T (2009) Surface modelling of soil pH. Geoderma 150, 113-119. doi:10.1016/j.geoderma.2009.01.020

Simuta-Champo R, Herrera-Zamarron GS (2010) Convergence analysis for Latin-hypercube lattice-sample selection strategies for 3D correlated random hydraulic conductivity fields. Geofisica Internacional 49, 131-140.

Sivapalan M (2003) IAHS decade on predictions in ungauged basins (PUB). Hydrological Sciences Journal 48, 857-880. doi: 10.1623/hysj.48.6.857.51421

Stenberg B, Viscarra Rossel RA, Mouazen AM, Wetterlind J (2010) Visible and near infrared spectroscopy in soil science. Advances in Agronomy 107, 163-215. doi: 10.1016/S0065-2113(10)07005-7

Sullivan DG, Shaw JN, Rickman D (2005) IKONOS imagery to estimate surface

soil property variability in two Alabama physiographies. Soil Science Society of America Journal 69, 1789-1798. doi: 10.2136/sssaj2005.0071

To J, Kay BD (2005) Variation in penetrometer resistance with soil properties: the contribution of effective stress and implications for pedotransfer functions. Geoderma 126, 261-276. doi: 10.1016/j.geoderma.2004.08.006

Tomasella J, Pachepsky YA, Crestana S, Rawls WJ (2003) Comparison of two techniques to develop pedotransfer functions for water retention. Soil Science Society of America Journal 67, 1085-1092, doi: 10.2136/sssaj2003.1085

van der Keur P, Iversen BV (2006) Uncertainty in soil physical data at river basin scale--a review. Hydrology and Earth System Sciences 10, 889-902. doi: 10.5194/hess-10-889-2006

van der Meer FD, de Jong SM, Bakker WH (2001) Basic physics of spectrometry. In 'Imaging spectrometry: basic principles and prospective applications. Remote sensing and digital image processing: Vol. 4'. (Eds FD van der Mecr, SM de Jong) pp. 12-13. (Springer: Berlin, Heidelberg)

van Genuchten MTh (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, 892-898. doi:10.2136/sssajl980.03615995004400050002x

van Genuchten MTh, Leij FJ, Yates SR (1991) The RETC code for quantifying the hydraulic functions of unsaturated soils, Version 1.0. EPA Report 600/2-91/065, U.S. Salinity Laboratory, USDA, ARS, Riverside, CA, USA.

Vereecken H, Maes J, Darius P (1989) Estimating the soil moisture retention characteristics from texture, bulk density and carbon content. Soil Science 148, 389-403. doi:10.1097/00010694198912000-00001

Webster R, Oliver MA (1992) Sample adequately to estimate variograms of soil properties. Journal of Soil Science 43, 177-192. doi: 10.1111/j.1365-2389.1992.tb00128.x

Wetterlind J, Stenberg B (2010) Near-infrared spectroscopy for within-field soil characterization: small local calibrations compared with national libraries spiked with local samples. European Journal of Soil Science 61, 823-843. doi:10.1111/j.1365-2389.2010.01283.x

Wosten JHM, van Genuchten MTh (1988) Using texture and other soil properties to predict the unsaturated soil hydraulic functions. Soil Science Society of America Journal 52, 1762-1770. doi: 10.2136/sssaj1988.03615995005200060045x

Wosten JHM, Lilly A, Nemes A, Le Bas C (1999) Development and use of a database of hydraulic properties of European soils. Geoderma 90, 169-185. doi: 10.1016/S0016-7061(98)00132-3

Wosten JHM, Pachepsky YA, Rawls WJ (2001) Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics. Journal of Hydrology 251, 123-150. doi: 10.1016/S0022-1694(01J00464-4

Wu C, Wu J, Luo Y, Zhang L, DeGloria SD (2009) Spatial prediction of soil organic matter content using cokriging with remotely sensed data. Soil Science Society of America Journal 73, 1202-1208. doi: 10.2136/sssaj2008.0045

Yates SR, Warrick AW (1987) Estimating soil water content using cokriging. Soil Science Society of America Journal 51, 23-30. doi: 10.2136/sssaj1987.03615995005100010005x

Ye M, Khaleel R, Schaap MG, Zhu J (2007) Simulation of field injection experiments in heterogeneous unsaturated media using cokriging and artificial neural network. Water Resources Research 43, W07413. doi: 10.1029/2006WR005030

Zhu Q, Lin HS (2010) Comparing ordinary kriging and regression kriging for soil properties in contrasting landscapes. Pedosphere 20, 594-606. doi: 10.1016/S1002-0160(10)60049-5

K. Liao (A), S. Xu (B,D), J. Wu (C), and Q. Zhu (A,D)

(A) Key Laboratory of Watershed Geographic Sciences, Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences, Nanjing 210008, China.

(B) Department of Environmental Science, Qingdao University, Qingdao 266071, China,

(C) Department of Hydrosciences, Nanjing University, Nanjing 210093, China.

(D) Corresponding authors. Email:;

Table 1. Statistical summary of basic soil properties
s.d., Standard deviation; CV, coefficient of variation

Variables            Minimum   Maximum   Mean     s.d.      CV

Bulk density          1.30      1.59     1.47     0.08     0.05
Organic matter (%)    0.68      2.65     1.36     0.38     0.28
Sand (%)              23.30     70.84    48.99    11.04    0.23
Silt (%)              12.02     50.69    33.88    8.37     0.25
Clay (%)              7.62      31.20    17.13    5.64     0.33

Table 2. Means and standard deviations of the fitted van
Genuchten parameters

Textural groups   Number     [[theta].sub.s]            [alpha]

Sandy loam          19     0.347 [+ or -] 0.061   0.030 [+ or -] 0.023
Silty loam          3      0.372 [+ or -] 0.020   0.014 [+ or -] 0.003
Loam                29     0.365 [+ or -] 0.056   0.038 [+ or -] 0.024
Sandy clay loam     3      0.350 [+ or -] 0.016   0.024 [+ or -] 0.010
Clay loam           4      0.376 [+ or -] 0.044   0.038 [+ or -] 0.016
All                 58     0.359 [+ or -] 0.054   0.034 [+ or -] 0.022

Textural groups               n                  [R.sup.2]

Sandy loam           1.302 [+ or -] 0.135    0.95 [+ or -] 0.03
Silty loam           1.270 [+ or -] 0.125    0.96 [+ or -] 0.03
Loam                 1.223 [+ or -] 0.072    0.95 [+ or -] 0.04
Sandy clay loam      1.315 [+ or -] 0.047    0.94 [+ or -] 0.05
Clay loam            1.128 [+ or -] 0.035    0.95 [+ or -] 0.05
All                  1.249 [+ or -] 0.108    0.95 [+ or -] 0.03

Table 3. Correlation matrix of the basic soil properties
via van Genuchten parameteras
* P < 0.05; ** P < 0.01

Variables           Sand        Silt        Clay

[[theta].sub.s]    -0.295 *     0.219       0.254
ln([alpha])        -0.061      -0.012       0.137
n                   0.489 **   -0.291 *    -0.524 **

Variables          Organic      Bulk
                    matter     density

[[theta].sub.s]     0.136      -0.653 **
ln([alpha])        -0.139      -0.306 *
n                  -0.468 **    0.192

Table 4. Pedotransfer functions for the van Genuchten parameters.
BD and OM denote bulk density and organic matter content,

Variables         Pedotransfer functions           [R.sup.2]

[[theta].sub.s]   1.034 - 0.46BD                     0.43
ln([alpha])       0.391 - 2.73BD                     0.10
n                 1.487 - 0.0075Clay - 0.08 lOM      0.34

Table 5. Correlation matrix of the basic soil properties via
digital numbers (DNs) of Landsat ETM imagery DN 1, DN 2, ...,
and DN 7 indicate the digital number of band 1, band 2, ...,
band 7 of Landsat ETM imagery. * P<0.05; ** P<0.01

                    DN 1        DN 2        DN 3

Sand               0.223       0.297 *     0.436 **
Silt              -0.173      -0.268 *    -0.398 **
Clay              -0.264 *    -0.366 **   -0.476 **
Organic matter    -0.475 **   -0.483 **   -0.435 **
Bulk density      -0.017      -0.032       0.102

                    DN 4        DN 5        DN 7

Sand              0.489 **     0.503 **    0.568 **
Silt             -0.412 **    -0.422 **   -0.453 **
Clay             -0.522 **    -0.545 **   -0.597 **
Organic matter   -0.316 **    -0.337 **   -0.276 *
Bulk density      0.104        0.065       0.082
COPYRIGHT 2014 CSIRO Publishing
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2014 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Liao, K.; Xu, S.; Wu, J.; Zhu, Q.
Publication:Soil Research
Article Type:Report
Geographic Code:9CHIN
Date:Aug 1, 2014
Previous Article:A combined equation to estimate the soil pore-water electrical conductivity: calibration with the WET and 5TE sensors.
Next Article:Using categorical soil structure information to improve soil water retention estimates of tropical delta soils.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters