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−1) 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.


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-19 th century, with the loss of global glacier ice accelerating during the 21 st century (Zemp et al., 2015). This change has the potential to affect over one 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 5 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 deglaciating catchments is recognized, but to date there has been little 10 direct hydrogeological investigation of groundwater-meltwater interactions (Levy et al., 2015, Vincent et al., 2018 with calls for more (Heckmann et al., 2016, Vincent et al., 2018. 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 15 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. 2018): e.g. subglacial groundwater behaviour (Sigurðsson, 1990;Boulton et al., 2001;Boulton et al., 2007aBoulton et al., , 2007b; groundwater flow in relict rock glaciers (Winkler et al., 20 2016;Harrington et al., 2018); and the behaviour of shallow (<3 m) groundwater in glacial outwash plains in Iceland (Robinson et al., 2008;Robinson, et al., 2009a;Robinson, et al., 2009b). 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 3D aquifer properties of a proglacial floodplain (referred to here as sandur) of the Virkisjökull glacier in SE Iceland, to 15 m depth, using geophysics, drilling and hydraulic conductivity testing; and provide 25 continuous time series data for groundwater, river stage/flow and precipitation over three years, with campaign sampling for stable isotopes. We explore the relationships between groundwater, glacier meltwater flows and precipitation, revealing seasonal and spatial hydrological patterns.

Study site
Virkisjökull is an outlet glacier of the Vatnajökull ice cap in SE Iceland (Figure 1), within the Virkisá river basin, which has a catchment area of ~32.5 km 2 to the confluence with the Svinafellsá river (MacDonald et al., 2016). Virkisjökull drains ice 5 steeply southwestwards from an elevation of >1800 m asl on the ice cap summit to <150 m asl 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 -1 (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 -1 in the main ablation zone (Flett, 2016). The equilibrium line altitude on Virkisjökull is approximately 1150 m asl (MacDonald et al., 2016). The glacier has been retreating since 1990 (Hannesdóttir et al., 2015), with a marked acceleration 10 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.
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 ( Figure 1). The river drains glacial meltwater and virtually all precipitation falling on Virkisjökull glacier, adjacent hillslopes 15 and proglacial moraines. It occupies a single channel across the upper sandur, separating into a number of distinct channels across the lower sandur ( Figure 1). The mean summer river flow over three years of continuous monitoring (2011)(2012)(2013)(2014) ranged from 5.3-7.9 m 3 s -1 ; and significant river flow occurred in winter (mean 1.6-2.4 m 3 s -1 ). 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 20 Virkisjökull sandur falls from 100 to 50 m asl with a surface gradient of 0.017 ( Figure 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 ( Figure 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 ( Figure 1).
The proglacial area has a maritime climate with cool summers (mean summer air temperature 8-12 °C) and mild winters (1 25 °C). Air temperature in the Virkisá basin is controlled mainly by altitude, with an average annual lapse rate of -5 °C km -1 (Flett, 2016;Mackay et al., 2018). Mean annual precipitation southwest 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 -1 , and can exceed 7000 mm a -1 on the ice cap summit (Guðmundsson, 2000). The proglacial area receives ~150 precipitation days per year, estimated from interpretation of three years of daily photographs (MacDonald et al., 2016), which also show that snow cover, even in winter, 30 rarely lasts for more than a week before melting. Potential evapotranspiration over the sandur was estimated at ~450 mm a -1 by Einarsson (1972) and actual evapotranspiration at 100-414 mm a -1 by Jónsdóttir (2008).

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 (Figure 1a). Sediment 5 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). 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 (Figure 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 10 calculated by the Laplace method (Reynolds et al., 1983) (Table S2) (Table S3); and (3) constant rate pumping tests of between 3.5 and 6 hours in each sandur piezometer, at rates of 0.5-1.8 l s -1 , and transmissivity estimated by the Jacob timedrawdown and Theis Recovery methods corrected for unconfined conditions (Kruseman and de Ridder, 1994). 15 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 minutes at 50-100 m lateral intervals and data processed to derive depth to bedrock assuming typical shear wave velocities of 400-600 m s -1 for Icelandic glacial sands and gravels (Bessason and Kaynia, 2002;Castellaro et al., 2005). These data were 20 interpreted with a previous seismic reflection survey in the area to infer sediment thickness and potential layering (Guðmundsson, 2002).

Groundwater, surface hydrology and precipitation monitoring
Monitoring of groundwater levels and temperature in sandur piezometers, at 15 minutes intervals, was undertaken from August 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 asl). 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. 5 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, aquifer thickness at the river gauge from the passive seismic interpretation, average measured groundwater level gradient of 0.016 and hydraulic conductivity from median of all measured values (n = 64). Uncertainty was calculated from the interquartile range of measured K and uncertainty in 10 aquifer thickness interpretation. The hydraulic conductivity of the deeper, unmeasured, sandur aquifer layer were 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 average aquifer porosity of 15% (Parrieux & Nicoud, 1990). 15

Groundwater isotopic sampling and analysis
Physico-chemical analysis and modelling were based on samples of groundwater from piezometers and springs collected  (Table S4), and of dissolved oxygen and redox potential (Eh), were made at the time of sampling. Samples for stable isotopes δ 18 O and δ 2 H were collected unfiltered in glass or Nalgene™ polyethelyne bottles and analyzed 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, 2016); measurement precision was ±0.1‰ for δ 18 O and ±1.0‰ for δ 2 H. Local precipitation stable isotope 25 composition and a local meteoric line were estimated from International Atomic Energy Agency station data for Reykjavik (IAEA/WMO, 2016), supported by estimates for southeast Iceland (Arnason, 1977) and for south 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-December 2014 (MacDonald et. al., 2016).
The high topographic gradient of the catchment, with large climatic differences between the upland glacial accumulation area 30 (>1800 m asl) and the lowland temperate proglacial area (<150 m asl), results in two easily distinguished isotopic compositions: (1) glacier meltwater; and (2) precipitation across the proglacial area. A binary mixing model for δ 2 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 δ 2 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 = (δ 2 HR -δ 2 HP) / (δ 2 HGW -δ 2 HR), where δ 2 HR is the composition of river water; δ 2 HP is the composition of local precipitation; and δ 2 HGW is the composition of sampled groundwater. Since most river recharge to the aquifer occurs during the summer months when river flow is high and 5 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.

Sandur structure and aquifer properties 10
The groundwater study catchment covers 6 km 2 and encompasses the sandur, adjacent hillslopes and moraines, and river outflow from the proglacial lake ( Figure 1). Geophysical evidence from the passive seismic and previous seismic reflection survey indicates that (Figure 2, Figure 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 ( Figure 2). All the sediment is of volcanic origin 15 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, 2002). Although not directly observed in the seismic data there is a possibility that at greater depth (>50 m) there exists more consolidated Pleistocene aged sediments which have been compacted by ice loading during earlier glaciations (Guðmundsson, 2002). Observations of bedrock from nearby exposures and two boreholes drilled 20 in bedrock reveal relatively massive and poorly fractured volcanic rock.  (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 25 of compaction due to burial and ice loading (Shmidt and McDonald, 1979), median hydraulic conductivity may have reduced to 15 m d -1 or at a worst case 6 m d -1 (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 2 d -1 ). 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 30 ground surface via perennial and ephemeral springs (Figure 2). A conservative estimate of the volume of groundwater stored in the full thickness of the aquifer is 51 ±15 million m 3 , approximately 1 -2 % of estimated ice volume in the glacier (Mackay et al. 2018).

Groundwater dynamics
Groundwater level elevation falls from upper to lower sandur, with a gradient of 0.018 across the upper and 0.013 across the lower sandur (Figure 2). In the upper sandur, closest to the glacier, groundwater levels adjacent to the glacial meltwater channel 5 are on average 1 m below river stage for most or all of the year (Figure 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 adjacent river stage 10 for much of the year, creating a piezometric gradient that drives visible groundwater discharge through seeps and springs to the river (Figure 2d Piezometers U1 and U2 (20 m 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 (Figure 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 (Figure 3a). During individual winter rain storms, groundwater levels in U2 rise higher than in U1, creating a 30 piezometric gradient towards the river (Figure 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 3 s -1 (IQ range 0.093-0.30 m 3 s -1 ), 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 3 s -1 (0.12-1.1 m 3 s -1 ) equivalent to 9.8% (3.6-22%) of mean annual river flow and 21% 5 (7.7-46%) of mean winter river flow.

Stable isotopes and temperature
Stable isotope composition (δ 2 H and δ 18 O) in groundwater from piezometers and springs was compared to that of glacier meltwater and local rainfall (Figure 4, Table S4). Previous studies have demonstrated that glacial meltwater and local rainfall on the proglacial area are easily distinguished using δ 2 H and δ 18 O due to the high elevation of the accumulation area 10 (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 -76.1 ± 2.6 δ 2 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 (Figure 4). Piezometers (U1, L1, 15 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. A binary mixing model developed for δ 2 H indicates the relative proportion of precipitation and glacier meltwater in 20 groundwater (Figure 4b, 4c) 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 (Figure 4c). Since the binary 25 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 (HCO3) are significantly lower in those piezometers strongly influenced by the river than those where precipitation influence is dominant ( Figure S1). River water temperature is relatively constant year-30 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).

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 5 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 & 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. 2018).

Groundwater flow
Our study shows that significant water can flow through a glacierized 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 15 is significant, 0.19 m 3 s -1 , equivalent to approximately 4.5% of mean annual river flow or 9.7% of mean winter river flow, with estimates 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 -1 ) 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 20 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). 25

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 30 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 ( Figures  5   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 10 soil development and vegetation cover. Recharge is likely to occur not only from direct precipitation on the sandur surface, but 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 Figure 4).
Looking forward, as the glacier continues to melt, the proglacial aquifer will continue to have a buffering effect on river flow. 15 High river flows will recharge the aquifer, whether caused by glacier icemelt, 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 20 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 deglaciating catchments.

Proglacial aquifers as strategic water resources 25
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 30 compared to surface water. If projected glacier losses and increased precipitation in glacierized catchments are realized (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.

Conclusions
Three years of investigations of groundwater and glacier meltwater at Virkisjökull, SE 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: 5 1. Direct measurements of aquifer characteristics show consistently high permeability (35 m d -1 n=64, IQR 25 -40 m d -1 ), and volume of groundwater storage (50 ±15 million m 3 ). The proglacial floodplain therefore forms a highly productive aquifer.
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-500 m from the river where glacier meltwater comprises > 25% and often >50% of recharge. 15 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.

20
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 characterizing 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. 25

Author contribution
BEOD managed the field sampling campaign and installation of piezometers and wrote early drafts of the manuscript. AMM oversaw the research and analysis and wrote the final draft of the paper. PW, MR, BEOD, JE, AB, AM and LJ undertook fieldwork and analysis of individual components of the research and WGD provide interpretation of the stable isotopes. All commented on final draft of the manuscript. 5

Competing interests
The authors declare that they have no conflict of interest.

Acknowledgments
This research was funded by the BGS-NERC Earth Hazards and Observatories Directorate and is published with the permission of the Executive Director of the British Geological Survey (NERC). We thank VatnajökulsÞjóðgarður for permission to install 10 monitoring equipment; Vatnsborun ehf for borehole drilling and installations; Icelandair for assistance with equipment transport; and Veðurstofa ĺslands and the people of Svinafell for research support. Tim Heaton at NERC Stable Isotope Facility undertook stable isotope analysis; James Sorensen carried out cluster analysis; and Craig Woodward at BGS assisted with diagrams (all BGS).