Groundwater–glacier meltwater interaction in proglacial aquifers

Groundwater plays a significant role in glacial hydrology and can buffer changes to the timing and magnitude of flows in meltwater rivers. However, proglacial aquifer characteristics or groundwater dynamics in glacial catchments are rarely studied directly. We provide direct evidence of proglacial groundwater storage, and quantify multi-year groundwater–meltwater dynamics, through detailed aquifer characterisation and intensive high-resolution monitoring of the proglacial system of a rapidly retreating glacier, Virkisjökull, in south-eastern Iceland. Proglacial unconsolidated glaciofluvial sediments comprise a highly permeable aquifer (25–40 m d⁻¹) in which groundwater flow in the shallowest 20–40 m of the aquifer is equivalent to 4.5 % (2.6 %–5.8 %) of mean river flow, and 9.7 % (5.8 %–12.3 %) of winter flow. Estimated annual groundwater flow through the entire aquifer thickness is 10 % (4 %–22 %) the magnitude of annual river flow. Groundwater in the aquifer is actively recharged by glacier meltwater and local precipitation, both rainfall and snowmelt, and strongly influenced by individual precipitation events. Local precipitation represents the highest proportion of recharge across the aquifer. However, significant glacial meltwater influence on groundwater within the aquifer occurs in a 50–500 m river zone within which there are complex groundwater–river exchanges. Stable isotopes, groundwater dynamics and temperature data demonstrate active recharge from river losses, especially in the summer melt season, with more than 25 % and often >50 % of groundwater in the near-river aquifer zone sourced from glacier meltwater. Proglacial aquifers such as these are common globally, and future changes in glacier coverage and precipitation are likely to increase the significance of groundwater storage within them. The scale of proglacial groundwater flow and storage has important implications for measuring meltwater flux, for predicting future river flows, and for providing strategic water supplies in de-glaciating catchments.

1 Introduction

A major challenge in modern hydrology is predicting changes in freshwater flows and storage resulting from glacier retreat in response to climate change (Jiménez Cisneros et al., 2014). Most glaciers worldwide have been in retreat since the mid-19th century, with the loss of global glacier ice accelerating during the 21st century (Zemp et al., 2015). This change has the potential to affect over 1 billion people who live in catchments where glacier melt contributes to river flow (Kundzewicz et al., 2008). Glacial retreat is expected to increase meltwater river flows until the mid-late 21st century (Jiménez Cisneros et al., 2014; Lutz et al., 2014). Longer term, as glacier ice loss continues, meltwater flows will decrease (Jiménez Cisneros et al., 2014). This lessening of the role of glaciers in regulating flows will change the nature of glacier-fed rivers and the importance of other water sources in glacier catchments: rainfall, snowmelt and groundwater. Predicted impacts include changes to the frequency and magnitude of flooding (Jiménez Cisneros et al., 2014); hydroelectric power production (Laghari, 2013); drinking water and irrigation (Kundzewicz et al., 2008); ecosystem functioning of catchments (Brown et al., 2007); and groundwater recharge (Taylor et al., 2013).

The role of groundwater storage in the hydrology of de-glaciating catchments is recognised, but to date there has been little direct hydrogeological investigation of groundwater–meltwater interactions (Levy et al., 2015; Vincent et al., 2019) with calls for more (Heckmann et al., 2016; Vincent et al., 2019). Indirect studies, inferred from river flow, indicate groundwater in Himalayan glacial catchments may be a significant source of delayed discharge to rivers (Andermann et al., 2012). Modelling of Himalayan catchments suggests that increased glacial melt this century will increase groundwater recharge from glacier runoff and the groundwater baseflow component in river flow (Immerzeel et al., 2013). Glacier meltwater rivers in Alaska can potentially lose half their annual flow to groundwater (Liljedahl et al., 2017). Groundwater can comprise 15 %–75 % of winter river flows in glacial catchments in the European Alps, Canadian Rockies, Peruvian Andes and Iceland (Malard et al., 1999; Hood et al., 2006; Bury et al., 2011; McKenzie et al., 2014; Baraer et al., 2015; MacDonald et al., 2016). Direct experimental studies of groundwater in glacial environments are rare (Vincent et al., 2019): e.g. subglacial groundwater behaviour (Sigurðsson, 1990; Boulton et al., 2001, 2007a, b); groundwater flow in relict rock glaciers (Winkler et al., 2016; Harrington et al., 2018); and the behaviour of shallow (<3 m) groundwater in glacial outwash plains in Iceland (Robinson et al., 2008, 2009a, b). The latter Icelandic studies demonstrated meltwater recharge to proglacial aquifers and linked retreating glaciers with declining groundwater levels (Levy et al., 2015).

In this study, we directly investigate the 3-D aquifer properties of a proglacial floodplain (referred to here as a sandur) of the Virkisjökull glacier in south-eastern Iceland, to 15 m depth, using geophysics, drilling and hydraulic conductivity testing; and provide continuous time series data for groundwater, river stage/flow and precipitation over 3 years, with campaign sampling for stable isotopes. We explore the relationships between groundwater, glacier meltwater flows and precipitation, revealing seasonal and spatial hydrological patterns.

Iceland provides an ideal observatory for studying groundwater in de-glaciating catchments. Ice melt from glaciers, which cover ∼11 % of Iceland, provides an estimated third of total river runoff (Björnson and Pálsson, 2008), but glacier retreat across Iceland (Sigurðsson et al., 2007) is forecast to produce significant changes in glacial catchment hydrology (Aðalgeirsdóttir et al., 2011). The British Geological Survey (BGS), in collaboration with Veðurstofa Íslands (the Icelandic Meteorological Office), has studied the Virkisjökull catchment since 2009, monitoring rapid glacier retreat (Bradwell et al., 2013), retreat mechanisms (Phillips et al., 2013, 2014), and researching glacial meltwater hydrology (MacDonald et al., 2016; Flett et al., 2017; Mackay et al., 2018).

Study site

Virkisjökull is an outlet glacier of the Vatnajökull ice cap in south-eastern Iceland (Fig. 1), within the Virkisá River basin, which has a catchment area of ∼32.5 km² to the confluence with the Svinafellsá River (MacDonald et al., 2016). Virkisjökull drains ice steeply south-westwards from an elevation of > 1800 m a.s.l. on the ice cap summit to <150 m a.s.l. at its terminus, with an average gradient of approximately 0.25. It has a high mass balance gradient, with net annual accumulation of more than 7 m w.e. a⁻¹ (metres of water equivalent per annum) at the ice cap summit (Guðmundsson, 2000) and net annual ice melt of more than 8 m w.e. a⁻¹ in the main ablation zone (Flett, 2016). The equilibrium line altitude on Virkisjökull is approximately 1150 m a.s.l. (MacDonald et al., 2016). The glacier has been retreating since 1990 (Hannesdóttir et al., 2015), with a marked acceleration in retreat rates since 2005 (Bradwell et al., 2013), during which time the glacier terminus has retreated by ∼1 km and there has been extensive surface lowering.

Figure 1 Virkisjökull study catchment. (a) Study area on Virkisjökull sandur, south-eastern Iceland, encompassing 6 km² groundwater catchment originating at a proglacial lake outlet. Hillshade model generated from LiDAR DEM ©Veðurstofa Íslands, 2010. (b) Piezometer M1 on the upper sandur near the catchment edge, showing established sandur vegetation and, in the middle distance, the area of moraines. (c) Virkisá River on the lower sandur in the summer melt season showing braided channels and an active, unvegetated sandur surface. Piezometers in the sandur are in three transects: Upper (U1–U2), Middle (M1–M3), and Lower (L1–L3).

The Virkisá River emerges from a small, shallow proglacial lake that has formed during the recent rapid deglaciation, and flows initially for 1 km over bedrock, flanked by moraines, and then for 4 km across the Virkisjökull sandur to the Svinafellsá River (Fig. 1). The river drains glacial meltwater and virtually all precipitation falling on Virkisjökull glacier, adjacent hillslopes and proglacial moraines. It occupies a single channel across the upper sandur, separating into a number of distinct channels across the lower sandur (Fig. 1). The mean summer river flow over 3 years of continuous monitoring (2011–2014) ranged from 5.3 to 7.9 m³ s⁻¹; and significant river flow occurred in winter (mean 1.6–2.4 m³ s⁻¹). Isotopic studies (MacDonald et al., 2016), validated by numerical modelling (Mackay et al., 2018), demonstrate that summer river flows are governed by glacier ice melt, and that winter flows are a combination of glacier meltwater, local precipitation and groundwater flow. The Virkisjökull sandur falls from 100 to 50 m a.s.l. with a surface gradient of 0.017 (Fig. 1). Over much of the sandur where river channels are actively migrating, there is little vegetation cover and no soil development. In more stable areas thin soils and more developed vegetation cover occur (Fig. 1). The groundwater catchment on the sandur associated with outflow from Virkisjökull has been estimated by using the surface water catchment identified from lidar and dGPS (Fig. 1).

The proglacial area has a maritime climate with cool summers (mean summer air temperature 8–12 ^∘C) and mild winters (1 ^∘C). Air temperature in the Virkisá basin is controlled mainly by altitude, with an average annual lapse rate of −5 ^∘C km⁻¹ (Flett, 2016; Mackay et al., 2018). Mean annual precipitation south-west of the Vatnajökull ice cap, including the Virkisjökull sandur, is ∼1800 mm; precipitation on the eastern side of the ice cap averages 3000 mm a⁻¹, and can exceed 7000 mm a⁻¹ on the ice cap summit (Guðmundsson, 2000). The proglacial area receives ∼150 precipitation days per year, estimated from interpretation of 3 years of daily photographs (MacDonald et al., 2016), which also show that snow cover, even in winter, rarely lasts for more than a week before melting. Potential evapotranspiration over the sandur was estimated at ∼450 mm a⁻¹ by Einarsson (1972) and actual evapotranspiration at 100–414 mm a⁻¹ by Jónsdóttir (2008).

2 Methodology

2.1 Aquifer characterisation

Eight boreholes were drilled into the sandur to 9–15 m depth during July and August 2012, in three transects approximately perpendicular to the river along a 3 km longitudinal reach in the upper, middle, and lower study catchment (Fig. 1a). Sediment samples collected during drilling were lithologically logged. The boreholes were installed as piezometers in September 2012, with 88 mm diameter uPVC plain casing to at least 5–12 m depth and a 3–6 m length of 0.5 mm slotted well screen below this (Table S1 in the Supplement). A further two boreholes were drilled into volcanic bedrock, to 5.5 and 13.75 m depth, between the glacier terminus and the upper edge of the sandur (Fig. 1a). Three methods were used to establish the physical aquifer properties of the sandur: (1) infiltration tests to 0.15 m depth at 20 locations, using a Guelph permeameter, and saturated hydraulic conductivity calculated by the Laplace method (Reynolds et al., 1983) (Table S2); (2) particle size analysis on 42 sandur sediment samples to 0.5 m depth, at 22 locations, and hydraulic conductivity estimated using a modified Hazen formula suitable for heterogeneous glacial deposits (MacDonald et al., 2012; Williams et al., 2019) (Table S3); and (3) constant rate pumping tests of between 3.5 and 6 h in each sandur piezometer, at rates of 0.5–1.8 L s⁻¹, and transmissivity estimated by the Jacob time-drawdown and Theis recovery methods corrected for unconfined conditions (Kruseman and de Ridder, 1994).

To measure aquifer thickness and depth to bedrock, two Tromino^® passive seismic surveys were undertaken transversely across the Virkisjökull sandur, and a third longitudinally down the Svinafellsandur aquifer 4.5 km to the west, using a single broad-band three-component seismometer with one vertical and two horizontal components. Measurements were recorded for 15 min at 50–100 m lateral intervals and data processed to derive depth to bedrock assuming typical shear wave velocities of 400–600 m s⁻¹ for Icelandic glacial sands and gravels (Bessason and Kaynia, 2002; Castellaro et al., 2005). These data were interpreted with a previous seismic reflection survey in the area to infer sediment thickness and potential layering (Guðmundsson et al., 2002).

2.2 Groundwater, surface hydrology, and precipitation monitoring

Monitoring of groundwater levels and temperature in sandur piezometers, at 15 min intervals, was undertaken from August 2012 to May 2015 (34 months) using In-Situ Inc. Rugged Troll 100 non-vented pressure transducers at 7–8.4 m depth (Ó Dochartaigh et al., 2019). Two In-Situ Rugged Barometer Trolls measured air pressure and temperature. River stage and discharge data are available for August 2012–May 2015 from an automatic stream gauge at Virkisá bridge (Fig. 1) with two water-level sensors, checked using daily photographs and continuous flow measurements from a radar mounted beneath a bridge (MacDonald et al., 2016). From April 2013 to March 2015 river stage and temperature were additionally monitored continuously every 15 min adjacent to piezometer U1 by an In-Situ Inc. Rugged Troll 100 pressure transducer (Fig. 1). Rainfall data and temperature for the proglacial area were measured at the closest of the three Automatic Weather Stations installed by BGS in the catchment (AWS1; 156 m a.s.l.). These weather stations were not equipped to measure snowfall, but daily photographs enabled periods of snowfall to be estimated. Long-term weather data from the Fagurhólsmýri weather station operated by the Icelandic Meteorological Office (IMO) approximately 12 km south of the study site, and national scale gridded products (Nawri et al., 2017), were used to check the plausibility of weather data measured on site.

Hierarchical cluster analysis of groundwater-level data was carried out on the entire dataset. Data were treated using the Standardized Groundwater level Index (Bloomfield and Marchant, 2013), which indicated the optimal number of clusters is four. Groundwater flow was estimated assuming a mean aquifer width of 1 km, an aquifer thickness at the river gauge from the passive seismic interpretation, an average measured groundwater-level gradient of 0.016, and hydraulic conductivity from the median of all measured values (n=64). Uncertainty was calculated from the interquartile range of measured K and uncertainty in aquifer thickness interpretation. The hydraulic conductivity of the deeper, unmeasured, sandur aquifer layer was estimated using the formula of MacDonald et al. (2012) taking into account a change in sediment state from very loose to loose and firm, which is likely to over-estimate the reduction in pore space due to loading (Schmidt and McDonald, 1979). The total volume of groundwater stored in the aquifer was estimated using a conservative estimate of the average aquifer porosity of 15 % (Parrieux and Nicoud, 1990).

2.3 Groundwater isotopic sampling and analysis

Physico-chemical analysis and modelling were based on samples of groundwater from piezometers and springs collected during three summer campaigns in September 2012, 2013, and 2014 and three winter (pre-melt) campaigns in January 2013, April 2013, and May 2014 (MacDonald et al., 2019). Groundwater sampling from piezometers was carried out after purging by low-flow pumping until stable readings were obtained for field-measured parameters. Field measurements of specific electrical conductance (SEC), temperature and bicarbonate alkalinity by titration pH (Table S4), and dissolved oxygen and redox potential (Eh) were made at the time of sampling. Samples for stable isotopes δ¹⁸O and δ²H were collected unfiltered in glass or Nalgene^™ polyethylene bottles and analysed at BGS laboratories by isotope ratio measurement on a VG-Micromass Optima mass spectrometer. Data are quoted in permil (‰) with respect to Vienna Standard Mean Ocean Water (VSMOW) (IAEA/WMO, 2019); measurement precision was ±0.1 ‰ for δ¹⁸O and ±1.0 ‰ for δ²H. Local precipitation stable isotope composition and a local meteoric line were estimated from International Atomic Energy Agency station data for Reykjavik (IAEA/WMO, 2019), supported by estimates for south-eastern Iceland (Arnason, 1977) and for southern Iceland (Sveinbjörnsdóttir et al., 1995), as described in MacDonald et al. (2016). The isotopic composition of the Virkisá River as it enters the sandur was established by sampling campaigns from September 2011 to December 2014 (MacDonald et al., 2016).

The high topographic gradient of the catchment, with large climatic differences between the upland glacial accumulation area (>1800 m a.s.l.) and the lowland temperate proglacial area (<150 m a.s.l.), results in two easily distinguished isotopic compositions: (1) glacier meltwater; and (2) precipitation across the proglacial area. A binary mixing model for δ²H was applied to investigate the relative contributions of local precipitation and of river water (which is dominated by glacier melt) to sandur groundwater, based on a two-component mixing equation. The end members applied for δ²H composition were −76.1 ‰ for river water (Table S4) and −58.5 ‰ for average annual local precipitation (MacDonald et al., 2016). The fraction of local precipitation in sandur groundwater (FGW) was calculated using the formula FGW $=({{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}_{\mathrm{R}}-{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}_{\mathrm{P}})/({{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}_{\mathrm{GW}}-{\mathit{\delta }}_{\mathrm{HR}}^{\mathrm{2}}$), where δ²H_R is the composition of river water; δ²H_P is the composition of local precipitation; and δ²H_GW is the composition of sampled groundwater. Since most river recharge to the aquifer occurs during the summer months when river flow is high and dominated by glacier melt, the impact of the small evolution in stable isotope composition down-river observed in winter due to groundwater baseflow (MacDonald et al., 2016) is insignificant, particularly when compared to the large difference between river flow and local precipitation isotopic composition.

3 Results

3.1 Sandur structure and aquifer properties

The groundwater study catchment covers 6 km² and encompasses the sandur, adjacent hillslopes and moraines, and river outflow from the proglacial lake (Fig. 1). Geophysical evidence from the passive seismic and previous seismic reflection survey indicates that (Figs. 2, S1) depth to bedrock increases from approximately 60–100 m in the upper sandur to 100–150 m in the lower sandur. The shallow aquifer material comprises loosely consolidated, moderately to poorly sorted, dominantly medium- to coarse-grained glaciofluvial sand, gravel, and cobbles (Fig. 2). All the sediment is of volcanic origin and has been transported and deposited by the Virkisá River. The deeper deposits are not exposed, but nearby seismic interpretation confirms that the material is generally uniform to >50 m, reflecting the similar sediment derivation and deposition mechanisms (Guðmundsson et al., 2002). Although not directly observed in the seismic data, there is a possibility that at greater depth (>50 m) there exist more consolidated Pleistocene aged sediments which have been compacted by ice loading during earlier glaciations (Guðmundsson et al., 2002). Observations of bedrock from nearby exposures and two boreholes drilled in bedrock reveal relatively massive and poorly fractured volcanic rock.

Figure 2 Geometry, geology, and hydrogeology of the sandur aquifer. (a) Hydraulic conductivity and summer groundwater-level contours. Other legend as Fig. 1. Hillshade model generated from LiDAR DEM © Veðurstofa Íslands, 2010; (b) schematic cross section of the hydrogeology, showing the locations of the piezometer transects, spring discharge area, and indicative groundwater flow lines; (c) geological section through the river bank showing heterogeneous glaciofluvial deposition; (d) perennial groundwater-fed stream on the lower sandur, associated with extensive growth of mosses and other aquatic vegetation; (e) groundwater discharge to the otherwise inactive river channel on the lower sandur, flowing to the active channel in the distance.

The sandur aquifer is highly permeable to at least 15 m depth, with a median hydraulic conductivity of 35 m d⁻¹ (IQ range 25–40 m d⁻¹) (Fig. 2a, Tables S2, S3). Transmissivity of the upper 15 m is 100–2500 m² d⁻¹ with a median value of 600 m² d⁻¹, consistent with hydraulic conductivity measurements (Table S1). The permeability of the deeper sandur aquifer was not directly measured. However, given the grain size distribution is the same as the shallow aquifer, and assuming the worst case of compaction due to burial and ice loading (Schmidt and McDonald, 1979), median hydraulic conductivity may have reduced to 15 m d⁻¹ or at a worst case 6 m d⁻¹ (MacDonald et al., 2012). By contrast, the underlying bedrock has low transmissivity, less than the lower limit from the experimental methods employed (transmissivity <0.25 m² d⁻¹). The sandur aquifer is unconfined. Depth to groundwater ranges from 0 to 4.4 m below ground level and maximum measured seasonal groundwater-level fluctuations are 1.0–3.6 m. From 1 km down-sandur from its upper edge, there is extensive groundwater discharge at the ground surface via perennial and ephemeral springs (Fig. 2). A conservative estimate of the volume of groundwater stored in the full thickness of the aquifer is 51 ±15 million m³, approximately 1 %–2 % of estimated ice volume in the glacier (Mackay et al., 2018).

3.2 Groundwater dynamics

Groundwater-level elevation falls from the upper to lower sandur, with a gradient of 0.018 across the upper and 0.013 across the lower sandur (Fig. 2). In the upper sandur, closest to the glacier, groundwater levels adjacent to the glacial meltwater channel are on average 1 m below river stage for most or all of the year (Fig. 3a, b), leading to a strong piezometric gradient away from the river to groundwater. Across the middle sandur, groundwater levels close to the river vary from 0.5 m below to 0.5 m above adjacent river stage, leading to complex river–groundwater interactions. Here, piezometric gradients are generally from river to aquifer in the summer melt season, when river flows are highest, and from aquifer to river in winter, driven by high winter precipitation and associated recharge. From 2 km down-sandur, groundwater levels are above the adjacent river stage for much of the year, creating a piezometric gradient that drives visible groundwater discharge through seeps and springs to the river (Fig. 2d) and ephemeral and perennial springs (Fig. 2e).

Figure 3 Groundwater levels, river stage, and precipitation. (a) Groundwater levels in the upper sandur during a 14 d dry period in summer (for the legend, see panel d). (b) Groundwater levels in the upper sandur during a 14 d rainy period in winter (for the legend, see panel d). Piezometer U1 (solid) is 20 m from the river; piezometer U2 (dashed) is 90 m from the river. (c) Dendogram obtained by hierarchical cluster analysis of groundwater-level data from sandur piezometers (piezometer locations in Fig. 1). The highest-level break shows two clusters representing piezometers where groundwater is influenced dominantly by local precipitation (U2, M1, M2, M3, L3) and piezometers where groundwater is influenced by the meltwater river (U1, L1, L2). Boxes show the four optimal sub-clusters indicated by data standardisation. (d) Monthly running mean of river stage (m), hourly precipitation (mm), and groundwater level (m), as a variation from the long-term average (LTA = 0).

Hierarchical cluster analysis of groundwater-level data indicates two patterns of groundwater-level fluctuation (Fig. 3c): one driven primarily by local precipitation; and the second driven partly by precipitation but also strongly influenced by river stage, especially in summer (Fig. 3d). Groundwater levels showing the first pattern (in piezometers U2, M1, M2, M3, and L3) fluctuate dominantly in response to individual precipitation events and longer-term precipitation patterns. The magnitude of groundwater-level fluctuations typically increases with distance from the river. Rainfall is higher than its long-term average throughout most of the winter and lower in summer, and this is generally reflected in the groundwater-level fluctuations (Fig. 3d). Groundwater levels showing the second pattern (in piezometers U1, L1, and L2) fluctuate in response to river stage as well as local precipitation, at seasonal (Fig. 3d) and also at event timescales (Fig. 3a). River stage is typically higher than its long-term average during peak summer melt, and groundwater levels in this group also remain close to or higher than their long-term average throughout the summer (Fig. 3d). The strongest response to river stage at a seasonal timescale is in piezometer L1, where groundwater levels during the 2013 summer melt season remained consistently higher than throughout the three winters from 2012 to 2014 (Fig. 3d). The strongest response at an event timescale is in piezometer U1, where groundwater levels show consistent diurnal fluctuations during the summer melt season that coincide with diurnal melt-controlled fluctuations in the river stage (Fig. 3a).

Piezometers U1 and U2 (20 and 90 m from the river, respectively) illustrate the relative impacts of summer glacier meltwater flows and large winter precipitation events on groundwater level–river stage gradients (Fig. 3). In summer, low precipitation and large glacier meltwater flows cause groundwater levels in U1 to rise above U2, creating a piezometric gradient away from the river (Fig. 3a). During individual winter rain storms, groundwater levels in U2 rise higher than in U1, creating a piezometric gradient towards the river (Fig. 3b) and driving baseflow to the river further downstream in the middle sandur.

Mean estimated annual groundwater flow through the shallow part of the aquifer calculated using Darcy's equation (20–40 m thick) is 0.19 m³ s⁻¹ (IQ range 0.093–0.30 m³ s⁻¹), equivalent to 4.5 % (2.7 %–5.8 %) of mean annual river flow and 9.7 % (5.8 %–12 %) of mean winter river flow. The relatively small seasonal variation in groundwater levels means there is no significant seasonal variation in estimated groundwater flow across the aquifer. Overall groundwater flow through the total depth of the sandur aquifer is estimated as 0.42 m³ s⁻¹ (0.12–1.1 m³ s⁻¹), equivalent to 9.8 % (3.6 %–22 %) of mean annual river flow and 21 % (7.7 %–46 %) of mean winter river flow.

3.3 Stable isotopes and temperature

Stable isotope composition (δ²H and δ¹⁸O) in groundwater from piezometers and springs was compared to that of glacier meltwater and local rainfall (Fig. 4, Table S4). Previous studies have demonstrated that glacial meltwater and local rainfall on the proglacial area are easily distinguished using δ²H and δ¹⁸O due to the high elevation of the accumulation area (MacDonald et al., 2016). Across individual piezometers, springs, and the river, variability between sampling campaigns was much less than variability between sites (Table S4). In particular, the river samples (taken as the river enters the sandur and therefore largely glacier meltwater) exhibited little seasonal variability $-\mathrm{76.1}\pm \mathrm{2.6}$ δ²H (n=19). Therefore, mean values from across the campaigns were used for analysis. Groundwater stable isotope compositions vary considerably across the sandur, spanning the range of compositions expected from glacier meltwater and local precipitation (Fig. 4). Piezometers (U1, L1, L2) identified from their hydrographs as most influenced by the river have isotopic compositions similar to river water, while piezometers whose hydrographs are influenced more by precipitation have a much wider range of isotopic composition, with U2 and M3 similar to local rainfall, and M2, M1, and L3 a mixture between local rainfall and river water. The springs showed a wide variety of compositions.

Figure 4 Stable isotope composition of waters and results of the binary mixing model of δ²H in groundwater. (a) Stable isotope composition of sandur groundwater in piezometers and springs. Individual piezometers labelled; piezometer locations in Fig. 1. Also shown are ranges in stable isotope composition of precipitation and river water (MacDonald et al., 2016). Plotted on the Local Meteoric Water Line (LMWL) calculated for Reykjavik. (b) Plot of the mean proportion of groundwater recharged from the river using binary mixing of the model of δ²H by perpendicular distance from the river and down-sandur. (c) Map of the mean proportion of groundwater estimated to be recharged from the river using a binary mixing model for δ²H. Hillshade model generated from LiDAR DEM © Veðurstofa Íslands, 2010.

A binary mixing model developed for δ²H indicates the relative proportion of precipitation and glacier meltwater in groundwater (Fig. 4b, c) and demonstrates a clear relationship with distance from the meltwater river. Within a zone extending up to 50 m from the river in the upper sandur, 130 m in the central sandur, and 500 m in the lower sandur, groundwater in piezometers generally comprises more than 50 % glacier meltwater. Shallower groundwater from springs within this river zone is more influenced by local precipitation, but still comprises more than 25 % glacier meltwater. Beyond this zone, groundwater from both piezometers and springs consistently comprises less than 25 % river water (Fig. 4c). Since the binary mixing model uses glacier meltwater as its endpoint, it is likely to be conservative in the proportion of river–groundwater interactions as it does not account for evolution of the river water stable isotope composition downstream due to groundwater baseflow. Selected hydrochemical tracers and water temperature also help distinguish these two zones (Table S4). Specific electrical conductance (SEC) and bicarbonate (HCO₃) are significantly lower in those piezometers strongly influenced by the river than those where precipitation influence is dominant (Fig. S1). River water temperature is relatively constant year-round at an average of 1.7 ^∘C, and mean annual groundwater temperature is lowest in piezometers close to the river and highest in those furthest from the river (Table S4).

4 Discussion

This study in Iceland shows that proglacial floodplains can form thick, highly permeable aquifers. By directly quantifying aquifer parameters and groundwater–glacier meltwater interaction, we have provided evidence of the significance of groundwater in proglacial hydrology. This has important implications for measuring glacial meltwater flux, for predicting future river flows and ecological impacts, and for water supplies in de-glaciating catchments. Similar thick proglacial glaciofluvial aquifers with high permeability and storage occur in other active glacial environments, e.g. elsewhere in Iceland (Robinson et al., 2008), the European Alps (Parrieux and Nicoud, 1990), and the Peruvian Andes (McKenzie et al., 2014), and with rapid deglaciation occurring globally proglacial aquifers are developing in many other locations, increasing the importance of characterising groundwater (Vincent et al., 2019).

4.1 Groundwater flow

Our study shows that significant water can flow through a glacierised catchment as groundwater, despite groundwater representing only a small proportion of the volume of water stored in glacial ice in the catchment. Reliable measurements of glacier meltwater are important for calibrating cryospheric-hydrological models (Bliss et al., 2014; Lutz et al., 2014; Mackay et al., 2018). The estimated volume of groundwater flow through the shallowest 20–40 m of the Virkisjökull proglacial aquifer is significant, 0.19 m³ s⁻¹, equivalent to approximately 4.5 % of mean annual river flow or 9.7 % of mean winter river flow, with estimates of 9.8 % and 21 %, respectively, if flow through the full thickness of the aquifer is considered. Other studies in Iceland have proposed that a similarly large proportion of meltwater (0.5–1 m w.e. a⁻¹) can flow through the groundwater system, either from sub-glacial or proglacial recharge (Sigurdsson, 1990; Hemmings et al., 2016); meltwater river losses to groundwater of up to 50 % have also been reported (Liljedahl et al., 2017). Measuring river flow in catchments with active glaciers is notoriously difficult, given the harsh conditions, actively changing river beds, and wide ranges in flows and sediment load. Measurements are therefore subject to high uncertainty. Here, we demonstrate that groundwater adds another source of uncertainty. Measurements of river flows that rely solely on river stage in the proglacial area are likely to underestimate total annual meltwater flows, with much higher relative errors at low flows. Similar potential underestimation in glacier melt estimations due to groundwater flow have recently been reported in South America (Saberi et al., 2019).

4.2 Meltwater–groundwater interaction

Groundwater–glacier meltwater interactions are controlled by relative differences in water levels between the river and the proglacial aquifer and vary both spatially, down the catchment, and seasonally. There is year-round active recharge of river water to the aquifer in the upper catchment, complex interaction in the middle of the sandur, and extensive groundwater baseflow to the river and springs across the lower catchment. Distinct patterns of groundwater–glacier meltwater dynamics are observed in groundwater-level fluctuations and in groundwater stable isotopic composition, temperature, and chemistry. In a zone extending up to 50–500 m from the river, the influence of the river on groundwater overshadows that of local precipitation. Here, recharge of glacier meltwater to the aquifer from river losses has a significant impact on the physical, chemical, and stable isotopic characteristics of groundwater in the proglacial aquifer. The aquifer provides additional water storage and groundwater discharges back to the river further downstream through a large number of springs and seeps (Figs. 1 and 4). This is consistent with other studies in glacier-dominated catchments, which inferred groundwater baseflow to rivers of 15 %–75 % (Malard et al., 1999; Hood et al., 2006; Bury et al., 2011;; McKenzie et al., 2014; MacDonald et al., 2016).

However, away from the river the aquifer is recharged dominantly from local precipitation. Active precipitation recharge to the aquifer is evident from groundwater stable isotopic composition and groundwater-level response to precipitation and reflects high annual precipitation (rainfall and snow), high aquifer permeability, and low evapotranspiration linked to limited soil development and vegetation cover. Recharge is likely to occur not only from direct precipitation on the sandur surface, but also from ephemeral streams draining from hillslopes and groundwater seepage from surrounding moraines. Groundwater discharge via springs and baseflow in the lower catchment supports surface water flows and local ecosystems and comprises groundwater derived mainly from local precipitation (Fig. 4).

Looking forward, as the glacier continues to melt, the proglacial aquifer will continue to have a buffering effect on river flow. High river flows will recharge the aquifer, whether caused by glacier ice melt, snowmelt or winter storms, as occurs in relic glacial outwash aquifers now in now-temperate areas (e.g. MacDonald et al., 2014), and will sustain springs, baseflow, and surface ecosystems further down the catchment. Local precipitation falling on the aquifer is likely to continue to be a major source of aquifer recharge and contribution to river baseflow in addition to the groundwater discharging from other glacial deposits emerging within the landscape (MacDonald et al., 2016). In upland areas in northern Europe where glaciofluvial deposits from past glaciations are present, detailed studies have demonstrated that groundwater often comprises more than 50 % of flow to river headwaters (Soulsby et al., 2005; Blumstock et al., 2015; Scheliga et al., 2017). Therefore, as glaciers continue to melt, groundwater baseflow is likely to become an increasingly important proportion of river flow in de-glaciating catchments.

4.3 Proglacial aquifers as strategic water resources

This study has demonstrated that the Virkisjökull sandur is a highly productive aquifer with regular recharge. Similar thick proglacial glaciofluvial aquifers occur throughout the world and are increasing in extent as glaciers recede, and are likely to also have the potential to sustain high-quality reliable water supplies. In formerly glaciated areas, these aquifers are often targeted for public water supply (e.g. Ó Dochartaigh et al., 2015) because of their ability to sustain high-yielding boreholes, their connectivity with rivers that provides additional recharge, and the generally high chemical quality of the groundwater compared to surface water. If projected glacier losses and increased precipitation in glacierised catchments are realised (Jiménez Cisneros et al., 2014), proglacial aquifers, recharged by local precipitation, represent a potentially significant store of high-quality water in regions around the world that currently rely on glacier melt for water supply.

5 Conclusions

Three years of investigations of groundwater and glacier meltwater at Virkisjökull, south-eastern Iceland, have enabled the aquifer parameters of the proglacial floodplain to be reliably characterised and seasonal groundwater–glacier meltwater dynamics to be quantified. The key findings from the research are the following.

1. Direct measurements of aquifer characteristics show consistently high permeability (35 m d⁻¹ n=64, IQR 25–40 m d⁻¹) and volume of groundwater storage (50±15 million m³). The proglacial floodplain therefore forms a highly productive aquifer.

2. Significant water flows as groundwater through the shallowest 20–40 m of the proglacial floodplain (0.19 m³ s⁻¹; IQ range 0.09–0.30 m³ s⁻¹), equivalent to 4.5 % of mean annual meltwater flow and 9.7 % of mean winter flow. If the full thickness of the aquifer is considered, then groundwater flows of 0.42 m³ s⁻¹ (IQ range 0.12–1.1 m³ s⁻¹) – equivalent to 9.8 % (3.6 %–22 %) of mean annual river flow and 21 % (7.7 %–46 %) of mean winter river flow – are possible.

3. Groundwater is recharged both from the glacial meltwater river and local precipitation falling on the aquifer or draining from nearby hillslopes. Glacier meltwater is particularly important in a zone from 50 to 500 m from the river, where glacier meltwater comprises >25 % and often >50 % of recharge.

4. There are complex but consistent river–groundwater interactions: in the upper sandur, closest to the glacier, the river loses to groundwater much of the year; in the middle sandur the river loses to the groundwater in the summer melt and gains from groundwater in the winter low flows; in the lower sandur groundwater provides baseflow to the river through springs and baseflow seeps.

Proglacial aquifers are common worldwide and increasing in extent with deglaciation. These findings, therefore, have wider implications for measuring glacier meltwater flux, for predicting future river flows, and for water supplies in de-glaciating catchments. Effectively understanding and characterising groundwater flows and storage in catchments with glaciers, and incorporating this in hydrological models, will strengthen our ability to predict and manage the hydrological and environmental impacts of accelerating glacier retreat.