STUDY GLACIERS

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 km² but areas range from 0.02–1.89 km². Only Sperry Glacier (∼0.87 km²) 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.

FIGURE 2.

Aerial photographs (National Agricultural Image Program imagery from August 2005) of the five instrumented glaciers with letters corresponding to those on Figure 1 (a) Grinnell, (b) Sexton, (c) Sperry, (d) Blackfoot, and (e) Pumpelly. Yellow circles show locations of the on-ice weather stations, red triangles indicate an ablation stake, dashed red line is the Continental Divide.

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.

TABLE 3

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

FIGURE 3.

Cross-sectional view of an on-ice weather station showing the tower structure, and placement of the instruments, solar panel, and data logger. Two horizontal cross bars form a “T” connecting the vertical legs above the glacier's surface; the second horizontal pole is normal to view and is depicted by circle on right side of the middle vertical pole.

TABLE 4

Data for stream gage sites.

MODELING SCHEME

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), T_c 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, R² = 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 km² 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 m², <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.

FIGURE 4.

(A) Time series plot of average daily temperatures during 2010 at the five on-glacier weather stations normalized by the mean elevation of these sites. (B) Normalized hourly temperatures at the five stations during the warmest day on record (8/26/10). (C) Normalized hourly temperatures during the coldest day on record (9/6/10). Dates are labeled as m/dd/yy.

FIGURE 5.

Comparison of daily cumulative melt (surface lowering converted to a water equivalent) measured with the sonic distance ranger at the five on-ice weather stations, June 27 to September 26, 2010. Dates are labeled as m/dd/yy.

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

FIGURE 6.

Sensitivity to melt coefficient values. Time series plot of modeled August glacier runoff for the year 2010 using the chosen coefficients (thick black line) and the 20 other simulations as thin colored lines.

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.

ERROR ASSESSEMENT

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 (km²) 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.

GLACIER RUNOFF

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 × 10⁶ m³ in 2009 and 21.92 × 10⁶ m³ in 2010, with the mean between both years at 25.15 × 10⁶. The mean meltwater volume produced each day was 9.2 × 10⁵ m³ in 2009 and 7.1 × 10⁵ m³ in 2010, a difference of 23%. The highest daily runoff was 1.31 × 10⁶ m³ in 2009, and 1.17 × 10⁶ m³ in 2010. The lowest daily runoff was 3.5 × 10⁵ m³ in 2009, but just 7 × 10⁴ m³ in 2010. The variability among individual glaciers was substantial. The Harrison Glacier in 2009 produced the highest monthly discharge with 3.0 × 10⁶ m³, while the lowest occurred on Gem Glacier in 2010 with 3.0 × 10⁴ m³. (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, R² 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.

FIGURE 7.

Colored circles represent the 2009 and 2010 average modeled volume of August meltwater runoff at each of the 39 named glaciers. The mean value per individual glacier was 6.4 × 10⁵ m³ with a standard deviation of 6.8 × 10⁵ m³, but the median was only 3.8 × 10⁵ m³.

BASIN CONTRIBUTIONS

The North Fork of the Flathead River basin is the largest basin in the study area (3960 km²) 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 km² 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 km² 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 km²), 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 × 10⁵ m³ 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 × 10⁶ m³) came from glacier-melt, then the difference between this non-precipitation-derived runoff and the modeled glacier-melt runoff (1.57 × 10⁶ m³) would be about 11%. Considering these factors, it is plausible that glacier melt likely produced ∼90% of August 2009 runoff.

TABLE 5

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 R_g = 0.537A + 1.24, where R_g 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 R² = 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 m³ mo^-1 during the 1960–1969 interval to 1.57 × 10^-6 m³ 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 × 10⁶ m³ mo^-1. This is still 29% less than average flows measured in the 1960s.

FIGURE 8.

Percent of total basin runoff that is potentially glacier derived (discharge percent) for Augusts 2009 and 2010 versus the percent glacier-covered area (basin percent). Note only August 2009 is included for Upper Grinnell Creek.

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 km² 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.