Accelerated retreat of coastal glaciers in the Western Prince William Sound, Alaska

ABSTRACT

Analyzing historical maps and Landsat imagery indicates that coastal glaciers in the western Prince William Sound (PWS) have retreated since the end of the Little Ice Age, exhibiting accelerated retreat after the mid-2000s. A multitemporal inventory of 43 glaciers was developed using historical field observations, topographic maps, and Landsat imagery. Area and length measurements are derived from digitized outlines, and center lines are derived from a semi-automatic, geographic information systems (GIS)-based algorithm. Land-based glaciers retreated at a peak rate of 48 m a⁻¹ from the mid-2000s to 2018, more than doubling the average rate of retreat (22 m a⁻¹) for the preceding fifty-year period. From ~1950 to 2018, the total area of land-based glaciers decreased by 228 km², with 36 percent of the glacier loss occurring after the mid-2000s. Simple upscaling of area and volume changes to unmeasured glaciers across the entire PWS resulted in an estimated aggregate glacier mass loss of 379 Gt, equivalent to a 1.047 mm rise in sea level from the 1950s to 2018. Tidewater glaciers respond asynchronously with differing periods of peak area and length loss and lower average rate of retreat compared to land-based glaciers. Glacier retreat correlates with increased summer and winter temperatures and decreased winter precipitation.

KEYWORDS:

• Alaska glaciers
• tidewater glaciers
• glacier fluctuation
• climate change
• rising sea level

Introduction

Glaciers play an essential role in the study of climate change, and advancing and retreating glaciers provide reliable indicators for climate change at decadal to centennial scales (Adalgeirsdóttir, Echelmeyer, and Harrison 1998; Wiles, Barclay, and Calkin 1999; Calkin, Wiles, and Barclay 2001; Santos et al. 2010; Qiang, Lei, and Chuan-Yan 2014). The recent decline in glacier mass worldwide has important implications for human habitats, including water resource disputes, desertification in areas where glaciers have vanished, and shipping hazards in areas with increased ice loss due to calving from tidewater glaciers (Intergovernmental Panel on Climate Change 2007; Y. Li and Li 2014; Qiang, Lei, and Chuan-Yan 2014). Of particular concern is the relationship between global warming, glacier loss, and sea level rise. Gulf of Alaska glaciers represent an estimated one-third of the global contribution to sea level rise from melting glaciers and ice caps (Arendt et al. 2002; Meier and Dyurgerov 2002). Glaciers in the western Chugach Mountains are among the largest contributors to the rising sea level in the Gulf of Alaska (Berthier et al. 2010).

Comprehensive glacier chronologies and repeat inventories over time for different regions are critical for assessing glacier changes due to the climate or other forces (Nuth et al. 2013; Y. Li and Li 2014; Kienholz et al. 2015). In this study, topographic maps and Landsat imagery were used to document the advance and retreat of land and lake-terminating (land-based) glaciers and coastal tidewater glaciers in the western Prince William Sound (PWS) in the northern Gulf of Alaska (Figure 1). The forty-three evaluated glaciers terminate on land or in fjords along the PWS and are divided into three subregions: the Sargent Icefield, Spencer–Blackstone Ice Complex, and the western Chugach Mountains.

Muir (1902) recorded the earliest observations of glaciers in the area in 1899. Glacier chronologies, extending back to the mid-1700s, exist for some glaciers in the study area (Field 1975; Sturm et al. 1991; Wiles, Barclay, and Calkin 1999; Barclay, Wiles, and Calkin 2003). However, the chronologies only record fluctuations up to the late twentieth century. The aim of this study is to understand how coastal glaciers in the study area fluctuated since the Little Ice Age (LIA) in response to climate change. To fulfill this aim, glacier chronologies, including extent and area, were developed for forty-three coastal glaciers larger than 10 km² (Table 1). The chronologies were based on analysis of topographic maps from ~1950 and Landsat imagery from 1973 and 1975 (1973–1975), 1986, 1994, 2004 and 2006 (2004–2006), and 2018. Further, temporal and spatial changes, including glacier length, area, and mass changes, were analyzed for land-based and tidewater glaciers. Morphometric variables, representing topographic and geometric glacier properties, were also evaluated as potential controls for glacier change. Finally, climate time series and trends were developed and compared with patterns of advance and retreat to understand glacier changes in response to climate change from the 1950s to 2018.

Table 1. General characteristics of study area glaciers

Region	GLIMS ID	Glacier name	Glacier type	Longitude (°)	Latitude (°)	Area (2018) (km^2)	Minimum elevation (m a.s.l.)	Maximum elevation (m a.s.l.)	Median elevation (m a.s.l.)
Chugach Mountains	G211267E61046N		Land-terminating	−148.73	61.05	11	594	1,524	1,041
	G211320E60963N		Land-terminating	−148.66	60.98	14	419	1,878	1,201
	G211348E60944N	Twentymile	Lake-terminating	−148.65	60.95	24	104	1,837	885
	G211353E61066N	Lake George	Lake-terminating	−148.65	61.06	84	143	2,164	1,049
	G211424E60880N	Billings	Land-terminating	−148.57	60.88	10	225	1,651	970
	G211481E60940N	Harriman	Tidewater	−148.53	60.94	50	6	1,808	709
	G211491E61046N	Surprise	Tidewater	−148.50	61.02	75	88	2,213	1,074
	G211615E61175N	Colony	Lake-terminating	−148.39	61.17	210	59	2,897	1,412
	G211805E61150N	Cascade	Tidewater	−148.19	61.16	15	192	2,456	1,463
	G211861E61338N	Knik	Lake-terminating	−148.14	61.34	412	32	3,631	1,743
	G211960E61143N	Coxe	Tidewater	−148.04	61.14	19	0	1,777	1,268
	G211952E61204N	Barry	Tidewater	−148.01	61.22	82	257	2,611	1,349
	G212501E61399N	Harvard	Tidewater	−147.64	61.32	450	0	4,010	1,776
	G212608E61223N	Meares	Tidewater	−147.39	61.23	115	4	2,541	1,449
	G212622E61302N	Yale	Tidewater	−147.34	61.32	129	364	3,599	2,168
Spencer–Blackstone Ice Complex	G211040E60441N		Land-terminating	−148.96	60.49	47	540	1,675	1,022
	G211075E60552N	Trail	Lake-terminating	−148.91	60.55	13	422	1,748	1,310
	G211088E60613N	Spencer	Lake-terminating	−148.87	60.63	74	67	1,810	898
	G211100E60420N	Wolverine	Land-terminating	−148.91	60.41	17	532	1,647	1,297
	G211128E60439N		Land-terminating	−148.87	60.44	17	226	1,654	1141
	G211194E60700N	Portage	Lake-terminating	−148.80	60.72	17	88	1,567	965
	G211241E60622N	Blackstone	Tidewater	−148.77	60.61	31	12	1,616	1,015
	G211244E60493N		Land-terminating	−148.73	60.55	14	245	1,737	910
	G211299E60617N	Beloit	Tidewater	−148.71	60.61	26	65	1,814	981
	G211320E60526N		Land-terminating	−148.68	60.55	12	351	1,824	952
	G211358E60580N	Taylor	Land-terminating	−148.66	60.58	12	95	1,766	934
	G211341E60606N	Cotterell	Land-terminating	−148.66	60.61	15	216	1,429	986
	G211468E60679N	Tebenkof	Land-terminating	−148.53	60.68	19	65	1,383	543
	G211232E60657N	Northland	Land-terminating	−148.77	60.66	21	109	1,377	735
Sargent Icefield	G211126E60158N	Ellsworth	Lake-terminating	−148.86	60.21	108	20	1,826	1,234
	G211143E60281N		Land-terminating	−148.85	60.29	18	540	1,565	1,277
	G211233E60341N		Land-terminating	−148.76	60.35	18	263	1,726	992
	G211266E60059N	Excelsior	Lake-terminating	−148.73	60.11	115	365	1,662	995
	G211299E60379N		Land-terminating	−148.70	60.38	12	223	1,591	985
	G211376E60282N	Chenega	Tidewater	−148.62	60.27	375	106	1,791	979
	G211386E60399N	Langdon	Land-terminating	−148.63	60.41	13	214	1,713	1,042
	G211420E60046N	Puget	Land-terminating	−148.58	60.07	13	415	1,391	1,012
	G211467E60175N	Tiger	Tidewater	−148.55	60.17	53	24	1,359	900
	G211493E60096N	Bainbridge	Land-terminating	−148.52	60.12	40	60	1,605	936
	G211561E60361N	Princeton	Land-terminating	−148.44	60.36	36	131	1,211	623
	G211617E60390N	Ultramarine	Lake-terminating	−148.38	60.39	16	157	797	525
	G211513E60431N	Nellie Juan	Tidewater	−148.49	60.43	41	226	1,383	657
	G211469E60427N	Falling	Land-terminating	−148.51	60.46	12	161	1,163	626

Note. Glaciers identified by glacier name (where available) and Global Land Ice Measurements from Space project identification numbers (GLIMS ID). Latitude and longitude represent glacier centroid. Elevation data based on 2018 outlines.

Setting

Glaciers have occupied fjords in the western PWS since the Pleistocene (Wiles, Barclay, and Calkin 1999). The LIA period of glacial advance began across southern Alaska around 1200 AD and continued through the early twentieth century (Porter 1986; Wiles, Barclay, and Calkin 1999). With the exception of some tidewater glaciers, the vast majority of glaciers in Alaska are now retreating (Molnia 2008). Berthier et al. (2018) applied a glacier mass balance approach to detect changes in Alaska’s Juneau and Stikine icefields and estimated rapid mass loss from 2000 to 2016. Arendt, Walsh, and Harrison (2009) used repeat airborne altimetry to study forty-six glaciers in Alaska and northwestern Canada and concluded that 75 percent of the glaciers were losing mass at an increasing rate since the mid-1990s, compared to an earlier period beginning in the 1950s to 1970s. The study attributed most glacier mass losses to increased summer (May–September) temperature, with likely contributions from winter (October–April) warming, resulting in increased thermal energy and precipitation as rain instead of snow.

Post (1975) described the tidewater glacier cycle as the pattern of slow advance followed by rapid collapse as glaciers pull back from shoals and into deeper water. Meier and Post (1987) studied Columbia Glacier, in the northern PWS, and concluded that “advance and retreat of tidewater glaciers are not related to climate fluctuations in a direct way” (9058) but considered tidewater cycles to be related more to basin and fjord geometry and glacier mass balance. McNabb and Hock (2014) developed a sixty-four-year chronology of glacier length for fifty Alaskan tidewater glaciers to understand the response of tidewater glaciers to climate change. The study relied on digitized glacier outlines derived from topographic maps and Landsat images from 1948 to 2012. The authors concluded that 60 percent of the glaciers are retreating and observed “step change” retreats “where glaciers quickly retreat from a stable position, usually located at a topographic constriction, and stabilize again at another topographic constriction” (McNabb and Hock 2014;165). Tidewater glaciers are influenced by oceanic forces, such as the buoyancy-driven circulation of warm water near the calving front (Motyka et al. 2011), and respond differently to climate change and often with asynchronous advance or retreat compared to nearby land-terminating glaciers (Post et al. 2011). Brinkerhoff, Truffer, and Aschwanden (2017) suggest that ice and sediment dynamics could explain the advance of some tidewater glaciers while most glaciers are retreating.

In the western PWS, the U.S. Geological Survey (USGS) maintains a glacier mass balance monitoring program at Wolverine Glacier, one of the glaciers included in this study. Wolverine Glacier is one of three North American glaciers selected by the USGS for long-term monitoring of glacier response to climate change. Results of long-term mass balance monitoring indicate that Wolverine Glacier is losing mass as a result of warmer summer temperatures (O’Neel et al. 2019). Repeat inventory studies conducted in regions throughout the world show similar trends of recent glacier retreat (Xiao et al. 2007; Błaszczyk, Jania, and Kolondra 2013; Bajracharya, Maharjan, and Shrestha 2014; Osipov and Osipova 2014; Paul and Mölg 2014; Qiang, Lei, and Chuan-Yan 2014; Wei et al. 2014; Winsvold, Andreassen, and Kienholz 2014).

The counterclockwise circulation of the Gulf of Alaska is responsible for the maritime climate in the study area and brings significant moisture into the coastal mountains along the PWS (Beamer et al. 2016). As a result, the region experiences heavy precipitation in the form of rain during the cool summers and snow in the winter. The combination of heavy moisture, cool temperatures, and mountainous topography supports the persistence of coastal glaciers. Mean temperatures for coastal communities in the study area range from −5°C in January to 14°C in July, with average annual snowfall as high as 8,280 mm at Whittier (Alaska Climate Research Center 2019).

Data and methods

Historical topographic maps of the study area provided the spatial data needed to extend glacier length change and area chronologies back to the 1950s. The twenty-one maps obtained for this study are available for download in a georeferenced format from the USGS (2019) allowing for use in geographical information systems (GIS) without further processing. The map publication dates span 1951 to 1960 and are produced at the 1:63,360 scale from aerial photographs acquired from 1948 to 1957 (~1950; Supplementary Table 1). Horizontal accuracy estimates for 1:63,360 scale USGS topographic maps are ~50 m, with some studies indicating much larger errors in Alaska (Adalgeirsdóttir, Echelmeyer, and Harrison 1998; Arendt et al. 2002, 2006; McNabb and Hock 2014).

Landsat images were accessed from an online service portal (Environmental Systems Research Institute 2019). The images are georeferenced and orthorectified by the USGS, allowing for direct integration into a GIS. The images, at 15- to 60-m resolution, provided the spatial data for the repeat measurement of glacier outlines spanning 1973 to 2018 (Supplementary Table 1). Previous studies provided LIA maximums for eight of the land-based glaciers analyzed in this study (Wiles, Barclay, and Calkin 1999; Barclay, Wiles, and Calkin 2003). The process of developing glacier center lines, necessary for the calculation of glacier length changes, required high-resolution digital elevation data provided by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) optical sensor. ASTER Global Digital Elevation Model data (ASTER GDEM.v2) is a product of Japan’s Ministry of Economy, Trade, and Industry and the National Aeronautics and Space Administration (30-m resolution, obtained from NASA/METI/AIST/Japan Spacesystems, and U.S./Japan ASTER Science Team 2011).

Annual climate time series were developed for the Cannery Creek, Cordova MK Smith Airport, Seward Airport, and Wolverine Glacier weather stations in the PWS region (Table 2). The climate data are grouped into May–September and October–April averaging periods to represent the summer ablation and winter accumulation seasons in Alaska (Arendt, Walsh, and Harrison 2009; Kienholz et al. 2015). Raw climate data were acquired from the National Oceanic and Atmospheric Administration (NOAA 2019a). Missing data are common in Alaska climate data sets, requiring gap filling in some cases (Supplementary Table 2). Techniques used to fill gaps in monthly records include filling in missing monthly data with daily data averages for the month, when available, and replacing missing months or months with incomplete data (less than twenty-five days of data) with the series average of that month (McAfee, Guentchev, and Eischeid 2013).

Table 2. Characteristics of study area weather stations

					Length of record
Weather station	Station identification	Altitude (m a.s.l.)	Longitude (°)	Latitude (°)	Temperature	Precipitation
Cannery Creek	USC00501240	25.9	−147.51	61.02	1979–2018	1979–2018
Cordova MK Smith Airport	USW00026410	9.4	−145.45	60.49	1950–2018	1950–2018
Seward Airport	USW00026438	6.7	−149.42	60.13	1950–2018	1950–2018
Wolverine Glacier	USGS15236895	990	−148.92	60.40	1967–2018	1967–2016

Mapping glacier outlines

Randolph Glacier Inventory (RGI Consortium 2017) outline shapefiles, based on imagery from 2006 (Sargent Icefield and Spencer–Blackstone Ice Complex) and 2009 (Chugach Mountain glaciers), provided the initial glacier boundaries. Manual adjustments of the digitized RGI glacier boundaries were based on the interpretation of topographic maps from the 1950s and Landsat images acquired in 1973–1975, 1986, 1994, 2004–2006, and 2018. GIS geoprocessing tools were used to calculate the area for each glacier outline and to measure changes over time (Supplementary Table 3). Glacier outlines are available from the Global Land Ice Measurements from Space (GLIMS) database (Raup et al. 2007; www.glims.org). Issues including the frequency of image collection and persistent cloud cover in the western PWS limited the availability of usable imagery. Several glaciers were excluded from the study due to the presence of debris cover, making it difficult to accuarately digitize glacier boundaries. Each study glacier is identified by a project identification number, including GLIMS and RGI identification numbers as well as glacier name, if available. All topographic maps and Landsat images were either acquired in or projected to UTM Zone 6N with WGS84 datum.

Table 3. Average annual retreat rates for the eight land-based study glaciers with estimated Little Ice Age maximum extents (LIAM date) based on tree-ring dating of moraines

		Length change rate
Glacier name	LIAM date	LIAM to 1950s (m a^−1)	1950s to 2018 (m a^−1)	1950s to 1973–1975 (m a^−1)	1973–1975 to 1986 (m a^−1)	1986 to 1994 (m a^−1)	1994 to 2006 (m a^−1)	2006 to 2018 (m a^−1)
Billings	1797	−9	−25	−64	−15	−1	−15	−8
Cotterell	1891	−25	−25	−17	−30	−33	−14	−35
Langdon	1737	−5	−11	−19	−4	−9	0	−10
Tebenkof	1891	−18	−25	−30	−11	−23	−13	−35
Wolverine	1713	−3	−18	−12	−8	−5	−9	−54
Ellsworth	1855	−12	−50	−44	−86	−20	13	−84
Excelsior	1797	−18	−135	−7	−119	−158	−168	−313
Ultramarine	1889	−17	−28	−9	−12	−34	−37	−65
Average		−14	−40	−25	−36	−35	−30	−76

Notes. Highlighted columns show average change rates before and after the 1950s.

Source. Wiles, Barclay, and Calkin (1999).

Volume–area scaling

Volume–area power law scaling was used to estimate changes in glacier volume. Bahr, Meier, and Peckham (1997) described the following simple relationship between area and volume:

(1) $V=c{S}^{\gamma }$(1)

where glacier volume V is estimated from glacier area S, c is a scaling constant in units of length raised to the power (3 − 2γ), and γ is a dimensionless scaling coefficient. Estimates of scaling parameters have been derived empirically from regional data (e.g., Macheret, Cherkasov, and Bobrova 1988; Chen and Ohmura 1990) and theoretical relationships (e.g., Bahr, Meier, and Peckham 1997; Van de Wal and Wild 2001). Among frequently cited volume–area scaling studies, c estimates vary widely by method and region, ranging from c = 0.12 m^3−2γ based on a global glacier model (Van de Wal and Wild 2001), to c = 0.2055 m^3−2γ from a statistical analysis applied to a large glacier area data set (Radić and Hock 2011). Here, lacking a known c estimate for the study area, the scaling constant c = 0.163 m^3−2γ was derived from the average of the upper and lower ends of range of published c estimates, using γ = 1.375. The scaling equation was applied to the entire population of study glaciers, rather than individual glaciers, to minimize error in volume estimation (Bahr, Pfeffer, and Kaser 2015). Changes in glacier mass, measured in gigatons, were calculated by multiplying ice volumes (km³) by the density of ice and dividing by the time interval.

Glacier mass changes and the equivalent rise in sea level (SLE) were calculated for the forty-three measured study glaciers and upscaled to all unmeasured glaciers throughout the PWS, following an approach described by Arendt et al. (2006). The unmeasured glaciers included all 2,828 PWS glaciers available in the GLIMS database (Figure 1), excluding the forty-three measured study glaciers. With the exception of Columbia Glacier, all of the unmeasured glaciers in the PWS were land-based glaciers. Shapefiles of unmeasured glacier outlines, from the mid- to late 2000s, were obtained from the GLIMS database and volume–area scaling was used to calculate an aggregate volume for the entire PWS. In order to estimate aggregate area and volume changes over time, measured glacier area change rates, obtained from analyzing the forty-three study glaciers, were applied to unmeasured glaciers backward to the 1950s and forward to 2018. Columbia Glacier was included in the aggregate volume estimate using a combination of reported volume and area estimates from 1957 to 2007 and extrapolation of measured glacier area and volume changes across 2007 to 2018. Columbia Glacier is the region’s largest contributor to sea level rise and has lost more than 130 km³ of volume since 1957 (Rasmussen et al. 2011; McNabb and Hock 2014). Measured area estimates for Columbia Glacier were obtained from the GLIMS database (2007 outline) and manual mapping of Landsat imagery (2018 outline).

Figure 1. Study area and glacier area loss between the 1950s and 2018

Glacier length

Glacier center lines were calculated using a semi-automatic, GIS-based algorithm (Le Bris and Paul 2013), resulting in a center line from the highest point of the glacier to the terminus. Glacier length changes were measured from the intersection of the center line with each glacier terminus, as defined by digitized glacier outlines. Repeat measurements for 1950s topographic maps and the Landsat images acquired in 1973–1975, 1986, 1994, 2004–2006, and 2018, resulted in a glacier length change chronology for each glacier (Supplementary Table 3). For a subset of eight glaciers, length change measurements extended back to the digitized LIA maximum terminus positions identified in previous studies (Table 3).

Morphometric and bathymetric analysis

Following an approach established by Y. Li and Li (2014), six derived morphometric variables representing topographic and geometric glacier properties were used to evaluate potential controls on glacier change. Table 4 includes the list of derived variables. Relative glacier area change (ΔA) was calculated to quantify glacier change from the 1950s to 2018. Glacier area is commonly used to represent change, because GIS allows for easy calculation of glacier area change across time periods (Y. Li and Li 2014). Relationships between ΔA and the six morphometric variables were quantified using multiple regression analysis and interpretation of correlation coefficients (R²) and p values.

Table 4. List of six derived topographic and geometric variables

Derived variable		Description	Data sources
Elevation mean		Average elevation of all DEM cells within the glacier outline.	ASTER GDEM
Aspect		Orientation of each DEM cell within a glacier outline as an azimuth (°) relative to 0° at due north	Randolph Glacier Inventory (based on imagery from 2006 to 2009; RGI Consortium 2017)
Slope		Average of slopes (°) of surfaces of each DEM cell within a glacier outline	Randolph Glacier Inventory (based on imagery from 2006 to 2009; RGI Consortium 2017)
Hypsometric interval		Difference between mean elevation and minimum elevation divided by relief (maximum elevation − minimum elevation) within the glacier outline (Pike and Wilson 1971; Y. Li and Li 2014)	ASTER GDEM
Shape index	Derived from $\sqrt{(4\pi \mathrm{x}A/P^2)}$, where A is the area and $P$ is the perimeter of a digitized glacier outline (Y. Li and Li 2014)		Landsat 7 ETM+ and Landsat 8 OLI/TIRS (July 2018)
Area (2018)		Glacier area based on digitized glacier outlines	Landsat 7 ETM+ and Landsat 8 OLI/TIRS (July 2018)

Note. GDEM = Global Digital Elevation Model; ETM+ = Enhanced Thematic Mapper Plus; OLI/TIR = Operational Land Imager/Thermal Infrared Sensor.

Meier and Post (1987) observed a relationship between rapid retreat rates and water depth in tidewater glaciers in Alaska. Fjord bathymetric data were obtained from NOAA (NOAA 2019b; https://maps.ngdc.noaa.gov/viewers/bathymetry/) to test this relationship in the study area. Fjord bathymetric data were derived from hydrographic surveys conducted from 1990 to 1997. Based on availability of depth soundings, mean water depths were calculated for fjord areas between the retreating terminus positions of eight tidewater glaciers, including multiple time periods for Beloit, Coxe, Nellie Juan, Tiger, and Yale glaciers (Table 5). Regression analysis was applied to assess the relationship between retreat rates and mean water depth and the potential influence of water depth on tidewater glacier fluctuation. Fjord width is another important controlling factor for tidewater glacier retreat and advance (Mercer 1961). To test this relationship, fjord widths were measured at the 2004–2006 and 2018 terminus of each tidewater glacier and compared to corresponding length changes. The 2004–2006 to 2018 interval represents the peak period of average tidewater glacier retreat in the study, including the rapid collapse of Barry, Nellie Juan, and Surprise glaciers.

Table 5. Comparison of retreat rates and mean water depth for eight tidewater glaciers in the study area

Glacier	Time period	Rate of retreat (m a^−1)	Mean water depth (m)
Barry–Coxe	1898–1957	83	78
Beloit	1950–1975	3	13
Beloit	1975–1986	45	26
Blackstone	1975–1986	6	14
Coxe	1910–1957	43	90
Nellie Juan	1950–1975	46	28
Nellie Juan	1975–1986	50	42
Nellie Juan	1986–1994	57	48
Surprise	1898–1975	30	63
Tiger	1950–1973	16	53
Tiger	1994–2006	31	53
Yale	1810–1910	49	72
Yale	1910–1957	6	32
Yale	1957–1975	135	82
Yale	1975–1986	179	133
Yale	1986–1994	44	80
Yale	1994–2004	53	88

Uncertainty

The sources of uncertainty relating to the use of remote sensing in glacier studies are well documented and include satellite image quality and resolution issues, georeferencing issues, manual editing errors (e.g., digitizing outlines), and the analyst’s degree of skill in interpreting features such as ice flow divides, snowlines, glacier boundaries, and debris cover (DeBeer and Sharp 2007; Racoviteanu et al. 2009; Paul et al. 2013, 2017; McNabb and Hock 2014; Osipov and Osipova 2014; Paul and Mölg 2014; Pfeffer et al. 2014; Qiang, Lei, and Chuan-Yan 2014; Kienholz et al. 2015). The following study design criteria were applied to minimize potential errors: (1) exclude glaciers with heavy debris cover at margins; (2) avoid satellite images with significant cloud cover or shadows along glacier margins; (3) select late summer (July–September) satellite images to minimize seasonal snow; and (4) obtain georeferenced and orthorectified maps and satellite images, allowing for direct integration into GIS, and check for geolocation errors.

Strategies recommended by Paul et al. (2017) were followed to assess the accuracy and precision of glacier outlines and lengths derived from satellite imagery, including the accuracy measures outline overlay, literature values, and multiple digitizing. Outline overlay involves a qualitative analysis of multiple glacier outlines over the source image and helps identify potential sources of error. This method, applied to six study glaciers, revealed differences in outlines in the range of 30–90 m (one to three pixels). The errors were in discrete areas of glacier margins where boundaries between snow and ice, or debris and ice, were difficult to assess.

Literature values provide an additional source of accuracy estimates for digitizing glacier outlines. Uncertainty estimates for manual glacier outline editing range from ~2 percent to 5 percent of the area measured (Paul et al. 2013; Bajracharya, Maharjan, and Shrestha 2014; Paul and Mölg 2014; Pfeffer et al. 2014). McNabb and Hock (2014) estimated a horizontal accuracy of 60 m for Landsat images, incorporating digitizing and georeferencing errors. Analyst precision testing performed by Fischer et al. (2014) resulted in a predicted precision of ±3.1 percent of the area for glaciers larger than 5.0 km². Previous studies show a decrease in estimates of uncertainty for glacier area measurement as glacier area increases (Fischer et al. 2014; MacHguth and Huss 2014; Pfeffer et al. 2014). Analyst precision testing in this study followed the methodology recommended by Paul et al. (2017), including the selection of six glaciers of different sizes (10–77 km²) and three iterations of manual outline digitizing for each glacier. The calculated standard deviation ranged from 0.16 to 0.42 km², and the percentage precision ranged from 0.7 to 2.0 percent. For comparison purposes it is important to note that this study avoided glaciers with debris cover along margins, whereas other studies may include all debris-covered glaciers.

Bahr (1997) and Bahr, Meier, and Peckham (1997) recognized log–log scaling relationships for glacier volume, area, and length. Establishing area–length scaling relationships allows for the comparison of differences between regions and methods of glacier center line development and provides a simple method for estimating glacier length from area measurements (Le Bris and Paul 2013; MacHguth and Huss 2014; Kienholz et al. 2015). Area–length scaling for the thirty land-based glaciers in this study resulted in the log–log relationship of L = 1.7836A^0.5416 (R² = 0.898), similar to the relationship established by Le Bris and Paul (2013), who developed the center line length calculation method applied in this study (Supplementary Figure 1).

The large uncertainty resulting from the volume–area power law scaling approach was estimated by calculating glacier volume using the upper (c = 0.2055 m^3−2γ; Radić and Hock 2011) and lower (c = 0.12 m^3−2γ; Van de Wal and Wild 2001) ends of the range of published c estimates. The uncertainty associated with the forward and backward extrapolation of area estimates for unmeasured glaciers was calculated as the percentage difference between measured and extrapolated average area change rates for the study glaciers from 1950–1957 to 2004–2006 and 2004–3006 to 2018. Average area change rates were derived from measured study glaciers larger than 10 km². This approach underestimates volume–area scaling results because relative area losses were greater for smaller measured glaciers, and over 98 percent of unmeasured glaciers were smaller than 10 km². The uncertainty was estimated as 3.8 percent for backward (2004–2006 to 1950–1957) and 2.9 percent for forward (2004–2006 to 2018) extrapolation.

Climate analysis

Climate time series were developed for four weather stations in the PWS region to assess the relationship between glacier changes and climate. Time series included temperature and precipitation for each station for annual, ablation season, and accumulation season averages. Detailed climate time series methodology is provided in the Supplementary Material. Previous studies suggest a relationship between abnormally warm summer sea surface temperatures and the onset of rapid glacier retreat (Schild and Hamilton 2013; McNabb and Hock 2014). To further investigate potential factors on tidewater glacier fluctuations, sea surface temperature (SST) data were obtained from NOAA (NOAA 2019c; https://tidesandcurrents.noaa.gov/stations.html) for the Cordova, Seward, and Valdez stations. Mean summer (June–August) SSTs were calculated for Valdez station from 1994 to 2018. Gaps in the water temperature record for Cordova and Seward limited the extent of the mean summer SST analysis to 2015.

Temperature in Alaska is influenced by climate signals including the Pacific Decadal Oscillation (PDO) and shifts in atmospheric circulation in the Arctic. The PDO refers to a cyclical pattern of SST and sea level pressure changes occurring in the North Pacific over periods of ten years or more. North Pacific SST anomalies define the PDO and shift between positive and negative phases. The positive PDO phase is associated with warmer and wetter conditions along the coast of Alaska (Arendt, Walsh, and Harrison 2009). The Arctic Oscillation (AO) refers to the changing strength of the polar vortex (Thompson and Wallace 1998), measured as an index based on sea level pressure anomalies. In the positive AO phase, strong circulation at the North Pole contains the spread of colder air to the midlatitudes, leading to cooler temperatures and increased winter precipitation in Alaska. PDO and AO data from 1950 to 2018 were acquired from NOAA (NOAA 2019d; https://www.ncdc.noaa.gov/teleconnections/) and used in regressions to quantify relationships with temperature and precipitation. Time series were developed to provide a visual comparison of glacier changes and significant shifts in PDO and AO index phases.

Results

Length change metrics document the horizontal retreat of land-based study glaciers from ~1950 to 2018 and a period of accelerated retreat from 2004–2006 to 2018. Land-based glaciers retreated by an average of 1.8 km (14 percent of total length) since the 1950s (Figure 2), with minimal spatial variation (~1 percent) across the three subregions. The average rate of retreat for all land-based glaciers more than doubled to 48 m a⁻¹ after 2004–2006 (Figure 3a, Table 6). Based on mapped LIA extents (Wiles, Barclay, and Calkin 1999), glacier retreat rates averaged 14 m a⁻¹ from the LIA maximum to the 1950s (Table 3). Average retreat rates for this subset of glaciers peaked after 2004–2006, consistent with average retreat rates for all other glaciers in this study. The combined area of all land-based glaciers decreased by 228 km² (17 percent loss) from ~1950 to 2018 (Figure 1), with 36 percent of the land-based glacier area loss occurring from 2004–2006 to 2018. The percentage of area loss was similar across the three subregions. The total glacierized area change rates and average rate of area loss for land-based glaciers peaked from 2004–2006 to 2018 (Figures 3b–3c). During this period, 80 percent of land-based glaciers experienced the highest annual area loss rate (Figure 4). The regression between glacier relative area change (1950s–2018) and the six derived morphometric variables, representing topographic and geometric glacier properties, indicates a weak relationship (R² = 0.17). Among the six variables, only glacier area showed a significant correlation (p = .03) with glacier relative area change (Table 7).

Table 6. Average length and area change metrics for land-based and tidewater glaciers

Metric		Glacier type	Unit	1950s to 1973–1975	1973–1975 to 1986	1986 to 1994	1994 to 2004–2006	2004–2006 to 2018
Length	Absolute change	Land-based	m	−451	−276	−199	−247	−601
		Tidewater	m	−175	−262	13	−172	−336
	Rate of change	Land-based	m a^−1	−20	−24	−25	−21	−48
		Tidewater	m a^−1	−8	−23	2	−15	−26
	Total cumulative change	Land-based	%	−4.2	−6.2	−7.7	−9.5	−14.1
		Tidewater	%	−0.6	−1.5	−1.4	−2.3	−4.3
Area	Absolute change	Land-based	km^2	−1.6	−1.1	−0.7	−1.4	−2.8
		Tidewater	km^2	−1.5	−2.0	−0.9	−1.2	−1.2
	Rate of change	Land-based	km^2 a^−1	−0.07	−0.09	−0.08	−0.12	−0.21
		Tidewater	km^2 a^−1	−0.08	−0.18	−0.11	−0.10	−0.09
	Total cumulative change	Land-based	%	−3.6	−6.2	−7.7	−11.3	−17.4
		Tidewater	%	−1.2	−2.5	−3.2	−4.6	−6.1

Table 7. Regression between glacier area change from 1950s to 2018 and six derived morphometric variables

	ΔArea
Explanatory variable	Estimate	p
Elevation mean	−0.01	.39
Aspect	−0.01	.51
Slope	0.60	.25
Hypsometric interval	51.96	.06
Shape index	−7.23	.53
Area	0.05	.03

Note. Bold value indicates significance at the 95 percent level (p < .05).

Mass of the forty-three measured study glaciers declined by 63 ± 10 Gt from the 1950s to 2018, equivalent to a sea level rise of 0.174 mm. Simple upscaling of area and volume changes to unmeasured glaciers across the PWS resulted in an estimated aggregate glacier mass loss for the PWS of 379 ± 38 Gt, equivalent to a 1.047 mm rise in sea level from the 1950s to 2018 (Table 8). The catastrophic retreat of Columbia Glacier accounted for over half (201 ± 11 Gt) of the mass loss for the PWS. Since the 1950s, over 37 percent of the aggregate PWS glacier mass loss occurred after 2004 with an equivalent sea level rise.

Table 8. Glacier mass (gigatons) and SLE changes from 1950–1957 to 2018

	n	Area (km^2)^a	Glacier mass change (Gt)		SLE rate (mm yr^−1)		SLE change (mm)
Measured glaciers	43	3,215	−63	±10	−0.003	±0.001	−0.174	±0.046
Unmeasured glaciers	2,827	4,651	−115	±18	−0.005	±0.001	−0.318	±0.084
Columbia Glacier	1	1,068	−201	±11	−0.009	±0.001	−0.555	±0.030
Total	2,871	8,934	−379	±38	−0.017	±0.003	−1.047	±0.160

Note. ^aArea refers to initial 1950–1957 estimates. SLE = sea level equivalent.

Tidewater glaciers behaved asynchronously compared to land-based glaciers, as indicated by differing periods of peak area and length loss and a much lower average rate of retreat. The combined area of study tidewater glaciers declined by 88 km² or 6.1 percent since the 1950s (Figure 1; Table 6), with peak losses occurring between 1973–1975 and 1986 (Figure 3b). Tidewater glaciers lost an average of 4.3 percent of total length, compared to a 14.1 percent loss for land-based glaciers since the 1950s (Figure 2). The average retreat rate for tidewater glaciers fluctuated across the study period, with highest rates occurring in 2004–2006 to 2018 and 1973–1975 to 1986 (Figure 3a). The period 1986 to 1994 represents the only study interval with positive average length change rates (2 m a⁻¹) for all tidewater glaciers, with four glaciers reaching peak advance rates (Figure 5). Nellie Juan Glacier provides an example of the tidewater glacier cycle with rapid retreat of over 5 km since the glacier withdrew from the terminal moraine shoal after 1935 (Barclay, Wiles, and Calkin 2003; Figure 6). Nellie Juan Glacier’s period of rapid retreat peaked at a rate of 124 m a⁻¹ from 2004–2006 to 2018. During the same period, Barry Glacier retreated over 2 km, into a widening fjord, from a relatively stable position at the head of the inlet shared with Cascade and Coxe glaciers. The glaciers of College Fjord highlight the asynchronous behavior of tidewater glaciers. Since the 1950s, Harvard Glacier advanced slowly, whereas the neighboring Yale Glacier retreated over 5 km (Figure 6).

Figure 2. Percentage glacier length and area change since the 1950s. Vertical bars represent calculated uncertainty (one standard deviation)

Figure 3. (a) Average annual glacier retreat rates; (b) average annual glacier area change rates; and (c) total glacierized area annual change rates. Change rates presented by glacier type (land-based n = 30; tidewater n = 13)

Figure 4. Period area change rates (bars) and glacier area in 2018 (dashed line) for land-based glaciers

Figure 5. Length change rates for individual tidewater glaciers by time period

Figure 6. (a) Chronology of retreat for Nellie Juan Glacier between 1935 and 2018. Terminus position for 1935 from photographs in Field (1975; as cited in Wiles, Barclay, and Calkin 1999) and (b) chronology of retreat for Harvard and Yale glaciers between 1810 and 2018. Yale Glacier 1810 terminus position from Field and others (1958; as cited in Mercer 1961) and Harvard Glacier 1910 terminus position from Tarr and Martin (1914; as cited in Sturm et al. 1991)

Climatology

Climate trend analysis included weather station precipitation and temperature data for the four study area weather stations (Supplementary Table 4). Trend analysis results indicate increasing average annual air temperatures in the study area, including an increase from 3.1°C to 4.7°C at Cordova from 1951 to 2017 (Figure 7a). Significant winter warming occurred at PWS maritime stations (+0.9°C) from 1988 to 2017 (Figure 7b). This increase is noteworthy because the average PWS maritime winter temperature is near freezing (−0.5°C). Increasing winter temperatures may expose larger areas of maritime glaciers to temperatures above freezing (Arendt, Walsh, and Harrison 2009). Average summer temperature increased from 5.4°C to 6.1°C from 1968 to 2017 at the higher elevation Wolverine Glacier weather station (Figure 7c). Annual winter precipitation decreased by an average of 319 mm (20 percent decline) at PWS maritime stations from 1976 to 2017 (Figure 7d). Analysis of summer water temperature data obtained from Valdez station indicate a 2.8°C rise from 2008 to 2016 (Figure 8b). Data from Cordova and Seward showed similar trends of increasing summer water temperatures beginning in the late 2000s.

Figure 7. Time series of (a) mean annual temperature at Cordova Airport station; (b) average winter temperature (PWS maritime stations); (c) mean summer temperature at Wolverine Glacier Station; and (d) average winter precipitation (PWS maritime stations). Dashed lines show statistically significant trends

Figure 8. Comparison between (a) annual summer (May–September) and winter (October–April) temperature and precipitation; (b) mean annual summer sea surface temperatures; (c) Pacific decadal oscillation (PDO); (d) Arctic Oscillation (AO); and (e) length change rates for land-based and tidewater glaciers

Figure 9. (a) Relationship between mean fjord water depth and retreat rates for tidewater glaciers and (b) comparison between tidewater glacier length and fjord width changes from 2004–2006 to 2018. Note, measurements based on early August 2004–2006 and late July 2018 imagery

Climate indices provide evidence of significant climate phase shifts in Alaska, with a period of positive PDO from 1977 to 1997 and a large shift to positive AO in 1989 (Arendt, Walsh, and Harrison 2009). Regression analysis was conducted to identify relationships between PDO, AO, air temperature, and precipitation from 1950 to 2018. No significant correlations were identified between PDO and AO. Regression between PDO and annual temperature at Cordova indicates a strong correlation from 1952 to 1989 (R² = 0.46; p < .0001); however, the correlation weakens from 1990 to 2017 (R² = 0.23; p = .012), suggesting a shift in forces responsible for recent warming in Alaska. Weak correlations were observed between winter PDO and precipitation at Seward (R² = 0.12; p = .026) from 1976 to 2017 and Wolverine Glacier (R² = 0.15; p = .034) from 1986 to 2015.

Discussion

The majority of coastal glaciers in the study area are retreating, with rapid mass loss occurring from the mid-2000s to 2018. These findings are consistent with reported glacier mass loss in Alaska (Arendt et al. 2006; Molnia 2008; Arendt, Walsh, and Harrison 2009; Berthier et al. 2010; Yang et al. 2020) and western North America (Menounos et al. 2019). O’Neel et al. (2019) measured an area loss of 1.5 km² for Wolverine Glacier from 1971 to 2018, which is similar to the findings in this study (1.4 km² over 1975–2018), and attributed mass loss to warmer summer temperatures. Both studies estimated the highest rates of area loss occurring from 2006 to 2018 at Wolverine Glacier.

Yang et al. (2020) compiled new inventories for glaciers in the Kenai Peninsula, including the Sargent Icefield and Spencer–Blackstone Ice Complex, and found a significant loss in total glacier area (~12 percent) from 1986 to 2016. Here, measured glaciers in the Sargent Icefield and Spencer–Blackstone Ice subregions lost ~8 percent of their total area from 1986 to 2018. The lower percentage of area loss compared to Yang et al. (2020) is attributed to the exclusion of glaciers smaller than 10 km². Other studies in Alaska (Meier and Dyurgerov 2002), the Northern Alps (Paul et al. 2004), and the Chinese Tien Shan (K. Li et al. 2011; Y. Li and Li 2014) suggest that glacier area as a key indicator of glacier loss, with smaller glaciers tending to have greater relative reductions.

Arendt et al. (2006) used airborne altimetry measurements and multiple extrapolation methods, including special considerations for tidewater glaciers, to estimate a regional contribution of 0.02 mm yr⁻¹ SLE and a net balance rate of −7.4 ± 1.1 km³ yr⁻¹ water equivalent for glaciers covering 9,300 km² of the western Chugach Mountains and Talkeetna Mountains from 1950–1957 to 2001–2004. Here, simple upscaling to a similar total aggregate glacier area of 9,600 km², including glaciers in the western Chugach Mountains, Sargent Icefield, and Spencer–Blackstone Ice Complex, results in an estimated rate of 0.013 mm yr⁻¹ SLE and a net balance rate of −4.6 ± 0.5 km³ yr⁻¹ water equivalent from 1950–1957 to 2004–2010.

The asynchronous glacier cycles observed among tidewater suggest the influence of factors other than climate, such as fjord geometry, water depth, mass balance, and sediment dynamics (Meier and Post 1987; Brinkerhoff, Truffer, and Aschwanden 2017). Warming waters in the PWS may have contributed to the increased average rate of tidewater glacier retreat observed since the mid-2000s (Figures 8b, Figure 8e), and the rapid collapse of Nellie Juan and Barry glaciers. Regionally, SST anomalies throughout the Gulf of Alaska increased by up to 3°C to 4°C since 2013 (Campbell 2018). Other studies from Alaska and Greenland stress the importance of ocean water temperatures as a controlling factor for submarine glacial melting and calving rates (Motyka et al. 2003, 2011; Post et al. 2011).

Earlier studies linked rapid retreat of glaciers in the Gulf of Alaska to water depth and emphasized the importance of nonclimatic controls on tidewater glacier fluctuations (Post 1975; Meier and Post 1987). Here, regression analysis between tidewater glacier retreat rate and mean water depth (Figure 9a) results in a significant positive relationship (R² = 0.58; p < .001). Water depth may explain the rapid retreat of Yale Glacier, occurring from 1973 to 1986. During this period, Yale Glacier retreated over the deepest waters in Yale Arm at a rate of 179 m a⁻¹ (Table 5). Further investigation into fjord geometry suggests that fjord width may be another controlling factor for tidewater glacier advance and retreat. From the mid-2000s to 2018, all four advancing tidewater glaciers moved into a narrowing fjord, whereas Nellie Juan and Barry glaciers experienced the most rapid retreat in fjords that were rapidly widening (Figure 9b). The relationships presented here suggest that factors other than climate are influencing tidewater glacier fluctuations. McNabb and Hock (2014) determined that tidewater glacier advance and retreat patterns in Alaska are highly variable and noted a lack of evidence for a coherent regional-scale retreat. Here, the lower average rates of tidewater glacier retreat and area loss compared to land-based glaciers are attributed to nonclimatic controls and the variability in tidewater glacier cycle position for the thirteen studied tidewater glaciers.

Increasing temperatures and decreasing winter precipitation explain overall glacier recession in the study area. The accelerated retreat after 2004–2006 is likely attributable to an unfavorable combination of temperature and precipitation changes that occurred before and after 2004–2006. Figure 8 shows the relationship between the average length change rates, temperature precipitation trends for summer and winter seasons, summer SSTs, and climate indices. Summer temperatures and SSTs increased from ~2012 to 2017. To further drive accelerated melting after 2004–2006, winter temperature increased sharply at all stations from ~2012 to 2016, and summer precipitation increased at all stations from ~2009 and ~2013. Preceding this period of increased summer precipitation, warming, and accelerated glacier retreat, a steep decline in winter season precipitation began around 2002 and continued to ~2010 (Cannery Creek and Seward).

The brief shifts to negative PDO and positive AO after 1988 coincide with a period of positive average length change rates for all tidewater glaciers (Figures 8c–Figure 8e). However, this relationship does not extend to the recent period of increased temperature and accelerated glacier retreat. After 2004, PDO was strongly negative with the exception of a short period of positive PDO from 2014 to 2016 (Figure 8c). Arendt, Walsh, and Harrison (2009) suggested that forces other than climate indices may be responsible for warming after 1999, when AO began to trend neutral and PDO shifted to the negative phase.

Conclusions

Analysis of topographic maps from the 1950s and Landsat imagery from 1973–1975, 1986, 1994, 2004–2006, and 2018 supported the development of repeat glacier inventories for forty-three glaciers in the western Prince William Sound, Alaska. Coastal glaciers in study area have been retreating since the end of the LIA, with a period of accelerated retreat and mass loss from the mid-2000s to 2018. Simple upscaling to the entire PWS results in an estimated aggregate glacier mass loss of 379 Gt, equivalent to a 1.047 mm rise in sea level from the 1950s to 2018. Glacier retreat correlates with climate trends including increased summer and winter temperatures and decreased winter precipitation over the last several decades. Sequencing of declining winter precipitation before 2004–2006, followed by increasing summer precipitation and summer and winter temperatures, may have forced accelerated melting after the mid-2000s. The lack of significant statistical relationships between topographic and geometric glacier properties and relative glacier area change further supports the suggestion that climate is the key driver for land-based glacier loss in the study area.

Tidewater glaciers behaved asynchronously compared to land-based glaciers, as indicated by differing periods of peak area and length loss and a much lower average rate of retreat (Figure 3). However, though some tidewater glaciers are advancing, the combined area of all tidewater glaciers in the study declined by approximately 6 percent since the 1950s. This overall decline in tidewater glacier area is likely driven by factors other than climate, such as fjord geometry, bathymetry, and warming summer SSTs. However, climate may have triggered cycles of retreat and indirectly contributed to the overall loss in tidewater glacier mass because warming waters in the PWS likely increased calving rates and the rapid collapse of Barry, Nellie Juan, and Yale glaciers. The tidewater glacier cycle and high variability in advance and retreat patterns for the thirteen study tidewater glaciers explains the lower average rates of tidewater glacier retreat and area loss compared to land-based glaciers.

Expanding the study to the eastern PWS and along the Gulf of Alaska may provide insights regarding the extent of recent accelerated retreat documented in this study. Incorporating remote sensing–based or physical measurements of glacier thinning and volume changes would provide a more accurate estimate of glacier response to climate change in the study area. Further glacier-specific studies are warranted to better understand the complex interaction of climate, fjord geometry, and tidewater glacier mass balance.

Acknowledgments

I thank Gregory Wiles for his suggestions and support during the research design process. I also thank Andrew Phillips for his constructive comments and helping to improve earlier drafts of this article. This work is not a product of the United States Government or the U.S. Environmental Protection Agency. The author has not carried out this work in any governmental capacity. The views expressed are his own and do not necessarily represent those of the United States or the U.S. Environmental Protection Agency.

Disclosure statement

No potential conflict of interest was reported by the author.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.