Published in Soil Sci. Soc. Am. J. 68:1526-1538 (2004).
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
DIVISION S-1—SOIL PHYSICS
Preferential Flow Revealed by Hydrologic Modeling Based on Predicted Hydraulic Properties
Frédéric Gérarda,*,
Mark Tinsleya,c and
K. Ulrich Mayerb
a INRA, Centre de Nancy, Unité Biogéochimie des Ecosystèmes Forestiers, Forêt d'Amance, 54280 Champenoux, France
b Univ. of British Colombia, Dep. of Earth and Ocean Sciences, 6339 Stores Road, V6T 1Z4 Vancouver, BC, Canada
c Present address: Chemistry School, Clark Hall, West Virginia Univ., Morgantown, WV 26506
* Corresponding author (gerard{at}nancy.inra.fr)
 |
ABSTRACT
|
|---|
Pedotransfer functions have shown a reasonable reliability and accuracy for predicting soil hydraulic properties. However, the contribution of the range of macroscopic features leading to preferential water flow is not readily taken into account. We modified and used the hydrological component of the reactive transport model MIN3P and the neural network-based code ROSETTA in an attempt to simulate 4 yr of daily measurements of the soil water content in a forest soil covered by Douglas-fir (Pseudotsuga menziessii Franco). A good fit of the mean measured water contents was obtained during periods of low soil moisture, while the model tended to overpredict water contents during periods of high soil moisture. This behavior is typical for the presence of significant preferential flow. Slightly better results were obtained by using predicted values of the saturated hydraulic conductivity, while the assumption of a water table located at shallow depth increased discrepancies. A good match was obtained by calibration of a simple preferential flow scheme, which was based on the assumption that the retention properties of the porous network control preferential flow. Accordingly, preferential flow seemed to initiate within the capillary pore domain. This causes a much greater sensitivity of the results to the position of the water table than with other schemes that consider pure gravity-driven flow in large macropores. Knowledge of the functional pore size is needed to ascertain the type of preferential flow scheme to be used.
Abbreviations: FV, functional variable • HC, hydraulic conductivity • MRE, mean residual error • PTF, pedotransfer function • REW, reserve of extractable water • RMSE, root mean square error • TDR, time-domain reflectometry • WRC, water retention curve
|
INTRODUCTION
|
|---|
KNOWLEDGE OF SOIL HYDRAULIC PROPERTIES is crucial for the modeling of flow in unsaturated soils. The accurate estimation of the water flux is of importance for environmental issues such as metal mobility and nutrient cycling. Two of the key properties of an unsaturated soil are the soil water retention curve (WRC) and the hydraulic conductivity (HC). Measurements of the soil HC and/or WRC are costly, time-consuming, and sometimes unreliable because of soil heterogeneity and experimental errors. Methods have been developed to estimate soil hydraulic properties from more easily measured soil properties. These methods involve the use of pedotransfer functions (PTFs). A number of PTFs can be found in the literature and can be classified according to the nature of the basic input soil properties and the method adopted to generate predictions (Wösten et al., 2001). Recent publications focus on comparing PTF predictions with independent data sets of hydraulic properties measured in the laboratory. Results of the comparisons vary; good agreement was found in several studies (Schaap and Leij, 2000; Cornelius et al., 2001; Rawls et al., 2001; Wagner et al., 2001) whereas significant discrepancies were observed by others (Chen and Payne, 2001; Pachepsky and Rawls, 2003; Soet and Stricker, 2003). The evaluation of the PTFs performance in various field situations has also been addressed in several publications (e.g., see review by Wösten et al., 2001). A set of on-site measured pressure heads and/or volumetric moisture contents at different soil depths are generally considered as functional variables (FV) to evaluate the predictive performance of the PTFs (e.g., Espino et al., 1996; Hack-ten Broeke and Hegmans, 1996; van Alphen et al., 2001; Nemes et al., 2003). Actual hydraulic properties were measured and served to pre-evaluate the predictive accuracy of PTFs (Nemes et al., 2003; Soet and Stricker, 2003) and, more frequently, to run numerical simulations in addition to those parameterized with PTF predictions (Wösten et al., 1990; Hack-ten Broeke and Hegmans, 1996; van Alphen et al., 2001). The accuracy of the predictions can vary appreciably according to the PTF used and the number of input data considered (e.g., van Alphen et al., 2001; Wösten et al., 2001; Nemes et al., 2003). Overall, however, the performance of PTFs for prediction can be described as reasonably good given the inherent complexity of the hydraulic behavior of soils at the field scale.
We would like to emphasize two general limitations in the past studies targeting the functional evaluation of PTFs. First, the widespread occurrence of preferential flow in soils has never been considered. Furthermore, functional evaluation has not been performed in a forest soil and there is ample evidence for fast preferential flow of water and solutes in such soils (e.g., Feyen et al. [1999] and references therein). A number of macroscopic features cause preferential flow to occur in soils in general, and in forest soils in particular. These include cracks produced by shrinkage of clay soils and biotic pores formed by the decay of dead plant roots and fauna (e.g., Larsson et al. [1999] and references therein). Second, in previous functional evaluation studies, the number of FV data points used was rather small compared with the duration of the simulated time series (typically 1–2 yr). Therefore, knowing that the accurate prediction of high water contents is difficult due to the presence of preferential flow channels (e.g., Simunek et al., 2003), we believe that the relative lack of FV data corresponding to these periods could preclude the detection of preferential flow while testing the performance of PTFs at the field scale. Furthermore, provisions are generally taken to avoid macroscopic heterogeneities in soil samples studied in the laboratory. In other words, samples often better reflect the hydraulic properties of the soil matrix rather than field-scale values that include a range of structural features potentially causing preferential flow to occur. This fundamental bias, often referred to as scaling effect, is well recognized in studies involving the use of PTFs at various spatial scales (e.g., Inskeep et al., 1996; Chen and Payne, 2001; Wagner et al., 2001). Researchers are attempting to overcome this issue (Lin et al., 1999; Jarvis et al., 2002; Pachepsky and Rawls, 2003). The PTFs proposed by Lin et al. (1999) gave reasonable predictions of some of the hydraulic properties controlling preferential flow, especially the HC at water saturation, but these PTFs also require a rather large number of predictors. Moreover, some of these, such as the macropore porosity, can be difficult to measure in the field. Jarvis et al. (2002) tentatively established a PTF for the prediction of the near-saturated soil HC involving a fewer number of predictors, but they obtained a poor match with measurements.
Concerns regarding forest decline and sustainable land management have motivated research on ecosystems evolution since the 1970s. Nowadays networks of forested sites are subjected to various degrees of monitoring of soil physical and chemical properties and many other ecological variables. One of these, the Vauxrenard field site (Rhône, France), was managed by INRA, FMN department (Institut National de la Recherche Agronomique, Forêts et Milieux Naturels), from April 1992 to January 2000. Among other instruments progressively implemented on-site, replicated time-domain reflectometry (TDR) probes have provided 4-yr records of hourly measurements of the soil volumetric water content at different soil depths. Clearly, this survey has led to an extraordinary database of FV data, which is thus particularly well suitable to a functional evaluation of PTFs in a forest soil with a special focus on preferential flow. The objectives of the present paper are two-fold: (i) to perform a functional evaluation of a selected PTF for a forest soil; (ii) to study whether significant preferential flow is occurring in this forested soil, and to determine the prediction errors made in terms of soil moisture when preferential flow is ignored.
|
MATERIALS AND METHODS
|
|---|
Study Field Site
The field site was situated on a gentle slope (<10°) nearby the summit of a small mountain (Montagne des Aiguillettes, Rhône, France) within the Beaujolais Hills. Site elevation is about 60 m at the foot of the hill and reaching an 842-m altitude at the top. A mean annual temperature of 7°C and a mean annual precipitation of roughly 1000 mm characterize the climate. The soil is of brown acidic type, classified as Typic Dystrochrept according to USDA (Soil Survey Staff, 1994), or as an Alocrisol (AFES, 1992). A 40-yr old douglas-fir forest covers the field site. The parent material from which the soil is derived consists of a volcanic tuff dating from the upper Visean. Regional metamorphism has led to its intensive recrystallization (Ezzaïm et al., 1999). Many other aspects of the site have been discussed in several studies, such as stand characteristics, nutrient dynamics, and silicate weathering mechanisms (e.g., Ranger et al., 1995; Marques and Ranger, 1997; Ranger et al., 2001; Gérard et al., 2003). Soil texture is loamy, organic matter is rather high (6–7% down to 20 cm), and most roots are located in the top 50 cm but some have been observed down to a depth of 100 cm (Ranger et al., 2001). Root diameter distribution according to soil depth is given in Table 1. Four soil horizons can be distinguished in a representative soil profile: the A11 (0- to 12-cm depth), the A12 (12- to 35-cm depth), the B (35- to 55-cm depth), and the C (55- to 120-cm depth).
Functional Variable Dataset and Evaluation Criteria
Soil water content were measured hourly at three soil depths (15, 30, and 60 cm) during 4 yr by means of TDR probes (Trase System, Soil Moisture Equipment Core, Goleta, CA). Six replicates per depth were installed to evaluate the influence of spatial heterogeneity. Data is available from January 1996 to December 1999, with some time gaps caused by the temporary breakdown of the monitoring system. Time domain reflectometry probes were internally calibrated and tested in free water and in dry sand for the maximum and minimum moisture values before installation. Onsite control tests have also been completed in June 1999 (Fig. 1) . Control values of the soil water content were measured nearby (approximately 0.5 m) in soil samples (gravimetric measurements) and in situ by means of a portable TDR probe (Delta-T Devices, Theta Probe, Thermo Electron Corp., Waltham, MA) installed in holes excavated with an Edelman auger. The auger diameter was marginally smaller than the probe to ensure good contact between the probe and the soil material. A statistical test (ANOVA, Unistat 5.0 software, Unistat Ltd., England) showed a lack of significance differences in the measured soil water contents.
View larger version (25K):
[in this window]
[in a new window]
|
Fig. 1. Comparison of on-site control measurements of the soil water content at different soil depths. Error bars correspond to the standard deviation. White bars: mean value measured with the trase-system Soil Moisure time domain reflectometry (TDR) probes; light-gray colored bars: mean value measured with the Theta Probe Delta-T Devices; dark-gray colored bars: gravimetric measurements.
|
|
Two statistical criteria were used for model evaluation purposes. The first is the root mean square error (RMSE):
 | [1] |
where Pi is the model prediction, Oi the observed value, and N is the number of observation.
This RMSE is a measure for the scatter around the observed mean water content between measured and simulated values. The mean residual error (MRE) was also used to evaluate the direction of the error:
 | [2] |
Values of the MRE close to zero indicate that measured and simulated moisture contents do not differ systematically from each other.
Model and Parameters
The hydrological component of the reactive-transport model MIN3P was used for the simulations. The MIN3P model is a coupled multicomponent reactive transport model (e.g., Mayer, 1999; Mayer et al., 2002). Richards equation is solved for pressure head using Eq. [3] to describe the relationship between pressure head and soil water content.
 | [3] |
where h is the pressure head, n and m are experimental exponents, S is the saturation given by the ratio 
/s, and Sr is the residual saturation given by the ratio r/s, where , r, and s are the volumetric water content, the residual water content, and the saturated water content, respectively. The van Genuchten–Mualem relation (Eq. [4]) is used for unsaturated hydraulic conductivity (van Genuchten, 1980).
 | [4] |
where krel is the relative hydraulic conductivity, L is an empirical pore tortuosity/connectivity parameter, and Se is the effective saturation given by:
 | [5] |
The unsaturated HC is calculated using:
 | [6] |
where K0 (m s–1) corresponds to the matching point at saturation. The empirical exponent L is commonly assumed to be 0.5, as recommended by Mualem (1976). However, different values (mostly negative ones) have also been widely used to improve the fit to experimental data under unsaturated conditions (e.g., Wösten et al., 1999; Schaap et al., 2001). This often corresponds to values of K0 one order of magnitude lower than the measured hydraulic conductivity at saturation, Ks (Schaap and Leij, 2000).
In MIN3P, the solution of Richards equation is obtained using integrated finite differences in space and implicit time weighting. The equations are linearized using a modified Newton's method and solved using a sparse iterative matrix solver (Vanderkwaak et al., 1997). For the present application the model has been extended in two ways; first, to account for the influence of evapotranspiration on the soil water content, and second to introduce functions that allow the simulation of preferential flow. Several numerical verification tests were conducted to ensure proper behavior of the model including single soil layer simulations and comparison with single permeability simulations. Also trends in soil water content were examined as various parameters of the system were varied to ensure that the model was behaving qualitatively reasonably.
The rate of water uptake and loss is calculated at each time step by adopting a method similar to the one used in the soil water module SWIF of the FORHYD model (Tiktak and Bouten, 1992; Tiktak and Bouten, 1994; Bouten, 1995). The mass balance equation for water in the forest ecosystem can be expressed as:
 | [7] |
where PET is the potential evapotranspiration (mm d–1), T is the actual tree transpiration rate (mm d–1), Eu corresponds to the sum of understory plus soil evaporation (mm d–1), and I is the amount of intercepted water by the tree canopy (mm). Intercepted water is evaporated at a certain rate, hence the fi factor is introduced (d–1). The correction of T by the fiI term serves to take account of the tree transpiration reduction by the evaporative demand of a wet canopy on days with precipitation.
Daily PET values (Penman), rainfall (P), and throughfall (Th) used in the simulations can be seen in Fig. 2
. Potential evapotranspiration and P were obtained from a nearby meteorological station (Météo France Co., Tarrare Station, Rhône). The daily amounts of throughfall were calculated with the process-oriented regression model of Vilette (1994), where parameters were obtained from the best fit between monthly throughfall measured onsite and monthly rainfall measured in a field located nearby. Consistent with the work of Vilette (1994), Eu was set to 8% of PET. This proportion was determined from the mean ratio between solar radiation above and below the canopy measured on-site. In agreement with most of the previous studies on water balance modeling in a forest soil, water uptake was arbitrarily limited to the topsoil layer. The generic value of fi = 0.2 from Granier et al. (1999) was used. Tiktak and Bouten (1994) used a similar value (fi = 0.141) for the simulation of the long-term water balance in a douglas-fir stand in the Netherlands.
View larger version (43K):
[in this window]
[in a new window]
|
Fig. 2. Variations in precipitation (P, dark solid line), throughfall (Th, square symbols), and potential evapotranspiration (PET, light solid line).
|
|
The actual tree transpiration is distributed in each ith discretization control volume according to the root length density of the jth soil layer and weighted by the local soil moisture content:
 | [8] |
where Qi,j (m3 s–1) is the water uptake rate by trees in the ith control volume of the jth soil layer, Vi,j is the volume of the control volume, Rj (m2 m–3) is root length density in the jth soil layer (i.e., the ratio of the root surface per volume of soil), Si,j is the water saturation in the ith control volume of the jth soil layer.
The decrease in the water uptake at low water content is also taken into account in the model. However, the use of an implicit numerical scheme in the model for solving mass balance equations did not allow the implementation of the common discontinuous reduction function. Instead, the continuous function used by Battaglia and Sands (1997) was implemented into MIN3P:
 | [9] |
where 
i,j is the dimensionless reduction factor, REWi,j is the reserve of extractable water, REW0 and p1 are fitting parameters. REWi,j is defined as:
 | [10] |
where Slimj is the water saturation at wilting point (i.e., no water uptake below this critical value), and Sfj is the water saturation at field capacity.
Measurements of the ratio between T and PET as a function of REW given by Granier et al. (1999) was fitted with Eq. [10] using the nonlinear fitting procedure of the Unistat 5.0 software (Unistat Ltd, England). We considered data corresponding to a leaf area index equal to 6 m2 m–2, which is close to the generic value proposed for coniferous trees in a recent review article (Breuer et al., 2003). The same leaf area index was used by Vilette (1994).
Values of the soil-related transpiration parameters; that is, R, the root length density, and REW, the reserve of extractable water, were obtained during calibration of the model (e.g., Musters and Bouten, 2000; Vrugt et al., 2001; Hupet et al., 2003). It is necessary to stress that parameter optimization is often required for practical reasons during model calibration. While the maximum rooting depth was fixed to 1.2 m based on measurements (see Table 1); in the present study R could not be easily measured in a 45-yr-old douglas-fir stand (i.e., trees of about 20 m high). Similar considerations can be made with respect to the definition of the reserve of extractable water, REW, except that the likely range of values for the parameters in Eq. [10] is better bounded. Typically, the water saturation at wilting point, Slim, is assumed as log|h(cm)|or pF = 4.2 in agronomy. Slightly greater values have been frequently reported depending on a variety of factors such as soil texture and plant species (e.g., Volaire and Lelièvre, 2001). The saturation at field capacity, Sf, defines the upper limit of available water for roots. This maximum again may vary according to a number of factors including plant species. In theory, Sf can range from close to saturation to some value greater than Slim.
Numerous techniques have been developed for the implementation of preferential flow (see Simunek et al. [2003] for a review). Briefly these techniques may be subdivided into two types; equilibrium and nonequilibrium formulations (i.e., instantaneous and slow transfer between macroporosity and the matrix, respectively). Both types of preferential flow schemes have been implemented in MIN3P. The equilibrium formulation involves a composite HC function to allow for divergences from the standard unsaturated HC behavior at higher water content (Mohanty et al., 1997). The calculated HC is an average property of the two pore systems. Below a certain pressure head, noted hpf, the standard van Genuchten–Mualem relationship is used to calculate the HC. Above that head the HC is assumed to increase exponentially according to:
 | [11] |
where 
is a scaling factor that determines the magnitude of preferential flow.
The nonequilibrium preferential flow scheme implemented in MIN3P is based on the work of Gerke and van Genuchten (1993a)(1993b). The water flux evolution is calculated by solving Richards equation in both the macropore and micropore system (using Eq. [3] and [4]), and a first-order process-based expression is used to calculate interstitial flow between micro and macropore regions. However, such a scheme requires an increased number of input parameters for simulating preferential flow. Since the present application would require calibration of poorly constrained preferential flow parameters to improve the fit between simulated and measured water contents, only the equilibrium formulation is considered here.
The PTF available in the ROSETTA software (Schaap et al., 2001) was used in this study to predict the soil hydraulic parameters. This PTF is based on an artificial neural network and constitutes one of the most recent PTFs that overall has shown reasonable predictions evaluation studies. For example, Nemes et al. (2003) recently showed that the functional performance of the ROSETTA model with different sets of input data was reasonably good for simulating soil moisture variations in the field. An advantage of neural network approaches resides in the fact that no a priori model concept is required to convert basic soil properties to hydraulic parameters. The ROSETTA model is flexible as its hierarchical structure enables the input of limited and more extended sets of predictors. A trend of improvement was observed with an increasing number of predictors (Nemes et al., 2003). This behavior has been frequently found in other evaluation studies of PTFs (Rawls et al., 2001; van Alphen et al., 2001; Wösten et al., 2001). It is unfortunate that soil hydraulic properties were never measured at the Vauxrenard field site, and are not measurable today because the storm on 27 Dec. 1999 destroyed the 45-yr-old douglas-fir stand (like millions of hectares throughout Western Europe) and the soil was deeply tilled by falling trees. Therefore, contrary to past studies only a blind-functional evaluation could be undertaken; that is, without any possible reference to measured hydraulic properties. However, this presented an opportunity to evaluate the reliability of a PTF already recognized for its predicting capacity in the context of great scarcity of information on soil properties. Available predictors for the present application are sand, silt and clay content, and bulk density data for each soil layer. Basic soil properties had been measured for various purposes in each soil layer at a number of locations randomly distributed throughout the site. For simplicity, we used the mean values of the basic soil layer properties (Table 2) as input for ROSETTA. Another interesting aspect of the ROSETTA model lies in the use of the bootstrap method to provide uncertainty estimates of the predicted hydraulic parameters, which can then be considered in the simulations. Also, the common relationship m = 1 – 1/n is used in ROSETTA and predicted values of both Ks (with L = 0.5) and K0 (with L 
0.5) are generated. Uncertainties exist in the reliability of either K0 or Ks values to model field water content. Besides the order of magnitude of difference typically encountered between these two soil parameters, Schaap et al. (2001) also outlined that using K0 < Ks leads to an untenable situation near saturation because the HC should be equal to Ks while Eq. [4] provides K0.
In addition to the lack of measured hydraulic parameters, the position of the water table was not determined. However, the local situation of the Vauxrenard field site (i.e., hill slope, crystalline bed rock presumably fractured, moderate soil thickness) suggests that if a perched aquifer were present or not beneath the field site it should be very deep. Moreover, the soil is recognized to be well aerated and hydromorphic features have never been observed in this soil to a depth of 2.5 m (i.e., the greatest depth reached by trenches excavated for various purposes). Accordingly, the total Fe concentration in soil solutions collected down to 1.20 m was always found to be very low, on the order of a few micromoles, as controlled by formation of low solubility Fe(III) hydroxides (Bourrié et al., 1999).
Modeling Procedure
All simulations are performed for a one-dimensional soil profile. The top 120 cm of soil are subdivided into 120 control volumes of equal size. Each soil horizon is defined in the simulation by a property layer having its own set of hydraulic and soil-dependent transpiration parameters with the upper boundary open to throughfall. In a first step, we use the outputs of ROSETTA to simulate daily variations in soil moisture during the 4-yr period of TDR measurements. The reliability of the two types of predictions for the K0 and L parameters was studied. Values of the soil-related transpiration parameters (i.e., R and REW) were manually optimized in each soil layer to obtain the best fit during the drought events of the January 1996 and December 1997 time series; that is, in the summer when their influence is greater due to the low water content and the maximum transpiration demand. Then, the calibrated values were kept constant for the simulation of the following 2-yr period of TDR measurements, that is, from January 1998 to December 1999. The influence of the position of the water table on the simulations was minimized at this stage by setting it to 100 m below the bottom of the soil profile (i.e., lower boundary condition h = –100 m). In a second step, the influence of prediction uncertainties and of the position of the water table is examined using the best simulation found previously. Lastly, the equilibrium preferential flow scheme and manual parameter estimation were used to improve the fit. In all the simulations, the steady-state solution was calculated to give the initial conditions in the soil column, with the mean value of PET, I, and Th been used to generate the steady- state solutions. Also, each simulation was started three months before January 1996. This initial period was ignored in data treatments to suppress any influence of the different initial conditions.
|
RESULTS AND DISCUSSION
|
|---|
Hydraulic parameters and uncertainties predicted with ROSETTA are presented in Table 3. As reported by Schaap and Leij (2000), the difference between K0 and Ks reached about one order of magnitude. Results of the simulations made by using the values of K0 (with L 0.5) and Ks (with L = 0.5) are illustrated in Fig. 3
. For comparison purposes, the same values of the soil-dependent evapotranspiration parameters were used in both cases. In each soil layer Sf was set to saturation and Slim corresponds to |h| = 160 m (pF = 4.2), which corresponds to the theoretical wilting point (see above). Results were found to be very similar. Whether Ks or K0 were used in MIN3P, results clearly exhibited an overestimation of the soil water content during periods of high soil moisture. Simulated water contents even exceeded the measured maximum as given by the mean measured soil water content plus the standard deviation. Also, calculated water contents were slightly less overestimated within the high soil moisture periods when using Ks (i.e., smaller MRE and RMSE—see Table 4). The value obtained from these statistic criteria decreased from 0.052 at the 15-cm depth to 0.028 at the 60-cm depth when K0 was considered, in comparison with 0.036 and 0.027, respectively, when using Ks. Consistently, the calculated RMSE follow the same trends with soil depth. By distinguishing the seasons on a calendar-basis, the MRE reaches its maximum in the winter and its minimum during the spring or in the summer. The maximum MRE value was 0.089 at the 15-cm depth with K0 in comparison to 0.066 when using Ks. In each season a tendency to decrease was found for the values of both statistical criteria with increasing soil depth. These results, and particularly the overpredictions of the water content during periods of high soil moisture and the such systematic decrease in discrepancy with increasing soil depth, strongly suggest the occurrence of preferential flow, which appears to be related to the root density distribution (see Table 1). Better results obtained by using Ks instead of K0, and the corresponding values of L, were likely caused by the much smaller value determined for K when K0 was considered and the soil water content was high. The same conclusions can be drawn to a lesser extent with respect to results obtained for the other seasons, notably in the summer and spring when the model becomes very sensitive toward the values of the evapotranspiration parameters. Readjusting the values taken by REW and R did not bring significant improvements. Improving the fit of the mean measured water content locally caused the predictions to become poorer elsewhere. Therefore, it appears that for the present application the Ks approach provides better results, and will be considered hereafter.
View larger version (53K):
[in this window]
[in a new window]
|
Fig. 3. Result of the simulations made with predicted hydraulic parameters and using the two types of prediction for hydraulic conductivity: Ks with L = 0.5 (cross), and K0 with L 0.5 (circles). Mean and spread of the soil moisture content as measured by means of six time domain reflectometry probe replicates are shown (thick line: mean; thin lines: mean ± sd).
|
|
View this table:
[in this window]
[in a new window]
|
Table 4. Statistical criteria at the 15-, 30-, and 60-cm depth for simulations using the two types of prediction for the hydraulic conductivity (Ks with L = 0.5, and K0 with L 0.5). Statistics are provided for the total simulation period (from January 1996 to December 1999) and for each season.
|
|
Predictions showed little variations in the simulated water contents within the range of uncertainty given by ROSETTA. An example at the 60-cm depth is shown in Fig. 4
, where the effect of uncertainties was large enough to warrant graphical presentation. The range of the uncertainty effects on the simulated water contents was established by testing all the combinations of the van Genuchten–Mualem parameters plus or minus uncertainties.
View larger version (35K):
[in this window]
[in a new window]
|
Fig. 4. Simulated water content at 60 cm by taking into account the prediction uncertainties given by ROSETTA (light curves) and using Ks and L = 0.5. Cross: simulation with the mean values of hydraulic parameters. Mean and spread of the measured soil moisture content are shown (thick line: mean; thin lines: mean ± sd).
|
|
The mean values of the predicted hydraulic parameters were used for testing the influence of a hypothetical water table. Given the uncertainty in the location of the water table and the potential for water table fluctuations (see above), several test simulations were conducted by setting a fixed water table at various depths during the wet periods (winter and fall), while a deep water table (i.e., water table set at the 100-m depth) was still assumed during the dry periods (spring and summer). Results showed that the impact of the water table becomes significant at the 60-cm depth when specified at a depth of about 3.5 m or shallower. As shown in Fig. 5
, discrepancies during the periods of high soil moisture become larger when the water table is closer to the surface while no significant change was observed during the periods of low water content (i.e., dry periods). Additional tests were performed, in which case the water table was kept at the same level during the entire course of the simulation. Under such conditions, simulated water contents markedly increased during dry periods as well, but could be compensated by readjusting the soil-dependent transpiration parameters. In the example shown in Fig. 6
, the water table was set to a 3-m depth, and a reasonable fit was obtained by setting Sf to the value given at |h| = 10 m (Slim unchanged, given by |h| = 160 m) in each soil layer. Values of the root length density needed to be readjusted as well. The overprediction during the periods of high soil moisture were even greater than previously calculated by minimizing the influence of the water table (i.e., position set at the 100-m depth).
View larger version (38K):
[in this window]
[in a new window]
|
Fig. 5. Simulated water content at 60 cm using Ks and L = 0.5 and by setting the water table during the wet periods (i.e., winter and fall) at the 2-m depth (open squares) and at the 3-m depth (filled squares). For comparison purposes, results calculated by neglecting the influence of the water table are given (cross). Mean and spread of the measured soil moisture content are shown (thick line: mean; thin lines: mean ± sd).
|
|
View larger version (35K):
[in this window]
[in a new window]
|
Fig. 6. Simulated water content at 60 cm using Ks and L = 0.5 and by setting the water table at the 3-m depth, without readjusting the soil-dependent transpiration parameters (filled squares), after readjustment (open diamonds). For comparison purpose, results calculated by neglecting the influence of the water table are given (cross). Mean and spread of the measured soil moisture content are shown (thick line: mean; thin lines: mean ± sd).
|
|
The preferential flow scheme proposed by Mohanty et al. (1997) was used in an attempt to correct for the overpredictions of calculated soil water content. As expected, accounting for preferential flow drastically improved predictions (Fig. 7)
. Statistics presented in Table 5 confirm this result. Corresponding fitted values of parameters in Eq. [11] are = 100 m d–1 in each horizon, and |hpf| = 8.3 and 4 m in the top three horizons and in the bottom soil horizon, respectively. The value of the HC associated with preferential flow, , is about two orders of magnitude larger than Ks predicted by ROSETTA. Such a difference was not surprising considering that, for example, Feyen et al. (1999) reported a decrease of travel time for a tracer of about three orders of magnitude in a forest soil in comparison to travel times estimated from saturated HC of the soil matrix. By application of the Jurin's law, values taken by hpf yield an equivalent pore diameter of about 4 µm at the 15- and 30-cm depths, and of 8 µm at the 60-cm depth, which are within the range of the capillary soil porous system. Such a low pore diameter limit can be viewed as evidence for the occurrence of a bimodal porosity distribution in this soil (Othmer et al., 1991), which was not considered by the ROSETTA PTF. However, such pore diameter values are theoretical since the complexity of the porous network is not considered (i.e., pore roughness, tortuosity, heterogeneity of the surrounding solids). This implies that the effective pore diameter corresponding hpf should be larger. Furthermore, we believe that the above limits are uncertain in the absence of physical measurements. These values may be caused by the type of preferential flow scheme used since the WRC predicted by the PTF is being considered together with the extra term accounting for preferential flow (see Eq. [11]). Any of the other modeling approaches for preferential flow also attempt to correct overestimations during periods of high soil moisture through parameter optimization. For example, the approach available in the MACRO model (e.g., Roulier and Jarvis, 2003) which considers preferential flow driven by gravity only; that is, occurring within macropores large enough to neglect the influence of capillary forces. Knowledge of the active pore diameters wherein preferential flow occurs could guide modelers in the selection of the preferential flow schemes (i.e., only gravity-driven flow in large macropores, gravity and capillary-driven flow in case of smaller macropores, or the need for the co-consideration of both formulations). In practice, comprehensive measurements of the soil hydraulic properties (i.e., WRC and HC function) especially at large water contents and at the field-scale should be considered.
View larger version (53K):
[in this window]
[in a new window]
|
Fig. 7. Simulated water content with preferential flow (circles), after calibration of the parameters involving in Eq. [11]. For comparison purposes, results calculated without preferential flow and using Ks and L = 0.5 are shown (cross). Mean and spread of the measured soil moisture content are shown (thick line: mean; thin lines: mean ± sd).
|
|
View this table:
[in this window]
[in a new window]
|
Table 5. Statistical criteria at the 15-, 30-, and 60-cm depth for the simulations performed including preferential flow. Statistics are provided for the total simulation period (from January 1996 to December 1999) and for each season.
|
|
The type of preferential flow scheme is also of relevance considering that the preferential flow schemes used here led to a much greater sensitivity of the results with respect to the position of the water table. The influence of the water table was minimized in the above calculations by setting it at the 100-m depth. By assuming a water table closer to the bottom of the simulated soil profile, capillary- driven upward preferential flow led to much larger discrepancies. This observation is illustrated in Fig. 8
, where the water table was set to a depth of 3 m during winter and fall, similar to the simulations shown in Fig. 5. It can be seen that the upward water movements became so significant that simulated water contents mirrored the stepwise adjustment of the water table. Clearly, such a large sensitivity to the water table level would not be obtained when using a gravity-driven preferential flow scheme. Future development of a more comprehensive scheme in MIN3P enabling the simulation of preferential flow in both macropores (i.e., gravity-driven flow) and in smaller pores (capillary-driven flow) would also be of advantage; as such processes are very likely to be present in a natural system.
View larger version (45K):
[in this window]
[in a new window]
|
Fig. 8. Influence of the position of the water table on the results of simulations performed with preferential flow at the 15-cm depth (Gray circles: water table set at the 3-m depth in winter and fall, and at the 100-m depth in summer and spring). For comparison purposes, results calculated by setting no influence of the water table are shown (open circles).
|
|
|
CONCLUSIONS
|
|---|
We presented an investigation of evapotranspiration and preferential water flow in the vadose zone of a forest soil using an extended version of the reactive-transport model MIN3P. With the exception for periods of high soil moisture, during which significant and systematic overpredictions were observed, the neural network-based PTF implemented in the ROSETTA model showed good predictive ability for simulating daily average values of the measured water content over a 4-yr period. Various sensitivity tests were performed and showed their inability to correct for such overestimates. As expected from these results and the literature, a good fit was obtained by considering the occurrence of preferential flow in this soil. Thus, it appears the ROSETTA PTF can give reasonable predictions of the hydraulic properties corresponding to the soil matrix of the present study field site. We further demonstrated the limited physical meaning of the preferential flow parameters obtained by calibration. These parameters are model-dependent in the absence of measured soil hydraulic properties at high water content. Process understanding can be enhanced by studying the corresponding size and nature of the active pores, as well as by means of comprehensive on-site measurements of the soil hydraulic properties. Results of such measurements should also help to determine the type of preferential flow formulation to be used in a model. Such findings should be of further practical relevance, because the use of preferential flow schemes that include capillary-controlled flow may lead to massive upward capillary-driven preferential flow.
|
ACKNOWLEDGMENTS
|
|---|
The authors gratefully acknowledge financial support for this cooperative research from the French Canada Research Foundation, the French Embassy in Ottawa, the INRA FMN department (Institut National de la Recherche Agronomique, Forêts et Milieux Naturels), and NSERC (Natural Sciences and Engineering Research Council of Canada). Many thanks also go to the GIP-ECOFOR Program for funding TDR instrumentation and measurements. The authors are grateful to D. Gelhaye for supervising the TDR system, and to the scientific manager of the field site, J. Ranger.
Received for publication August 20, 2003.
|
REFERENCES
|
|---|
- AFES. 1992. Principaux sols d'Europe. Référentiel pédologique Association Française pour l'Etude des Sols. Institut National de Recherches Agronomique, Paris.
- Battaglia, M., and P. Sands. 1997. Modelling site productivity of Eucalyptus Globulus in response to climatic and site factors. Aust. J. Plant Physiol. 24:831–850.
- Bourrié, G., F. Trolard, J.-M.R. Génin, A. Jaffrezic, V. Maître, and M. Abdelmoula. 1999. Iron control by equilibria between hydroxy-green rusts and solutions in hydromorphic soils. Geochim. Cosmochim. Acta 63:3417–3427.
- Bouten, W. 1995. Soil water dynamics of the Solling spruce stand, calculated with the FORHYD simulation package. Ecol. Modell. 83:67–75.
- Breuer, L., K. Eckhardt, and H.-G. Frede. 2003. Plant parameter values for models in temperature climates. Ecol. Modell. 169:237–293.
- Chen, C., and W.A. Payne. 2001. Measured and modeled unsaturated hydraulic conductivity of a Walla Walla silt loam. Soil Sci. Soc. Am. J. 65:1385–1391.[Abstract/Free Full Text]
- Cornelius, W.M., J. Ronsyn, M. Van Meirvenne, and R. Hartmann. 2001. Evaluation of Pedotransfer functions for predicting soil moisture retention curve. Soil Sci. Soc. Am. J. 65:638–648.[Abstract/Free Full Text]
- Espino, A., D. Mallants, M. Vanclooster, and J. Feyen. 1996. Cautionary notes on the use of pedotransfer functions for estimating soil hydraulic properties. Agric. Water Manage. 29:235–253.
- Ezzaïm, A., M.-P. Turpault, and J. Ranger. 1999. Quantification of weathering processes in an acid brown soil developed from tuff (Beaujolais, France). Part I. Formation of weathered rind. Geoderma 87:137–154.
- Feyen, H., H. Wunderli, H. Wydler, and A. Papritz. 1999. A tracer experiment to study flow paths of water in a forest soil. J. Hydrol. (Amsterdam) 225:155–167.
- Gérard, F., J. Ranger, C. Ménétrier, and P. Bonnaud. 2003. Silicate weathering mechanisms determined using soil solutions held at high matric potential. Chem. Geol. 202:443–460.
- Gerke, H.H., and M.Th. van Genuchten. 1993a. Evaluation of a first order water transfer term for variably saturated dual-porosity models. Water Resour. Res. 29:1225–1238.
- Gerke, H.H., and M.Th. van Genuchten. 1993b. A dual porosity model for simulating the preferential movement of water and solutes in structure porous media. Water Resour. Res. 29:305–319.
- Granier, A., N. Bréda, P. Biron, and S. Vilette. 1999. A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands. Ecol. Modell. 116:269–283.
- Hack-ten Broeke, M.J.D., and J.H.B.M. Hegmans. 1996. Use of soil physical characteristics from laboratory measurements or standard series for modelling unsaturated water flow. Agric. Water Manage. 29:201–213.
- Hupet, F., S. Lambot, R.A. Feddes, J.C. van Dam, and M. Vanclooster. 2003. Estimation of root water uptake parameters by inverse modeling with soil water content data. Water Resour. Res. 39:1312–1328.
- Inskeep, W.P., J.M. Wraith, J.P. Wilson, R.D. Snyder, R.E. Macur, and H.M. Gaber. 1996. Input parameter and model resolution effects on predictions of solute transport. J. Environ. Qual. 25:453–462.[Abstract/Free Full Text]
- Jarvis, N.J., L. Zavattaro, K. Rajkai, W.D. Reynolds, P.-A. Olsen, M. McGechan, M. Mecke, B. Mohanty, P.B. Leeds-Harrison, and B. Jacques. 2002. Indirect estimation of near-saturated hydraulic conductivity from readily available soil formation. Geoderma 108:1–17.[ISI]
- Larsson, M.H., N.J. Jarvis, G. Torstensson, and R. Kasteel. 1999. Quantifying the impact of preferential flow on solute transport to tile drains in a sandy field soil. J. Hydrol. (Amsterdam) 215:116–134.
- Lin, H.S., K.J. McInnes, L.P. Wilding, and C.T. Hallmark. 1999. Effect of soil morphology on hydraulic properties: II. Hydraulic pedotransfer functions. Soil Sci. Soc. Am. J. 63:955–961.[Abstract/Free Full Text]
- Marques, R., and J. Ranger. 1997. Nutrient dynamics in a chronosequence of Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) stands on the Beaujolais Mounts (France). 1: Qualitative approach. For. Ecol. Manage. 91:255–277.
- Mayer, K.U. 1999. A numerical model for multicomponent reactive transport in variably saturated porous media. University of Waterloo, Waterloo, Canada.
- Mayer, K.U., E.O. Frind, and D.W. Blowes. 2002. A numerical model for the investigation of reactive transport in variably saturated media using a generalized formulation for kinetically controlled reactions. Water Resour. Res. 38:1174–1194.
- Mohanty, B.P., R.S. Bowman, J.M.H. Hendrick, and M.Th. van Genuchten. 1997. New piecewise-continuous hydraulic functions for modelling preferential flow in an intermittent flood-irrigated field. Water Resour. Res. 33:2049–2063.
- Mualem, Y. 1976. A new model predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
- Musters, P.A.D., and W. Bouten. 2000. A method for identifying optimum strategies of measuring soil water contents for calibrating a root water uptake model. J. Hydrol. (Amsterdam) 227:273–286.
- Nemes, A., M.G. Schaap, and J.H.M. Wösten. 2003. Functional evaluation of pedotransfer functions derived from different scales of data collection. Soil Sci. Soc. Am. J. 67:1093–1102.[Abstract/Free Full Text]
- Othmer, H., B. Diekkrüger, and M. Kutilek. 1991. Bimodal porosity and unsaturated hydraulic conductivity. Soil Sci. 152:139–149.
- Pachepsky, Y.A., and W.J. Rawls. 2003. Soil structure and pedotransfer functions. Eur. J. Soil Sci. 54:443–451.
- Ranger, J., R. Marques, M. Colin-Belgrand, and N. Flammang. 1995. The dynamics of biomass and nutrient accumulation in a Douglas-fir (Pseudotsuga menziesii Franco) stand studied using a chronosequence approach. For. Ecol. Manage. 72:167–183.
- Ranger, J., R. Marques, and J.-H. Jussy. 2001. Forest dynamics during stand development assessed by lysimeter and centrifuge solutions. For. Ecol. Manage. 144:129–145.
- Rawls, W.J., Y. Pachepsky, and M.H. Shen. 2001. Testing soil water retention estimation with the MUUF pedotransfer model using data from the southern United States. J. Hydrol. (Amsterdam) 251:177–185.
- Roulier, S., and N. Jarvis. 2003. Analysis of inverse procedures for estimating parameters controlling macropore flow and solute transport in the dual-permeability model MACRO. Vadose Zone J. 2:349–357.[Abstract/Free Full Text]
- Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated hydraulic conductivity with the Mualem–van Genuchten. Soil Sci. Soc. Am. J. 64:843–851.[Abstract/Free Full Text]
- Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. (Amsterdam) 251:163–176.
- Simunek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. (Amsterdam) 272:14–35.
- Soet, M., and J.N.M. Stricker. 2003. Functional behaviour of pedotransfer functions in soil water flow simulation. Hydrol. Process. 17:1659–1670.
- Soil Survey Staff. 1994. Keys for soil taxonomy. 6th ed. USDA-SCS, Washington, DC.
- Tiktak, A., and W. Bouten. 1992. Modelling soil water dynamics in a forested ecosystem. III: Model description and evaluation of discretization. Hydrol. Process. 6:455–465.
- Tiktak, A., and W. Bouten. 1994. Soil water dynamics and long-term water balances of a Douglas fir stand in the Netherlands. J. Hydrol. (Amsterdam) 156:265–283.
- van Alphen, B.J., H.W.G. Booltink, and J. Bouma. 2001. Combining pedotransfer functions with physical measurements to improve the estimation of soil hydraulic properties. Geoderma 103:133–147.
- van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
- Vanderkwaak, J.E., P.A. Forsyth, K.T.B. MacQuarrie, and E.A. Sudicky. 1997. WatSolv- sparse matrix interative solver, user's guide for version 2.16. University of Waterloo, Waterloo, Canada.
- Vilette, S. 1994. Etablissement du bilan hydrique sur une chronosequence de peuplement de Douglas (Pseudotsuga menziesii (Mirb.) Franco). (In French.) ENESA-ENITA/INRA, Dijon.
- Volaire, F., and F. Lelièvre. 2001. Drought survival in Dactylis glomerata and Festuca arundinacea under similar rooting conditions in tubes. Plant Soil 229:225–234.
- Vrugt, J.A., J.W. Hopmans, and J. Simunek. 2001. Calibration of a two-dimensional root water uptake model. Soil Sci. Soc. Am. J. 65:1027–1037.[Abstract/Free Full Text]
- Wagner, B., V.R. Tarnawski, V. Hennings, U. Muller, G. Wessolek, and R. Plagge. 2001. Evaluation of pedo-transfer functions for unsaturated soil hydraulic conductivity using an independent data set. Geoderma 102:275–297.
- Wösten, J.H.M., C.H.J.E. Schuren, J. Bouma, and A. Stein. 1990. Functional sensitivity analysis of four methods to generate soil hydraulic functions. Soil Sci. Soc. Am. J. 54:832–836.[Abstract/Free Full Text]
- Wösten, J.H.M., A. Lilly, A. Nemes, and C. Le Bas. 1999. Development and use of a database of hydraulic properties of European soils. Geoderma 90:169–185.[ISI]
- Wösten, J.H.M., Ya. A. Pachepsky, and W.J. Rawls. 2001. Pedotransfer functions: Bridging the gap between available basic soil data and missing soil hydraulic characteristics. J. Hydrol. (Amsterdam) 251:123–150.
TOP
ABSTRACT
INTRODUCTION
