Proglacial groundwater storage dynamics under climate change and glacier retreat

Proglacial aquifers are an important water store in glacierised mountain catchments that supplement meltwater‐fed river flows and support freshwater ecosystems. Climate change and glacier retreat will perturb water storage in these aquifers, yet the climate‐glacier‐groundwater response cascade has rarely been studied and remains poorly understood. This study implements an integrated modelling approach that combines distributed glacio‐hydrological and groundwater models with climate change projections to evaluate the evolution of groundwater storage dynamics and surface‐groundwater exchanges in a temperate, glacierised catchment in Iceland. Focused infiltration along the meltwater‐fed Virkisá River channel is found to be an important source of groundwater recharge and is projected to provide 14%–20% of total groundwater recharge by the 2080s. The simulations highlight a mechanism by which glacier retreat could inhibit river recharge in the future due to the loss of diurnal melt cycling in the runoff hydrograph. However, the evolution of proglacial groundwater level dynamics show considerable resilience to changes in river recharge and, instead, are driven by changes in the magnitude and seasonal timing of diffuse recharge from year‐round rainfall. The majority of scenarios simulate an overall reduction in groundwater levels with a maximum 30‐day average groundwater level reduction of 1 m. The simulations replicate observational studies of baseflow to the river, where up to 15% of the 30‐day average river flow comes from groundwater outside of the melt season. This is forecast to reduce to 3%–8% by the 2080s due to increased contributions from rainfall and meltwater runoff. During the melt season, groundwater will continue to contribute 1%–3% of river flow despite significant reductions in meltwater runoff inputs. Therefore it is concluded that, in the proglacial region, groundwater will continue to provide only limited buffering of river flows as the glacier retreats.


| INTRODUCTION
Groundwater is increasingly being recognised as an important component of water cycling in the foreland areas of glacierised mountain catchments and could become strategically more important as water supply from rainfall and meltwater become less reliable under 21st century climate change (Ó Dochartaigh et al., 2019;Taylor, 2013;Vincent, Violette, & Aðalgeirsdóttir, 2019). Proglacial aquifers offer an accessible source of fresh water to downstream communities year-round (e.g., for drinking and irrigation, Stefania et al., 2018) and provide a steady supply of baseflow to glacier-fed rivers, helping to supplement low flows outside of the melt season (Jódar et al., 2017;MacDonald et al., 2016;Wilson, Williams, Kayastha, & Racoviteanu, 2016). This baseflow input can be an important regulator of glacier-fed river physiochemistry, providing more favourable habitat conditions for aquatic fauna (Crossman, Bradley, David, & Milner, 2012;Hotaling, Finn, Joseph Giersch, Weisrock, & Jacobsen, 2017;Khamis, Brown, Hannah, & Milner, 2016).
Recent studies indicate that climate change could cause longterm decline of proglacial aquifer storage due to reductions in diffuse recharge, from rainfall and snow melt (Somers, Mckenzie, Mark, & Lagos, 2019), and in focused recharge from glacier-fed river channels (Liljedahl, Gädeke, O'Neel, Gatesman, & Douglas, 2017). The dynamic behaviour of proglacial aquifer storage dynamics (e.g., intraannual variability) will also change in response to regional precipitation variability and meltwater runoff dynamics (Allen, Whitfield, & Werner, 2010). Meltwater runoff dynamics will be most important for aquifers hydraulically connected to meltwater river channels, where groundwater dynamics are closely coupled to river stage variability. These systems may exhibit groundwater storage behaviour that is more strongly tied to climatic conditions in the upstream glacierised region away from the aquifer. The trajectory of proglacial groundwater storage dynamics under climate change is, therefore, likely to be complex: driven in part by perturbations in diffuse recharge, and in part by perturbations in meltwater runoff dynamics as a consequence of glacier retreat. One must consider linkages along this climate-glacier-groundwater response cascade in order to understand proglacial groundwater vulnerability to climate change.
Yet, this cascade remains poorly understood and is likely to be highly non-linear and time-dependant (La Frenierre & Mark, 2014;Vincent et al., 2019).
Numerical models of glacio-hydrology and groundwater flow provide a means to investigate the climate-glacier-groundwater response cascade. To date, only one study has integrated the required process models to investigate this (Somers et al., 2019). They forced a distributed GSFLOW surface-groundwater model of the Shullcas Watershed in Peru with runoff simulations from a glacier melt model. They used this integrated modelling approach to project 21st century changes in the flow regime of the Shullcas River and showed that baseflow from groundwater could buffer reductions in dry-season runoff due to glacier retreat.
Other studies have integrated models to investigate groundwater storage dynamics in mountain catchments without glaciers. For example, Okkonen and Kløve (2011) linked a hydrological, land surface and distributed MODFLOW groundwater model to study changes in groundwater storage seasonality in an esker aquifer in Finland due to climate change. Their results showed an earlier rise and fall of groundwater storage due to a projected shift in seasonal peak diffuse recharge from spring to winter. Sridhar, Billah, and Hildreth (2018) integrated the Variable Infiltration Capacity hydrological model with MODFLOW to simulate climate change impacts on water table dynamics in an unconfined fractured basalt aquifer in the north-west of the United States. They also showed that future trajectories in water table elevation were strongly controlled by changes in diffuse recharge inputs. Huntington and Niswonger (2012) used GSFLOW to investigate the mechanism for observed reductions in summer streamflow in three mountain watersheds in California and Nevada. In contrast to other studies, they found that the trajectory of groundwater level dynamics were driven by changes in runoff patterns from neighbouring mountains rather than diffuse recharge (note, Scibek, Allen, Cannon, & Whitfield, 2007, also found this). More specifically, they found that climate warming induced an earlier recession in mountain spring snow melt runoff, inducing earlier drainage of the alluvial aquifer to the stream, and subsequently, a fall in baseflow to the stream in the summer.
The findings of Huntington and Niswonger (2012) raises the question of how groundwater storage dynamics in glacierised catchments might respond to changes in meltwater runoff due to glacier retreat. This question is particularly pertinent, given that glaciers around the world are retreating rapidly (Zemp et al., 2019). A number of field observation studies have demonstrated the sensitivity of proglacial groundwater storage dynamics to glacier meltwater runoff variability, especially where aquifers are in strong connection with meltwater river channels (Dragon, Marciniak, Szpikowski, Szpikowska, & Wawrzyniak, 2015;Levy, Robinson, Krause, Waller, & Weatherill, 2015;Ó Dochartaigh et al., 2019;Robinson, Fairchild, & Arrowsmith, 2009). Yet, our understanding of how proglacial groundwater storage dynamics will respond to 21st century climate change and prolonged glacier retreat remains limited. This is in part, due to the scarcity of high quality observation data of local climate, glacier mass balance, meltwater flows and aquifer storage and characterisation (Vincent et al., 2019).
This study uses an integrated modelling approach consisting of a glacio-hydrological model and a distributed Newton-Raphson formulation of MODFLOW (MODFLOW-NWT) to simulate groundwater storage dynamics and surface-groundwater interactions in the foreland region of the well-characterised and monitored Virkisá River basin in Iceland. Through driving the models with climate change projections, they are used to investigate: (1) twenty-first century changes in groundwater storage dynamics in the proglacial aquifer in response to climate change and glacier retreat; and (2) the principal drivers of these changes including diffuse recharge and meltwater runoff dynamics.

| Study site
The Virkisá River basin covers an area of 22 km 2 on the western side of the ice-capped Öraefajökull stratovolcano in south-east Iceland (Figure 1(a)). Approximately 60% of basin is covered by Virkisjökull and Falljökull, which together comprise a twin-lobed outlet glacier of the Öraefajökul ice cap (herein referred to as Virkisjökull). Between 1988 and 2011 Virkisjökull lost 0.3 km 3 of ice and retreated 0.5 km. Continued retreat over the last decade has formed a small proglacial lake at the terminus which forms the headwater of the Virkisá River, the only drainage pathway for melt and rainfall runoff from the basin. The river flows through an 800 m bedrock-controlled section flanked on either side by push moraines and then over a sandur floodplain which forms a shallow, unconfined and gently sloping (average surface gradient = 0.017) aquifer. This is made up of loosely consolidated, moderately to poorly sorted, dominantly medium-to F I G U R E 1 Study site including topographical map of Virkisá river basin and proglacial sandur groundwater catchment with instrumentation (a); conceptual model of water flow along cross section from glacier to groundwater adapted from Ó Dochartaigh et al. (2019) (b); and photograph of sandur floodplain, meltwater channels and Virkisjökull taken from the southern boundary of the groundwater catchment (c) coarse-grained glaciofluvial sand, gravel and cobbles which have been deposited from actively shifting meltwater streams and frequent (c. 5 per century) jökulhlaups (Robinson, Fairchild, & Russell, 2008). Geophysical evidence from Tromino® passive seismic surveys indicate the deposits are up to 150 m thick within several kilometres of the lake outlet (Ó Dochartaigh et al., 2019). The Virkisá groundwater catchment forms part of the world's largest sandur, Skeiðarársandur.
The sandur aquifer is recharged by diffuse rainfall and snow melt infiltration as well as focused infiltration of mountain runoff along the Virkisá River channel (Figure 1(a)). Shallow borehole investigations have shown that the groundwater flows downstream, approximately parallel to the river (Ó Dochartaigh, MacDonald, Wilson, & Bonsor, 2012). However, stable isotope data have shown that there is a zone of meltwater-groundwater mixing within 500 m of river channel where more than 25% of groundwater is derived from mountain meltwater runoff (Ó Dochartaigh et al., 2019). River-aquifer water exchanges in both directions are promoted by the shallow water table which is typically between 0 and 4.4 m below ground level. Transmissivity has been estimated to be between 100-2500 m 2 d −1 with a median value of 600 m 2 d −1 from pumping tests conducted at eight sandur piezometers (green markers in Figure 1 The local maritime climate is characterised by cool summers (10 C on average near the glacier terminus) and mild winters (1 C on average) with an average temperature lapse rate of −0.44 C 100 m −1 (Flett, 2016). There is a significant lateral precipitation gradient due to the prevailing north-easterly winds and orographic effects.

| Integrated modelling framework
The modelling framework consists of the three principal components to represent the climate-glacier-groundwater response cascade including historical and future driving climate data, a distributed conceptual glacio-hydrological model and a distributed physically-based groundwater model. The climate data provide boundary conditions for the glacio-hydrological model which simulates runoff from the glacierised region to the Virkisá River and diffuse recharge from rainfall and snow melt over the sandur. These simulations provide boundary conditions for the groundwater model of the sandur aquifer. The climate data and glacio-hydrological model are documented in detail by Mackay et al. (2018Mackay et al. ( , 2019 and accordingly, are only briefly outlined here. 2.2.1 | Historical and future climate data Mackay et al. (2018) compiled what are considered to be the most reliable record of hourly near-surface air temperature, precipitation and incident solar radiation (ISR) data for the Virkisá River basin. All data span 1981-2016 inclusive. A continuous temperature record for the weather station at the glacier terminus (Figure 1(a)) was constructed by combining in-situ weather station measurements with records from the Icelandic Meteorological Office weather station network. A combination of in-situ ISR measurements and a statistical model for converting temperature to ISR was used to generate a continuous ISR time series for the catchment. For precipitation, state-ofthe-art 2.5 km gridded reanalysis precipitation data were used (Nawri et al., 2017) after bias-correction against in-situ rainfall records using equidistant quantile mapping (Li, Sheffield, & Wood, 2010).

| GHM++ glacio-hydrological model
The distributed GHM++ glacio-hydrological model code developed by Mackay (2020)  3 years of seasonal ablation stake data; and remotely-sensed ice and snow cover data. GHM++ has not been used to provide river discharge or diffuse recharge boundary conditions to a groundwater model before, but was deemed suitable for this application as it has shown to be efficient at simulating the Virkisá River flow dynamics and it uses a soil water balance routine to calculate diffuse recharge which has been applied extensively in temperate environments (Jackson, Wang, Pachocka, Mackay, & Bloomfield, 2016;Mackay, Jackson, & Wang, 2014;Mansour, Wang, Whiteman, & Hughes, 2018) and compares favourably to physicallybased models at the field scale where interception losses are small (Sorensen, Finch, Ireson, & Jackson, 2014).

| MODLFOW-NWT groundwater model
The distributed, physically-based MODFLOW-NWT groundwater modelling code (Niswonger, Panday, & Ibaraki, 2011) was used to simulate groundwater flow through the sandur and water exchanges with the Virkisá River and springs. MODLFOW-NWT uses a Newton-Raphson solver which has been designed specifically for simulating non-linear unconfined groundwater flow problems. The model domain was divided into a Cartesian grid with 50 m horizontal and 10 m vertical resolution. The model time step was set to hourly. These discretisations were selected as a compromise between model complexity and runtime requirements, to ensure consistency with GHM++ and to ensure that the aquifer geometry could be adequately represented. The upstream model boundary (B1 in Figure 2(a)) coincides with known outcropping of the bedrock 800 m downstream of the lake outlet. Pumping tests conducted by Ó Dochartaigh et al. (2012) have shown that the permeability of the volcanic bedrock material is negligible, and so no flow model nodes were specified at this boundary. No flow conditions were also imposed along the boundaries to the east and west of the river (B2 and B3 in Figure 2 (a)), given that groundwater predominantly flows parallel to the river.
At the downstream boundary (B4 in Figure 2(a)), specified flux boundary conditions were prescribed by assuming that the hydraulic gradient is equal to the topographic slope, which is reasonable given that the water table resides close to the ground surface year-round.
The thickness of the sandur model was set according to the inter- were also imposed at the aquifer-bedrock interface.
The Virkisá River was represented by a single rectangular channel running between the upstream and downstream model boundaries.
F I G U R E 2 Groundwater model lateral (a) and vertical (b) extent and model used to define river width variations (c). Aerial image taken on 7 August, 2012,Source:DigitalGlobe (Vivid -Iceland), created using ArcGIS. Copyright © Esri. All rights reserved This is contrary to observed channel braiding in the lower sandur, but was deemed a necessary simplification, given that sediment reworking causes regular (sub-annual) changes to channel morphology and the deactivation and reactivation of channels (Marren, 2005) which would be practically infeasible to simulate. The braided channel system promotes meltwater exchanges with the aquifer. To mimic this, the width of the river channel was varied in the groundwater model according to the degree of braiding downstream. Here, the total width of all river braids were measured at 50 m intervals along the main river channel (blue dots in Figure 2(c)). A polynomial was then fit to these data points (black dash line in Figure 2(c)) to parameterise the river channel width in the groundwater model. The streamflow-routing package (SFR1) (Prudic, Konikow, & Banta, 2004) was used to route flow along the river channel. It calculates surface-groundwater exchanges across the river bed using Darcy's Law.
No direct measurements of the river bed topography exist. However, numerous sections of the river have been waded during field excursions where water depth is typically between 0.2 and 1 m with an estimated average water depth of 0.5 m. Accordingly, the river bed bottom elevation was set for each river section by subtracting 0.5 m from the 5 m airborne lidar digital elevation model in Figure 1(a).
In addition to the river, an extensive network of ephemeral and perennial springs exist, particularly at the downstream end of the sandur. These springs are thought to form a significant drainage output and so were included in the model using the MODFLOW drain package. Here, drains were placed at the land surface elevation at each x and y coordinate covering the lateral extent of the model domain.

| Groundwater model calibration
Eight unknown model parameters relating to aquifer and river hydraulic properties needed to be identified (Table 1). Depositional processes on alluvial aquifers bring about heterogeneity in aquifer permeability and porosity (Chen, Song, & Wang, 2010;Neton, Dorsch, Olson, & Young, 1994). Pumping tests at the sandur piezometers indicate that transmissivity is high (median 600 m 2 d −1 ), but do not provide evidence of spatial variability in aquifer permeability or storage (Ó Dochartaigh et al., 2012). Surface permeability measurements indicate no significant spatial variability (Ó Dochartaigh et al., 2019). Accordingly, the hydraulic properties were uniformly parameterised over the sandur.
A number of preliminary steady-state simulations were used to explore parameter sensitivity and refine their behavioural ranges.
Groundwater level simulations showed to be insensitive to the prescribed river bed thickness parameter. This was set to an estimated thickness of 0.5 m.
Given that the sandur is unconfined, the specific storage was also fixed to a representative value for unconsolidated coarse sand aquifers (Domenico & Schwartz, 1990). The river bed hydraulic conductivity and drain conductance were both set using the mean of the surface permeability measurements undertaken by Ó Dochartaigh et al. (2019). Finally, the Manning's roughness coefficient was set to 0.05 based on the "normal" value for mountain streams with beds made up of cobbles and large boulders (Chow, 1959).
The horizontal hydraulic conductivity, vertical anisotropy and specific yield were calibrated against hourly groundwater level time series from the sandur piezometers which span August 2012-December 2017 inclusive. Groundwater level fluctuations at one of the piezometers are known to be affected by a small ephemeral surface water channel which causes discrete focused groundwater recharge in the immediate vicinity of the piezometer. Accordingly, only data from seven piezometers were used for model calibration (see Figure 2(a)).
The parameters were calibrated from a 5000-run Monte Carlo procedure with parameter sets drawn from uniform distributions using Sobol sampling (Brately & Fox, 1988 hereafter referred to as scenarios, were selected to represent a range of potential future outcomes. While this approach precludes a full examination of model projection uncertainty, using a handful of scenarios crucially allows a more process-oriented analysis of the individual scenarios to be undertaken to better-capture the nuanced feedbacks between climate, glacio-hydrology and groundwater that might be lost when averaging simulations over large ensembles (Knutti, Furrer, Tebaldi, Cermak, & Meehl, 2010). This is particularly beneficial for evaluating water cycling in unconfined aquifers given their non-linear response to boundary forcing (Cayar & Kavvas, 2009).
Given the emphasis on a small number of scenarios, it was important to adopt a robust scenario selection procedure that captured a range of future glacio-hydrological and climatic behaviours important to proglacial groundwater dynamics. Here, the kmeans clustering (MacQueen, 1967) approach was used, which finds natural groupings in multivariate data that are distinct from one another and has been used in the past for climate projection scenario selection (Cannon, 2015;Wilcke & Bärring, 2016). For this study, clustering was based on 11 quantitative signatures of climatic and glacio-hydrological behaviour which were calculated for each of the 94,080 possible scenarios ( Table 2). The MATLAB® kmeans function was used which uses an iterative approach to determine the optimal location of k centroids. Given the potential to reach local minima in the solution space during the optimisation procedure, 100 replicates were undertaken using random centroid locations initialised using the kmeans++ algorithm.
It was decided to focus the projections on the time slice centred on the 2080s (2073-2097) given that Mackay et al. (2019) found this period to have the largest spread in projected evolution of climate and glacio-hydrological characteristics.
3 | RESULTS proglacial fans and sandur aquifers (Parriaux & Nicoud, 1990). The model captures the seasonality of peizometer observations as well as the timing of individual, event-scale peaks (Figure 3(a-g)). It also shows some deficiencies including a systematic underestimation of summer levels at U2 and M1 and overestimation at M2 and L3. At L1 and L2 the amplitude of seasonality is overestimated. Even so, the model captures the spatial distribution of mean depth to water table accurately, obtaining an R 2 of 0.84 (Figure 3(h)). It is also relatively consistent across the calibration and evaluation periods where the average RMSE only rises from 0.39 to 0.40 (Table 3). Figure 4 shows where the signatures for the five selected scenarios (coloured from orange to blue to indicate the degree of glacier retreat)

| Twenty-first century scenario selection
F I G U R E 3 Observed (blue) and simulated (yellow) groundwater level time series at the seven piezometers (a-g) and mean observed and simulated depth to water table (h). Note, dashed black line in time series indicates divide between calibration (left) and evaluation (right) periods lie on the distribution of all possible future scenarios. The selected scenarios cover between 38% (spring river discharge) and 76% (summer river discharge) of the ensemble. The extremes are not always represented in the selected scenarios, but on average, they cover 61% of the full range of possible future scenarios.
Each scenario has been assigned a code with the format G*-Q*-P* where G, Q and P represent glacier coverage, mean annual river discharge and mean annual sandur precipitation respectively and * is either 1, 2 or 3 representing low, moderate and high respectively.
These bandings were defined based on which tercile of the full distribution they fall into. Table 4 provides a summary of each scenario.
Average near-surface air temperature is projected to rise by 0.7-3.1 C between the historical reference period  and the 2080s. Glacier coverage (ΔG) is simulated to reduce by 3.7-8.1 km 2 . Glacier retreat is greatest for the warmest scenarios. All scenarios show a reduction in mean river discharge (ΔQ) for the 2080s. The largest reductions are observed when glacier retreat is greatest with the exception of the G3-Q1-P1 scenario in which glacier retreat is relatively small, but the reduction in river flow is large. This is because it is the driest scenario in terms of sandur precipitation. All scenarios except for G2-Q2-P3 show a reduction in total precipitation over the sandur relative to the reference period (ΔP).

| Groundwater storage dynamics
Daily average groundwater level and diffuse recharge simulations are presented in Figure 5. Note, these are shown as an average over the    Note: Each have been assigned a code with the format G*-Q*-P* where G, Q and P represent glacier coverage, mean annual river discharge and mean annual sandur precipitation respectively and * is either 1, 2 or 3 representing low, moderate or high respectively. Also shown are the changes in nearsurface air temperature (ΔT), glacier coverage (ΔG), upstream river discharge (ΔQ) and total precipitation over the sandur (ΔP).

| Baseflow dynamics
F I G U R E 5 Intra-annual distribution plots of simulated groundwater level (a) and diffuse recharge (b) over the upper sandur, averaged over all scenarios for the reference period (1991-2015, dashed black line) and shown for individual scenarios for the 2080s (coloured lines). Also shown are the changes in groundwater level (c) and diffuse recharge (d) between the reference and 2080s periods between September and March and the smallest baseflow fluxes occur during the warmer months between May and August. The reference period simulations indicate that groundwater contributes up to 15% of the total sandur runoff in the winter months. This is within the range estimated from stable isotope data (MacDonald et al., 2016).
Note, the exact contributions in the winter months depend on the reference scenario as shown by the reference run range (grey band in Figure 5(b)).
All of the scenarios show a small reduction in baseflow in the winter months (Figure 6(d)), which partly explains the simulated reduction in groundwater contribution to sandur runoff in these months ( Figure 6(e)). For the three warmest scenarios (orange, yellow and green), this is compounded by the simulated rise in runoff from the Virkisá River basin during the winter months as more precipitation falls as rainfall (Figure 6(c, f)).
During the melt season (May to September), the direction of change in baseflow fluxes by the 2080s depends on the scenario, but all scenarios consistently simulate 1% increase in proportional contribution of groundwater to sandur runoff due to the simulated reduction in runoff from the Virkisá River basin (Figure 6(f)). During the melt season, groundwater drainage only makes up 1%-2% of sandur runoff in the reference and future simulations.

| River recharge dynamics
The seasonality of river recharge is consistent across the reference and future scenarios where river recharge fluxes are highest in the summer months and lowest in the winter months (Figure 7(a)). Over the reference period, the 30-day average river recharge contributes up to 39% of total recharge to the groundwater catchment in the melt season. The reference simulations indicate that overall, river recharge contributes 13%-17% of total recharge to the groundwater catchment. Note, the range reported here is due to differences in simulated river recharge between the reference runs (grey band in Figure 7(a)), particularly the G2-Q2-P3 scenario (indicated with arrow in Figure 7 (a)) which showed a much higher peak river recharge during the summer months (up to 54 m 3 h −1 higher than the other scenarios).
The seasonality of simulated changes in river recharge (Figure 7 (c)) closely follow the relative simulated changes in length of losing river (Figure 7(b)) which is largely controlled by changes in groundwater level. All simulations for the 2080s project increases in river recharge between November and March where the length of losing river is projected to increase. For the spring and summer months, the wettest scenarios (yellow and blue: higher diffuse recharge, higher groundwater level and lower length of losing river) simulate a reduction in spring and summer recharge while the driest scenarios (orange, green and cyan: lower diffuse recharge, lower groundwater level and higher length of losing river) show an increase in river recharge for the majority of the year. The relative wetness (P in scenario code) is more influential on the simulated change in river recharge, than the magnitude of river flow (Q in scenario code). The projections indicate that the contribution of river recharge to total recharge in the groundwater catchment will increase to between 14% and 20% by the 2080s.
However, this is not true for the G2-Q2-P3 scenario (yellow) which shows a reduction in July river recharge of >50 m 3 h −1 which is F I G U R E 6 Intra-annual distribution plots of simulated baseflow (a), proportion of sandur runoff from groundwater (b), and runoff from the Virkisá River basin (c), averaged over all scenarios for the reference period (1991-2015, dashed black line with range shown by grey band) and shown for individual scenarios for the 2080s (coloured lines). Also shown are the corresponding changes between the reference and 2080s periods (d-f) not reflected in the simulated change in losing river length. Instead, this reduction is associated with a large reduction in the specific river recharge that is, the rate of recharge per unit area of losing river bed (Figure 7(d)).
To investigate this, Figure 8 in a larger downward head gradient (Figure 8(b)) and an enhanced river channel infiltration flux (Figure 8(c)). It is the loss of this diurnal melt signal for the G2-Q2-P3 scenario that explains the large reduction in projected river recharge during the melt season for the 2080s.

| DISCUSSION
From the results of this study, four key findings are identified and discussed below in the context of groundwater resources in glacierised basins, methodological limitations and future research needs.
F I G U R E 7 Intra-annual distribution plots of simulated river recharge (a). Note, the dashed yellow line indicated with an arrow shows the reference simulation for the G2-Q2-P3 scenario. Also shown are the changes in length of loosing river (b), river recharge (c) and specific river recharge (d) between the reference and 2080s periods F I G U R E 8 Simulated river stage (a), river stage minus groundwater head (b) and river channel infiltration (c) time series for a single MODFLOW-NWT river node 1 km downstream of the upstream model boundary for 3 days during July 2013 for the reference runs of all scenarios 4.1 | The Virkisá River is a significant source of proglacial groundwater recharge The contribution of meltwater channel infiltration to groundwater recharge in proglacial aquifers has rarely been quantified (see Liljedahl et al., 2017;Ó Dochartaigh et al., 2019, for two notable exceptions) and, to our knowledge, only the study of Somers et al. (2019) has calculated the relative significance of this (2% of total recharge). The integrated modelling approach used in this study estimates that between 1991 and 2015, the Virkisá River contributed 13%-17% of total recharge to the groundwater catchment. The river is, therefore, a significant source of groundwater recharge in the Virkisá groundwater catchment. Furthermore, the relative contribution of river channel infiltration to total groundwater recharge is projected to rise to 14%-

| Glacier retreat could inhibit river recharge
All five scenarios evaluated in this study simulated reductions in mean runoff to the Virkisá River in the 2080s due to glacier retreat. However, the corresponding groundwater model simulations did not provide any clear evidence that changes in mean runoff inhibits future river recharge. The two principal mechanisms by which a reduction in mean river flow could inhibit river recharge is through lowering the downward head gradient at the river-aquifer interface and through reducing the area of active river bed in hydraulic connection with the aquifer. The simplified river geometry adopted in this study (single rectangular channel) means that only the former could be investigated. However, the geometry of the Virkisá River is highly dynamic with frequent activation and deactivation of channels which could lead to significant changes in active river bed area. It is, therefore, important to highlight that the apparent insensitivity of river recharge to projected reductions in river flow could be partly explained by the limitations of the river model.
Even without the ability to fully evaluate the impact of reduced river runoff on river recharge, this study did show that the flattening of diurnal flow variations in the runoff hydrograph due to the loss of daily glacier melt cycling could inhibit river recharge. The diurnal melt cycle can induce a downward head gradient at the river-aquifer interface that drives river recharge inputs to the aquifer. For one of the scenarios, the loss of this cycling reduced summer river recharge fluxes by up to 29%. This is the first study to demonstrate this mechanism of glacier retreat impact on groundwater recharge. It is hypothesised here, that this mechanism is relevant to the vulnerability of other proglacial aquifers and therefore requires further investigation.
Crucial to these investigations will be implementing monitoring and numerical modelling approaches that consider sub-daily (ideally hourly) fluctuations in key hydrological variables.

| Groundwater storage dynamics are resilient to changes in river recharge
The simulations showed that groundwater storage dynamics in the Virkisá groundwater catchment are driven by changes in diffuse recharge and are largely resilient to changes in river recharge. The drier scenarios with reduced sandur precipitation consistently simulated lower groundwater levels in the 2080s while shifts in the timing of groundwater storage seasonality were controlled by changes in the magnitude of diffuse recharge during the main recharge season. This finding is perhaps to be expected, given that the sandur is situated in a temperate environment and receives precipitation year-round which makes up the majority of groundwater recharge. However, numerical modelling studies in semi-arid mountain aquifers have found groundwater storage dynamics to be most sensitive to river flow dynamics (Allen, Mackie, & Wei, 2004;Huntington & Niswonger, 2012;Scibek et al., 2007). The primary drivers of proglacial groundwater storage dynamics are, therefore, likely to be strongly influenced by key environmental factors such as the prevailing climate. The findings from this study are most relevant to glacierised catchments situated in temperate environments.
It is also important to highlight that the simulations presented in this study represent an areal average over the groundwater catchment and therefore do not capture some of the spatial variability that might be expected. Ó Dochartaigh et al. (2019) showed that the groundwater level records from the piezometers closest to the river (U1, M1 and L1 in Figure 2) show a strong signal resembling the river stage dynamics. Accordingly, regions closer to the river are expected to be more influenced by changes in river flow dynamics.

| Groundwater will provide limited buffering of proglacial river runoff under climate change
The reference period simulations estimated that groundwater has contributed up to 15% of 30-day average flow in the Virkisá River outside of melt season. This is at the lower end of the estimate based on stable isotope analysis of water samples (MacDonald et al., 2016), but still demonstrates the historical importance of groundwater in buffering river flow outside of the melt season. In the 2080s, groundwater is projected to contribute up to 3%-8% of 30-day average river flow outside of the melt season due to simulated lowering of the water table and increased meltwater and rainfall-runoff in the winter months. During the melt season, groundwater will continue to contribute 1%-3% of river flow despite significant reductions in meltwater runoff inputs. Groundwater will, therefore, continue to provide limited buffering of river flows throughout the year, although, will contribute proportionally less to river flows outside of the melt season. The proportional reduction in groundwater contribution to river flow in the proglacial region could have implications for the ecological health of the river (Khamis et al., 2016). However, it should be emphasised, that these results are only relevant for the proglacial region considered in this study. If, for example, the length of the groundwater catchment were extended to include additional downstream regions of the Virkisá River, one would expect the contribution and buffering from groundwater to increase.

| CONCLUSIONS
Through implementing an integrated numerical modelling framework to represent the climate-glacier-groundwater response cascade in a well-characterised catchment, this study has provided a rare analysis The findings indicate that groundwater will continue to provide limited buffering of river flows over the 21st century, but that this buffering is sensitive to changes in groundwater storage dynamics, which are mainly driven by future patterns of diffuse recharge. Our understanding of water cycling in proglacial aquifers remains limited and it is likely that groundwater storage dynamics for catchments with different climate regimes will respond to climate change and glacier retreat differently. It is imperative that future research of these systems is undertaken in other parts of the world, particularly given the potential socio-economic and environmental importance of proglacial aquifers. Integrated modelling approaches and high quality observation data should underpin this research.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. The convergence ratio only reaches 57% for κ ( Figure A1(c)) indicating that it is less identifiable than K h . Even so, the parameter still converges toward a value >0.5 and the calibrated κ is 0.98 indicating an aquifer lithology that is approximately isotropic. The calibrated S y is 0.15, but the convergence ratio only reaches 4% ( Figure A1(d)) indicating that this parameter is not identifiable. This is likely because inter-piezometer groundwater level variability is much larger than variations in groundwater levels at individual piezometers. Accordingly, the calibration of hydraulic model parameters to all of the groundwater level time series data simultaneously likely renders S y insensitive. To investigate this further, an additional experiment was undertaken whereby the model was calibrated to each piezometer individually. It was found that, except for the U1 and U2 piezometers, which are situated close to the river and therefore controlled primarily by river dynamics, S y showed to be more sensitive with the convergence ratio ranging from 27% to 99% with a mean of 73% and the calibrated S y ranging from 11% to 25% with a mean of 15%. Given that the S y obtained when taking the average of the calibration values at all sensitive piezometers is in agreement with the value obtained from the initial Monte Carlo calibration and that this is in the representative porosity range (10%-20%) for proglacial fans and sandur aquifers (Parriaux & Nicoud, 1990), an S y of 0.15 was deemed justifiable. Figure A1(e-h) show equivalent convergence plots for the mean simulated state variables as the RMSE threshold is reduced. All state variables give a convergence ratio ≥89% indicating that the model behaviour is well constrained by the calibration procedure.
F I G U R E A 1 Frequency histogram (blue bars) and cumulative frequency diagram (black line) of RMSE scores obtained from successful Monte Carlo calibration runs (a). Also shown is the range (light blue lines), mean (dashed blue line) and convergence ratio (yellow line) of the calibration parameters (b-d) and model state variables (e-h) for all model runs under each point of the cumulative frequency diagram