- Email: [email protected]
Consequences of preferential flow in cracking clay soils for contamination-risk of shallow aquifers
Journal of Environmental Management (1995) 43, 359-373
Consequences of Preferential Flow in Cracking Clay Soils for Contamination-risk of Shallow Aquifers K. Oostindie and J. J. B. Bronswijk DLO Winand Staring Centre for Integrated Land, Soil and Water Research (SC-DLO), P.O. Box 125, 6700 AC Wageningen, Netherlands Received 3 March 1994
A simple method is presented to assess, on a regional scale, the contaminationrisk of shallow aquifers covered by cracking clay soils, with special emphasis on preferential flow of contaminants through shrinkage cracks. A water extraction area of several hundred hectares is divided into units with homogeneous soil types and hydrological conditions. For each unit, a onedimensional computer simulation model is used to compute residence times in the unsaturated zone of the cracking clay soils. Analytical equations are used to estimate residence times in the saturated zone. Finally, a classification system is proposed to summarize and visualize the model computations in a simple way. The developed method is illustrated with a case study from the Netherlands. It is known how contamination-risk maps of water extraction areas can be obtained. Policies to prevent contamination of the aquifer can be focused on the most vulnerable areas as indicated on such a map. The contamination-risk maps, obtained using the present method, deviate sharply from the traditional, more or less concentric protection zones, computed with saturated groundwater flow models.
Keywords: aquifer contamination, clay soil, shrinkage crack, preferential flow.
1. Introduction A saturated clay has a low hydraulic conductivity. As a result, the residence time of water flowing through saturated clay layers is long. Because of these long residence times and the generally high adsorption capacity of clay, contaminants will often be adsorbed, or transformed by chemical or microbiological processes, before passing a clay layer. Therefore, saturated clay soils overlying an aquifer effectively protect that aquifer against contamination originating from the soil surface. When saturated clay soils dry, however, their volume decreases with the result that there is subsidence of the soil surface and the appearance of shrinkage cracks. Water and solutes can be transported rapidly to the subsoil and the groundwater by preferential flow via these cracks (e.g. Beven and Germann, 1982; Bouma, 1984). In shallow aquifers, preferential flow in a cracked clay soil may lead to very short residence times for solutes from the soil surface to the aquifer. 359 0301-4797/95/040359+ 15 $08.00/0
© 1995 Academic Press Limited
360
Preferential flow and aquifer protection
At many locations in the Netherlands, pumping stations withdraw water from aquifers for drinking water. In the water extraction areas surrounding these pumping stations, special regulations are enforced to prevent contamination of the aquifer. Water extraction areas are divided into zones, depending on the estimated average travel times of water from the soil surface to the extraction point of water from the aquifer. Generally, one-year, 10-year and 25-year-zones are distinguished. In the one-year-zones, very strict regulations, for instance with respect to the use of pesticides by farmers, are enforced. For the 10- and 25-year-zones, less severe regulations apply. Several water extraction areas are located in the riverine area in the central part of the Netherlands. These water extraction areas have a similar hydrogeological profile: a Holocene heavy clay soil overlies an aquifer from which drinking water is extracted. In some cases, the clay deposits are rather thin (less than one metre). In such situations, preferential flow of solutes via shrinkage cracks to the subsoil and the aquifer may be a serious hazard for the water quality in the aquifer. To enable policy makers to develop scientifically based regulations for land use in the vicinity of drinking water extraction points in clay soil areas, it is essential to be aware of the risk of preferential flow on a regional scale, and its consequences for residence times of water and solutes in the clay soils overlying the aquifer. The transport of water through cracking clay soils is, however, a complex process. Qualitatively, it is determined by the following factors (Bronswijk, 1991): •
•
•
•
•
•
Soil properties The more severely the clay cracks, the more water will infiltrate into the cracks, and the greater the preferential flow will be. The higher the hydraulic conductivity of the soil matrix, the shorter the residence times will be. Precipitation surplus When the precipitation surplus is high, in general more water will finally reach the aquifer. Precipitation intensity In the case of showers with a high precipitation intensity, preferential flow will be high. Drainage In an intensively tile-drained soil, a considerable amount of the infiltrating water will flow via drains to surface waters, and less water will reach the aquifer. On the other hand, intensive drainage results in deeper water-tables and deeper cracks, leading to shorter residence times. Upward and downward seepage to or from the aquifer Downward seepage enhances downward flow and decreases residence time from soil surface to aquifer. For upward seepage, the opposite is true. The thickness of the clay layer Aquifers covered by thin clay layers are more vulnerable than those covered by thicker ones.
To compute residence times of water in the covering clay layers, it is necessary to consider all these factors together with their complex interactions, and the use of computer simulation models is indispensable. Because solute transport from soil surface to the aquifer is mainly convective, only the transport of water needs to be considered for a first estimation of residence times and contamination-risk. The objective of this paper is to present a simple method to assess the contaminationrisk of shallow aquifers covered by cracking clay soils on a regional scale, with special
K. Oostindie and J. J. B. Bronswijk
361
emphasis on preferential flow through shrinkage cracks. A water extraction area of a few hundred hectares is divided into units with homogeneous soil types and hydrological conditions. For each unit, a one-dimensional computer simulation model is used to compute residence times in the unsaturated zone of the cracking clay soils. Analytical equations are used to estimate residence times in the saturated zone. Finally, a classification system is proposed to summarize and visualize the model computations in a simple way. The developed method will be illustrated with a case study from the Netherlands. 2. Materials and methods
2.1.
BASIC PRINCIPLES OF THE METHOD
The developed method consists of six steps: 1. The boundaries of the area under consideration are defined. In the present study, we have used the 25-year-zones around the pumping stations. These are the ones where the residence time of water towards the extraction filter of the pumping station is less than 25-years, as computed with steady state groundwater flow models. 2. A soil survey is conducted on a 1:50000 scale to identify the soil types in the considered area. Subsequently, soil physical properties of the various soil types are collected. 3. The hydrological and climatic data, required for the model simulations, are collected. Areas with homogeneous hydrological conditions, such as drain depth, are grouped into hydrological units. For climatic data, an average year, a wet year, and a dry year, with respect to precipitation and evaporation, are selected. 4. Combination of the soil map and the hydrological map yields a simulation map, in which each unit is considered homogeneous with respect to: --soil profile and soil physical properties; ----drainage conditions; ----depth of the aquifer; --hydraulic head in the aquifer. 5. For each simulation unit, residence times are then computed for an average year, a wet year, and a dry year, using the collected input data in computer simulation models and analytical equations. 6. Finally, the computed residence times for the average, the dry and the wet year, are classified, and a contamination-risk map of the considered area is constructed. 2.2.
COMPUTING RESIDENCE TIMES IN THE UNSATURATED ZONE
Figure 1 shows the schematization as applied in the present research. Water transport in the unsaturated zone is calculated with the FLOCR computer simulation model (Bronswijk, 1988; Oostindie and Bronswijk, 1992). FLOCR computes unsaturated vertical water flow in cracking clay soils, taking into account swelling and shrinkage processes and preferential flow through cracks. The soil is divided into soil matrix and shrinkage cracks. Precipitation is divided between matrix and cracks. Water flow in the soil matrix is very slow. Water infiltrating into the cracks is assumed to reach the bottom of these cracks instantaneously, without lateral infiltration into crack walls. Due to limited contact area between water and crack walls, and weak diffusivity forces,
Preferential flow and aquifer protection
362 Unsaturated zone
Watertable
Saturated zone
Clay
ii ................................................................
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii i iiiiiiiiiiiiiiiiiiiiiiiii Figure 1. Schematization of the soil profile in riverine water extraction areas.
Variable crack volume
I E
F
Subsidence o
lI 2 3 4
27 28
29 3O ......
Matrix-crack system at time t Matrix-crack system at time t + At
Figure 2. Principles of the FLOCR model, for computing water balance, cracking and subsidence in clay soils./, infiltration rate in soil matrix (m s-I); let, part of total crack infiltration caused by rainfall intensity exceeding maximum infiltration rate of soil matrix (m s-~); Ic2, part of total crack infiltration caused by rainfall directly into the cracks (m s-t); E, actual evapotranspiration (m s-~); V, Darcy flux between two nodal points (m s-a); V,aL,,drain discharge (m s-I).
lateral infiltration has been f o u n d to be only a m i n o r fraction o f the soil water balance ( H o o g m o e d and B o u m a , 1980; Bronswijk, 1988; Booltink and B o u m a , 1993). The principles o f F L O C R are pictured schematically in Figure 2. F L O C R produces o u t p u t per one-fifth o f a day.
K. Oostindie and J. J. B. Bronswijk
363
Water in the cracks and water in the matrix eventually reach the groundwater. There, some of the water flows, via the subsoil, towards drains and ditches, and some will reach the aquifer. For water flowing through the cracks, the residence time in the unsaturated zone is set at zero days: to=0
(1)
where t~ is the residence time of water in the cracks (d). For water inside the soil matrix, the residence time in the unsaturated zone is approximated by (e.g. Van Lanen, 1984):
tu= ~Ps
(2)
where tu is the residence time of water in the matrix (d), Tu is the average thickness of the unsaturated zone (cm), ® is the volumetric water content, and Ps is the average precipitation surplus (cm d-t). Equation 2 yields only a crude estimation of the residence time in the unsaturated zone of the soil matrix. Due to upward and downward movements of the water in the root zone, and to subsequent precipitation and evapotranspiration, the residence time in the soil matrix, as computed using Equation 2, must be considered as average. This rather crude approximation is justifiable because of the much lower residence times of matrix flow in comparison with crack flow in clay soils. As a result, contamination risk is determined much more by crack flow than by matrix flow. 2.3.
COMPUTING RESIDENCE TIMES IN THE SATURATED ZONE
The Holocene layer in the majority of locations in the considered riverine area is a few metres thick, whereas water-tables are around l'50m or less. Therefore, in most situations, after reaching the groundwater table, a water particle still remains for some time in the clay layer, before reaching the aquifer. Because the clay is saturated here, the residence times will generally be high. The residence time in the saturated zone is calculated using the equation: t s - TsOs
(3)
v
in which ts is the residence time in the saturated zone (d), Ts is the thickness of the saturated zone (cm), v is the water flux through the saturated clay soil (cm d-l), and ®s is the saturated volumetric water content of the clay soil. The water flux (v) is the vertical flux between the water-table and the bottom of the clay layer which is calculated using Darcy's equation: v= -Ks AH
Az
(4)
where Ks is the saturated hydraulic conductivity of the 'clay layer (cm d-~), AH is the
Preferential flow and aquifer protection
364
50
500
40
400
I Crack infiltration I-1 Matrix infiltration
E E ,- 30
~E 300
0
g
.c_ 200 .x
-~- 20 U
10
100
0
0 100
250 Residencetime (d)
1000
Figure 3. Example o f a computed residence time distribution. Class width is 25 d.
difference between the water-table and hydraulic head of the aquifer (cm), and Az is the thickness of the saturated clay layer (cm). Finally, the total residence time of water in the clay layer, overlying the aquifer, is calculated using one of the following equations (see also Figure 1): 1.
For water in shrinkage cracks penetrating into the aquifer (Equation 1): t=0d
2.
For water in shrinkage cracks not penetrating into the aquifer (Equation 3): t =t,
3.
For water in the soil matrix (Equations (2) and (3)): t=tu+ts
2.4.
CLASSIFICATION OF RESIDENCE TIMES
By combining the FLOCR model and Equations (2), (3) and (4), a residence time is computed for each volume of water infiltrating at the soil surface. For a whole year, this computation yields a residence time distribution. An example is presented in Figure 3 for a hypothetical situation where the aquifer starts at a depth of 180 cm below the soil surface. For this situation, FLOCR, using the climatic conditions of 1973, has computed a total annual net infiltration (infiltration minus evapotranspiration) of 299 mm, subdivided in a matrix infiltration of 185 mm, and a crack infiltration of up to l l 4 r n m The 185mm of water infiltrating into the soil matrix has an average residence time of about 1000 d in the clay layer, before reaching the aquifer. Of the 114 rnm of water infiltrating into the cracks, 24 mm has a residence time of between 76 and 100 d in the clay layer, 39 mm has a residence time of between 101 and 125 d, 33 mm has a residence time of between 126 and 150 d, and so on.
K. Oostindie and J. J. B. Bronswilk
365
TABLE1. Classification of residence time distributions, based upon quantities of water reaching the aquifer and residence times of water in the covering clay layer. Quantity of water reaching the aquifer (mm a -I) >75 25-75 10-25 0-10
0-50 d
51-400 d
>400 d
1 1 2 3
3 4 4 5
5 5 5 5
The classes range from very vulnerable to contamination (class 1) to safe (class 5).
The contamination risk of the aquifer below this covering clay layer is now determined by the quantity of water that reaches the aquifer, and by the velocity with which the water reaches the aquifer. Therefore, a risk assessment based on the computed results presented in Figure 3 should consider both quantities of water reaching the aquifer and residence times inside the covering clay layer. To this end, the classification system of Table 1 has been developed. The five distinguished classes vary from very vulnerable (class 1) to safe (class 5). According to Table 1, the residence time distribution of Figure 3 should be assigned to class 3 because more than 75 mm water (in this case 114 mm) has a residence time of between 50 and 400 d in the clay layer. Using Table 1, the risk of aquifer contamination can be classified for each simulation unit, and made visible on a contamination-risk map. It should be emphasized that Table 1 is a subjective interpretation of the objective simulation results shown in Figure 3. An alternative classification system to the one developed and given in Table 1, would yield different categories. Although the objectivity of the model results is lost, a quick and easy interpretation of the data as presented in Figure 3 requires a classification system such as presented in Table 1. 2.5. LOCATION, SOILS AND MODEL INPUT PARAMETERS The developed method was applied to a number of water extraction areas in the Netherlands (Oostindie et al., 1990). In this article, the risk-assessment of the area of Druten is presented. The location of the Druten area is shown in Figure 4. The considered area was about 350 ha. Some 3 000 000 m 3 drinking water is extracted from the aquifer annually. In the Druten area, two soil units were distinguished: a heavy clay soil (Rn44C/ Rn47C on the 1:50 000 soil map of the Netherlands) and a clay loam soil (Rd90A/ Rn95A). The water retention curves and hydraulic conductivity curves for these soil types were taken from WOsten et al. (1987) and are shown in Figure 5. The depth of the aquifer was also determined in a field survey. The clay cover above the aquifer varied between zero and five metres. For upper boundary conditions, climatic data from the average year 1973, the wet year 1974, and the dry year 1976, were used. Precipitation and potential evapotranspiration are given in Table 2. 2.6. MODELVALIDATION Validating methods and models as applied in this study is extremely difficult. Such
366
Preferential flow and aquifer protection
Legend .,
1-year protection zone
0
2S-year protection zone
........
Provincial boundary
m.m
National border
0
10
20
I
I
I
30 km I
........."..
#.
~-°.T-° ,.,,'
.- .,-,., |
"'g cu,.~bo4
~
Verschuer~ ~
(~lmmerloo
~'~._.~, k;.~ ..,-'fx'-
Figure 4. Location of the investigated water extraction area of Druten.
K. Oostindie and J. J. B. Bronswilk
367
10 6
105 % "*...
lo 4 °
"0 t--
*o.o % °-°° % °.°o° % "..o.
lO 3 ~
a_
% °'°" ~N'°°°
lo 2
101
I 0.6
10 0 0.1
0.2
0.3
0.4
0.5
I 0.7
Volumetric water content
I""...
Topsoil Rn95A
1~ ~ k ~ A "7 "0
101 F ~ " , ~ . k
/
E
10 0 --
..............
TopsoilRn47C
------
SubsoilRn95A
"~'x "'.'~
t; "0 C
8
?.1
lO-,
"5 "1-
...... "-~--
10 -2
""-~. ..... o
t
I
10. 3 lO 0
lo ~
lo z
to ~
104
Pressure head (cm)
Figure 5. Water retention (a) and hydraulic conductivity (b) curves of the heavy clay soils in the Druten area.
validation would require long-term monitoring of water transport in the future on numerous locations, or historical data on water transport in the past. Neither options are realistic in practice. However, in previous studies, the applied FLOCR model has been validated with good results (e.g. Bronswijk, 1988). In the present study, we used existing data of minimum and maximum water-tables in the Druten area to validate our model results.
Preferential flow and aquifer protection
368
TABLE2. Climatic data used in the model computations of residence times for an average, a wet and a dry year Year
Precipitation (mm)
Potential evaporation (mm)
1993 (average) 1974 (wet) 1976 (dry)
779 993 536
559 581 646
3. Results and discussions
3.1.
WATER BALANCE AND RESIDENCE TIMES
Although two soil types have been distinguished, we will concentrate in the presentation of the results on the heavy clay soil area, because cracking and preferential flow are most pronounced in this soil type. The water balance of the heavy clay soils covering the shaUow aquifer, as computed with the FLOCR model, is pictured in Figure 6. In an average climatic year (1973), actual evapotranspiration is about 440mm and precipitation about 779mm. Of these 779mm, some 13% (104ram) flows via the shrinkage cracks as preferential flow to the subsoil, whereas the remaining 87% (675 ram) infiltrates into the soil matrix. Subtracting the actual evapotranspiration from these 675 mm gives a net matrix infiltration of about 235 mm. Because the grasslands in this area are poorly drained, the drain discharge is negligible. The net flux of water into the aquifer in an average climatic year is 296 mm per year. For a wet and a dry year the percentages of preferential flow of the total precipitation are 13% (126mm) and 14% (73 mm), respectively. The quantities of water flowing into the aquifer are 428 and 189 mm per year, for a wet and a dry year, respectively. For a simulation map unit with a heavy clay soil and an aquifer depth of 180 cm below the soil surface, the residence time distributions for the average, wet, and dry year, are shown in Figure 7. The residence times for water infiltrating into the soil matrix are about 800 d. For water flowing into the shrinkage cracks, the residence times are between 75 and 275 d. Such residence time distributions were derived for each simulation map unit. These residence time distributions have been classified using the classification system of Table 1. The present example would fall into class 3. A classification was obtained, using the procedure outlined above, for each of the simulation map units. The resulting contamination-risk map is shown in Figure 8. The heavily shaded areas are the most vulnerable to aquifer contamination. The light dotted areas are the safest. The former 1-, 10- and 25-year protection zones, computed with steady state groundwater flow models, were more or less concentric, and centred around the pumping station with increasing radius. In these earlier computations, horizontal solute transport in the aquifer towards the pump constituted the main part of the solute travel path. New provincial policy is aimed at minimizing pollution of the aquifer from which water is extracted, with lesser emphasis on the lateral distance to the pump. As a result, vertical solute transport in the unsaturated zone becomes much more important. Therefore, the new protection zones show a completely different picture. As explained in the introduction, .several processes determine the residence time of solutes in the unsaturated zone at a specific site, and therefore the contamination risk of the shallow aquifer at that site.
369
K. Oostindie and J. J. B. Bronswijk
104
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
675
.
.
Clay O~296
G 0
Average(1973)
73
126 ~
_
780
.
.
.
.
.
.
87
.
~
428
0
Wet (1974)
iiiiii+iiiii!iii
Dry (1976)
Figure 6. Computed water balance of the heavy clay soils in the Druten area for an average, a wet and a dry year. All quantities in mm per year. a, actual evapotranspiration; b, crack infiltration; c, matrix infiltration; d, runoff; e, flux to groundwater; f, flux to aquifer; g, flux to drains/ditches; h, change in storage.
3.2. MODEL VALIDATION
The quality of the model computations depends on a number of factors, such as the quality of the model principles and the quality of the input data. For residence time predictions, validation of model predictions using field data is very difficult, because long measurement periods would be required. However, it is possible to use other model output data, and compare them with available field data, in order to gain insight into the quality of the model computations. For the present study, we have used watertables to compare model predictions with field data. A water-table class map was available for the Druten area. On such maps, which are widely used in the Netherlands to characterize the hydrological conditions in a region, units are distinguished on the basis of the mean highest water-table (MHW)
370
Preferential flow and aquifer protection 50
500
ii
40
400
E
II
Crack infiltration
r-I
Matrix infiltration
E
E 300
~- 30 0
-g 217
--
200
O
c
"x
U
100
10
o
m
" 100
"
"
"
"
m
::
25O
o 800
50
500
40
400
30
300
20
200
10
100
0
0 100
25O
800
50
C
40
--
-
-
5OO
400
30
--
300
20
--
200
10
,,1! !
100
25O
100
0
800
Residence time (d)
Figure 7. Computed residence time distributions of water in the heavy clay soil area with an aquifer depth of 180era. a, average year (1973); b, wet year (1974); c, dry year (1976).
371
K. Oostindie and J. J. B. Bronswijk •'
':
:-
+:7.'7": :~........
"'~G= ~ - - : :
'~ '
'~
abr<-,;.~' •
_
• -', •
; -'
]::7]t-Oi'ilff{:~:
,~i"
',',~:,":-'
=~
'
.L~l
,.~-
'-"
"
I..
.
i
II -~
"-~.
~.
:
! 7":i..:%cI: :
".' I~'- ' '
~l.c-..'.'_-~-:~
-.
i
',~li'.~:q"7,
-
.
"~-
:.,
"
"-
:..~"
-?.,. '
;
- ,
........
le
'?-':=:
,
ill'""
~
,I';~
~'
/~'}i~'.~{l~ ~ I
:
"'
..... ~
'~
-6 7
";..
.
.
-,
r.
• . ".~.-
!-~.
. '
~7
..
t~
.
-
,:.',~'
."
.
,-'
'. .
..... i- . ,'<*
.
.
, ....
'
. ..'
-'
..
, '
*t:~,?.
Neersteind.'.
Safe (Class 5)to
,~:
.....
~..
.
~# . . . . " . . .
contamination
.
. ,";A~'~'~:.2Horssen"._i_~ . , . ; . , ' ~ { . ~ f . . L ~ i : ..0 ' ..;~-:,.il; q,,.;":;~-:",.:L ~•. . . .
(Class 3) ~
(Class 2) Very vulnerable
(Class 1)to
contamination
)
~
,
•
Laagveld ...................
,~_~ I~..~. oo.~.ze ~ :, ' " -
1-year protection
risk
-..
•
:..
' : ' : ~ . . ;, .=:. ,. L,:::,_ ,..., ~ ......... :c~. ~ _ L ~;L~,.,,L2/~.. ,=~-4::~-----.. ~................ '.-----. '- ,.; ....
b ~' ,:
~.
54.
Her
:
,
.
...-<~.,~.,
'
;
.
"
• .....
"~"
'
.'.#~)
~.
~ '
..: . . . . , . . . . . .
.
.
~>~
..'
:
. "''~
"" "'"
~
-'L~I
-';,
'
.~ -...., ~,~..<,,, . . . . . . . .
,
...... 2 . ~ 1 ~ t
, ..... .,,,.-,
~ ....... '
~'
~,',
":.~-'~
{'7~"
., I
,.., :
~
.~ ~..,'~,'.
.......
•
' ",
ri~!,~"
".'.;Ig'.lll
:
Horssep
.,<,
t#
'l~.£~!~"-';:"TiT;;"r";J','L,,.,"/`,
'"
"
t,.
i -i~,;.
'~"."~
Druten-:..--,-,F-_,-
'~
"'
-'
i-~l.l, , ,
zone
lO-year
protection
zone
25-year
protection
zone
risk
Figure 8. Computed assessment of contamination risk of the shallow aquifer in the Druten area.
Preferential flow and aquifer protection
372
TABLE3. Calculated and observed water-tables (cm below soil surface) in the Druten area Mapping unit
RN95A RN47C
Water-table class
VIo IIIb
Calculated
Observed
MHW
MLW
MHW
MLW
58 25
166 109
40-80 25--40
120-180 80-120
MHW, mean highest water-table; MLW, mean lowest water table.
and the mean lowest water-table (MLW). For instance, class 1 implies a mean lowest water-table of less than 50 cm below the soil surface (i.e. very wet conditions), whereas class 7 implies a mean highest water-table deeper than 80 cm and a mean lowest watertable deeper than 120 cm (i.e. very dry conditions). The model computations yield daily water-tables for all the map units. Computed MHWs and MLWs were derived from these computed water-table time series. A comparison between computed and observed MLW and MHW is given in Table 3. All the computed water-tables fall in the correct water-table class, which supports the validity of the applied model. 4. Conclusions
The developed method enables a quick and relatively easy assessment of the spatial distribution of residence times of water between soil surface and shallow aquifers in a clay soil area, with special attention to preferential flow through shrinkage cracks. In regions where drinking water is extracted from such shallow aquifers, the method can be applied to assess the contamination-risk of the aquifers, and to locate zones where residence times for water in the soil are short, and therefore the contamination-risk of the underlying aquifer is high. Measures to prevent aquifer contamination can be focused on those areas, for instance by restricting the use of pesticides and fertilizers. The present method assumes that any contamination of an aquifer from which water is extracted for drinking should be avoided. As a result, processes inside the unsaturated zone of the soil, including swelling and shrinkage and rapid transport of solutes by preferential flow through shrinkage cracks, play a major role. The contamination-risk maps, obtained using the present method, strongly deviate from traditional concentric protection zones computed with saturated groundwater flow models, which were aimed more at computing saturated horizontal transport in the aquifer to the pump. We would like to thank D. J. Groot-Obbink for preparing the soil map which is the basis of the present research project. The project was partially funded by the authorities of the Province of Gelderland, the Netherlands. References Beven, K. and Germann, P. (1982). Macropores and water flow in soils. Water Resources Research 18, 1311-1325. Booltink, H. W. G. and Bouma, J. (1993). Sensitivity analysis on processes affecting bypass flow. Hydrological Processes 3, 33-43. Bouma, J. (1984). Using soil morphology to develop measurements methods and simulation techniques for water movement in heavy clay soils. I L R I Publication 37, 298-315.
K. Oostindie and J. J. B. Bronswijk
373
Bronswijk, J. J. B. (1988). Modeling of water balance, cracking and subsidence of clay soils. Journal of Hydrology 97, 199-212. Bronswijk, J. J. B. (1991). Magnitude, modeling and significance of swelling and shrinkage processes in clay soils. PhD Thesis. Wageningen Agricultural University, Wageningen, the Netherlands. Hoogmoed, W. B. and Bouma, J. (1980). A simulation model for predicting infiltration into cracked clay soil. Soil Science Society of America Journal 44, 458-461. Oostindie, K. and Bronswijk, J. J. B. (1992). FLOCR--a simulation model for the calculation of water balance, cracking and surface subsidence of clay soils. Report 47, DLO-Winand Staring Centre, Wageningen, the Netherlands. Oostindie, K., Groot-Obbink, D. J. and Bronswijk, J. J. B. (1990). Watertransport door afdekkende kleilagen boven een aantal waterwinningen in Gelderland Rapport 95, DLO-Winand Staring Centre, Wageningen, the Netherlands. Van Lanen, H. A. J. (1984). Verblijftijd van water in de onverzadigde zone van zandgronden in gebieden met diepe grondwaterspiegels. 1-120 17, 9-16. W6sten, J. H. M., Bannink, M. H. and Beuving, J. (1987). Waterretentie- en doorlatendheidskarakteristieken van boven- en ondergronden in Nederland: de Staringreeks; herz. versie. Stiboka, Rapport 1932/ICW, Rapport 18, Wageningen (Netherlands).
Consequences of Preferential Flow in Cracking Clay Soils for Contamination-risk of Shallow Aquifers K. Oostindie and J. J. B. Bronswijk DLO Winand Staring Centre for Integrated Land, Soil and Water Research (SC-DLO), P.O. Box 125, 6700 AC Wageningen, Netherlands Received 3 March 1994
A simple method is presented to assess, on a regional scale, the contaminationrisk of shallow aquifers covered by cracking clay soils, with special emphasis on preferential flow of contaminants through shrinkage cracks. A water extraction area of several hundred hectares is divided into units with homogeneous soil types and hydrological conditions. For each unit, a onedimensional computer simulation model is used to compute residence times in the unsaturated zone of the cracking clay soils. Analytical equations are used to estimate residence times in the saturated zone. Finally, a classification system is proposed to summarize and visualize the model computations in a simple way. The developed method is illustrated with a case study from the Netherlands. It is known how contamination-risk maps of water extraction areas can be obtained. Policies to prevent contamination of the aquifer can be focused on the most vulnerable areas as indicated on such a map. The contamination-risk maps, obtained using the present method, deviate sharply from the traditional, more or less concentric protection zones, computed with saturated groundwater flow models.
Keywords: aquifer contamination, clay soil, shrinkage crack, preferential flow.
1. Introduction A saturated clay has a low hydraulic conductivity. As a result, the residence time of water flowing through saturated clay layers is long. Because of these long residence times and the generally high adsorption capacity of clay, contaminants will often be adsorbed, or transformed by chemical or microbiological processes, before passing a clay layer. Therefore, saturated clay soils overlying an aquifer effectively protect that aquifer against contamination originating from the soil surface. When saturated clay soils dry, however, their volume decreases with the result that there is subsidence of the soil surface and the appearance of shrinkage cracks. Water and solutes can be transported rapidly to the subsoil and the groundwater by preferential flow via these cracks (e.g. Beven and Germann, 1982; Bouma, 1984). In shallow aquifers, preferential flow in a cracked clay soil may lead to very short residence times for solutes from the soil surface to the aquifer. 359 0301-4797/95/040359+ 15 $08.00/0
© 1995 Academic Press Limited
360
Preferential flow and aquifer protection
At many locations in the Netherlands, pumping stations withdraw water from aquifers for drinking water. In the water extraction areas surrounding these pumping stations, special regulations are enforced to prevent contamination of the aquifer. Water extraction areas are divided into zones, depending on the estimated average travel times of water from the soil surface to the extraction point of water from the aquifer. Generally, one-year, 10-year and 25-year-zones are distinguished. In the one-year-zones, very strict regulations, for instance with respect to the use of pesticides by farmers, are enforced. For the 10- and 25-year-zones, less severe regulations apply. Several water extraction areas are located in the riverine area in the central part of the Netherlands. These water extraction areas have a similar hydrogeological profile: a Holocene heavy clay soil overlies an aquifer from which drinking water is extracted. In some cases, the clay deposits are rather thin (less than one metre). In such situations, preferential flow of solutes via shrinkage cracks to the subsoil and the aquifer may be a serious hazard for the water quality in the aquifer. To enable policy makers to develop scientifically based regulations for land use in the vicinity of drinking water extraction points in clay soil areas, it is essential to be aware of the risk of preferential flow on a regional scale, and its consequences for residence times of water and solutes in the clay soils overlying the aquifer. The transport of water through cracking clay soils is, however, a complex process. Qualitatively, it is determined by the following factors (Bronswijk, 1991): •
•
•
•
•
•
Soil properties The more severely the clay cracks, the more water will infiltrate into the cracks, and the greater the preferential flow will be. The higher the hydraulic conductivity of the soil matrix, the shorter the residence times will be. Precipitation surplus When the precipitation surplus is high, in general more water will finally reach the aquifer. Precipitation intensity In the case of showers with a high precipitation intensity, preferential flow will be high. Drainage In an intensively tile-drained soil, a considerable amount of the infiltrating water will flow via drains to surface waters, and less water will reach the aquifer. On the other hand, intensive drainage results in deeper water-tables and deeper cracks, leading to shorter residence times. Upward and downward seepage to or from the aquifer Downward seepage enhances downward flow and decreases residence time from soil surface to aquifer. For upward seepage, the opposite is true. The thickness of the clay layer Aquifers covered by thin clay layers are more vulnerable than those covered by thicker ones.
To compute residence times of water in the covering clay layers, it is necessary to consider all these factors together with their complex interactions, and the use of computer simulation models is indispensable. Because solute transport from soil surface to the aquifer is mainly convective, only the transport of water needs to be considered for a first estimation of residence times and contamination-risk. The objective of this paper is to present a simple method to assess the contaminationrisk of shallow aquifers covered by cracking clay soils on a regional scale, with special
K. Oostindie and J. J. B. Bronswijk
361
emphasis on preferential flow through shrinkage cracks. A water extraction area of a few hundred hectares is divided into units with homogeneous soil types and hydrological conditions. For each unit, a one-dimensional computer simulation model is used to compute residence times in the unsaturated zone of the cracking clay soils. Analytical equations are used to estimate residence times in the saturated zone. Finally, a classification system is proposed to summarize and visualize the model computations in a simple way. The developed method will be illustrated with a case study from the Netherlands. 2. Materials and methods
2.1.
BASIC PRINCIPLES OF THE METHOD
The developed method consists of six steps: 1. The boundaries of the area under consideration are defined. In the present study, we have used the 25-year-zones around the pumping stations. These are the ones where the residence time of water towards the extraction filter of the pumping station is less than 25-years, as computed with steady state groundwater flow models. 2. A soil survey is conducted on a 1:50000 scale to identify the soil types in the considered area. Subsequently, soil physical properties of the various soil types are collected. 3. The hydrological and climatic data, required for the model simulations, are collected. Areas with homogeneous hydrological conditions, such as drain depth, are grouped into hydrological units. For climatic data, an average year, a wet year, and a dry year, with respect to precipitation and evaporation, are selected. 4. Combination of the soil map and the hydrological map yields a simulation map, in which each unit is considered homogeneous with respect to: --soil profile and soil physical properties; ----drainage conditions; ----depth of the aquifer; --hydraulic head in the aquifer. 5. For each simulation unit, residence times are then computed for an average year, a wet year, and a dry year, using the collected input data in computer simulation models and analytical equations. 6. Finally, the computed residence times for the average, the dry and the wet year, are classified, and a contamination-risk map of the considered area is constructed. 2.2.
COMPUTING RESIDENCE TIMES IN THE UNSATURATED ZONE
Figure 1 shows the schematization as applied in the present research. Water transport in the unsaturated zone is calculated with the FLOCR computer simulation model (Bronswijk, 1988; Oostindie and Bronswijk, 1992). FLOCR computes unsaturated vertical water flow in cracking clay soils, taking into account swelling and shrinkage processes and preferential flow through cracks. The soil is divided into soil matrix and shrinkage cracks. Precipitation is divided between matrix and cracks. Water flow in the soil matrix is very slow. Water infiltrating into the cracks is assumed to reach the bottom of these cracks instantaneously, without lateral infiltration into crack walls. Due to limited contact area between water and crack walls, and weak diffusivity forces,
Preferential flow and aquifer protection
362 Unsaturated zone
Watertable
Saturated zone
Clay
ii ................................................................
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii i iiiiiiiiiiiiiiiiiiiiiiiii Figure 1. Schematization of the soil profile in riverine water extraction areas.
Variable crack volume
I E
F
Subsidence o
lI 2 3 4
27 28
29 3O ......
Matrix-crack system at time t Matrix-crack system at time t + At
Figure 2. Principles of the FLOCR model, for computing water balance, cracking and subsidence in clay soils./, infiltration rate in soil matrix (m s-I); let, part of total crack infiltration caused by rainfall intensity exceeding maximum infiltration rate of soil matrix (m s-~); Ic2, part of total crack infiltration caused by rainfall directly into the cracks (m s-t); E, actual evapotranspiration (m s-~); V, Darcy flux between two nodal points (m s-a); V,aL,,drain discharge (m s-I).
lateral infiltration has been f o u n d to be only a m i n o r fraction o f the soil water balance ( H o o g m o e d and B o u m a , 1980; Bronswijk, 1988; Booltink and B o u m a , 1993). The principles o f F L O C R are pictured schematically in Figure 2. F L O C R produces o u t p u t per one-fifth o f a day.
K. Oostindie and J. J. B. Bronswijk
363
Water in the cracks and water in the matrix eventually reach the groundwater. There, some of the water flows, via the subsoil, towards drains and ditches, and some will reach the aquifer. For water flowing through the cracks, the residence time in the unsaturated zone is set at zero days: to=0
(1)
where t~ is the residence time of water in the cracks (d). For water inside the soil matrix, the residence time in the unsaturated zone is approximated by (e.g. Van Lanen, 1984):
tu= ~Ps
(2)
where tu is the residence time of water in the matrix (d), Tu is the average thickness of the unsaturated zone (cm), ® is the volumetric water content, and Ps is the average precipitation surplus (cm d-t). Equation 2 yields only a crude estimation of the residence time in the unsaturated zone of the soil matrix. Due to upward and downward movements of the water in the root zone, and to subsequent precipitation and evapotranspiration, the residence time in the soil matrix, as computed using Equation 2, must be considered as average. This rather crude approximation is justifiable because of the much lower residence times of matrix flow in comparison with crack flow in clay soils. As a result, contamination risk is determined much more by crack flow than by matrix flow. 2.3.
COMPUTING RESIDENCE TIMES IN THE SATURATED ZONE
The Holocene layer in the majority of locations in the considered riverine area is a few metres thick, whereas water-tables are around l'50m or less. Therefore, in most situations, after reaching the groundwater table, a water particle still remains for some time in the clay layer, before reaching the aquifer. Because the clay is saturated here, the residence times will generally be high. The residence time in the saturated zone is calculated using the equation: t s - TsOs
(3)
v
in which ts is the residence time in the saturated zone (d), Ts is the thickness of the saturated zone (cm), v is the water flux through the saturated clay soil (cm d-l), and ®s is the saturated volumetric water content of the clay soil. The water flux (v) is the vertical flux between the water-table and the bottom of the clay layer which is calculated using Darcy's equation: v= -Ks AH
Az
(4)
where Ks is the saturated hydraulic conductivity of the 'clay layer (cm d-~), AH is the
Preferential flow and aquifer protection
364
50
500
40
400
I Crack infiltration I-1 Matrix infiltration
E E ,- 30
~E 300
0
g
.c_ 200 .x
-~- 20 U
10
100
0
0 100
250 Residencetime (d)
1000
Figure 3. Example o f a computed residence time distribution. Class width is 25 d.
difference between the water-table and hydraulic head of the aquifer (cm), and Az is the thickness of the saturated clay layer (cm). Finally, the total residence time of water in the clay layer, overlying the aquifer, is calculated using one of the following equations (see also Figure 1): 1.
For water in shrinkage cracks penetrating into the aquifer (Equation 1): t=0d
2.
For water in shrinkage cracks not penetrating into the aquifer (Equation 3): t =t,
3.
For water in the soil matrix (Equations (2) and (3)): t=tu+ts
2.4.
CLASSIFICATION OF RESIDENCE TIMES
By combining the FLOCR model and Equations (2), (3) and (4), a residence time is computed for each volume of water infiltrating at the soil surface. For a whole year, this computation yields a residence time distribution. An example is presented in Figure 3 for a hypothetical situation where the aquifer starts at a depth of 180 cm below the soil surface. For this situation, FLOCR, using the climatic conditions of 1973, has computed a total annual net infiltration (infiltration minus evapotranspiration) of 299 mm, subdivided in a matrix infiltration of 185 mm, and a crack infiltration of up to l l 4 r n m The 185mm of water infiltrating into the soil matrix has an average residence time of about 1000 d in the clay layer, before reaching the aquifer. Of the 114 rnm of water infiltrating into the cracks, 24 mm has a residence time of between 76 and 100 d in the clay layer, 39 mm has a residence time of between 101 and 125 d, 33 mm has a residence time of between 126 and 150 d, and so on.
K. Oostindie and J. J. B. Bronswilk
365
TABLE1. Classification of residence time distributions, based upon quantities of water reaching the aquifer and residence times of water in the covering clay layer. Quantity of water reaching the aquifer (mm a -I) >75 25-75 10-25 0-10
0-50 d
51-400 d
>400 d
1 1 2 3
3 4 4 5
5 5 5 5
The classes range from very vulnerable to contamination (class 1) to safe (class 5).
The contamination risk of the aquifer below this covering clay layer is now determined by the quantity of water that reaches the aquifer, and by the velocity with which the water reaches the aquifer. Therefore, a risk assessment based on the computed results presented in Figure 3 should consider both quantities of water reaching the aquifer and residence times inside the covering clay layer. To this end, the classification system of Table 1 has been developed. The five distinguished classes vary from very vulnerable (class 1) to safe (class 5). According to Table 1, the residence time distribution of Figure 3 should be assigned to class 3 because more than 75 mm water (in this case 114 mm) has a residence time of between 50 and 400 d in the clay layer. Using Table 1, the risk of aquifer contamination can be classified for each simulation unit, and made visible on a contamination-risk map. It should be emphasized that Table 1 is a subjective interpretation of the objective simulation results shown in Figure 3. An alternative classification system to the one developed and given in Table 1, would yield different categories. Although the objectivity of the model results is lost, a quick and easy interpretation of the data as presented in Figure 3 requires a classification system such as presented in Table 1. 2.5. LOCATION, SOILS AND MODEL INPUT PARAMETERS The developed method was applied to a number of water extraction areas in the Netherlands (Oostindie et al., 1990). In this article, the risk-assessment of the area of Druten is presented. The location of the Druten area is shown in Figure 4. The considered area was about 350 ha. Some 3 000 000 m 3 drinking water is extracted from the aquifer annually. In the Druten area, two soil units were distinguished: a heavy clay soil (Rn44C/ Rn47C on the 1:50 000 soil map of the Netherlands) and a clay loam soil (Rd90A/ Rn95A). The water retention curves and hydraulic conductivity curves for these soil types were taken from WOsten et al. (1987) and are shown in Figure 5. The depth of the aquifer was also determined in a field survey. The clay cover above the aquifer varied between zero and five metres. For upper boundary conditions, climatic data from the average year 1973, the wet year 1974, and the dry year 1976, were used. Precipitation and potential evapotranspiration are given in Table 2. 2.6. MODELVALIDATION Validating methods and models as applied in this study is extremely difficult. Such
366
Preferential flow and aquifer protection
Legend .,
1-year protection zone
0
2S-year protection zone
........
Provincial boundary
m.m
National border
0
10
20
I
I
I
30 km I
........."..
#.
~-°.T-° ,.,,'
.- .,-,., |
"'g cu,.~bo4
~
Verschuer~ ~
(~lmmerloo
~'~._.~, k;.~ ..,-'fx'-
Figure 4. Location of the investigated water extraction area of Druten.
K. Oostindie and J. J. B. Bronswilk
367
10 6
105 % "*...
lo 4 °
"0 t--
*o.o % °-°° % °.°o° % "..o.
lO 3 ~
a_
% °'°" ~N'°°°
lo 2
101
I 0.6
10 0 0.1
0.2
0.3
0.4
0.5
I 0.7
Volumetric water content
I""...
Topsoil Rn95A
1~ ~ k ~ A "7 "0
101 F ~ " , ~ . k
/
E
10 0 --
..............
TopsoilRn47C
------
SubsoilRn95A
"~'x "'.'~
t; "0 C
8
?.1
lO-,
"5 "1-
...... "-~--
10 -2
""-~. ..... o
t
I
10. 3 lO 0
lo ~
lo z
to ~
104
Pressure head (cm)
Figure 5. Water retention (a) and hydraulic conductivity (b) curves of the heavy clay soils in the Druten area.
validation would require long-term monitoring of water transport in the future on numerous locations, or historical data on water transport in the past. Neither options are realistic in practice. However, in previous studies, the applied FLOCR model has been validated with good results (e.g. Bronswijk, 1988). In the present study, we used existing data of minimum and maximum water-tables in the Druten area to validate our model results.
Preferential flow and aquifer protection
368
TABLE2. Climatic data used in the model computations of residence times for an average, a wet and a dry year Year
Precipitation (mm)
Potential evaporation (mm)
1993 (average) 1974 (wet) 1976 (dry)
779 993 536
559 581 646
3. Results and discussions
3.1.
WATER BALANCE AND RESIDENCE TIMES
Although two soil types have been distinguished, we will concentrate in the presentation of the results on the heavy clay soil area, because cracking and preferential flow are most pronounced in this soil type. The water balance of the heavy clay soils covering the shaUow aquifer, as computed with the FLOCR model, is pictured in Figure 6. In an average climatic year (1973), actual evapotranspiration is about 440mm and precipitation about 779mm. Of these 779mm, some 13% (104ram) flows via the shrinkage cracks as preferential flow to the subsoil, whereas the remaining 87% (675 ram) infiltrates into the soil matrix. Subtracting the actual evapotranspiration from these 675 mm gives a net matrix infiltration of about 235 mm. Because the grasslands in this area are poorly drained, the drain discharge is negligible. The net flux of water into the aquifer in an average climatic year is 296 mm per year. For a wet and a dry year the percentages of preferential flow of the total precipitation are 13% (126mm) and 14% (73 mm), respectively. The quantities of water flowing into the aquifer are 428 and 189 mm per year, for a wet and a dry year, respectively. For a simulation map unit with a heavy clay soil and an aquifer depth of 180 cm below the soil surface, the residence time distributions for the average, wet, and dry year, are shown in Figure 7. The residence times for water infiltrating into the soil matrix are about 800 d. For water flowing into the shrinkage cracks, the residence times are between 75 and 275 d. Such residence time distributions were derived for each simulation map unit. These residence time distributions have been classified using the classification system of Table 1. The present example would fall into class 3. A classification was obtained, using the procedure outlined above, for each of the simulation map units. The resulting contamination-risk map is shown in Figure 8. The heavily shaded areas are the most vulnerable to aquifer contamination. The light dotted areas are the safest. The former 1-, 10- and 25-year protection zones, computed with steady state groundwater flow models, were more or less concentric, and centred around the pumping station with increasing radius. In these earlier computations, horizontal solute transport in the aquifer towards the pump constituted the main part of the solute travel path. New provincial policy is aimed at minimizing pollution of the aquifer from which water is extracted, with lesser emphasis on the lateral distance to the pump. As a result, vertical solute transport in the unsaturated zone becomes much more important. Therefore, the new protection zones show a completely different picture. As explained in the introduction, .several processes determine the residence time of solutes in the unsaturated zone at a specific site, and therefore the contamination risk of the shallow aquifer at that site.
369
K. Oostindie and J. J. B. Bronswijk
104
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
675
.
.
Clay O~296
G 0
Average(1973)
73
126 ~
_
780
.
.
.
.
.
.
87
.
~
428
0
Wet (1974)
iiiiii+iiiii!iii
Dry (1976)
Figure 6. Computed water balance of the heavy clay soils in the Druten area for an average, a wet and a dry year. All quantities in mm per year. a, actual evapotranspiration; b, crack infiltration; c, matrix infiltration; d, runoff; e, flux to groundwater; f, flux to aquifer; g, flux to drains/ditches; h, change in storage.
3.2. MODEL VALIDATION
The quality of the model computations depends on a number of factors, such as the quality of the model principles and the quality of the input data. For residence time predictions, validation of model predictions using field data is very difficult, because long measurement periods would be required. However, it is possible to use other model output data, and compare them with available field data, in order to gain insight into the quality of the model computations. For the present study, we have used watertables to compare model predictions with field data. A water-table class map was available for the Druten area. On such maps, which are widely used in the Netherlands to characterize the hydrological conditions in a region, units are distinguished on the basis of the mean highest water-table (MHW)
370
Preferential flow and aquifer protection 50
500
ii
40
400
E
II
Crack infiltration
r-I
Matrix infiltration
E
E 300
~- 30 0
-g 217
--
200
O
c
"x
U
100
10
o
m
" 100
"
"
"
"
m
::
25O
o 800
50
500
40
400
30
300
20
200
10
100
0
0 100
25O
800
50
C
40
--
-
-
5OO
400
30
--
300
20
--
200
10
,,1! !
100
25O
100
0
800
Residence time (d)
Figure 7. Computed residence time distributions of water in the heavy clay soil area with an aquifer depth of 180era. a, average year (1973); b, wet year (1974); c, dry year (1976).
371
K. Oostindie and J. J. B. Bronswijk •'
':
:-
+:7.'7": :~........
"'~G= ~ - - : :
'~ '
'~
abr<-,;.~' •
_
• -', •
; -'
]::7]t-Oi'ilff{:~:
,~i"
',',~:,":-'
=~
'
.L~l
,.~-
'-"
"
I..
.
i
II -~
"-~.
~.
:
! 7":i..:%cI: :
".' I~'- ' '
~l.c-..'.'_-~-:~
-.
i
',~li'.~:q"7,
-
.
"~-
:.,
"
"-
:..~"
-?.,. '
;
- ,
........
le
'?-':=:
,
ill'""
~
,I';~
~'
/~'}i~'.~{l~ ~ I
:
"'
..... ~
'~
-6 7
";..
.
.
-,
r.
• . ".~.-
!-~.
. '
~7
..
t~
.
-
,:.',~'
."
.
,-'
'. .
..... i- . ,'<*
.
.
, ....
'
. ..'
-'
..
, '
*t:~,?.
Neersteind.'.
Safe (Class 5)to
,~:
.....
~..
.
~# . . . . " . . .
contamination
.
. ,";A~'~'~:.2Horssen"._i_~ . , . ; . , ' ~ { . ~ f . . L ~ i : ..0 ' ..;~-:,.il; q,,.;":;~-:",.:L ~•. . . .
(Class 3) ~
(Class 2) Very vulnerable
(Class 1)to
contamination
)
~
,
•
Laagveld ...................
,~_~ I~..~. oo.~.ze ~ :, ' " -
1-year protection
risk
-..
•
:..
' : ' : ~ . . ;, .=:. ,. L,:::,_ ,..., ~ ......... :c~. ~ _ L ~;L~,.,,L2/~.. ,=~-4::~-----.. ~................ '.-----. '- ,.; ....
b ~' ,:
~.
54.
Her
:
,
.
...-<~.,~.,
'
;
.
"
• .....
"~"
'
.'.#~)
~.
~ '
..: . . . . , . . . . . .
.
.
~>~
..'
:
. "''~
"" "'"
~
-'L~I
-';,
'
.~ -...., ~,~..<,,, . . . . . . . .
,
...... 2 . ~ 1 ~ t
, ..... .,,,.-,
~ ....... '
~'
~,',
":.~-'~
{'7~"
., I
,.., :
~
.~ ~..,'~,'.
.......
•
' ",
ri~!,~"
".'.;Ig'.lll
:
Horssep
.,<,
t#
'l~.£~!~"-';:"TiT;;"r";J','L,,.,"/`,
'"
"
t,.
i -i~,;.
'~"."~
Druten-:..--,-,F-_,-
'~
"'
-'
i-~l.l, , ,
zone
lO-year
protection
zone
25-year
protection
zone
risk
Figure 8. Computed assessment of contamination risk of the shallow aquifer in the Druten area.
Preferential flow and aquifer protection
372
TABLE3. Calculated and observed water-tables (cm below soil surface) in the Druten area Mapping unit
RN95A RN47C
Water-table class
VIo IIIb
Calculated
Observed
MHW
MLW
MHW
MLW
58 25
166 109
40-80 25--40
120-180 80-120
MHW, mean highest water-table; MLW, mean lowest water table.
and the mean lowest water-table (MLW). For instance, class 1 implies a mean lowest water-table of less than 50 cm below the soil surface (i.e. very wet conditions), whereas class 7 implies a mean highest water-table deeper than 80 cm and a mean lowest watertable deeper than 120 cm (i.e. very dry conditions). The model computations yield daily water-tables for all the map units. Computed MHWs and MLWs were derived from these computed water-table time series. A comparison between computed and observed MLW and MHW is given in Table 3. All the computed water-tables fall in the correct water-table class, which supports the validity of the applied model. 4. Conclusions
The developed method enables a quick and relatively easy assessment of the spatial distribution of residence times of water between soil surface and shallow aquifers in a clay soil area, with special attention to preferential flow through shrinkage cracks. In regions where drinking water is extracted from such shallow aquifers, the method can be applied to assess the contamination-risk of the aquifers, and to locate zones where residence times for water in the soil are short, and therefore the contamination-risk of the underlying aquifer is high. Measures to prevent aquifer contamination can be focused on those areas, for instance by restricting the use of pesticides and fertilizers. The present method assumes that any contamination of an aquifer from which water is extracted for drinking should be avoided. As a result, processes inside the unsaturated zone of the soil, including swelling and shrinkage and rapid transport of solutes by preferential flow through shrinkage cracks, play a major role. The contamination-risk maps, obtained using the present method, strongly deviate from traditional concentric protection zones computed with saturated groundwater flow models, which were aimed more at computing saturated horizontal transport in the aquifer to the pump. We would like to thank D. J. Groot-Obbink for preparing the soil map which is the basis of the present research project. The project was partially funded by the authorities of the Province of Gelderland, the Netherlands. References Beven, K. and Germann, P. (1982). Macropores and water flow in soils. Water Resources Research 18, 1311-1325. Booltink, H. W. G. and Bouma, J. (1993). Sensitivity analysis on processes affecting bypass flow. Hydrological Processes 3, 33-43. Bouma, J. (1984). Using soil morphology to develop measurements methods and simulation techniques for water movement in heavy clay soils. I L R I Publication 37, 298-315.
K. Oostindie and J. J. B. Bronswijk
373
Bronswijk, J. J. B. (1988). Modeling of water balance, cracking and subsidence of clay soils. Journal of Hydrology 97, 199-212. Bronswijk, J. J. B. (1991). Magnitude, modeling and significance of swelling and shrinkage processes in clay soils. PhD Thesis. Wageningen Agricultural University, Wageningen, the Netherlands. Hoogmoed, W. B. and Bouma, J. (1980). A simulation model for predicting infiltration into cracked clay soil. Soil Science Society of America Journal 44, 458-461. Oostindie, K. and Bronswijk, J. J. B. (1992). FLOCR--a simulation model for the calculation of water balance, cracking and surface subsidence of clay soils. Report 47, DLO-Winand Staring Centre, Wageningen, the Netherlands. Oostindie, K., Groot-Obbink, D. J. and Bronswijk, J. J. B. (1990). Watertransport door afdekkende kleilagen boven een aantal waterwinningen in Gelderland Rapport 95, DLO-Winand Staring Centre, Wageningen, the Netherlands. Van Lanen, H. A. J. (1984). Verblijftijd van water in de onverzadigde zone van zandgronden in gebieden met diepe grondwaterspiegels. 1-120 17, 9-16. W6sten, J. H. M., Bannink, M. H. and Beuving, J. (1987). Waterretentie- en doorlatendheidskarakteristieken van boven- en ondergronden in Nederland: de Staringreeks; herz. versie. Stiboka, Rapport 1932/ICW, Rapport 18, Wageningen (Netherlands).