Surface Elevation Change of Glaciers Along the Coast of Prudhoe Land, Northwestern Greenland From 1985 to 2018

Outlet glaciers in Greenland exhibited rapid frontal retreat, thinning, and flow acceleration in the 21st century. In contrast to a number of studies in the last two decades, long‐term glacier changes are not well known, hence limiting the understanding of ongoing glacier retreat. In this study, we present surface elevation changes of 16 glaciers along the coast of Prudhoe Land, northwestern Greenland, using digital elevation models between 1985 and 2018. All the glaciers experienced surface lowering at a mean rate of −0.55 ± 0.24 m a−1 over the period of study. Surface elevation change has shown a clear shift from 0.14 ± 0.17 m a−1 in 1985–2001 to −1.31 ± 0.20 m a−1 in 2001–2018. Thinning after 2001 was more significant at the glaciers terminating in Inglefield Bredning, represented by the most rapid thinning of Tracy and Farquhar Glaciers at −3.91 ± 0.13 and −2.91 ± 0.18 m a−1, respectively. Glaciers terminating in shallower fjords directly connected to Baffin Bay showed a thinning rate 40% slower than those in the Inglefield Bredning region. Warming trends in atmospheric and ocean temperature since the late 1990s are most likely triggers of the glacier regime shift around 2001. The heterogeneity in the glacier change is attributed to the glacier geometry and fjord bathymetry of each individual glacier, since glaciers terminating in deep fjords are subjected to greater acceleration and are more affected by deep ocean warming.

To quantify glacier mass change over a large area on a decadal temporal scale, remote sensing is a suitable and commonly employed approach. Near the coast of Greenland, repeated laser altimetry by NASA's Airborne Topographic Mapper (ATM) indicated thinning of the ice during the 1990s at a rate exceeding 1 m a −1 (Krabill et al., 2000(Krabill et al., , 2004. Based on the combination of ATM and subsequent ICESat altimetry data sets, Csatho et al. (2014) revealed that 48% of the ice sheet mass loss from 1993 to 2012 was due to dynamic thinning (vertical straining due to accelerated ice flow). The laser altimetry provides a precise measure, but the incomplete data between satellite tracks limits spatial resolution of elevation analysis over a relatively small area. As an alternative method for measuring glacier elevation change in Greenland, digital elevation models (DEMs) derived from satellite or aerial images have been applied in recent years. By means of DEM differencing, mass loss of outlet glaciers has been evaluated in northeastern Greenland (Khan et al., 2014), on Melville Bay in northwestern Greenland  and along the entire margins of the Greenland Ice Sheet (Kjeldsen et al., 2015). The results indicate that glacier mass loss has increased over the Greenland Ice Sheet since the 2000s. More recently, surface elevation change of 1,526 peripheral glaciers and ice caps in west-central Greenland has been quantified as −0.5 ± 0.2 m a −1 for the period from 1985 to 2012 (Huber et al., 2020). However, the timing of the mass loss and its magnitude are highly heterogeneous across Greenland and for each glacier (e.g., Khan et al., 2010;Moon et al., 2012). Thus, it is crucial to perform a detailed regional analysis by taking climatic and geometrical settings into consideration. Marine-terminating glaciers exhibit a particularly high degree of spatial variability, even within a relatively small region (Porter et al., 2018). For example, 37 marine-terminating glaciers in central east Greenland showed a significantly large variation in surface elevation change, ranging from −156 to +7 m between 2000 and 2010. The authors concluded that coastal ocean heat transport was primarily responsible for the glacier change (Walsh et al., 2012). Therefore, comparison of glacier changes under different settings in terms of fjord geometry, ocean temperature, circulation, and stratification, which control ice-ocean interactions (Straneo et al., 2010), may help us better understand the mechanisms controlling glacier change.
Tidewater glaciers in northwestern Greenland have shown rapid thinning Pritchard et al., 2009), significant acceleration (Moon et al., 2012), and persistent terminus retreat (Bunce et al., 2018) since the 2000s, and are important contributors to the mass loss of the Greenland Ice Sheet (Mouginot et al., 2019). However, the previous studies focused on Melville Bay, located in the southern part of northwestern Greenland. More recently, Sugiyama (2018, 2020) reported variations in the ice front positions and flow speed of marine-terminating glaciers along the coast of Prudhoe Land, which is situated 300 km north of Melville Bay (Figure 1). Synchronized glacier retreat and acceleration began around 2000 , approximately the same time as the rapid change of the glaciers in Melville Bay (Bunce et al., 2018;Moon et al., 2012). Surface elevation change has been reported for some of the glaciers along the coast of Prudhoe Land, e.g., Bowdoin and Tugto Glaciers , Tracy and Heilprin Glaciers (Porter et al., 2014;Willis et al., 2018), which highlights a significant heterogeneity in the magnitude and timing of the glacier change. Despite the effort of the previous studies, observations of surface elevation change in the region are limited to the individual glaciers. Moreover, elevation change over the 20th century has not previously been reported in the Prudhoe Land region.
In this study, we used the recently released AeroDEM (Korsgaard et al., 2016) and ASTER-VA DEM (Fujisada et al., 2005;Hirano et al., 2003) to assess the surface elevation changes of all the glaciers along the coast of Prudhoe Land. We analyzed spatial variations in the elevation change from 1985 to 2018, as well as the change in the trend between two subperiods (1985-2001 and 2001-2018). Each individual glacier was assessed, as well as the glacier groups subdivided into those terminating in Inglefield Bredning and those facing Baffin Bay ( Figure 1). Subsequently, we discuss possible drivers of recent rapid glacier change, as well as controls of spatial and temporal patterns of the glacier mass loss.

Study Site
We studied 16 outlet glaciers of the Greenland Ice Sheet located along the coast of Prudhoe Land in northwestern Greenland (Figure 1). The glaciers are distributed across the area of 77.45-78.03°N and 65.95-72.03°W. Of all the glaciers studied, only Tugto Glacier terminates on land; the other 15 are marine-terminating glaciers. Glaciers in the eastern part (Heilprin, Tracy, Farquhar, Melville, Sharp, Hart, Hubbard, and Bowdoin) feed into Inglefield Bredning, a fjord with a width of approximately 20 km and a length of 80 km. Glaciers located in the western part of the study site (Sun, Verhoeff, Meehan, Morris Jesup, Diebitsch, Clements, and Bamse) feed into relatively small fjords, which have short and direct connections to Baffin Bay.
The marine-terminating outlet glaciers along the coast of Prudhoe Land showed frontal retreat between the 1980s and 2014. Most of the glaciers began retreating after 2000, with substantial acceleration observed at those terminating in Inglefield Bredning. From 2000 to 2014, Tracy, Heilprin, Farquhar, and Bowdoin Glaciers retreated at a rate greater than 80 m a −1 . Tracy, Melville, Farquhar, and Heilprin Glaciers accelerated at a rate in excess of 10 m a −2 . Glaciers in the region showed seasonal ice speed variations with an initial speedup between late May and early June, a peak speed between late June and early July, and subsequent decrease from July to September, while the timing and magnitude of the speedup varies among these glaciers .
The largest two glaciers in this area, Heilprin and Tracy Glaciers located in the easternmost part of Inglefield Bredning (Figure 1), have a width of 5.7 and 5 km, respectively (Hill et al., 2017). They are separated by only 20 km in distance but have shown significantly different patterns of frontal variations. Tracy Glacier has retreated by ∼15 km since the 19th century, whereas the retreat of Heilprin Glacier during the same period was only several kilometers (Dawes & van As, 2010). From 1980 to 2014, the frontal displacement rates of Tracy and Heilprin Glaciers were −200 m a −1 and −56 m a −1 , respectively . Although both glaciers have shown acceleration between 2000 and 2014, Tracy Glacier (51 m a −2 ) accelerated almost four times faster than Heilprin Glacier (13 m a −2 ) . Contrasting patterns were also observed in the change in glacier thickness between 2011 and 2012, i.e., Tracy Glacier showed a thinning rate double (−12 m a −1 ) that of Heilprin Glacier (−6 m a −1 ) (Porter et al., 2014). More recently, from 2016 to 2017, Tracy Glacier thinned by 9.9 m near the terminus while the thinning at Heilprin Glacier was only 1.9 m (Willis et al., 2018). Previous studies attribute the more rapid retreat of Tracy Glacier to the difference in the depth of the fjord. Tracy Glacier terminates in deeper water than Heilprin Glacier (>600 m at Tracy Glacier and ∼350 m at Heilprin Glacier), therefore the ice is plausibly more exposed to warmer deep water and subjected to rapid submarine melting (Porter et al., 2014(Porter et al., , 2018Willis et al., 2018). Bowdoin Glacier is a 3-km-wide glacier terminating in Bowdoin Fjord, a 20-km-long fjord connected to Inglefield Bredning ( Figure 1). This glacier thinned at a rate of −4.1 m a −1 between 2007 and 2010, which was roughly 50% faster than the adjacent land-terminating Tugto Glacier (−2.8 m a −1 ) (Sugiyama et al., 2015;Tsutaki et al., 2016). This result highlights the importance of ice dynamics and ice-ocean interactions in the recent thinning of marine-terminating glaciers in this region.

AeroDEM
AeroDEM is a product derived from aerial photographs taken between 1978 and 1987 (Korsgaard et al., 2016). This DEM covers the entire margins and surrounding region of the Greenland Ice Sheet with a spatial resolution of 25 m. The accuracies in horizontal and vertical directions are reported as <10 and <6 m, respectively, while the precision is greater than 4 m (Korsgaard et al., 2016). This product has been used in previous studies for evaluating mass loss along the margins of the Greenland Ice Sheet (Felikson et al., 2017;Khan et al., 2014;Kjaer et al., 2012;Kjeldsen et al., 2015;Mouginot et al., 2019) and peripheral glaciers (Abermann et al., 2020;Huber et al., 2020;Marcer et al., 2017;von Albedyll et al., 2018). The part of AeroDEM used in this research was generated from aerial photographs taken in 1985.

ASTER DEMs
The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) is an imaging instrument, developed by Japan and the USA, onboard the Terra satellite launched in December 1999. Advanced Industrial Science and Technology in Japan distributes value-added ASTER data (ASTER-VA), which contain DEMs derived from stereo pair ASTER images (https://gbank.gsj.jp/madas/). Each ASTER-VA DEM covers an area of 60 × 60 km 2 , with a horizontal resolution of 30 m (Fujisada et al., 2005;Hirano et al., 2003). Fujisada et al. (2005) reported that the vertical accuracy of ASTER-VA DEM is 20 m with 95% confidence. We used 13 ASTER-VA DEMs between 2001-2003 and 2016-2018. To select DEMs acquired late in the ablation season with minimal influence of cloud and snow cover, corresponding ASTER visible and near-infrared nadir band (vnir3n) images were inspected. Dates and IDs of the utilized DEM are summarized in Table 1.

GIMP DEM
Greenland Ice Mapping Project (GIMP) DEM was derived from high resolution panchromatic stereoscopic imagery collected between 2009 and 2015 by GeoEye-1 and WorldView 1-3 . This DEM, with a spatial resolution of 30 m, is available through the National Snow and Ice Data Center  (https://nsidc.org/data/NSIDC-0715/versions/1). In this study, we used GIMP DEM as a reference for co-registration of other DEMs as described in the Method section (Section 3.2.1).

Glacier Mask
Surface elevation changes were evaluated over glacier areas defined by the GIMP Land Ice and Ocean Classification Mask . The glacier boundaries were revised in this study by inspecting ASTER vnir3n images taken in the late ablation season in 2010. DEM differencing was performed over areas below 850 m a.s.l., which are covered by the AeroDEM. Glacial basin boundaries were defined by a previous study (Mouginot et al., 2019) and glacier areas were subdivided into elevation bands based on the GIMP DEM. The total glacier area analyzed in this study was 1,624.54 km 2 , as measured in 2010.

Meteorological and Oceanic Data
Atmospheric temperature recorded at a meteorological station in Qaanaaq (77.48°N, 69.38°W; 16 m a.s.l.), the nearest station to the studied glaciers, is available for the period from 1996 to the present. To extend the investigation into an earlier period, we also used daily mean air temperatures from 1985-2018 at Thule Airbase, (76.53°N, 68.75°W; 77 m a.s.l.) situated 108 km south of Qaanaaq. These temperature data are accessible via the website of the National Oceanic and Atmospheric Administration. The mean summer temperatures (June, July, and August) at Thule Airbase are well correlated with those at Qaanaaq (r 2 = 0.83, p < 0.001). Thus, we assumed that the atmospheric temperature at Thule Airbase and Qaanaaq had a similar variation pattern (Figure 7a), which enabled us to investigate temperature change over the entire study period. In addition to this data obtained at low elevation areas, we analyzed the temperature measured at SIGMA-B, (77.54°N, 69.07°W; 944 m a.s.l.) located on Qaanaaq Ice Cap, where an automatic weather station has been operated since 2012 (Aoki et al., 2014).
We utilized monthly mean ocean temperatures during 1991-2019 from TOPAZ4 Arctic Ocean reanalysis data (Sakov et al., 2012) distributed by the Copernicus Marine Environment Monitoring Service. As a measure of thermal conditions of Polar Surface Water (Bevan et al., 2019), we extracted 5 m potential temperature from TOPAZ4 data and calculated summer mean values (June-August) within an area of 77-78°N and 66-72°W. To assess possible warming in subsurface Atlantic water, we adopted a reconstruction of ocean thermal forcing (depth-averaged temperature above freezing point) from 1992 to 2017 near Inglefield Bredning (76.8-77.3°N, 71-73°W) . The data set of thermal forcing is widely used to study the warming of Atlantic water and its role in the recent glacier retreat in Greenland Wood et al., 2018Wood et al., , 2021.

Corrections of DEM Biases
To derive surface elevation changes from multi-temporal DEMs, co-registration of the DEMs is necessary to remove potential offsets in horizontal and vertical directions (Nuth & Kääb, 2011). All DEMs were projected to WGS 1984 UTM zone 19N before co-registration. AeroDEM was resampled to a resolution of 30 m by cubic convolution to obtain the same cell size as ASTER and GIMP DEMs. Subsequently, we sampled all DEMs on to the same grid, based on the GIMP DEM.
We employed a three-dimensional co-registration method, as proposed by Nuth and Kääb (2011). The co-registration was performed on terrain with slopes smaller than 45° (Berthier et al., 2016) because the performance of DEMs derived from optical stereo-pair images on steep terrain is poor (Toutin, 2002). Horizontal and vertical offsets of AeroDEM were corrected using GIMP DEM as a reference, excluding ice, snow, water surfaces, steep terrain (>45°), and pixels where elevation differs from the reference DEM by more than 50 m (Berthier et al., 2016). As for the ASTER DEMs, along/cross track biases introduced by satellite acquisition geometry were corrected, based on the method proposed by Nuth and Kääb (2011), which was followed by horizontal and elevation-dependent corrections (Gardelle et al., 2013;Nuth & Kääb, 2011). Vertical offsets of the DEMs over the ice-free areas before and after the corrections are shown in Supporting Information S1 ( Figure S1). Outliers in the data were excluded, based on previously published elevation change records within the research region. The greatest glacier elevation change reported in our study area was −9.9 m a −1 in Tracy Glacier between 2016 and 2017 (Willis et al., 2018). We used these previous estimates as a guide, but set a more conservative threshold because the study area and period are not fully consistent with our research. We also expect surges and other short-term rapid glacier changes, which may not be captured if the tolerance is small. Thus, we excluded pixel values as outliers when the difference relative to the GIMP DEM fell outside the range −20 (twice the previously reported largest thinning rate) to +10 m a −1 (set arbitrarily to capture possible surges). Over the glacier mask analyzed in this study, only 1.5% of the pixels were excluded. The outliers were carefully inspected so that large elevation changes outside of this range due to glacier dynamics such as surging were also included in the analysis.

DEM Differencing
We separated the DEMs into three timeframes: (a) the 1985 Aerial DEM, (b) ASTER DEMs acquired from 2001-2003, and (c) 2016-2018, hereafter referred to as T0, T1, and T2, respectively. Elevation change (Δh) and its mean rate (Δh/Δt) were obtained by differencing DEMs for three periods (T0-T2, T0-T1, and T1-T2). All available Δh/Δt values were mosaicked to produce a Δh/Δt map for each period, covering all 16 glacier surfaces below 850 m a.s.l. As a measure of the uncertainty in the elevation change, Δh/Δt on ice-free areas near the glaciers were analyzed in each period and presented as histograms in the Supporting Information S1 ( Figure S2).
Elevation changes were analyzed for every 50 m altitude bin, based on the GIMP DEM. Within each elevation band, we calculated a mean Δh/Δt after excluding values deviating from the average by more than twice the standard deviation (Berthier et al., 2004;Gardelle et al., 2013), with 5% of pixels excluded in this process. Because elevation change was unavailable near the fronts of rapidly retreating glaciers, we analyzed the lowest elevation bin (0-50 m) only if available data covered more than 10% of the area. Banded structures were observed in the elevation change on Heilprin and Tracy Glaciers (Figure 2), which we attribute to satellite jitter (Girod et al., 2017) and the signals being within the range of uncertainty.

Uncertainty Analysis
To estimate the uncertainty in the elevation change, we considered the spatial correlation of the Δh (Gardelle et al., 2013). The estimator was the standard error of the mean (e) defined as where E Δh/Δt represents the standard deviation of the elevation change rate and N eff represents the effective number of the observations on the ice-free area. Considering the spatial correlation (Rolstad et al., 2009), the N eff is given by the total number of observations (N tot ) where R is the pixel size (30 m) and d is the spatial autocorrelation distance (713 m), determined by the mean of Moran's I autocorrelation index on elevation differences in the ice-free area (Gardelle et al., 2013). Uncertainties within Inglefield, Baffin Bay and for all of the studied glaciers (e region ) were computed by weighting uncertainties obtained in individual glaciers included in each region (e i ) for corresponding glacier areas (A i ).

Elevation Change Over the Study Area
Across most of the glacier areas, surface elevation showed a significant decrease during the entire study period from T0 to T2 (Figure 2a). The mean rate of surface elevation change over the glaciers studied during T0-T2 was −0.55 ± 0.24 m a −1 . All the elevation bands (0-850 m) experienced surface lowering, and the magnitude of the rate decreased in higher elevation (Figure 3a). The most rapid change (−3.08 m a −1 ) was observed near the glacier termini (elevation band 0-50 m), and the smallest change (−0.14 m a −1 ) was found in the highest elevation band 800-850 m.
Contrasting elevation change patterns were observed for the two subperiods, that is, surface lowering accelerated from T0-T1 to T1-T2 over the entire elevation range (Figure 3b). Elevation change was nearly zero or slightly positive from T0 to T1 (Figure 2b), while strongly negative elevation change was observed from T1 to T2 (Figure 2c). The mean elevation change rate for T0-T1 was 0.14 ± 0.17 m a −1 with an upglacier increasing trend (Figure 3b). In contrast to the first period, clearly negative elevation change was detected during the second subperiod (T1-T2). The mean rate over the study area was −1.31 ± 0.20 m a −1 and the most rapid change was observed near the frontal areas (−5.47 m a −1 at 0-50 m) (Figure 3b). Even in the elevation band 800-850 m, the glaciers experienced surface lowering at a rate of −0.57 m a −1 .
Surface lowering was slightly enhanced at the elevation band 550-600 m (Figure 3). This trend is due to the altitudinal distribution of Tracy Glacier, the most rapidly thinning and the second largest among the studied glaciers. The surface area of this glacier accounts for 8.4% of the study area at 550-600 m, while it only accounts for 3.5% and 7.9% in the adjacent 500-550 and 600-650 m bands, respectively. Therefore, rapid

Individual Glaciers
The elevation change showed significant variations among the glaciers ( Table 2). Most of the glaciers exhibited no obvious change or only a slight thickening during the period T0-T1. Obvious surface lowering was observed only in the frontal areas of Tracy, Farquhar, Sharp, and Sun Glaciers (Figure 2b). The most significant surface lowering during this period was found below 600 m a.s.l. in Tracy and Farquhar Glaciers (Figures 4b and 4c), which formed a common glacier tongue until it disintegrated in 2002 . Nevertheless, even in these regions, the magnitude of thinning in the period T0-T1 was less than 2 m a −1 . Spatial patterns of the elevation change in T1-T2 are significantly different from those in the previous period. All the glaciers experienced surface lowering at a rate substantially different for each glacier (Figure 2c and Table 2). The thinning was particularly rapid at Tracy and Farquhar Glaciers, where the mean rates of the elevation change were −3.91 ± 0.13 and −2.91 ± 0.18 m a −1 , respectively. The rates near the calving front of these two glaciers reached −9 m a −1 (Figures 4b and 4c). Heilprin Glacier thinned at a rate of −0.51 ± 0.13 m a −1 , which was ∼60% slower than the mean of the glaciers studied. The magnitude increased toward the calving front, but did not exceed −1.5 m a −1 (Figure 4a). Thinning rates of Morris Jesup and Verhoeff Glaciers (−0.08 ± 0.35 and −0.14 ± 0.14 m a −1 ) were the two lowest among the glaciers studied during T1-T2 (Table 2). Interestingly, Morris Jesup Glacier showed contrasting altitudinal distribution patterns between the two subperiods ( Figure 5d). During T0-T1, the glacier thickened over the area below 550 m a.s.l., whereas thinning was observed above 550 m a.s.l. The thinning was more significant at a higher elevation, being enhanced from −0.27 m a −1 at 500-550 m a.s.l. to −1.12 m a −1 at 800-850 m a.s.l. During T1-T2, the glacier thinned and thickened below and above 600 m a.s.l., respectively. No clear altitudinal variation was found at Verhoeff Glacier during T0-T1 and T1-T2 (Figure 5b), where elevation change over the last three decades was the smallest among the glaciers studied. The only land-terminating glacier, Tugto Glacier, showed a relatively small elevation change at rates of 0.23 ± 0.14 m a −1 for T0-T1 and −0.83 ± 0.14 m a −1 for T1-T2. The latter corresponds to approximately 60% of the mean of the glaciers studied. No apparent dependence on altitude was observed in either subperiod for this glacier (Figure 5h).

Regional Variation
The magnitude of thinning and its temporal patterns vary between the Baffin Bay and Inglefield Bredning regions (Figure 2). To investigate regional patterns, thinning rates of glaciers terminating in Inglefield Bredning (Heilprin, Farquhar, Melville, Sharp, Hart, Hubbard, and Bowdoin) are compared with those terminating in Baffin Bay (Sun, Verhoeff, Meehan, Morris Jesup, Diebitsch, Clements, and Bamse) ( Figure 6). Tracy Glacier was excluded from this analysis to avoid the influence of its much higher thinning rate.   In the period T1-T2, elevation changes were negative over the entire elevation range in the two regions ( Figure 6b). The magnitude of thinning was ∼70% greater in the Inglefield Bredning region (−1.07 m a −1 ) than in the Baffin Bay region (−0.64 m a −1 ). The thinning was most accelerated near the front of the glaciers in the Inglefield Bredning region, which is represented by the change in the 0-50 m a.s.l. band from −0.31 m a −1 (T0-T1) to −3.69 m a −1 (T1-T2) ( Figure 6). In the Baffin Bay region, the most significant change from 0.36 m a −1 (T0-T1) to −1.54 m a −1 (T1-T2) occurred at 100-150 m a.s.l.

Atmosphere and Ocean Temperature
Summer

Comparison With Previous Studies
Here, we compare our results with previously published estimates for glacier mass loss and elevation change in Greenland and adjacent areas. We converted the observed volume change to mass change using the ice density of 910 kg m −3 . The mass loss over the studied areas from T0-T2 was −0.81 ± 0.44 Gt a −1 , which corresponds to 0.5% of the total mass loss of the Greenland Ice Sheet (−150 Gt a −1 ) from 1992 to 2018 (Shepherd et al., 2020). After 2001, the mass loss increased to −1.94 ± 0.15 Gt a −1 during T1-T2, accounting for 0.8% of the Greenland Ice Sheet mass change for 2005-2015 (Shepherd et al., 2020). On the western side of Baffin Bay, thinning rates of glaciers and ice caps on the Queen Elizabeth Islands during 2005/06-2012/14 were more than three times greater than those during 1995-2000. In this region, marine-terminating glaciers showed greater acceleration than those terminating on land, suggesting that an increase in ocean heat flux plays an important role in addition to atmospheric warming (Mortimer et al., 2018).
Elevation change of Tracy Glacier from 2003 to 2007 was reported as −5.46 m a −1 at 600 m a.s.l. (Pritchard et al., 2009), which is higher than the rate we obtained for T1-T2 (−3.77 m a −1 at 550-600 m a.s.l.). The difference is potentially due to the shorter sampling period or the limited number of ICESat sampling sites analyzed in the previous study. More recent research showed thinning within 6 km from the front (approximately 0-400 m a.s.l.) at a rate of −9.9 m a −1 from 2016-2017, which is ∼30% greater than our result in Note. The area of each glacier and glacier group is also shown in the table.

Table 2 Mean Elevation Change (m a −1 ) During the Three Periods (T0-T2, T0-T1, and T1-T2) on Each Glacier as Well as Glacier Groups in Inglefield Bredning (Numbers in Brackets Indicate the Results Calculated Without Tracy Glacier), Baffin Bay Region, and All Glaciers Studied
T1-T2 (−7.40 m a −1 ) (Willis et al., 2018). Porter et al. (2014) reported that the thinning rate in the ablation zone doubled from −6 m a −1 in 2002-2010 to −12 m a −1 in 2011-2012 (Porter et al., 2014). Therefore, mass loss of Tracy Glacier appears to have increased in the 21st century.

Driving Mechanism of the Elevation Change
In general, thinning of the glaciers accelerated from T0-T1 to T1-T2. However, the timing and magnitude of the acceleration showed substantial variation from glacier to glacier. Here, we discuss air temperature, ocean temperature, ice dynamics, and fjord bathymetry as possible drivers of the spatiotemporal variability in the thinning rate.

Atmospheric Warming
A substantial rise in air temperature has been observed in Greenland since the mid-1990s (Box et al., 2009;Mernild et al., 2011). The warming atmospheric conditions have enhanced surface melting (van den Broeke et al., 2009), resulting in a mass loss of −76 Gt a −1 over the Greenland Ice Sheet from 1992 to 2008 (Shepherd et al., 2020). Glacier thinning in our study area increased from 0.14 ± 0.17 m a −1 in T0-T1 to −1.31 ± 0.20 m a −1 in T1-T2. Temperatures at Thule Airbase and in Qaanaaq showed warming trends during these periods (Figure 7a), suggesting the influence of atmospheric warming on the observed glacier changes. Temperature sensitivity of the mass loss is calculated as −0.61 m a −1 K −1 from the elevation change rate from T0-T2 and temperature increase during the same period. This number is slightly greater than −0.48 m a −1 K −1 , as reported in the Canadian Arctic Archipelago (Gardner et al., 2011). The sensitivity is more than double those estimated for Arctic glaciers and ice caps, based on only surface mass balance (Woul & Hock, 2005), indicating that the observed mass loss was not entirely due to atmospheric warming.
In addition to the long-term warming trend, unprecedented melt events due to exceptional atmospheric circulation patterns have been reported in Greenland (e.g., Hanna et al., 2014;Nghiem et al., 2012;Tedesco & Fettweis, 2020). The large amount of meltwater generated during melt events either forms supraglacial lakes and streams or is delivered into the ice-bed interface through moulins and crevasses (Chu, 2014). The supraglacial lakes and streams enhance ablation because of their lower albedo relative to the surrounding ice . Moreover, the increased input of meltwater into an inefficient subglacial hydrological system elevates basal water pressure and enhances basal sliding, resulting in glacier acceleration and dynamic thinning (Bartholomew et al., 2010;Zwally et al., 2002). Furthermore, the increased volume of subglacial discharge potentially enhances submarine melting because the discharge activates fjord circulations and facilitates more efficient oceanic heat transport to the glacier (Motyka et al., 2013;.

Oceanic Forcing
Oceanic forcing is recognized as a key control on Greenlandic outlet glaciers. Ocean warming drives not only glacier retreat (Fahrner et al., 2021;Wood et al., 2021), but also glacier acceleration (Holland et al., 2008;Howat et al., 2008) and thinning (Felikson et al., 2017;Thomas et al., 2009) through ice-ocean interaction at the glacier front. Intrusion of warming Atlantic waters is considered to be the driver of recently enhanced submarine melting of Greenland's outlet glaciers along the western coast . More rapid retreat, acceleration, and thinning were observed at glaciers located in deep fjords, which are more exposed to deep water warming (Rignot et al., 2010;Wood et al., 2018Wood et al., , 2021. Wood et al. (2021) investigated the influence of recent ocean warming on submarine melting at some of the glaciers in our study site. During the ∼2°C increase in the ocean thermal forcing from 1998 to 2007 (Figure 7b), an approximately fourfold increase was estimated for the upper bound of the melt rates at Tracy, Heilprin, and Bowdoin Glaciers (see Table S1 in Wood et al., 2021). These glaciers terminate in deep fjords in Inglefield Bredning (150-400 m deep), whereas no significant change was reported for those terminating in shallower fjords in the Baffin Bay region (e.g., Diebitsch, Verhoeff, and Morris Jesup glaciers) . Greater influence of ocean warming is a possible interpretation for the regional variations in the glacier thinning rate, since enhanced submarine melting causes persistent glacier retreat, acceleration, and dynamic thinning as observed at the glaciers in Inglefield Bredning .
In addition to warming in the deep layer, the TOPAZ4 data showed an increase in near-surface water temperatures during the period 1996-2012 (Figure 7b), suggesting its influence on ice front melting and calving. Near-surface ocean warming enhances melting at or below the waterline, resulting in an unstable ice front and an increase in calving rates (Benn et al., 2007). In some cases, enhanced calving can trigger initial retreat of a glacier from a bedrock bump, which leads to the destabilization of the glacier terminus, ice-flow acceleration, and dynamic thinning (Felikson et al., 2017;Porter et al., 2018).

Ice Dynamics
Recently observed rapid thinning of Greenlandic marine-terminating glaciers cannot be attributed to increasing surface melt alone. Rather, a significant portion of the thinning is due to ice dynamics driven by atmospheric and oceanic forcing. Dynamic thinning associated with ice acceleration has been observed extensively in the coastal areas of Greenland (e.g., Abdalati et al., 2001;Krabill et al., 2004;Thomas et al., 2009), and most prominently in the southeastern and northwestern regions (Csatho et al., 2014;Pritchard et al., 2009).
The magnitude of the glacier thinning obtained in this study is substantially greater than those reported from peripheral glaciers and ice caps in the region, most of which are terminating on land. Elevation changes of glaciers and ice caps in northwestern Greenland were reported to be −0.54 m a −1 in the period (Gardner et al., 2013 and −0.60 m a −1 in the period (Bolch et al., 2013. The thinning rate of six ice caps in the Qaanaaq region was −1.1 ± 0.1 m a −1 for the period 2006-2010 (Saito et al., 2016). These thinning rates are less than the mean rate obtained in this study for the period 2001-2018 (−1.31 ± 0.20 m a −1 ), indicating the significance of dynamic thinning at marine-terminating glaciers in northwestern Greenland.
To investigate a link between the thinning, acceleration, and retreat in the Prudhoe Land region, we used ice speeds and frontal displacement rates of the studied glaciers reported for the period between 2000 and 2014 . As represented by Tracy Glacier, rapidly thinning glaciers are characterized by greater retreat and acceleration (Figure 8). The correlation between the elevation change and acceleration is significant (r 2 = 0.70, p < 0.001) (Figure 8a). A relationship is still evident even if Tracy Glacier is excluded from the analysis, although the correlation coefficient decreases (r 2 = 0.39, p < 0.05). Correlation is also observed between the elevation change and retreat rate (r 2 = 0.57, p < 0.005), but it is insignificant when Tracy Glacier is excluded (r 2 = 0.12, p = 0.22) (Figure 8b). Observations at Bowdoin Glacier suggested such a link between glacier mass loss and ice dynamics. Sugiyama et al. (2015) observed a two-fold acceleration near the glacier front from 1999 to 2003, which coincided with a small glacier retreat. The acceleration interrupted relatively stable glacier conditions since the 1980s, leading to rapid thinning and retreat since 2008. The thinning rate of Bowdoin Glacier from 2007 to 2010 was approximately 50% greater than that of the adjacent land-terminating Tugto Glacier. The rapid thinning was thus attributed to longitudinal stretching caused by acceleration . Similar connections between acceleration and thinning were also suggested for other glaciers in Inglefield Bredning, i.e., Heilprin, Tracy, Farquhar, and Diebitsch Glaciers, based on data obtained between 2000 and 2014 . Thinning and acceleration are often associated with glacier front retreat. Our data show that the magnitude of acceleration has a stronger correlation with surface elevation change rate than does the magnitude of retreat (Figure 8). Therefore, ice acceleration is more strongly linked with the thinning rates of the glaciers in the Prudhoe Land region during the period T1-T2.
The intriguing elevation change pattern observed at Morris Jesup Glacier is another example of glacier change driven by ice dynamics. Mouginot et al. (2019) reported a more than two-fold increase in the solid ice discharge of this glacier from 1975 to 1990 (from 0.15 to 0.35 Gt a −1 ), which corresponds to thinning (thickening) in the upper (lower) reaches during T0-T1 ( Figure 5d). Presumably, rapid ice transport resulted in the thickening downglacier, which was followed by a 5 m a −2 glacier deceleration  and thickening (thinning) in the upper (lower) reaches during T1-T2 (Figure 5d). This observation indicates a long-lived surge event similar to that which has been documented at Hagen Brae (Solgaard et al., 2020) and Storstrømmen .
Glaciers terminating in relatively narrow fjords (1.2-2.5 km wide), including Bamse, Cements, Meehan, Sun, Hubbard, and Hart Glaciers, showed a stable or slightly decreasing trend in the ice speed ( Figure 8a) , while our data indicate ice thinning during T1-T2 (Figures 4 and 5). It is likely that ice dynamics of these glaciers were relatively insensitive to frontal retreat and decrease in basal drag after thinning because driving stress is more supported by lateral drag (Cuffey and Paterson., 2010). These glaciers decelerated under a greater influence of reduction in the driving stress, while the glaciers thinned by increasingly negative mass balance rather than changes in glacier dynamics.

Fjord Bathymetry
Fjord bathymetry is critically important for ice dynamics near the glacier front (e.g., Enderlin et al., 2013;. It also affects frontal ablation because fjord depth controls the access of the relatively warm deep ocean water to the ice front (Holland et al., 2008;Porter et al., 2014). Among the glaciers studied, the changes observed in the neighboring Tracy and Heilprin Glaciers have varied greatly in terms of their front position, flow speed, and surface elevation since the 2000s (Moon et al., 2012;Rignot & Kanagaratnam, 2006). Tracy Glacier has experienced more rapid retreat, speedup, and surface lowering than Heilprin Glacier (Porter et al., 2014;Pritchard et al., 2009;Willis et al., 2018). Porter et al. (2014) suggested that Tracy Glacier is more vulnerable to changes in the ocean because of its deeper grounding line (>600 m), compared to Heilprin Glacier (∼350 m), causing its ice front to be more exposed to warm water at depth. To assess the influence of fjord bathymetry on glacier thinning, we investigated the basal topography from the BedMachine v3 product   (Figure 9). The fjords on the Baffin Bay side are shallower than those belonging to Inglefield Bredning, which generally exceed several hundred meters in depth. Ocean temperature data from Bowdoin Fjord indicated that a depth below ∼200 m is occupied by relatively warm Atlantic water (Kanna et al., 2018;Ohashi et al., 2020). Therefore, more rapid glacier thinning in the Inglefield Bredning region is most likely affected by greater exposure to oceanic heat.
This situation contrasts to that on the western side of Baffin Bay. Cook et al. (2019) studied recent frontal variations of marine-terminating glaciers in the Canadian Arctic Archipelago, namely Queen Elizabeth, Baffin, and Bylot Islands. They concluded that elevated atmospheric temperature was the primary driver of glacier change, and that the influence of deep ocean water temperature was insignificant because the fjords are shallower than those along the western coast of Greenland and deep ocean water was not able to access the glacier front (Cook et al., 2019).
In addition to the role of fjord depth in ocean heat transfer, glaciers terminating in deep fjords are more subject to the influence of buoyancy when ice thins close to flotation. Although the glaciers along Inglefield Bredning are grounded, Heilprin, Tracy, Farquhar, and Bowdoin Glaciers are near flotation condition Sugiyama et al., 2015). Therefore, even a small perturbation in ice thickness or front position has a significant impact on the stress balance and stability of the terminus region (Benn et al., 2007). Due to the ice's exposure to heat in the deep ocean, as well as the greater susceptibility in the stress balance, we hypothesize that deep glacier bed geometries have an effect on the rapid thinning that has been observed.

Conclusions
In this study, we measured the glacier surface elevation change of 16 outlet glaciers along the coast of Prudhoe Land in northwestern Greenland. The measurements were taken from 1985 to 2018 using AS-TER-VA DEM and the recently released AeroDEM. All the glaciers studied experienced surface lowering between 1985 and 2018 at a mean rate of −0.55 ± 0.24 m a −1 . The elevation change greatly accelerated after the 2000s. We detected a slight thickening (0.14 ± 0.17 m a −1 ) prior to 2001, whereas substantial thinning (−1.31 ± 0.20 m a −1 ) was observed during the period 2001-2018. Of the glaciers studied, Tracy and Farquhar Glaciers located in Inglefield Bredning thinned most rapidly, at a rate exceeding −9 m a −1 in the period 2001-2018. Summer air temperatures have shown a warming trend from 1996 to 2019, indicating that enhanced surface melting is an important driver of glacier mass loss. Nevertheless, the acceleration of thinning that has been observed cannot be attributed to atmospheric warming alone. Deep ocean temperatures showed a warming trend between 1998 and 2007, which was potentially the driving force behind the rapid thinning and retreat of the glaciers located in deep fjords after the year 2001. Relatively deep fjords increase the exposure of the ice front to deep warm ocean water, which might have amplified the impact of deep ocean warming on the glaciers terminating in Inglefield Bredning. Near-surface ocean temperature showed a similar warming trend to that of the atmosphere, suggesting a possible influence on the glacier change. There was a significant correlation between the rate of change in glacier elevation and flow acceleration, suggesting the importance of dynamic thinning in the surface lowering of the glaciers studied. We also identified signs of a surge event on the elevation change pattern of Morris Jesup Glacier. Therefore, we conclude that the recent thinning of the glaciers along the Prudhoe Land region was generally controlled by atmospheric and oceanic conditions. In addition to the general trend, the response of individual glaciers to the atmospheric/ oceanic forcing was modified by ice dynamics and fjord bathymetry.