Glacier recession since the Little Ice Age: Implications for water storage in a Rocky Mountain landscape

ABSTRACT

Glaciers have significant influence on hydrology, vegetation, and wildlife in mountainous regions, and are receding globally. To quantify the impacts of sustained glacier loss, we mapped a complete set of glacier areas from the Little Ice Age (LIA) using very high-resolution satellite imagery (30 cm) within Glacier National Park (GNP), Montana, a region that encompasses 4098.81 km² in the northwestern United States. We measured glacier change across the park using LIA glacier area as a baseline and then estimated change in glacier area and volume over time. An estimated 146 glaciers existed within the current boundaries of GNP during the LIA. By 2005, only fifty-one (35 percent) persisted. Nearly 90 percent of LIA glaciers had lost more than 50 percent of their area by 2005. This decrease in glacier area equates to an estimated ice volume loss of 1.52 km³, or 1.37 km³ of water storage, roughly equivalent to 90 percent of Lake McDonald, the largest lake in the park. Understanding rates of deglaciation and implications for water storage and use can assist local resource managers and downstream communities in planning for change.

KEYWORDS:

• Glacier
• moraine
• recession
• Little Ice Age
• climate change
• imagery

Introduction

Neoglaciation in the northern Rocky Mountains, USA, began 5,300–7,000 years BP (Carrara 1989; Munroe et al. 2012) and reached its peak at the end of the Little Ice Age (LIA). The LIA was an anomalous cold period that lasted from approximately 1300 to 1850 in North America, and was likely triggered by volcanism (Miller et al. 2012), leading to cold summers and ice growth. Northern Hemisphere cold summers were maintained by sea-ice/ocean feedbacks, aided by an insolation minimum. During the LIA, fluvial and drought periods, as determined by tree-ring analyses, alternated in multi-decadal cyclicity (Pederson et al. 2004). A sustained fluvial episode of ~50 years occurred near the end of the LIA that increased alpine glacier mass more than previous periods. Glacial advance was robust and left distinct and measurable moraines that indicated the maximum LIA glacier extents. Moraines are continuous arcuate rock deposits with distinct crests that form at the ice edge of active glaciers. In Glacier National Park (GNP), Montana, the LIA advance overrode most moraines left by the earlier, smaller glacier advances that occurred since the onset of neoglaciation (Carrara 1989); only a few moraines from earlier advances remain. Since the furthest glacial advance of the LIA, glaciers in GNP have retreated rapidly (Dyson 1948; Carrara and McGimsey 1981; Key, Fagre, and Menicke 2002; Hall and Fagre 2003; Fountain 2007; Fagre et al. 2017). Initial glacier retreat has been attributed to recovery from the LIA, but Marzeion et al. (2014) provide documentation of increasing anthropogenic influence on glacier retreat in recent decades. Carrara and McGimsey (1981) estimated that GNP, established in 1910, had lost more than two-thirds of its glaciers by 1980. During the forty-nine years between 1966 and 2015, Fagre et al. (2017) documented a mean reduction in glacier area of 40 percent. Like many other mountain landscapes in western North America, exploration, documentation and mapping of mountain features, including glaciers, did not occur in GNP until early in the twentieth century. A complete map of the park based on aerial photography was not produced until 1968, underscoring the limited availability of geospatially accurate glacial information prior to that time.

A longer record of the magnitude of glacier retreat helps scale the ongoing loss of ice, and provides important historical context for landscape managers. Such information is critical to understand potential past variations in hydrological dynamics across the region (Fountain and Tangborn 1985; Jansson, Hock, and Schneider 2003), and provides a contrast to current mountain water storage as ice (Fagre et al. 1997; Brittain and Milner 2001; Clark, Harper, and Fagre 2015). Threatened and rare alpine species (e.g., Bull Trout, Salvelinus confluentus; Meltwater Stonefly, Lednia tumana; Ptarmigan, Lagopus muta) and local human communities have depended upon the mountains as the “water towers” of the Northern Rockies for millennia. Understanding glacier contribution to past water supplies provides a basis for assessing future changes in water storage with continued glacier loss.

Modern regional glacier inventories have been developed for several parts of the world, including a recent inventory of the American West (Fountain, Glenn, and Basagic 2017), which includes the Lewis and Livingston Ranges in GNP. However, inventories typically capture glacier extent at one point in time and do not address glacier loss over time, which is important for understanding climate drivers, water balance, and sea level rise. Carrara (1987, 1989) and Carrara and McGimsey (1981, 1988) mapped moraines for some glaciers in the park and also identified more than 150 potential moraines parkwide using aerial photography, but the moraines had not been mapped to modern accuracy standards using GIS and remote sensing technologies. Changes in glacier area have been established for subsets of glaciers over various spans of time in GNP; however, until recently, no comprehensive mapping of glacier area at the end of the LIA existed for GNP that would provide a longer record of change.

In this article, we first analyze a recent inventory of the maximum extent of glaciers in GNP at the end of the LIA that used high-resolution satellite imagery to map moraines (Fagre and Martin-Mikle 2018). We compare this to a glacier inventory from 2005 (Fagre et al. 2017; Fagre, Martin-Mikle, and Clark 2019) to determine the magnitude of ice loss by estimating glacier areas and volumes from the LIA and modern glaciers. For a subset of thirty-nine glaciers, we examined a slightly longer period ending in 2015 (Fagre et al. 2017).

Methods

Study area

GNP, located in the northwest corner of Montana adjacent to the Canadian border (Figure 1), is dominated by the Lewis and Livingston ranges, which run northwest to southeast. The Continental Divide bisects the park, dividing GNP into eastern and western sides that differ in climate and vegetation. Elevation ranges from 947 m at the confluence of the North Fork and Middle Fork of the Flathead River in west central GNP to 3194 m at Mount Cleveland in northcentral GNP. GNP covers 4098.81 km² of protected land area, and is designated as a World Heritage Site, International Peace Park, and Biosphere Reserve. As a relatively undisturbed natural system, it is a bellwether for climate change (Gonzalez et al. 2018). GNP climate is transitional between colder, drier continental and warmer, wetter maritime climate regimes, and its snow avalanche climatology has been characterized as intermountain (Mock and Birkeland 2000). The intermountain climate regime is most common between 2000 and 3500 m, the elevation zone for GNP glaciers, and is classified by annual snowfall >560 cm and mean winter temperatures >−7°C. Finklin (1986) reported 24 h mean summer temperatures between 7.5–9°C, and mean July/August precipitation of 135–150 mm, in the 2000–3500 m zone along the Continental Divide. Most of the thirty-one modern named glaciers in GNP (including two located in the adjacent Flathead National Forest) are cirque glaciers (e.g., Figure 2), with some described more specifically as niche glaciers (Cogley et al. 2011). No valley glaciers are present.

Figure 1. Glacier National Park is located in northwestern Montana, U.S., and is bordered by Canada to the north. Little Ice Age glaciers (Fagre and Martin-Mikle 2018), including two southwest of the park in the Flathead National Forest, are shown in red and are located in the higher elevations (black) of the park in the Lewis and Livingston mountain ranges. Elevation data source: National Elevation Dataset, 2002; U.S. Geological Survey web site.

Figure 2. Swiftcurrent Glacier, 27 August 2005, National Agricultural Imagery Program. Swiftcurrent Glacier in Glacier National Park is an example of the cirque glaciers that are characteristic of the glaciers present in the park today. The Little Ice Age extent of Swiftcurrent Glacier is denoted by the red line.

We consistently analyzed the dataset using two different size thresholds. Following convention, and to maintain consistency with previously published glacier area data for GNP, we applied a minimum area threshold of 0.1 km² to define a glacier (Paul et al. 2004; Schiefer, Menounos, and Wheate 2008). We also applied a 0.01 km² threshold to the data, the recently suggested international standard for small alpine glaciers (Paul et al. 2009; Fountain, Glenn, and Basagic 2017). Using the smaller threshold added three additional glaciers to the modern glacier inventory by Fagre et al. (2017). In this article, all LIA and modern glacier areas were greater than the 0.01 km² threshold and are henceforth referred to as “all glaciers.”

Data sources

LIA glacier area

The LIA glacier area dataset (Fagre and Martin-Mikle 2018) was created by digitizing all glacier moraines visible in Digital Globe’s pansharpened, multispectral WorldView satellite imagery collected during the years 2014, 2015, and 2016 at the 1:2000 scale. All WorldView imagery had a spatial resolution of 30 cm and was obtained in late summer (July through September) to minimize the presence of seasonal snow. WorldView imagery did not have complete coverage across GNP and so was supplemented with 2009 imagery from the Department of Homeland Security (1 m spatial resolution) and 2005 imagery from the National Agricultural Imagery Program (NAIP) (1 m spatial resolution) (Table 1). Oblique terrestrial and aerial photography, including an aerial collection of photographs of glacier moraines shot by John Scurlock in 2009, was referenced to help delineate moraines in cases where seasonal snow or shadows made image interpretation difficult (Portland State University 2009). Historic photographs, dating back to 1887, were used to constrain moraine position based on visible glacier extent during the late nineteenth and early twentieth centuries (U.S. Geological Survey n.d.). In locations where imagery was poor in quality, we used Google Earth to view archived imagery at oblique perspectives. Together, consulting these complementary imagery sources ensured the most accurate reconstruction of LIA glacier extent possible.

Table 1. Imagery-based data sources used to derive LIA glacier area data.

Data source	Acquisition date	Spatial resolution	Coverage
WorldView Satellite Imagery	2014–2016	30 cm	Mountainous regions of GNP
Border Imagery, Department of Homeland Security (DOHS)	2009	1 m	Northern two-thirds of GNP
National Agricultural Imagery Program (NAIP)	2005	1 m	Entirety of GNP

Potential inaccuracies inherent to WorldView and NAIP imagery were not addressed, except by excluding images that were visually distorted or heavily shaded. WorldView, NAIP, and other products downloaded from web sources were already corrected for spatial accuracy and were used as is; i.e., no additional geo- or ortho-rectification was conducted. Because off-nadir angles are reported for all WorldView images, we were able to select for images with the lowest off-nadir angle possible to minimize spatial error.

Digitized LIA glacier areas based on moraines, by definition, reflect LIA glaciers that were active. These digitized areas do not represent a comprehensive accounting of all permanent snow, ice, or periglacial (i.e., rock glaciers) features in the park that do not leave moraines. Moraines within GNP that pre-date the LIA were mapped by Carrara (1987) and Carrara and McGimsey (1988) and are not included in this analysis. The criteria used for identifying the former size and location of glaciers were (1) the presence of a terminal moraine; (2) the presence of lateral moraines; and (3) whether the identified glacier area exceeded 0.01 km². A potential glacier was not omitted if it failed to meet both of the first two criteria (i.e., a LIA glacier that terminated at a cliff with well-defined lateral moraines but no discernible terminal moraine).

The primary investigator digitized each LIA glacier margin by tracing the cirque headwalls and visible moraines in a GIS. Two reviewers independently examined the imagery and each polygon of the LIA glacier margins for positional accuracy. For 58 percent of the glaciers, the original margin trace was confirmed by the two reviewers. For 29 percent of the LIA glaciers, there were fractions of the margin that were poorly defined or were deeply shadowed (2.87 percent of the total 451.83 km of digitized margin length for all 146 glaciers combined) and reviewers differed by ±1 m in positioning the margin. These differences were further examined utilizing additional imagery and historical photography for clarification and then averaged for the final location of the margin. In <17 percent of the LIA glaciers (1.40 percent of the total 451.83 km of digitized margin length for all 146 glaciers combined), positional differences between reviewers were ±2.5 m for fractions of the margins, primarily where moraine crests were rounded rather than sharply ridged. The area uncertainty due to interpretation differences in margin position for all 146 LIA glaciers combined was ±0.04%. Eighty-five glacier margins were confirmed by reviewers; sixty-one glacier margins, where differences occurred, had ±2.65 percent uncertainty. The mean uncertainty for 146 individual glaciers was ±0.24 percent.

Another source of uncertainty occurred as a result of the digitization process. This uncertainty was estimated by using the ground width of the mapped line to buffer each polygon ±0.5 m. The resulting uncertainty for the total LIA area of $59.57\mathrm{k}{\mathrm{m}}^2\mathrm{w}\mathrm{a}\mathrm{s}{}_{-0.43\mathrm{\%}}^{+0.32\mathrm{\%}}$. The mean uncertainty for individual glaciers was ±0.83 percent with a range of 0.15 percent to 2.45 percent.

Modern glacier area

Modern area data for a subset of thirty-nine GNP glaciers from the years 2005 and 2015 were obtained from Fagre et al. (2017) and are based on NAIP and WorldView imagery, respectively. We used additional data on perennial snow and ice developed from 2005 NAIP imagery (Fagre, Martin-Mikle, and Clark 2019) that included a subset of ice bodies characterized as glaciers by displaying crevassing and banding indicative of active glacial ice and being >0.01 km². We merged these two datasets to have complete coverage of 2005 glaciers for comparison to the 146 LIA glaciers.

Analysis

After digitizing moraines and obtaining areas for all LIA and modern glaciers, volume change was estimated (R Studio 2015) using a simple area-to-volume scaling relationship (Chen and Ohmura 1990; Fountain  et al. 2017).

The scaling relation is expressed as:

$V=\alpha A^{\beta }$

where α, β are constants and can be either empirically or theoretically derived (Chen and Ohmura 1990; Bahr, Meier, and Peckham 1997; Basagic and Fountain 2011). Chen and Ohmura (1990) empirically derived nine sets of regionally applicable parameters, four sets of which include data from alpine glaciers located in the western United States. The four relevant parameter sets are: “Cascade and other areas,” “Cascade, small glaciers,” “Alps, Cascade, and other areas,” and “Alps, Cascade, and Svalbard etc.” We estimated volume using each of the four sets of constants (Figure 3). Basagic and Fountain (2011), Dick (2013), and Fountain, Glenn, and Basagic (2017) found that the number and volume of glaciers in the western United States were best estimated using “Cascades, small glaciers” parameters from Chen and Ohmura (1990), where $\alpha =21.346$, and $\beta =1.145$. While we illustrate differences between the four parameter sets, we focus on results derived from the “Cascade, small glaciers” set of constants. Due to the deviation of individual glaciers from the regression, Granshaw and Fountain (2006) estimate that uncertainty at an individual glacier level can be large (± 50 percent). Because there is only one glacier in our study area with volume data for verification, we provide an estimate of total volume across GNP only and do not include an uncertainty estimate, with the caveat that volume estimates should be regarded as approximate (Basagic and Fountain 2011).

Figure 3. Estimated areas and volumes for 146 Little Ice Age glaciers based on four sets of constants developed by Chen and Ohmura (1990) used in area-volume scaling. Fountain and Basagic (2011) and Fountain, Glenn, and Basagic (2017) found “Cascade, small glaciers” to best fit the western United States and, as such, we employ this set of constants as well. Granshaw and Fountain (2006) estimate that uncertainty can be as much as ±50 percent for individual glaciers, but when calculating total volume on a regional scale, can be smaller as a result of errors offsetting one another. Following Fountain and Basagic (2011), these volumes should be regarded as approximate.

Results

There were 146 LIA glaciers identified in GNP, including two in the adjacent Flathead National Forest (Table 2). Eighty-two LIA glaciers were >0.1 km². Of these larger glaciers, thirty-four remained >0.1 km² in 2005. Ninety-five LIA glaciers (65%) no longer contained any active glacial ice by 2005. Nearly 90 percent (131/146) of LIA glaciers lost >50 percent of their area by 2005 (Figure 4).

Table 2. Little Ice Age to 2005 number, area, and volume estimates of glaciers in GNP, including two in the adjacent Flathead National Forest. In the data presented here, all glacier areas are >0.01 km².

	LIA	2005	Percent decline
Number of all glaciers	146	51	65.07
Number of glaciers ≥0.1 km^2	82	34	58.54
Total area of all glaciers (km^2)	59.57	16.40	72.48
Total area of glaciers ≥0.1 km^2 (km^2)	56.32	15.56	72.37
Total volume of all glaciers (km^3)	2.03	0.51	75.10
Total volume of glaciers ≥0.1 km^2 (km^3)	1.96	0.49	75.12

Note. Granshaw and Fountain (2006) estimate that uncertainty in volume estimates can be as much as ±50 percent for individual glaciers, but when calculating total volume on a regional scale, can be smaller as a result of errors offsetting one another. Following Fountain and Basagic (2011), these volumes should be regarded as approximate.

Figure 4. Percent of LIA glacier area decrease (black) from the LIA to 2005.

LIA glaciers decreased from 59.57 km² to 16.40 km² (−72.48 percent). The mean LIA area was 0.41 km² (s.d. = 0.79) with the largest at 4.97 km² (Blackfoot Glacier). By 2005, the mean area of active glacial ice (> 0.01 km²) within LIA glacier margins was 0.32 km² (s.d. = 0.39), with the largest at 1.67 km² (Harrison Glacier). Ninety-five glaciers (43.17 km²) disappeared (Figure 4). Additionally, the largest glaciers (>1.91 km²) lost 42.98–79.97 percent (Figures 5 and 6). The volume of LIA glaciers was estimated to be 2.03 km³ during the LIA (Table 2, including uncertainty caveat). The estimated volume of these same glaciers in 2005 was 0.51 km³ (−74.88 percent).

Figure 5. Distribution, Little Ice Age (LIA) relative sizes, and percent area decrease from the LIA to 2005 of all glaciers in Glacier National Park (and two southwest of the park in the Flathead National Forest). LIA glaciers with high percent area decrease are shown in yellow, and those with low percent area decrease are shown in dark blue.

Figure 6. The relationship between GNP LIA glacier area and percent decrease from the LIA to 2005 demonstrates that variability in percent decrease exhibited by small glaciers is large, whereas all larger glaciers (>1.91 km²) have lost between 42.98 percent and 79.97 percent.

For a subset of thirty-nine LIA glaciers in GNP, including two in the adjacent Flathead National Forest, glacier area data was available for 2015 (Fagre et al. 2017) (Table 3). Of the thirty-nine LIA glaciers, thirty-one contained glacial ice in 2015. Twenty-eight were >0.1 km² in 2015, including two glaciers located in the adjacent Flathead National Forest. Seventy-seven percent (30/39) of the subset lost >50 percent of their LIA area by 2015. The thirty-nine glaciers decreased from 44.78 km² during the LIA to 13.73 km² in 2015 (−69.34 percent). In 2015, the mean area of remaining glaciers (>0.01 km²) within LIA margins for this subset was 0.44 km² (s.d. = 0.41), with the largest at 1.66 km². Total estimated volume of the thirty-nine glaciers was 0.43 km³ in 2015 (Table 3, with uncertainty caveat), a decrease from 1.64 km³ (−73.78 percent) during the LIA.

Table 3. Little Ice Age, 2005 and 2015 number, area, and volume estimates for subset of thirty-nine glaciers in GNP, including two in the adjacent Flathead National Forest. In the data presented here, all glacier areas are >0.01 km².

	LIA	2005	2015	2005 percent decline	2015 percent decline	2005–2015 percent decline
Number of all glaciers	39	31	31	20.51	20.51	0.00
Number of glaciers ≥0.1 km^2	37	29	28	21.62	24.32	3.45
Total area of all glaciers (km^2)	44.78	14.75	13.73	67.06	69.35	6.92
Total area of glaciers ≥0.1 km^2 (km^2)	44.68	14.67	13.55	67.17	69.66	7.63
Total volume of all glaciers (km^3)	1.64	0.47	0.43	71.34	73.72	8.51
Total volume of glaciers ≥0.1 km^2 (km^3)	1.64	0.46	0.43	71.95	73.92	6.52

Note. Granshaw and Fountain (2006) estimate that uncertainty in volume estimates can be as much as ±50 percent for individual glaciers, but when calculating total volume on a regional scale, can be smaller as a result of errors offsetting one another. Following Fountain and Basagic (2011), these volumes should be regarded as approximate.

Discussion

A conservative estimate of 146 LIA glaciers in what is now GNP (including two in the adjacent Flathead National Forest) is based on very high-resolution satellite imagery analysis of moraines (Fagre and Martin-Mikle 2018). We acknowledge that an unknown number of small LIA glaciers could have existed but did not leave discernable geomorphological evidence. Our estimate is comparable to Carrara (1987), who identified more than 150 potential moraines within GNP. Our dataset is the first to map all discernable moraines in GNP at fine scale (1:2,000) using satellite imagery and digital technologies.

Within GNP and the adjacent Flathead National Forest, we documented a −72.48 percent change in LIA glacier area by 2005. These values are broadly consistent with other regional estimates of glacier area loss that occurred for shorter spans and/or for fewer glaciers within the Lewis and Livingston Ranges, GNP. For instance Key, Fagre, and Menicke (2002) reported a –73 percent change for eleven selected glaciers from end of LIA to 1993; Hall and Fagre (2003), based on work by Carrara and McGimsey (1988), reported a –65 percent change for glaciers in the Jackson-Gunsight area from end of LIA to 1979; and Fountain (2007) described a –66 percent change for glaciers in Northwest Montana from end of LIA to present, the majority of which are in GNP. Glacier loss in GNP is occurring more rapidly than in other regions of the western United States over similar time periods; e.g., Sierra Nevada, CA, −55 percent from 1903–2004 (Basagic and Fountain 2011); and Wind River Range, WY, −47 percent from 1900–2006 (DeVisser and Fountain 2015). While Clark, Harper, and Fagre (2015) demonstrated that most glaciated basins in GNP do not currently contribute significantly to stream flow, our estimate of glacier area and volumetric loss over time suggests that they likely did in the past. The percent change in glacier area from LIA to 2005 ranged from little change at Gem Glacier (−6.75 percent) to −89.24 percent at Thunderbird Glacier (Table 4).

Table 4. Little Ice Age and 2005 glacier areas (m²) and corresponding percent decrease of glacier area for a subset of thirty-nine glaciers in Glacier National Park and in the adjacent Flathead National Forest.

Glacier name	LIA area	2005 area	Percent decline from LIA to 2005
Baby Glacier	136133.57	*	100
Boulder Glacier	829577.99	*	100
North Swiftcurrent Glacier	227306.21	*	100
Herbst Glacier	343244.32	*	100
Harris Glacier	79286.24	*	100
Red Eagle	477070.9	*	100
Shepard Glacier	497244.16	*	100
Hudson Glacier	132164.98	*	100
Thunderbird Glacier	1075097.51	115652	89.24
Whitecrow Glacier	741530.93	110866	85.05
Two Ocean Glacier	1076919.16	187043	82.63
Sperry Glacier	3793322.12	883968	76.7
Dixon Glacier	671001.22	156474	76.68
Agassiz Glacier	4246977.56	1020710	75.97
Jackson Glacier	3094285.19	794333	74.33
Swiftcurrent Glacier	594231.7	175401	70.48
Grinnell Glacier	1976494.97	602469	69.52
Stanton Glacier	950202.24	298733	68.56
Blackfoot Glacier	4970064.26	1608800	67.63
Ipasha Glacier	594065.79	193386	67.45
Kintla Glacier	2780343.02	919674	66.92
Carter Glacier	708091.34	234275	66.91
Chaney Glacier	1078922.83	358892	66.74
Grant Glacier	842943.74	284506	66.25
Lupfer Glacier	171200.25	63826.7	62.72
Miche Wabun	263262.3	104614	60.26
Vulture Glacier	811528.41	329091	59.45
Pumpelly Glacier	1978996.13	898616	54.59
Harrison Glacier	3470484.89	1666770	51.97
Logan Glacier	694335.54	350642	49.5
Siyeh Glacier	378010.17	208857	44.75
Sexton Glacier	528179.2	298084	43.56
Rainbow Glacier	1910417.47	1089260	42.98
Weasel Collar Glacier	809813.09	505863	37.53
Ahern Glacier	765116.39	508211	33.58
Old Sun Glacier	499662.38	340581	31.84
Salamander Glacier	250056.48	172592	30.98
Piegan Glacier	311173.34	247162	20.57
Gem Glacier	24364.74	22719.1	6.75

* Values are below the 0.01 km² threshold considered to be glaciers.

⁺ Previously published 2005 area values may differ slightly from 2005 area values in this table, which represent areas within Little Ice Age boundaries.

We estimated the number of LIA glaciers that contained >0.01 km² of glacial ice in 2005 to be fifty-one. This is in contrast to the Fountain, Glenn, and Basagic (2017) estimate of 151 glaciers in ~1966 in the GNP area. Our two studies are not directly comparable; in addition to the large difference in dates, the Fountain, Glenn, and Basagic (2017) regional-scale analysis was based on 1:24,000-scale USGS maps created from aerial photography collected between 1955–1966, whereas we estimated the number of LIA glaciers with existing glacial ice based on a combination of 30 cm WorldView imagery from 2014–2016 and 1 m 2005 NAIP imagery. Additionally, although Fountain et al. defined glaciers based on a 0.01 km² minimum threshold, they differentiated between perennial snow and ice and glaciers based on a basal shear stress threshold, where we utilized visual indicators of active glacial ice such as banding and crevassing.

Volume estimates derived from area-volume scaling suggest an ice volume loss of about 1.52 km³ (Table 2). Assuming an ice density of 900 kg m³, this translates to a loss of ~ 1.37 Gt of water across the GNP system between the LIA and 2005. This is roughly equivalent to 90 percent of Lake McDonald, the largest lake in the Park. The largest six glaciers (LIA total = 0.903 km³) accounted for approximately the same amount of glaciated volume currently held in the entire state of Montana (Fountain, Glenn, and Basagic 2017).

It is apparent that all glaciers have not retreated at the same rate (Figure 6). Area loss is likely highly dependent on local topographic variables (e.g., elevation, aspect, slope, and area), micro-climates, and shading. While the largest glaciers tend to have lost more area, as glaciers retreat and become smaller, loss rates become extremely variable. The ability to more accurately model future water storage and shifts in glacial meltwater contribution to streamflow would allow researchers and managers to better grasp the implications of continued glacier loss to communities and ecosystems. This study documents a glacier water storage transition for one specific mountain landscape, but similar transitions are likely to become more common on a global basis.

Acknowledgments

We would like to thank Lisa McKeon for assistance in developing the LIA dataset. We would also like to thank Caitlyn Florentine and Andrew Fountain for providing significant improvements to an earlier draft, and we would also like to thank multiple reviewers for their helpful comments, which greatly improved this paper. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the U.S. Geological Survey Land Resources Land Change Science Program.