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 km2 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 km3, or 1.37 km3 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.


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 (1987Carrara ( , 1989 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 thirtynine glaciers, we examined a slightly longer period ending in 2015 (Fagre et al. 2017).

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 2 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.
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 2 to define a glacier Schiefer, Menounos, and Wheate 2008). We also applied a 0.01 km 2 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 2 threshold and are henceforth referred to as "all glaciers."

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.
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 geoor 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 2 . 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:57km 2 was þ0:32% À0:43% . 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 2 . 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: 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). 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 α ¼ 21:346, and β ¼ 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).
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 2 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 2 during the LIA to 13.73 km 2 in 2015 (−69.34 percent). In 2015, the mean area of remaining glaciers (>0.01 km 2 ) within LIA margins for this subset was 0.44 km 2 (s.d. = 0.41), with the largest at 1.66 km 2 . Total estimated volume of the thirty-nine glaciers was 0.43 km 3 in 2015 (Table 3, with uncertainty caveat), a decrease from 1.64 km 3 (−73.78 percent) during the LIA.

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  Chen and Ohmura (1990) used in area-volume scaling. Fountain andBasagic (2011) andFountain, Glenn, andBasagic (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.  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.
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(Basagic and Fountain 2011and 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).
We estimated the number of LIA glaciers that contained >0.01 km 2 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, 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 2 ) have lost between 42.98 percent and 79.97 percent.  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.
although Fountain et al. defined glaciers based on a 0.01 km 2 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 3 (Table 2). Assuming an ice density of 900 kg m 3 , 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 3 ) 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.