1. Introduction
The hydrology of Vatnajökull is a subject of long-standing scientific interest and of practical importance for Iceland. It has attracted the attention of scientists studying jökulhlaups, or outburst floods (e.g. Reference BjörnssonBjörnsson, 1974, Reference Björnsson1976, Reference Björnsson1978,Reference Björnsson1988, Reference Björnsson1992; Reference NyeNye,1976; Reference Spring and HutterSpring and Hutter, 1981; Reference FowlerFowler, 1999; Reference Roberts, Russell, Tweed and KnudsenRoberts and others, 2000; Reference Björnsson, Pálsson, Gudmundsson and FlowersBjörnsson and others, 2001; Reference JóhannessonJóhannesson, 2002), glacier-bed mechanics (e.g. Reference Boulton and HindmarshBoulton and Hindmarsh, 1987) and fast-flowing ice (e.g. Reference BjörnssonBjörnsson, 1998; Reference Fuller and MurrayFuller and Murray, 2002). Vatnajökull feeds most of Iceland’s large rivers (Fig. 1) and is the source of the most destructive jökulhlaups in recent decades (e.g. Gudmundsson and others,1997).
Fig. 1. Location map of Vatnajökull and its major outlet rivers. Ice divides are shown to delineate major outlet glaciers.
Studies of Vatnajökull’s regional-scale hydrology over the last 15 years (e.g. Reference BjörnssonBjörnsson, 1982, Reference Björnsson1986a, Reference Björnsson1988) have relied on detailed topographic data and static subglacial fluid-potential calculations using the method of Reference ShreveShreve (1972, Reference Shreve1985). This “hydraulic-potential method” has been used to delineate subglacial drainage basins and predict flood routing from subglacial lakes in western Vatnajökull (e.g. Reference BjörnssonBjörnsson, 1978, Reference Björnsson1986a). While the hydraulic-potential method captures the first-order geometric controls on subglacial drainage, it takes no account of spatial and temporal variations in water supply. Furthermore, it does not include the potentially important effect of a permeable subglacial substrate. To construct a dynamic description of Vatnajökull hydrology, we apply recent advances in theoretical modelling along with Vatnajökull field data. In particular, we adapt the glaciohydraulic model of Reference Flowers and ClarkeFlowers and Clarke (2002) to Vatnajökull. Digital elevation models (DEMs) of the ice surface and bed (Reference BjörnssonBjörnsson, 1986b, Reference Björnsson1988; Science Institute, University of Iceland, unpublished data) define the glacier geometry, and the surface forcing is derived from 1999/2000 measured mass balance (Reference Pálsson, Björnsson, Eydal and HaraldssonPálsson and others, 2001). Subsurface hydrogeologic properties are estimated based on the geological map of Iceland (Reference Jóhannesson, Sæmundsson and JakobssonJóhannesson and others, 1990) and a review paper of Reference SigurðssonSigurðsson (1990).
The objectives of this study are (1) to compare the modelled subglacial drainage structure of Vatnajökull to that obtained by the hydraulic potential method, (2) to quantify the contribution of geothermally derived meltwater to glacier runoff, (3) to estimate glacial discharge through subsurface aquifers and (4) to constrain the glacial contributions to major rivers originating at Vatnajökull. This study forms the basis for an investigation of Vatnajökull geometry and hydrology in response to future climate, undertaken with both ice-dynamical and hydrological models.
2. Methods
We consider a two-layer coupled model for water flow at the glacier bed and in an underlying aquifer. Surface and basal water are sources to the subglacial drainage system, and water exchange between the subglacial and ground-water systems acts as a source or sink depending on the direction of water flow. Each system is treated as a vertically integrated layer governed by a local water balance (Reference Flowers and ClarkeFlowers and Clarke, 2002). In this section we briefly outline the model theory. Superscript s is reserved for variables associated with subglacial drainage, and superscript a for those associated with the aquifer.
2.1. Subglacial drainage
The ablation area of Vatnajökull drains largely through conduit networks. On rare occasions, most notably during jökulhlaups from Grímsvötn beneath Skeiðararjökull, these channels extend tens of kilometres into the centre of the ice cap. However, successful hydrological predictions have been made assuming subglacial water pressure is equivalent to the ice overburden pressure (e.g. Reference BjörnssonBjörnsson, 1982, Reference Björnsson1986a), suggesting that Vatnajökull drainage cannot be wholly conduit-dominated. More precisely, even in a conduit system the subglacial water pressure on a regional scale is likely to be higher than the pressure in the conduit itself, as conduits occupya small areal fraction of any given region and have a finite spatial influence (e.g. Reference Hubbard, Sharp, Willis, Nielsen and SmartHubbard and others, 1995). Where they exist, conduits would carry most of the discharge, but the basal hydraulic potential field would be jointly controlled by conduit flow and flow through the system feeding the conduits.
Following standard glaciological theory, we expect conduit networks to develop where and when there is a sufficient supply of surface water (see Reference PatersonPaterson, 1994, ch. 6, for a review). Hence the area serviced by conduits expands and contracts seasonally. Changes in the extent of the conduit network are mirrored by changes in a slower “distributed” system. These morphological transitions in the subglacial drainage system also occur in conjunction with surge onset and termination (e.g. Reference KambKamb and others, 1985; Reference BjörnssonBjörnsson, 1998). Such complexities present a substantial modelling challenge, especially as a conduit theory that is readily imported into continuum models has not yet been published.
One end-member formulation of basal drainage would be to construct a conduit network and solve the equations of conduit growth and decay in response to water input (e.g. Reference Arnold, Richards, Willis and SharpArnold and others, 1998). This approach requires a priori specification of the drainage-system structure — which is substantially more difficult for an entire ice cap than for an individual valley glacier — leaving only the temporal (rather than spatial) evolution of the drainage system undetermined. Whereas a continuum approach would allow some crude approximation of morphological transitions in the basal drainage system (from fast to slow, and vice versa), the conduit-based approach cannot emulate the slow distributed drainage that takes place over most of the ice cap unless it is coupled to a distributed model. The choice between these two end-member approaches depends on the goals of the study. For determining mid- to late-melt season discharge hydrographs, the conduit approach is superior despite the quantity of prior information required and its associated uncertainty. For investigating seasonal drainage-system transitions, objectively determining basal water routing, and relating basal hydraulics to dynamics with quantities such as effective pressure, the more flexible continuum approach is preferable.
Treatment of the subglacial system outlined below assumes that drainage takes place in a macroporous sediment horizon at the glacier bed (Reference Flowers and ClarkeFlowers and Clarke, 2002). This morphology lends itself to continuum mathematics, while enabling a wide range of hydrological behaviour. Changes in water flux are accommodated by adjustments in porosity, as is possible for dilatant materials such as glacial till (e.g. Reference ClarkeClarke, 1987). The glacier bed is assumed to be at the pressure-melting point. We introduce the notion of fast (conduitlike) and slow (distributed) drainage systems (Reference Raymond, Benedict, Harrison, Echelmeyer and SturmRaymond and others, 1995) by allowing subglacial hydraulic conductivity, the controlling parameter in this formulation of subglacial transport, to vary with subglacial water-sheet thickness.
For an incompressible fluid, conservation of mass leads to the subglacial water balance
where xj
is the horizontal spatial coordinate, h
s (x, y, t) is the subglacial water thickness, 
is the water flux, 
is the basal source term that includes melting due to geothermal heat and glacier sliding friction, 
is the surface water infiltration rate to the bed, and 
represents water exchange with the underlying aquifer. In this study, we do not have sufficient knowledge of basal sliding rates to include the frictional component of 
, so
is strictly a function of the geothermal heat flux Q
G, ice density ρ
I = 910 kg m−3 and latent heat of fusion of ice L = 3.34 × 105 J kg−1. We treat surface water as a direct source to the glacier bed, 
, rather than explicitly modelling its transport through the ice. This introduces substantial computational savings, and is justified insofar as water delivery to the bed is efficient. For a model with gridcells ≥ 1 km2 applied to the temperate and low-sloping ice cap Vatnajökull, efficient delivery of surface water through moulins and crevasses is a reasonable assumption below the equilibrium line. Surface water that originates above the equilibrium-line altitude (ELA) is routed to the appropriate model gridcell below the ELA in a predefined supraglacial drainage network. This network is derived using a steepest-descent algorithm that connects each interior pixel to an outlet point (e.g. Reference Zevenbergen and ThorneZevenbergen and Thorne, 1987; Reference TarbotonTarboton,1997). A similar strategy of surface water routing was used successfully by Reference Arnold, Richards, Willis and SharpArnold and others (1998) in a multicomponent hydrology model of Haut Glacier d’Arolla, Switzerland.
As for Darcian flow, we write the vertically integrated water flux in Equation (1) as
with a homogeneous and isotropic hydraulic conductivity K s and fluid potential ψ s = p s + ρ w g z B, where p s is water pressure, ρ w = 1000 kg m−3, g = 9.81 m s−2 and z B is bed elevation. Mathematical closure requires a relationship for p s as a function of h s.We use
(Reference FlowersFlowers, 2000), where p
I = ρ
I
g h
I with ice thickness h
I = z
S – z
B, and 
is the critical thickness of the water sheet (porosity × layer thickness) corresponding to buoyancy pressure (p
s = p
I ). We allow hydraulic conductivity to fluctuate in space and time as a function of h
s according to
where K = 1 m s−1 is introduced to non-dimensionalize the argument of the logarithm, k
a modulates the abruptness of the transition from 
to 
and k
b determines its position (Reference Flowers and ClarkeFlowers and Clarke, 2002).
2.2. Subsurface ground-water flow
The balance equation for the aquifer analogous to Equation (1) is
(Reference Flowers and ClarkeFlowers and Clarke, 2002). For a depth-independent aquifer porosity n a, the water thickness in the aquifer h a is defined as n a (z w – z L), where z w and z L are the elevations of the saturated horizon and the lower boundary of the aquifer, respectively. Fluid compressibility effects are included in the variable water density ρ a which obeys the equation of state
with β = 5.04 × 10−10 Pa−1 and p a the water pressure in the aquifer. Ground-water flux is written analogously to Equation (3) as
with fluid potential ψ
a = p
a + ρ
w
g z
L. Density ρ
w is used in Equation (8), rather than ρ
a, because Q
a appears as the argument of a spatial derivative in Equation (6) and we assume that 
.
Equation (6) governs both saturated and unsaturated flow in the aquifer. When the aquifer is unsaturated, the water table is a free surface, so (∂ρ a /∂t) = 0. Saturation implies that (∂h a /∂t) ≈ 0. For the unsaturated case,
with n a and d a the aquifer porosity and thickness, respectively. For the saturated case,
where α a is the aquifer compressibility and Equation (10) is derived for vertical stresses on an aquifer (see Reference Freeze and CherryFreeze and Cherry, 1979, p. 57) due to changes in water content (Reference Flowers and ClarkeFlowers and Clarke, 2002).
2.3. Water exchange
To avoid the complexities of a three-dimensional model, we parameterize the exchange between the subglacial sheet and underlying aquifer φ s:a. Cast in terms of dependent variables h s and h a, φ s:a evolves simultaneously with both systems. For saturated and unsaturated conditions in the aquifer, respectively,
where K t and d t are the vertical conductivity and thickness of the debris layer (aquitard) separating the basal hydraulic system from the aquifer. This is analogousto a system of interbedded aquifers and aquitards where flow is horizontal in the aquifers and vertical in the aquitards (Reference Bredehoeft and PinderBredehoeft and Pinder, 1970). The quantity ρ w g d t represents the driving potential due to elevation differences between the glacier bed and the top surface of the aquifer. For the purpose of calculating φ s:a, fluid pressure at the top surface of the aquifer is required. Thus, in the saturated case, the first term in Equation (10) is dropped and p a = (h a – n a d a)/α a d a. In the unsaturated case (Equation (9)), the top surface of the aquifer is unpressurized, so p a = 0. Both positive and negative values of φs:a are permitted, corresponding to aquifer infiltration and exfiltration, respectively.
2.4. Model numerics
For Vatnajökull the zonal (Lx ) and meridional (Ly ) limits of the model domain are approximately 18°20′–15°20′ W and 63°55′–65°00′ N, respectively. This area is discretized into nx × ny = 150 × 107 Cartesian gridcells, resulting in grids of size Δx = Δy = 1 km. Equations (1) and (6) are solved simultaneously on identical grids for 2 × nx × ny = 32100 unknown variables. Boundary conditions on the subglacial system (p s = 0) are imposed at the ice margin. Boundary conditions on the ground-water system are imposed distal from the ice cap at the grid edge. We assume the aquifer is well drained, hence h a = 0 at the grid boundary. Physical and numerical parameters are listed in Table 1.
Table 1. Physical constants and numerical parameters
Conventional finite-difference approximations are used to discretize the governing equations, with second-order centred differences in space and forward differences in time. We use a staggered numerical grid (see Reference PatankarPatankar, 1980) where scalar quantities (e.g. h
s) are evaluated at cell-centred nodes and vectors (e.g. 
) are evaluated across cell interfaces. Gradients of scalars then apply to cell interfaces and divergence of vectors to cell-centred nodes, which leads to a more realistic representation of the flux field and thus a more stable numerical scheme. Dependent variables h
s and h
a are represented as an equal blend of implicit (future) and explicit (present) values in all terms but the time derivatives. This partitioning is known as the Crank–Nicolson scheme, and serves to reduce the order of the error in the time discretization (Reference FletcherFletcher, 1991). A Newton–Krylov procedure is used to solve the resulting sparse system of equations. The problem is linearized in the outer iteration, and the linear system is solved with a preconditioned biconjugate gradient algorithm (Reference Press, Teukolsky, Vetterling and FlanneryPress and others, 1992, p.77–82).
3. Inputs and Implementation
Model initialization requires DEMs of the ice surface and bed, distributions of surface and basal melt rates and a variety of model parameters that define the vertical dimensions and properties of each system. Below we outline the sources and preprocessing of these inputs.
3.1. Ice-cap geometry
Radio-echo sounding data collected on Vatnajökull since 1976 have been processed and merged into high-resolution (200 × 200 m) DEMs of the glacier surface and bed (Reference BjörnssonBjörnsson, 1986b, Reference Björnsson1988; Reference Björnsson and EinarssonBjörnsson and Einarsson, 1991; Science Institute, University of Iceland, unpublished data). For this study we generate more computationally convenient DEMs by averaging each separate 5 × 5 pixel cluster into a single value. Figure 2 plots the resulting representation of Vatnajökull topography.
Fig. 2. DEMs (1 × 1 km pixels). (a) Glacier surface. Squares denote locations of automatic weather stations. (b) Glacier bed topography. Subglacial geothermal areas are labelled (KF, Kverkfjöll; SC, Skaftá Cauldrons; GV, Grímsvötn; PF, Pálsfjall). Dashed line partitions permeable (west) and impermeable (east) regions of the glacier bed.
The model domain is partitioned along a geological boundary (Fig. 2b) roughly separating permeable Quaternary basalts to the west from nearly impermeable Tertiary deposits to the east (Reference Jóhannesson, Sæmundsson and JakobssonJóhannesson and others, 1990). The permeability contrast between these formations is several up to eight orders of magnitude (personal communication from F. Sigurðsson, 2002). The ground-water equation (6) is only solved in the western (permeable) sector, where we assume the aquifer to be of uniform thickness and subparallel to the glacier bed.
3.2. Geothermal heat flux
Geothermal heat flux is a potentially important source of meltwater beneathVatnajökull. A heat-flux gradient striking northwest–southeast surrounds the ice cap, with the highest background fluxes (0.20–0.25 W m−2) overlying the active rift zone that divides Iceland (Reference Flóvenz and SæmundssonFlóvenz and Sæmundsson, 1993). Much higher localized values of heat flux (up to ∼50 W m−2) are found coincident with individual central volcanoes. In the absence of a subglacial heat-flux map, a uniform background value of 0.18 W m−2 is assigned to the eastern sector of Vatnajökull. In the western sector, subsurface hydrothermal circulation is believed to be sufficiently vigorous to prevent heat from reaching the base of the ice (personal communication from O. G. Flóvenz, 2001). Active geothermal areas beneath western Vatnajökull (Fig. 2b) are assigned values in W m−2 as follows(Reference BjörnssonBjörnsson,1988): Kverkfjoll,30; Skafta Cauldrons, 50; Grímsvötn/Grímsfjall, 50; Pálsfjall, 10. Geothermal heat flux Q G is related to the basal melt rate by Equation (2).
3.3. Mass balance, 1999–2000
Digital models of Vatnajökull winter (b w) and summer (b s ) balances for the glacial year 1999/2000 have been constructed to match the (200 × 200 m) DEM grid (Reference Pálsson, Björnsson, Eydal and HaraldssonPálsson and others, 2001). The net balance b n = b w + b s of the ice cap was −0.85 m for 1999/2000. Digital balance models with 150 × 107 pixels are created in the manner described for the DEMs.
From these balance distributions we compute the mean summer 
and annual 
surface melt rates as
with t s = 5 months (May–September) and t a = 1 year. Negative signs are used in Equations (12) and (13) to generate positive melt rates. These melt rates are minimum estimates as they do not include snow deposited and melted within the same measurement interval (e.g. May–September or October–April).
3.4. Surface-melt time series
In the absence of spatially dense measurements of surface melt rate, we utilize observed temperature records from summer 2000 to construct a time series of the surface melt rate for each model gridcell. The degree-day method (e.g. Reference Braithwaite, Olesen and OerlemansBraithwaite and Olesen, 1989), which assumes a proportionality between positive temperature and ablation, provides a ready means of estimating melt rates from air-temperature measurements. Provided we know the ablation totals of snow and ice and can construct a spatially variable temperature record, we need only assume a proportionality between the degree-day factors (DDFs) for snow and ice in order to estimate surface melt rates in each gridcell.
Six automatic weather stations (AWSs), ranging in elevation from 755 to 1724 m (see Fig. 2a), provide 2 m air temperatures for May–September 2000. We use splines to map these records onto a common time axis, and a linear regression to compute a lapse rate at each hourly time-step. All lapse rates are then averaged to obtain a single value for the whole interval. Six sea-level temperature time series are constructed by applying this calculated lapse rate (4.5°C km−1) to the original data. Figure 3 shows the result of averaging these time series after they have been adjusted to sea level (using the calculated lapse rate) to obtain a single synthetic sea-level temperature record. For modelling purposes, this reconstruction describes the temporal variation of air temperature. Spatial variations are introduced through the calculated lapse rate and the distribution of glacier surface elevations. Correspondence between observed and computed temperatures for the six AWS locations is reasonable, with standard deviations of 0.8–1.7°C.
Fig. 3. Synthetic sea-level temperature record, May–September 2000, constructed from five AWS air-temperature time series (station locations in Fig. 2a).
To compute time-dependent surface melt rates, we determine DDF values that produce the correct snow and ice ablation totals (A snow and A ice ) in each gridcell. For cells with A ice (i,j) = 0,
and for cells with A snow(i, j) = 0,
with Δt T = 3600 s and T (i, j, t) the reconstructed temperature vector of length N. A snow and A ice are computed as functions of the measured mass balance as
A snow is a minimum estimate of the total snow ablation since it does not include snow deposited and melted within the same mass-balance measurement interval.
For gridcells with A snow(i, j) > 0 and A ice (i, j) > 0, we take DDFsnow(i, j) = λ DDFice (i, j) with λ < 1. Assuming that snow melts to completion in any cell (i, j) before ice ablation begins,
where M is a time index that marks the transition from snow ablation to ice ablation. Setting DDFsnow(i, j) = λ DDFice(i, j) yields an equation in M:
We choose λ = 0.75 to match the highest value calculated from the literature (Reference Laumann and ReehLaumann and Reeh, 1993), to reflect the relatively mild climate (and hence warm snowfall) of Vatnajökull. For Sátujökull, an outlet of Hofsjökull ice cap, interior Iceland, Reference Jóhannesson, Sigurdsson, Laumann and KennettJóhannesson and others (1995) report values leading to λ = 0.73. M is then determined for each gridcell such that Equation (20) is satisfied, and DDFs are then back-calculated from Equations (18) and (19). These values are used in conjunction with the reconstructed sealevel temperature curve and calculated lapse rate to estimate the surface melt rate as a function of space and time.
4. Results and Discussion
We present steady-state simulations for the isolated subglacial system, followed by simulations that include geothermal heat flux, ground-water drainage and time-dependent surface ablation. Physical model parameters are listed in Table 2.
Table 2. Physical model parameters
4.1. Steady-state subglacial hydrology
The 1999/2000 mean annual surface-melt distribution is used to force steady-state simulations of Vatnajökull subglacial drainage. Meltwater above the ELA is routed over the glacier surface, as described earlier, before it is deposited at the glacier bed in the ablation area.
4.1.1. Subglacial flow field
The modelled subglacial flow field is plotted in Figure 4a (
, Equation (1)). Intensified subglacial drainage is suggested in basins with vigorous surface melt and large ablation areas, such as Brúarjökull and Skeiðararjökull. These two outlets, along with Dyngjujökull and Breiðamerkurjokull, are Vatnajökull’s most hydraulically active. Modelled flow convergence near the ice margin in Figure 4a highlights possible outlet locations for sub-basin-scale drainage networks. Nearly all of these implied outlets can be associated with actual river tributaries, except in areas with recent surface lava flows. Figure 4b is the hydraulic-potential gradient vector field with basal water pressure assumed to be equal to the ice overburden pressure: ∇ψ = ρ
I
g∇z
S + (ρ
w – ρ
I)g∇z
B, where ρ
I and ρ
w are the densities of ice and water, and z
S and z
B are glacier surface and bed elevations. This vector field is used to determine the subglacial catchment structure of Vatnajökull from which Reference Pálsson, Björnsson, Eydal and HaraldssonPálsson and others (2001) estimate glacier runoff to various rivers.While hydraulic-potential analysis is useful for identifying water divides, comparison between Figure 4a and b emphasizes the importance of a realistic basal water-pressure distribution (and hence a consideration of water supply) in imaging the large-scale subglacial drainage structure.
Fig. 4. Subglacial drainage structure. Interior contour is the equilibrium line. (a) Dynamic model steady-state. (b) Static model (hydraulic-potential method).Vectors are parallel and proportional to the fluid potential gradient.
4.1.2. Subglacial effective pressure
Effective pressure is defined as the difference between the ice overburden pressure and basal water pressure, p E = p I – p s. Because it affects basal shear stress and sediment shear strength, it is the hydrological metric most often used in laws for glacier sliding and bed deformation (e.g. Reference KambKamb, 1970; Reference Iverson, Baker, Hooke, Hanson and JanssonIverson and others, 1999; Reference KavanaughKavanaugh, 2000; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000). Figure 5 shows the simulated steady-state effective pressure for Vatnajökull.
Fig. 5. Simulated steady-state effective pressure pE = pI – ps expressed in metres of water.The pE = 0 contour is shown as a dashed line for the simulation forced by the mean annual surface melt rate, and as a solid line for a simulation forced by the summer surface melt rate.
Drainage basins of Brúarjökull, Skeiðararjökull and Breiðamerkurjokull contain significant areas of low effective pressure that expand in response to the increased summer melt rate. A common feature of these three outlets is low basal and surface topography, and, in the cases of Skeiðararjökull and Breiðamerkurjökull, deep bedrock trenches in the ablation area. The area enclosed by the dashed contour resembles the extent of Vatnajökull affected by surges of Brúarjökull, Skeiðararjökull and Breiðamerkurjökull (Reference ThorarinssonThórarinsson, 1969; Science Institute, University of Iceland, unpublished data). Areas affected by surging are smaller than the total ice catchment basins, especially in the case of Breiðamerkurjokull which only experiences surges in its eastern sector.
This is the first plausible and spatially variable subglacial effective-pressure distribution computed for Vatnajökull. Standard fluid potential calculations assume p s = p I and result in a uniform distribution p E = 0. Assuming p s ∝ p I but p s < p I yields an effective-pressure distribution that reflects the ice thickness and bears little resemblance to Figure 5.
4.1.3. Glacier runoff
Table 3 compares annual glacial discharge contributions to various outlet rivers as estimated with static and dynamic models. The static model (Reference Pálsson, Björnsson, Eydal and HaraldssonPálsson and others, 2001) uses measured mass balance and subglacial catchment structure, as determined by the hydraulic-potential method, to compute glacier runoff.The dynamic model, as described in this paper, makes no assumptions about subglacial catchment structure. For each river basin, the glacier runoff estimated with the dynamic model is the sum of simulated discharge from each ice-marginal gridcell.
Table 3. Comparison of 1999/2000 glacier runoff to selected outlet rivers as estimated with static (Reference Pálsson, Björnsson, Eydal and HaraldssonPálsson and others, 2001) and dynamic models
Agreement between statically and dynamically modelled discharges is good for Jökulsá á Breiðamerkursandi,Tungnaisá, Kaldakvísl andJökulsá í Fljótsdal. The Sylgja catchment area is poorly resolved with1 km2 gridcells, hence the dynamically modelled discharge is less reliable than for the larger outlets.
The discrepancy between statically and dynamically modelled discharges for Jökulsá á Brú and Kreppa/Kverká is strictly a function of catchment divides. Considered together, the discharge sum of the two rivers is approximately the same in both cases.
4.1.4. Contributions from basal melt
By including calculated basal melt rates in the simulations, we can evaluate the influence of geothermal heat flux on drainage under Vatnajökull, and comment on the provenance of outlet river discharge. This is relevant to hydrochemical studies of Vatnajökull (e.g. Reference Björnsson and KristmannsdóttirBjörnsson and Kristmannsdóttir, 1984), as water derived from the bed has a chemical signature distinct from that of surface water. Table 4 summarizes steady-state subglacial discharge to major outlet rivers as predicted by simulations with and without basal melt. Although geothermally derived meltwater is known to drain episodically in jökulhlaups from Vatnajökull, these simulations assume continuous meltwater drainage and hence highlight all possible routes for basal water. At least one study asserts that more basal meltwater drains from Vatnajökull than is accounted for in jökulhlaups (personal communication from F. Sigurðsson, 2001), raising questions as to the quantity and drainage routes of this additional water.
Table 4. Simulated glacial discharge contributions to selected outlet rivers with and without basal melt (1999/2000). Jökulhlaup physics is not included in the model, so simulations with basal melt represent continuous, rather than episodic, drainage of geothermally derived water
Although simulated discharge perturbations caused by introducing basal melt cannot be directly equated to geothermally derived discharge, they give some indication of how basal melt might be distributed among the outlets. To compute true provenance fractions would require a tracer transport model. The mean annual discharge due to surface melt from Vatnajökull is calculated to be ∼570 m3 s−1 for 1999/2000. Reference Björnsson, Pálsson, Guðmundsson and HaraldssonBjörnsson and others (1998) estimate the discharge for a zero net-balance year to be 486 m3 s−1.
For both annual and summer simulations, Jökulsá á Fjöllum experiences the largest absolute change in discharge, ∼8 m3 s−1, when basal melt is included. Tungnaá, Sylgja, Kaldakvísl and Jökulsá a Fjöllum all register annual differences over 10%. All rivers that appear to be appreciably affected by basal melt tap into geothermal systems beneath western Vatnajökull. Rivers draining the eastern sector (Jökulsá á Breiðamerkursandi, Jökulsá á Brú, Jökulsá í Fljótsdal) discharge negligible quantities of basal melt. Hence the background heat flux in the east is not particularly important for Vatnajökull hydrology, even on annual time-scales.
Figure 6 shows perturbations to subglacial hydraulichead profiles caused by geothermal heat sources. The effect of meltwater production at Grímsvötn (Fig. 6a) is localized over the geothermal area, except for a small difference in hydraulic head on upper Skeiðararjökull. Since we are not modelling subglacial lakes, the elevation of the ice over Grímsvötn is misrepresented in the profile, and should be higher by an amount roughly equivalent to the lake depth. A similar localization of basal meltwater is shown in the transverse profile (Fig. 6b) which intersects the Skaftá Cauldrons and north Grímsvötn.
Fig. 6. Profiles of hydraulic head with (solid line) and without (dashed line) geothermal heat sources for a steady-state simulation forced by the mean annual surface melt rate. Solid and dashed lines overlap except above geothermal areas. Glacier surface and bed topography are shown (bold lines) along with hydrostatic head (dotted line). (a) North–south transect through Grímsvötn (GV) and down Skeiðararjökull. (b) East–west transect through Skaftá Cauldrons (SC) and north Grímsvötn (GV).
4.2. Steady-state coupled hydrology
Springs in the vicinity of Vatnajökull, particularly north of Dyngjujökull (feeding Jökulsá á Fjöllum), attest to the presence of an active ground-water system (Reference SigurðssonSigurðsson, 1990). Glacially derived ground-water has been detected tens of kilometres from the ice margin, and ground-water recharge by Iceland’s four largest ice caps has been estimated at ∼130–220 m3 s−1 (Reference SigurðssonSigurðsson, 1990). Using the coupled subglacial–subsurface drainage model, we make independent estimates of the glacial contribution to ground-water drainage beneath western Vatnajökull.
A lack of complete hydrogeologic information for the area requires that we make the simplest possible assumptions. The aquifer is taken to be of uniform thickness and parallel to the glacier bed, and the hydraulic conductivity homogeneous and isotropic. Unless otherwise noted, simulations use parameters inTable 2 and are forced by the 1999/2000 surface melt rate.
Figure 7 shows steady-state distributions of water exchange with the aquifer (Equation (3)) and the corresponding ground-water flow field. The pattern of exchange in Figure 7a persists through simulations with a variety of parameter values. It is characterized by aquifer recharge beneath most of the glacier ablation area, and ground-water expulsion near the glacier margin and in the periglacial area. Bright areas are suggestive of springs and are to be found along the northern margin of Dyngjujökull, west of Köldukvíslarjökull and to the southwest of Tungnaárjökull and Síðujokull. The most vigorous predicted ground-water expulsion coincides with the river basins Skaftá and Hverfisfljót, but the latter is partially an artifact of the prescribed aquifer boundary. Compared to a map of known spring areas (see Reference SigurðssonSigurðsson, 1990, fig. 2) our results show artesian conditions over relatively short distances from the glacier margin. This is a result of modelling ground-water transport in homogeneous isotropic horizontal lava beds and neglecting the highly conductive anisotropic fissures.
Fig. 7. Ground-water reference equilibrium forced by the mean annual surface melt rate. (a) Modelled vertical water exchange with the glacier bed (ϕs:a, Equation (11)). Positive values represent aquifer recharge. (b) Modelled horizontal ground-water flow field.Vector lengths are proportional to discharge magnitude through the aquifer and are oriented parallel to the flow field.
In Figure 7a, the gray-scale limits (−5 to 10 m s−1) correspond to annual recharge rates of −1.6 to 3.2 m a−1, but the maximum modelled recharge rate is 7.6 m a−1. For comparison, Reference SigurðssonSigurðsson (1990) argues that the subglacial geology may permit recharge rates of 3–4 m a−1, but that rapid subglacial water transport in ice conduits probably limits these to 1.5–2 m a−1. Since our continuum model assumes extensive contact between basal water and the substrate, we regard the simulated recharge rates as maximum estimates. Integrating recharge rate over ice-cap area, we obtain a simulated total volumetric recharge rate beneath the glacier of 118 m3 s−1. This is consistent with Reference SigurðssonSigurðsson ’s (1990) estimate of 50–100 m3 s−1 of volumetric recharge in the basins of Tungnaá and Jökulsá á Fjöllum alone. The simulated total volumetric expulsion rate from the aquifer to the surface outside of the ice cap is 40 m3 s−1.
Despite the prescription of an isotropic aquifer conductivity, the simulated ground-water flow field (Fig. 7b) has a dominant northeast–southwest orientation roughly aligned with the geologic strike of the substrata. Modelled water transport in the aquifer is intensified beneath Dyngjujökull in the north and Tungnaárjökull and Síðujökull in the southwest, with topographically driven channelization occurring distal from the glacier. Modelled flow directions in the aquifer are similar to those at the glacier bed, but exhibit less small-scale spatial variability. For hydraulic conductivities ≤10−4 m s−1 in the aquifer, the model predicts ground-water transport to be negligible compared to subglacial drainage.
For different prescribed values of aquifer thickness d a, the spatial extent of ground-water expulsion to the surface would be expanded (for thin aquifers) or reduced (for thick aquifers). Increasing n a d a to 50 m would reduce the bright areas surrounding western Vatnajökull in Figure 7b to a small isolated patch along the southern margin (in front of Síðujökull). The ground-water flow field as shown in Figure 7b would appear more topographically channelized if the aquifer were thicker. Linear increases in the transport capacity of the aquifer (through increases in n a, d a or K a) result in a quadratic decline of the fraction of total runoff discharged from the ice–bed interface rather than through the subsurface. The sensitivity of individual drainage basins to changes in n a d a appears to be directly controlled by the subglacial catchment fractions underlain by permeable formations.
Simulated subglacial (Q s) and subsurface (Q a) discharges to the major outlet rivers are recorded in Table 5 for conservative and maximum estimates of aquifer properties. In the maximum model, five of twelve river basins are dominantly ground-water-fed (Q a /(Q a + Q s) > 0.5). In the conservative scenario, all western Vatnajökull river basins, from Jökulsá á Fjöllum in the north to Hverfisfljót in the south, have subsurface discharge fractions ≥0.15. Glacierwide, the aquifer transports ∼7% and ∼30%, respectively, of the total discharge in the conservative and maximum model simulations.
Table 5. Simulated 1999/2000 subglacial and ground-water discharge to various rivers for two ground-water scenarios: conservative (na da = 25 m, Ka = 10−3 m s−1) and maximum (na da = 50 m, Ka = 10−2 m s−1). Results represent a steady state forced by the mean annual surface melt rate
Reference SigurðssonSigurðsson (1990) illustrates the effect of ground-water drainage beneath Vatnajökull by comparing Jökulsá á Fjöllum and Jökulsá á Brù, whose drainage basins are situated on permeable and impermeable formations, respectively. He notes that both rivers carry approximately 150 m3 s−1 when leaving the glacially influenced ground-water basins, but that the winter base flow in Jökulsá á Brù is ∼5–10 m3 s−1 compared to ∼50 m3 s−1 in Jökulsá á Fjöllum. Our reference model simulations for 2000 predict summer runoff rates from Jökulsá á Fjöllum and Jökulsá á Brù to be 154 and 181 m3 s−1, respectively, and annual runoff rates to be approximately equal around 90 m3 s−1. The simulated mean annual ground-water discharge, which can be compared to Reference SigurðssonSigurðsson ’s (1990) winter base flow estimate, is ∼5 m3 s−1 for Jökulsá á Brù and ∼60 m3 s−1 for Jökulsá á Fjöllum.
These simulations provide only a rough estimate of the large-scale drainage budget of Vatnajökull. Ground-water drainage in the vicinity of Vatnajökull is undoubtedly affected by spatial variations in hydrogeologic properties, which have been entirely neglected in this analysis. The most important future development of this aspect of the modelling is to introduce an anisotropic hydraulic conductivity tensor; in the limit that no ground-water transport is permitted perpendicular to the geologic strike, subsurface catchment structure would be profoundly affected.
4.3. Time-dependent hydrology
Realistic melt-season representations of Vatnajökull hydrology require time-evolving forcing and response. Making use of the temperature record in Figure 3 and DDFs computed in Equations (18) and (19), we can simulate spatial evolution of the subglacial drainage system (Fig. 8). For simplicity, this simulation does not include ground-water drainage or meltwater produced by geothermal heat sources beneath western Vatnajökull, and assumes that water originating above the ELA does not reach the bed.
Fig. 8. Snapshots of simulated subglacial drainage structure. Snowline is contoured in bold. (a) 31 May 2000 (day 152). (b) 30 July 2000 (day 212).
In all time-dependent simulations, we find large changes in the Brúarjökull catchment area in response to surface meltwater input. In early summer, the Brúarjökull subglacial catchment area extends ∼50 km interior from the northern glacier margin (Fig. 8a). As the snowline retreats, the fluid potential gradient created by incoming meltwater is sufficient to override the gentle elevation potential gradient. Consequently, water above the snowline is driven uphill and southward toward Breiðamerkurjökull and the hydraulic divide migrates in time (Fig. 8b). While this may be an artifact of oversimplifying the basal hydrology, we cannot dismiss it entirely. We do not see any decisive indication of such divide migration in the western half of Vatnajökull where it could have serious implications for jokulhlaup routing. Surface elevation gradients are more severe in the west, so that meltwater-induced divide perturbations are less likely.
4.3.1. Discharge hydrographs, May–September 2000
Introducing a time-dependent surface melt rate permits the first estimates of seasonally varying glacial discharge contributions to Vatnajökull’s outlet rivers. Figure 9 shows the mean daily measured discharge in five rivers (data of the National Power Company of Iceland) compared to the simulated glacial discharge. This crude approach to surface hydrology enables us to bracket a range of possibilities, while obviating the need for a complicated physical model.
Fig. 9. Simulated glacial discharge (fine lines) compared to measured hydrographs (bold lines) for five outlet rivers, May–September 2000.The dashed lines represent end-member simulations, one assuming that all snowmelt and ice melt reaches the bed (upper curve) and one assuming that only ice melt reaches the bed (lower curve). Fine lines represent a compromise simulation in which snowmelt above the ELA is stored on the glacier surface. Data provided by H. H. Haraldsson, National Power Company of Iceland. (a) Tungnaá; (b) Jökulsá á Fjöllum; (c) Kreppa/Kverká; (d) Jökulsá á Brù; (e) Jökulsá í Fljótsdal.
Gauge measurements include significant non-glacial runoff and are variously affected by ground-water drainage. Hence, an unknown fraction of the total basin discharge is unaccounted for in the measurements (except in the case of Jökulsá í Fljótsdal). This uncertainty is probably largest for Jökulsá á Fjöllum. Calculations of glacial runoff by Reference Björnsson, Pálsson, Guðmundsson and HaraldssonBjörnsson and others (1998) for 1994/95 yielded mean values for Jökulsá á Fjöllum and Kaldakvísl that were higher than the gauge measurements.
Measured hydrographs for the smaller rivers Tungnaá, Kreppa/Kverká and Jökulsá í Fljótsdal have early-season snowmelt peaks (before day 150) that are comparable to or larger than peaks due to glacier ablation. In contrast, hydrographs for the large rivers Jökulsá á Fjöllum and Jökulsá á Brú are dominated by late-season discharge. Except in the river Tungnaá, the largest amplitude fluctuations in both simulated and observed hydrographs occur on seasonal time-scales, followed by fluctuations on 10–20 day time-scales. Fluctuations on the diurnal time-scale are by far the smallest for both simulated and observed records (personal communication from H. H. Haraldsson, 2001), though it is not shown in the measurements. For all but Jökulsá á Brá, the glacial component of basin discharge generally increases from early to late melt season. The variability structure in the simulated records is similar to that in the observations, with most simulated peaks appearing to correlate with peaks in the observations. This is not necessarily expected, considering the simplistic method used to construct the surface forcing and the fact that summer precipitation is not included in the model.
Hydrographs simulated using the minimum model (lower dashed line in each panel of Fig. 9) show little early-season glacial discharge because we assume snowmelt does not contribute to runoff. Maximum and compromise models agree well in the early season, because they differ only in their treatment of snow above the ELA. Minimum and compromise models converge in the late season, when most of the runoff is derived from ice ablation. The discrepancy between end-member simulations for each river reveals the sensitivity of the respective drainage basins to snow hydrology, with Jökulsá á Fjöllum being the most sensitive to snowmelt routing, and Jökulsá á Brù the least. This sensitivity is a direct reflection of glacier surface hypsometry, with more sensitive basins having a larger fraction of their total area at high elevation.
To gauge the effect of ground-water drainage on glacier discharge as presented in Figure 9, we briefly reintroduce the coupled model. The simulated glacier discharge fraction routed through the aquifer (Q a /(Q a + Q s)) is plotted in Figure 10 for ground-water maximum (n a d a = 50 m, K a = 10−2 m s−1) and conservative (n a d a = 25 m, K a = 10−3 m s−1) models. The quantity Q a /(Q a + Q s) is predicted by the model to be highest in the early season and lowest during the peak melt season between days 180 and 250. The prominent dip preceding day 140 represents the arrival of snowmelt. In the conservative case (lower curves, Fig. 10), ground-water discharge comprises >60% of the early-season glacier runoff in Jökulsá á Fjöllum and Jökulsá á Brù, and less than ∼10% during the peak melt season for all basins. With the ground-water maximum model, Q a /(Q a + Q s) is generally higher and more variable, especially in Jökulsá á Fjöllum where ground-water comprises 25–100% of the discharge. Determination of this ground-water component is prerequisite to a rigorous quantification of glacially derived discharge in rivers fed by Vatnajökull. This is a challenging task given the potentially continuous exchange of surface and subsurface water in these drainage basins.
Fig. 10. Simulated fraction of glacial runoff routed through the aquifer, May–September 2000, for maximum (upper curves) and conservative (lower curves) ground-water models. (a) Tungnaá; (b) Jökulsá á Fjöllum; (c) Kreppa/Kverká; (d) Jökulsá á Brù.
5. Summary and Conclusions
A simple two-layer physical model of drainage under Vatnajökull has facilitated new insight into its basal hydrology and water budget. Our results complement existing studies of Vatnajökull by providing a physical link between the surface mass-balance distribution and glacier runoff.
Steady-state simulations of Vatnajökull suggest a summer 2000 mean discharge of ∼1350 m3 s−1, compared to a 1999/2000 annual mean of ∼570 m3 s−1. Geothermally derived meltwater is predicted by the model to account for approximately 5% of the mean annual glacier discharge. Modelled discharge provenance of Vatnajökull’s major outlet rivers is roughly consistent with prior estimates based on mass balance and a statically derived subglacial catchment structure.While this reflects the similarity in subglacial flow directions implied by the two methods, the overall basal drainage structures appear markedly different because the distribution of surface meltwater is accounted for in the physical model.
Brúarjökull, Skeiðararjökull and Breiðamerkurjökull appear to be the most hydrologically active drainage basins. All three glaciers occupy low-lying areas of the bed and have gently sloping surfaces. Simulations predict an early and ultimately extensive snowline retreat on these glaciers, with correspondingly high basal water fluxes. Low subglacial effective pressures are expected in these basins, at least prior to the seasonal development of conduit networks, suggesting favourable conditions for glacier sliding and sediment deformation. In steady-state simulations forced by 2000 mean summer melt rates, the zero effective-pressure contour resembles the perimeter of surge-affected areas in these three glacier basins. Dyngjujökull (feeding the river Jökulsá á Fjöllum) appears to be the most hydrologically complex drainage basin, exhibiting the highest sensitivities to geothermally derived meltwater, subsurface ground-water loss and surface snowmelt routing.
We estimate that the ground-water system beneath western Vatnajökull evacuates up to ∼30% of glacier discharge annually. In drainage basins feeding the rivers Hverfisfljót, Skaftá, Tungnaá, Kaldakvísl and Jökulsá á Fjöllum, ground-water may account for 70–80% of the annual discharge. Although hydraulic anisotropy is neglected in this study, regional-scale ground-water transport is nevertheless aligned with the southwest–northeast strike of subglacial geologic formations. Using reference model parameters, we calculate the mean annual ground-water discharge from Vatnajökull for the balance year 2000 to be ∼145 m3 s−1, which falls within Reference SigurðssonSigurðsson’s (1990) estimate of 130–220 m3 s−1 for the total ground-water runoff from Iceland’s four largest ice caps.
Our results highlight possible topics of interest for future field investigations of Vatnajökull, including (1) discharge with chemical signatures indicative of basally derived meltwater in rivers other than Skaftá and Skeiðará, (2) indications of ground-water emergence beneath glacier margins in westernVatnajökull, and (3) evidence for subglacial hydraulic divide migration between Brúarjökull and Breiðamerkurjokull during the melt season. Hydraulic divide stability in western Vatnajökull, where jökulhlaups initiate, remains an open question. To address this issue and its implications for jokulhlaup routing requires a model that includes a better subgrid representation of drainage-system evolution and a proper treatment of subglacial lakes.
Acknowledgements
G. Flowers gratefully acknowledges the Science Institute, University of Iceland, for its hospitality during her postdoctoral tenure there. She was supported by the U.S. National Science Foundation under the International Research Fellow Awards Program (award No. 0000-425). Mass-balance data collection was sponsored by the National Power Company and the Road Authority of Iceland, in conjunction with the European Union (Framework IV — Environment and Climate, TEMBA 1996–97, ICEMASS 1998–2000). Radio-echo data have been collected over many years, with funding from the National Power Company, Icelandic Road Authority, Icelandic Research Council and Research Fund of the University of Iceland. H. H. Haraldsson of the National Power Company of Iceland provided all the river discharge data in this paper. We wish to acknowledge the field volunteers, especially the members of the Icelandic Glaciological Society. The manuscript benefited from the careful reviews ofT. Schuler, A.G. Fountain and J. S.Walder.
