Groundwater Buffers Decreasing Glacier Melt in an Andean Watershed—But Not Forever

Accelerating mountain glacier recession in a warming climate threatens the sustainability of mountain water resources. The extent to which groundwater will provide resilience to these water resources is unknown, in part due to a lack of data and poorly understood interactions between groundwater and surface water. Here we address this knowledge gap by linking climate, glaciers, surface water, and groundwater into an integrated model of the Shullcas Watershed, Peru, in the tropical Andes, the region experiencing the most rapid mountain‐glacier retreat on Earth. For a range of climate scenarios, our model projects that glaciers will disappear by 2100. The loss of glacial meltwater will be buffered by relatively consistent groundwater discharge, which only receives minor recharge (~2%) from glacier melt. However, increasing temperature and associated evapotranspiration, alongside potential decreases in precipitation, will decrease groundwater recharge and streamflow, particularly for the RCP 8.5 emission scenario.


Introduction
Mountains sustain water resources by gradually and consistently supplying snow and ice melt for downstream communities. However, accelerating glacier recession in a warming climate threatens the sustainability of these natural water towers (Barnett et al., 2005;Viviroli et al., 2007).
Glacier retreat has accelerated worldwide since the 1980s (Kaser et al., 2006;Paul et al., 2004), with tropical alpine glaciers experiencing the most rapid areal recession (Bradley et al., 2006;Paul et al., 2004;Rabatel et al., 2013;Vuille et al., 2018). The Peruvian Andes, which host 70% of the world's tropical glaciers, have had coverage decrease by over 40% since the 1970s (Autoridad Nacional del Agua, 2014). Communities and industries in the Andes and on the arid Pacific coast are vulnerable to climate-driven intensification of dry season (May-September) water shortages , when only 20% of annual precipitation occurs in the Peruvian Andes (Bury et al., 2013).
Groundwater, previously often ignored in mountain hydrology (Bales et al., 2006;Huss et al., 2008;Juen et al., 2007;Ragettli et al., 2014;Voeckler et al., 2014), is now understood to be a significant contributor to streamflow generation in the Andes Glas et al., 2018;Gordon et al., 2015;Saberi et al., 2019;Somers et al., 2016) and in mountain systems globally (Andermann et al., 2012;Frisbee et al., 2011; ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2019GL084730
Key Points: • Groundwater accounts for a large fraction of streamflow and only receives minor (~2%) recharge from glaciers in the study catchment in Peru • As meltwater decreases, groundwater provides consistent discharge in the near term (~30 years), becoming a larger fraction of streamflow • In the long term (>60 years), groundwater storage and discharge decrease in response to higher evapotranspiration and lower precipitation Supporting Information: • Supporting Information Figure S1 • Figure S1 • Figure S2 • Figure S3 • Figure S4 Hood & Hayashi, 2015; Liu et al., 2004;McClymont et al., 2010). While groundwater supplies may prove more resilient than snow and glaciers to climate change (Lemieux et al., 2008;Tague et al., 2008), our knowledge of future impacts is limited by oversimplification of groundwater processes in hydrological models due computational and/or observational limitations in mountain regions (Juen et al., 2007;Ragettli et al., 2014). While there have been some implementations of integrated surface-subsurface hydrologic models in glacierized watersheds (Lemieux et al., 2008;Saberi et al., 2019), it is unclear how results from these few studies generalize to other settings or extrapolate into the future as glaciers retreat.
Here we directly address this research gap by developing a modeling approach that integrates distributed groundwater flow and recharge with surface-water and glacier-melt processes. We constrain our model with extensive field observations from 1998 to 2018. Using our integrated modeling approach, we first elucidate the role of mountain groundwater in the proglacial Andes by characterizing dominant streamflow contributors and flow paths, as well as quantifying glacier-groundwater interactions. Secondly, we apply climate projections to the integrated model to assess the capacity of groundwater to buffer the hydrological impacts of climate change and mountain glacier recession.

Study Area
The Shullcas Watershed (area = 170 km 2 ), central Peru, is a typical proglacial watershed (Figure 1), and is comprised of glacierized peaks, steep alpine grasslands, and flatter valley-bottom wetlands. It is underlain by glacio-fluvial sediments above metamorphic rock of the Huaytapallana Complex in the north and weakly metamorphosed sedimentary rock of the Mitu, Cabanilla, and Pucara Groups in the south (Chew et al., 2016). The steep alpine grasslands are dominated by Jarva ichu grass, while the flat alpine wetlands (known as bofedales) host cushion-forming vegetation such as Disticha muscoides and other plants including Plantaeog rigida, and Jarva ichu (Maldonado Fonkén, 2014;Salvador et al., 2014). The Shullcas River receives meltwater from the glaciers of the Cordillera Huaytapallana (Figure 2a), which have retreated by 44% between 1984(López-Moreno et al., 2014. The Shullcas River is the primary water source for the city of Huancayo (population: 365 K; Instituto Nacional de Estadística y Informática, 2012), which In the future, air temperature will continue to increase, while changes in precipitation are more uncertain. In response, glacier melt will decrease and eventually disappear, ET will increase, and groundwater storage will decrease (i.e., lower water table), all leading to lower stream discharge for all modeled climate scenarios.
experiences frequent water shortages during the dry season (Mark et al., 2017). The social and economic importance of this watershed and availability of climate, glacier, hydrologic, and hydrogeological data make the Shullcas Watershed an ideal test site to use a system approach to assess the impact of climate change on Andean hydrology.

Data Collection and Processing
We obtained daily maximum and minimum air temperature and precipitation, from the National Meteorology and Hydrology Service of Peru (SENAMHI) for two meteorological stations (Huaytapallana, 4,684 m above sea level, and Shullcas, 3,839 m above sea level; Figure 2a) within the watershed from 2008 to 2018 (SENAMHI, 2018). Data from two additional weather stations (Huayao and Santa Ana) located within 20 km of the watershed and from NASA's Tropical Rainfall Measurement Mission (TRMM version 3B42 RR; Huffman & Bolvin, 2015) were used to extend temperature and precipitation inputs, respectively, back to 1998 and fill data gaps using linear regression to the primary stations.
We obtained monthly streamflow data at the watershed outflow (Location 6; Figure 2a) from 1985 to 2009 from the Peruvian National Water Authority (ANA). In June 2015, we installed five stream-discharge gauges and two lake-level gauges. In July 2016, we installed an additional stream gauging station and seven shallow (<2.5 m) groundwater wells in two clusters ( Figure 2a).

Integrated Hydrological Model Setup
We integrated GSFLOW v2.0.0 (Markstrom et al., 2008), a coupled groundwater-surface-water model, with a glacier melt module to simulate the complete hydrological system of the Shullcas Watershed. GSFLOW combines MODFLOW (Langevin et al., 2017), a finite-difference groundwater flow model, and PRMS (Markstrom et al., 2015), a distributed-parameter, process-based, surface water model. GSFLOW has a daily time step, and we prepared our GSFLOW input files using GRASS-GSFLOW version 1.0.0 (Ng et al., 2018). Input parameters were based on field measurements where possible, and literature values otherwise (see supporting information and Table S1). For our purposes, we define the shallow unsaturated zone as the top~0.5 m of the subsurface, which is controlled by PRMS within GSFLOW. The MODFLOW domain above the water table is defined as the deep unsaturated zone and the MODFLOW domain below the water table is the saturated zone. We discretized the surface water component of the model into 45 hydrologic response units (HRUs) classified as either glacier/rock, grassland, or grassland/wetland based on a supervised maximum likelihood classification of Landsat satellite imagery (after Mather, 1985). We used the ASTER Global Digital Elevation Model with 30-m resolution (Tachikawa et al., 2011) to define the watershed boundaries and topography. Precipitation and temperature were spatially distributed based on elevation lapse rates calculated from the two meteorological stations within the watershed. Given the lack of humidity and wind speed data, potential evapotranspiration (PET) was calculated using a Priestley-Taylor approach (Liu et al., 2017;Priestley & Taylor, 1972). Actual evapotranspiration (ET) was calculated in GSFLOW according to how much water was available in the capillary reservoir of the soil zone. After ET, excess water in the soil zone could go to interflow, based on fast and slow interflow coefficients, and to the deeper subsurface MODFLOW domain (See Markstrom et al., 2008).
The subsurface was discretized horizontally into 2,850, 242 by 246-m grid cells. We impose an upper layer thickness of 20 m based on road cut observations that exposed over 10 m of poorly sorted soils. Beneath this, the fractured bedrock layer extends an additional 180 m. Below 200-m depth, we assume that bedrock is impermeable. We placed three constant-head MODFLOW cells downstream of the watershed outlet point to allow for some groundwater outflow from the modeled domain and prevent forced groundwater discharge to the stream network.
We used a semidistributed, modified temperature-index model to simulate glacier melt and mass balance for the Huaytapallana glaciers (after Quick & Pipes, 1977; see supporting information S3 for more details). Temperature and precipitation were calculated for each 100-m elevation band per time step. Melt, m, for each elevation band, b, was calculated as where MF s is the melt factor of a given season s (dry or rainy), T is the average daily temperature (°C), and T crit is the minimum temperature where melt is possible (−1°C in this case). We calibrated seasonal melt factors to match glacial lake outflow for Chuspicocha Lake from 2015 to 2018. At each time step, t, the mass balance of the glacier, Mass, is calculated as where is the density of water (1,000 kg/m 3 ), Sn is the total daily snowfall (m water equivalent), Sub is the daily sublimation (m water equivalent), and A is the area of the elevation band (m 2 ).
We ran the integrated model over two time periods. From 1998 to 2018, we drove the model with meteorological observations from within the Shullcas Watershed for model calibration and analysis of hydrograph components. During this period, we also ran the model with and without glacier melt to determine glacier melt contribution to groundwater discharge. From 1961 to 2099, we drove the model using downscaled GCM outputs for hydrologic projections.

Model Calibration and Sensitivity Analysis
We calibrated selected input parameters (using >100 model runs) to minimize an error metric (equation (3)) that combines multiple data sources to ensure that the model adequately captures relevant processes across the modeled domain. Twenty key parameters were varied within a range of possible values (Table S1). The root-mean-square error (RMSE) was calculated for each available streamflow record and included in an aggregate error, E Agg : where E [1998][1999][2000][2001][2002][2003][2004][2005][2006][2007][2008][2009] is the RMSE of the monthly discharge record from ANA from 1998-2009 and E 1-6 are the calculated RMSE for each of the six daily discharge stations between 2015 and 2018.
We used the Nash-Sutcliffe Efficiency (E N-S ) to express model fit, where 10.1029/2019GL084730

Geophysical Research Letters
SOMERS ET AL. 13,019 Q obs and Q sim are observed and simulated stream discharge (m 3 /s), respectively, at any time step i.
Lastly, we used sensitivity analysis to establish error bounds based on four parameters that we identified as relatively uncertain given our inability to measure them across the modeled domain: (1) vertical and horizontal hydraulic conductivity of bedrock, (2) Priestly-Taylor α coefficients for PET, (3) interflow coefficients, and (4) glacier melt factors.
Given the computational intensity of the model, we use a computationally efficient, one-at-a-time sensitivity analysis (Zheng & Bennett, 2002) to assess the uncertainty of the results (see supporting information).

Climate Projections
Climate projections were generated from the dynamic Eta regional climate model (Chou,  Daily temperatures were corrected using the additive (classical) bias correction technique (Bordoy & Burlando, 2013;Hawkins et al., 2012) while daily precipitation raw data were corrected in two steps. First, linear scaling (Smitha et al., 2018), also called multiplicative shift (Ines & Hansen, 2006), bias correction were used, with a three-month sliding window precipitation. Second, the quantile mapping bias correction technique (Ines & Hansen, 2006) was applied following the procedure by Grillakis et al. (2013). The biascorrected RCM results produce comparable glacier melt amounts to the local meteorological data.

Model Calibration and Sensitivity Analysis
Our calibrated, integrated hydrological model adequately represents field observations where E N-S for stream discharge ranged from −0.46 to 0.74 ( Figure S1 and Table S2), with four of the six gauging stations having E N-S greater than 0.5. A hydrological model is often considered to adequately represent observed data if E N-S > 0.5 (Moriasi et al., 2007). Furthermore, temporally averaged modeled water-table depth closely matches observations from our two groundwater well clusters with root-mean-square errors (RMSE) of 0.20 and 1.50 m for the Chuspicocha and Huacracocha areas, respectively, from 2016 to 2018. The model underestimates the observed range in water table fluctuation, in part due to the large grid cells which may not accurately represent the point locations of the wells which are all located in wetlands close to the river for accessibility. Future work should seek to expand groundwater head data, particularly in highlands where piezometer and well data are extremely scarce.
Dry and rainy season glacier melt factors of 2.8 and 3.0 mm d −1°C−1 , respectively, yielded the best fit to observed meltwater production. The relatively low melt factors are consistent with the dry, high-radiation setting where more of the available energy is partitioned to sublimation (Fernàndez & Mark, 2016;Hock, 2003).
The climatic setting of the tropical Andes lends itself well to the calibration of a joint surface-water and groundwater model. Because very little precipitation falls during the dry season, the baseflow recession pattern and dry-season streamflow serve as ready calibration targets for the groundwater-flow parameters. Meanwhile, the rainy season hydrograph serves as a good calibration target for the interflow parameters.
Sensitivity analysis indicates that model results are most sensitive to the Priestly-Taylor α coefficient, followed by the hydraulic conductivity of the fractured bedrock layer, the interflow coefficient, and the glacier melt factors (Table S3). Simulated streamflow has an uncertainty of ±23% for a 95% confidence interval (see supporting information).

Groundwater Dominates Dry Season Streamflow
From 1998 to 2018, simulated historical groundwater discharge from the saturated zone to the stream contributes an average of 37% (ranging from −4.4 to 93%) of outlet streamflow annually and 72% during the dry season ( Figure 3a). During the rainy season, interflow (shallow, preferential flow through the unsaturated zone) dominates (60%) streamflow. These results complement work by others (e.g., Andermann et al., 2012;Baraer et al., 2015;Hood & Hayashi, 2015;Pohl et al., 2015;Saberi et al., 2019;Somers et al., 2016) to demonstrate that groundwater is an important contributor to streamflow in glacierized mountain catchments.
Conversely, our calibrated model simulations indicate that glacier melt contributes only 8.0% of Shullcas River outflow discharge from 1998 to 2018. To validate this result, we compare the model outputs against hydrochemical mixing analysis for streamflows from 10-22 July 2014 (Crumley, 2015). Our simulated 11% glacier-melt contribution to stream discharge for this specific time period lies within the geochemically estimated 8.9-17% streamflow contribution.

Glacier Melt Contributes Little to Groundwater Recharge
Only a small amount (2.3%) of saturated zone groundwater discharge to the Shullcas River was originally sourced from glaciers that covered 2.4% of the watershed area in 2015. The low contribution of glacier melt to groundwater suggests that overall, groundwater recharge is resilient to a reduction in glacier meltwater input. To our knowledge, the only other quantitative estimate of glacier melt contribution to groundwater comes from a headwater catchment on Volcán Chimborazo in Ecuador, where 15% of groundwater discharge was found to originate from glaciers that cover 34% of the watershed area (Saberi et al., 2019). These results suggest a useful relationship between glacier coverage as a proportion of watershed area, and glacier-sourced groundwater discharge, similar to the relationship proposed by Baraer et al. (2015) between glacierized area and glacier melt contribution to streamflow. This relationship likely varies with climate, bedrock lithology, and topography (Markovich et al., 2016;Saberi et al., 2019;Tague et al., 2008), and requires more estimates in different basins to expound.
The simulated groundwater table (Figure 2b) remains close to the ground surface in valley bottoms and near streams, but is located as deep as 176 m below the ground surface in highland and steep areas (similar spatial trends to Ofterdinger et al., 2014), well below the 20-m-thick unconsolidated surficial aquifer. Therefore, much of the modeled groundwater flow toward the streams passes through the bedrock aquifer, as suggested in other mountain hydrological studies (Andermann et al., 2012;Liu et al., 2004;Tague et al., 2008). The longer and slower flow paths through bedrock allow groundwater to contribute to streamflow year round.

Unsaturated Zone Storage Fluctuates Seasonally While Saturated Groundwater Storage Follows Decadal Trends
Subsurface water storage fluctuates seasonally, increasing during the rainy season and decreasing during the dry season. For example, in 2016-2017 the range in cumulative annual storage change is largest in the deep unsaturated zone (0.0158 km 3 ), followed by the shallow unsaturated zone (0.0095 km 3 ) and the saturated zone (0.0060 km 3 ; Figure 3b). For comparison, the nearby city of Huancayo requires approximately 0.023 km 3 of water per year, based on the 175 L per person-day estimated by the United Nations (United Nations Development Report, 2006).
The timing of changes in storage is important for water resource management. The unsaturated zone is quickly depleted after the end of the rainy season, whereas the saturated zone continues to drain and supply water to streams throughout the dry season, during maximum water stress. The negative net change in saturated storage in a particular year (Figure 3b) illustrates the decadal time scale on which groundwater fluctuates and therefore its ability to buffer interannual climatic variability (Cuthbert et al., 2019), even in a mountain environment where hydraulic gradients are high.

Rapid Projected Climate Warming Causes Disappearance of Huaytapallana Glaciers
Observed temperatures at the Huayao Meteorological Station from 1975 to 2005 increased at a rate of 0.42°C per decade. The averaged projected rates of warming from the Eta RCM over the twenty-first century are 0.46 and 0.88°C per decade for RCP 4.5 and 8.5 emission scenarios, respectively (Figure 4a). The high projected rates of warming are consistent with elevation-dependent warming (Mountain Research Initiative EDW Working Group, 2015;Seth et al., 2010;Urrutia & Vuille, 2009;Yarleque et al., 2018). For the RCP 4.5 scenario, minimum and maximum temperatures begin to stabilize around 2070. For RCP 4.5, little change is projected in precipitation and the GCMs disagree on the direction of change, whereas RCP 8.5 projects more consistent negative trends in precipitation.
All emission scenarios and GCMs run through our glacier melt module project the complete disappearance of the Huaytapallana glaciers by the end of the twenty-first century (Figure 4a), with the slowest melt scenario projecting disappearance (glacier area <100 m 2 ) in 2085. Our glacier projections are comparable to equilibrium-line-elevation projections for the Quelccaya ice cap (Yarleque et al., 2018), approximately 400 km southeast of our study area, and with a similar summit elevation (5,680 versus 5,557 m for the Huaytapallana glaciers). Although local climatic conditions vary, approximately 2,330 glaciers in the Peruvian Andes have summit elevations of less than 5,500 m (as of 2014, calculated from ANA (2014)) and therefore face similar outlooks.
Lacking historical humidity and solar radiation records onsite, our temperature-index glacier module provides a reasonable estimate of glacier melt (Fernàndez & Mark, 2016) with the available temperature data and compares well (R 2 = 0.67) with observed historical glacier areas by López-Moreno et al. (2014) (Figure 4a). However, the glacier module does not include edge-effect feedback, in which longwave radiation from rock surrounding the glacier plays an increasingly important role as the glacier retreats (Aubry-Wake et al., 2015), or elevation feedback, in which glacier thinning lowers the elevation of the glacier surface and is subjected to higher temperatures (Weertman, 1961). These unaccounted feedback make our glacier melt projections conservative estimates of mass loss through time.

Groundwater Temporarily Buffers Loss of Glacier Melt Under Future Climatic Conditions
As glaciers recede, meltwater production initially increases to a maximum, known as "peak water," before dropping off as the glaciers continue to retreat (Gleick & Palaniappan, 2010). Approximately half of the world's glacierized drainage basins have passed peak water (Huss & Hock, 2018), including >70% of studied proglacial watersheds in the Cordillera Blanca, Peru (Baraer et al., 2012). Our model simulates multiple peaks in meltwater production, the last of which occurs when glacier areal retreat is fastest, around 2013, and then steadily drops off until it becomes zero when the glacier completely disappears (Figure 4b).
While the projected decrease in glacier melt in turn decreases stream discharge, it is partially mitigated by groundwater discharge to the stream. During the dry season, glacier melt steadily decreases after peak water. Groundwater discharge to the stream, on the other hand, remains relatively constant (Figure 4b), thereby maintaining a buffer to stream flow during the water-stressed dry season. However, as temperature increases, so does ET, which reduces groundwater recharge. Furthermore, changes in the frequency of extreme events and rain-snow fractions are implicitly included in the climate projections and also affect the amount of groundwater recharge. The effect of decreasing recharge is observed in the gradual decreasing trend in dry season groundwater discharge (Figure 4b). This effect is more severe in the RCP 8.5 scenario where temperature continues to increase through 2099, and where there is a severe reduction in precipitation for some GCMs.
Although both RCP 4.5 and 8.5 predict complete loss of glacier meltwater by the end of the twenty-first century, there is a marked difference in streamflow between the scenarios. RCP 4.5 (RCP 8.5) projects 20% (50%) lower streamflow and 18% (35%) lower groundwater discharge year round due to both lower precipitation and higher ET (Table S4 and Figure 4). Both future scenarios project that a larger proportion of streamflow will come from groundwater. These projections do not account for changes in ET as vegetation colonizes the deglaciating area. However, the additional ET from the~2% glacierized area of the watershed is assumed to be relatively minor considering the cold, steep, and rocky environment.
All of the modeled climate scenarios project a decrease in water table elevation and groundwater storage (Table S5). On average, the groundwater table drops by 8.8 and 13.4 m for RCP 4.5 and 8.5, respectively, over the twenty-first century. The greatest changes in water table elevation occur beneath the Huaytapallana glaciers, where the loss of meltwater decreases groundwater recharge, and beneath the steepest hillslopes,  Figure 1a). Peak glacier runoff occurs in 2013. (c) Change in modeled mean monthly streamflow between the first and last decade of the twenty-first century. (d) Change in groundwater discharge between the first and last decade of the twenty-first century. (e) The proportion of streamflow that comes from groundwater discharge, where a value of 1 means that all streamflow is sourced from groundwater discharge, demonstrating the seasonality of groundwater importance annually.
where hydraulic gradient is greatest. The least amount of change occurs near the stream network and in the high flat areas of the watershed. Similar spatial trends in groundwater table sensitivity to recharge have been suggested in the Alps (Ofterdinger et al., 2014).

Conclusions
Though historically considered a minor component of mountain streamflow (Liu et al., 2004), our results demonstrate the current and future importance of groundwater in a mountain hydrologic system. Groundwater dominates (72%) dry season streamflow, including substantial contribution from fractured bedrock aquifers, but receives little recharge from glacier melt (2%), making groundwater storage resilient to decreases in glacier melt. In the near term (~30 years), consistent groundwater discharge mediates the loss of glacier melt. But, in the long term (>60 years), increasing ET and decreasing precipitation (for some GCMs) decrease groundwater recharge and therefore stream discharge, particularly for the high carbon emission scenarios (Figure 1). The relatively small (~0-2 m) projected drop in water table beneath flat areas may lend some resilience to high-altitude wetlands located in these areas. However, Andean wetland ecosystems, considered biodiversity hot spots, may be quite sensitive to changes in water table (Maldonado Fonkén, 2014).
More broadly, our results demonstrate that mountain water resource response to climate change is shaped not only by glacier change but also by surface and-particularly-groundwater hydrological processes. These findings provide important context to large-scale glaciological studies (e.g., Huss & Hock, 2018). Future efforts to model mountain hydrology under climate change should include a realistic groundwater component.