- Email: [email protected]
A distributed physical approach for surface-subsurface water transport modeling in agricultural watersheds
Journal &fdOgy Journal of Hydrology 203 (1997) 79-92
A distributed physical approach for surface-subsurface water transport modeling in agricultural watersheds F. Bouraoui*, Luboratoire
G. Vachaud, R. Haverkamp,
d’Etude des Transferts en Hydrologic
et Environnement
38041,
(LTHEYINPG.
B. Normand UJF, CNRS IJMR 5564).
BP 53X
Grenoble, France
Received 13 December 1995; accepted 20 June 1997
Abstract Surface cover and soil type have a major influence upon groundwater recharge and groundwater quality in agricultural watersheds. However, several hydrological models focus on simulating groundwater recharge without including the influence of agricultural practices and soil characteristics. In this study, ANSWERS, a distributed parameters surface nonpoint source model has been modified to include the simulation of water transport in the vadose and saturated zones. This model takes into account the spatial and temporal variability of crop cover and management practices, and the spatial variability of soil type and rainfall distribution. It is physically based and uses parameters that can be easily determined from readily available soil and plant information. It has been validated at multiple scales: local scale, field scale and watershed scale. At the local and field scale, it predicts accurately drainage below the root zone and evapotranspiration on different type of soil cover. At the watershed scale, it reproduces well the piezometric levels and trends of variation. 0 1997 Elsevier Science B.V. Keywords: Model; Physical approach; Simulation; Surface; Subsurface; Water transport; Agricultural watersheds
1. Introduction Management of surface and groundwater quality has emerged as an environmental priority, particularly in agricultural watersheds, where conjunctive use of water for agriculture and human needs, represents a major problem. Modeling is the most cost effective means of determining the impact of management alternatives on water resources. In the past, modeling efforts have either focused on developing surface water or groundwater models. Only recently, more emphasis was given to the interactions between the saturated and unsaturated zones. Many of the comprehensive hydrologic models were developed for use at small scales, whereas the need of application is mostly * Corresponding author. 0022-1694/97/$17.00 0 1997 Elsevier Science B.V. All rights reserved PII SOO22- 1694(97)00085- 1
at the watershed scale. An integrated approach is thus necessary for a better understanding of the dynamic interactions between the different processes involved. Three essential aspects have to be addressed in the study of the transport processes through the soil to the groundwater: ??
??
.
modeling of the whole watershed as an entire system; characterization of the spatial and temporal variability of the soil characteristics and management type; and influence of vegetation and agricultural practices on the quality and the quantity of the groundwater recharge.
In this sense, it appears that a continuous distributed approach is probably the most appropriate. The
use of distributed models allows the determination of the response at each point within the watershed, and the identification of the critical areas such as the regions of a watershed which are likely to contribute the most to the recharge and/or to the contamination of groundwater. For a complete review of available watershed models the reader is referred to Singh (1995). During the last decade, two general tendencies have emerged in distributed modeling. One consists of a it uses the fundamental mechanistic approach; equations of conservation of mass and energy at the large scale, a typical example being the SHE model (Abbott et al., 1986). This approach can become quickly limiting because of overparametrization and the lack of readily available data. Furthermore, because of the discrepancy between the scales of the measurement and the model grid size, lumping has to be done (Beven, 1989). The other type of approach uses capacity storage-type models instead. This type of models tries to characterize the average response of each cell rather than to capture in detail the physics and the variations that might occur within the cell. To illustrate the efficiency of such an approach, we want to present in this article a time continuous, spatially distributed parameters model developed recently. At this stage, only the hydrodynamic aspect will be considered. The model simulates the entire hydrologic cycle. Simple procedures based on soil properties and characteristics are used to determine the required parameters. Special emphasis has been given to the influence of coarse fragments on the physical and hydrodynamic properties of soils. The model was first used at the grid scale to determine its sensitivity to different management practices and its ability to simulate soil water transport. It was then used at the watershed scale, with limited calibration, with an appropriate spatially and temporarily distributed management scenario and climatic input.
2. Model description 2.1. Introduction The following modeling approach is based on the concept of ANSWERS (Area1 Nonpoint Source
Flow Path
Soil Horizons
Piazometric Head
Fig. I. Representation ANSWERS model.
of the hydrologic
cycle as simulated
by the
Watershed Environment Response Simulation; Bouraoui and Dillaha, 1994; Bouraoui and Dillaha, 1996): a distributed parameters, nonpoint source pollution surface model for long term simulation of infiltration, runoff, sediment transport and nutrient (nitrogen and phosphorus) transport and transformations. The original work was mostly focused on prediction (both qualitative and quantitative) of runoff, sediment transport and drainage below the root zone. In this study, the hydrodynamic part of the model has been modified to take into account the interaction between vadose zone and aquifer by including multiple layers. A conceptualization of the watershed representation is given Fig. 1. The watershed is discretized in volume units with square surface areas; one hectare is usually suggested, however a larger size can be considered for homogeneous conditions. Each block is characterized by topographic, soils hydrodynamic properties and crops characteristics. Vertically the blocks are discretized in different layers to account
F. Bouraoui et al./Joumal of Hydrology 203 (I 997) 79-92
<
Precipitation >--f--i
Evapotranspiration
81
1
Yes lnitiallze rimuktlon time t=o 1 -
layers 1,2 and 3
6
Net Pmcipitation
1
Groundwater discharge
infiltration/Runoff
Update plezometric
i
t=t+NAt (At = 60 s)
DAY = DAY +1
I
.
No
Is the simulation
Output summary
Fig. 2. Flow chart for the ANSWERS
for water movement through the soil profile from the surface to the aquifer. The soil profile is divided into three vertical layers, two representing the root zone and a third one extending from beyond the root zone to the water table. A fourth layer of variable thickness is used to represent the saturated zone. As illustrated by the flow chart (Fig. 2) the following flow phenomena are considered: infiltration, surface flow, drainage, evapotranspiration and groundwater flow. During a rainfall event, infiltration starts once the interception volume, which is -plant
model.
dependent, is satisfied. Infiltration is computed based on the Green-Ampt model. At the end of the rainstorm, soil water redistribution is determined based on the assumption of a gravity flow (head gradient equal unity) under unsaturated conditions. The evapotranspiration, the recharge of the aquifer, which is the result of water percolating out of the vadose zone, and the groundwater discharge are computed daily. A daily water balance is maintained. More details will be given about the drainage and groundwater components because they have been modified from the original ANSWERS model.
3. Model overview At the surface, the net rainfall (i.e. rainfall minus interception) is subject to infiltration and runoff for each time step and each block. Infiltration starts once the interception volume is filled. A maximum intercepted volume, dependent on the plant type, is interpolated linearly according to the plant stage of growth. Plant growth is represented by a varying leaf area index (LAI) and root zone. For each crop cover, tabulated values of the LA1 (Knisel, 1980) are used as input to the model for 10 stages of the plant growth. A linear interpolation is made on a daily base for each stage. The root depth, which determines the extent of the root zone, is computed using a sinus function related to the number of days after planting and the number of days to reach maturity (Borg and Williams, 1986). The intercepted water is then returned to the atmosphere through evaporation. Evapotranspiration is determined based on the approach of Ritchie (1972). The potential evapotranspiration is computed with the Priestly-Taylor equation (Priestly and Taylor, 1972) using the net solar radiation and mean air temperature. The daily value of LA1 is used to partition the potential evapotranspiration into potential soil evaporation and potential plant transpiration, the latter term being limited by the amount of water in the root zone. Soil evaporation is assumed to take place in two different stages. In the first stage (constant rate stage), soil evaporation is just limited by the amount of energy available and proceeds at the potential rate. Beyond an upper limit depending of the soil water transmission characteristics the second stage begins. During this stage, the falling rate stage, soil evaporation rate decreases as a square root function of time. For any given day, the sum of plant transpiration and soil evaporation cannot exceed potential evapotranspiration. Infiltration is simulated using the physically-based Green-Ampt equation (Green and Ampt, 1911). It assumes (i) a step water retention function h(8) (wetting front capillary potential) describing the relation between volumetric water content, 19(cm’/cm”), and soil water pressure, h (cm); and (ii) a one-point hydraulic conductivity function K, (cm/h) where K, is the hydraulic conductivity for volumetric water content at natural saturation 8,. The infiltration process is represented as a wetting front traveling down
the soil protile with 0 = 8, and K = K, behind rhe wetting front and 0 = 0,) and K = 0 ahead of the wetting front. The saturated hydraulic conductivity and the wetting front capillary potential are determined from textural soil properties using statistical regression functions (pedotransfer functions) given by Rawls and Brakensiek ( 1989). Once the infiltration is determined, rainfall excess is routed to the watershed outlet. A modified discharge-water depth equation of the Manning type is used for both overland and channel flow routing. Every square element of the discretized watershed acts as an overland flow plane with a user specified slope and slope direction. Channel elements (described in terms of their slope, width and roughness coefficient) collect flow from overland flow elements and route the runoff to the watershed outlet. The drainage in the unsaturated zone is determined based on Darcy’s approach on the assumption of a gravity flow under unsaturated conditions. Travel time of percolating water through the soil matrix is regulated by the unsaturated hydraulic conductivity. This conductivity K(0) is modeled using the Brooks and Corey ( 1964) equation: / a\‘i (1) where 11represents the conductivity shape parameter, 0, (cmj/cm”), the water content at natural saturation. both determined from pedotransfer functions (Rawls and Brakensiek, 1989). Travel time [ 7T, (h)] through a particular soil layer is computed using a linear storage equation:
where d, is the thickness of the specific layer (cm). Drainage DR (cm) during a specific time step, At (h) is determined using an exponential function (Williams et al., 1985):
1 The water percolating out of the last unsaturated layer (Fig. 1) is added to the amount of water present in the groundwater cell and becomes available for lateral saturated flow. Because of the high infiltrating properties of the soils in the study area, it is assumed that there is no
F. Bouraoui et al./Joumal
Table 1 Major soil parameters
of Hydrology 203 (1997)
83
79-92
used in the model Layer O-30 cm
Bulk density (g/cm’) Clay fraction (%)** Sand fraction (%)** Saturated conductivity (cm/h) Wetting front capillary potential, ht (cm) Saturation water content, 6, (%)*** Conductivity shape parameter, n Volume of coarse fragments, VP (%)**** Fine fraction. or (%)*
1.339 17.5 41.0 1.19 26.1 90 19.6 25.2 60.0
Layer 30- 120 cm 1.316
16.4 51.7 3.44 27.3 90 13.4 49.4 33.7
* percentage of the total weight. ** percentage of total weight of fine particles. *** percentage of total porosity. **** percentage of total volume.
lateral subsurface flow from one cell to another in the unsaturated zone. Thus, the unique recharge of an aquifer cell is the excess water draining from the corresponding overlying cells. The draining water is added to the groundwater volume present in each aquifer cell. The total groundwater volume is then divided by the specific storage to obtain the piezometric head for each cell (Hi). An average daily value of aquifer discharge, Q (m3/ day), is computed at the downstream limit using a Darcy-type equation, such as: Q-m:
(4)
where K represents the conductivity of the aquifer material (m/day), M/AZ is the average hydraulic gradient (m/m) between the upstream and downstream limit of the aquifer, and A is the wetted cross section of the aquifer (m*). The resulting total outflow obtained from Eq. (4) is distributed to each cell and a new piezometric head is determined for each cell.
4. Experimental layout and data description 4.1. The site The site selected for the calibration and validation is ‘La C&e St Andre’, 60 km northeast of Grenoble (southeast of France). The choice of this site is related
to important threads concerning the qualitative and quantitative degradation of the groundwater resources with the development of intensive agricultural practices. It should be noted that the groundwater is the only drinking water supply for the local collectivities; its evolution is a source of conflict between farmers and cities. For this reason, an intensive, pluriannual, interdisciplinary study was initiated in 1991 in order to optimize agricultural practices, and to develop sustainable management schemes. Special attention was given to the characterization of the relationships between fertilization, irrigation and crop production, in particular to the water and nutrient balance in the soil-plant-atmosphere continuum (Kengni et al., 1994). Another important objective of the study is to predict the evolution with time of the groundwater resource. Nitrate is the major nonpoint source groundwater pollutant with nearly one half of the wells reaching, or overpassing the European limit of drinkability (50 ppm). 4.2. The scales of investigation Three different scales have been intensively investigated: local, field and watershed. At the local scale, measurements were carried out on an undisturbed lysimeter (1.2 m diameter, 1.5 m height), with bare soil or grass coverage, in order to obtain direct measurements of drainage and evapotranspiration and to characterize in detail the hydrodynamic transport properties in the soil. The lysimeter was
Table 2 Monthly rainfall and potential
evapotranspiration
1991
February March April May June July August September October November December
for the study site I992
I993
19V4
Rain
PET
Rain
PET
Rain
PET
Rain
18.1 26.6 118.4 51.6 21.6 92.8 49.6 21.9 94.4 133.6 53.3 93
18.3 3 I .6 49.1 17.3 106.8 103.4 163.2 153.5 92. I 41.4 23.5 20
25.8 37.4 75.9 57.9 96.9 120.6 99.7 51.5 93.3 182.1 114.1 64.1
18.4 29. I 50 77.9 106.7 88.6 127.8 141.6 76 28. I 19.2 21.4
2X.8 9.5 8.7 110.2 104.1 86.6 129.7 36.2 258.3 250.7 26.9 83.4
26.6 24 59.9 73.7 105.9 126.8 129.6 122.8 64.5 41.5 23.5 24.4
126.4 26.9 15.3 126.5 88.2 45.4 I7 37.9 192.8 80.4 107.9 30.5
25.7 61.5 55.4 95.1 129.4 170. I 156.5 ?4.3 49..3 33 ._’ -_ 14.9
780.9
880.2
823.2
8V5.2
886
1019.3
instrumented with neutron probe access tubes, tensiometers and soil temperature gauges, and direct outflow was automatically read. At the field scale, experiments were carried out on two types of soil cover: maize and bare soil. Measuring devices identical to those installed in the lysimeter were set up on two plots, one for each treatment, to determine the water balance in the root zone (Kengni et al., 1994). At the watershed scale, three different sets of data are looked upon: the soil, the soil cover (crop cover and management practices) and the aquifer (boundaries, transmissivities, and piezometric level and water quality variations). The watershed is about 320 km2 and the elevation varies from 480 to 250 m. The major part of the watershed (about 80%) is a flat plain with very permeable, fertile, unconsolidated soil which is surrounded by hills, with an average slope of 7.7%. Hills are covered with grass and forest; the plain is cultivated either with irrigated maize and tobacco (20% of the area) or with dry farming crops, such as sunflower and corn. In the plain, the arable layer is shallow (average 1 m) and rich in organic matter. There is almost no hydrographic network due to the high infiltrability of the arable layer which consists of a sandy loam becoming progressively coarser with increasing depths. Roughly speaking, the soil is characterized by two different layers (Table 1).
784.8
1133.1
PET 2 I .6
The groundwater holding layer consists of a very coarse glacial material extending from a depth of 1 m to about 20-30 m, with a water table aquifer. Due to the presence of coarse fragments, it was not possible to install any measuring device to characterize transport in the unsaturated part of this glacial material. The aquifer is characterized by very high transmissivities (range from 5 x lo-* to 10-l m*/s), high velocities (between 5 and 10 m/day). Consequently, the average discharge is approximately 4 m3/s. The effective porosity of the aquifer is about lo%, and its depth fluctuates around 5 m during a year, with a typical range of average depth of 10 - 15 m below the soil surface. The changes of piezometric levels may be very fast during the rain season. The average annual rainfall over 30 y is around 1000 mm, with two characteristics: a very high interannual variability and an annual pattern with the most important rain events always in September-October, and another rain season, weaker and more uncertain, in April-May (Table 2). Mean evapotranspiration computed by the Penman-Montieth approach is around 850 mm/y. 4.3. Available
data
In the lysimeter and on the bare and corn plots, soil water pressure, and volumetric soil water content were measured at different depths in the upper soil
F. Bouraoui et al./Journal
layer (thickness 1 m) during the entire crop season, at a time step varying from 30 min to one week. Flux draining out of the lysimeter were measured continuously, whereas on the field sites, the fluxes below the root zone were estimated on a daily basis by the use of Darcy’s law following a method described previously (Kengni et al., 1994), and resting on time course measurement of soil water content, soil water pressure, and direct estimation of the K(8) relationship. No measurement of the flux and soil transport parameters could be done in the underlying layer because of the presence of coarse gravel and pebbles. The evaporation on the bare soil plot and the evapotranspiration on the grass lysimeter and maize plot were determined by mass balance as the difference between the estimated drainage and the changes in soil water storage. All relevant micrometeorological data (rainfall, radiation, pan evaporation, etc.) were measured, at a time step of 30 min, during the entire experiment. Finally, at the watershed scale, the actual availability of data concerns the following items: in term of the aquifer, an intensive geophysical series of measurements (seismic and electric) were made 10 years ago to determine the geometry of the bedrock; in terms of soil, a survey was conducted by the Regional Chamber of Agriculture resulting in a database including texture, organic matter, and density; concerning the soil coverage and management, a map was made by the Institut Supkrieur d’Agriculture Rh6ne Alpes (ISARA, Lyon) following a systematic survey of farming practices. Continuous measurements of the piezometric level (on the basis of weekly values) are available for three years for five different wells.
of Hydrology 203 (1997) 79-92
X5
measured data mean is a better predictor than the predicted data). The mean and standard deviation are also reported. Care should be given, however, when interpreting these values since the measurement periods can vary from plot to plot (no crosscomparison is possible for the different types of cover). 5.2. Grid scale validation The model was first used at a grid scale to obtain a validation between the model output and the measurement of water flux for different surface cover. This internal validation at the grid scale is an essential step before conducting simulations at the watershed scale in order to ensure the proper representation by the model of the different components of the water balance for different surface covers. It was assumed that within the grid, the soil and its cover have the same characteristics as those determined on the measurement sites. Three different applications were carried out: first on the bare soil lysimeter (later planted with grass) with comparison between simulated and measured drainage at the lysimeter outflow and real evapotranspiration; then on the bare and corn plots with comparison between simulated and estimated drainage. On the lysimeter about 18 months of measured drainage was available. The lysimeter was installed in September 1992 and was kept bare until June 1993 after which it ‘was planted with grass. With 100
g
E I0 E
5. Validation 5.1. Statistical
6 zii analysis
1
s 0
5 The accuracy of the model predictions were analyzed using a correlation coefficient which gives a measure of linear association between the predicted and measured data, and a coefficient of efficiency which represents the fraction of total variance of the measured data explained by the predicted data (a negative value of this coefficient indicates that the
0.1
E
u
I” 0.01 0.15
0.25
0.35
0.45
Water content cm3km3 Fig. 3. Hydraulic conductivity curve for the lysimeter.
Table 3 Measured
and simulated
R*
Cover type
Year
Bare soil
91-92 92-93 93-94 91-92 92-93 93-94 92-93
Corn
Lysimeter
actual evapotranspiration
640 909
II00 640 909 II00 1394
(AET, mm) and drainage
(D, mm) for different
surface cover
PET*
AET,
AET,
AETJAET,
D,
D\
DJD
726 595 514 726 595 574 962
320 392 379 495 514 548 609
355 345 330 580 524 513 663
1.109 0.880 0.87 I 1.172 I.019 0.936 1.089
516 572 753 342 467 551 908
so4 560 762 280 478 559 955
0.977 0.979 I.(112 0.819 1.024 1.015
,n
I .os1
-
* R: rainfall
(mm). ** PET: potential evapotranspiration m measured. s simulated.
(mm).
regard to the physical properties of the soil, the soil profile was discretized in two horizons: O-O.3 m and 0.3- 1.2 m. The values of K, and hr are computed from the pedotransfer functions developed by Rawls and Brakensiek (1989) and are given in Table 2, along with the other parameters necessary to run the model. The saturated hydraulic conductivity and the porosity are adjusted for the presence of coarse fragments according to the procedure described in Appendix A. A comparison between the values predicted by the pedotransfer function and the measured K(0) is given in Fig. 3. It must be noted that due to specificities of the ‘zero flux’ method, it is not possible to obtain measured values of K(8) in the wet range. The results for cumulative model predictions for drainage and real evapotranspiration (AET),
Table 4 Statistical
without any calibration are very satisfactory (Fig. 4). Total drainage and AET are well reproduced (Table 3) and are within 10% of the measured data. The performance of the model is quite acceptable considering that no adjustments were made (Table 4). Identical soil characteristics were used in the simulation of the corn and the bare plots. Fig. 5 represents the measured and the simulated cumulated drainage, and the estimated and simulated AET, respectively, for the bare plot for three years of data. The predictions of drainage and evaporation were good for the three years of simulation. Similar results are also obtained on the corn plot (Fig. 6). For both surface covers, the cumulative predictions are within 20% of the measured data (Table 3). The predictions of drainage are better than those of AET and this
analysis for all simulations
Actual evapotranspiration Baresoil Corn lysimeter Bare soil Corn Lysimeter Well
11.48 17.28 16.46 Drainage 19.27 15.12 15.67 Piezometric 342.29
* correlation coefficient. ** coefficient of efficiency m measured. s simulated.
10.70 17.96 17.94
8.52 13.98 1I .50
5.80 9.88 9.54
0.667 0.741 0.834
0.439 0.553 0.683
23.14 14.61 14.88
26.82 26.82 27.82
25.85 25.85 30.48
0.933 0.916 0.881
0.873 0.837 0.731
342.46
I.12
I .08
0.957
0.892
head
.
0001 I
..
'(pa1qIz,Q.m03
doJ3 aql Bu!%ueqs l~ql alexpu! satyeA paJnseam ~(ai!pm ‘SWJ% 'pas mq)
aaql
~03 AlaleJnm? apnb
J1AO3 JO SadB luaJaJJ!p pauuoyad
11
-sacylxzJd
JO sadlClaaJylaqlJojL~vJoJ wql ah.qmp
qmu
a~ slInsaJaqL -clappaJnwaur aq
_103~awq LIalenbape
lwld aql 30 uopmy!~dLu~sJaAo aql 01 anp aq hxu
lCJaAsf Iapow aq3 leqlmoqs sllnsaJsno!AaJdaqJ
pampoJdaJ IapouIaql ‘aseqLIyaaM B ug -1uauodwos
luama%euk?wpu~JaAo~doJ3u~ sa%.mq301 aA!suodsaJ
.Ja]aur!s/([ aq) 104 13~
pm a%eu!wp pawpaJd SA pamsea~
‘p +?!,I
e&u!wp pew~nw!s S&“!EJp
pamsew
??
I i ooz Z6J@HI -t----f 3 - ooz -8 ii 3
oat
3 d s
18
F. Bouraoui et d./Journulof Hwirology 203 (IYY7) 70-02
88
600 1
-
. measured AET simulated AET
4oot
E-
s fk e .e
. 600
b
600
1
-
measured drainage simulated drainage
D
1000
Fig. 6. Measured
vs predicted
cover from maize to bare conditions, results into an increase of drainage of 44%, 14% and 28% for the years 1991, 1992, 1993, respectively. The model estimated these increases at 33%, 19% and 26%. Finally, the model predicted well the fact that most drainage does not occur during the growing season of the plants, even during irrigation, but during the winter, when the soil is bare and the values of AET are very low. 5.3. Watershed
validation
The next step was the application of the model to the entire plain to test its ability to predict the recharge of the aquifer as influenced by the surface cover management and soil depth to the aquifer. The soil was discretized into three layers: two layers in the root zone, and a third one which extends from below the root zone to the top of the aquifer. It is assumed that all the water leaving the root zone is transferred in its entirety to the aquifer. The surrounding boundary conditions for the aquifer were assumed to be a no-flow boundary condition as suggested from previous hydrogeological studies: the unique recharge comes from atmospheric water. Fig. 7(a) illustrates the elevation map of the underlying impervious layer. This
drainage
and AET for the corn plot
map of the substratum was interpolated using well profile information. Fig. 7(b) provides the piezometric map for the whole watershed for 7 January 1991. The piezometric map was kriged using 78 data points spread throughout the region; additional details about the kriging procedure can be found in Ruy (1992). As mentioned earlier, three years of continuous measurements of water level at five different wells, one of them (Feuges) located near the study site in the heart of the irrigated area, are available. The vegetation parameters, and the water input (rain or rain + irrigation) were distributed over the plain according to Fig. 7(c). The first six months of the data set (Feuges well) were used for calibrating the transfer parameter in the unsaturated zone extending from below the root zone to the water table (third layer). The calibration consisted in adjusting the hydraulic conductivity until the predicted and measured value of the piezometric head at the Feuges well were in agreement in terms of timing of the rise of the water level. The transfer parameters in the root zone (first and second layers) were determined from the pedotransfer functions developed by Rawls and Brakensiek (1989). The model was then run for 30 months in a purely predictive way. The prediction for the water level for the
x9
F. Bouraoui et alJJourna1 of Hydrology 203 (1997) 79-92
a) IMPERVIOUS
LAYER
ELEVATION
(m)
Feuges well
b) PIEZOMETRIC
c) LAND SURFACE
Fig. 7. Substratum depth (a), piezometric (above sea level).
HEAD
AND CROP
(m)
COVER
head (b), land surface and crop cover for the plain of Bievre. Isovalue lines represent elevation in m
‘Feuges’ well is given in Fig. 8. The piezometric head is expressed in meters above sea level, the soil surface elevation is 356.10 m. The model reproduced well the variations of the water table level (Table 4), in particular the high level of water during winter and early spring when the soil is completely bare (Fig. 8). The model also predicted that hardly any significant rise of the aquifer level will occur during the summer. This is mainly due to the presence of maize, a highly transpirating plant which limits the amount of deep percolation. Thus, during these periods of no deep percolation, there
will be no risk of nitrate leaching. This has already been noticed by Kengni et al. (1994). The model calibration leads to an estimate of travel time through the vadose zone of about three weeks for heavy rain periods. This value indicates how quick the aquifer responds to the excess rainfall. Finally, an increase of the water table level from the beginning to the end of the simulation period (3 y) by more than 4 m is properly predicted. This indicates that the model is very responsive to large variation in climatic input. The results for the Feuges well were confirmed by the four additional wells (Fig. 9).
90
measured 340
339 _ _.~_ Nov.90
Fig. 8. Predicted
-
1
~~_ ~~~ em- _~_ Jun-91
YS measured
Deco1
piezometric
~~
Jul-92
~
Jan-93
-
simulated
~-. Mar-94
Aug-93
level for the ‘Feugea’ well (soil surface elevation
6. Summary and conclusions
layers, and excess water from the root zone is routed to the aquifer. A Darcian-type equation was used to estimate groundwater discharge. Local measurements on a bare soil lysimeter, a bare and a corn plot were used to validate the model in the root zone over three years. Since the model is intended to be used for predictive purposes on watersheds where measured information is limited or lacking, indirect methods of estimation of soil transfer properties (pedotransfer function) and tabulated values of plant growth were
A continuous, distributed nonpoint source model for estimating groundwater recharge as affected by the spatial distribution of soil property characteristics, spatial and temporal variability of management practices, was developed. The model simulates infiltration by the Green-Ampt method. Ritchie’s approach is used to estimate plant transpiration and soil evaporation separately. The soil is discritized into different
- well . well 0 well 0 well . well
1
2 3 4 6
-+-_ ._~~ 336
340
344
346
Measured piezometrlc head (m)
Fig. 9. Predicted
350.10 m)
vs measured piezometric
level for the five validation
wells
91
F. Bouraoui et al.Nournal of Hydrology 203 (1997) 79-92
privileged over the measured characteristics. The information obtained at the local scale was then used to run the model on a large watershed. The model performed quite well in reproducing the piezometric head values and variations at several wells. However, it is important to note that the research described above showed that a simple model, with little input data required, can be a powerful tool for predicting and managing natural water resources. It may be added that the proper functioning of the model certainly results from a good agreement between its main hypothesis (infiltration governed by Green and Ampt assumption) and the type of soil characterizing the watershed. Questions remain open about the efficiency of this type of model to predict groundwater recharge in clayey material.
Acknowledgements This study was realized with the joint support of a European Union Program (EV5V-CT94-0484, ‘Evaluation of the Effect of Climatic Variations on the Recharge of Aquifers in Southern European Catchments’, DGXII) and of the Centre National pour la Recherche Scientifique, CNRS, Paris, Programme Environnement, Vie et SociCtC.
Appendix
A Adjustment
for Coarse Fragments
Most of the measurement techniques give an estimate of the whole-soil volumetric water content (e,,, cm ‘/cm 3). The representative x/olume for which these water content values are determined may include not only fine soil but also coarse fragments and stones, when present (as in our case study). The terms ‘coarse fragments’ and ‘stones’ as used in this paper concern mineral fragments more than 0.2 cm in diameter regardless of their specific size and shapes (Soil Survey Laboratory Staff, 1992). The whole-soil bulk density, pap (g/cm’), applies to the total sample volume [VT (cm’)] for which 8,, is determined. The usual bulk density, pf (g/cm’), refers to the volume of fine soil, Vf (cm3), and the specific density pc (g/cm’) to the volume fraction of coarse fragments, V, (cm”) which are linked through:
Writing
the conservation
M,+Mf=M,
(Al)
(A2)
where M,, Mr and MT represent the mass of coarse fragments, the mass of fine material and the total mass, respectively. Eq. (Al) can be rewritten as: @
Mf
(A3)
p_+-=vT ‘ Pf
The combination of Eq. (A2) and Eq. (A3) allows the calculation of the fine soil bulk density as: Pf =
‘yf Mw --___
(A4)
(l-af)
&i,MT
PC
where lyf is the fraction by weight of the fine material and A4w is the mass of soil water of the whole-soil sample determined through the process of weighing, drying and reweighing of the sample. Taking pc equal to 2.65 g/cm” and pr being known from Eq. (A4), the volume fraction of coarse fragments, Vi’ defined as the ratio V, and VT becomes: (1 -at.)
v,“=
l+cr. t
(AS)
p__,( Pt >
The volumetric tine soil water content and/or soil porosity is then adjusted by a correction factor equal to (1 - VP) and the saturated hydraulic conductivity is adjusted by cxf. Additional details about this correction procedure can be found in Haverkamp et al. (1997).
References Abbott,
M.B.,
Bathurst,
J.C..
Cunge.
J.A..
Connell,
P.E..
Rasmussen, .I.. 1986. An introduction to the European Hydrological System (Systbme Hydrologique structure of a physically-based,
Europeen.
‘SHE’,
2:
distributed modeling syctem).
Journal of Hydrology, 87, 6 I-77. Borg. H., Williams,
D.W.,
1986. Depth development of roots with
time: An empirical description, Transactions of the ASAE. (I).
29
194-197.
Beven,
K.J.,
1989.
Changing
ideas in hydrology:
physically based models, Journal of Hydrology. Bouraoui, F., Dillaha. T.A.,
1994. ANSWERS-2000.
the case of I OS, I S7- 172. Continuous
simulation version. ASAE Paper No. 94-2 I SS. presented at the ASAE
v, + v, = vr
of mass gives:
USA.
International
Winter
Meeting,
Kansas City. Miwourt.
Y2
I,‘. Rouruoui
rt ul.Nournul
Bouraoui. F.. Dillaha, T.A.. 1996. ANSWERS-2000: Runoff and sediment transport model. Journal of Env,ironmental Engineering, ASCE, 122 (6), 493-502. Brooks, R.H., Corey, C.T., 1964. Hydraulic properties of porous media. Hydrology Paper 3. Colorado State University, Fort Collins, USA. Green, W.H., Ampt, G., I91 I. Studies of soil physics, Part 1. The flow of air and water through soils, Journal of Agricultural Science, 4, l-24. Haverkamp, R., Arrue, J.L., Ross, P.J., Dalton. F.N., Lopez. Ma.V.. Bouraoui, F., Garcia. R., 1997. Influence of stones and porous coarse fragments on water content measurements and bulk density at field scale. Soil Science Society of America Journal (submitted). Kengni, L., Vachaud, G., Thony, J.L., Laty. R., Garino. B.. Casabianca, H., Jame. P., Viscogliosi, R., 1994. Field measurements of water and nitrogen losses under irrigated maize. Journal of Hydrology, 162, 33-46. Knisel, W.G., 1980. CREAMS: A Field Scale Model for Chemical, Runoff, and Erosion from Agricultural Management Systems. Conservation Research Report No.26, USDA-SEA.
of Hvdrolo,cv 203 (lYY7) 7Y--92 Priestly. C.H.B, Taylor, R.J., 1972. On the assessment 01 rurtacl, heat flux and evaporation using large-scale parametric. Monthly Weather Review 100, 8 I-92. Rawls. W.J.. Brakenstek, D.L.. 1989. Estimation of soil water retention and hydraulic properties. In: Morel-Seytoux. H.J. (Ed.). Unsaturated flow in Hydrologic Modeling: Theory and Practice. NATO ASI Series c: Mathematical and Physical Sciences, Vol.275. Kluwer Academic, Dordrecht. Holland. Ruy. S.. 1992. Etude de la nappe phreatique de la plaine de Ia Bievre. DEA. Laboratoire d’Etude des Transferts en Hydrologie et Environnement. Grenoble, France. Ritchie, J.T.. 1972. A model for predicting evapotranspiratton from a row crop with incomplete cover, Water Resources Research. 8 (51, 1204-1213. Singh, V.J.. 1995. Computer models of watershed hydrology. Water Resources Publications, Highlands Ranch, CO. Soil Survey Laboratory Staff, lY92. Soil survey laboratory methods manual. Soil Surv. Invest. Rep. 42 Version 2.0 IJSDA-SCS U.S. Gov. Print. Oftice. Washington, D.C., USA. Williams. J.R.. Nicks, A.D.. Arnold. J.G.. 1985. Simulator of water resources in rural ba.:ins. J. Hydr. Div., ASCE, I I I (6),970-9X6.
A distributed physical approach for surface-subsurface water transport modeling in agricultural watersheds F. Bouraoui*, Luboratoire
G. Vachaud, R. Haverkamp,
d’Etude des Transferts en Hydrologic
et Environnement
38041,
(LTHEYINPG.
B. Normand UJF, CNRS IJMR 5564).
BP 53X
Grenoble, France
Received 13 December 1995; accepted 20 June 1997
Abstract Surface cover and soil type have a major influence upon groundwater recharge and groundwater quality in agricultural watersheds. However, several hydrological models focus on simulating groundwater recharge without including the influence of agricultural practices and soil characteristics. In this study, ANSWERS, a distributed parameters surface nonpoint source model has been modified to include the simulation of water transport in the vadose and saturated zones. This model takes into account the spatial and temporal variability of crop cover and management practices, and the spatial variability of soil type and rainfall distribution. It is physically based and uses parameters that can be easily determined from readily available soil and plant information. It has been validated at multiple scales: local scale, field scale and watershed scale. At the local and field scale, it predicts accurately drainage below the root zone and evapotranspiration on different type of soil cover. At the watershed scale, it reproduces well the piezometric levels and trends of variation. 0 1997 Elsevier Science B.V. Keywords: Model; Physical approach; Simulation; Surface; Subsurface; Water transport; Agricultural watersheds
1. Introduction Management of surface and groundwater quality has emerged as an environmental priority, particularly in agricultural watersheds, where conjunctive use of water for agriculture and human needs, represents a major problem. Modeling is the most cost effective means of determining the impact of management alternatives on water resources. In the past, modeling efforts have either focused on developing surface water or groundwater models. Only recently, more emphasis was given to the interactions between the saturated and unsaturated zones. Many of the comprehensive hydrologic models were developed for use at small scales, whereas the need of application is mostly * Corresponding author. 0022-1694/97/$17.00 0 1997 Elsevier Science B.V. All rights reserved PII SOO22- 1694(97)00085- 1
at the watershed scale. An integrated approach is thus necessary for a better understanding of the dynamic interactions between the different processes involved. Three essential aspects have to be addressed in the study of the transport processes through the soil to the groundwater: ??
??
.
modeling of the whole watershed as an entire system; characterization of the spatial and temporal variability of the soil characteristics and management type; and influence of vegetation and agricultural practices on the quality and the quantity of the groundwater recharge.
In this sense, it appears that a continuous distributed approach is probably the most appropriate. The
use of distributed models allows the determination of the response at each point within the watershed, and the identification of the critical areas such as the regions of a watershed which are likely to contribute the most to the recharge and/or to the contamination of groundwater. For a complete review of available watershed models the reader is referred to Singh (1995). During the last decade, two general tendencies have emerged in distributed modeling. One consists of a it uses the fundamental mechanistic approach; equations of conservation of mass and energy at the large scale, a typical example being the SHE model (Abbott et al., 1986). This approach can become quickly limiting because of overparametrization and the lack of readily available data. Furthermore, because of the discrepancy between the scales of the measurement and the model grid size, lumping has to be done (Beven, 1989). The other type of approach uses capacity storage-type models instead. This type of models tries to characterize the average response of each cell rather than to capture in detail the physics and the variations that might occur within the cell. To illustrate the efficiency of such an approach, we want to present in this article a time continuous, spatially distributed parameters model developed recently. At this stage, only the hydrodynamic aspect will be considered. The model simulates the entire hydrologic cycle. Simple procedures based on soil properties and characteristics are used to determine the required parameters. Special emphasis has been given to the influence of coarse fragments on the physical and hydrodynamic properties of soils. The model was first used at the grid scale to determine its sensitivity to different management practices and its ability to simulate soil water transport. It was then used at the watershed scale, with limited calibration, with an appropriate spatially and temporarily distributed management scenario and climatic input.
2. Model description 2.1. Introduction The following modeling approach is based on the concept of ANSWERS (Area1 Nonpoint Source
Flow Path
Soil Horizons
Piazometric Head
Fig. I. Representation ANSWERS model.
of the hydrologic
cycle as simulated
by the
Watershed Environment Response Simulation; Bouraoui and Dillaha, 1994; Bouraoui and Dillaha, 1996): a distributed parameters, nonpoint source pollution surface model for long term simulation of infiltration, runoff, sediment transport and nutrient (nitrogen and phosphorus) transport and transformations. The original work was mostly focused on prediction (both qualitative and quantitative) of runoff, sediment transport and drainage below the root zone. In this study, the hydrodynamic part of the model has been modified to take into account the interaction between vadose zone and aquifer by including multiple layers. A conceptualization of the watershed representation is given Fig. 1. The watershed is discretized in volume units with square surface areas; one hectare is usually suggested, however a larger size can be considered for homogeneous conditions. Each block is characterized by topographic, soils hydrodynamic properties and crops characteristics. Vertically the blocks are discretized in different layers to account
F. Bouraoui et al./Joumal of Hydrology 203 (I 997) 79-92
<
Precipitation >--f--i
Evapotranspiration
81
1
Yes lnitiallze rimuktlon time t=o 1 -
layers 1,2 and 3
6
Net Pmcipitation
1
Groundwater discharge
infiltration/Runoff
Update plezometric
i
t=t+NAt (At = 60 s)
DAY = DAY +1
I
.
No
Is the simulation
Output summary
Fig. 2. Flow chart for the ANSWERS
for water movement through the soil profile from the surface to the aquifer. The soil profile is divided into three vertical layers, two representing the root zone and a third one extending from beyond the root zone to the water table. A fourth layer of variable thickness is used to represent the saturated zone. As illustrated by the flow chart (Fig. 2) the following flow phenomena are considered: infiltration, surface flow, drainage, evapotranspiration and groundwater flow. During a rainfall event, infiltration starts once the interception volume, which is -plant
model.
dependent, is satisfied. Infiltration is computed based on the Green-Ampt model. At the end of the rainstorm, soil water redistribution is determined based on the assumption of a gravity flow (head gradient equal unity) under unsaturated conditions. The evapotranspiration, the recharge of the aquifer, which is the result of water percolating out of the vadose zone, and the groundwater discharge are computed daily. A daily water balance is maintained. More details will be given about the drainage and groundwater components because they have been modified from the original ANSWERS model.
3. Model overview At the surface, the net rainfall (i.e. rainfall minus interception) is subject to infiltration and runoff for each time step and each block. Infiltration starts once the interception volume is filled. A maximum intercepted volume, dependent on the plant type, is interpolated linearly according to the plant stage of growth. Plant growth is represented by a varying leaf area index (LAI) and root zone. For each crop cover, tabulated values of the LA1 (Knisel, 1980) are used as input to the model for 10 stages of the plant growth. A linear interpolation is made on a daily base for each stage. The root depth, which determines the extent of the root zone, is computed using a sinus function related to the number of days after planting and the number of days to reach maturity (Borg and Williams, 1986). The intercepted water is then returned to the atmosphere through evaporation. Evapotranspiration is determined based on the approach of Ritchie (1972). The potential evapotranspiration is computed with the Priestly-Taylor equation (Priestly and Taylor, 1972) using the net solar radiation and mean air temperature. The daily value of LA1 is used to partition the potential evapotranspiration into potential soil evaporation and potential plant transpiration, the latter term being limited by the amount of water in the root zone. Soil evaporation is assumed to take place in two different stages. In the first stage (constant rate stage), soil evaporation is just limited by the amount of energy available and proceeds at the potential rate. Beyond an upper limit depending of the soil water transmission characteristics the second stage begins. During this stage, the falling rate stage, soil evaporation rate decreases as a square root function of time. For any given day, the sum of plant transpiration and soil evaporation cannot exceed potential evapotranspiration. Infiltration is simulated using the physically-based Green-Ampt equation (Green and Ampt, 1911). It assumes (i) a step water retention function h(8) (wetting front capillary potential) describing the relation between volumetric water content, 19(cm’/cm”), and soil water pressure, h (cm); and (ii) a one-point hydraulic conductivity function K, (cm/h) where K, is the hydraulic conductivity for volumetric water content at natural saturation 8,. The infiltration process is represented as a wetting front traveling down
the soil protile with 0 = 8, and K = K, behind rhe wetting front and 0 = 0,) and K = 0 ahead of the wetting front. The saturated hydraulic conductivity and the wetting front capillary potential are determined from textural soil properties using statistical regression functions (pedotransfer functions) given by Rawls and Brakensiek ( 1989). Once the infiltration is determined, rainfall excess is routed to the watershed outlet. A modified discharge-water depth equation of the Manning type is used for both overland and channel flow routing. Every square element of the discretized watershed acts as an overland flow plane with a user specified slope and slope direction. Channel elements (described in terms of their slope, width and roughness coefficient) collect flow from overland flow elements and route the runoff to the watershed outlet. The drainage in the unsaturated zone is determined based on Darcy’s approach on the assumption of a gravity flow under unsaturated conditions. Travel time of percolating water through the soil matrix is regulated by the unsaturated hydraulic conductivity. This conductivity K(0) is modeled using the Brooks and Corey ( 1964) equation: / a\‘i (1) where 11represents the conductivity shape parameter, 0, (cmj/cm”), the water content at natural saturation. both determined from pedotransfer functions (Rawls and Brakensiek, 1989). Travel time [ 7T, (h)] through a particular soil layer is computed using a linear storage equation:
where d, is the thickness of the specific layer (cm). Drainage DR (cm) during a specific time step, At (h) is determined using an exponential function (Williams et al., 1985):
1 The water percolating out of the last unsaturated layer (Fig. 1) is added to the amount of water present in the groundwater cell and becomes available for lateral saturated flow. Because of the high infiltrating properties of the soils in the study area, it is assumed that there is no
F. Bouraoui et al./Joumal
Table 1 Major soil parameters
of Hydrology 203 (1997)
83
79-92
used in the model Layer O-30 cm
Bulk density (g/cm’) Clay fraction (%)** Sand fraction (%)** Saturated conductivity (cm/h) Wetting front capillary potential, ht (cm) Saturation water content, 6, (%)*** Conductivity shape parameter, n Volume of coarse fragments, VP (%)**** Fine fraction. or (%)*
1.339 17.5 41.0 1.19 26.1 90 19.6 25.2 60.0
Layer 30- 120 cm 1.316
16.4 51.7 3.44 27.3 90 13.4 49.4 33.7
* percentage of the total weight. ** percentage of total weight of fine particles. *** percentage of total porosity. **** percentage of total volume.
lateral subsurface flow from one cell to another in the unsaturated zone. Thus, the unique recharge of an aquifer cell is the excess water draining from the corresponding overlying cells. The draining water is added to the groundwater volume present in each aquifer cell. The total groundwater volume is then divided by the specific storage to obtain the piezometric head for each cell (Hi). An average daily value of aquifer discharge, Q (m3/ day), is computed at the downstream limit using a Darcy-type equation, such as: Q-m:
(4)
where K represents the conductivity of the aquifer material (m/day), M/AZ is the average hydraulic gradient (m/m) between the upstream and downstream limit of the aquifer, and A is the wetted cross section of the aquifer (m*). The resulting total outflow obtained from Eq. (4) is distributed to each cell and a new piezometric head is determined for each cell.
4. Experimental layout and data description 4.1. The site The site selected for the calibration and validation is ‘La C&e St Andre’, 60 km northeast of Grenoble (southeast of France). The choice of this site is related
to important threads concerning the qualitative and quantitative degradation of the groundwater resources with the development of intensive agricultural practices. It should be noted that the groundwater is the only drinking water supply for the local collectivities; its evolution is a source of conflict between farmers and cities. For this reason, an intensive, pluriannual, interdisciplinary study was initiated in 1991 in order to optimize agricultural practices, and to develop sustainable management schemes. Special attention was given to the characterization of the relationships between fertilization, irrigation and crop production, in particular to the water and nutrient balance in the soil-plant-atmosphere continuum (Kengni et al., 1994). Another important objective of the study is to predict the evolution with time of the groundwater resource. Nitrate is the major nonpoint source groundwater pollutant with nearly one half of the wells reaching, or overpassing the European limit of drinkability (50 ppm). 4.2. The scales of investigation Three different scales have been intensively investigated: local, field and watershed. At the local scale, measurements were carried out on an undisturbed lysimeter (1.2 m diameter, 1.5 m height), with bare soil or grass coverage, in order to obtain direct measurements of drainage and evapotranspiration and to characterize in detail the hydrodynamic transport properties in the soil. The lysimeter was
Table 2 Monthly rainfall and potential
evapotranspiration
1991
February March April May June July August September October November December
for the study site I992
I993
19V4
Rain
PET
Rain
PET
Rain
PET
Rain
18.1 26.6 118.4 51.6 21.6 92.8 49.6 21.9 94.4 133.6 53.3 93
18.3 3 I .6 49.1 17.3 106.8 103.4 163.2 153.5 92. I 41.4 23.5 20
25.8 37.4 75.9 57.9 96.9 120.6 99.7 51.5 93.3 182.1 114.1 64.1
18.4 29. I 50 77.9 106.7 88.6 127.8 141.6 76 28. I 19.2 21.4
2X.8 9.5 8.7 110.2 104.1 86.6 129.7 36.2 258.3 250.7 26.9 83.4
26.6 24 59.9 73.7 105.9 126.8 129.6 122.8 64.5 41.5 23.5 24.4
126.4 26.9 15.3 126.5 88.2 45.4 I7 37.9 192.8 80.4 107.9 30.5
25.7 61.5 55.4 95.1 129.4 170. I 156.5 ?4.3 49..3 33 ._’ -_ 14.9
780.9
880.2
823.2
8V5.2
886
1019.3
instrumented with neutron probe access tubes, tensiometers and soil temperature gauges, and direct outflow was automatically read. At the field scale, experiments were carried out on two types of soil cover: maize and bare soil. Measuring devices identical to those installed in the lysimeter were set up on two plots, one for each treatment, to determine the water balance in the root zone (Kengni et al., 1994). At the watershed scale, three different sets of data are looked upon: the soil, the soil cover (crop cover and management practices) and the aquifer (boundaries, transmissivities, and piezometric level and water quality variations). The watershed is about 320 km2 and the elevation varies from 480 to 250 m. The major part of the watershed (about 80%) is a flat plain with very permeable, fertile, unconsolidated soil which is surrounded by hills, with an average slope of 7.7%. Hills are covered with grass and forest; the plain is cultivated either with irrigated maize and tobacco (20% of the area) or with dry farming crops, such as sunflower and corn. In the plain, the arable layer is shallow (average 1 m) and rich in organic matter. There is almost no hydrographic network due to the high infiltrability of the arable layer which consists of a sandy loam becoming progressively coarser with increasing depths. Roughly speaking, the soil is characterized by two different layers (Table 1).
784.8
1133.1
PET 2 I .6
The groundwater holding layer consists of a very coarse glacial material extending from a depth of 1 m to about 20-30 m, with a water table aquifer. Due to the presence of coarse fragments, it was not possible to install any measuring device to characterize transport in the unsaturated part of this glacial material. The aquifer is characterized by very high transmissivities (range from 5 x lo-* to 10-l m*/s), high velocities (between 5 and 10 m/day). Consequently, the average discharge is approximately 4 m3/s. The effective porosity of the aquifer is about lo%, and its depth fluctuates around 5 m during a year, with a typical range of average depth of 10 - 15 m below the soil surface. The changes of piezometric levels may be very fast during the rain season. The average annual rainfall over 30 y is around 1000 mm, with two characteristics: a very high interannual variability and an annual pattern with the most important rain events always in September-October, and another rain season, weaker and more uncertain, in April-May (Table 2). Mean evapotranspiration computed by the Penman-Montieth approach is around 850 mm/y. 4.3. Available
data
In the lysimeter and on the bare and corn plots, soil water pressure, and volumetric soil water content were measured at different depths in the upper soil
F. Bouraoui et al./Journal
layer (thickness 1 m) during the entire crop season, at a time step varying from 30 min to one week. Flux draining out of the lysimeter were measured continuously, whereas on the field sites, the fluxes below the root zone were estimated on a daily basis by the use of Darcy’s law following a method described previously (Kengni et al., 1994), and resting on time course measurement of soil water content, soil water pressure, and direct estimation of the K(8) relationship. No measurement of the flux and soil transport parameters could be done in the underlying layer because of the presence of coarse gravel and pebbles. The evaporation on the bare soil plot and the evapotranspiration on the grass lysimeter and maize plot were determined by mass balance as the difference between the estimated drainage and the changes in soil water storage. All relevant micrometeorological data (rainfall, radiation, pan evaporation, etc.) were measured, at a time step of 30 min, during the entire experiment. Finally, at the watershed scale, the actual availability of data concerns the following items: in term of the aquifer, an intensive geophysical series of measurements (seismic and electric) were made 10 years ago to determine the geometry of the bedrock; in terms of soil, a survey was conducted by the Regional Chamber of Agriculture resulting in a database including texture, organic matter, and density; concerning the soil coverage and management, a map was made by the Institut Supkrieur d’Agriculture Rh6ne Alpes (ISARA, Lyon) following a systematic survey of farming practices. Continuous measurements of the piezometric level (on the basis of weekly values) are available for three years for five different wells.
of Hydrology 203 (1997) 79-92
X5
measured data mean is a better predictor than the predicted data). The mean and standard deviation are also reported. Care should be given, however, when interpreting these values since the measurement periods can vary from plot to plot (no crosscomparison is possible for the different types of cover). 5.2. Grid scale validation The model was first used at a grid scale to obtain a validation between the model output and the measurement of water flux for different surface cover. This internal validation at the grid scale is an essential step before conducting simulations at the watershed scale in order to ensure the proper representation by the model of the different components of the water balance for different surface covers. It was assumed that within the grid, the soil and its cover have the same characteristics as those determined on the measurement sites. Three different applications were carried out: first on the bare soil lysimeter (later planted with grass) with comparison between simulated and measured drainage at the lysimeter outflow and real evapotranspiration; then on the bare and corn plots with comparison between simulated and estimated drainage. On the lysimeter about 18 months of measured drainage was available. The lysimeter was installed in September 1992 and was kept bare until June 1993 after which it ‘was planted with grass. With 100
g
E I0 E
5. Validation 5.1. Statistical
6 zii analysis
1
s 0
5 The accuracy of the model predictions were analyzed using a correlation coefficient which gives a measure of linear association between the predicted and measured data, and a coefficient of efficiency which represents the fraction of total variance of the measured data explained by the predicted data (a negative value of this coefficient indicates that the
0.1
E
u
I” 0.01 0.15
0.25
0.35
0.45
Water content cm3km3 Fig. 3. Hydraulic conductivity curve for the lysimeter.
Table 3 Measured
and simulated
R*
Cover type
Year
Bare soil
91-92 92-93 93-94 91-92 92-93 93-94 92-93
Corn
Lysimeter
actual evapotranspiration
640 909
II00 640 909 II00 1394
(AET, mm) and drainage
(D, mm) for different
surface cover
PET*
AET,
AET,
AETJAET,
D,
D\
DJD
726 595 514 726 595 574 962
320 392 379 495 514 548 609
355 345 330 580 524 513 663
1.109 0.880 0.87 I 1.172 I.019 0.936 1.089
516 572 753 342 467 551 908
so4 560 762 280 478 559 955
0.977 0.979 I.(112 0.819 1.024 1.015
,n
I .os1
-
* R: rainfall
(mm). ** PET: potential evapotranspiration m measured. s simulated.
(mm).
regard to the physical properties of the soil, the soil profile was discretized in two horizons: O-O.3 m and 0.3- 1.2 m. The values of K, and hr are computed from the pedotransfer functions developed by Rawls and Brakensiek (1989) and are given in Table 2, along with the other parameters necessary to run the model. The saturated hydraulic conductivity and the porosity are adjusted for the presence of coarse fragments according to the procedure described in Appendix A. A comparison between the values predicted by the pedotransfer function and the measured K(0) is given in Fig. 3. It must be noted that due to specificities of the ‘zero flux’ method, it is not possible to obtain measured values of K(8) in the wet range. The results for cumulative model predictions for drainage and real evapotranspiration (AET),
Table 4 Statistical
without any calibration are very satisfactory (Fig. 4). Total drainage and AET are well reproduced (Table 3) and are within 10% of the measured data. The performance of the model is quite acceptable considering that no adjustments were made (Table 4). Identical soil characteristics were used in the simulation of the corn and the bare plots. Fig. 5 represents the measured and the simulated cumulated drainage, and the estimated and simulated AET, respectively, for the bare plot for three years of data. The predictions of drainage and evaporation were good for the three years of simulation. Similar results are also obtained on the corn plot (Fig. 6). For both surface covers, the cumulative predictions are within 20% of the measured data (Table 3). The predictions of drainage are better than those of AET and this
analysis for all simulations
Actual evapotranspiration Baresoil Corn lysimeter Bare soil Corn Lysimeter Well
11.48 17.28 16.46 Drainage 19.27 15.12 15.67 Piezometric 342.29
* correlation coefficient. ** coefficient of efficiency m measured. s simulated.
10.70 17.96 17.94
8.52 13.98 1I .50
5.80 9.88 9.54
0.667 0.741 0.834
0.439 0.553 0.683
23.14 14.61 14.88
26.82 26.82 27.82
25.85 25.85 30.48
0.933 0.916 0.881
0.873 0.837 0.731
342.46
I.12
I .08
0.957
0.892
head
.
0001 I
..
'(pa1qIz,Q.m03
doJ3 aql Bu!%ueqs l~ql alexpu! satyeA paJnseam ~(ai!pm ‘SWJ% 'pas mq)
aaql
~03 AlaleJnm? apnb
J1AO3 JO SadB luaJaJJ!p pauuoyad
11
-sacylxzJd
JO sadlClaaJylaqlJojL~vJoJ wql ah.qmp
qmu
a~ slInsaJaqL -clappaJnwaur aq
_103~awq LIalenbape
lwld aql 30 uopmy!~dLu~sJaAo aql 01 anp aq hxu
lCJaAsf Iapow aq3 leqlmoqs sllnsaJsno!AaJdaqJ
pampoJdaJ IapouIaql ‘aseqLIyaaM B ug -1uauodwos
luama%euk?wpu~JaAo~doJ3u~ sa%.mq301 aA!suodsaJ
.Ja]aur!s/([ aq) 104 13~
pm a%eu!wp pawpaJd SA pamsea~
‘p +?!,I
e&u!wp pew~nw!s S&“!EJp
pamsew
??
I i ooz Z6J@HI -t----f 3 - ooz -8 ii 3
oat
3 d s
18
F. Bouraoui et d./Journulof Hwirology 203 (IYY7) 70-02
88
600 1
-
. measured AET simulated AET
4oot
E-
s fk e .e
. 600
b
600
1
-
measured drainage simulated drainage
D
1000
Fig. 6. Measured
vs predicted
cover from maize to bare conditions, results into an increase of drainage of 44%, 14% and 28% for the years 1991, 1992, 1993, respectively. The model estimated these increases at 33%, 19% and 26%. Finally, the model predicted well the fact that most drainage does not occur during the growing season of the plants, even during irrigation, but during the winter, when the soil is bare and the values of AET are very low. 5.3. Watershed
validation
The next step was the application of the model to the entire plain to test its ability to predict the recharge of the aquifer as influenced by the surface cover management and soil depth to the aquifer. The soil was discretized into three layers: two layers in the root zone, and a third one which extends from below the root zone to the top of the aquifer. It is assumed that all the water leaving the root zone is transferred in its entirety to the aquifer. The surrounding boundary conditions for the aquifer were assumed to be a no-flow boundary condition as suggested from previous hydrogeological studies: the unique recharge comes from atmospheric water. Fig. 7(a) illustrates the elevation map of the underlying impervious layer. This
drainage
and AET for the corn plot
map of the substratum was interpolated using well profile information. Fig. 7(b) provides the piezometric map for the whole watershed for 7 January 1991. The piezometric map was kriged using 78 data points spread throughout the region; additional details about the kriging procedure can be found in Ruy (1992). As mentioned earlier, three years of continuous measurements of water level at five different wells, one of them (Feuges) located near the study site in the heart of the irrigated area, are available. The vegetation parameters, and the water input (rain or rain + irrigation) were distributed over the plain according to Fig. 7(c). The first six months of the data set (Feuges well) were used for calibrating the transfer parameter in the unsaturated zone extending from below the root zone to the water table (third layer). The calibration consisted in adjusting the hydraulic conductivity until the predicted and measured value of the piezometric head at the Feuges well were in agreement in terms of timing of the rise of the water level. The transfer parameters in the root zone (first and second layers) were determined from the pedotransfer functions developed by Rawls and Brakensiek (1989). The model was then run for 30 months in a purely predictive way. The prediction for the water level for the
x9
F. Bouraoui et alJJourna1 of Hydrology 203 (1997) 79-92
a) IMPERVIOUS
LAYER
ELEVATION
(m)
Feuges well
b) PIEZOMETRIC
c) LAND SURFACE
Fig. 7. Substratum depth (a), piezometric (above sea level).
HEAD
AND CROP
(m)
COVER
head (b), land surface and crop cover for the plain of Bievre. Isovalue lines represent elevation in m
‘Feuges’ well is given in Fig. 8. The piezometric head is expressed in meters above sea level, the soil surface elevation is 356.10 m. The model reproduced well the variations of the water table level (Table 4), in particular the high level of water during winter and early spring when the soil is completely bare (Fig. 8). The model also predicted that hardly any significant rise of the aquifer level will occur during the summer. This is mainly due to the presence of maize, a highly transpirating plant which limits the amount of deep percolation. Thus, during these periods of no deep percolation, there
will be no risk of nitrate leaching. This has already been noticed by Kengni et al. (1994). The model calibration leads to an estimate of travel time through the vadose zone of about three weeks for heavy rain periods. This value indicates how quick the aquifer responds to the excess rainfall. Finally, an increase of the water table level from the beginning to the end of the simulation period (3 y) by more than 4 m is properly predicted. This indicates that the model is very responsive to large variation in climatic input. The results for the Feuges well were confirmed by the four additional wells (Fig. 9).
90
measured 340
339 _ _.~_ Nov.90
Fig. 8. Predicted
-
1
~~_ ~~~ em- _~_ Jun-91
YS measured
Deco1
piezometric
~~
Jul-92
~
Jan-93
-
simulated
~-. Mar-94
Aug-93
level for the ‘Feugea’ well (soil surface elevation
6. Summary and conclusions
layers, and excess water from the root zone is routed to the aquifer. A Darcian-type equation was used to estimate groundwater discharge. Local measurements on a bare soil lysimeter, a bare and a corn plot were used to validate the model in the root zone over three years. Since the model is intended to be used for predictive purposes on watersheds where measured information is limited or lacking, indirect methods of estimation of soil transfer properties (pedotransfer function) and tabulated values of plant growth were
A continuous, distributed nonpoint source model for estimating groundwater recharge as affected by the spatial distribution of soil property characteristics, spatial and temporal variability of management practices, was developed. The model simulates infiltration by the Green-Ampt method. Ritchie’s approach is used to estimate plant transpiration and soil evaporation separately. The soil is discritized into different
- well . well 0 well 0 well . well
1
2 3 4 6
-+-_ ._~~ 336
340
344
346
Measured piezometrlc head (m)
Fig. 9. Predicted
350.10 m)
vs measured piezometric
level for the five validation
wells
91
F. Bouraoui et al.Nournal of Hydrology 203 (1997) 79-92
privileged over the measured characteristics. The information obtained at the local scale was then used to run the model on a large watershed. The model performed quite well in reproducing the piezometric head values and variations at several wells. However, it is important to note that the research described above showed that a simple model, with little input data required, can be a powerful tool for predicting and managing natural water resources. It may be added that the proper functioning of the model certainly results from a good agreement between its main hypothesis (infiltration governed by Green and Ampt assumption) and the type of soil characterizing the watershed. Questions remain open about the efficiency of this type of model to predict groundwater recharge in clayey material.
Acknowledgements This study was realized with the joint support of a European Union Program (EV5V-CT94-0484, ‘Evaluation of the Effect of Climatic Variations on the Recharge of Aquifers in Southern European Catchments’, DGXII) and of the Centre National pour la Recherche Scientifique, CNRS, Paris, Programme Environnement, Vie et SociCtC.
Appendix
A Adjustment
for Coarse Fragments
Most of the measurement techniques give an estimate of the whole-soil volumetric water content (e,,, cm ‘/cm 3). The representative x/olume for which these water content values are determined may include not only fine soil but also coarse fragments and stones, when present (as in our case study). The terms ‘coarse fragments’ and ‘stones’ as used in this paper concern mineral fragments more than 0.2 cm in diameter regardless of their specific size and shapes (Soil Survey Laboratory Staff, 1992). The whole-soil bulk density, pap (g/cm’), applies to the total sample volume [VT (cm’)] for which 8,, is determined. The usual bulk density, pf (g/cm’), refers to the volume of fine soil, Vf (cm3), and the specific density pc (g/cm’) to the volume fraction of coarse fragments, V, (cm”) which are linked through:
Writing
the conservation
M,+Mf=M,
(Al)
(A2)
where M,, Mr and MT represent the mass of coarse fragments, the mass of fine material and the total mass, respectively. Eq. (Al) can be rewritten as: @
Mf
(A3)
p_+-=vT ‘ Pf
The combination of Eq. (A2) and Eq. (A3) allows the calculation of the fine soil bulk density as: Pf =
‘yf Mw --___
(A4)
(l-af)
&i,MT
PC
where lyf is the fraction by weight of the fine material and A4w is the mass of soil water of the whole-soil sample determined through the process of weighing, drying and reweighing of the sample. Taking pc equal to 2.65 g/cm” and pr being known from Eq. (A4), the volume fraction of coarse fragments, Vi’ defined as the ratio V, and VT becomes: (1 -at.)
v,“=
l+cr. t
(AS)
p__,( Pt >
The volumetric tine soil water content and/or soil porosity is then adjusted by a correction factor equal to (1 - VP) and the saturated hydraulic conductivity is adjusted by cxf. Additional details about this correction procedure can be found in Haverkamp et al. (1997).
References Abbott,
M.B.,
Bathurst,
J.C..
Cunge.
J.A..
Connell,
P.E..
Rasmussen, .I.. 1986. An introduction to the European Hydrological System (Systbme Hydrologique structure of a physically-based,
Europeen.
‘SHE’,
2:
distributed modeling syctem).
Journal of Hydrology, 87, 6 I-77. Borg. H., Williams,
D.W.,
1986. Depth development of roots with
time: An empirical description, Transactions of the ASAE. (I).
29
194-197.
Beven,
K.J.,
1989.
Changing
ideas in hydrology:
physically based models, Journal of Hydrology. Bouraoui, F., Dillaha. T.A.,
1994. ANSWERS-2000.
the case of I OS, I S7- 172. Continuous
simulation version. ASAE Paper No. 94-2 I SS. presented at the ASAE
v, + v, = vr
of mass gives:
USA.
International
Winter
Meeting,
Kansas City. Miwourt.
Y2
I,‘. Rouruoui
rt ul.Nournul
Bouraoui. F.. Dillaha, T.A.. 1996. ANSWERS-2000: Runoff and sediment transport model. Journal of Env,ironmental Engineering, ASCE, 122 (6), 493-502. Brooks, R.H., Corey, C.T., 1964. Hydraulic properties of porous media. Hydrology Paper 3. Colorado State University, Fort Collins, USA. Green, W.H., Ampt, G., I91 I. Studies of soil physics, Part 1. The flow of air and water through soils, Journal of Agricultural Science, 4, l-24. Haverkamp, R., Arrue, J.L., Ross, P.J., Dalton. F.N., Lopez. Ma.V.. Bouraoui, F., Garcia. R., 1997. Influence of stones and porous coarse fragments on water content measurements and bulk density at field scale. Soil Science Society of America Journal (submitted). Kengni, L., Vachaud, G., Thony, J.L., Laty. R., Garino. B.. Casabianca, H., Jame. P., Viscogliosi, R., 1994. Field measurements of water and nitrogen losses under irrigated maize. Journal of Hydrology, 162, 33-46. Knisel, W.G., 1980. CREAMS: A Field Scale Model for Chemical, Runoff, and Erosion from Agricultural Management Systems. Conservation Research Report No.26, USDA-SEA.
of Hvdrolo,cv 203 (lYY7) 7Y--92 Priestly. C.H.B, Taylor, R.J., 1972. On the assessment 01 rurtacl, heat flux and evaporation using large-scale parametric. Monthly Weather Review 100, 8 I-92. Rawls. W.J.. Brakenstek, D.L.. 1989. Estimation of soil water retention and hydraulic properties. In: Morel-Seytoux. H.J. (Ed.). Unsaturated flow in Hydrologic Modeling: Theory and Practice. NATO ASI Series c: Mathematical and Physical Sciences, Vol.275. Kluwer Academic, Dordrecht. Holland. Ruy. S.. 1992. Etude de la nappe phreatique de la plaine de Ia Bievre. DEA. Laboratoire d’Etude des Transferts en Hydrologie et Environnement. Grenoble, France. Ritchie, J.T.. 1972. A model for predicting evapotranspiratton from a row crop with incomplete cover, Water Resources Research. 8 (51, 1204-1213. Singh, V.J.. 1995. Computer models of watershed hydrology. Water Resources Publications, Highlands Ranch, CO. Soil Survey Laboratory Staff, lY92. Soil survey laboratory methods manual. Soil Surv. Invest. Rep. 42 Version 2.0 IJSDA-SCS U.S. Gov. Print. Oftice. Washington, D.C., USA. Williams. J.R.. Nicks, A.D.. Arnold. J.G.. 1985. Simulator of water resources in rural ba.:ins. J. Hydr. Div., ASCE, I I I (6),970-9X6.