Freely available online through the author-supported open access option.

Abstract

We revisited the Vereecken database, which has been used to derive pedotransfer functions (PTFs) to estimate the soil hydraulic parameters of Belgian soils. We developed new PTFs based on the Mualem–van Genuchten model, constraining m = 1 − 1/n and using fewer parameters. The goodness-of-fit was similar to the one originally obtained by Vereecken. We used a one-step procedure that allows direct quantification of the correlation matrix and the uncertainties of the estimated parameter values. The coefficients of the new PTFs were estimated using a global search algorithm and they were validated against independent data. The PTFs have a wider range of applicability since: (i) they allow the use of the closed-form solution of the unsaturated hydraulic conductivity in the Mualem–van Genuchten model; and (ii) they consider the effect of macroporosity. We determined that the hydraulic conductivity measured close to saturation could not be estimated based on the available estimators; however, the hydraulic conductivity in the matrix domain was predicted with high accuracy.

A new set of pedotransfer functions, based on a closed-form model, was developed from the Vereecken data set. The performance of these functions was as good as the existing published pedotransfer functions. Using the matrix saturated conductivity Ko instead of saturated hydraulic conductivity Ksat improved the shape of the predicted hydraulic conductivity curve.

Whatever the scale, soils are intrinsically heterogeneous and this heterogeneity controls their hydraulic behavior (Vogel and Roth, 2003). At the landscape scale, improving and developing methods to assess the hydraulic properties of the representative soil types in a quick and exhaustive manner is a crucial topic in soil science research. Besides the development of new experimental techniques to determine soil hydraulic properties, much effort has been devoted to building PTFs. This term, first introduced by Bouma (1989), covers mathematical functions predicting soil properties from more readily available data (e.g., texture, soil organic C, and bulk density). The main assumption underlying most PTFs is that textural properties dominate the hydraulic behavior of soils (Lilly and Lin, 2004).

There are some limitations to the PTFs presently available in the literature, however. First, most PTFs do not incorporate any structural information of the soil. It is well known that structural properties such as aggregation and macroporosity influence the soil hydraulic behavior. Structural information is especially crucial in the wet range of the moisture retention characteristic and hydraulic conductivity function. This is well illustrated by the low accuracy and the uncertainty of the predictions of saturated hydraulic conductivity, Ksat, from PTFs based on textural properties, bulk density, and soil organic matter. Nowadays, information about the soil structure is becoming increasingly available and is incorporated in predictive models (Pachepsky et al., 2006). Second, the reliability of PTFs used in a geographic region different from the one for which they were originally developed is often limited due to geologic, geohydrologic, climatic, and land use factors. Third, in the last two decades, simple nonlinear least squares optimization techniques have been the main tools to estimate hydraulic parameters and to derive PTFs. Current computer capabilities allow the use of global search algorithms, which can improve the estimation of hydraulic parameters and the development of PTFs.

Based on these considerations, we revisited the database of Vereecken, which contains data on Belgian soils. It was first used to derive PTFs for predicting soil hydraulic parameters of the van Genuchten model for the moisture retention characteristic and of the Gardner model for the hydraulic conductivity from texture, bulk density, and organic C content (Vereecken et al., 1989, 1990). These PTFs have been reported by several researchers to be among the most accurate they evaluated (e.g., Tietje and Tapkenhinrichs, 1993; Romano and Santini, 1997; Cornelis et al., 2001; Mermoud and Xu, 2006). They are not based on a closed-form model of the soil hydraulic properties, however, and this formulation is preferred when simulating transfer in the vadose zone.

The main objectives of this study were: (i) to explore the efficiency of different optimization methods, such as global search algorithms and a one-step method to derive PTFs, compared with the sequential estimation approach, performed with a local optimization algorithm used by Vereecken et al. (1989, 1990); (ii) to investigate the existence and value of a threshold pressure head that distinguishes macroporous from matrix flow in hydraulic conductivity; and (iii) to provide PTFs for the Vereecken data set based on a closed-form expression of the soil hydraulic functions.

Materials and Methods

Data Sets

The Vereecken data set was used to develop the new PTFs (Vereecken, 1988). The soil profiles were chosen to be representative of the main soil types of northern Belgium. Figure 1 shows the major textural classes of the sampled soil profiles. A total of 39 soil profiles were sampled and described according to the Belgian soil taxonomy. Soil cores for determining the moisture retention curve (MRC) and the hydraulic conductivity curve (HCC) were taken in the middle of each of the 182 horizons. Once the hydraulic measurements were done, the soil cores were used to determine the textural composition, organic C content, and bulk density. The granulometric composition was measured using dry sieving to separate the sand fractions and the pipette method to discern the finer particles. Nine fractions were thus determined that correspond to the Belgian classification scheme. The organic C content was determined by means of the Walkley and Black (1934) method. The bulk density was measured on the 100-cm3 soil cores used for the MRC, which were dried for 24 h at 105°C. The ranges across which the available properties varied are presented in Table 1 . Soil water content was determined at five matrix suctions below 100 cm using a sand-box apparatus and at higher suctions with low- and high-pressure chambers. For the HCC, two different methods were used: the crust method (Bouma et al., 1983) for the wet range and the hot air method (Arya et al., 1975) for the dry range. Part of the raw data could not be retrieved so that the final data set contained only 136 HCC and 166 MRC.

To evaluate the reliability of the new PTFs, we also took advantage of an independent data set containing moisture retention data collected in 48 horizons from 24 soil profiles sampled in the forests of Flanders (North Belgium). It was developed in the framework of assessing the predictive quality and usefulness of the Belgian soil map and historical forest soil profile data for soil mapping purposes (Cornelis et al., 2005). Each horizon was sampled twice using Kopecky rings of 100 cm3. Cylinders with stones, charcoal, or roots >2 mm in diameter were rejected and the soil was resampled in the same horizon. The particle-size distributions were determined on disturbed samples using a Coulter LS200 laser diffractometer (Beckman Coulter, Fullerton, CA). Organic matter content was determined using the Walkley and Black (1934) method. It ranged from 2.3 to 130.0 g kg−1. Bulk densities were measured by weighing the 100-cm3 undisturbed soil samples at 10 kPa and subtracting the corresponding mass of water measured on a 25-cm3 subsample. They varied from 0.76 to 1.78 Mg m−3.

For each soil ring, water content was measured at nine soil-matrix potentials to obtain the MRC. For the pressure potentials ranging from 1 to 10 kPa, the sand-box apparatus (Eijkelkamp Agri-Equipment, Giesbeek, the Netherlands) was used. For higher potentials, subsamples were used in pressure chambers (Soilmoisture Equipment, Santa Barbara, CA). Undisturbed 20-cm3 soil cores were used for measuring the water content at 20 and 33 kPa, whereas two disturbed subsamples served for water content determination at 1500 kPa.

Unfortunately, this database contained no hydraulic conductivity measurements and was therefore not suitable to evaluate the PTFs predicting the parameters of the conductivity model. As the extrapolation capacity of PTFs is generally very bad, the test set was reduced to contain only the measurements for which the basic properties fell within the range of the training set. Therefore only 39 soil horizons were available to assess the reliability of the MRC PTFs. The ranges of their properties are also presented in Table 1.

Texture fractions (sand, 50–2000 μm; silt, 2–50 μm; clay, <2 μm), bulk density, and organic C content were used as predictors for developing the new PTFs.

Models

Our major objective was to explore Vereecken's data set using new optimization methods and to assess the possibility of better estimating Ksat to improve the PTFs of Vereecken using fewer parameters. Several models were fitted on the data set. Table 2 displays their equations. For the retention curve only, several forms of the van Genuchten (1980) model were tested (shape parameter m fitted [VG0], m = 1 [VG1], and m = 1 − 1/n [VG2]), but also the models proposed by Kosugi (1996), Brooks and Corey (1964), and Ippisch et al. (2006), further referred to as KS, BC, and IP, respectively. The model referred to here as IP is actually the modified VG model proposed by Vogel et al. (2001), rewritten by Ippisch et al. (2006). Concerning the hydraulic conductivity, the model of Mualem (1976), referred to as Mu, was coupled with several of the models used for the MRC. In its general form, it can be written (van Genuchten, 1980) as 
\[K(S){=}K_{\mathrm{o}}S^{\mathrm{{\lambda}}}\left[\frac{{{\int}_{0}^{S}}\frac{1}{h(S)}\mathrm{d}S}{{{\int}_{0}^{1}}\frac{1}{h(S)}\mathrm{d}S}\right]^{2}\]
[1]
where K(S) is the unsaturated hydraulic conductivity [L T−1] as a function of the relative saturation S; Ko is a hydraulic conductivity acting as a matching point [L T−1] at suction head h = 0 [L], and λ accounts for the tortuosity and the connectivity of the porous network. The model of Gardner (1958), as used by Vereecken et al. (1990), was also tested. The equations are also displayed in Table 2.
We slightly modified the retention models by considering the moisture variation along the sample height, as suggested by Peters and Durner (2006). Rather than estimating the pressure head corresponding to the measured water content at half the sample height, its whole length was used. The measured water content was assumed to be the sample average water content, which varies. This can lead to more accurate results, especially near saturation where the water content may greatly vary along short pressure intervals for coarse or macroporous soils. The water content is now expressed as 
\[\mathrm{{\bar{{\theta}}}}{=}\frac{1}{L}{{\int}_{h_{\mathrm{l}}}^{h_{\mathrm{u}}}}\mathrm{{\hat{{\theta}}}}(\mathrm{B},h)\mathrm{d}h\]
[2]
where θ is the water content, the overbar represents the mean of the actual water content in the sample of total height L, hl and hu are the pressure heads at the bottom and the top of the sample, respectively, and B is the vector of parameters of the chosen model.
As MuVG in Table 2 was intended to simulate the hydraulic conductivity in the soil matrix domain, a distinction was made between Ko and Ksat. The former acts as a matching point at h = 0, without considering the macropores. The latter is the effective saturated conductivity. To take macroporosity into account, multimodal models (e.g., Priesack and Durner, 2006) could have been fitted on the measurements, but the number of parameters would have been too large with regard to the number of data. Hence, we opted for the introduction of a threshold value ht, from which macropores do not influence the HCC. Since too few measurements were available in the wet range, ht could not be estimated by optimization, but we fixed it at −6 cm, where our first unsaturated conductivity measurements occurred. This value is also in accord with the review of Jarvis (2007), who wrote that macropores can be defined as pores having an equivalent cylindrical diameter larger than 0.3 to 0.5 mm, which is equivalent to a water entry pressure head of −6 to −10 cm, according to the Young–Laplace equation: 
\[h{=}\frac{2\mathrm{{\gamma}{\,}cos}(\mathrm{{\theta}})}{\mathrm{{\rho}}gR}\]
where h is the pressure head, γ is the surface tension, θ is the contact angle between the two fluids, ρ is the fluid density, g is the gravitational acceleration, and R is the radius of the capillary.

Pedotransfer Functions

Method of Building Pedotransfer Functions

The classical sequential method to build PTFs consists in (i) fitting the hydraulic models on the measurements to derive the models' parameters, and (ii) building regression equations from the fitted parameters, using for example a stepwise method to choose the relevant variables entering the PTFs. Although widely used, this method has several drawbacks: (i) PTF coefficients are obtained from uncertain hydraulic parameters but uncertainty on these hydraulic parameters is usually not taken into account; (ii) the final uncertainty of the PTF coefficients is therefore hard to assess; and (iii) the forward model is largely overparameterized. We therefore opted for another method, proposed by Scheinost et al. (1997), which optimizes the PTFs' coefficients in one step. Empirical relations between each of the parameters of the hydraulic model and the potential predictors (i.e., variables entering the PTFs) must be set a priori. The coefficients of the PTFs are then directly fitted to the retention and conductivity observations through the hydraulic model (Eq. VG and MuVG in Table 2) using nonlinear regression. We further improved this method by allowing a stepwise optimization process to keep only significant regression coefficients. With this stepwise extended method, uncertainty in the final PTF coefficients can be fully characterized.

The main drawback of the method is that the equations of the relationships between the parameters and their predictors have to be set a priori. They must also be kept relatively simple to maintain the number of coefficients low enough for the optimization. Table 3 presents the equations chosen in this study for the MuVG model. For the retention parameters, we based our a priori relations on the equations retained by Vereecken for his published PTFs (Vereecken et al., 1989). Regarding the Mualem–van Genuchten model, we entered into the equations all our available predictors without any transformation. The coefficients of these equations were optimized through a stepwise procedure that minimizes the objective function presented below. The stepwise procedure drops the less relevant predictors. In case some coefficients are not significantly different from zero (P > 0.05), the optimization is performed again without the coefficient with the highest P value, and so on until all the coefficients included in the model have a P value <0.05. The results from a stepwise procedure depend on the chosen path to remove the nonsignificant coefficients. Even if it does not necessarily lead to the best performance, removing the coefficient with the largest P value is a common practice and seems the most objective choice.

As parameter Ksat is not included in this set of PTFs, another regression equation was developed to predict it, using a stepwise regression procedure.

Objective Functions

The objective functions for the MRC and HCC are classical sums of squares, based on the direct models presented in Table 2. Note that for HCC, observations for |h| < |ht| were not included in the optimization process, since KoKsat for soils with macropores. The value of ht was chosen as 6 cm of suction, which falls in the range of 6 to 10 cm proposed by Jarvis (2007) and beneath which we did not have other measurement data than the saturated conductivity itself.

For the coupled optimization, the objective functions of both the MRC and the HCC were combined to result in 
\begin{eqnarray*}&&{\Phi}(\mathrm{B}){=}{{\sum}_{i{=}1}^{n_{1}}}v_{i}\left[\mathrm{{\bar{{\theta}}}}_{i}{-}\frac{1}{L}{{\int}_{h_{\mathrm{l}_{i}}}^{h_{\mathrm{u}_{i}}}}\mathrm{{\hat{{\theta}}}}(\mathrm{B}_{\mathrm{MRC}},h_{i})\mathrm{d}h\right]^{2}\\&&{+}{{\sum}_{j{=}1}^{n_{2}}}w_{j}\left\{\mathrm{ln}(K_{j}){-}\mathrm{ln}\left[{\hat{K}}(\mathrm{B}_{\mathrm{HCC}},h_{j})\right]\right\}^{2}\end{eqnarray*}
[3]
where vi and wj are used to balance the two parts of the function. A different method was applied to compute the weights for the choice of the model to be applied to the data set and for building the PTFs.

In the first case, the weights that were used are the inverse of the variance of the model errors when the MRC or HCC models were directly fitted on the data.

In the second case, the PTFs were first fitted once separately for the MRC and for the HCC parameters, and the variances of the errors for the whole data set for the MRC, σθ2, and the HCC, σK2, were computed. The weights were also computed as the inverses of these variances. The same weights were hence applied for all the MRC data on the one hand and for all the HCC data on the other hand. The weights have here no statistical significance but are just a way to balance the two parts of the objective function. Applying the same weight to all the data of one kind would prevent some data that would have been badly fitted in the first optimization to be condemned to never improve.

Optimization Algorithm

The optimizations were performed using a global multilevel coordinate search (GMCS) combined with a local Nelder–Mead Simplex (NMS) algorithm, as proposed by Lambot et al. (2002). The GMCS is an algorithm written in Matlab (The MathWorks, Natick, MA) for bound constrained global optimization using function values only, based on a multilevel coordinate search that balances global and local searches. The local search is done via sequential quadratic programming (Huyer and Neumaier, 1999). In the approach used in this study, a global search was first performed by omitting the local search algorithm to avoid unnecessary function evaluations. Then a local search with NMS was realized on the parameter set provided as output from the GMCS calculations.

Accuracy and Reliability

The accuracy and reliability of the newly developed PTFs were evaluated using the criteria presented in Table 4 .

Uncertainty

The uncertainties were evaluated during the PTF calibration process. Assuming that measurement errors are normally distributed and that the sum of squares is locally linear in the parameters, we approximated the variance–covariance matrix of the coefficients C as suggested by Kool and Parker (1988): 
\[\mathbf{\mathrm{C}}{\sim}\frac{\mathbf{\mathrm{e{^\prime}}}_{\mathrm{w}}\mathbf{\mathrm{e}}_{\mathrm{w}}[\mathbf{\mathrm{J{^\prime}}}_{\mathrm{w}}\mathbf{\mathrm{J}}_{\mathrm{w}}]^{{-}1}}{N{-}p}\]
[4]
where J is the Jacobian matrix, i.e., the jth column of J is the partial derivative of the simulated data with respect to the jth coefficient; e is the vector of differences between the measurements and the simulated data; and index w means that the matrix has been multiplied by the Cholesky decomposition of the weight matrix. This is the diagonal matrix formed by the elements of the vector of the weights. The covariance matrix obtained in this way allowed us to evaluate confidence intervals for the regression coefficients as well as their correlation matrix. The statistic t =
\(\mathrm{{\hat{{\beta}}}}_{i}/{\surd}(\mathbf{\mathrm{C}}_{ii}\)
) and follows a Student statistical distribution with Np degrees of freedom, and its associated P value gives an indication of the statistical significance of the estimation of coefficient βi.

Results and Discussion

Water Retention and Hydraulic Conductivity Models

To develop the new PTFs on the basis of a relevant model for the description of the soil hydraulic functions, we fitted several models on the training data set, presented in Table 2. The parameters of the models were optimized using the global algorithm GMCS. Their performances are presented in Table 5 .

Comparing the performances of the van Genuchten model with parameter m = 1 (VG1) and 1 − 1/n (MuVG) shows that VG1 is better. Comparing the RMSE distributions of the two models showed that their mean RMSEs are significantly different (P = 0.0171). This corroborates the findings of Vereecken et al. (1989). Figure 2 presents an example of the difference in the shapes of the retention curve. The VG1 model fits the decrease in water content slightly better. The predictions of the two models mainly differ in the wet and dry range. Besides, near saturation, the curves do not behave the same way. With m = 1, n can be <1, which leads to a decreasing slope rather than to the sill observed with m = 1 − 1/n and n > 1. The retention model of Ippisch et al. (2006), here referred to as IP, is a rewriting of the modified VG introduced by Vogel et al. (2001). It introduces an air-entry point in the VG model that does not lead to great changes in retention but can improve the behavior of the simulated HCC near saturation, mostly when n < 2. When the air-entry point, he, is 0, the model reduces to MuVG; however, it has an extra parameter that must be optimized. Concerning the combination of the MRC and the HCC, Table 5 shows that, if VG1 proved to better simulate the MRC than VG2, the Gardner model performed worse than MuVG. The MuKS model performed the worst. The MuIP model performed the best, but it has an extra parameter compared with MuVG and its RMSE for the MRC is larger than that of MuVG. Besides, MuVG is already implemented in many models. Hence, we chose the classic MuVG for the development of the new PTFs.

Pedotransfer Function Coefficients

Table 6 presents the predictors and their coefficients retained for each MuVG parameter in the new PTFs. Parameter θr does not appear in the table; it is not significantly different from 0, and was therefore excluded from the new PTFs. This has already been observed in various PTFs (e.g., Wosten et al., 1999, Zacharias and Wessolek, 2007). The saturated water content, θs, was predicted with clay and bulk density. As bulk density (BD) is directly related to the soil porosity F by the relation F = 1 − BD/ρs (where ρs is the density of the solid particles of the soil), it is expected that the coefficient multiplying bulk density will be negative (Scheinost et al., 1997) as soils with large porosity tend to retain more water near saturation. The predictors for parameter α* are sand, clay, and organic C contents. The presence of clay and sand variables as predictors reflects the dependency of α on the geometric mean particle diameter (Scheinost et al., 1997), as the two predictors together fully characterize the soil texture. Parameter n*, which is linked to the shape of the particle-size distribution, is predicted using clay and sand contents as well as the square of the latter. Parameter K* depends on sand content, bulk density, and organic C content. Parameter λ, which reflects the tortuosity and the connectivity of the pores, depends on the clay and sand contents.

Pedotransfer Function Performance

From the 25 coefficients proposed in Table 3, the stepwise procedure retains only 18 (see Table 6).

We compared the new PTFs (referred to as WEY) with the old ones (VER) and with two other published PTFs: Hypres (Wosten et al., 1999) and Rosetta (Schaap et al., 2001). Hypres (HYP) was developed using multiple regressions on a European database, whereas Rosetta (ROS) is based on artificial neural networks built from soils in temperate to subtropical climates in North America and Europe.

Moisture Retention Curve

To evaluate the accuracy and the reliability of the PTFs (see Table 7 ) for the MRC, we used a training set and a test set. Using the training set, the performance of the new PTFs (WEY) was almost as good as the old one (VER). Figure 3a shows that the spread of the simulated values plotted against the observations is quite narrow, except for a few points that were overestimated. The mean relative error (MRE) is 40% larger than the value observed for the old PTFs (VER). This means that the errors made on small water content values were larger with the new PTFs (WEY). Nevertheless, the means of the relative errors averaged on each soil horizon are not significantly different (P = 0.0522). The results produced by HYP and ROS are comparable to those obtained with VER and WEY, although a little worse, with ROS being worse than HYP. This illustrates the very specific character of hydraulic properties as related to soil-climate conditions. The HYP PTFs developed from European soils performed better than ROS, which encompasses European and U.S. soils. But HYP, as a pan-European set of functions, was less accurate than VER or WEY, developed for Belgian soils.

Testing the performances of the PTFs with an independent data set, we observed larger deviations compared with the training set. Figure 3b illustrates for WEY and VER that the scattering around the straight 1:1 line is much broader than for the training set. Moreover, the high water content values were systematically underestimated. There are three main explanations for this. First, forest soils are typically highly structured and this is not directly included in the PTFs. Information about the soil structure (such as the size and shape of the soil aggregates or the pedological development of the profile) could be integrated into PTFs, but it is mostly stored as qualitative and subjective data. A second reason may be the difference in the textural characterizations used for both data sets. Textural properties in the training set were determined using the classical combination of sieving and decantation, whereas a laser method was used for the test set. A third reason may be the overparameterization of the model. The comparison between the distribution of the relative errors of VER and WEY PTFs shows that their means are highly significantly different (P ≤ 0.001). It is also worth noting that the performances of HYP and ROS on the test set were also far worse than on the training set. This pleads for the first reason proposed above since none of these PTFs really consider the soil structure.

Hydraulic Conductivity Curve

The performances of the PTFs in predicting the hydraulic conductivity are presented in Table 8 . Most of the indicators show better performance of the new PTFs (WEY). Only the MRE is large compared with the one observed with the old PTFs (VER). Figure 4 displays the scattering of the data predicted with WEY and VER.

The results obtained with HYP and ROS are also presented in Table 8. Relative to the results of the MRC alone, their performances are worse. Depending on the indicator, ROS performed better than HYP, or vice versa. The results obtained for the MRC by the large-scale PTFs were comparable to the local PTFs. This is not true for the HCC, where the need to take local specifications into account is more crucial.

Unfortunately, no test data set containing independent hydraulic conductivity measurements was available in this study to evaluate the reliability of these PTFs. Note that the measured saturated conductivity was not accounted for in these PTFs.

Table 9 presents the coupled performance indicators obtained with the different PTFs on the training set. The MRC and the HCC indicators were balanced using the inverse of the variances of the measurements. The WEY performed the best except for the mean absolute error, but only slightly better than VER. The ROS and HYP PTFs were both worse.

Matching-Point vs. Saturated Hydraulic Conductivities

The Mualem–van Genuchten model (MuVG) only describes hydraulic conductivity in the matrix flow domain. Therefore we did not include hydraulic conductivity measurements at and near saturation (|h| ≤ |ht|) in the optimization process. Predictions of the estimated Ksat using the same variables as in the PTFs did not produce satisfactory results, even when using various transformations of the variables (maximum obtained R2 = 0.25).

Figure 5 shows the three first components of the principal component analysis realized on ln(Ksat) and the predictors used for the PTFs (sand content, clay content, organic C content, and bulk density). These three components explain 83.58% of the total variance. The ln(Ksat) mainly contributes to the third component, whereas the different predictors are more represented by the first and second components. This means that these predictors are not suitable to predict the saturated conductivity. Two reasons could explain why Ksat cannot be predicted by basic soil properties. First, it is well known that it is mainly affected by macropores that are not directly linked to the textural properties of the soils. Therefore other variables that are more related to soil structure (e.g., profile development, horizon, land use, the presence of roots) could be added for developing new PTFs for Ksat. Second, the measurement support or the number of samples was not adequate to characterize Ksat. Conductivity measurements are indeed highly scale dependent and sensitive to the measurement method, especially close to saturation (Javaux and Vanclooster, 2006; Mallants et al., 1998). As PTFs are usually used for large-scale studies, excluding any nonsignificant very localized effect (namely the location of one or several macropores) is justified. There is a need to develop new methods to characterize effective large-scale Ksat based on other indicators, however, such as soil surface structural properties, agricultural practices, or crops.

The distribution of the ratio ln(Ksat)/ln(Ko) by soil type can be used to assess the presence of soil macroporosity. Figure 6 presents a box plot of the ratio for each textural class. The dispersion of the ratio tends to increase for soils having a decreasing content of sand (sand < silt loam < silt).

The number of observations available for the other classes is too low to draw conclusions with respect to this. Plotting sand content against the ratio, however, as shown in Fig. 7 , supports the finding that the ratio ln(Ksat)/ln(Ko) shows a larger spread for soils with a low content of sand. Nothing can be concluded about the clay content, for which the sampling does not encompass the whole textural range. This suggests that fine-textured soils are more prone to macroporous flow; however, measurements in fine-textured soils are more prone to errors from side effects.

Pedotransfer Function Uncertainties and Correlations

The one-step inversion method allows a direct quantification of the PTF uncertainties and the correlation matrix of their coefficients. Contrary to the sequential method, the final uncertainty encompasses at the same time data set experimental errors, model errors, and uncertainty due to fitting PTFs to the data set. The values of the relevant PTF coefficients, their confidence intervals, and their significance are presented in Table 6. These uncertainties, along with the correlation matrix shown in Table 10 , can be used to perform uncertainty analysis for specific flow problems.

The correlation matrix of the hydraulic parameters predicted with the PTFs (see Table 11 ) reveals that α* and K* on the one hand and n* and λ on the other hand are very strongly correlated.

The regression between parameters ln(Ko) and ln(α), obtained from fitted HCC and MRC functions (see above) produced the relation ln(Ko) = 8.94 + 1.78 ln(α) with an R2 of 0.59. The value of 1/α corresponds to the inflection point of the MRC described by the VG model. We can relate it to Ko using a set of physics laws. First, Poiseuille's law 
\[\mathrm{{\phi}}{=}\frac{\mathrm{d}V}{\mathrm{d}t}{=}v\mathrm{{\pi}}R^{2}{=}\frac{\mathrm{{\pi}}R^{4}}{8\mathrm{{\eta}}}\left(\frac{{-}{\Delta}P}{{\Delta}x}\right){=}\frac{\mathrm{{\pi}}R^{4}}{8\mathrm{{\eta}}}\frac{{\Delta}P}{L}\]
explicitly relates the volumetric laminar stationary flow of an incompressible uniform viscous liquid through a cylindrical tube with constant circular cross-section (ϕ) to the internal radius of the tube (R) at power 4 (V is a volume of the liquid poured [m3], t is time [s], v is the mean fluid velocity along the length of the tube [m s−1], x is a distance in the direction of flow [m], ΔP is the pressure difference between the two ends [Pa], η is the dynamic fluid viscosity [Pa s], and L is the total length of the tube in the x direction [m]). By the Young–Laplace equation, 1/α is inversely proportional to the average pore radius. Combining these two laws with Darcy's law, 
\[q{=}\frac{\mathrm{{\phi}}}{A}{=}K\frac{{\Delta}H}{L}\]
where q is the flux density, φ is the discharge, A is the unit cross-sectional area, K is the proportionality factor called hydraulic conductivity, and ΔH/L is the hydraulic gradient, we can conclude that ln(Ko) is linearly linked with ln(α), with a slope of 2, which is close to the value of 1.78 obtained in the regression above. The intercept should be equal to ln[(ρg/2γcosθ)2(|ΔP|/8ηL)], which is 12.35 in the units used here. Again, this value is not far from 8.94. Note that this strong correlation does not exist between α and Ksat, as experimental Ksat is mainly influenced by the macropore network (R2 = 0.0101).

Both λ and n influence the general slope of the HCC and MRC. It is therefore expected that the same type of information will affect these parameters and that they therefore will be correlated. These strong correlations suggest that the parameters specific to the HCC in the MuVG model could be predicted using the optimized parameters of the MRC modeled with VG2. This interesting suggestion deserves further investigation. Rosetta (Schaap et al., 2001) has already proposed to predict Ko and λ from the VG model parameters optimized on measurement data on the basis of the work of Schaap and Leij (2000).

Conclusions

New PTFs were developed for Belgium based on the Vereecken data set. The new PTFs, developed using the coupled model of Mualem–van Genuchten (van Genuchten, 1980), performed as good as the ones published by Vereecken et al. (1989, 1990), based on the models of van Genuchten (1980) with m = 1 − 1/n and Gardner (1958), and they used fewer coefficients (18 vs. 26), fewer parameters (six vs. seven) and were more accurate for K(h). The one-step inversion method coupled with a stepwise optimization procedure led to a consistent parameter set for the PTFs. In addition, it allowed the direct quantification of the correlation matrix and the uncertainties on the PTFs. The strong correlations observed between some parameters can be explained by physical reasoning. The possibility of predicting the HCC parameters from the MRC parameters was suggested but needs to be further investigated. Besides, we showed that close to saturation, conductivity measurements cannot be predicted based on any of the available variables. Instead, the new PTFs predict the saturated hydraulic conductivity for the matrix Ko with high accuracy. The resulting PTFs therefore predict the hydraulic conductivity in the soil matrix and not in the macropores. There is therefore a need for developing alternative methods to asses relevant Ksat values based on other predictors, experiments, or support scales. Improved Ksat could then be linked to the hydraulic conductivity at the ht value. The function between h = 0 and h = ht could be chosen based on experimental data sets. We advise that more systematic measurements should be performed in the range 0 to 15 cm to further investigate the optimal threshold value ht.

Further research could focus on the improvement of the prediction if any information about the soil structure was introduced into the model. As shown here, predictions using texture, bulk density, and organic C content failed to predict the hydraulic characteristics of highly structured soils or the hydraulic conductivity near saturation. Structure-related information could consist of pedologic data, gathered during soil surveys (Lin et al., 2006a,b; Pachepsky et al., 2006). Information about land use and crop management, which strongly affect the soil surface, could also be taken into account (Jarvis, 2007).

The research of M. Weynants is funded by the Fonds pour la Formation à la Recherche dans l'Industrie et dans l'Agriculture (F.R.I.A.). The test set used in this study was set up in the framework of Project 00/05 of the VLINA research program (1995–2001), which was funded by the Ministry of the Flemish Community and conducted by Ghent University and the Institute for Forestry and Game Management.