Introduction
Mountain glaciers around the world are important contributors to streamflow and therefore act as a water resource for a wide range of uses, including irrigation, drinking water supply and hydroelectricity (Reference ØstremØstrem, 1966; Reference MeierMeier, 1969; Reference Marston, Pochop, Kerr, Varuska and VeryzerMarston and others, 1991; Reference Yao, Wang, Liu, Pu, Shen and LuYao and others, 2004; Reference Bradley, Vuille, Diaz and VergaraBradley and others, 2006; Reference MooreMoore and others, 2009). Mountain glaciers are also sensitive to climate change and are melting rapidly, therefore drastically reducing their volume (Reference Lemke and SolomonLemke and others, 2007) and contributing to sea-level rise (e.g. Reference GardnerGardner and others, 2013). They may also contribute to potential reduction in streamflow availability for water resource use once they have melted completely (e.g. Reference Mark and SeltzerMark and Seltzer, 2003; Reference VuilleVuille and others, 2008; Reference XuXu and others, 2009; Reference Kaser, Grosshauser and MarzeionKaser and others, 2010), particularly for late summer in arid regions. In the contiguous USA, rapid glacial melting linked to climate change has been noted across several mountain ranges, including the North Cascades and within Glacier National Park, Montana (Reference Hall and FagreHall and Fagre, 2003; Reference PeltoPelto, 2006). With projected increases in temperature as a direct consequence of climate change, many mountain glaciers are expected to melt completely by 2100 (Reference Lemke and SolomonLemke and others, 2007; Reference Nesje, Bakke, Dahl, Lie and MatthewsNesje and others, 2008).
In our study region, the Wind River Range of Wyoming, USA, the temperature at an elevation of 4000 m increased by 3.5°C between the mid-1960s and the early 1990s, as determined from d18O evolution on an extracted ice core from Upper Fremont Glacier (southern part of the range) (Reference NaftzNaftz and others, 2002). Another study suggests that Wind River Range glaciers will disappear as early as 2021 (Reference Marston, Pochop, Kerr, Varuska and VeryzerMarston and others, 1991). While that study is over 20 years old, it is the only published study known to us that uses surface elevation change along with in situ measurements of ice depth to estimate the demise of glaciers in the Wind River Range.
In contrast to the projected rapid melting of many mountain glaciers, small cirque glaciers are often protected by large headwalls that allow shadows to protect against melting, and differential snow accumulation through drifting (Reference KuhnKuhn, 1995). In addition, climatic conditions such as snowfall are not necessarily homogeneous in mountainous regions, leading glaciers to respond differently from one basin to the next (Reference Brown, Harper and HumphreyBrown and others, 2010). Owing to these particular conditions some smaller mountain glaciers may not disappear as quickly as previously anticipated. For example, glaciers in Glacier National Park are projected to melt completely by ~ 2030, based on historical retreat and climate analysis (Reference Hall and FagreHall and Fagre, 2003). However, Sperry Glacier was recently estimated to remain through 2080, as determined by modeling changes in mass balance under various climate change scenarios (Reference Brown, Harper and HumphreyBrown and others, 2010).
Regardless of how rapidly glaciers are melting, glacial meltwater remains a constant water resource for many populations surrounding glaciated mountain ranges. Measurements of glacier melt rates and ice volume are the two key variables needed to calculate estimates of how long a glacier will survive. To calculate ice volume, it is necessary to determine the depth of the ice and create an up-to-date ice surface elevation grid, which allows for creation of a three-dimensional (3-D) model of the glacier. Ice-penetrating radar is a widely established technique for measuring ice thickness of mountain glaciers, ice caps or ice sheets (e.g. Reference Bogorodsky, Bentley and GudmandsenBogorodsky and others, 1985; Reference Nolan, Motyka, Echelmeyer and TrabantNolan and others, 1995; Reference Matsuoka, Maeno, Uratsuka, Fujita, Furukawa and WatanabeMatsuoka and others, 2002; Reference Bauder, Funk and GudmundssonBauder and others, 2003; Reference AzamAzam and others, 2012). While airborne ice-penetrating radar is widely used on ice sheets (Reference GogineniGogineni and others, 2001), ice caps (Reference Dowdeswell, Benham, Gorman, Burgess and SharpDowdeswell and others, 2004) and large mountain glaciers (Reference Conway, Smith, Vaswani, Matsuoka, Rignot and ClausConway and others, 2009), it remains a field-based method when applied to small glaciers in rugged terrain where steep rock walls and narrow valleys cause multiple reflections using airborne radar systems.
Recent technological improvements have made ground-based ice-penetrating radar systems lighter (Reference Mingo and FlowersMingo and Flowers, 2010), making it easier to conduct mass-balance and ice depth measurements on remote glaciers. However, field research is still time-consuming and difficult because of environmental challenges. In some cases, wilderness restrictions that prohibit transportation by helicopter, such as those within US Forest Service wilderness areas (Reference ZahniserZahniser, 1964), add to the difficulty in accessing glaciers. Researchers have to face these challenges to study the glaciers of the Wind River Range, where only a handful of field studies on mass balance and/or ice depth have been published over the past few decades (e.g. Reference Marston, Pochop, Kerr, Varuska and VeryzerMarston and others, 1991; Reference Naftz and SmithNaftz and Smith, 1993; Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). Despite these challenges, we present a study focusing on the mass balance, ice volume and potential remaining lifespan of Continental Glacier in the Wind River Range.
Study Area and Previous Research
Located in west-central Wyoming, the Wind River Range contains the largest concentration of glacial mass (~ 63 glaciers) in the Rocky Mountains of the USA (Reference Cheesbrough, Edmunds, Tootle, Kerr and PochopCheesbrough and others, 2009). This range is oriented northwest to southeast, with the Continental Divide running the length of the range (Fig. 1). The region is characterized as a dry continental climate, with orographically lifted Pacific air masses providing ~ 100–130 cm of precipitation annually, and an average snowfall of 500 cm at the highest elevation of the range (Reference Pochop, Marston, Kerr and VaruskaPochop and others, 1989). Field observations of air temperature on Continental Glacier during August 2012 ranged between 1°C at night and 158C during the day, with maximum extremes of 20°C. Snowmelt onset in the range has occurred earlier in the season since the mid-20th century (Reference Stewart, Cayan and DettingerStewart and others, 2005; Reference Lundquist, Stewart, Dettinger, Cayan and WagnerLundquist and others, 2009) and the snowpack has generally disappeared below 3000 m by the end of June during the 2000s (Reference Hall, Foster, DiGirolamo and RiggsHall and others, 2012).
Fig. 1. Map of Continental Glacier with a hillshade background (inset shows location of Continental Glacier in relation to other Wind River glaciers). High-accuracy (<10 cm) GPS point locations are represented by white circles. Ice radar transects are shown as black lines. Thin black lines represent trend lines parallel to glacier flow along which interpolated points were created for the generation of the surface DEM (based on the GPS points) and bedrock DEM (based on the ice depth). Thin dashed line represents snowline altitude during fieldwork in August 2012. Terminal debris for the upper section is indicated by the speckled pattern to the east of the glacier. Elevation contours are thin gray lines (contour interval of 20 m with index contours labeled every 100 m). Glacier outline (digitized from 2012 US National Agriculture Imagery Program (NAIP) imagery) is represented by a thick black line.
The majority of glaciers are situated near the summits of the range, at elevations between 3500 and 4200 m, with 77% located on the east side of the Continental Divide (Reference Marston, Pochop, Kerr, Varuska and VeryzerMarston and others, 1991). Consequently, most of the glacial meltwater drains into the Wind River, which flows into the Missouri–Mississippi River system (Fig. 1). The Wind River Indian Reservation is also situated directly to the east of the Wind River Range and its population commonly uses snow and glacial meltwater for a variety of purposes, including crop irrigation and supply to fisheries for hunting and fishing (Reference Flanagan and LaituriFlanagan and Laituri, 2004).
While the existence of glaciers in the Wind River Range has been documented by European pioneers since the late 19th century (Reference HaydenHayden, 1878; Reference Wentworth and DeloWentworth and Delo, 1931), field research on the glaciers did not begin until 1950 (Reference MeierMeier, 1951). Since this initial work, only a few field measurements on a few glaciers have been published, mostly due to the challenges posed by the rugged conditions and remote access (Reference Marston, Pochop, Kerr, Varuska and VeryzerMarston and others, 1991; Reference Naftz and SmithNaftz and Smith, 1993; Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). Most Wind River Glacier studies conducted over the last decade have been based either on remote sensing (Reference Cheesbrough, Edmunds, Tootle, Kerr and PochopCheesbrough and others, 2009; Reference Thompson, Tootle, Kerr, Sivanpillai and PochopThompson and others, 2011) or hydrology, using streamflow and isotope data collected downstream of the glaciers (Reference Cable, Ogle and WilliamsCable and others, 2011).
Continental Glacier (43°19′48′ N, 109°41′27′ W;~ 3800 m a.s.l.) is located on the northern end of the Wind River Range. Glacier meltwater drains into West Torrey Creek (a tributary of Torrey Creek) and eventually into the Wind River near Dubois, WY (Fig. 1). The glacier is oriented to the east-northeast and has a north–south dimension of 4 km and an east–west dimension of 800 m (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). The slope of the glacier is relatively low (~ 7–14°), with the exception of a steep terminus (408) at the north end of the glacier. There are a number of mountain peaks along the western head of the glacier, but there are no steep headwalls that could cast shadows onto the glacier. Although topographic maps list Continental Glacier as one complete glacier, it is actually divided into two sections: a larger section of 1.76 km2 at higher elevations between 3700 and 3985 m and a smaller section of 0.74 km2 at lower elevations between 3400 and 3785 m. The two sections appear to respond differently to the environment (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013).
While the glacier appears to have little to no ice flow (as indicated by the lack of movement of debris or crevasses observed on aerial imagery since 1956), using only topo-graphic evidence the glacier would appear to flow from south to north, where the lower section flows into a steep canyon. However, further inspection shows that there are large terminal debris deposits directly to the east of the upper section, indicating a general flow direction from west-southwest to east-northeast (Fig. 1). This observation is consistent with flowlines (thin debris lines exposed by melting that show the deformation of ice through past movement), which are visible in the bare ice on aerial imagery and in the field. Therefore, the glacier is actually wider than it is long. The main outlet stream of the upper section is located at the northeast portion of the glacier (Fig. 1). From this location, the stream flows across the eastern edge of the lower section, where it drains into West Torrey Creek.
Previous research on Continental Glacier was conducted using remote-sensing techniques to determine surface elevation and volume changes (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). This was achieved by differencing two digital elevation models (DEMs): a DEM derived from a United States Geological Survey (USGS) topographic map created from aerial photos in 1966, and an Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)-derived DEM from 2006. Results from this DEM differencing showed that the upper section of the glacier was thinning (0.51 ± 0.19 m a–1), whereas the lower section was essentially stagnant (+0.06 ± 0.19 m a–1) (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). In addition, the annual snowmelt pattern of the upper section observed on satellite images from 1972, 1995, 2004 and 2011 appears to begin in the center of the glacier and progress downward to lower elevations. In conjunction, the areas of greatest thinning, as noted by the 1966–2006 DEM differencing, occur within these areas of early ablation onset (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). During August 2011, fieldwork was conducted on Continental Glacier, but only a few GPS surface elevations and a single ice thickness measurement were obtained.
Owing to the lack of field data from 2011, a second year of fieldwork was conducted to gather more extensive glacier surface elevation and ice depth data to better constrain the current ice volume of Continental Glacier. In August 2012, an extensive GPS survey was conducted to produce a complete up-to-date DEM of the glacier surface. A survey of ice depth transects was also conducted and combined with the surface DEM to create a 3-D model of the glacier. This 3-D model was used in conjunction with calculated surface elevation change rates to determine the remaining lifespan of the glacier. As a result of the larger size and more dynamic conditions of the upper section (i.e. more rapid melt rate), as noted by the 1966–2006 DEM differencing results (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013), this research focuses on the upper section of Continental Glacier.
Data and Methodology
To calculate the remaining ice volume of Continental Glacier, two measurements are required: (1) ice depth collected with an ice-penetrating radar system and (2) surface elevation collected with a differential GPS unit. The surface elevation data were required for creation of the 3-D model of glacier volume and to compare with older surface elevations (particularly from a 1966 topographic map) of the glacier, to produce a map of surface elevation change over time. Following collection of the ice-depth and surface elevation datasets, two DEMs (one for the bed and one for the glacier surface) were produced with a grid size similar to that of the older elevation datasets for consistent comparison.
2012 radar and GPS data collection
Fieldwork was conducted on Continental Glacier during 10–13 August 2012. At that time ~ 60% of the upper section of the glacier was still covered with snow/firn, while the remaining 40% was bare ice (38%) or exposed rock debris (2%). Snow depth data were not collected in 2012; however, the locations where snow depth data (averaging 2 m) were collected in 2011 were melted to bare ice in August 2012. Therefore we assume the remaining snow depth from the previous winter was minimal at the time of the 2012 fieldwork.
A portable impulse radar system from Blue Systems Integration Ltd (Reference Mingo and FlowersMingo and Flowers, 2010) was used for common-offset surveys to measure ice thickness along continuous transects across the glacier (Fig. 1). The radar system includes a monopulse radar transmitter (Reference Narod and ClarkeNarod and Clarke, 1994), two antennas, and a digitizer that integrates data acquisition and storage, and provides real-time visual display of ice depth. For the antennas, a nominal center frequency of 10 MHz is used for sounding the glacier; assuming the resolution is limited to one-quarter of the wavelength, the vertical resolution is ~ 4 m (Reference Mingo and FlowersMingo and Flowers, 2010). A sampling rate between 3 and 5 s allows depth at point locations to be measured every 3 m while traversing the glacier at low speed. Horizontal locations of each depth were recorded by a low-accuracy GPS unit (accuracy <15 m) incorporated into the receiver.
The radar survey resulted in six continuous transects of the glacier, with efforts to conduct transects both perpendicular and parallel to glacier flow (Fig. 1). The nonuniform coverage was due to the presence of crevasses in the accumulation area, making much of it inaccessible. While the upper section of Continental Glacier was the focus of this study, one transect was also completed across the lower section.
Static surface elevation point data were collected using high-accuracy TrimbleTM GeoXH and TrimbleTM GeoXT handheld GPS units. Both GPS units have a horizontal accuracy <10 cm and a vertical accuracy <1 m after differential corrections are applied using nearby base stations (Dubois and Moose, WY, at a distance of 20 and 85 km, respectively). For the surface elevation data, four parallel traverses were made across the glacier with the two GPS units spaced ~ 250 m apart with points collected every ~ 100 m, yielding a total of 142 points. The GPS points were based on the World Geodetic System 1984 ellipsoidal elevation (WGS84) horizontal datum, with vertical measurements based on the Earth Gravitational Model 1996 (EGM96) geoid. As with the ice-penetrating radar transects, certain portions in the accumulation area of the glacier were avoided because of safety concerns. The GeoXH GPS unit was also used to collect sample points along the ice-penetrating radar transects to adjust the ice radar GPS elevations and estimate their associated accuracy (<1 m).
Historical elevation data
A USGS topographic map (1: 24 000) of Continental Glacier produced from September 1966 aerial photographic stereo-pair images was used to create a DEM for comparison with the 2012 surface elevation data. The DEM was generated from a series of elevation points digitized from the contours (every 50 m) on the topographic map (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). Each digitized contour was then interpolated using an ordinary kriging method (Reference KitanidisKitanidis, 1997). A Gaussian model with 15 lags and a 90 m lag size was the best fit for the semivariogram in the kriging process. The final DEM was produced with a spatial resolution of 1 arcsec (30 m × 22 m) and was readjusted to the WGS84 horizontal datum, with elevation converted from feet to meters (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). Uncertainty within the topographic map is assumed to be within half the vertical distance between contour intervals, which in this case is 6.1 m. The absolute accuracy of the DEM was analyzed by comparing DEM elevations with GPS elevations collected in the field on bedrock around the edge of the glacier during August 2011, which indicated an offset of 2.9 ± 3.4 m, with the 1966 DEM being higher than the GPS elevations (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). Comparing the 1966 DEM with six GPS points collected off-glacier for this study during August 2012, the offset between the vertically unadjusted DEM and the GPS points also indicated that the DEM was higher by 2.8 ± 2.3 m. Owing to the consistency in absolute vertical accuracy analysis between VanLooy and others, (2013) and this study, the 1966 DEM used by VanLooy and others (2013) was adjusted for the 2.9 m offset and used for the current study.
2012 surface DEM generation
The differentially corrected (TrimbleTM) GPS surface elevations collected in 2012 were used in conjunction with edge-of-glacier bedrock elevations determined from the adjusted 1966 DEM to interpolate a surface DEM of the glacier for 2012. A preliminary DEM interpolation of these points showed that the limited number of GPS points collected in the field resulted in a poorly interpolated DEM because the shape of the glacier surface is not covered adequately by these few points. However, the continuous transects from the lower-accuracy radar GPS show that the glacier generally has a parabolic surface shape parallel to the glacier flow direction, which was confirmed visually through field observations (Fig. 2). Therefore, to produce more elevation points and generate a more accurate DEM representing the surface shape of the glacier, parabolic trend lines were fitted along the high-accuracy GPS points parallel to the flow direction of the glacier (Fig. 1). Each trend line was spaced ~ 100 m apart, as designated by the frequency of the GPS point collection. With the edge-of-glacier bedrock elevations as ‘tie points’ for the trend lines, new elevation points were interpolated along the trend lines every ~ 30 m. On the eastern edge of the glacier, the bedrock beyond the terminus debris was used as the ‘tie points’ due to the possible instability of the debris. However, trend-line elevations were not interpolated at the southern end of the glacier due to the surface geometry.
Fig. 2. Parabolic trend lines fitting surface elevations (black triangles) of the adjusted ice radar GPS locations and ice thickness (black squares) from the ice radar transect parallel to glacier flow. Division between the accumulation and ablation areas is an approximation based on the late-summer snowline at the time of the fieldwork. The black line on the photograph (view looking north) shows the approximate location of the transect (also shown in Fig. 1).
At its southern end, the glacier rises to its highest elevation and flattens along a ridgeline. This topography is not well modeled by the parabolic shape used for the rest of the glacier; therefore, it was necessary to use only the collected GPS points for the generation of the DEM at this location. However, a limited number of high-accuracy TrimbleTM GPS points were collected at this location and therefore the lower-accuracy elevations from the radar-embedded GPS were used as supplementary data. To obtain a higher vertical accuracy, the low-accuracy ice radar GPS and high-accuracy TrimbleTM GPS points were compared where overlap occurred. Overall across the glacier, 16 points overlapped (within 3 m horizontally), resulting in an offset of 0.31 ± 2.62 m, and the ice radar GPS was higher in elevation. This offset was subtracted from each radar GPS point and the points were used only where there was a lack of high-accuracy TrimbleTM GPS points. Overall, 18% of the elevation points came from the radar GPS.
Ordinary kriging was used to interpolate the new surface points, resulting in a new 2012 DEM with a spatial resolution of 1 arcsec (30 m × 22 m) and a root-mean-square error (RMSE) of 1.88 m. The 2012 DEM was also compared with the remaining high-accuracy GPS points (e.g. points collected while following the ice radar transects and not included in the DEM interpolation process) as another measure of accuracy. The results of this analysis indicated that the 2012 kriged DEM was interpolated, on average, 2.52 ± 3.31 m too high for the absolute elevation. The 2012 DEM was then adjusted by subtracting this value from the entire DEM.
2012 bedrock DEM generation
The bedrock DEM beneath the ice surface was generated using the same method as for the glacier surface DEM. The same trend-line locations parallel to glacier flow were used in the interpolation process for ice thickness in areas where ice radar measurements were lacking. However, uncertainty in the interpolated ice thickness points will increase along trend lines where only one radar transect is present. Therefore, to prevent introducing more uncertainty in the final bedrock DEM, ice thickness was interpolated only where two or more ice radar transects crossed perpendicular to the trend lines. Mapping bedrock elevations first involved extracting ice thickness at each transect point and then subtracting the ice thickness from the newly created 2012 surface DEM. Radar two-way travel time was converted to ice thickness assuming constant radar electromagnetic waves travel into the ice at 1.68 × 108 m s–1 (Reference Mingo and FlowersMingo and Flowers, 2010). The bedrock elevations were used with the edge-of-glacier elevations for interpolating bedrock elevations in locations where ice radar transects were not available.
While verifying glacier bed shape is difficult because it is not possible to confirm glacier morphology visually, radar profiles of the bed following the transects that run parallel to glacier flow generally indicate a concave (or parabolic) shape, consistent with a typical cirque glacier morphology (Figs 2 and 3). The concave pattern of the glacier bed continued to appear as the other ice thicknesses were plotted along their corresponding trend lines. However, the concave pattern flattens out as it nears the terminus (within 250 m) where ice depth becomes relatively small (~ 20–40 m).
Fig. 3. Parabolic trend lines fitting surface elevations (black triangles) of the adjusted ice radar GPS and ice depth (black squares) along a transect parallel to glacier flow. The two asterisks indicate the surface elevation and depth (on the graph) and location (in Fig. 1) of the August 2011 ice radar measurement.
Bedrock elevations located away from the radar profiles were interpolated using the best-fitting trend line of the known bedrock and edge-of-glacier elevation points. The majority of these were third-order polynomial equations. Interpolated points were spaced 30 m apart along the flowlines and each flowline was separated by ~ 100 m. The interpolated bedrock points along with known bedrock elevations determined from the original radar profile depths were used in ordinary kriging (RMSE 4.38 m) to produce a final bedrock DEM with a spatial resolution matching the two surface DEMs (1 arcsec; 30 m × 22 m). The single ice thickness measurement collected with the stationary radar system in August 2011 allowed for a comparison with the interpolated thicknesses from the 2012 data. The August 2011 measurement located near the center of the glacier indicated a depth of 100 ± 4 m, which is in agreement with the interpolated line for the 2012 ice thickness along that particular flowline (Fig. 3). This result confirms that the interpolation method is appropriate for creating the bedrock DEM.
Analytical process and uncertainty analysis
Following the geodetic method (Reference Fountain, Krimmel and TrabantFountain and others, 1997), the 2012 DEM of glacier surface was differenced from the 1966 DEM (using the Raster Calculator Tool in Environmental Systems Research Institute’s (ESRI) ArcGISTM), providing surface elevation change values at each cell in the DEM. These cell values were multiplied by the area of the individual cells (in this case 30 m × 22 m in the Universal Transverse Mercator (UTM) coordinate system) and summed to produce a final volume change between 1966 and 2012. The resulting volume change was divided by the area of the glacier (1.76 km2 as determined from the 1966 topographic map) to calculate the final surface elevation change. Similarly, the remaining volume of the glacier was determined by differencing the 2012 bedrock DEM from the 2012 glacier surface DEM. The total volume of the glacier was then calculated by multiplying the cell size (30 m × 22 m) by the elevation difference results for each cell, and finally summing these results.
There are five sources of uncertainty in the resulting surface elevation change and total volume calculations, all previously described: (1) DEM offset analysis as compared with the GPS points (δa = 3.31 m); (2) RMSE from the interpolation process (δb = 1.88 m); (3) topographic map uncertainty (δc = 6.10 m); (4) vertical uncertainties for the high-accuracy GPS (δ d = 0.10 m); and (5) vertical uncertainties for the radar GPS (δe = 2.62 m). To calculate the total uncertainty (δ q ) for either surface elevation change or total remaining glacier volume, the uncertainty components were summed in quadrature (Reference TaylorTaylor, 1982) using
The surface elevation change uncertainty was then divided by the number of years (46) between the two datasets for an uncertainty value of ± 0.17 m a–1.
Uncertainty in the ice volume includes all of the uncertainty components from the 2012 DEM (i.e. δa, δb , δc, δd and δe ) along with the vertical uncertainty in the ice radar returns (± 4 m). This uncertainty calculation, which was also summed in quadrature, is first expressed in meters (as all components of potential uncertainty are deduced from vertical offsets), but was then converted to a volume uncertainty by multiplying by the area of the glacier (1.76 km2). The result has a total uncertainty in ice volume of ± 10.8 ± 106 m3.
Results and Discussion
Surface-Elevation Changes
The difference calculation between the 1966 and 2012 DEMs resulted in an average surface lowering of 13.8 ± 7.8 m (0.30 ± 0.17 m a–1) for the upper section of Continental Glacier, with rates varying across the glacier between +0.32 and –0.98 m a–1. These estimates of glacier melting are in contrast to those for other Wind River Range glaciers, which are thinning at faster rates (i.e. Dinwoody: 0.88 m w.e. a–1 between 1958 and 1983; Fremont: 0.93 m a–1 in the early 1990s) (Reference Marston, Pochop, Kerr, Varuska and VeryzerMarston and others, 1991; Reference Naftz and SmithNaftz and Smith, 1993). This emphasizes the point that glacial changes can vary from basin to basin, even within the same mountain range (Reference Brown, Harper and HumphreyBrown and others, 2010).
The greatest glacial thinning on the upper section of Continental Glacier (~ 1 m a–1) occurred in three patches located on the central and southern portions of the glacier (Fig. 4). This thinning pattern is interesting as the southern thinning patches occur at high elevations on the glacier. Previous analysis of aerial photographs and satellite imagery over several decades shows that summer ablation consistently begins at three locations near the center and southern end of the glacier, as indicated by visible bare ice versus snow cover on the imagery (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013). This study confirms the presence of the thinning patches and their relationship with the locations of annual ablation onset across the glacier (Fig. 4).
Fig. 4. Surface elevation change of the upper section of Continental Glacier between 1966 and 2012. Lines indicate the location of ablation onset for 1956 (thick white line), 1972 (thick white dashed line), 1995 (thick gray line), 2004 (thick black dashed line) and 2011 (thin black line).
A surface elevation change calculation was carried out between the 2006 ASTER DEM and the 2012 DEM. This was conducted using the same method as for the 1966 and 2012 DEM comparison; however, it was discovered through the resulting difference image that the 2006 ASTER DEM had significant elevation errors in the southern thinning patches. While the original 1966–2006 DEM difference showed plausible thinning rates in the southern thinning patches of 2.9 m a–1 at the extreme (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013), the surface-elevation changes between 2006 and 2012 at the same locations indicated substantial thickening of up to 12 m a–1. This indicated that the 2006 ASTER-derived DEM was the source of the error as these locations were represented by elevations that were too low compared with the 2012 high-accuracy GPS elevations collected in the field, and it is highly unlikely that the glacier thickened by 12 m a–1 between 2006 and 2012. Therefore, the 2006 DEM was discarded in the surface elevation change analysis in this study, and we emphasize the need for caution when using ASTER-derived DEMs in glacier surface change calculations (Reference Nuth and Ka¨a¨bNuth and Kӓӓb, 2011). For these reasons, the current study provides more reliable average glacier thinning rates for the upper section of Continental Glacier for 1966– 2012 (0.30 ± 0.17 m a–1) than those for 1966–2006 (0.51 ± 0.19 m a–1) reported by VanLooy and others (2013).
Glacier area change analysis
Aerial photographs taken by the United States Department of Agriculture (USDA) for the National Agriculture Imagery Program (NAIP) along with Landsat multispectral scanner (MSS) and thematic mapper (TM) imagery were obtained for an analysis of total glacier area change. The NAIP imagery is typically collected near the end of a hydrological year (end of summer) when most of the previous winter’s snow is expected to have melted from the glacier. At this time, the equilibrium line of the glacier is theoretically equivalent to the snowline. Also, the Landsat images for August and September provided a more accurate representation of glacier area, unimpeded by snow cover. However, for six of the ten years of imagery (aerial photograph: 15 September 1956; MSS image: 6 August 1972; TM image: 20 August 1995; TM image: 4 August 2004; aerial photograph: July 2009; and ground photograph: 6 August 2011) snow from the previous winter still covered much of the glacier. The only areas during these years in which bare ice was visible were the same thinning patch locations noted in the surface elevation change difference image (Fig. 4). Two of the years (aerial photographs: mid-September 2006 and late July 2012) showed substantial melting of the previous winter’s snowpack as approximately half of the glacier was exposed as bare ice. The imagery from the other two years (aerial photographs: September 1980 and 1983) showed fresh snowfall, making it difficult to distinguish the extent of the summer ablation areas.
An analysis of glacier area changes was conducted by digitizing the glacier outline over time (1956–2012) from the aerial photographs (1 m resolution) and Landsat imagery (30 m resolution). Minimal changes occurred around most of the glacier’s edge (reduced in size by <3% over the 56 year period), and most of the change occurred near the highest point of the glacier at the southern end where topography is broad and flat. Analysis of rock debris at the terminus and along other portions of the glacier shows negligible change. The only consistent change in debris was along a linear pattern stretching from the northeast portion of the glacier. This area of debris appears to grow substantially (extending 465 m) between 1994 and 2012 (Fig. 5). However, the debris area (including visible large boulders) extends but preserves its original shape over this time period, suggesting that it is being exposed as the glacier thins. The lack of substantial change in glacier area along with this ‘growing’ debris line (despite the lack of boulder or other debris movement) suggests that Continental Glacier is thinning vertically without shrinking horizontally, acting more as a remnant ice field or stagnant glacier.
Fig. 5. NAIP imagery of the exposed debris evolution caused by progressive glacier melting along the northeast portion of the glacier: (a) 1994, (b) 2006 and (c) 2012.
Glacier bed morphology and total remaining ice volume
The bedrock beneath the glacier is composed of two cirque depressions with a steep headwall located generally to the east side of the glacier (Fig. 6a). A south-to-north transect of the glacier bed shows that the bedrock has a profile similar to a cyclopean staircase in which the two cirques step down from higher to lower elevations (Fig. 6c). This profile fits the pattern of other previously glaciated cirques surrounding Continental Glacier and throughout the Wind River Range, visible on topographic maps.
Fig. 6. (a, b) Map of (a) bedrock elevation and (b) ice depth from glacier surface. (c) Transect profile highlighting cyclopean staircase and glacier surface. Estimated tarn lake areas are noted by dashed lines in (a).
The glacier thickness follows a similar pattern to the shape of the base of the glacier (e.g. two deep spots in proximity to the two cirques) (Fig. 6b). The greatest depth of the glacier was ~ 184 m in the northern cirque and ~ 150 m in the southern cirque. Beyond the two cirques, the glacier is relatively thin: 33% of the glacier is <20 m deep, located primarily around the periphery, and 26% is between 20 and 40 m deep. The remaining 41% of the glacier includes the two cirques between 40 and 184 m. An analysis was conducted on the bedrock DEM to determine the possible formation of tarn lakes if the glacier were to melt completely and the two cirques were to fill with water (Fig. 6a). In this event, two lakes would form, covering areas of 0.16 km2 (for the northern cirque) and 0.03 km2 (for the southern cirque). While the lakes would be only a fraction of the glacier in terms of potential water storage, they would still act as reservoirs for snowmelt and rainfall for downstream water use.
The total ice volume of the upper section of Continental Glacier in August 2012 is 72.1 × 106 ± 10.8 × 106 m3, representing ~ 64.9 × 106 ± 9.7 × 106 m3 meltwater for potential water resource use. Given the known volume of the glacier and the surface elevation lowering rates between 1966 and 2012, an estimation of the expected remaining lifespan of the upper section of Continental Glacier was calculated. The glacier volume at each cell (of the DEM difference) was divided by the average surface elevation change rates for the corresponding cells across the entire glacier. This resulted in a glacier map of the number of years until the ice thickness reaches zero (Fig. 7). Given the amount of glacial ice and the relatively slow thinning rate, the glacier is expected to reduce in volume by 24% over the next 50 years, using this simple thinning rate model. Over the next 100 years, the glacier is estimated to reduce in volume by 43% and will bifurcate into multiple ice masses. The remaining glacial mass is estimated to disappear in ~ 300–400 years.
Fig. 7. Lifespan estimation for Continental Glacier. Areas with missing values were due to either positive surface elevation change values (only negative values were used for determining remaining glacial lifespan) or errors within the bedrock and/or surface DEMs leading to inverted ice depth, particularly in locations where the ice thickness was small (near zero), such as at the southern end of the glacier.
There are a number of assumptions with this model. First, only the surface elevation change values indicating thinning were used, as any positive surface elevation changes would indicate unlikely continued accumulation. However, thickening composed only 10% of the glacier and is present around the periphery where the majority of elevation changes between the 1966 and 2012 DEMs are small and within the uncertainty. Second, these estimates do not take into account dynamic changes in melting rates as the glacier shrinks. Finally, the estimates do not take into account any projected climatic change over the next 400 years. Instead, the surface elevation change rate (directly related to climatic conditions) used for this analysis is the average rate measured between 1966 and 2012. It is possible that a more recent surface elevation change rate (e.g. between 2006 and 2012) could lead to different results (likely an acceleration of the melting rate). However, even if the melting rates of the glacier were to triple, our estimates indicate it would still take 100–134 years before the glacier melted completely.
This study also highlights the difficulties in making regional glacier lifespan estimates for the entire mountain range, since differential melting is observed both within an individual glacier and between different glaciers. For example, Reference Marston, Pochop, Kerr, Varuska and VeryzerMarston and others (1991) found that Dinwoody Glacier had 69 × 106 m3 w.e. remaining, and that it was melting at a rate of ~ 0.88 m w.e. a–1. Given these values (along with the area of the glacier, 2.91 km2), a direct calculation of the remaining lifespan of the glacier, assuming a constant melt rate, is 27 years. Applying the same method to Continental Glacier, it will take 123 years to melt completely (using the 46 year average melt rate of 0.30 m a–1). This represents even less time if our projected melt rates were tripled using differential melt rates across the glacier (134 years). However, the lifespan results using either method (averaged melt rates versus variable melt rates) for Continental Glacier are substantially longer than the extrapolated projection by Reference Marston, Pochop, Kerr, Varuska and VeryzerMarston and others (1991) of 30 years for all Wind River Glaciers. This study emphasizes that glacier mass-balance variability in relation to ice thickness at coincident locations needs to be taken into account when attempting to project remaining glacier lifespan.
2012 observations of the lower section
While most of this study focuses on the upper section of Continental Glacier, due to its larger size and dynamic characteristics compared with the lower section analyzed by VanLooy and others (2013), some data were collected and a few observations were made in 2012 concerning the lower section. One ice radar transect and high-accuracy GPS survey was conducted across the lower section (Fig. 1). The ice thickness measurements indicate that the lower section is generally shallow, with an average thickness of 17 m, ranging between 0 and 30 m. Considering this average depth over the area of the lower section (0.74 km2) the amount of remaining glacial mass for this section is estimated to be 12.6 × 106 ± 10.8 × 106 m3 (11.3 × 106 ± 9.7 × 106 m3 w.e.). However, it should be noted that this is a first-order estimate of remaining glacial volume as a lack of ice depth data did not allow for detailed analysis.
Although only six GPS points were collected across the lower section, these points were differenced with the 1966 DEM to determine an average elevation change for the point locations. The elevations between 1966 and 2012 showed no significant change as these locations on the lower section of the glacier actually thickened by 0.06 ± 0.19 m a–1.
Visual observations showed that the only portion of the lower section that had bare ice was located at the terminus, where the topography is steeper (≥ 40°) at the lowest elevations of the glacier (~ 3370 m). The remainder of the lower section (91%) is at elevations between 3600 and 3800 m, with similar slopes as the upper section (8–14°), and was still snow-covered above 3540 m in August 2012. This is in contrast to the upper section, which had ~ 38% bare ice. The contrast in surface elevation change and snow-cover characteristics between the upper and lower sections raises questions because these two sections are changing at different rates despite their proximity, similar slopes and same general aspect and the fact that they are located within the same basin.
Conclusions
Despite being represented as one glacier, Continental Glacier should actually be viewed as two separate glaciers responding differently to the environment. Owing to its substantially larger size and significant changes noted in a previous study (Reference Van Looy, Forster, Barta and TurrinVanLooy and others, 2013), the upper section of Continental Glacier was chosen for detailed analysis of surface elevation changes and remaining glacial volume in order to estimate the remaining time it will exist as a water resource. We observed 72.1 × 106 ± 10.8 × 106 m3 remaining glacial volume (~ 64.9 × 106 ± 9.7 × 106 m3 w.e.). The glacier is thinning at an average rate of 0.30 ± 0.17 m a–1 (24.3 × 106 m3 loss of glacial volume since 1966), with surface changes ranging from –0.98 to +0.32 m a–1. In addition, the glacier area reduced by <3% from 1956 to 2012. This observation highlights the importance of combining surface elevation changes with area changes to represent total glacier changes. If the rates of surface-elevation change on Continental Glacier were to persist, it is estimated that 43% of the glacier would melt out over the next 100 years and the glacier would disappear completely within 300–400 years. Although these estimates are based on linear average thinning between 1966 and 2012 and do not take into account future climate projections, even if glacial melting rates were to triple it is estimated that it would take >130 years for the glacier to melt completely.
Acknowledgements
This work was partially supported by NASA grant NNX10AH20G and by a Faculty Seed Money Award from the Office of the Vice President for Research at the University of North akota. Clément Miège is supported by a NASA ESS fellowship. We thank Sergey Molodtsov, Tatiana Molodtsova and David Buttz for assistance in the field. We also thank the National Forest Service for permission to conduct the fieldwork within the Bridger and Fitzpatrick Wilderness Areas, as well as Blucher Creek Outfitters and Derrick Sellergren for support and assistance in traveling to the field site. Finally, we thank two reviewers for comments which improved the quality of the manuscript.
