VZJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 23 August 2007
Published in Vadose Zone J 6:638-650 (2007)
DOI: 10.2136/vzj2006.0077
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Sheikh, V.
Right arrow Articles by van Loon, E. E.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Sheikh, V.
Right arrow Articles by van Loon, E. E.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Sheikh, V.
Right arrow Articles by van Loon, E. E.
Related Collections
Right arrow Watershed and Landscape Processes
Right arrow Effective Parameters
Right arrow Water Flow Models

ORIGINAL RESEARCH

Comparing Performance and Parameterization of a One-Dimensional Unsaturated Zone Model across Scales

Vahedberdi Sheikha and E. Emiel van Loonb,*

a Erosion and Soil & Water Conservation Group, Wageningen Univ., Nieuwe Kanaal 11, 6709 PA Wageningen, The Netherlands, and Gorgan Univ. of Agricultural Sciences & Natural Resources, Gorgan, Iran
b Inst. for Biodiversity and Ecosystem Dynamics, Univ. of Amsterdam, Nieuwe Achtergracht 166, 1018 WV Amsterdam, The Netherlands
* Corresponding author (vanloon{at}uva.nl).

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.


Received 29 May 2006.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 Materials and Methods
 Results
 Discussion and Conclusions
 REFERENCES
 
The utility of an unsaturated zone soil moisture model is not only its ability to describe the soil moisture dynamics at a given point but also the possibility to generalize the results to larger areas. In this study we investigated the predictive performance of a one-dimensional unsaturated zone soil moisture model when applied at point, field, response unit, and catchment scales, using detailed field observations from a 0.42-km2 catchment in the Netherlands. Our main question was how model parameterization and model performance could be compared across these scales. We considered two different calibration–validation schemes and three performance statistics. In all cases we applied the same Levenberg–Marquardt optimization scheme. Differences between calibration–validation schemes (interpolation vs. extrapolation) were surprisingly small. Using one particular model parameterization across the various aggregation levels, the optimal Mualem–van Genuchten parameters for a coarser aggregation level can be derived from an underlying level by simple arithmetic averaging. The different performance indices (RMSE, index of agreement, and Nash–Sutcliffe coefficient) were highly variable between observation locations and for different aggregation levels. Overall, the indices were more favorable at higher aggregation levels, and in correspondence with errors reported in comparable studies. The unsaturated-zone model did not, however, provide satisfactory predictions of independent flux observations, in this case daily catchment discharge. Moreover we did not succeed in deriving a meta-model to scale model performance indices with aggregation level. Our case study therefore supports the view that multiscale calibration studies that use both state and flux observations are required to compare results from unsaturated zone models at different aggregation levels.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
Materials and Methods
 Results
 Discussion and Conclusions
 REFERENCES
 
Soil moisture content in the top 1 to 2 m of the earth's surface controls the success of agriculture and regulates partitioning of precipitation into surface runoff, evapotranspiration, and drainage into deeper groundwater storage. Therefore, it is considered a key parameter in meteorological or hydrological studies (Walker, 1999). Despite its importance, it is not monitored regularly. Because it is highly variable both spatially and temporally, it is not cost effective to measure soil moisture routinely using conventional methods. One possibility to cover spatial variability of soil moisture across large areas is the application of remote sensing techniques; however, the availability of satellite images at a low temporal frequency and at a relatively coarse spatial resolution (while only covering a very shallow surface layer) still limits their applicability. The efforts to enhance satellite sensors and extract better information are ongoing.

Such efforts greatly benefit from well-designed in situ soil moisture measurements, especially to calculate accurate spatially averaged estimates (Walker, 1999; Western and Grayson, 1998). Dynamical unsaturated zone models are generally believed to be appropriate tools to obtain such averages (Giacomelli et al., 1995; Teuling and Troch, 2005), while statistical interpolation techniques, when used in isolation, are only of limited value for soil moisture monitoring (e.g., Western and Grayson, 1998; Western et al., 2004). During the past few decades, a wealth of experimental studies have been done on infiltration and water movement in soil profiles (Scotter et al., 1982; Devaurs and Gifford, 1984; De Roo and Riezbos, 1992; Bagarello and Iovino, 2003). These results promoted considerable progress in the conceptual understanding and mathematical description of water flow and solute transport processes in the unsaturated zone. This is illustrated by the variety of analytical (Ahuja et al., 1981; Hurley and Pantelis, 1985; Stagnitti et al., 1992; Selim, 1987) and numerical (Beven, 1981; Rogers and Selim, 1989; Jackson, 1992; Simunek et al., 1992) models that have been developed.

Nearly all of the currently available water flow and solute transport simulation models (e.g., SWAP, SWMS-2D, WAVES, HYDRUS) are based on the numerical solution of Richards' equation with different schemes. Models of this type have been shown to be useful tools for extrapolating information from a limited number of field experiments to different soil, crop, land management, and climatic conditions (e.g., Yu et al., 2000). Apart from obvious omissions from these models (e.g., effects of temperature, freezing, or macropore flow), more intricate problems of scale also remain (Klemes, 1983; Blöschl and Sivapalan, 1995; Harter and Hopmans, 2004; Lin et al., 2006). It has also been outlined in several studies that a Richards' equation model is not an appropriate parameterization at every resolution (Downer and Ogden, 2004; Van Dam and Feddes, 2000). Yet the questions remain at exactly which scales it does provide a suitable model, and how model performance should be evaluated in this context. The latter question comprises the choice of a calibration–validation scheme as well as a performance criterion. This study was meant to contribute in answering these questions by demonstrating, with a case study, the effect of different calibration–validation schemes, and how model error criteria may vary for different levels of aggregation.


   Materials and Methods
TOP
 ABSTRACT
 INTRODUCTION
 Materials and Methods
Results
 Discussion and Conclusions
 REFERENCES
 
Study Area
The data used in this study originate from the Catsop catchment. The catchment is situated in the hilly loess region of South Limburg in the Netherlands (inside the box 50.93–50.94° N, 5.77–5.79° E; Fig. 1 ). It has an area of 0.42 km2, and is almost entirely used for agricultural purposes. Elevation ranges from 79 to 112 m above sea level, and the terrain is undulating with gently to moderately sloping topography. About 86% of the slopes have a gradient of ≤10% and 3.5% of the slopes are steeper than 15%. The catchment has a dominant flow direction to the west, the direction of the River Meuse. During major storms, 3 to 30% of the total rainfall reaches the catchment outlet (De Roo 1993). The soils of the catchment have been formed on aeolian deposits (Ice age, 0.01–0.2 million yr ago) and are typical of the loess areas in northwestern Europe. The depth of these deposits in the Catsop catchment exceeds 5 m.


Figure 1
View larger version (20K):
[in this window]
[in a new window]

  		FIG. 1. Geographical position of study area and location of the measurements (Gr = grass, WW = winter wheat, CST = conservation tillage, CVT = conventional tillage, YM = yellow mustard, NO = 2-yr-old orchard, and OO = 5-yr-old orchard).

 
Within this small catchment, four land use types were distinguished: arable (79.5%), orchard (7.9%), grassland (11.8%), and infrastructure (0.8%). Figure 2 presents the land use pattern of the Catsop catchment during the winter season of 2003–2004. Arable land is cultivated mainly with winter wheat (Triticum aestivum L.), spring barley (Hordeum vulgare L.), sugarbeet (Beta vulgaris L. ssp. vulgaris), potato (Solanum tuberosum L.), and yellow mustard (Sinapsis alba L.). Yellow mustard is a second crop, functioning as an erosion prevention measure. It is planted after harvesting the cereals at the end of the summer, and later in the season chopped into small pieces and spread on the land surface as green manure and protection cover. Grasslands are used in two ways. One field is grazed by cows from April until October (<0.5 cows ha–1) and the other fields are harvested by machinery. During the last 5 yr, the area in orchard has increased from nil to about 8%. The only available infrastructure within the catchment are a few tarred roads and one ditch near the outlet, which runs parallel to the main road (shown in Fig. 1). Table 1 lists some of the properties of each land use.


Figure 2
View larger version (77K):
[in this window]
[in a new window]

  		FIG. 2. Land use in the Catsop catchment during the winter season of 2003–2004.

 

View this table:
[in this window]
[in a new window]

  		TABLE 1. General information on the measurement locations.

 
The climate is temperate humid, with a mean annual precipitation of approximately 740 mm. Precipitation occurs mainly as rainfall and is evenly distributed throughout the year; however, the rainfall patterns in winter and summer are different. Summer events are shorter and more intense, while winter events are on average longer and less intense.

Data Collection
Soil moisture was monitored in 15 tubes at 1-m depth with a time domain reflectometry system (Trime-FM, IMKO Micromodultechnik, Ettlingen, Germany) once every week for the period of November 2003 to September 2004. Measurements were taken for five depth intervals: 0 to 20, 20 to 40, 40 to 60, 60 to 80, and 80 to 100 cm from the surface downward. The measurement technique integrates soil moisture content throughout each entire 20-cm layer. Due to large errors in the measurements of the 80- to 100-cm layer, its values were not used in this study. The locations of these tubes and other measurement equipment are shown in Fig. 1. The topographic characteristics of each observation location are given in Table 1. For each land use type and management practice, at least two tubes were installed. Two tipping bucket rain gauges and one small stand-alone weather station were installed as well. All observed meteorological variables did agree well with the observations at the Beek weather station (50°55' N, 5°47' E), at a distance of <2 km from the study area. The evapotranspiration rate has been calculated with the Penman–Monteith equation using the daily meteorological data of the Beek station. Discharge was measured at the catchment outlet using a partial flume with a capacity of 950 L s–1 and a stilling well with a vertical float recorder. This station was initially installed by Wageningen University in 1987. Later the water board "Roer en Overmaas" took over the responsibility of this station. It has been used sporadically during the last two decades. For the purpose of this study, discharge measurements were collected at 10-min intervals during the period November 2003 to September 2004. Although the discharge observations may be relevant to questions about (multiscale) unsaturated zone soil moisture modeling, these data were only used to illustrate the (in)adequacy of models to predict runoff when calibrated on soil moisture observations only.

Modeling
The Soil–Water–Atmosphere–Plant (SWAP) model was used to describe soil water movement in the unsaturated zone of our study catchment. The SWAP model is the successor of the agrohydrological model SWATRE (Feddes et al., 1978; Kroes and Van Dam, 2003). It assumes that the main flow processes occur in the vertical direction only (hence it is a one-dimensional model) and considers soil moisture movement in close interaction with crop growth. The SWAP model can simulate crop growth with either a simple or a detailed crop growth model. In the simple crop growth model (used in this study) the measured leaf area index, crop height, and rooting depth are prescribed as a function of crop development stage, which either is controlled by a temperature sum or is linear in time. The simple module does not simulate the effect of water or salt stress. The detailed crop growth model comprises an implementation of WOFOST (Hijmans et al., 1994). This detailed model does include the effect of water and salt stress on crop growth. For more detailed information, see Kroes and Van Dam (2003) and Van Dam (2000).

Soil Water Flow
Transient soil water flow is simulated by the well-known Richards' equation.

Formula 1[1]
where C = d{theta}/dh is the soil water capacity (L –1), h is the soil water pressure head (L), K is the hydraulic conductivity (L T –1), Sa is the root water uptake rate (T–1), and z is the vertical coordinate (L) (positive upward).

The SWAP model solves Eq. [1] numerically using an implicit finite difference scheme as described by Van Dam and Feddes (2000). The numerical solution of Eq. [1] is subject to specified initial and boundary conditions and soil hydraulic functions, which relate {theta}, h, and K. The Mualem–van Genuchten equations (Mualem 1976; van Genuchten, 1980) are implemented to relate {theta}, h, and K:

Formula 2[2]

Formula 3[3]
where {theta}res is the residual water content (L3 L –3), {theta}sat is the saturated water content (L3 L –3), Se is the relative saturation (dimensionless),{alpha} (L –1) and n (dimensionless) are empirical shape factors, Ksat is the saturated hydraulic conductivity (L T –1), and {lambda} (dimensionless) is an empirical coefficient.

The upper boundary conditions are controlled by the rates of potential evapotranspiration, irrigation, precipitation, and interception. Potential evapotranspiration (ETp) is estimated by the Penman–Monteith equation, using daily weather data of solar radiation, air humidity, wind speed, and air temperature as well as crop characteristics such as minimum resistance, reflectance (albedo), and crop height (Allen et al., 1998).

When soil surface is partly covered by a crop, ETp is partitioned into soil evaporation rate, Ea (L T –1), and potential transpiration rate, Tp (L T –1). This partitioning is achieved through the leaf area index, LAI (L2 L –2), or soil cover, Sc (dimensionless), as a function of crop development stage (Goudriaan, 1997; Belmans et al., 1983):

Formula 4[4]

Formula 5[5]

Formula 6[6]
where ETp is the potential evapotranspiration rate of the wet crop as calculated with the Penman–Monteith equation, neglecting crop aerodynamic resistance, kgr is the extinction coefficient for global solar radiation (0.39 for common crops), and Pi is the interception rate (L), which is estimated with the Von Hoyningen-Hüne (1981) formula:

Formula 7[7]
where Pgross is gross precipitation (L), a is an empirical coefficient (L), and b is the soil cover fraction (~LAI/3; dimensionless). Generally for ordinary crops, it is assumed that a = 0.25.

Under wet soil conditions, actual soil evaporation E (L T –1) is governed by the atmospheric demand, and equals Ep. When soil becomes drier, the actual evaporation decreases to a rate that is quantified with Darcy's law as

Formula 8[8]
where K1/2 is the average hydraulic conductivity (L T –1) between the soil surface and the first node, hatm is the soil water pressure head (L) in equilibrium with the air relative humidity, h1 is the soil water pressure head (L) of the first node, and z1 is the soil depth (L) at the first node.

The SWAP model also has an option to choose the empirical evaporation functions of Black et al. (1969) or Boesten and Stroosnijder (1986) to avoid the serious limitations of the Emax model given in Eq. [8]. In this study, Eq. [8] has been used. It is still not clear to what extent the soil hydraulic functions, which usually represent a top layer of a few decimeters, are valid for the top few centimeters of a soil, because the top few centimeters of soil is subject to splashing rain, dry crust formation, root extension and various cultivation practices (Van Dam, 2000). The SWAP model determines Ea by taking the minimum value of Ep, Emax, or an actual evaporation rate calculated based on the empirical functions.

The lower boundary condition is defined by the fluxes at the bottom of the soil profile. Since groundwater level at the catchment is deeper than 2 m for all the measurement locations, we applied the free-drainage condition to calculate the bottom flux qbot (L T –1):

Formula 9[9]

Parameter Estimation and Model Validation
Generally speaking, models of soil water flow are highly sensitive to soil hydraulic characteristics as parameterized via the soil hydraulic retention {theta}(h) and hydraulic conductivity, K({theta}) functions (Van Dam, 2000; Mertens et al., 2005). In SWAP, the Mualem–van Genuchten equations (Eq. [2–3]Go) are used to describe the soil hydraulic functions. In total, there are six parameters ({theta}res, {theta}sat, {alpha}, n, {lambda}, and Ksat) in these equations (the parameter m is calculated by 1 – 1/n). Of these parameters {theta}sat and Ksat have a clear physical meaning and can in theory be measured directly. Due to practical difficulties in the measurement of Ksat, however, its value is usually derived through calibration. We have divided the profile into three layers: 0 to 20, 20 to 40, and 40 to 80 cm (from the soil surface downward). This division was made based partially on our field observations, and partially on the basis of previous studies of soil properties in the Catsop catchment (Stolte et al., 1994). During augering to install TRIME access tubes, at most of the augering locations a compacted layer was discovered at a depth of between 15 and 50 cm in the soil profile. A penetrometer resistance survey with an electronic penetrologger confirmed the presence of this compacted layer throughout the catchment. Stolte et al. (1994) noticed a considerably lower hydraulic conductivity in the 20- to 40-cm layer, compared with the 0- to 20-cm layer. The {theta}res is usually assigned a value near zero. In our study, it was given a fixed value of 0.01 for all three layers. The empirical shape parameter {lambda} was based on the values found by Ritsema et al. (1996): a value of 0 was used for the top two soil layers and –1.17 for the third layer. Both {theta}res and {lambda} show low sensitivities in relation to soil water flow (Van Dam 2000; Singh 2005). The shape parameters {alpha} and n can be estimated by fitting measured water retention data or can be derived via calibration of a flow model to observed water content and pressure head data. We chose to derive {alpha} and n via calibration. In summary, there are three uncertain parameters per soil layer that have to be derived by calibration: {alpha}, n, and Ksat (nine parameters in total). For our system, this large number of parameters leads to an ill-posed inverse problem (Kool and Parker, 1988). Thus we decided to fix three of the nine parameters. We followed two strategies to decrease the number of parameters to be optimized. The first was a preliminary optimization procedure with nine parameters to find less sensitive soil hydraulic parameters, and subsequently continued optimization with six parameters. The second strategy was to optimize with different combinations of two of three parameters ({alpha}, n, and Ksat) for each soil layer. The initial values for the nine parameters were taken from Ritsema et al. (1996) and Stolte et al. (1994). Parameters of the first layer have been extracted from measurements reported in Stolte et al.(1994). The parameter values for the second and third layers were taken from Ritsema et al. (1996). They derived the Mualem–van Genuchten parameters on the basis of soil profile descriptions.

For calibration, we used the PEST calibration software (Doherty, 2002), model-independent optimization software that requires the user to provide a communication link between the specified model and PEST. We used a communication link that is schematically shown in Fig. 3 . The PEST software uses the Levenberg–Marquardt algorithm to find an optimal parameter set that minimizes the differences between model results and observations, using the following objective function {Phi}(b):

Formula 10[10]
where b is the vector with fitting parameters, yj*(ti) denotes the observation of type j at time ti, yj(b,ti) is the corresponding model prediction, and wj is the weighting factor. In the case of random observation errors only, according to maximum likelihood the weighting factor wj should be equal to the standard deviation of the observation error of observation type j. Since we have only one observation type in our optimization process, the objective function decreases to

Formula 11[11]
where {theta}obs(ti) is the observed soil moisture content at time ti, N is the number of observations, and {theta}sim(b,ti) is the simulated values of soil moisture using an array with parameter values b at time ti. In comparing the observed with simulated soil moisture values, the observations for the 0- to 20- and 20- to 40-cm layers match the layers used in the model. The observations for the 40- to 60- and 60- to 80-cm layers were averaged to match the 40- to 80-cm layer in the model.


Figure 3
View larger version (15K):
[in this window]
[in a new window]

  		FIG. 3. Communication between SWAP and PEST during parameter calibration.

 
To evaluate model performance and uniqueness of the optimized parameters set, we used two calibration–validation schemes. In the first scheme, we used the first part of the available data for every tube (all observations before 5 February) for calibration and the second part of the data (all data after 5 February) for validation. This calibration scheme is similar to that applied by Singh (2005). In the second scheme, we used observations, ranked by date, with uneven rank numbers for calibration and the remaining ones for validation. It is generally believed that the first is a more appropriate test for a model that leads to higher prediction errors for the validation set (McCuen and Snyder, 1986). Given that all model states occur more or less random in time (at least at a monthly time scale), the first represents an extrapolation scheme, whereas the second is an interpolation scheme; we refer to them as such. For both calibration–validation schemes, the initial parameter values were randomly chosen from a hypercube around the prior parameter set to evaluate whether the results were unique.

Evaluation of Simulation Quality
There are different criteria for evaluation of a model performance. A review of some recent soil moisture simulation literature (Hymer et al., 2000; Mapfumo et al., 2004; Wegehenkel, 2005) indicates that four criteria have been used most frequently. These are the Pearson product-moment correlation coefficient (R), root mean squared error (RMSE), the index of agreement (IA) (Willmott, 1981), and the Nash–Sutcliffe (NS) coefficient (Nash and Sutcliffe, 1970). For evaluation of SWAP results, however, usually the RMSE has been used (Singh, 2005; Crescimanno and Garofalo, 2005). We evaluated all these criteria. The equations for the evaluation criteria are as follows:

Formula 12[12]

Formula 13[13]

Formula 14[14]

Formula 15[15]
where Formula 15sim and Formula 15obs are the mean values of the simulated and measured values, respectively. The range of IA is between 0 and 1, and the closer it is to 1, the better is the fit between measured and simulated model outputs. On the basis of the above criteria we would reject a model if NS < 0.5, NS < 0, or if RMSE > 0.1 (the last value is an approximate upper bound on the basis of a literature review, where similar data and models were used; see below).

In addition to considering only soil moisture to evaluate model quality, we considered the adequacy of the SWAP model to predict daily catchment discharge (while catchment discharge data was not used in any way for model calibration). This adequacy was measured by RMSE and compared with the RMSE when predicting discharge by daily rainfall depth only. The choice to use discharge at a daily rather than a more detailed time scale was made because SWAP lacks a routing module. In SWAP, runoff from a land unit is assumed to go straight to the catchment channel, and runoff for the total catchment is calculated by a simple summation of all land units.

The Relation between Level of Aggregation and Model Performance
After evaluation of the model performance for each individual tube (coded as AggL0), the performance of the model was evaluated at three different aggregation levels. In the first aggregation level (AggL1), the tubes located in each field were put in the same group. There were, at most, two tubes in each field. For the second aggregation level (AggL2), the tubes with the same land use were put in one group. We assumed three different land use types: cropland, grassland, and orchard. At the third aggregation level (AggL3), all 15 tubes were put in one group. Values were aggregated separately for each layer. At Aggregation Levels 1, 2, and 3, the layer-specific values used for model calibration were the average values of the point measurements within each aggregation level. Thus, the depth-specific average across an aggregation level is a substitute for larger scale measurements, e.g., remote sensing data with a resolution equal to the larger (aggregated) level.


   Results
TOP
 ABSTRACT
 INTRODUCTION
 Materials and Methods
 Results
Discussion and Conclusions
 REFERENCES
 
General
The many calibration–validation runs that were made in this study led to a large number of time series of predicted vs. observed soil moisture. Typical results are shown in Fig. 4 to 7. Figure 4 shows results for the unaggregated data. Only the values for the upper soil layer (0–20 cm) of the four tubes in the orchards are shown. Figure 5 shows the predictions per field (AggL1) along with the unaggregated observations. Figure 6 shows the predictions per land use (orchard in this case), and Fig. 7 shows the predictions for the entire catchment. In this case, the extrapolation scheme was applied (i.e., the period till 5 February was used for calibration and the remaining period was used for validation). As expected, the model error is slightly smaller in the calibration period than in the validation period. Also, the model error increases with aggregation level. The graphs clearly indicate that it will be very hard if not impossible to relate predictions at AggL2 and AggL3 to point observations. Figures 4 to 7 are merely used as an illustration of the results and the complexity of the multivariate output (different time instants, different soil layers, different aggregation levels). The full complexity is not shown because several land uses and other soil layers are omitted. With regard to general simulation quality, it is interesting to note that the model results are unbiased at all aggregation levels when considered for periods of a week or longer.


Figure 4
View larger version (21K):
[in this window]
[in a new window]

  		FIG. 4. Comparison of observed and predicted soil moisture in the top layer (0–20 cm) for tubes within orchards for unaggregated observations and predictions (Aggregation Level 0). The graphs in the left column show time series of observations (symbols) and predictions (lines). The vertical line in each time series separates calibration (left side) and validation (right side) period. The graphs in the right column are scatter graphs of observation vs. predictions. Similar graphs can be made for the second and third layer of the model; the overall pattern in (dis)agreement between observation and prediction is the same for all three layers. The labels NO up, NO down, OO up, and OO down refer to the tubes depicted in Fig. 1.

 

Figure 7
View larger version (11K):
[in this window]
[in a new window]

  		FIG. 7. Comparison of observed and predicted soil moisture in the top layer (0–20 cm). Observations are unaggregated (so these are the same as those in Fig. 4–6GoGo) and predictions are for Aggregation Level 3 (catchment average). The graph at the right shows unaggregated observations vs. aggregated predictions.

 

Figure 5
View larger version (14K):
[in this window]
[in a new window]

  		FIG. 5. Comparison of observed and predicted soil moisture in the top layer (0–20 cm) for the new orchard (NO) and old orchard (OO). Observations are unaggregated (so the observations for, e.g., OO up and OO down are the same as those in Fig. 4) and predictions are for Aggregation Level 1. The graphs in the right column are scatter graphs of unaggregated observations vs. aggregated predictions.

 

Figure 6
View larger version (11K):
[in this window]
[in a new window]

  		FIG. 6. Comparison of observed and predicted soil moisture in the top layer (0–20 cm). Observations are unaggregated (so these are the same as those in Fig. 4 and 5) and predictions are for Aggregation Level 2 (here, orchard). The graph at the right shows unaggregated observations vs. aggregated predictions.

 
Calibration and Parameter Uniqueness
As stated above, two strategies were evaluated to decrease the number of parameters for optimization. Both strategies showed that optimization with {alpha} and n for each soil layer gives the best results. For the first strategy, the least sensitive parameter in at least one of the soil layers was Ksat in 70% of the cases. Our initial hypothesis was that this low sensitivity of SWAP to Ksat could be attributed to the use of daily mean rainfall intensity as an input to the model. A test with rainfall data at a fine resolution, however, led to the same (low) sensitivities for Ksat, and hence we found no supporting evidence for this hypothesis. Also De Roo et al. (1996) investigated the sensitivity to K({theta}) for the Catsop catchment (using the Lisem model) and also noted the relatively low sensitivities of model output to Ksat (a 60% model change due to a 20% change in Ksat).

In the second strategy, the objective function ({Phi}, Eq. [10]) was minimized for the calibration period while only changing combinations of two parameters. The optimization process with the combination of {alpha} and Ksat for each soil layer led to the worst result. These two parameters also appeared to be very strongly correlated. Probably an approach using scale factors (which are related to the hydraulic parameters) as in Rockhold et al. (1996) and Zhang et al. (2004) would have led to better results. With our present calibration schemes, the best results were obtained for the combination of {alpha} and n. The results for the combination of n and Ksat were for most of the tubes comparable with the results of combining {alpha} and n. But the standard deviation of Ksat was generally much higher than that of {alpha}. Also, the relative sensitivity of Ksat was lower. Therefore we calibrated six parameters: {alpha}1, n1(the first layer), {alpha}2, n2 (the second layer), {alpha}3 and n3 (the third layer). The values of Ksat were obtained from Ritsema et al. (1996) and Stolte et al. (1994). When calibrating with these six parameters, up to 20% uncorrelated Gaussian noise to the initial parameter values (by random selection within a hypercube) always led to the same optimal parameter sets. Adding a 5% uncorrelated Gaussian noise to the soil moisture observations and the rainfall data led to only minor changes in the parameter values. Hence our choice to represent the optimal parameter sets as a unique set rather than a parameter range seems to be justified. Apparently the information content of the observations in different soil layers, together with the objective function (Eq. [11]), led to a unique optimum. Notwithstanding this result, we consider the optimal parameter sets as a realization from a joint probability distribution rather than some fixed value. The reason why we nonetheless represent the parameter sets as unique values in this study is to limit the complexity of our analysis.

The Effect of Different Validation Schemes
The rationale for evaluating two different validation schemes (viz. extrapolation and interpolation) was to answer the question to what degree a validation scheme influences model performance. The two schemes chosen in this study are in fact the most extreme examples for extrapolation and interpolation. We therefore expected that the results would generalize to intermediate cases. Both schemes used the same initial conditions, parameter values, and optimization routine. The resulting parameter sets and performance statistics are indeed different for the two schemes (see Tables 2 and 3). In the interpolation scheme, {Phi} and R are more similar between calibration and validation than in the extrapolation scheme. In fact, for interpolation there is no significant difference in {Phi} and R between calibration and validation (evaluated with a paired t test), whereas for extrapolation there is a significant difference (P < 0.05). Also, the parameter values across the soil layers showed a more regular pattern for interpolation than for extrapolation. The model performance for the extrapolation scheme (considering validation {Phi} and R) is somewhat better, however, than for interpolation (significant at the 0.05 level, according to a paired t test). The almost equal performance for calibration–validation via interpolation or extrapolation is somewhat counterintuitive. It is generally assumed that interpolation should lead to higher model performance than extrapolation, especially in a situation where some boundary conditions (like temperature and rainfall pattern in our system) and system states (hydraulic properties of the top layers) change considerably across the time scale of a few months (see Ritsema et al., 1996). Our explanation for this result in our study is that the observation frequency for the interpolation scheme (observation times >10 d apart) is too low to capture the essential dynamics of a dry-down cycle. With higher observation frequencies, we would expect better performance indices for the interpolation scheme than the extrapolation scheme.


View this table:
[in this window]
[in a new window]

  		TABLE 2. Soil hydraulic parameter values (saturated water content {theta}sat, empirical shape factors {alpha} and n, empirical coefficient {lambda}, and saturated hydraulic conductivity Ksat, see Eq. [2] and [3]) after optimization according to the extrapolation calibration–validation procedure. Objective function {Phi} is the optimization criterion to be minimized (Eq. [10] and [11]), and R is the correlation coefficient between observed and predicted soil moisture (Eq. [12]).

 

View this table:
[in this window]
[in a new window]

  		TABLE 3. Soil hydraulic parameter values (saturated water content {theta}sat, empirical shape factors {alpha} and n, empirical coefficient {lambda}, and saturated hydraulic conductivity Ksat, see Eq. [2] and [3]) after optimization according to the interpolation calibration–validation procedure. Objective function {Phi} is the optimization criterion to be minimized (Eq. [10] and [11]), and R is the correlation coefficient between observed and predicted soil moisture (Eq. [12]).

 
Parameter Values across Scales
Very little structure is seen in the parameter values when moving from point to catchment scales. When moving from lower to higher aggregation levels, the parameter values change in a more or less random fashion, whereby the ordering over layers is not related between aggregation levels. The change of parameter values across scales is shown in Fig. 8 . The difference between parameter values from a lower to a higher aggregation level is in all cases insignificant (tested with pairwise Mann–Whitney U tests). Stated differently, the mean parameter values obtained at a lower aggregation level were in all cases quite close to the values obtained at a higher aggregation level. This result is similar to that in studies by Pachepsky and Rawls (1999) and Zhu and Mohanty (2002).


Figure 8
View larger version (12K):
[in this window]
[in a new window]

  		FIG. 8. The values of Mualem–van Genuchten parameters {alpha} and n per layer for different aggregation levels. Each dot is a parameter value for a particular land unit, each + represents the mean of the parameter values per aggregation level. The dots are jittered per aggregation level to distinguish the individual parameter values. The parameter values for Aggregation Level 0 correspond with the values in Table 2.

 
Relations between the Different Performance Indices
Our study confirms that using only one index to evaluate model performance does not give a comprehensive view (Willmott, 1981; Legates and McCabe, 1999; McCuen and Snyder, 1986; Krause et al., 2005). This is illustrated by the apparent absence of a relation between R and {Phi} (see Tables 2 and 3) as well as IA, NS, and RMSE (Table 4). In Table 4, the statistics for each layer, the mean of four layers, and the mean of four layers and all observation times are given per observation location. Especially the Nash–Sutcliffe coefficient sometimes deviates from the other statistics (see, e.g., the indices for land uses OO.u and CVT.u in Table 4). Considering the three evaluation criteria together, it is quite difficult to judge about the performance of the SWAP model for individual layers. With respect to the IA values alone, however, the model gives consistently better results for the first layer. We attribute this to the fact that in deeper soil layers there are larger influences from lateral flows than in the top layer (while SWAP is a one-dimensional model that ignores any lateral flows). When considering soil layers separately, IA and NS are somewhat related, whereas RMSE is unrelated to the other two statistics. Strikingly, the relation between IA and NS becomes much stronger on averaging, whereas a relation with RMSE remains absent. We did not find any information about relations between performance indices under the influence of aggregation level in the hydrological literature.


View this table:
[in this window]
[in a new window]

  		TABLE 4. The values of several error criteria of soil moisture prediction (RMSE, index of agreement IA, and Nash–Sutcliffe NS, see Eq. [13–15]GoGo) per measurement location and layer, using the extrapolation calibration–validation scheme. Although our model distinguishes three layers, we distinguish four layers when considering the error criteria because there are separate observations for Layers 3 (40–60 cm) and 4 (60–80 cm).

 
Considering all three evaluation criteria, the quality of the model predictions would not be satisfactory when considering the separate layers; in 73% of the cases, either the IA or the NS coefficient is <0.50. When considering the averaged data, however, model predictions would be judged sufficient. Apparently the evaluation criteria are only meaningful at a specific scale.

When considering the RMSE alone, it seems to agree well with the results of previous studies. The SWAP study by Singh (2005) reported an RMSE that varied between 0.016 and 0.033 for different layers to 3 m down. Crescimanno and Garofalo (2005), using Mualem–van Genuchten equations (Mualem, 1976; van Genuchten, 1980), reported RMSE ranges of 0.037 to 0.101 and 0.035 to 0.078 at soil depths of 30 and 45 cm, respectively. Mertens et al. (2005), using soil moisture measurements at 25 locations at three different depths (at the surface and 30- and 60-cm depths) on an 80- by 20-m hillslope, reported RMSE values ranging from 0.041 to 0.089 when applying the MIKE-SHE model. Heathman et al. (2003), comparing the simulation results of the RZWQM model with soil moisture data from different depths up to 60 cm, found that RMSE ranged from 0.016 to 0.097. In a validation study of the THESEUS model in Germany, Wegehenkel (2005) obtained RMSE ranges of 0.030 to 0.043, 0.028 to 0.030, and 0.026 to 0.049 for Layers 1 (0–30 cm), 2 (30–60 cm), and 3 (60–90 cm), respectively. While the current study resulted in RMSE ranges of 0.015 to 0.045, 0.016 to 0.040, 0.012 to 0.050, and 0.011 to 0.051 for Layers 1 (0–20 cm), 2 (20–40 cm), 3 (40–60 cm), and 4 (60–80 cm), respectively.

The Relation between Level of Aggregation and Ability to Predict Soil Moisture
After evaluation of the model performance for each individual tube (AggL0), the evaluation was repeated at three different aggregation levels. In the first aggregation level (AggL1), the tubes located in each field were put in the same group. There were at maximum two tubes in each field. For the second aggregation level (AggL2), the tubes with the same land use were put in one group. We assumed three different land use types: cropland, grassland, and orchard. In the third aggregation level (AggL3), all tubes were put in one group. For every group, the observed soil moisture content of each layer is the average of the corresponding layer of the included tubes in the group. For the comparison across different aggregation levels, we used only RMSE as a performance measure (NS and IA showed behavior similar to RMSE for aggregation).

Like the optimization process for individual tubes, the initial values of the parameters per aggregation level were obtained from the literature according to soil type and land management practices. Since the tubes in each group of AggL1 have the same soil type and land management practices, the initial values of the parameters were in this case the same as for the model calibration with individual tubes. In AggL2, both the soil type and land management practices were different only for the tubes included in the cropland group. For this group, the calibration process started with the initial parameter values that had been used for the tube CST.u at AggL1 because its parameter values were intermediate to all other tubes in this group. Evapotranspiration differences between crop types did not play a role since the study period data covers the winter (transpiration is negligible and evaporation is the same for all tubes in this group).

In AggL3, the contemporary measured soil moisture content of all tubes for each layer were averaged and used as observation values. Initial values of the parameters were obtained by averaging the initial values of the parameters for each individual tube. For calculating the evapotranspiration rate in this aggregation level, grass was used as the representative land use type. The results of the different aggregation levels are presented in Tables 5 and 6. Table 5 shows the RMSE of the simulation results when they are compared with the aggregated soil moisture values for each group. Table 6 presents the RMSE values of the simulated moisture content for each group when compared with the measured ones of each individual tube in the related group. As expected, the RMSE values for each group (on the basis of the average of observed soil moisture for each group) are lower than the RMSE values for each individual tube (see also Fig. 4–7GoGoGo). The ranges of RMSE values for the three layers were 0.009 to 0.037, 0.006 to 0.025, and 0.009 to 0.025 m3 m–3 for AggL1, AggL2, and AggL3, respectively. These results are quite comparable with those found in other studies. For instance, Crow et al. (2005) found that the RMSE of spatially averaged predicted soil moisture content decreased from 0.043 to 0.032 cm3 cm–3 when compared with spatially averaged measured soil moisture contents at the field and regional scales, respectively (measuring and modeling soil moisture in the 0–6-cm surface layer). On the other hand, when comparing measured soil moisture of each individual tube with the simulated soil moisture of the pertinent group at different aggregation levels indicates that, by increasing the aggregation level, the RMSE increases, too (Table 6). The ranges of the RMSE for all layers are 0.011 to 0.051, 0.015 to 0.12, 0.018 to 0.174, and 0.015 to 0.181 cm3 cm–3 for AggL0, AggL1, AggL2, and AggL3, respectively. The degree of change in the RMSE is not equal for all soil layers. Layers 1 and 4 show relatively large changes in the RMSE values. Layer 2 shows the least changes. We think that the large changes in the first layer are related to land use type and land management practices, while the large changes in the fourth layer are related to the topography. To illustrate this last point, the RMSE values of the CST.d and CST.u tubes for AggL1 and AggL0 can be compared with the RMSE of the YM.d and YM.u tubes. Topography of CST.d and CST.u (which have the same land use and land management practices) differs a lot while for the YM.d and YM.u tubes it is comparable (see Table 1). Table 1 shows that the catchment area of CST.d is four times larger than that of CST.u. In this case, we think the higher values of soil moisture content in the deeper layers of CST.d are due to reinfiltration of surface runoff. The difference in wetness between CST.d and CST.u explains why the RMSE values of CST.d and CST.u (especially for Layers 3 and 4) are quite high for AggL1 in comparison to AggL0, while for YM.d and YM.u the RSME values are comparable. Beyond these qualitative considerations, we were unable to relate RMSE (or NS or IA) to its value at a lower (or higher) aggregation level and other explanatory variables with regression techniques. The explanatory variables used were those listed in Table 1 (land use, soil, slope steepness, upslope length, upslope area) and Table 2 ({theta}sat, {alpha}, n, {lambda}, and Ksat—all per layer, as well as averaged across the three soil layers—calibration {Phi}, calibration R, validation {Phi}, and validation R). Both generalized linear regression (McCullagh and Nelder, 1989) as well as generalized additive models (Hastie and Tibshirani, 1990) were used, while evaluating all combinations of the explanatory variables up to four dimensions.


View this table:
[in this window]
[in a new window]

  		TABLE 5. The RMSE of soil moisture predictions at different aggregation levels and per layer, using the extrapolation calibration–validation scheme. The measured values were averaged before calculating the RMSE.

 

View this table:
[in this window]
[in a new window]

  		TABLE 6. The RMSE of soil moisture predictions at different aggregation levels (AggL) and per layer, using the extrapolation calibration–validation scheme. To calculate RMSE, the point scale measurements were used.

 
The Relation between Level of Aggregation and Ability to Predict Discharge
After evaluation of the model performance at different aggregation levels with respect to soil moisture only, performance was also evaluated with respect to daily discharge. Figure 9 shows the comparison of observed vs. predicted discharge for the four aggregation levels.


Figure 9
View larger version (17K):
[in this window]
[in a new window]

  		FIG. 9. Discharge as predicted by the SWAP models for the different aggregation levels (AggL0, AggL1, AggL2, and AggL3) against the observed daily discharge. SWAP does slightly overpredict discharge for low observed discharge and underpredicts severely for high observed discharge.

 
Figure 9 shows that there is a strong effect of aggregation (lower aggregation levels tend to underpredict discharge more, while higher aggregation levels tend to overpredict, especially for lower observed discharge). The model also tends to overpredict the number of runoff events. To place the performance of the SWAP model in context, a few regression models, using only rainfall depth and event duration, were also fitted (using only the calibration data and evaluating for the validation data). Without exception, these had a much better performance (RMSE <2000). This comparison makes clear that SWAP, when calibrated with only soil moisture observations, does not provide informative discharge predictions.


   Discussion and Conclusions
TOP
 ABSTRACT
 INTRODUCTION
 Materials and Methods
 Results
 Discussion and Conclusions
REFERENCES
 
In the hydrologic literature, one-dimensional Richards' models are applied at all levels of aggregation that are encountered in this study, that is, point scale (Mertens et al., 2005), field scale (Feddes et al., 1978; Ritsema et al., 1996), and more inhomogeneous response units (Heathman et al., 2003; Crow et al., 2005). With these studies at various scales but a similar system conceptualization in mind, the question presents itself: Why is it (still) so hard to compare results between these studies at different aggregation levels? To contribute to answering this intriguing question, we have investigated model parameterization and performance for a research catchment at various levels of aggregation.

Despite the theoretical impossibility to derive effective parameters for the nonlinear Richards' equation by simple aggregation or disaggregation procedures—at least when considering fluxes (Kim and Stricker, 1996), many studies have ignored this fact and investigated the derivation of Mualem–van Genuchten parameters by both aggregation and disaggregation (see, e.g., Pachepsky and Rawls 1999; Zhu et al., 2006). These approaches have been stimulated by the availability of parametric and nonparametric pedotransfer functions, which in most cases do not explicitly refer to an aggregation level for their application (e.g., Wösten and van Genuchten, 1988; Pachepsky et al., 1996; Schaap et al., 1998).

In brief, our study shows that for our case study, large-scale measurements of soil moisture may be used to satisfactorily calibrate one-dimensional unsaturated model moisture predictions based on Richards' equation (as good as calibrating at a smaller scale). It also shows that the large-scale (calibrated) {alpha} and n parameters are statistically identical to those obtained based on averaging calibrated (or otherwise obtained) parameters representative of a much smaller (measurement) scale and averaged spatially. In addition, our study shows that the effect of a different calibration–validation scheme was small. The effect of using different error criteria was, however, significant (but unpredictable). This applies as long as only soil moisture distribution is predicted. This study tells nothing about the much more important ability to obtain upscaled fluxes from small-scale measurements of either soil moisture, fluxes, or model parameters (Ksat, {alpha}, and n). We illustrate this limitation by showing the limited ability of the SWAP model to predict catchment discharge (the only flux directly observed in our experiment). Also, this study does not tell anything about the relevance of the soil hydraulic parameters in relation to other important hydrological controls (such as terrain shape, land management, and vegetation), nor does it address the effect of selecting certain profile characteristics as valid for complete mapping units. A full sensitivity study would be useful for this.

The values for the Mualem–van Genuchten parameters {alpha} and n at higher aggregation levels are very close to the arithmetic mean of these parameters at lower aggregation levels (Ksat is omitted from this analysis because it was constrained at an early stage in the calibration procedure), see Fig. 8. These results are in line with the recommendations by Zhu and Mohanty (2002). The approach in our study also differs in a few ways from the aforementioned "effective parameter" studies. In the first place, our study lacks any geostatistical component, whereas most effective parameter studies consider spatially correlated fields of soil properties. Using this information about spatially heterogeneous soils, effective parameter studies normally focus on matching the steady-state vertical flow, whereas our study does not match an integrated flux but rather a set of distributed soil moisture observations. Finally, as stated above, effective parameter studies assume some kind of analytical or statistical relation between soil hydraulic parameters at different levels of aggregation. On the basis of this distribution, effective values for hydraulic parameters are directly derived at another aggregation level, i.e., without recalibrating a model. In our study, these hydraulic parameters are derived by calibration at each distinct aggregation level.

Evaluation of model performance at multiple scales is rare. Our study shows, however, that studies of this type are in fact required to enable the intercomparison between results from studies at different aggregation levels. Some frequently used performance measures (RMSE, the Nash–Sutcliffe coefficient, and the index of agreement) show a dependence on aggregation level, albeit a dependence that we were not able to fully explain or characterize in a simple manner. With higher aggregation levels, the prediction error becomes smaller, but not in an equal manner for the various indices. Hence, the relation between the various indices changes with aggregation level.

Usually, the validation scheme, as well as the parameter calibration procedure, is quite distinct between different modeling studies. This also complicates a comparison of results from different studies. To investigate the relative importance of these factors, we varied the validation scheme (i.e., interpolation vs. extrapolation) and the parameter optimization scheme (i.e., the procedure to select the subset of parameters to calibrate). The results led to the conclusion that both the validation and the parameter optimization scheme do not significantly influence the results of model parameter values and performance at different aggregation levels.

By calibrating the SWAP model on soil moisture data only, it is not possible to predict catchment discharge very well. The use of discharge data in addition to distributed soil moisture observations in calibration is an extension of the current study that deserves further investigation. In calibration, an additional integrative observation like discharge acts as a regularizing constant across aggregation levels. In the research catchment of this study, discharge was measured only at one location. Hence the discharge record to be matched is the same for each aggregation level, and parameter values would therefore become more similar between aggregation levels. A more extensive data set with discharge observations at various locations in the stream network could lead to a more complex behavior.


   ACKNOWLEDGMENTS
 
We acknowledge MSRT (Ministry of Science and Research and Technology of Iran) for financial support of this study. Jos van Dam is kindly acknowledged for helping with the SWAP-PEST interface. We are grateful to Piet Peters and Wim Spaan for their help in collecting the field data.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 Materials and Methods
 Results
 Discussion and Conclusions
 REFERENCES
 





This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Sheikh, V.
Right arrow Articles by van Loon, E. E.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Sheikh, V.
Right arrow Articles by van Loon, E. E.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Sheikh, V.
Right arrow Articles by van Loon, E. E.
Related Collections
Right arrow Watershed and Landscape Processes
Right arrow Effective Parameters
Right arrow Water Flow Models

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome