The second largest concentration of glaciers in the U.S. Rocky Mountains is located in Glacier National Park (GNP), Montana. The total glacier-covered area in this region decreased by ∼35% over the past 50 years, which has raised substantial concern about the loss of the water derived from glaciers during the summer. We used an innovative weather station design to collect in situ measurements on five remote glaciers, which are used to parameterize a regional glacier melt model. This model offered a first-order estimate of the summer meltwater production by glaciers. We find, during the normally dry month of August, glaciers in the region produce approximately 25 × 106 m3 of potential runoff. We then estimated the glacier runoff component in five gaged streams sourced from GNP basins containing glaciers. Glacier-melt contributions range from 5% in a basin only 0.12% glacierized to >90% in a basin 28.5% glacierized. Glacier loss would likely lead to lower discharges and warmer temperatures in streams draining basins >20% glacier-covered. Lower flows could even be expected in streams draining basins as little as 1.4% glacierized if glaciers were to disappear.
Mountain glaciers in North America are reservoirs of snow and ice that store winter precipitation and release meltwater during summer months. Water from snow and ice melt can partially compensate for the lack of precipitation during a hot, dry year (Rothlisberger and Lang, 1987) and help maintain base flows and cool water temperatures in streams that are critical for aquatic organisms (Jacobsen et al., 2012). Glaciers can moderate the effects of multi-year droughts and, during hot summers, add significantly to seasonal discharge, particularly during late summer in the western United States (Meier and Tangborn, 1961; Fountain and Tangborn, 1985). Previous studies demonstrate that even very small glaciers can have disproportionately large effects on basin runoff. For example, Nolin and others (2010) concluded that glacier-derived runoff made up 41%–73% of the late summer discharge in the Upper Middle Fork of the Hood River, Oregon, which drains a basin only 6.6% glacierized. Huss (2011) calculated glacier-melt contributions to rivers draining the European Alps and found that even in basins that were only 0.20% and 0.06% glacierized, glacier runoff supplied 6.6% and 2.8% of the total August discharge, respectively. In the Canadian Rockies, glacier-melt can account for as much as 27% of the total runoff in basins as little as 1% glacierized (Comeau et al., 2009). Hence, glaciers can play a key role in the hydrology of the basins they occupy by modulating the timing and volume of runoff.
The ∼5000 km2 landscape within and bordering Glacier National Park (GNP) in northwest Montana contains the second largest concentration of glaciers in the U.S. Rocky Mountains (Fig. 1). A total of 39 named glaciers cover 17.2 km2 (Table 1). Cold water issuing from glacier outlet streams is considered integral to this region's ecosystems, especially during the late summer months (Pederson et al., 2010). The meltwater stonefly (Lednia tumana) is endemic to GNP and lives almost entirely in stream reaches that are within 500 m of melting snow and ice. It is now a candidate species for listing under the Endangered Species Act due to the predicted warming climate and subsequent loss of glaciers and permanent snowfields (Muhlfeld et al., 2011). Bull trout (Salvelinus confluentus), another native species in this region, require colder water compared to other North American salmonids (Selong et al., 2001). Due to projected impacts of a warming climate along with other stressors such as habitat degradation/fragmentation and invasive species (Jones et al., 2013), bull trout have already been listed as “threatened” under the Endangered Species Act (USFWS, 2010). Humans also rely on this region's glacier meltwater. Irrigation in central Montana consumes 50%–90% of the discharge from the Saint Mary River (USBR, 2012) that drains a basin containing two-thirds of the total glacier area in GNP east of the Continental Divide. Glaciers west of the Continental Divide in GNP lie within the headwaters of the Columbia River watershed, a basin containing 31 hydroelectric facilities.
The retreat of glaciers in GNP has come to epitomize the impacts of a warming climate on the landscape and the hydrologic cycle of the western United States. Recession of glaciers over the past 100 years has been well documented here (Dyson, 1948; Johnson, 1980; Carrara and McGimsey, 1981) with rates of retreat being higher than in other U.S. mountain ranges (Fountain, 2007). With continued climate warming, glaciers are expected to continue shrinking or even disappear. A geospatial model projection suggests the disappearance of five glaciers in GNP by 2030 (Hall and Fagre, 2003), while another process-based model of one glacier suggests current conditions could cause elimination by about 2080 (Brown et al., 2010). Substantial habitat loss for native cold-water aquatic species and possible extinctions are projected for the GNP region if glaciers were to disappear (Pederson et al., 2010; Muhlfeld et al., 2011). Other North American studies (Moore et al., 2009) have shown that glacier retreat would result in less water for hydroelectric power generation and increased competition for water used to irrigate crops.
Despite previous research describing the possible consequences of diminishing glaciers, no study has quantified glacier melt runoff in GNP explicitly, leaving uncertain the significance of this region's glaciers on hydrological and ecological processes. This uncertainty is in part due to an almost complete lack of observational data from the remote and difficult to access glaciers and their associated outlet streams. Mass balance data are available for just one glacier in the study area, so quantifying region-wide glacier runoff from glaciological techniques (Ostrem and Brugman, 1991; Kaser et al., 2003) is not possible. There is only one gaged glacier outlet-stream, but it is rarely operational, with the last period of record spanning 2004–2009. This precludes determining glacier-derived runoff by hydrological methods (e.g., Pellicciotti et al., 2010). These problems are common to other glacierized areas worldwide, so many studies model glacier melt as a proxy for runoff (Hock, 2005; Comeau et al., 2009; Huss, 2011; Racoviteanu et al., 2013) to name a few. The complete water balance in a basin has many components including, but not limited to: groundwater flow, evaporation, transpiration, and precipitation. Data for these components is also lacking or non-existent and difficult to measure across GNP. Therefore, a method was needed to make a simple, first-order estimate of the glacier melt component in streams draining glacierized basins.
List of August runoff and specific discharges for the 39 study glaciers.
In this paper, we develop a glacier-melt model parameterized by strategic in situ glaciological measurements taken on five selected glaciers. The model then quantified glacier meltwater production across GNP. From these results we present necessary estimates of glacier-derived runoff during the month of August in five gaged streams and also determine a relationship between these estimates and the percentage of glacier-covered area in each stream's catchment. The discussion addresses both the strengths and limitations of our methods and the prospect of using them in similar data-sparse and difficult to access regions. We also interpret the significance of our results in light of the current trend of receding glaciers.
Geographic data for weather stations.
We model melt on 39 named glaciers, which include 37 within GNP, and 2 more, located 10 km from the park's western boundary (Fig. 1; Table 1). These named glaciers are the largest in the study area and contain the majority of the snow and ice in this region persisting through the summer months and from year-to-year. During any given year, there are an undetermined number of unnamed much smaller snowfields and ice masses persisting until autumn. These are not included in this study.
Most glaciers in GNP are wider than they are long relative to the flow direction, and are classified as cirque and niche glaciers (Key et al., 2002). Meltwater typically issues from many small outlets and often discharges over cliffs. The glaciers are located primarily above 2000 m elevation and near the Continental Divide in steep east- to northeast-facing cirques. They likely accumulate much of their snow mass from avalanching and wind-transport (Graf, 1976; Kuhn, 1995; Allen, 1998). The 2005 average size of the glaciers is 0.44 km2 but areas range from 0.02–1.89 km2. Only Sperry Glacier (∼0.87 km2) has a surface mass balance monitoring program where glacier-wide (Cogley et al., 2011) balances are calculated using glaciological techniques (Ostrem and Brugman, 1991). The program is run by the U.S. Geological Survey (USGS) and began in 2005. Snow depth and density measurements along with seven ablation stakes are used to estimate seasonal and annual balances each year.
STUDY TIME INTERVAL
We ran the glacier melt model for months July, August, and September, but focus our results on the month of August since this is the month when glacier runoff has potential to be most critical to the hydrograph in GNP. The seasonal snowpack on non-glacier-covered regions can linger well into July above 2000 m (Gillan et al., 2010), with peak snowmelt runoff in late June, but is typically completely gone by August. August is the second driest and second warmest month of the year, so precipitation contributes little to stream flow (Finklin, 1986). By September, mean temperatures decline and precipitation increases, both changing by 30% from the mean August values of 16 °C and 4.5 cm (Finklin, 1986) and glacier melt shows a marked decline. We modeled two years, 2009 and 2010. Long-term temperature records from Flattop Mountain located inside GNP and Kalispell Airport located 25 km southwest of GNP (Table 2) show 2009 was warmer than average and 2010 was cooler. However, both years were within 1 standard deviation (SD) of the mean, and the 2009/2010 average was within 1% of the mean for the entire period of record. Given this temperature record, we assume that August glacier melt was high in 2009 and low in 2010. USGS ablation measures from Sperry Glacier (USGS unpublished data) support this assumption as ablation was above average in 2009 and below average in 2010. We then averaged the model results from these two contrasting years and reasoned that these mean values would represent meltwater production during an August with average temperatures.
IN SITU MEASUREMENTS
We installed on-ice meteorological stations on five different glaciers: Blackfoot, Grinnell, Pumpelly, Sexton, and Sperry (Figs. 1 and 2). The five glaciers were chosen to sample the spectrum of elevations, aspects, and topographic shading of the 39 glaciers (Table 3). The remote glaciers could only be accessed on foot, requiring that instrumentation be lightweight and simplified. We built a tower on each glacier consisting of three vertical support legs with two horizontal cross bars (Fig. 3). All materials needed to be backpackable, so the tower structure was cut into seventeen 1.50-m-long sections of 0.025 m diameter aluminum pipe. We built each leg by fastening five sections together in a telescoping fashion using roll pins and a 0.15-m-long aluminum rod at the joint. The legs were then placed vertically into 7 m deep holes that were drilled into the glacier in a tripod pattern using a backpackable steam drill. Two horizontal cross-bars, placed in a “T” shape, connected the three legs above the glacier's surface. All instruments were then fastened to these cross-bars. Power was supplied by a battery and recharging solar panel.
Measurements were made at 10 minute intervals and then stored in the data logger as hourly averages. A 0.1 °C resolution thermistor inside a radiation shield measured temperature, and a pyranometer measured incoming solar radiation with 1%–4% accuracy at 300–1100 nm. Melt rates were inferred from measurements of glacier surface lowering. This was obtained using a sonic distance ranger aimed at the glacier. The height loss measured by the sonic ranger was multiplied by the density of snow or ice and provided a water-equivalent depth for the mass lost. As the glacier's surface lowered beneath the station, the horizontal cross-bars could be dropped down the vertical legs without unfastening the instruments in order to keep the increasingly top-heavy structure from falling over. Fifteen ablation stakes were also installed (Fig. 2) to infer melt rates from surface lowering on variable multi-day intervals at locations different from the weather stations. Stations were in place from early June through mid- to late September, depending on the glacier. We visited each glacier on 3–5 week rotations to lower the instruments on their support poles, and to measure surface lowering at each ablation stake.
Geographic data for on-glacier weather stations with melt rates*.
To measure snow density, we dug four snow pits, one each on Sexton, Pumpelly, and Grinnell Glaciers during mid-July, 2010, and an additional snow pit on Sexton Glacier in mid-July, 2011. Snow pits were located within 10 m of the weather stations, and depths ranged from 1.0 m to 1.2 m. Density measurements were made at 0.10 m intervals. We also referenced data from the USGS Sperry Glacier mass balance program and other studies (LaChapelle, 1954; Pelto, 1996; Krimmel, 1998; Miller and Pelto, 1999) to compare results.
In 2011, we measured water temperature at 15 minute intervals with a HOBO water level logger in an outlet stream within 50 m of the terminus of Sexton Glacier starting on 4 August and ending on 8 October. Other sources of water temperature in streams came from the USGS stream gage on the North Fork of the Flathead River (Fig. 1; Table 4) where records date from 1997-present, and a recent study done in the North Fork of the Flathead basin (D'Angelo and Muhlfeld, 2013).
Data for stream gage sites.
Our goal was to obtain a two-year snapshot spanning most of the summer ablation season for the years 2009 and 2010 revealing the meltwater produced by melting of snow, firn, and ice. To do this, we quantified glacier-derived runoff by simulating surface melt on the glaciers as delineated on a digital topographic domain obtained from the USGS National Elevation Dataset (NED). A 30 m digital elevation model (DEM) was downloaded from NED then resampled to 60 m resolution for the sake of computational efficiency. The margin of each glacier was determined by hand-digitizing polygons using aerial photographs obtained from the National Agricultural Image Program (NAIP). These photos were taken from 21–31 August 2005, after the previous winter snow cover had disappeared from the landscape, thus clearly revealing each glacier's margin. Some outlines were field checked for accuracy by USGS employees during the late summer, years 2006–2009.
We simulated daily melt using a regionally distributed temperature-index model modified to include solar forcing. On large spatial and temporal scales, empirical models such as this often perform just as well as more sophisticated physically based energy balance models, especially in regions similar to GNP where input data is lacking (Ohmura, 2001; Hock, 2003). Compared to traditional temperature-index and/or degree-day models, our model includes the addition of a separate solar radiation term such that
Here, M is melt (m d-1), a is an empirically derived temperature coefficient ([m d-1] [°C-1]-1), T is a temperature index representing mean daily air temperature (°C), Tc is the critical temperature for melt (0 °C), SR is a solar radiation index (W m-2) representing the cumulative effect of short-wave solar radiation, b is an empirically derived radiation coefficient ([m d-1] [W m-2]-1). The solar radiation component improves model performance because incoming solar radiation is a major source of melt energy treated independently of air temperature (Hock, 1999; Pellicciotti et al., 2005).
Mean daily air temperatures (T) for each grid cell were distributed across the model domain using data from local weather stations (Fig. 1; Table 2) and an inverse-distance weighting interpolation adjusted with a locally calculated lapse rate of 0.0074 °C m-1 (Dodson and Marks, 1997; Gillan et al., 2010). Inputs to our simulation of 2010 include data from the five on-glacier weather stations. The radiation index for each grid cell was obtained by calculating the potential, direct solar radiation received for each hour in every grid cell (Hock, 1999). These calculations account for northing and easting coordinates, time of year and day, slope angle and aspect, as well as shading effects due to the surrounding topography. Hourly values were then summed for each day. This value, a clear sky potential radiation index, was then reduced each day to account for cloud cover by using solar radiation data from local weather stations and calculating a daily scaling factor. In 2010, this scaling factor was derived from solar radiation measurements collected at five on-glacier weather stations. For example, on 15 August 2010 the station on Blackfoot Glacier received 77% of the potential radiation, Grinnell 60%, Pumpelly 73%, Sexton 79%, and Sperry 95%. We then averaged these percentages, and so the scaling factor for this day was 0.77. Next, every potential radiation value in each grid cell was multiplied by 0.77, generating the radiation indices in each grid cell (SR in Equation 1) for 15 August 2010. It is important to note that the solar radiation term does not give the equation an energy balance component. Rather, it is an attempt to separately account for the influence of solar radiation on glacier melt using an empirical relationship.
Temperature and solar radiation, and melt data from the on-glacier weather stations were used to derive the coefficients α and β in Equation 1 iteratively by optimizing for the lowest sum of square errors such that
Here, Φ is the sum of the squared errors between the daily melt, measured in meters of water equivalent, on the ith day (melti) and the modeled melt. For modeled melt, tempi is the measured mean daily air temperature on the ith day, and solari is the solar radiation index for the ith day derived from measurements. The optimal values for the month of August were 0.0038 m d-1 C-1 for α and 3.93 × 10-6 (m d-1) (W m-2)-1 for β. We specifically performed regressions of average daily air temperature and the daily solar radiation index taken from the five stations. This revealed that the two variables are not collinear. R2 values from the five weather stations ranged from 0.28 to 0.53. When all data from the stations were combined into one set and then regressed, R2 = 0.32.
We did not use different coefficients for snow, firn, or ice to represent differing albedo and heat transfer processes. We acknowledge that some other studies do acquire separate coefficients, or melt factors, for these surfaces (Hock, 2005; Shea et al., 2009). However, we had no reliable method to discern what portions of the glacier-covered area across GNP was partitioned into snow, firn, or ice on a daily basis. To account for the differences in melting due to varying surfaces, we used different coefficients on a monthly time step. We combined all measurements from the glacier weather stations to obtain coefficients that best-fit each month of the melt season: July, August, and September. In this way the resulting coefficients better represent the melting of mostly snow-covered glaciers in July and a mix of snow and ice surfaces that were melting in August and September.
We applied Equation 1 at daily time steps to the 60 m digital landscape. It is important to note that for the purposes of our study, we include melt derived from both the seasonal snow component and the firn and ice components of glacier surface ablation as glacier-derived runoff. Equation 1 yielded a modeled surface meltwater depth that was then multiplied by the area of each grid cell to obtain a volume. Summation of these grid cells for each glacier yielded a glacier-wide volume of potential daily runoff.
COMPARISON TO GAGED STREAMFLOWS
Our aim was to provide a baseline point of reference for the potential relevance of glacier-derived runoff. The model provides a daily volume of meltwater generation for each glacier, but does not account for any sinks to this meltwater such as evaporative losses and losses to groundwater flow. We also do not calculate a lag time for meltwater to flow downstream, through lakes, and/or through the groundwater system to reach the stream gages. Hence, our results are estimates that represent the maximum possible runoff derived from glaciers. This effectively sets an upper limit on glacier contributions to streamflow.
These contributions to the total basin runoff were assessed by comparing the modeled glacier runoff with total basin runoff measured by USGS stream gages. For each gaged basin, daily runoff from the glaciers within it was summed for August and divided by the total basin runoff measured at the stream gage.
The study area is divided into five gaged watersheds (Fig. 1; Table 4): North Fork Flathead River, Middle Fork Flathead River, Saint Mary River, Swiftcurrent Creek, and Upper Grinnell Creek. Two other watersheds, the Belly and Waterton River basins, are 1.0% and 0.5% glacierized, respectively, but these basins lack recent stream gage data and so were unusable in this study. The North and Middle Forks of Flathead River are the two largest watersheds in our study. Located west of the Continental Divide they also have the smallest percentages of glacier-covered area. East of the Continental Divide, the Saint Mary watershed is the third largest basin, and water drains northeast to the South Saskatchewan River and eventually to Hudson Bay. The Swiftcurrent Creek catchment is a sub-basin of the Saint Mary watershed.
The Upper Grinnell Creek catchment is further nested within the Swiftcurrent Creek basin and has the smallest catchment area and the highest percentage of glacierized area, currently at 28.5%. A stream gage sits within 500 m of Grinnell Glacier, but a 0.33 km2 pro-glacial lake sits between it and the terminus. Only two periods of stream discharge data exist, one spanning 1959–1971 and the other 2004–2009. Thus, we were unable to make a comparison between modeled glacier runoff to total basin runoff for August 2010. Despite these problems, this site is valuable because it is the only gage inside the study area with discharge data from a glacier outlet stream.
During August, the remaining seasonal snow not within the glacierized area reduces to less than 25,000 m2, <4% of the total snow/ice-covered area in the basin. There are no other streams feeding the catchment because it sits on the Continental Divide. Precipitation in August is quite low. Data from the Many Glacier SNOTEL site (records 1978-present), the closest weather station and located 6 km from Grinnell Glacier, show annual precipitation ranges from 0.72 m to 1.81 m with a mean of 1.18 m with a standard deviation (SD) of 0.23 m. Values for August range from 0.00 m to 0.14 m with an average of 0.05 m (SD 0.04 m) meaning less than 5% of the annual accumulated precipitation typically occurs in August.
There is a high probability that glacier-melt is the dominant source of discharge in this stream during the month of August. The surface of the basin upstream from the gage has very little soil and consists mostly of exposed bedrock. Although it has never been studied explicitly in this basin we assume groundwater inputs and losses in the catchment to be negligible. It is also possible that some of the August discharge is sourced from rainfall and snowmelt that occurred earlier in the summer and is stored in the small pro-glacial lake.
IN SITU DATA
The five on-glacier stations provided a new and comprehensive data set of meteorological conditions and surface melt on daily time steps. Removal of elevation effects by normalizing average daily temperatures to the average elevation of the stations (2287 m) and using the local lapse rate indicates that mean daily temperatures vary little between the stations (Fig. 4, part a). However, variability among stations is greater on an hourly time scale (e.g., Fig. 4, parts b and c) and can change by as much as 15 °C between days.
The mean melt rate during 27 June to 26 September 2010 was 0.035 m d-1 of water equivalent (w.e.) (SD 0.016 m d-1 w.e.) (Table 3). The cumulative daily melt from the five stations during this period (Fig. 5) showed spatial and temporal variability. From late June to late July, the mean melt rate of the five stations was 0.043 m d-1 w.e., and declined to 0.037 m d-1 w.e. during late July to late August. Mean daily melt was reduced to 0.017 m d-1 w.e. from late August to late September if only melting days are considered. If all days during this interval are included in the mean value, then the glaciers actually accumulated water during a few late summer snowfalls by an average rate of 0.003 m d-1. The sonic on Pumpelly Glacier only worked from 10 August to 11 September, so we used an average rate from the other four stations of 0.040 m d-1 w.e. to generate the Pumpelly curve up to 10 August.
Melt rates estimated from ablation measurements at the 15 stakes generally agree with results from the sonic rangers, although comparisons are limited by the coarse time resolution of the stake measurements. Stake averages from late June to mid-July showed the melt rate was 0.036 m d-1 w.e., rising to 0.044 m d-1 w.e. during mid-July to mid-August, and then declined to 0.017 m d-1 w.e. from mid-August to late September. The summer season average (27 June to 26 September) using the combined data was 0.033 m d-1 w.e.
We found that snow density varied little with increasing depth in the snowpack, nor did it change greatly across the landscape. For example, at our highest elevation site on Pumpelly Glacier (2525 m a.s.l.) the average bulk density on 15 July 2010 was 562 kg m-3 (SD 51 kg m-3). Meanwhile, on Grinnell Glacier, our lowest elevation site at 2058 m a.s.l. and 20 km north of Pumpelly Glacier, the average bulk density on 17 July 2010 was 554 kg m-3 (SD 30 kg m-3). On Sexton Glacier, density measurements were made on 10 July 2010 and on 11 July 2011, and we found average values of 544 kg m-3 (SD 48 kg m-3) and 541 kg m-3 (SD 22 kg m-3), respectively. Given the low variability found in our density measurements, we calculated an average bulk snow density of 550 kg m-3 (SD 36 kg m-3) from the four snow pits and then used this value for all subsequent calculations of snow water equivalents.
The low variability found in our results is consistent with the USGS measurements made on Sperry Glacier during June, years 2005–2014. Here the mean density was 530 kg m-3 (SD 62 kg m-3) from snow pits ranging in elevation from 2360 to 2525 m. Glacier mass balance programs in both Washington and Alaska have also found that snow density varies little across the landscape in a late summer snowpack that has been sitting for several weeks at melt conditions (LaChapelle, 1954; Pelto, 1996; Krimmel, 1998; Miller and Pelto, 1999). Research on Storglaciären, Sweden, found that using the bulk mean density obtained from snow pits was a reliable method for determining snow-derived water equivalents and that other methods, such as modeling density with depth using various functions did not yield significantly different results (Jansson, 1999).
Water temperature in the Sexton Glacier outlet stream for the month of August 2011 averaged 0.43 °C with a SD of 0.19 °C. Temperatures ranged from 1.11 °C to 0.23 °C. August water temperature at the stream gage in the North Fork of the Flathead River averages 15.0 °C with a SD of 1.5 °C for the period of record.
To investigate the sensitivity of our results to our optimized coefficients αand β, we ran the model using a range of coefficients. Values for αand βwere adjusted at 2% increments from -10% to +10%, and the model was run 20 times with random pairings of values (Fig. 6). The model is more sensitive to changes in α than to changes in β. For example, a 10% increase in α with a 10% decrease in βgives an overall 3% increase in August runoff. When both coefficients were increased or decreased by 10%, daily meltwater production from all glaciers increased/decreased by 10%.
In grid cells containing the on-ice weather stations, we examined the sensitivity of model results to the model's inputs. We compared the original temperature and solar radiation indices to another set of indices; however, these were derived only from the gridding process. First, we omitted data from the five on-ice stations from the gridding scheme and then interpolated new temperature and solar grids. We then compared the new modeled indices from cells with weather stations to the original, measurement-derived indices. Comparisons of both temperature and solar radiation indices showed the residuals were normally distributed. For 53% of days the measured average daily temperatures were warmer than the gridded values. The average difference was 2.0 °C (SD 1.4 °C). Temperatures were colder on the remaining 47% of days by an average of 1.7 °C (SD 1.7 °C). The gridded solar radiation index differed from the measurement-derived solar radiation index by an average of 500 W (SD 400 W). Again, periods of more/less solar radiation were split almost equally, with 51% of the days receiving more radiation than modeled and 49% of the days receiving less.
The uncertainties in temperature and radiation inputs can combine to produce a total input uncertainty. If modeled temperature and solar radiation typically differ from actual conditions by 2 °C and 500 W, this translates to a difference of ±0.01 m d-1 in meltwater equivalent on each grid cell. If this difference is assumed constant across the study area, the total glacier runoff would change by 24% of the 2010 total and 19% of the 2009 total, both of which are substantial. However, this estimate assumes both the temperature and radiation indices are concurrently and uniformly too high or too low every day, which has a low probability of occurring. Our data show that over the entire month, periods of positive and negative differences between model input and observations are almost equal, largely creating a “cancellation effect” where melt is over- and undercomputed about equally. Hence the calculated estimates of uncertainty above are an upper limit and the realized total input uncertainty is likely much less.
We compared our in situ measurements of glacier melt at each ablation stake to the modeled melt from grid cells containing the stakes. We acknowledge this comparison suffers from point-to-cell scaling issues. In 2010, the modeled melt is within 5% of measured values at 9 of the 15 grid cells. In 6 of these cells, modeled melt was within 1%–3% of measured melt, and at one cell there was no difference between the two. The largest error occurred in a cell where the model over-computed melt by 37%. The model over-computed melt in 10 cells and the mean difference was 7% and the median was 3%. In 2009, the average difference at seven stakes was just 1%, but the range was large where the model over computed melt at one cell by 22% and under computed by 28% in another.
In 2009 and 2010, the USGS calculated a glacier-wide summer balance for Sperry Glacier (USGS unpublished data) with an ablation season for both of these years that closely matched the model's time interval of this study. Comparison of modeled glacier-wide meltwater loss to the USGS glacier-wide summer balance suggests model performance improves significantly at larger spatial and temporal scales. In 2009, the modeled July–September meltwater equivalent lost was 4.00 m, while the USGS glacier-wide summer balance was 3.90 m, a difference of 3%. In 2010, the modeled total melt was 2.80 m water equivalent and the USGS balance was 2.95 m, a 5% difference.
The error assessments demonstrate the strengths and limitations of our modeling approach. The largest potential for error is introduced by the temperature and solar radiation indices (±19%–24%), with the model coefficients adding additional potential error (±10%). Thus, daily melt at a single grid cell likely has high associated error, perhaps as much as ±34%. However, as is typical of temperature index models, model performance improves at larger spatial scales (km2) and longer time steps (weeks), with absolute errors down to ±7% and as low as ±3% when cancellation effects are taken into account. For these reasons we believe our model is sufficient to quantify a first-order estimate of glacier-derived runoff across the GNP region.
More runoff from glaciers occurred in August 2009 than in August 2010 (Tables 1 and 5). The glacier-derived runoff from all 39 glaciers was 28.37 × 106 m3 in 2009 and 21.92 × 106 m3 in 2010, with the mean between both years at 25.15 × 106. The mean meltwater volume produced each day was 9.2 × 105 m3 in 2009 and 7.1 × 105 m3 in 2010, a difference of 23%. The highest daily runoff was 1.31 × 106 m3 in 2009, and 1.17 × 106 m3 in 2010. The lowest daily runoff was 3.5 × 105 m3 in 2009, but just 7 × 104 m3 in 2010. The variability among individual glaciers was substantial. The Harrison Glacier in 2009 produced the highest monthly discharge with 3.0 × 106 m3, while the lowest occurred on Gem Glacier in 2010 with 3.0 × 104 m3. (Table 1; Fig. 7).
The August specific discharge (monthly glacier runoff volume divided by glacier area) averaged over all 39 glaciers was 1.68 m in 2009 and 1.33 m in 2010 (Table 1). The standard deviation between glaciers was 0.19 m in 2009 versus 0.17 m in 2010. The maximum specific discharges of individual glaciers were 2.16 m during 2009 and 1.84 m during 2010. Both of these occurred on Lupfer Glacier, the lowest elevation glacier in the study area (mean altitude 1940 m). Minimum values were 1.26 m during 2009 and 1.00 m during 2010, both on Hudson Glacier, which interestingly is not the highest elevation glacier, but the 11th lowest elevation glacier (mean altitude of 2223 m). A weak correlation, R2 values of 0.55 for 2009 and 0.60 for 2010, between mean monthly specific discharge and mean glacier altitude among the 39 study glaciers reveals that specific discharge declined with increasing altitude at a rate of 0.002 m m-1.
Comparison of glaciers with similar elevations but different aspects reveals that glaciers with more southern exposures did not consistently produce more water than glaciers with more northerly ones. For example, Harrison Glacier (mean elevation 2494 m) faces primarily southeast, and the mean monthly specific discharge is 11% greater than Kintla Glacier (mean 2540 m) but 2% lower than Sperry Glacier (mean 2460 m), both of which face predominantly north. A more dominant control on melt appears to be topographic shading. For example, Hudson Glacier occupies a steep cirque with a 200–300 m headwall often keeping it shaded, and produces little melt despite its relatively low elevation.
The North Fork of the Flathead River basin is the largest basin in the study area (3960 km2) and has the lowest percent of glacierized area at 0.12% (Table 4). A stream gage is located about 55 km downstream from the nearest glacier. If all meltwater generated by glaciers within the basin were to flow by the gage, glacier-derived runoff would account for 6.3% of total August discharge in 2009 and 5.0% of the total August discharge in 2010 (Table 5). The next largest catchment, the Middle Fork of the Flathead at 2922 km2 and 0.16% glacierized, has a gage located 33 km downstream from the nearest glacier. Glacier runoff could account for up to 7.5% of total August flows in the Middle Fork during 2009 and 6.3% in 2010.
East of the Continental Divide, the Saint Mary River watershed has a stream gage located 31 km downstream from the closest glacier. This is the third largest basin at 715 km2 and it is 0.70% glacierized. Here, as much as 12.7% of August discharge could originate as glacier melt during a hot and dry summer like 2009, and 8.9% during a cooler summer like 2010. In the Swiftcurrent Creek watershed (80 km2), 1.4% of the catchment was covered by ice and the closest glacier is 9 km upstream from the gage. The potential contribution of glacier runoff was greater, ranging from 27.7% in 2009 to 23.4% in 2010.
On Upper Grinnell Creek, the modeled glacier runoff during August 2009 exceeded the measured total basin runoff, with a relative contribution of 102.6%. The value of 102.6% helps to reveal the magnitude of the likely sources of errors in our methods, which include but are not limited to (1) unaccounted for losses of meltwater such as to the ground or the atmosphere, (2) errors in the melt modeling calculations, (3) an unaccounted for lag-time for glacier meltwater to move through the pro-glacial lake, and (4) unknown precipitation inputs to basin runoff. Records from the Many Glacier SNOTEL indicate 0.04 m of rain fell during August 2009, which was drier than normal. There are no precipitation records from Grinnell Glacier basin during 2009, but assuming the nearby SNOTEL value for the basin and applying it uniformly across the catchment would yield 1.20 × 105 m3 of runoff at the stream gage. This would equate to 7.8% of the August 2009 total discharge. If one accepts this precipitation-runoff statistic and assumes that the remaining discharge (1.41 × 106 m3) came from glacier-melt, then the difference between this non-precipitation-derived runoff and the modeled glacier-melt runoff (1.57 × 106 m3) would be about 11%. Considering these factors, it is plausible that glacier melt likely produced ∼90% of August 2009 runoff.
Percent of glacier-derived runoff within the five gaged watersheds for the month of August.
The percentage of streamflow originating from glacier melt is well approximated by the function Rg = 0.537A + 1.24, where Rg is the log of the percentage of total basin runoff originating from glaciers and A is the log of the percentage of the basin glacierized (Fig. 8). This relationship yielded an R2 = 0.95 and is statistically significant (p < 0.001) using a t-test with 95% confidence intervals.
STRENGTHS AND LIMITATIONS
Our on-glacier weather stations (Fig. 3) were especially effective in GNP. The simple, lightweight, and backpackable design allowed us to access the remote terrain where glaciers exist. Data from these stations created a unique opportunity to calibrate the model specifically to conditions in GNP. Our study demonstrates that a simple, data-driven melt model can still be effective enough to produce estimates of glacier-derived runoff at the first order level. In GNP and similarly remote areas, terms in the water balance equation such as groundwater flow, evapotranspiration, and precipitation are often unknown, and estimates of these would likely come with high uncertainties. Instead of trying to approximate all these terms and back-calculate glacier-derived runoff from stream discharge data, we argue that it was more advantageous to quantify glacier melt rates and that this delivered a more reliable result.
The close agreement between model results, ablation stake measures, glacier-wide balances on Sperry Glacier, and August 2009 discharge at the Grinnell Glacier outlet stream support our estimates of glacier-derived runoff. These comparisons, combined with our error analysis help to validate our methods and suggest that they could be applied to other large and inaccessible glacierized regions where glaciological, meteorological, and hydrological data are lacking. We acknowledge the errors in our model and the shortcomings in our methods in that we do not actually account for how glacier meltwater moves through the glacier and hydrologic systems. Since the possible sinks to glacier runoff are not quantified across GNP, the percentages reported above are surely over-estimates of glacier-derived runoff. The results from Upper Grinnell Creek help to demonstrate this. However, it is important to note that these percentages are intended only to establish upper limits, meaning that glacier-derived runoff in GNP is not likely to exceed these values. Knowing this upper limit provides a new and necessary quantitative understanding of how glaciers affect streams sourced in GNP. Ultimately our work presents a baseline against which more robust future work can be compared.
GLACIERS AND STREAMFLOW
Our results allow for several simple comparisons among the five gaged basins. The impact of glaciers on August stream discharge and temperature diminishes as the upstream area becomes larger and progressively less glacierized. The role of melting glaciers is clearly substantial in the high-elevation streams close to the Continental Divide and the glaciers (Fig. 7), such as in Upper Grinnell Creek. In this stream, the majority of discharge is likely sourced and kept cold by glacier melt through the typically dry and hot month of August.
Another example that illustrates the importance of glacier-melt to Upper Grinnell Creek discharges can be drawn from a comparison between stream gage data and glacier extent measurements from the 1960s and current data. Glaciers covered almost 45% of the catchment area during the 1960s (Johnson, 1980) compared to the current measure of 28.5%. The mean August stream discharge at the gage decreased by 31% from 2.27 × 10-6 m3 mo-1 during the 1960–1969 interval to 1.57 × 10-6 m3 mo-1 during the 2004–2009 interval. While changes in groundwater flow and runoff sourced from precipitation cannot be ruled out as partial causes of the reduction, it is more likely that glacier recession is responsible for the majority of the streamflow reduction. Precipitation data (Johnson, 1980) from nearby Many Glacier show that August rainfall averaged about 0.01 m greater during the 1960s than during 2004–2009. However, the fact that it rained more in the 1960s does not fully account for the greater stream discharges during this time. If one was to assume average August precipitation during the 2004–2009 period was equal to that in the 1960s, average August stream discharge from 2004–2009 would still only increase to 1.60 × 106 m3 mo-1. This is still 29% less than average flows measured in the 1960s.
The streamflow reductions found on Upper Grinnell Creek are consistent with another large-scale study in neighboring British Columbia where discharges in glacier-fed rivers declined as glacierized areas decreased (Stahl and Moore, 2006). Moore and Demuth (2001) also found a reduction in glacier outlet streamflow corresponding to recession of the Place Glacier in British Columbia. These results imply that glacier recession will likely have the greatest impact on high-elevation, headwater streams draining GNP basins with a percentage of glacierized area exceeding 20%. However, in streams draining catchments where glacierized area is much smaller, such as in Swiftcurrent Creek, our results indicate glacier-melt is still an important source of stream discharge during August. Given our estimates, it appears streamflows would still decline due to glacier recession in catchments with percentages of glacierized area numbering in the single digits.
As one moves downstream away from the Continental Divide and the glaciers, the upstream watershed area increases rapidly and the percentage of the catchment covered with glaciers diminishes (Table 4). Thus, the glacier-derived component to discharge in the low elevation rivers, such as the North and Middle Forks of the Flathead, is much less. In these rivers, the coefficient of variation of August discharge is about five times greater than the percent of discharge that could potentially be glacier-derived (Tables 4 and 5). This indicates that other controls on the basin water balance are more dominant than glaciers and that a loss of glacier runoff would not likely have a major effect on August flows. For example, these rivers would not necessarily experience record low flows only because glaciers have disappeared.
The effect of glacier runoff on water temperatures in the North Fork of the Flathead River can be examined by setting up a basic calorimetry problem using the mean August water temperature and discharges measured at the gage and the modeled glacier runoff. We calculate that even if a 5% glacier melt component was added to the site of the North Fork gage at 0 °C (i.e., no warming occurred in transit), the total river discharge would be cooled by no more than 0.8 °C. Water temperature measurements on two different glacier-fed North Fork tributary streams show average August water temperatures had already warmed to 7–8 °C in as little as 4–6 km from the termini of glaciers (D'Angelo and Muhlfeld, 2013). Furthermore, measurements from this same study, taken in four glacier-fed streams near their confluences with the North Fork River, showed average August water temperatures had already reached 16–17 °C. Thus we conclude that glacier runoff has little influence on water temperature in the main body of the North Fork River when assuming all glacier-melt moves through the catchment as surface water.
Other research (Saar and Manga, 2003; Earman et al., 2006; MacDonald et al., 2014) suggests snowmelt-derived water is a significant source of groundwater recharge in mountainous regions. It is possible or even likely that glacier-melt in GNP could indirectly affect stream volumes and temperatures through the groundwater system. Indeed, the interactions between groundwater and surface water in the North and Middle Forks of the Flathead are complex (Hauer et al., 2007), and more research is needed to understand the effects that glacier meltwater has on these alluvial-bed rivers.
One important relationship between glaciers and humans may exist in the Saint Mary River. Here much of the Saint Mary River discharge is delivered into the Milk River via a canal during the months of April through October. This water is used to irrigate approximately 570 km2 of cropland in central Montana (USBR, 2012). From 2000 to 2010, the Saint Mary Canal removed on average 69% of the Saint Mary River discharge during the month of August. While glaciers do not appear to be a dominant source of water in the Saint Mary at the gage, which is 2 km upstream from the canal intake, the estimated percentage of water sourced from glaciers is not trivial considering the irrigation demands on this river.
Our work illustrates the varying effects glaciers have on water availability in five different streams sourced in GNP. Our use of an innovative on-glacier weather station allowed us to access five remote glaciers. This in turn supplied a novel data set from which we could calibrate a glacier melt model specifically to local conditions. Based on comparisons to independent measurements of glacier surface melt and glacier outlet stream runoff, we argue that this model is an effective and reliable tool for estimating glacier-derived runoff at the first-order level. We suggest our methods could be applied to other glacierized regions similar to GNP where glaciological, meteorological, and hydrological data are sparse and/or difficult to obtain.
Our model results estimate that the 39 small named glaciers in the GNP region currently contribute on the order of 25 × 106 m3 of water to August stream flows. The percentage of stream discharge potentially derived from glaciers diminished as a function of the percentage of the basin glacierized upstream from the gage (Fig. 8). The potential contribution to August streamflows varies from 5%–6% in basins as little as 0.12% glacierized to 23%–28% in basins 1.4% glacierized (Table 5). Glacier-melt was probably responsible for more than 90% of August 2009 runoff in a basin 28.5% glacierized. A complete loss of the region's glaciers would reduce August flows and likely lead to warmer water temperatures in the high-elevation stream reaches that drain basins >20% glacierized. Even in streams sourced from basins as little as 1.4% glacierized, such as Swiftcurrent Creek, continued glacier recession will likely reduce late-summer stream flows. Agricultural demands combined with diminishing glacier-melt runoff could increase competition for water from the Saint Mary River. The complete hydrologic role glaciers have on the North and Middle Forks of the Flathead remains uncertain due to the complex processes controlling discharge and water temperature in these two rivers. Even so, our estimates suggest that glaciers are not a dominant source of water to these two rivers during August.
The authors would like to thank the Montana Water Center, the American Alpine Club, the Geological Society of America, the Patrick McDonough Foundation, and the Inland Northwest Research Alliance their financial support. We also thank multiple reviewers for their helpful comments, which significantly improved this paper. This work is also a contribution of the Western Mountain Initiative of the U.S. Geological Survey. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Clark would especially like to thank the many volunteers who worked so hard to haul weather station equipment across some of Glacier Park's most rugged terrain.