The effects of Pleistocene glaciations on the geohydrological system of Northwest Europe

Journal

of Hydrology Journal of Hydrology 195 (1997) 137-159

ELSEVIER

The effects of Pleistocene glaciations on the geohydrological system of Northwest Europe F . H . A . v a n W e e r t a'*, K. v a n Gijssel b, A. L e i j n s e a'c, G.S. B o u l t o n d =National Institute of Public Health and Environmental Protection, P.O. Box 1, 3720 BA Bilthoven. Netherlands bGeological Survey of The Netherlands, Haarlem, Netherlands CDelmrtment of WaterResources, Agricultural University Wageningen, Wageningen, Netherlands dDepartment of Geology and Geophysics, University of Edinburgh, Edinburgh, UK

Received 31 May 1995, revised 13 March 1996; accepted 2 August 1996

Abstract A large-scale hydrological model study is carried out to investigate the geohydrological responses to glacial climate conditions in Northwest Europe. The vertically integrated groundwater model is based on a supraregional hydrogeological model of the Cenozoic and Mesozoic subsurface in Northwest Europe. Three different layers are distinguished above the relatively impervious base of Palaeozoic and Precambrian rocks. Boundary conditions are inferred from indicative palaeoenvironmentai reconstructions of the last three glacial cycles to simulate groundwater flow related to ice sheet expansions into theNorthwest European lowlands. Recharge of the groundwater system due to basal glacial melting is deduced from ice sheet model simulations driven by a transient climate function. Results of the large-scale model study show relatively high groundwater velocities and pressures in the subgiacial areas and the ice-marginal permafrost areas. Extreme high velocities may develop when the ice sheet has advanced to the southern margins of the upper Plio/Pleistocene aquifer in Northwest Europe. Drainage of the highly pressurized groundwater system mainly occurs in proglacial ice-dammed lakes, ice-marginal seas and zones of discontinuous permafrost, including river valleys.

1. I n t r o d u c t i o n In Europe several geological formations are considered as possible hosts for the disposal of high-level radioactive waste: rock salt, clay, shale and granite. To be able to carry out a safety assessment study of such a disposal site, one must understand the processes that * Conesponding author at: Department of Water Management, Environmental mui SanitaryEngineering, Delft University of Technology, PO Box 5048, 2600 GA, Delft, Netherlands. 0022-1694/97/$17.00 © 1997 Published by Elsevier Science B.V. All rights reserved PH S0022-1694(96)03248-9

138

F.H.A. van Weert et alJJournal o f Hydrology 195 (1997) 137-159

play a role and how these processes may change in future. Groundwater flow plays an important role in all scenarios where radionuclides are released from a repository in a geological formation. It may undermine the integrity of the disposal site and forms a means of transport of radionuclides to the biosphere (International Atomic Energy Agency, IAEA, 1981). Therefore, it is extremely important to consider the geohydrological system in long-term risk assessment studies. It is to he expected that long-term climatic changes like glaciations will greatly affect the NW European groundwater system. Recent developments in mathematical modelling, and the improved knowledge of the European geology and climatic changes, allow us to model some geohydrological effects of past and possible future glaciations. The present model studyis conducted with the objective to obtain general knowledge of the behaviour of the Northwest European geohydrological system under Pleistocene glacial conditions. It is based on knowledge of Pleistocene palaeociimatic changes. The model simulates groundwater pressures and velocities and focuses on the geographical indication of areas where the integrity of possible disposal sites might be undermined by geohydrological processes. In the next section a short overview is given on possible effects of glaciations on the geohydrological system. Next, we discuss the hydrogeological model of NW Europe which forms the geometrical framework of the groundwater model. Further, the concepts, assumptions and mathematics of the two-dimensional vertically integrated groundwater model are discussed. Palaeogeographical reconstructions referring to various stages in the Pleistocene glacial periods are described. The boundary conditions which are deduced from these reconstructions are discussed. Next, the simulation results regarding hydraulic heads, groundwater velocities and recharge/discharge patterns are shown and discussed for the given reconstructions. Finally, conclusions based on this model study are drawn.

2. Glaciations and geohydrologicai response It is generally accepted that Pleistocene glaciations were initiated by global cooling of the atmosphere. Large amounts of ocean water were stored in large sheets of continental ice, and consequently the sea level dropped globally, usually more than 100 m. Most glaciers and continental ice sheets melt basally because of thermal friction due to basal sliding and the geothermal heat flux. If the ice sheet rests on a lithified, relatively impermeable bed, drainage of basal meltwater may occur in episodic basal sheet flows and/or through a subgiacial system of tunnels in the ice-marginal zone. Relicts of these latter drainage pathways are the esker systems found in many parts of Fennoscandia. The basal water sheets may act as sliding planes causing large surges of ice lobes (Shoemaker, 1992). Transport of the subgiacial meltwater by groundwater flow through the underlying bed may he negligible because of the very low permeability of the underlying rock. For an ice sheet overlying unconsolidated sediments it is assumed that when the groundwater pressure is smaller then the overburden pressure, all the subgiacial meltwater infiltrates into the sediments (Boulton et al., 1987). In case of high meltwater production and/ or small aquifer thickness and/or small permeability values, the groundwater pressure may become very high, equal to the overburden pressure. As a result, the effective

F.H.A. van Weert et al./Journal of Hydrology 195 (1997) 137-159

139

pressure in the porous medium may be very small and the shear strength of the sediments is reduced. Such conditions favour deformation, thrusting and erosion of the sediments (Bonlton and Jones, 1979; Shoemaker, 1986; Mooers, 1990). Deep subglacial channel valleys may have been formed in the weakened sediments to discharge the excess meltwater (Boulton et al., 1993). Relicts of these channel valleys are the sedimentary channel fills found for instance in the northern parts of Germany. During the glaciations in the Elsterian stage and the subsequent Saalian stage ice sheets from Scandinavia and the British Isles overrode the extensive sedimentary aquifer systems of Northwest Europe. The area overridden was subjected to total landscape remodelling, with old river systems destroyed or buried (Gibbard, 1988). At the ice margin, major river valleys were dammed all across the region during the Pleistocene glaciations. The German rivers (Elbe, Weser, Rhine) were deflected westwards and changed their courses continually. Large temporary lakes were formed along the ice margin. Gibbard (1988) concluded that a massive ice-dammed lake in the southern North Sea had developed during the Elsterian glaciation. The main Northwest European rivers discharged into this lake. Overspill of this lake almost certainly initiated the Straits of Dover. In the area in front of the ice sheets the meltwater released from the ice together with large amounts of sediments had to be carried away by proglacial braided rivers (sandr). Groundwater movement in permafrost zones was restricted by the presence of both perennially and seasonally frozen ground. The permafrost, being nearly impermeable, acted as a confining layer (French, 1976). This reduced the aquifer thickness and hence the transmissivity. An important factor influencing the groundwater discharge was the absence or discontinuity of permafrost along the rivers and underneath seas and lakes. Thus, these rivers, seas and lakes did not only act as outlets to discharge surface runoff, but they probably also played an important role in draining the highly pressurized groundwater system. Obviously, thick ice sheets caused immense overburden pressures in the underlying sediments. Under high hydraulic diffusivity conditions this may have resulted in high effective pressures and consequently possible compaction. The compaction should have had consequences both for the porosity and permeability. Ice sheet loading is also assumed to have affected halokinesis, i.e. the autonomous isostatic movement of salt formations (Prij and Benneker, 1989). The low density of rock salt compared with the density of the sedimentary overburden causes gravitational instability. This instability is influenced by strong horizontal lithostatic pressure gradients caused by differential glacial loading. Another result of ice sheet loading is the isostatic reaction of the lithosphere. Because of the extra ice load, the earth's crust is thought to depress and flex. The Pleistocene depressions lead to significant topographic changes and forced river courses to migrate.

3. The hydrogeological model The geology in the basin area is strongly affected by salt tectonics and block faulting (Ziegler, 1991). Thickness variations are mainly related to these tectonic processes. During several stages halokinesis of Rotliegend and Zechstein evaporates occurred in the

140

F.H.A. van Weert et alYJournal of Hydrology 195 (1997) 137-159

northern parts of The Netherlands, Lower Saxony and the southern part of SchleswigHolstein. Locally, salt diapirs even displaced Tertiary and Quaternary deposits and reached up regionally to less than 100 m below the present ground surface. Consequently, many fault systems developed in the overlying strata. Mostly, halokinesis was tectonically influenced and related to phases of regional uplift and gx~bsidence. During Late Jurassic and Early Cretaceous subsidence of the Lower Saxony Basin and the Central Netherlands Basin took place. In the Upper Cretaceous and Palaeocene these basins were uplifted again and became inverted. In Western and Central Europe the accumulation of chalky carbonates ended at the Danian/Thanetian transition. It was followed by a clastic dominated depositional regime. In the Cenozoic a large new sedimentary basin developed, the North Sea Basin, of which the southern part lies in the model area. Sea level fluctuations strongly influenced the sedimentation pattern in the North Sea Basin, particularly in its shallower marginal parts during the Tertiary. Furthermore, the upper (Cenozoic) part of the model sequence is influenced by Quaternary glacial processes. The ice-sheet influenced groundwater flow model is based on a supraregional hydrogeological model extending from Great Britain to the German-Polish border and from southern Scandinavia to the German highlands (Fig. 1). On this large scale it is assumed that the N-W European hydrogeological system can be modelled adequately by subdividing

Hydrogoological rasp of m o ~ l area Ihnit1

m

l hit3

Fig. 1. Hydrogeologicalmap of modelarea. A subdivisioninto threeunits: Unit I consistsof the relativelyhighly permeable,unconsolidatedQuaternaryand (Upper)Tertiarydesposits,Unit 2 consistsof the semi-perviousclays of the (Lower)Tertiaryand Units 3 consist of the consolidated, low permeable and semi-perviousMesosoic depsosits.The basementconsistsof relativelyimperviousdepositsof Palaeozoicand Precamhrianage, including the Permianrock salt. Furthermore,the locationof profileline I-II is shown.

F.H.A. van Weert et al./Journal of Hydrology 195 (1997) 137-159

141

it into three units, two aquifers separated by a thick aquitard. This means that several regionally significant layers had to be excluded. The units are chosen according to lithological and hydrogeological properties (porosity, permeability) and the spatial distribution and extension of the strata. The division is based on the lithostratigraphic correlation charts of the Mesozoic and Cenozoic sedimentary records in the southern North Sea and northern Central Europe (Ziegler, 1991). However, the hydrogeological boundaries do not always coincide with the lithostratigraphical boundaries, because of lateral variations in sedimentary facies of the stratigraphical units. When the hydrogeological boundaries cannot be easily detected lithostratigraphical boundaries are used as the boundaries between the different units. Concerning the distribution and thickness of the units the following criteria are used: •

•

Only structures with a minimal lateral extent of 15 x 15 km 2 are considered, with the exception of some salt diapirs; consequently, relatively small structures like the salt-related fault systems and the Elsterian channels are neglected. The units have a minimal thickness of 50 m.

Based on this approach, the following three hydrogeological units have been distinguished: 1. the unconsolidated Quaternary and (Upper) Tertiary, sandy deposits--locally (Lower) Tertiary and Cretaceous sandy deposits are included; 2. the semi-pervious clay beds of (Lower) Tertiary age; 3. the consolidated Mesozoic deposits consisting of Upper Cretaceous (and Danian) chalk (permeable when fissured) and low permeable and semi-pervious Lower Cretaceous, Jurassic and Triassic deposits. The base of the model consists of the relatively impervious deposits of Palaeozoic and Precambrian age, including the rock salt of the Perinian Zechstein Group. This division is consistent with the International Hydrogeological Map of Europe (Karrenberg et al., 1973) and reports on the groundwater resources of the European Community (Commission of the European Communities, CEC, 1982). The values of the intrinsic permeability and the effective porosities of the different units are derived from data reported in literature (Karrenberg et al., 1973; Johannsen, 1980). The lithological characteristics of the various units are given in Table 1. Contour maps of the top of each hydrogeological unit were compiled from various geological and structural review maps. The distribution of the thicknesses of the various hydrogeological units derived from these maps are shown in Fig. 2.

4. The two-dimensional vertically integrated groundwater model Because of the large areal extent of the basin compared with its thickness and based on the objectives of this study the system is modelled as a two-dimensional vertically integrated layered system. In this type of model the geohydrological units are represented by layers with spatially varying transmissivities. In the layers representing aquifers the groundwater flow is assumed to be horizontal only, while groundwater flow in layers representing aquitards is strictly vertical. Boundary conditions are given at the side

Quarternary Holocene Pleistocene Tertiary Pliocene Miocene Miocene Oiigocene Eocene Palaeocene Mesozoic Upper Cretaceous Lower Cretaceous Jurassic Triassic Paleozoic Permian (+pre-Permian)

1

Base

3

2

Stratigraphy

Unit

Evaporites

Chalk, marl, claystone, shale, sandstone, evaporites

Clay, some sandy intervals

Fine to coarse sand, locally clay

Lithology

Table I Hydrogeological characteristics of the units in the supraregional model

1-30

0.30

Impermeable

1 x 10 -5

< 1 (low)

0.30

0.30

20

200-50000 (high)

0.35 •

0.1 (overlaid)

I (outcrop)

Hydr. cond. (m day -I)

Permeability (roD)

Effective porosity (-)

::t

.:e

to

F.H.A. van Weert et al3.1ournal of Hydrology 195 (1997) 137-159

T

~

(me~) o~~ o ~ c ~

.nits UnJtl

]

O
m u m°,.

~t2

I

0,<500

m l~.C

m..,_

Unit3

m~_ m.,_ m..~

0

l~O

Im

Fig. 2. Isopaches of various hydrogcological units.

~m0ow.=

143

144

F.H.A,

van Weert et alJJournal of Hydrology 195 (1997) 137-159

boundaries of the layers. Recharge and discharge at the top of the layered system are described by a groundwater head dependent flux, representing different processes. The governing equation for groundwater flow in porous media follows from the principle of conservation of mass. It is assumed that on this scale the equation of conservation of momentum is given by Darcy's Law: q=nv= -K.Vh

(1)

where q is the specific discharge (L T-I), n is the porosity (-), v is the actual velocity (L T-I), K =pgkll~ the hydraulic conductivity tensor (L T-I), p the liquid density (M L-J), g the gravity constant (L T-2), k the intrinsic permeability tensor (L2), /~ the liquid viscosity (M L -1 T-l), V the gradient operator (L-l), h = (p/pg + z) is the piezometric head (L), p the pressure in the liquid phase (M L -l T -2) and z the vertical coordinate (L). It is assumed that the density of groundwater is not influenced by dissolved salt concentrations and depends on the pressure p only. Substitution of Darcy's Law and the equation of state for the liquid density in the mass conservation equation gives:

So ~t = V.(KVh)+ R

(2)

where So = pg(~S + n~) is the storage coefficient (L-I),/~s and/~ are the compressibilities of the solid grains and water, respectively (L T 2 M -I) and R is a sink and source term (T -1) that accounts for mass exchange by freezing and thawing. As has been stated before, possible compaction of the sediments due to ice-sheet loading may affect the hydraulic properties of the sediments. In this study, however, we will neglect these effects and will assume that both So and K are independent of the groundwater head h. Consequently, results of the model study should be interpreted as indications of what will happen under glacial conditions. Eq. (2) can be vertically integrated over the thickness of each aquifer. Under the assumption that the groundwater head h is independent of z (Dupuit-assumption), and that the principal axes of the conductivity tensor coincide with the coordinate axes, the following two-dimensional flow equation is obtained:

sOh= L : T dh~ d__:T' Oh~ Ot Ox ~, x Ox,] + Oy ~, y Oy] + qt - qb + Ra

(3)

where S is the storativity (-), T~ and Tr are the components of the transmissivity tensor in x- and y-directions (L 2 T-l), respectively, qt and qb are the vertical fluxes at top and bottom of the aquifer, respectively, (L T -l) (positive downward) and R a is the sink/source term R integrated over the thickness of the aquifer (L T-I). Since thawing and freezing only occurs in the top aquifer, this term is zero in the bottom aquifer. The storativity S is the storage coefficient So integrated over the thickness of the aquifer, and the transmissivities are the integrated values of the hydraulic conductivities. For the top aquifer (Unit 1) qt accounts for the groundwater recharge of the aquifer through basal meltwater influx, meteoric recharge and interactions of groundwater with surface water bodies (seepage) while qb is the flux through the aquitard (Unit 2). For the bottom aquifer (Unit 3), qt is the flux through the aquitard, while qb is identically zero because the base of the system is assumed to be impervious. Only steady-state situations will be considered in this study. If the sink/source term Ra

F.H.A. van Weert et alJJournal of Hydrology 195 (1997) 137-159

145

for the top aquifer is included in the top flux qt, Eq. (3) reduces to: ax \

ax]

= qb-- q~

(4)

The vertical movement of groundwater through aquitard is dependent on the hydraulic heads above and below the aquitard:

ht-h3 hi-h3 q~l=Kz'clay bclay- Cclay

(5)

where qaq is the leakage through the aquitard, positive downward (L T-l), K ~ y is the vertical component of the hydraulic conductivity of the aquitard (L T-l), hi and h3 are the hydraulic heads in the top and bottom aquifer, respectively, (L), bc~yis the thickness of the aquitard (L) and Cclayis the hydraulic resistance of the aquitard (T). Eq. (5) defines qb for Unit 1 and qt for Unit 3. The recharge and discharge of the top aquifer (q0 is described as a continuous function composed of linear segments in terms of the groundwater head in the top aquifer. These functions may vary spatially and allow the mathematical description of different types of recharge/discharge processes. An approximate solution of Eq. (4), together with Eq. (5) and the prevailing relation for the recharge/discharge of the top aquifer is obtained by a finite element model. Once the hydraulic head distribution is known, the groundwater velocities can be calculated from Darcy's Law, Eq. (1). Based on the objectives of the study and the available data, and restricted by the conslTaints of numerical modelling, simplifying assumptions have been made for the present model. These assumptions pertain to the hydraulic characteristics of the units and the recharge/discharge of the top aquifer. The hydraulic conductivities of the various hydrogeological units are assumed to be homogeneous and iso¢opic. For the lower aquifer, two homogeneous zones are distinguished: the outcrop with a relatively high hydraulic conductivity, and the part overlaid by the aquitard with relatively low hydraulic conductivity. The transmissivides of the aquifers are approximated by the product of hydraulic conductivity and aquifer thickness. Similarly, the hydraulic resistance of the aquitard is proportional to the thickness of the aquitard. In the continuous and discontinuous permafrost zones, part of the top aquifer is permanently frozen. In these frozen parts horizontal groundwater flow is negligible. Therefore, the aquifer thickness available for groundwater flow, and thus the transmissivity, is reduced in these zones. The lateral extent of the various units differ. For instance, the lateral extent of the upper aquifer is smaller than the extent of the underlying aquitard and lower aquifer. To define a three-layer system everywhere in the model domain, small transmissivity values are assigned to the upper aquifer in those zones where the hydrogeological unit is missing. In this way zones are created where hardly any borizontal flow in the aquifer occur. In areas where the aquitard is missing, a small value of the hydraulic resistance is assigned, thus creating a direct communication between the top and bottom aquifer or between the recharge and the bottom aquifer. The geohydrological response of the system is not significantly altered by this approach.

146

F.H.A. van green et alTJournal of Hydrology 195 (1997) 137-159

To describe the different groundwater recharge/discharge processes for the top aquifer four different zones are distinguished: (1) zones covered by ice sheets, (2) permafrost zones, (3) zones covered by surface waters (seas, lakes, large rivers), (4) no-permafrost zones. 1. Groundwater recharge under the ice sheet is by subglaciai meltwater infiltration. In this model the process of meltwater infiltration is based on the effective pressure principle. The infiltration is constant and equal to the subglacial meltwater production rate if the groundwater pressure is less than the overburden pressure (effective pressure >0). It is assumed that if this condition is not satisfied, the ice sheet will be lifted and/or deep channels will be formed in the sediments, to discharge the excess meltwater and to drain the overpressurized groundwater. However, in the steady-state groundwater model presented here these failure processes are not incorporated and the groundwater pressure is allowed to exceed the overburden.pressure. In this condition the groundwater is assumed to discharge subglacially with a maximal rate equal to the meltwater production rate (Fig. 3). The model as presented here is formulated in terms of

i Orounclwat~ head ( m r ) in upper aquifer

_ +1 meter

i

Fig. 3. Subglacialrecharge/dischargefunction in terms of the groundwaterhead in the upper aquifer. When the groundwaterhead is less then the hydraulic potentialequivalenceof the overlyingice sheetgroundwaterrecharge will take place by means of meltwater infiltration. When the groundwater pressure exceeds this hydraulic potential equivalence groundwater will discharge subglacially. The maximal infiltafion rate is equal to the meltwater production rate.

F.H.A. van Weert et alJJournal of Hydrology 195 (1997) 137-159

147

hydraulic heads rather than pressures. The condition of a positive effective pressure is satisfied when the groundwater head is less than the hydraulic potential equivalent of the overlying ice sheet (i.e. approximately 0.9 times the ice sheet thickness). The transition of recharge to discharge is not discrete but defined over a hydraulic head difference of 2 m. Such a transition is required to prevent numerical problems in the solution of the governing equation. In this study, supraglacial meltwater does not contribute to the basal meltwater infiltration, i.e. only subglacially produced meltwater is considered. Typical values of basal meltwater production are 0-15 mm year-l. 2. In the continuous permafrost zone, no recharge or discharge of groundwater takes place because the aquifer is totally frozen at the top. In the discontinuous permafrost zone, only discharge will occur during the advancing stage of a glacier. Recharge is assumed to be negligible due to constant freezing of surface water and excess precipitation. In the deglaciation stage, however, excess precipitation, melting surface ice and melting permafrost will contribute to the groundwater recharge in this zone. In the discontinuous permafrost zones the recharge and discharge processes are linear functions of the groundwater head in the upper aquifer. The resistance to recharge/discharge is dependent on the permeability of the aquifer, the thickness of the permafrost zone and the area available for recharge. 3. In zones covered with large bodies of running or standing surface water (rivers, seas and lakes) the underground will not be frozen permanently, and infiltration/discharge of the top aquifer is possible. The rate of infiltration/drainage is proportional to the groundwater head in the aquifer. Infiltration (surface water entering the groundwater system) will occur if the groundwater head is lower than the surface water level. Drainage (groundwater seeping into the surface water bodies) occurs if the groundwater head exceeds the surface water level. In general, the hydraulic resistances between the surface water bodies and the top aquifer are small, such that the groundwater head practically equals the surface water level. 4. In the non-permafrost zones, meteoric recharge takes place. Annual precipitation rates in Northern Europe were low in the glacial periods: in The Netherlands approximately 200-500 mm in contrast to the mean annual precipitation of 600-900 mm during interglacial periods (Zagwijn et al., 1992). In this model, the meteoric recharge is independent of the groundwater head. The amount of precipitation actually recharging the groundwater system during the cold periods varies spatially from 0 to 50 mm year-l. During the warmer, interglacial periods meteoric replenishment varies with precipitation excess, topography, slope and the permeability of the subsurface. It varies from 0 to 250 mm year-I. Apart from the meteoric recharge, a groundwater head dependent recharge/discharge is defined in these zones. This component of the recharge/discharge accounts for small local surface water systems which as such are not included in the model.

S. PalaeogeograpMeal reconstructions and boundary conditions

The impact of Pleistocene glaciations on the groundwater system is studied with the aid of an ice sheet model that simulates distributions of ice sheet thicknesses and basal melting

148

F.H.A. van Weert et al.LJoumal of Hydrology

195 (1997)

137-159

in space and time. The numerical ice sheet model is based on a mathematical description of mass and heat balances for the ice sheet (Boulton and Payne, 1993). It is driven by a transient climate function which is derived from the marine oxygen isotope Stage 6, correlated to the Saalian glaciation (Imbrie et al., 19&l), and constrained by geological data pertaining ice sheet extent, sea level changes and permafrost development. To study the effects of various stages of the Pleistocene glaciations on the geohydrology, a general classification of palaeogeographical situations has been made. These palaeogeographical situations are studied and reconstructed from litho-, morpho- and biostratigraphical data from the European continent. Furthermore, the resulting data on basal melt rates, ice sheet thicknesses and the ice sheet extent from the ice sheet model are used as input in the groundwater model for these different palaeogeographical situations. Six general situations have been distinguished:

rates

A-periods of warm climate conditions associated with high sea levels in the North Sea; B 1-periods of cold climate conditions associated with the presence of an ice sheet’ in Fennoscandia, low sea levels in the North Sea and discontinuous permafrost zones in the inland areas; B2-periods of cold climate conditions associated with the presence of an ice sheet in Fennoscandia, low sea levels in the North Sea and continuous permafrost zones in the inland areas; Cl-periods of cold climate conditions associated with the expansion of the confluencing Scottish and Fennoscandian ice sheets into lowland Northwest Europe; C2-periods of cold climate conditions associated with the expansion of the separated Scottish and FeMoscandian ice sheets into lowland Northwest Europe; C3-periods of cold climate conditions associated with the expansion of the Fennoscandian ice sheet into lowland Northwest Europe. Type A defines the interglacial (present day) situation. Type B 1 relates to deglaciation periods, while Type B2 relates to periods with advancing ice sheets. The Elsterian glaciation is mainly characterized by a Cl situation, the Saalian glaciation by Cl, C2 and C3 situations and the Weichselian glaciation by a C2 situation. Extensive areas of extraglacial discontinuous and/or continuous permafrost exist during the Type C 1, C2 and C3 situations. Fig. 4 shows the palaeogeographical map for each of the general situations. These maps contain data on the distribution of land, sea and ice, the occurrence of continuous, discontinuous and no-permafrost zones and the position of large river systems such as Thames, Scheldt, Meuse, Rhine, Weser and Elbe. Contours of the ground surface are given with an interval of 500 m. The figure also shows basal meltwater production rates. The topography in the different reconstructions is taken relative to the current sea level. In the reconstructions, sea level differences are considered. Isostatic reaction of the lithosphere because of glacial loading is also taken into account. Areas covered by ice will have subsided whereas areas more distant from the ice sheet may have uplifted in a glacial ‘foreland bulge’. After deglaciation a rebound effect has taken place. These isostatic responses are included in the contours and the land-sea distribution of the palaeogeographical reconstructions.

149

F.H.A. van Weert et al./Journal of Hydrology 195 (1997) 137-159

m

Imllmt

I

~mm(mm~)

J

TolpolImlPkY (m)

Fig. 4. Palaeogeographicalsituations representing various stages in the Pleistocene glaciations and the present day situation. One of the most important elements of the groundwater model are the boundary conditions. These conditions have to be applied to the sides of the model layers. They are based on the geometry of the NW European hydrogeological model and on the palaeogeographical reconstructions. Because both aquifers taper off in northern, eastern and southern directions the boundary conditions for these sides are assumed to be no-flow. No groundwater is entering or leaving the aquifers through these sides. Along the

150

F.H.A. van Weert et al./Journal o f Hydrology 195 (1997) 137-159

northwestern side the aquifers have considerable thicknesses. If this boundary is covered by sea, constant head boundary conditions are applied. In case of an ice cover, either a constant head or a no-flow boundary condition is used, depending on the ice sheet geometry. With the set of glacial and interglacial situations it is possible to construct time sequences covering the last three Pleistocene glaciations: Elsterian, Saalian and Weichselian. A time sequence is shown in Fig. 5. No consensus exists whether the Elsterian took place in isotopic Stage 12 or 10. Here, the Elsterian is in the isotopic Stage 10. Each of the palaeogeographical situations is assumed to represent different stadia in the Pleistocene record. For each of these situations a steady-state simulation of the groundwater system is carried out. Quasi time-dependent results are then obtained by creating a correct sequence of steady state simulations. 6. Simulation results and discussion To get some idea of the behaviour of the geohydrological system under extreme conditions (i.e. conditions which are totally different from the present day situation), a simulation was carried out for an extreme glacial situation. This is a Type C 1 situation with a maximum extent of the confluent Scottish and Fennoscandian ice sheets. The distributions of ice, land, sea and the permafrost zones for this situation are shown in Fig. 4. The distribution of the calculated steady-state groundwater heads in the upper and lower aquifer almost follows the same pattern, with some regional differences. These distributions are shown in Figs. 6 and 7. In large areas underneath the ice sheet the groundwater head in the upper aquifer is smaller than the hydraulic potential equivalent of the ice sheet thickness, and the total amount of basal meltwater will infiltrate in the aquifer. Hydraulic heads underneath the eastern lobe of the ice sheet are in general larger than underneath the western lobe, because larger amounts of meltwater have to infiltrate in a thinner aquifer. In parts of the aquifer close to the margin of the glacier, the hydraulic head is equal to, or is even slightly larger than the equivalent ice sheet potential. Because of the close relation between the hydraulic head and the ice sheet thickness, large hydraulic gradients exist in these areas. Drainage of the groundwater in proximal lakes, rivers and discontinuous permafrost zones increases these hydraulic gradients. Large gradients also occur where meltwater is forced into a thinning aquifer. In the lower aquifer, a relatively large hydraulic gradient is observed where the hydraulic conductivity changes from 1 to 0.1 m day -1. The distribution of groundwater recharge/discharge at the top of the model is shown in Fig. 8. In most areas underneath the ice sheet the subglacial produced meltwater infiltrates completely in the underlying sediments. In other subglacial areas the basal meltwater infiltrates only partly or even subglacial seepage of groundwater occurs. The hydraulic head approximately equals or even exceeds the equivalent ice sheet hydraulic potential in these areas. Such areas are located e.g. in northern Germany because the upper aquifer is very thin, in the current North Sea area because of converging groundwater flow from the western and eastern ice lobes and, east of Denmark because the meltwater has to infiltrate into the semi-pervious clay layer. These areas suggest zones where part of the subglacial meltwater has to be discharged

F.H.A. van Weert et aL/,lournalof Hydrology 195 (1997) 137-159

151

by channels or tunnels and subglacial lakes may be formed. The areas of possible tunnelling and channelling calculated by the model do not conflict with geological observations. Deep sediment fiUcd subglacial channels of the Elstcrian are very common in Northern Germany. Many eskers are found east of Denmark on the outcrops of the clayey aquitard. The effective pressure in these model areas may be very small. It is assumed that these areas correspond with zones favourable for shear failures and channel forming. Such areas

Global Chronology

Magneto-

stratigraphy

Ocean

Northwest Europe

Isotope stratigraphy

Chrono- I Climatostratigmphy

and pa.~u~gq~gn~cal sequence (oplion 1)

¢hro~ I

0

.~

K

benthordc,

6~eo (%0)

I.mlC~NE \

3

LATE

ehke

.~:,

\ .=

\

5

d__e

/

5

w,

8

\

/ 'a

MIDDLE

HOLSTEINIAN

<-

iO

ELSTIERIAN

\\ Blwa III

~

11

CROMERIAN COMPLEX

Palaeogeog~

sequence:

A: cool and win'm-temperate dirnate condi~ons relative high sea l e v ~ n North Se~, forest vegetation B: col0 climate condi~ns, ice sheet in Fennoscandia: 1) tltsconlnuous permaltost, 2) continuous permafrost J

C: cold climate conditions, ice sheet advance into NW European Iowkmds: 1) Scandina~an ice sheet In Noah SM area, 2) ~ e Scandinavian and Scottish Ice sheets In Norlh Sea area, 3) confluent Scandinavian and Scolesh ice sheets in Norlh Sea area

Fig. 5. Late Pleistocene chrono- and climatostratigraphic time series (based on isotopic variations) and the resulting palaeogeographical situation. The Elsterian is dated in isotopic Stage 10. The Isotope Stratigraphic data for the ocean (OPD677) were taken from Shacideton et al. (1990).

152

F.H.A. van Weert et al./Journal of Hydrology 195 (1997) 137-159

Hydraulic head in upper aquifer (meters)

-1oo 1(]o 300

5o0

7oo 9oo 11oo

?

0

]m

~0

NOlCm

Fig. 6. Hydraulic heads in the upper aquifer for the present day situation and various stages of the Pleistocene glaciations.

153

F.H.A. van Weert et alJJournal of Hydrology 195 (1997) 137-159

Hydraulichead in lower aqui~ (meters)

-lOO lOO 3o0

~o 700 9oo

$it~tiooC2 11oo

0

Im

sm

mini

Fig. 7. Hydraulic heads in the lower aquifer for the present day situation and various stages of the Pleistocene glaciations.

154

F.H.A. van Weert et alYJournal of Hydrology 195 (1997) 137-159

Vertical fluxes through top surface (millimeters per

year)

:, 100 15-100

0-15

[----]

0

0-15 15-100

100-1000

> 1000

S ~ m ~ C3

Fig. 8. Recharge/discharge distribution at the top surface for the present day situation and various stages of the Pleistocene glaciations.

F.H.A. van Ween et al./Jourml

of Hydrology 195 (1997) 137-159

155

observed in eastern parts of The Netherlands where thrusted ice-pushed ridges am found and in northern Germany where deep sediment-filled subglacial channels are found. Model results do not contradict these observations. The extensive system of rivers, the proglacial lakes in the south and the North Sea in the north mainly drain the highly pressurized groundwater system. The maximal groundwater discharge occurs directly at the snout of the glacier. Discharge rates in lakes and seas can be as large as 2000 mm year-‘. This value is corresponds to the value calculated by Boulton et al. (1993). He suggested a discharge rate of 10 000 mm year-’ at the Saalian glacier margin in The Netherlands. Drainage in the discontinuous permafrost zone is usually small because the hydraulic resistance is large. Without the draining facilities of the rivers, lakes and seas, subglacial hydraulic heads would be considerably larger, resulting in larger areas favourable for channelling, tunnelling and sediment deformation. In most parts of the lower Tertiary clay aquitard (Unit 2) vertical fluxes are small. The difference in hydraulic heads between the upper and lower aquifer is about 100 m at maximum. With an average hydraulic resistance of about lo6 years, the vertical flux amounts to 0.1 mm year-‘. Locally these fluxes can be much larger when the clayey aquitard is thin. Large vertical fluxes from bottom to top aquifer up to 100 mm year-’ have been found in the northern regions of the aquitard. The main objective of this study was to obtain information about the magnitude and direction of groundwater velocities under glacial conditions. At present, the velocities in the upper aquifer are in the range of a few to tens of metres per year. The magnitudes of the velocities in the lower aquifer are usually one or two orders of magnitude smaller because the hydraulic conductivities in this aquifer are one or two orders smaller. During extreme glacial conditions however, these velocities are much higher. Fig. 9 shows the actual horizontal groundwater velocities (m year-‘) in the upper aquifer as calculated for the present day situation and the glacial situations. According to this study, groundwater velocities in the upper aquifer of the NE Netherlands will be ten times larger under extreme glacial conditions than at present. Recent studies (Gostrom et al., 1993) concluded that groundwater velocities in the northern Netherlands during the Late Saalian glaciation were three times larger than present day velocities. Boulton et al. (1993) concluded that for the Netherlands lower aquifer groundwater velocities were one to two orders higher in glacial times compared with present. According to this groundwater model study extreme high groundwater velocities (up to hundreds of metres per year) may develop in the upper aquifer in northern Germany when the ice sheet has advanced to the southern edge of the aquifer. Possibly three different glacial periods occurred during the Pleistocene. Each of these glacial periods can be characterized by a sequence of palaeographical situations. The sequences for the different glacial periods are very similar: a Type A situation (interglacial), a Type B2 situation (advance of ice sheet), one or more Type C situations (ice covering NW Europe), a Type Bl situation (retreating ice sheet) and finally a Type A situation (interglacial). Differences between the different glacial periods are mainly in the combination of Type C situations and the duration of the different stadia. The starting point for each of the glacial cycles is an interglacial situation (Type A situation) which is very similar to the present day situation. In the early stages of the glaciation (Type B2 situation) the ice sheet does not yet cover the lowland areas. Because are

156

F.H.A. van Weert et al./Journal of Hydrology 195 (1997) 137-159

Horizontal groundwater velocities (meters per year) in upper aquifer Smmion 132

lhwmt I~y

0-10 10-50

50-100

>

$immtiuaC3

~

]

'

100

"

~mlgm Fig. 9. Horizontal groundwater velocities in the upper aquifer for the present day situation and various stages of the Pleistocene glaciations.

157

F.H.A. van Weert et alJJournal of Hydrology 195 (1997) 137-159

of a reduction in the precipitation and permafrost confinement, the groundwater recharge will decrease. Groundwater drainage decreases rapidly in this period due to the permanent and continuous freezing of the soil. Groundwater flow is still in the northwest direction towards the lowered North Sea. When the ice advances over the lowland area (Type C situation), the groundwater velocities will increase mainly due to a decrease of the drainage possibilities. The groundwater recharge in the extreme glacial situation (mainly subglacial meltwater infiltration) is relatively small compared with present values (10 to 20%). The aquifer becomes fully confined due to the extensive permafrost and ice cover of the area. Groundwater velocities are much larger than during the interglacial period in most parts of the model area. The groundwater flow pattern is totally changed. Groundwater flow underneath the ice sheet is directed along the glacial pressure gradient which is in a southern direction. In the most extreme case, the confluent ice sheets cut off the normal drainage pathways for the surface water and groundwater to the North Sea in the northwestern part of the model area. When the ice sheets retreat (Type B 1 situation), the groundwater velocity decreases gradually. During these thawing periods the permafrost decays and becomes discontinuous, allowing drainage of the pressurized groundwater system. The thawing of the massive ice sheets and permafrost results in large amounts of surface water. Recharge will be about 30-70% of the recharge during interglacial periods. Simulations have been carried out for all six characteristic situations. A time sequence for the period considered is created by replacing different stadia in the Pleistocene sequence by their representative characteristic palaeogeographic situation. Results of one of the six groundwater flow simulations are assumed to be representative for that period. Based on that approach, the magnitude of the hydraulic head in the lower aquifer along profile I-II (see Fig. 1) is given as a function of time. The results are shown in Fig. 10. It Hydraulic head (m) in lower aquifer along profile line I-.11 -50 . 0 o100 100. 200 200300 300400 400500 500750 750 - 1000 1000 . 1500

i

l

A

-loo [

FI

I [ _.I-.

~l r ~

-200

.... Q

"00

I i=

0

470000

-5OO

940(~0

Distance (m) along profile I-II Fig. 10. The temporalchangeof the horizontalgroundwatervelocitiesin the loweraquiferalongprofileI-II.

158

F.H.A. van Weert et alJJournal of Hydrology 195 (1997) 137-159

appears that the periods of high hydraulic heads are relatively short and last for 15000 years at maximum. The figure shows evidently that the hydraulic gradients and consequently the groundwater velocities are increased in the glacial periods.

7. Conclusions Based on the results of the modelling study to the Pleistocene glacial effects on the NW European geohydrological system the following conclusions can be drawn. During the Pleistocene glaciations subglacial meltwater infiltrates totally in the aquifer system in most parts of the ice covered area. In some parts of the model area the groundwater head equals or exceeds the equivalent ice sheet's hydraulic potential. Here, only part of the meltwater infiltrates or even subglacial groundwater discharge takes place. It is assumed that in the actual situation the excess meltwater discharges through a drainage system of channels, tunnels and lakes at the ice-bed interface in these areas. There is yet no consensus about how these processes work but discharge may occur by catastrophic events. Subglacial channels are found in the field which scour into the underlying sediments for hundreds of metres. Groundwater pressures in these areas may be very high, resulting in very small effective pressures in the sediments. Under such conditions brittle and shear deformation of the medium can occur, triggered by lateral lithostatic pressure gradients due to glacial loading. According to the results of this modelling study, these extreme hydrogeological responses are restricted to the upper aquifer in the northern parts of The Netherlands and Germany. These ar¢as correspond to regions where many sediment filled Pleistocene channel valleys and eskers are found. Under glacial conditions, the highly pressurized groundwater system is mainly drained by the system of proglacial lakes, rivers and zones of discontinuous permafrost. Without this drainage, the subglacial groundwater head would be much higher, resulting in larger areas where channelling, tunnelling, sediment deformation and erosion can occur. The groundwater velocities are relatively large under glacial conditions. When the glacier snout overlies an aquifer that tapers off, extreme high velocities ate found. Within the Pleistocene sequence, periods of high groundwater velocities are relatively short. However, the effect of high groundwater velocities may be substantial, since these will influence the integrity of disposal sites in geological formations like rock salt and clay. Furthermore, transport distances of radionuclides escaped from a repository during a period of high groundwater velocities may be considerable. The groundwater flow pattern, i.e. the direction of flow, under glacial conditions may be very different from the present day flow pattern.

Acknowledgements This study is carried out for the EC project: "Simulation of the effects of long-term climatic change on groundwater flow and the safety of geological disposal sites". The project is in the frame of the 4th CEC R&D Programme on Management and Storage of Radioactive Waste (1990-1994), Part 1A, Task 4: "Disposal of Radioactive Waste" CEC contract NO. FI2W-CT90-(X)46 (SSMA).

F.H.A. van Weert et al./Journal of Hydrology 195 (1997) 137-159

159

References Bouiton, G.S. and Jones, A.S., 1979. Stability of temperate ice caps and ice sheets resting on a deformable sediments. J. GlacioL, 24(90): 29-43. Bouiton, G.S. and Payne, A., 1993. Simulation of the European ice sheet through the last glacial cycle and prediction of future glaciation. SKB Technical Report 93-14, Swedish Nuclear Fuel and Waste Management Co., Stockholm, Sweden. Boulmn, G.S. and Hindmarsh, R.C.A., 1987. Sediment deformation beneath glaciers: rbeology and geological consequences. J. Geophys. Res., 92(B9): 9059-9082. Boulton, G.S., Slot, T., Blessing, K., Glasbergen, P., Leijnse, T. and van Cdjssel, K., 1993. Deep circulation of groundwater in ovefl~essured subglacial aquifers and its geological consequences. Quaternary Sci. Rev., 12: 739-745. Commission of the European Communities (CEC), 1982. Groundwater resources of the European Community. Synthetical Report, Directorate General for the Environment, Consmuex Protection and Nuclear Safety, Schafer, Hannover, Germany. Freuch, H.M., 1976. The Periglacial Environment. Longnum Inc., New York. Gibbard, P.L., 1988. The history of the great northwest European rivers during the past three million years. Phil. Trans. R. Soc. Lond., B318: 559-602. Imbrie, J., Hays, J.D., Mclntyre, A., Morley, D.J., Pisias, N.G., Prell, W. and ShacHeton, N.J., 1984. The orbital theory of Pleistocene climate; support from a revised chronology of the marine 6 tSO record. In: A. Berger, J.D. Imbrie, G. Kukla and B. Salznmn (Editors), Milanicovitch and Climate, D. Reidel, Hingham, MS, pp. 269-305. International Atomic Energy Agency (IAEA), 1981. Safety assessment for the underground disposal of radioactive waste. R~. 56, IAE, Vienna. Johmmsen, A., 1980. Hydrogeologie yon Schleswig-Holstein. Geologisch Jnhrbuch, Reihe C, Heft 28, 583 pp. Karsenberg, H., Deutloff, O. and Von Stempel, C., 1973. Report commission for hydrogeological maps in the AIH. BGR/Unesco, Hannover/Paris. Mouers, H.D., 1990. Ice-marginal thrusting of drift and bedrock: thermal regime, subglecial aquifers, and glacial surges. Can. J. Earth Sci., 27: 849-862. Oostmm, M., van Gijssel, K. and Ziffi W., 1993. Modelling subrosion and groundwater flow in the vicinity of the Zuidwending salt diapir on the basis of (palaeo)hydrological conditions for the Northeastern Netherlands. The SESAM Project, RIVM Report 715205004, National Institute of Public Health and Environmental Protection, Bilthoven, The Netherlands. Prij, J. and Benneker, P.BJ.M., 1989. Invloed van een tijdelijk temporamor daling en landljs-bedekking op de gesteentetemperatuur en -dmk. ECN R a p ~ No. ECN-87-113. Shackleton, N.J., Berger, A. and Peltier, W.R., 1990. An alternative astronomical calibration of the lower Pleistocene time scale based on ODP Site 677. Trans. R. Soc. Edinb.: Earth Sci., 81: 251-261. Shoemaker, E.M., 1986. Subglacial hydrology for an ice sheet resting on a deformable aquifer. J. Glaciol., 32(10): 20-30.

Shoemaker, E.M., 1992. Subglacial floods and the origin of low-relief ice-sheet lobes. J. Glaciol., 38(128): 105112. Zagwijn, W.H., Ruegg, G.H.J. and Cleveringa, P., 1992. Reconstructie paleohydrologische randvonrwaarden in bet Noord-Nederlandse OPLA gebied. RGD Report 30.107HRB4, Dutch Geological Survey, Haarlem, The Netherlands. Ziegier, P.A., 1991. Geological Atlas of Western and Central Europe. Shell/Elsevier, 's-Gravenhage, The Netherlands.