An alternative deterministic method for the spatial interpolation of water retention parameters

There are two approaches available for mapping water retention parameters over the study area using a spatial interpolation method. (1) Retention models can be first fitted to retention curves available at sampling locations prior to interpolating model parameters over the study area (the FI approach). (2) Retention data points can first be interpolated over the study area before retention model parameters are fitted (the IF approach). The current study compares the performance of these two approaches in representing the spatial distribution of water retention curves. Standard geostatistical interpolation methods, i.e., ordinary kriging and indicator kriging, were used. The data used in this study were obtained from the Las Cruces trench site database, which contains water retention data for 448 soil samples. Three standard water retention models, i.e., Brooks and Corey (BC), van Genuchten (VG), and Kosugi (KSG), were considered. For each model, standard validation procedures, i.e., leave-oneout cross-validation and split-sample methods were used to estimate the uncertainty of the parameters at each sampling location, allowing for the computation of prediction errors (mean absolute error and mean error). The results show that the IF approach significantly lowered mean absolute errors for the VG model, while also reducing them moderately for the KSG and BC models. In addition, the IF approach resulted in less bias than the FI approach, except when the BC model was used in the split-sample approach. Overall, IF outperforms FI for all three retention models in describing the spatial distribution of retention parameters. Correspondence to: H. Saito (hiros@cc.tuat.ac.jp)


Introduction
The predictions of soil moisture distributions in the vadose zone or estimates of contaminant arrival time to groundwater rely heavily on robust estimates of the spatial distribution of soil hydraulic parameters.Spatial interpolation techniques have been used to estimate the unknown values of these parameters at unsampled locations from available observations.It has been widely accepted that soil hydraulic parameters are spatially correlated to a greater or lesser extent.Because of this, techniques that take such information into account must be used for spatial interpolation.Among many available techniques, including the inverse distance method or the linear interpolation method, only a least-square interpolation technique called kriging accounts for spatial correlations between variables.Kriging is now commonly used for mapping soil physical, chemical, and/or hydraulic properties.Kriging estimates not only the values of an attribute at unsampled locations, but also their uncertainties in terms of an error variance known as kriging variance, when the underlying geostatistical model is correct (Goovaerts, 1997).
Unlike water flow in saturated systems, predictions of variably-saturated water flow in soils depend not only on knowledge of saturated hydraulic conductivities, but also on knowledge of the water retention and unsaturated hydraulic conductivity characteristics that are usually described by various functional forms.Although experimental data for both retention and hydraulic conductivity functions are required for efficient parameterization, in many cases, only retention data is available.Retention curves are relatively easy to collect when in comparison with unsaturated hydraulic conductivities, which are generally difficult and time-consuming to Published by Copernicus Publications on behalf of the European Geosciences Union.measure.Typically, various pore size distribution models (e.g., Mualem, 1976;Burdine, 1953) are then used to predict unsaturated hydraulic conductivity functions from the retention data.In this study, only water retention curve data will be analyzed, since we have assumed that only this information is available, while information about unsaturated hydraulic conductivities is lacking.The water retention parameters for a particular soil are usually obtained by first measuring a series of pressure head and water content data pairs from core samples in a laboratory and then by fitting the constructed discrete curve using a simple analytical model, such as the well-known van Genuchten model (van Genuchten, 1980).Water flow in soils is then usually obtained using a numerical model that simulates variably-saturated water flow, and uses analytical models of soil hydraulic properties.Predictions of water flow are thus greatly improved when estimated parameters for a given soil hydraulic model adequately represent the highly nonlinear relationships between water contents, pressure heads and hydraulic conductivities.Although accurate estimation of water retention curves (WRC) or their model parameters is rarely a goal, it is one of the most crucial steps in modeling water flow and solute transport in the vadose zone.
It is well known that soil water retention parameters are not only soil-type dependent, but also spatial-location dependent.There are two possible approaches to map water retention parameters over the study area using a spatial interpolation method.(1) Retention models are first fitted to experimental retention curves available at sampling locations, and then the obtained model parameters are interpolated over the study area.(2) Retention data can first be interpolated over the study area, and then the retention model is fitted to this interpolated data.The former approach is called the FI (f itfirst and interpolate-later) approach, and the latter is referred to as the IF (interpolate-first and f it-later) approach (Fig. 1).When non-linear processes are involved in either step, results differ depending upon the approach chosen (Heuvelink and Pebesma, 1999;Leterme et al., 2007).
The FI approach treats the estimated water retention parameters as meaningful physical parameters.However, since some water retention parameters are derived in a purely em-pirical way and do not have a physical meaning, the evaluation of these parameters has to be done with great care, or the results may be misleading and meaningless (Romano, 2004).However, there are still a number of studies that rely solely on water retention parameters, because one can dramatically reduce the number of parameters that need to be analyzed.For example, Zhu and Mohanty (2002) averaged the widely used van Genuchten "parameters" to simulate large-scale infiltration and evaporation processes.In Oliveira et al. (2006), stochastic fields of "parameters" were generated using conditioning parameters to simulate variably-saturated water flow in soils.In some studies, the number of parameters is even reduced to one or two by using the scaling method (e.g., Shouse et al., 1995;Oliveira et al., 2006).The question then arises whether or not we can consider these derived (i.e., averaged or randomly generated) parameters to be equivalent to experimentally measured soil physical properties, such as volumetric water contents.In other words, to map water retention parameters, which approach, FI or IF, should we take?
Although there are a great number of studies that have analyzed and/or modeled the spatial distribution of soil hydraulic properties, including water retention model parameters (e.g., Wendroth et al., 2006;Unlu et al., 1990), we do not know of any studies that would rigorously investigate this question.Tomasella et al. (2003) compared two approaches to map retention parameters, in a manner similar to our study.In one approach, called the point-based approach, water retention data (the relationship between water contents and water potentials) were predicted using pedotransfer functions (PTF) prior to fitting retention models, while in the other approach, called the parametric approach, retention model parameters were directly predicted using PTFs.By using a comprehensive soil water retention database of Brazilian soils, they found that the point-based approach outperformed the parametric approach.Sinowski et al. (1997) also compared two approaches to map water retention parameters using PTFs.In one approach, soil properties that were used as input parameters of PTFs were interpolated before applying PTFs.In the other approach, water retention parameters obtained through PTFs were interpolated directly.Their comparison showed that the former approach outperformed the latter approach.However, the use of PTFs in predicting retention data introduces additional uncertainty (Schaap et al., 2001) that may be larger than uncertainties due to the use of fitting and interpolation techniques.
This study therefore aims to quantify how well the estimated soil hydraulic model parameters reproduce the actual water retention curve at a given location using one of two approaches discussed above (i.e., IF and FI approaches).The two approaches are evaluated in terms of prediction errors using standard validation techniques.

Water retention data and models
The soil hydraulic data set (i.e., water retention curves) used in this study is available through the Las Cruces Trench Site database (Wierenga et al., 1989).The database was generated as part of a comprehensive field study conducted in southern New Mexico, near Las Cruces, for validating and testing numerical models of water flow and solute transport in the unsaturated zone (Wierenga et al., 1991).A 24.6-m long by 6.0-m deep trench wall was excavated, and a total of 450 undisturbed samples, 7.6 cm internal diameter and 7.6 cm long, were taken from nine layers, 50 samples were taken from each layer, for every 0.5 m.Soil samples (undisturbed and disturbed) were then analyzed in the laboratory for soil properties, such as bulk density, saturated hydraulic conductivity, and soil water retention curve.Soil water retention curves were determined at 448 sampling locations, u i (i=1, 2. . .448) (Fig. 2).At each location, the water contents, θ u i ; h j , were measured at eleven pressure heads, h j (j=1, 2. . .11) of 0, −10, −20, −40, −80, −120, −200, −300, −1000, −5000, and −15 000 cm H 2 O.While undisturbed soil cores were used for the wet range (−300 cm to 0 cm), disturbed soil samples were used with a standard pressure plate apparatus for the dry range (−15 000 cm to −1000 cm).More details regarding the experimental procedures can be found in Wierenga et al. (1989).The soil characterization and infiltration experiments at this site have been analyzed in a number of studies, including recent articles by Rockhold et al. (1996), Oliveira et al. (2006), andTwarakavi et al. (2008).
Figure 2 shows the location map of measured saturated water contents when the pressure head is equal to 0 cm or h 1 at the site.The saturated water content data are, as expected, spatially heterogeneous with mean, maximum, and minimum values of 0.322, 0.529, and 0.218, respectively.Although not shown here, water contents measured for other pressure heads are also spatially heterogeneous.Figure 3 shows sample histograms of water contents for eleven pressure heads.The mean of each distribution decreases as the pressure head decreases, as expected.Not all histograms show highly skewed or asymmetric distributions, indicating that no data transformation is required.Figure 4

depicts the experimental
water retention curves obtained for seven different depths (z) at x=10.75 m.There are some discontinuities between the water contents observed for pressure heads of −300 cm, and those observed for −1000 cm, where water contents increase as pressure heads decrease.This is likely due to the different measurement procedures used; undisturbed soil cores were used for water content measurements from −10 to −300 cm pressure heads, while disturbed soil cores were used for pressure heads ranging from −1000 to −15 000 cm (Hills et al., 1993).Such inconsistency should have an impact on the results of model fitting, as all data are equally weighted.Investigating the impacts of this inconsistency on final parameter estimates by FI and IF approaches requires rigorous analysis of how errors are propagated between different approaches.Not surprisingly, all seven curves are quite different.Although this variation in the vertical direction is expected to be much more pronounced than in the horizontal direction, variations in the horizontal direction also exist and cannot be ignored, as can been seen from Fig. 2.This confirms the importance of analyzing the spatial distribution of water retention curves at this site.Among many available retention curve models, the Brooks and Corey (BC) model (Brooks and Corey, 1964), the van Genuchten (VG) model (van Genuchten, 1980), and the Kosugi (KSG) model based on the lognormal pore size distribution (Kosugi, 1996) are the most widely used analytical expressions for representing the dependence of the water content on the capillary pressure head for unimodal pore systems.Detailed discussions of the commonly used soil hydraulic models can be found in Leij et al. (1997).Analytical expressions of the three closed-form models used in this study are listed in Table 1.The parameter S e is an effective saturation defined as follows: where θ s and θ r are the saturated and residual water contents (m 3 m −3 ), respectively.As can be seen from  empirical) parameters, leading to a total of four parameters representing the retention curve models.
Application of these water retention functions requires reliable estimates of the parameters for a soil of interest.While θ s can be obtained independently, the other retention parameters usually have to be estimated indirectly, by fitting the analytical function to the experimental water retention data using an optimization approach, such as a non-linear least-squares minimization approach (e.g., implemented in the RETC code).In most optimization procedures, performance is improved if the initial estimates are close to the "true" solution.In other words, if initial estimates are too far from the "true" solutions, the solution may not converge to the global minimum, but to the local minimum during the optimization process.Seki ( 2007) recently developed a program code that uses a full-automatic procedure to estimate soil water retention model parameters.The program automatically selects initial estimates based on observations, so that the user does not have to make his/her selections.De- 1.00 2.00 3.00 .000 .040 .080 .120  tails of this approach will not be discussed in this paper, but it should be noted that the full-automatic approach worked very well in most cases when applied to a number of water retention data from the UNSODA database (Leij et al., 1996).Such an approach is especially attractive when a great number of retention curves need to be fitted simultaneously.All three models were fitted to 448 retention data points using the full-automatic procedure to obtain water retention parameters at all sampling locations.Figure 5 shows histograms of retention parameters, including θ s and θ r .Similar to the water content histograms shown in Fig. 3, θ s and θ r histograms show no significant skewness and are more-or-less symmetric, except for θ r of the BC and VG models, where about 12 and 2.5% of data have a value of 0, respectively.As for the shape parameters (e.g., α v and n v ), most are positively skewed.The parameter h b even shows a bimodal distribution.

Geostatistical interpolation theory
To estimate an unknown value of a given soil property at unsampled locations, a geostatistical least-square interpolation technique, known as kriging, was used in this study.Consider the problem of estimating the value of a soil attribute z (e.g., water content or soil hydraulic parameters) at an unsampled location u, where u is a vector of spatial coordinates.The available information consists of values of the variable z at N locations u i , i=1, 2. . .N. All univariate kriging estimates are variants of the general regression estimate z * (u) defined as: where λ α (u) is the weight assigned to datum z(u α ) and m(u) is the trend component of the spatially varying attribute (Goovaerts, 1997) given neighborhood or window W (u) centered on u, while the influence of those farther away are discarded (Saito and Goovaerts, 2000).One of the most common kriging estimators is ordinary kriging (or OK), which estimates an unknown value as a linear combination of neighboring observations: In OK, unlike the constant value used in simple kriging, the mean (or trend) at each estimation location (i.e., local mean, m(u)) is implicitly re-estimated.OK weightsλ OK α are determined so as to minimize the error or estimation variance σ 2 (u)= Var{Z * (u) − Z(u)} under the constraint of an unbiased estimate (Goovaerts, 1997).These weights are obtained by solving the system of linear equations known as the ordinary kriging system: where µ OK (u) is the Lagrange parameter, γ (u α -u β ) is the semivariogram between observations at u α and u β , and γ (u α -u) is the semivariogram between the datum location u α and the location being estimated, u.The semivariogram γ (h) models the variability between observations separated by a vector h.The only information required by the system (Eq.4) are the semivariogram values, γ , for different separation distances.These values are readily derived from the semivariogram model fitted to experimental values (i.e., linear model of regionalization): (5) where N(h) is the number of data pairs for a given separation vector h.The choice of the model is limited to functions that ensure a positive definite covariance function matrix of the left-hand-side of the kriging system (Eq.4).Spatial correlations often vary with direction, and such a case requires one to compute semivariograms in different orientations to fit anisotropic (direction-dependent) semivariogram models.
In this study, differences in horizontal and vertical variations observed in retention data were taken into account by considering anisotropic semivariograms.Details of model fitting can be found in Deutsch and Journel (1998), Goovaerts (1997), andKitanidis (1997).In addition to an estimate of the unknown z value, OK provides an error variance that is computed as When histograms show positively skewed distributions and/or are not unimodal, OK might not be the best approach to estimate values at unsampled locations (Saito and Goovaerts, 2000).In such case, indicator kriging (IK), in which all data are transformed into either 0 or 1 depending upon the exceedence of a given threshold (Journel, 1983), can be used.Indicator kriging allows one to estimate the probabilities of exceeding a series of thresholds at unsampled locations, leading to the construction of conditional cumulative distribution functions (ccdf) of the variable of interest.Once the ccdf is obtained, an optimal estimate and the uncertainty associated with the estimated value can be quantified by taking the mean (E-type estimate) and variance of the ccdf, respectively.
The conditional probabilities at u for a series of thresholds z k discretizing a range of variation of z are estimated as a linear combination of n(u) surrounding indicator transformation of data i(u α ; z k ): where Kriging weights λ α (u; z k ) are then computed from the indicator kriging system, similar to the OK system (Eq.4), where all experimental semivariograms γ (u) are replaced by indicator semivariograms γ I (u).A major drawback of IK is that in order to construct ccdfs at any location u, K indicator semivariograms, usually 9 to 16 (Deutsch and Journel, 1998), have to be computed, and K indicator kriging systems have to be solved.Consequently, IK is much more computationally demanding than OK.It is common to use nine deciles of the sample distribution as thresholds.More details on indicator kriging can be found in Goovaerts (1999).

Prediction performances: FI vs. IF
To compare the FI and IF approaches in terms of reproduction of WRC at i locations u i , S e (h; u i ), two standard validation techniques were used: a leave-one-out cross-validation, in which one observation at a time is temporarily removed from the dataset and re-estimated from remaining data, and a split-sample method, where the dataset is divided into two non-overlapping observation and validation sets (Fig. 6).The number of observation locations in each set were 199 (N p ) and 249 (N v ), respectively.In the split-sample method, the observation set was used to estimate values at all locations in the validation set so that estimation errors could be calculated explicitly.Figure 6 depicts the location map of both observation and validation sets used in the split-sample method.
For complete full validation, different split scenarios had to be tried.In this study, only the scenario shown in Fig. 6 was used.This scenario was considered because it is the most likely and realistic scenario where samples were taken layerby-layer, but not randomly.In addition, the main purpose of using the split-sample method was to investigate the impact of the sample size.Note that, in this study, the semivariogram models used in kriging were those fitted to semivariograms for all 448 observations.Details about the FI and IF approaches used in the study are given below.

FI approach
1. Parameters for three water retention functions (BC, VG, and KSG) were first obtained for all 448 water retention curves S e (h; u i ) using the automatic fitting procedure (Seki, 2007).
2. For each soil hydraulic parameter, unknown values were estimated at all 448 locations for the leaveone-out cross-validation, and at 249 locations for the split-sample method.Following two interpolation approaches were compared: (a) Ordinary kriging (OK) was used for all parameters.
(b) While OK was used for θ s and θ r , indicator kriging (IK) was used for the shape parameters.
The former is referred to as the FI OK approach, while the latter will called the FI IK approach for the remainder of the paper.
3. For each model and approach, a discrepancy between observed water retention curves S e (h; u i ) and those calculated from estimated parameters S e (h; u i ) was obtained at each location by computing differences in water contents at the eleven corresponding pressure heads.Mean absolute error (MAE) and mean error (bias ME) are used in this study in the following way: where N is 448 for the leave-one-out cross-validation and 249 for the split-sample method, θ u i ; h j is the water content at the location u i calculated from estimated parameters for the pressure head h j , and θ u i ; h j is the observed water content at the location u i for the pressure head h j .

IF approach
1. Unknown water content values corresponding to eleven pressure heads were estimated using OK at 448 locations for the leave-one-out cross-validation and 249 locations for the split-sample method.
2. Using the automatic fitting procedure (Seki, 2007), parameters for three WRC models (i.e., BC, VG, and KSG) were obtained at each estimated location.To account for estimation errors when fitting retention models, each h-θ pair was weighted based on the reciprocal of the kriging variance after all kriging variances were standardized by the total sills of the semivariograms.More weights are therefore given to those data points for which the error variances were small.
3. For each model, the discrepancy between observed water retention curves and those calculated from parameters obtained in the previous step was obtained by computing MAE and ME, similar to what was done for the FI approach.
4. The impact of the number of data pairs used to construct WRC when estimating parameters was investigated by repeating 1 through 3, using only the following six pressure heads: 0, −20, −80, −200, −1000, and −15 000 cm.This approach is referred to as IF6 in the remainder of the paper, while the approach considering all eleven pressure heads is referred to as IF11.The number of data pairs used in the IF approach is arbitrary.The optimum number depends on many factors, including soil type.In this study, instead of finding the optimum number of data pairs for this particular data set, we have investigated how prediction errors change when the number of data pairs is reduced to about half of the original number.

Results and discussions
The results are divided into two sections.The first section summarizes experimental semivariograms computed from water retention data and model parameters.In the second section, different approaches (IF and FI) are compared in terms of describing the spatial distribution of soil hydraulic parameters.pressure heads, h 1 − h 11 , with fitted geometric anisotropy models are depicted in Fig. 8.All these semivariograms were calculated using 448 observations.Except for α v , there is clear anisotropy in all semivariograms with major spatial continuity observed in the horizontal direction.While horizontal semivariograms, are generally well structured, in most cases vertical semivariograms fluctuate a lot and are not smooth.There are not a sufficient number of pairs in the vertical direction to obtain well structured semivariograms.Existence of soil horizons would also not result in spatial correlation of soil variables in the vertical direction at the scale of observations.All experimental horizontal semivariograms were fitted using either exponential or spherical models, and they all display a clear nugget effect.Ranges and sill values vary depending on the variable.For example, as expected, semivariograms for the saturated (θ s ) and residual (θ r ) volumetric water contents are similar for all three models (Fig. 7, top two rows).However, other parameters have semivariograms of different shapes, which was expected as the spatial distribution of these parameters varies.Geometric anisotropy models, in which the same nugget and sill values are used for different directions, were fitted for vertical semivariograms.Nugget effects may be underestimated using this approach, since the number of data pairs in the vertical direction is not sufficient.Using a geometric anisotropy model instead of a zonal anisotropy model, which can be used to account for difference in spatial correlation between vertical and horizontal directions, is justified and needed since kriging systems requires the same nugget effect in different directions (Gringarten and Deutsch, 2001).In practice, unless there are conclusive physical explanations for phenomena, it is common to use geometric anisotropy models.

Semivariograms
As for the semivariograms of water contents at particular pressure heads (Fig. 8), the ranges of water contents near saturation in the horizontal direction (i.e., a major range) are generally larger than those for drier conditions.This suggests that the spatial continuity of larger pores is more profound than that of smaller pores.When pressure heads are smaller than −40 cm, the shape of the semivariograms becomes almost identical.From the Laplace equation for capillary rise, the radius of a pore corresponding to a pressure head of −40 cm can be calculated to be about 0.035 mm (assuming that the contact angle is equal to 0 • ).Therefore, it can In the IF approach, the KSG model was fitted to the constructed water retention curve (red).In the FI approach, KSG model parameters were kriged using ordinary kriging to obtain the FI-KSG curve.be expected that spatial distributions of pores smaller than 0.035 mm are similar in this soil.

Prediction performance
Using both the FI and IF approaches, water retention curve model parameters can be obtained at any location.Figure 9 shows differences between the FI and IF approaches in estimating a water retention curve using the KSG model at a selected location, (x, z)=(7.25,2.16), where observed retention data (squares) are available.Due to the exactitude property of kriging, data at this location were not included in kriging (i.e., the leave-one-out cross-validation).While the blue line is obtained with parameters predicted by the FI OK approach, the red line is calculated with parameters obtained by the IF11 approach.In the FI OK approach, model parameters were estimated directly from surrounding conditioning parameters using OK.In the IF11 approach, the water retention curve was first constructed by interpolating water content values corresponding to eleven pressure heads, h 1 − h 11 .
Constructed water retention data are shown in Fig. 9 using red symbols.The KSG model (the solid red line in Fig. 9) was then fitted to the constructed retention data using the aforementioned non-linear least squared method to obtain model parameters.While the model obtained using the FI OK approach could not reproduce the S-shape trend observed in retention data well, the IF11 approach could capture the trend of WRC reasonably well.This confirms that, depending on the chosen approach, resulting WRC models can be significantly different.In the following section, both approaches were compared in a more comprehensive manner using the leave-one-out cross-validation and the split-sample methods.
Figure 10 depicts the mean absolute prediction errors (MAE) calculated for all four approaches, i.e., FI OK , FI IK , IF11, and IF6, using the leave-one-out cross-validation and the split-sample procedures for the three WRC models.In general, both IF approaches resulted in smaller MAE than the FI approaches, regardless of which WRC model was used.Decreases in prediction errors were much larger when the VG model was used, mainly because the prediction performance of the VG model for the FI approach was much worse than that for the other two models.At one location, (x, z)=(12.75,3.21), the fitted parameter α v (=28.16cm −1 ) was three orders of magnitude greater than in the rest of the domain, where the mean of fitted α v is 0.12 cm −1 and the median is 0.04 cm −1 (see the histogram of α v in Fig. 5).Although the fit was acceptable at this particular location (Fig. 11), the large α v value (an outlier) affected the estimation of α v values at surrounding locations, leading to a relatively poor performance of the FI OK approach for the VG model.Indicator kriging was expected to reduce the impact of outliers in the FI approach, as all data are transformed first to binary data, either 0 or 1, depending upon exceedence of the threshold value.However, when E-type estimates were computed, the observed extreme value was taken as a possible maximum value.Because of that, the indicator approach could not improve the prediction.The problem associated with the outlier can be avoided if the IF approach is used.For the other two WRC models, BC and KSG, average prediction errors decreased by about 10-15% when the IF approach  was used in comparison with the prediction errors for the FI approach, where there were no obvious outliers in initially evaluated parameters.
For the FI OK approach, mean absolute errors vary depending on the WRC model used.On the other hand, prediction errors do not depend upon the model used for the IF approach.This is not a surprising result since three models were fitted to the same reconstructed WRCs in the IF approach.Therefore, the prediction errors calculated for the IF approach are mainly the summary of model fitting performance, but do not directly reflect the accuracy of interpolation as those for the FI approach do.
With regard to MAE, there is no distinct trend in whether the IF11 approach is better than the IF6 approach is better.While IF6 outperforms IF11 when the VG model was used for both the leave-one-out cross validation and the splitsample method, IF11 had slightly smaller mean absolute errors or almost the same errors when compared to IF6 for the BC and KSG models.This indicates that reducing the number of θ-h pairs by about half does not necessarily worsen prediction performances, at least for this particular data set, as long as the IF approach is adopted.This result is quite important, since the main reason for preferring the FI approach over the IF approach in most studies is that the number of variables one needs to analyze can be significantly reduced.The number of parameters required in the WRC models used in this study is four, which is significantly smaller than the original eleven θ-h pairs used to construct water retention curves.To use all pairs in the IF approach, one needs to perform kriging eleven times, which requires significantly more effort than performing it only four times in the FI approach.However, it was just shown here that when the number of data pairs is reduced to six, the IF approach is still better than the FI approach.With a little extra work, we can improve the description of the spatial distribution of water retention curves using the IF approach.This may not be true when soils with multimodal pore systems (e.g., aggregated soils) are considered.Such soils may require models that can distinguish between different pore systems, such as the Durner model (Durner, 1994).To obtain parameters for such models, we need to have more θ-h pairs than required by traditional unimodal models.An appropriate number of pairs thus varies depending upon the soil type and can be assessed using the IF approach.This is another advantage of the IF approach over the FI approach.
As for the mean error, there is no obvious trend of overor under-estimation, except for the VG model where mean errors are always positive (meaning that predicted water contents are, on average, greater than observed values).When comparing different approaches, ME is much smaller (i.e., has less bias) for the IF approaches than for the FI approaches, except when the BC model was considered in the split-sample method.Like that for MAE, this result confirms the superiority of the IF approach over the FI approach.While for the BC and KSG models, IF11 resulted in less bias than IF6, for the VG model, IF6 had much smaller biases than IF11 for both the leave-one-out cross-validation and the split-sample method.This indicates again that reducing the number of h-θ pairs from eleven to six in the IF approach does not necessarily worsen the overall performance.
Figure 13 shows spatial distributions of locations where IF11 outperformed FI OK (open squares), and vice versa For the VG model especially, the IF approaches resulted in smaller prediction errors at more than 80% of the sampling locations.Although the difference between IF11 and IF6 is small, IF11 was slightly better than IF6 for BC and KSG, while IF6 was better than IF11 for VG.Overall, these percentages confirm the superiority of the IF approach over the FI approach in predicting WRCs at any given location.

Conclusions
This study compared the performance of two interpolation approaches, a commonly used fit-first interpolate-later (FI) approach and a proposed interpolate-first fit-later (IF) approach, to estimate the spatial distribution of water retention curve parameters.In the former approach, the retention curve model parameters are first evaluated at sampled locations from observed retention data, and then these parameters are interpolated directly for locations without measurements.In the latter approach, retention data are first interpolated for locations without measurements, and then retention curve model parameters are fitted to these estimated retention curves.Three retention models, Brooks and Corey, van Genuchten, and Kosugi models were considered.As for spatial interpolation, a least-squared interpolation technique, known as kriging, was used.For variables that do not show asymmetric distributions, ordinary kriging was used, while for those with asymmetric distributions, indicator kriging was used as well.Standard validation techniques (the leaveone-out cross-validation and the split-sample methods) were used to compare the two approaches in terms of reproducing retention curves.
Both validation results showed that, regardless of the retention model selected, the IF approach had lower prediction errors (mean absolute errors) when compared to the FI approach.Improvements were especially substantial for the VG model because an extreme value in initially estimated parameters led to the poor estimate of that parameter at surrounding locations in the FI approach.There is always a risk of including parameter outliers in the estimation process as long as the FI approach is used.In the IF approach, such a problem rarely occurs.In addition, the IF approach performed better than the FI approach even when the number of data pairs used to construct retention curves was reduced from eleven to six.This makes the IF approach comparable to the FI approach in terms of workload, because the number of variables used in the FI approach is four, while that in the IF approach would be only six.The FI approach has been preferred mainly because one can reduce the computational time and workload.However, as this study shows, with a little bit of extra work, one can achieve a much better estimation of the spatial distribution of water retention curves by using the IF approach.In the future, the IF approach should be applied to estimate unsaturated hydraulic conductivity as well.
Overall, this study shows the superiority of the IF approach over the FI approach in evaluating the spatial distribution of water retention curves, which is the initial but critical step in modeling variably-saturated water flow and solute transport in field-scale heterogeneous soils.This study shows that interpolating retention model parameters, such n and α in the VG model, should be done with extreme care.In future studies, the actual impact of the FI approach on water flow and solute transport still needs to be investigated.

FitFig. 1 .
Fig. 1.A difference between the fit-first interpolate-later (FI) and interpolate-first fit-later (IF) approaches used in this study.

Fig. 5 .
Fig. 5. Sample histograms of water retention curve model parameters for BC (left), VG (middle), and KSG (right) models.The saturated and residual water contents are in the top two rows, two shape parameters are in the bottom two rows.

Figure 7 Fig. 7 .
Figure7shows the experimental semivariograms of WRC model parameters with fitted geometric anisotropy models.Semivariograms of water contents corresponding to eleven 31

Table 1 .
Water retention, S e , models used in this study.

Table 2 .
Percentages (%) of locations where mean absolute errors (MAE) for the IF approach were smaller than those for the FI approach in the leave-one-out cross-validation.closedsquares), in terms of MAE in the leave-one-out crossvalidation.There is no specific trend observed, such as a cluster of open or closed squares.This means that the prediction of water retention curves was spatially unbiased.Proportions of open squares for both IF11 and IF6 (spatial distributions are not shown) are summarized in Table2for all cases.For all models, the IF approaches outperformed the FI approaches at more than 60% of the sampling locations. (