Change Points Detected in Decadal and Seasonal Trends of Outlet Glacier Terminus Positions across West Greenland

: We investigated the change in terminus position between 1985 and 2015 of 17 marine-terminating glaciers that drain into Disko and Uummannaq Bays, West Greenland, by manually digitizing over 5000 individual frontal positions from over 1200 Landsat images. We ﬁnd that 15 of 17 glacier termini retreated over the study period, with ~80% of this retreat occurring since 2000. Increased frequency of Landsat observations since 2000 allowed for further investigation of the seasonal variability in terminus position. We identiﬁed 10 actively retreating glaciers based on a signiﬁcant positive relationship between glaciers with cumulative retreat > 300 m since 2000 and their average annual amplitude (seasonal range) in terminus position. Finally, using the Detecting Breakpoints and Estimating Segments in Trend (DBEST) program, we investigated whether the 2000–2015 trends in terminus position were explained by the occurrence of change points (signiﬁcant trend transitions). Based on the change point analysis, we found that nine of 10 glaciers identiﬁed as actively retreating also underwent two or three periods of change, during which their terminus positions were characterized by increases in cumulative retreat. Previous literature suggests potential relationships between our identiﬁed change dates with anomalous ocean conditions, such as low sea ice concentration and high sea surface temperatures, and our change durations with individual fjord geometry.


Introduction
Over recent decades, the Greenland Ice Sheet has undergone increasing rates of mass loss, a large component of which is the result of retreat and speed up of marine-terminating outlet glaciers that have been subject to both increasing atmospheric and oceanic temperatures. Previous work has indicated high variability in the duration and magnitude of outlet glacier terminus responses to environmental perturbations [1][2][3][4]. Specifically, studies have found a rapid onset of terminus retreat in response to warming air temperatures, yet a much slower advance response to a cooling climate [5,6]. In addition to terminus retreat, recent observations of other coupled glacier dynamics, such as thinning [7], velocity [8,9], and calving, are also changing in response to long-term climate trends. Challenges arise in determining a primary control on glacier terminus change, as the magnitude of retreat observed at the terminus is determined by not only the atmospheric and oceanic temperatures, but also physical constraints unique to each glacier, such as ice mélange buttressing and glacial runoff [10], underlying bedrock topography [4], and fjord bathymetry [11]. These factors unique to each glacier have often led to asynchronous terminus retreat (or advance) trends of glaciers in the same region in response to climate changes [2,3]. This temporally non-linear glacier behavior is seen at glaciers throughout the world, for example, marine-terminating glaciers in Novaya Zemlya underwent dramatic retreat between 2000 and 2013, corresponding with high air temperatures and low sea ice concentrations, yet subsequently retreat rates decreased between 2013 and 2015 [12].
Monitoring terminus change is one of the simplest ways to track individual outlet glacier behavior and can inform the relationship between changes at the calving front, subglacial hydrology, and the inland propagation of such disturbances on glacier flow [5]. Many studies have utilized remotely sensed imagery to monitor outlet glacier terminus change of glaciers in West Greenland, over different spatial and temporal resolutions [1,2,4,6,9,[13][14][15][16]. Previous remote sensing studies across different regions of Greenland have shown spatial variability in seasonal and interannual glacier terminus behavior [13,15], indicating no simple linkage to air or sea temperatures, but rather, a strong control by geometric factors, such as the seasonal range in terminus position being proportional to individual glacier width [13]. In Uummannaq Bay, specifically, the spring ice mélange clearing date, which has also been linked to fjord geometry [11], was correlated with seasonal terminus variability [1]. While the magnitude and duration of individual terminus retreat may largely be controlled by geometric factors, such as overdeepening behind the terminus, onset of retreat at glaciers across central West Greenland was near synchronous beginning in 1998 ± 3 [4].
A tipping point in a system is defined as a non-linear response to a gradual forcing [17]. As such, the response of terminus retreat (advance) of outlet glaciers to global temperature change provides an example of an environmental system that may be vulnerable to reaching a tipping point. Specifically, modelling studies suggest that the slow response rate of marine-terminating glaciers to recent atmospheric and oceanic forcings points to not only a continuation of current glacier change observations, but also a commitment to even greater changes [18]. Understanding both the short-term (seasonal) and long-term (decadal) glacier terminus variability prior to past transitions in their behavioral trends can help elucidate how current terminus retreat (advance) responds to recent atmospheric and oceanic conditions and how this response may be indicative of approaching a tipping point in the future. Previous work found that critical slowing down occurred prior to the crossing of a tipping point in the paleoclimate system [19]. In particular, this critical slowing down, or rate of rebound to a previous state, was measured by decreases in the rate of change within the system coincident with increases in short-term autocorrelation [19]. Previous work has investigated tipping points in sea ice [20,21], paleoclimate records [19], the Amazon rainforest [22], and Atlantic thermohaline circulation [22]. However, direct application of these types of analyses applied to glacier terminus front behavior are currently limited.
Recent outlet glacier terminus studies with higher temporal resolution datasets have begun to broach this topic. For example, Catania et al. [4] identified the onset of "progressive terminus retreat" by requiring that terminus retreat quantities within a given summer be greater than the previous two summers. Yet, those authors also point to difficulties in this method for determining retreat onset timing at glaciers with large seasonal variability but small decadal retreat [4]. Other statistical techniques have proven robust for identification of change points based on significant breaks in time series of terminus positions [12,23,24].
To systematically address the types of challenges that arise in identifying the occurrence of change points in trends, in this study we apply a method for identification of change points in outlet glacier terminus positions developed from remotely sensed imagery that uses a segmentation algorithm to determine gradual and abrupt trend transitions. We first present a seasonal and decadal inventory of terminus positions for outlet glaciers in Disko and Uummannaq Bays, West Greenland, from 1985 to 2015. We utilized Landsat 5-8 imagery to investigate the relationship between the annual amplitude in glacier terminus position and the 31 year trends in cumulative terminus retreat (advance) of 17 marine-terminating glaciers. The temporal resolution of our terminus observations (14) per year since 2000 allowed us to test for the existence of change points in the trends of glacier front position using the Detecting Breakpoints and Estimating Segments in Trend (DBEST) program in the R software [25]. We then assessed the relationship between the annual amplitude in position and the terminus retreat (advance) behavior in the context of the identified change points in order to understand glacier fluctuation during the different phases.

Study Area
Disko and Uummannaq Bays are located in West Greenland ( Figure 1). Uummannaq Bay outlets into Baffin Bay and is bounded by the Svartenhuk Peninsula to the north, and Nuussuaq Peninsula to the south. Disko Bay is south of Nuussuaq Peninsula and opens into Davis Strait. Disko Bay is also home to Sermeq Kujalleq (Jakobshavn Isbrae), which is the fastest-moving glacier in the world and alone drains~7% of the ice sheet [8], while Uummannaq Bay includes Kangilliup Sermia (Rink Isbrae) and Sermeq Kujalleq (Store Glestcher), which are the sixth and fifth fastest-flowing glaciers in Greenland [26]. Disko and Uummannaq Bays include 19 marine-terminating outlet glaciers-eight and 11, respectively-which we will refer to by the official Greenlandic name [27] and/or a reference letter we defined for the purpose of this study (Table A1).
Recent field campaigns across the study region have revealed a wide range in glacial fjord bathymetry and ocean temperatures at depth, including influences from both Polar (<2.5 • C) and Atlantic (>2.5 • C) waters [11]. Sea ice in Baffin Bay and Davis Strait has been declining during recent decades and the increased heat flux associated with the subsequent increase in open water has been directly linked to enhanced ice sheet surface melt across West Greenland [28].

Data
Approximately 1200 Level 1T Landsat images (30 m; 16 day repeat cycle) from missions 5, 7 and 8, generally available from March through October, over the years 1985-2015, were downloaded from United States Geologic Survey Earth Explorer (http://earthexplorer.usgs.gov/). Level 1T data were previously radiometrically and geometrically corrected by the USGS. Landsat imagery was not selected based on cloud cover criteria, as the percentage of overall scene cloud cover may not be indicative of individual glacier terminus visibility. Imagery also did not require atmospheric correction because spectral information was not needed for terminus visibility. Gaps exist in the Landsat record during the 1990s when sensors were turned off for regions outside of the United States, as well as annually in winter months (DJF) owing to insufficient sunlight at high latitudes. We found Landsat image co-registration errors easily detectable when each glacier's box (see Section 2.3, below) was overlaid in a scene and the fewer than 20 images in which this occurred were not included in our analysis. An uncommon source of error was the failure of the scan line corrector on Landsat-7 beginning in 2003; however, the black banding associated with the scan line failure did not frequently interfere with terminus visibility, e.g., (1) the terminus fell in between two bands, leaving it fully visible; or (2) bands that fell perpendicular to the terminus did not substantially inhibit the visibility of a majority of the terminus [1,6,9,16]. Some glaciers lie at the edge of multiple Landsat path/row combinations and owing to scene overlap that increases towards higher latitudes, these glaciers have more frequent observations than the 16 day Landsat repeat cycle [29]. Therefore, the number of visible glacier termini varied within individual images and the total number of observations of each glacier terminus ranged from 234 to 352 over the 1985-2015 period (Table A1).

Retreat Calculation
All visible glacier terminus positions within each Landsat scene were digitized manually ( Figure  2) by a single user in order to limit potential biases produced by multiple digitizers [31]. User error in manual terminus digitization was estimated by repeat digitization of two glacier termini 10 times each within a winter Landsat-8 image. We chose a winter image with ice mélange present at the glacier fronts, as this should represent a near upper bound of potential errors associated with the manual digitization technique. Ultimately, we found digitized termini to be consistent within 20 m, less than the spatial resolution of the Landsat imagery. Glacier termini that were unable to be differentiated from ice mélange, typically in early spring images, were not digitized. This most frequently occurred at Sermeq Silardleq (f), leading to ~25% fewer observations (~150) of this glacier terminus relative to the average (>200) since 2000 (Table A1). Although infrequent, in cases of the Landsat-7 scan line failure (Section 2.2), we linearly interpolated across missing lines that crossed the

Retreat Calculation
All visible glacier terminus positions within each Landsat scene were digitized manually ( Figure 2) by a single user in order to limit potential biases produced by multiple digitizers [31]. User error in manual terminus digitization was estimated by repeat digitization of two glacier termini 10 times each within a winter Landsat-8 image. We chose a winter image with ice mélange present at the glacier fronts, as this should represent a near upper bound of potential errors associated with the manual digitization technique. Ultimately, we found digitized termini to be consistent within 20 m, less than the spatial resolution of the Landsat imagery. Glacier termini that were unable to be differentiated from ice mélange, typically in early spring images, were not digitized. This most frequently occurred at Sermeq Silardleq (f), leading to~25% fewer observations (~150) of this glacier terminus relative to the average (>200) since 2000 (Table A1). Although infrequent, in cases of the Landsat-7 scan line failure (Section 2.2), we linearly interpolated across missing lines that crossed the glacier terminus. Glaciers determined to have a substantial portion of their terminus covered by the Landsat-7 banding were not digitized. glacier terminus. Glaciers determined to have a substantial portion of their terminus covered by the Landsat-7 banding were not digitized. To calculate retreat (advance) over the 1985-2015 period, we utilized the box method [32]. The box method accounts for unequal change along the terminus and allows for the measured retreat (advance) of glaciers of varying sizes to be comparable [1,15,32,33]. With this method, approximate sides of the glacier/fjord and an arbitrary upflow reference line remain consistent over the study period, while the fourth side of the box varies as the digitized terminus for a given image date [30]. Retreat was then calculated by subtracting the change in the box area at subsequent time periods and dividing by the average glacier width. We calculated a cumulative retreat (advance) as a running sum from a prescribed zero assigned at the first terminus observation of each glacier in 1985. To have a study inclusive of all outlet glaciers draining into the bays, an alternative method, e.g., the method described in [33], to the box method would have been required in order to measure cumulative retreat (advance) of Salliarutsip Sermia (Ingia Isbrae) in Uummannaq Bay and Sermeq Kujalleq (Jakobshavn Isbrae) in Disko Bay and as such, they were not included in this study ( Figure 1).
The total number of glaciers included in this study was 17-seven that drain into Disko Bay and 10 into Uummannaq Bay (Table A1). With the abovementioned considerations, the total number of terminus observations for the 17 glaciers focused on in this study was greater than 5000, averaging 297 ± 35 (one standard deviation) per glacier over the entire 1985-2015 period. Approximately two- To calculate retreat (advance) over the 1985-2015 period, we utilized the box method [32]. The box method accounts for unequal change along the terminus and allows for the measured retreat (advance) of glaciers of varying sizes to be comparable [1,15,32,33]. With this method, approximate sides of the glacier/fjord and an arbitrary upflow reference line remain consistent over the study period, while the fourth side of the box varies as the digitized terminus for a given image date [30]. Retreat was then calculated by subtracting the change in the box area at subsequent time periods and dividing by the average glacier width. We calculated a cumulative retreat (advance) as a running sum from a prescribed zero assigned at the first terminus observation of each glacier in 1985. To have a study inclusive of all outlet glaciers draining into the bays, an alternative method, e.g., the method described in [33], to the box method would have been required in order to measure cumulative retreat (advance) of Salliarutsip Sermia (Ingia Isbrae) in Uummannaq Bay and Sermeq Kujalleq (Jakobshavn Isbrae) in Disko Bay and as such, they were not included in this study ( Figure 1).
The total number of glaciers included in this study was 17-seven that drain into Disko Bay and 10 into Uummannaq Bay (Table A1). With the abovementioned considerations, the total number of terminus observations for the 17 glaciers focused on in this study was greater than 5000, averaging 297 ± 35 (one standard deviation) per glacier over the entire 1985-2015 period. Approximately  (Table A1).

Change Point Analysis
To detect potential change points in our data, we used the change detection algorithm in the Detecting Breakpoints and Estimating Segments in Trend (DBEST) program version 1.8 in the R software package. This program uses a segmentation algorithm to detect changes in both long-term trends, as well as identify abrupt events [25].
Input data into DBEST can be with or without any cyclical (seasonal) component. Cyclical data require regularly spaced intervals in time in order to precisely extract an underlying trend from the seasonality. Therefore, we used an Akima spline method of interpolation on the Landsat observations for each glacier over the 2000-2015 period to create a daily time series. The Akima spline method gives better fits to curves with a rapidly changing second derivative and performs better than linear methods on coarse or unequally spaced time intervals, resulting in a lower maximum error [34]. We acknowledge that glacier termini are known to be highly variable, and the interpolation was not performed to assume daily glacier positions, but solely for the ability to ensure all of our observed data points were included in the change point analysis. As such, our largest temporal interpolation occurs over the winter season, between the last observation of a given year and the first observation of the following year. Interpolated data values should have little effect on the change point results because the change detection algorithm in DBEST is primarily based on local peak and valley points where there is a significant change in the direction of the underlying trend. We tested interpolation to a coarser monthly temporal resolution, but this resulted in dampening the peak and valley points of interest. To minimize our assumptions of glacier behavior through interpolation, we limited the change point analysis to the 2000-2015 period owing to the higher frequency of glacier terminus observations during this time interval (Table A1), which gave us greater confidence in the daily interpolation of our time series. Additionally, this 16 year period was when the majority of terminus position change was observed.
The settings for the change detection algorithm in the DBEST program (Table 1) include data type (cyclical or non-cyclical), seasonality (if cyclical, over what time period), first and second level shift thresholds (tests for abruptness of a change), duration (temporal period over which a change in the mean needs to persist to be considered significant), statistical significance level alpha (in this case, 0.05), distance threshold (statistical boundary between outlier and normal residuals of an obtained fit), and the change magnitude threshold (test for the existence of a large enough change). These DBEST thresholds and parameters were uniquely determined for each glacier terminus dataset, discussed in detail in Section 2.5. First, DBEST tests for the occurrence of discontinuities-in this case of glacier terminus position-by analyzing the absolute differences between consecutive data point pairs and comparing this to the first level shift threshold. If the difference is greater than the first level shift threshold, then it tests whether or not this difference caused a significant modification to the mean glacier position and persisted over the duration threshold. If the mean glacier position before and after this identified discontinuity is greater than the second level shift threshold, DBEST considers this a level shift point. DBEST then repeats this process for all data points, sorts them into descending order based on the absolute value of the glacier position difference, and tests if the spacing between a data point and an identified level shift point is at least the duration threshold. Using the Seasonal-Trend decomposition based on Loess (STL) [35], a given glacier time series is then decomposed into trend (low frequency variation), seasonal variation (variation at or near the seasonal frequency), and remainder (remaining variation) components [25] ( Figure A1) separately for each portion of the time series divided by the identified level shift points. The remainders (unexplained variance) are generally very small and not the focus of attention in DBEST whose main purpose is detecting changes in the trend component. The trend component of the time series is then segmented using a peak/valley detector function and a method which draws a straight line through detected peak/valley points and compares perpendicular distances to the non-peak and non-valley points between them with a distance-threshold (discussed more in the following paragraph and Section 2.5). If the distance is greater than this threshold, these points are added to the set of detected peak/valley points and level shift points, all of which are called turning points. Detected turning points are then fit to the glacier position trend using a piecewise linear modelling and those turning points that minimize the Bayesian Information Criterion (BIC) [36] are considered breakpoints. Here, we use the change detection algorithm of DBEST which tests if the trend has segments with variation lower than the threshold value set for the magnitude of a change (Table 1), and as such, identifies a final set of breakpoints that are larger than this magnitude. The results of the algorithm include the break date, or start and end time between breakpoints; the change duration, or the temporal period over which this change occurred; the change value, or the amount of change that occurred over this time that has to be larger than the change magnitude threshold; the change type, whether the change is abrupt or gradual; and the change significance, based on statistical alpha level. Because we used the change detection algorithm within DBEST, we refer to our final set of identified breakpoints in our results as "change points." Detailed explanations for the following parameter choices are discussed in detail in the following Section 2.5. The individual glacier daily terminus position time series data were set as cyclical type over a 365 day seasonality period ( Table 1). The flexible level shift setting allowed for the change detected by DBEST to be either abrupt or gradual depending on whether the temporal duration of the change was one time-step (day) or longer. In general, choosing larger level shift thresholds and a longer duration-threshold is more cautious, because smaller settings may lead to false detections of seasonal variations [25]. The first level shift threshold was set to 2.675× standard deviations of the difference between actual Landsat observations of each glacier. The change magnitude threshold and second level-shift threshold were set to a statistically defined boundary between outlier and normal variations, corresponding to 2.675× standard deviations of the annual amplitudes of each time series for the individual glacier termini determined from the seasonal component results of the STL decomposition. Therefore, if a detected change was quick (between two consecutive observations/days) and large enough (2.675× standard deviations of the difference between Landsat observations) to shift the mean over the user-set duration (2 × 365 days) before and after the point, it was characterized as an abrupt change, otherwise it was considered a gradual change. The distance threshold was set to be a default that is derived internally to each glacier's cumulative terminus retreat dataset, discussed further in the following error analysis section. However, through empirical testing, the distance threshold was set to a smaller value (200 m) for four glaciers (a, f, h, and n) in order to more accurately capture the underlying trend of the data and estimate the break date (Section 2.5).

DBEST Parameter Rationale and Error Estimation
In Jamali et al. [25], the reliability of DBEST was assessed using 500 iterations of simulated Normalized Difference Vegetation Index (NDVI) data with combinations of different noise levels, site locations, and change point types (e.g., abrupt or gradual). This analysis revealed that DBEST's algorithm for change point detection was statistically significant and valid from the fit quality standpoint based on BIC [36] for optimal model selection. Identification of breakpoints was found to be largely dependent on the distance threshold parameter and noise level (inter-observation variability) of the data, or in our case, daily variability in terminus position. In general, the greater the noise, the lower the precision. However, accuracies improve and root mean square errors (RMSEs) decrease for time (days) as the distance threshold increases, but RMSEs increase for magnitude (glacier terminus change units, m) as the distance threshold increases. In order to minimize RMSEs in both time and magnitude, the default calculation of the distance threshold is based on a statistical boundary, 3×RMSE between outlier and normal residuals of an obtained fit. For example, a distance threshold of zero would equate to all data points being deemed change points and, subsequently a small change magnitude (m). A distance threshold equal to the maximum residual, or maximum absolute difference (m) between data values, would likely lead to missed detection of change points. The default 3×RMSE is calculated as significant at the α~0.01 level, leading to 99% of significant breakpoints captured and the potential for 1% of breakpoints to be either false or missed detections. This means that data points with outlier residuals will most likely be detected as breakpoints in trend by fulfilling the distance threshold criterion. Hence, it is unlikely that noise, or interpolated values, lead to false detection of these such points because they normally correspond to small residuals. In an updated but not yet published version of DBEST, error in timing and magnitude of identified change points are calculated based on duration and fit coefficient (or slope) and its error estimated in the piecewise linear modelling. Here, error estimates range between~1 and 10 days and in magnitude, between~0.4 and 7 m (Table A2). As above, choosing larger level shift thresholds is more cautious, and the first level shift should be less than or equal to the second level shift. Here, we tested the first level shift set as equal to the second level shift (at 2.675× standard deviations of the annual amplitudes). This ultimately led to the same glaciers having change points identified as presented in our results below, but with fewer break dates identified over a similar change duration. For glaciers where visual inspection suggests the potential of abrupt changes (such as a, d, e, f, m, and n), we tested how low the first level shift threshold setting needed to be in order to identify such changes. From this, we determined that the identification of abrupt change would require a first level shift setting between 2 and 15 m, depending on the glacier, which we suggest is not a large change in terminus position, even within a single day, to be considered an unlikely event. We settled on the first level shift threshold as the 2.675× standard deviations of the difference between actual Landsat observations of each glacier and a second level shift threshold of 2.675× standard deviations of the annual amplitudes. This created a first level shift that was less than the second level shift for all glaciers except Kangilliup Sermia (b), likely owing to its strong seasonal variability. In the case of Kangilliup Sermia (b), these level shift thresholds were within 8 m based on the above rules, and we chose to set them both equal at the lesser of the two values. Final settings for first level shift range from 60 to 492 m, and second level shift range from 75 to 1022 m. Ultimately, in this study, we find that all change points were of the gradual type, rather than abrupt.

Annual Amplitude in Terminus Position
We estimated the annual amplitude in terminus position as the difference between maximum and minimum points of the seasonal component for each year generated by the STL decomposition in DBEST. In a given year, we typically observed a trend of retreat that persisted from late spring to early fall, followed by a return to a trend of advance which was sustained into the first observations of the following calendar year. Therefore, we do not sample the exact farthest retreated/advanced position because of the limitations of the interpolation of our observations, the Landsat repeat cycle, and the highly variable glacier terminus on a daily time scale [1,13,15,32,37]. While not exact, our observations Remote Sens. 2020, 12, 3651 9 of 25 capture the season during which the majority of annual terminus change occurs and offer insights into the broader patterns of the seasonal variability at individual glacier termini. Given the lack of Landsat image availability for this region, particularly in the 1990s, we calculated the annual amplitude in terminus position for only the interpolated 2000-2015 period input into DBEST, owing to more frequent Landsat availability and a greater number of observations of all 17 glaciers.

Cumulative Terminus Change
Using the box method, we calculated the cumulative retreat over the 1985-2015 period for all 17 glaciers to average 1213 ± 1172 m ( Table 2). For glaciers that drain into Uummannaq Bay, the average was 1264 ± 1437 m, ranging from an advance of 131 m at Sermeq Kujalleq (j) to a retreat just under 4000 m at Sermeq Silarleq (f) (Figure 3a). In Disko Bay, the average cumulative retreat was 1139 ± 885 m, with the greatest retreat of nearly 2400 m at Kangilerngata Sermia (m) and the least of 50 m at Sermeq Avannarleq (k) (Figure 3b). Overall, we found that 15 of the 17 glaciers retreated over the entire 1985-2015 period and the majority (~80%) of this retreat occurred over the 2000-2015 period (Figure 3a

Annual Amplitude in Terminus Position
We observed a signal of advance and retreat at all termini during each year of the study period (Figure 4b). Over the 2000-2015 period for which we had more frequent measurements for all 16 years, we calculated an average amplitude in terminus position of 184 ± 19 m in Uummannaq Bay,

Annual Amplitude in Terminus Position
We observed a signal of advance and retreat at all termini during each year of the study period (Figure 4b). Over the 2000-2015 period for which we had more frequent measurements for all 16

Change Points Analysis
A change point (or breakpoint), in DBEST, was defined as a point in time at which the direction of a glacier's terminus retreat (advance) trend started to change significantly. We chose to use only the 2000-2015 year period of glacier terminus observations because of the higher frequency that allowed for Akima spline interpolation into daily spaced time intervals as input into the DBEST software. Our cumulative terminus position results also support previous studies that have shown this 16 year period as that over which the majority of retreat occurred (e.g., 2; 4). Change start dates that were calculated at the first observation of a glacier's time series (e.g., 2000-2002 at Perlerfiup Sermia (e), Figure 6a) were not included, as this was likely a consequence of the daily interpolation performed.
We detected the gradual type of change points at 10 glaciers-five in Uummannaq Bay (a, d, e, f, and h) (Figure 6a) and five within Disko Bay (m, n, o, p, and q) (Figure 6b). Kangerluarsuup Sermia (d) underwent three change durations, while the remainder of the glaciers (a, e, f, h, m, n, o, p, and q) underwent two. Most change magnitudes (m) were negative, or differences in terminus position between identified change points indicated increases in the rate of retreat, ranging between a change value of -98 m at Sermilik (h) and -2780 m at Sermeq Silarleq (f). We also observed a positive change value of 85 m at Sermilik (h), indicating changes in the terminus advance trend. However, we note this change was within 5 m of the change magnitude threshold (e.g., 2.675× standard deviations of the annual amplitude).

Change Points Analysis
A change point (or breakpoint), in DBEST, was defined as a point in time at which the direction of a glacier's terminus retreat (advance) trend started to change significantly. We chose to use only the 2000-2015 year period of glacier terminus observations because of the higher frequency that allowed for Akima spline interpolation into daily spaced time intervals as input into the DBEST software. Our cumulative terminus position results also support previous studies that have shown this 16 year period as that over which the majority of retreat occurred (e.g., 2; 4). Change start dates that were calculated at the first observation of a glacier's time series (e.g., 2000-2002 at Perlerfiup Sermia (e), Figure 6a) were not included, as this was likely a consequence of the daily interpolation performed.
We detected the gradual type of change points at 10 glaciers-five in Uummannaq Bay (a, d, e, f, and h) (Figure 6a) and five within Disko Bay (m, n, o, p, and q) (Figure 6b). Kangerluarsuup Sermia (d) underwent three change durations, while the remainder of the glaciers (a, e, f, h, m, n, o, p, and q) underwent two. Most change magnitudes (m) were negative, or differences in terminus position between identified change points indicated increases in the rate of retreat, ranging between a change value of -98 m at Sermilik (h) and −2780 m at Sermeq Silarleq (f). We also observed a positive change value of 85 m at Sermilik (h), indicating changes in the terminus advance trend. However, we note this change was within 5 m of the change magnitude threshold (e.g., 2.675× standard deviations of the annual amplitude).  We observed that the majority of glaciers (a, d, e, f, m, n, o, p, and q) which had undergone change points had also cumulatively retreated >300 m. As such, two groupings of general glacier behavior are revealed based on a 300 m threshold in cumulative retreat over the 2000-2015 period, whereby glaciers with cumulative retreat larger than this threshold showed a significant correlation (R = 0.95; p < 0.05) with average annual amplitude in terminus position (Figure 7). For glacier termini with cumulative retreat <300 m, we do not find a significant relationship with average annual amplitude in terminus position. We therefore defined these groupings as those glaciers which underwent "active retreat" (>300 m) and those that were stable, or had no discernable trend, (<300 m) between 2000 and 2015. Sermilik (h) that underwent change in the terminus advance signal (positive change in trend) remained stable with a cumulative retreated <300 m.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 26 We observed that the majority of glaciers (a, d, e, f, m, n, o, p, and q) which had undergone change points had also cumulatively retreated >300 m. As such, two groupings of general glacier behavior are revealed based on a 300 m threshold in cumulative retreat over the 2000-2015 period, whereby glaciers with cumulative retreat larger than this threshold showed a significant correlation (R = 0.95; p < 0.05) with average annual amplitude in terminus position (Figure 7). For glacier termini with cumulative retreat <300 m, we do not find a significant relationship with average annual amplitude in terminus position. We therefore defined these groupings as those glaciers which underwent "active retreat" (>300 m) and those that were stable, or had no discernable trend, (<300 m) between 2000 and 2015. Sermilik (h) that underwent change in the terminus advance signal (positive change in trend) remained stable with a cumulative retreated <300 m.  (Table A1). Active retreat (hollow circles) are those glaciers that have retreated >300 m (gray dashed line), while stable glaciers (black circles) are those which have advanced or retreated <300 m.
To garner a better understanding for the timing and magnitude of the retreat occurring at the actively retreating glaciers determined by greater than 300 m of cumulative retreat in Figure 7, we further investigated the annual amplitude behavior of the nine glaciers (a, d, e, f, m, n, o, p, and q). As defined by the DBEST change point outputs, we compared the total cumulative retreat to the average annual amplitude in terminus position specifically throughout the defined duration of change (e.g., between the red points connected by red lines, Figure 6a;b) with those during the periods of more typical glacier terminus behavior (e.g., outside of red points, Figure 6a;b). At all nine glaciers, we found that the majority of the cumulative retreat occurred during the duration of change defined by the change points (Figure 8), indicating a good performance of the DBEST algorithm. We found that glacier Umiammakku Sermiat (a) had an average annual amplitude in terminus positions 3× greater throughout the duration of change as opposed to outside of this period, while the average annual amplitude of the other glaciers (d, e, f, m, n, o, p, and q) remained nearly consistent, ± 50-100 m, whether within the DBEST defined change duration or not.  To garner a better understanding for the timing and magnitude of the retreat occurring at the actively retreating glaciers determined by greater than 300 m of cumulative retreat in Figure 7, we further investigated the annual amplitude behavior of the nine glaciers (a, d, e, f, m, n, o, p, and q). As defined by the DBEST change point outputs, we compared the total cumulative retreat to the average annual amplitude in terminus position specifically throughout the defined duration of change (e.g., between the red points connected by red lines, Figure 6a;b) with those during the periods of more typical glacier terminus behavior (e.g., outside of red points, Figure 6a;b). At all nine glaciers, we found that the majority of the cumulative retreat occurred during the duration of change defined by the change points (Figure 8), indicating a good performance of the DBEST algorithm. We found that glacier Umiammakku Sermiat (a) had an average annual amplitude in terminus positions 3× greater throughout the duration of change as opposed to outside of this period, while the average annual amplitude of the other glaciers (d, e, f, m, n, o, p, and q) remained nearly consistent, ± 50-100 m, whether within the DBEST defined change duration or not. Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 26 Figure 8. Letters refer to individual glaciers (Table A1). For the actively retreating (>300 m) glaciers defined in Figure 7, the cumulative terminus change (m) and average annual amplitude in terminus position (m) within change durations (hollow circles) and outside change durations (black circles), as defined by the start date and end dates of change durations in DBEST (e.g., red squares in Figure 6a;b).

Cumulative Terminus Retreat and Annual Amplitude in Terminus Position
In this study, we present a dataset of 17 marine-terminating glacier frontal positions that drain into Disko and Uummannaq Bays, West Greenland, and provided an analysis of the relationship between their seasonal variability and frontal position change over the last three decades. While we observed that the majority (15 of 17) of the glaciers retreated over the 30 year record and the majority (~80%) of this retreat occurred since 2000 (Table 2), the duration and magnitude of overall retreat was highly variable among the individual glaciers (Figure 3a;b) despite experiencing similar regional increases in atmospheric and sea surface temperatures.
We found general agreement with the long-term terminus retreat (advance) behavior of our study glaciers compared with previous literature, noting some differences likely attributable to retreat calculation methodologies [33], study period [32], and temporal/spatial resolution of the remote sensing platform utilized [1]. In general, there has been an increase in the magnitude of glacier retreat in this area of Greenland from the period 1992-2000 to 2000-2006, with retreat onset correlated with atmospheric warming since 2000 [32]. Jensen et al. [16] investigated area changes at 42 glacier termini, including two glaciers in Disko Bay and 10 in Uummannaq Bay, and found that the majority (35 of 42) of their study glaciers retreated between 1999 and 2013, while the other seven glaciers, including Sermeq Kujalleq (j) studied here, remained stable. Similarly, we find agreement with previous studies who calculated that the terminus position of Sermeq Kujalleq (j) had not changed since 1948 [24], and gained mass since 1985 [7]. Consistent with the hypothesis of perturbations at the terminus propagating thinning up-glacier because of changes in the stress regime [1,7,38], the three glaciers we measured as having undergone the largest terminus retreat, Umiammakku Sermiat (a) and Sermeq Silarleq (f) in Uummannaq Bay, and Kangilerngata Sermia (m) in Disko Bay, also accounted for ~84% of the regional mass loss between 1985 and 2015 [7]. Changes in glacier flow  (Table A1). For the actively retreating (>300 m) glaciers defined in Figure 7, the cumulative terminus change (m) and average annual amplitude in terminus position (m) within change durations (hollow circles) and outside change durations (black circles), as defined by the start date and end dates of change durations in DBEST (e.g., red squares in Figure 6a;b).

Cumulative Terminus Retreat and Annual Amplitude in Terminus Position
In this study, we present a dataset of 17 marine-terminating glacier frontal positions that drain into Disko and Uummannaq Bays, West Greenland, and provided an analysis of the relationship between their seasonal variability and frontal position change over the last three decades. While we observed that the majority (15 of 17) of the glaciers retreated over the 30 year record and the majority (~80%) of this retreat occurred since 2000 (Table 2), the duration and magnitude of overall retreat was highly variable among the individual glaciers (Figure 3a;b) despite experiencing similar regional increases in atmospheric and sea surface temperatures.
We found general agreement with the long-term terminus retreat (advance) behavior of our study glaciers compared with previous literature, noting some differences likely attributable to retreat calculation methodologies [33], study period [32], and temporal/spatial resolution of the remote sensing platform utilized [1]. In general, there has been an increase in the magnitude of glacier retreat in this area of Greenland from the period 1992-2000 to 2000-2006, with retreat onset correlated with atmospheric warming since 2000 [32]. Jensen et al. [16] investigated area changes at 42 glacier termini, including two glaciers in Disko Bay and 10 in Uummannaq Bay, and found that the majority (35 of 42) of their study glaciers retreated between 1999 and 2013, while the other seven glaciers, including Sermeq Kujalleq (j) studied here, remained stable. Similarly, we find agreement with previous studies who calculated that the terminus position of Sermeq Kujalleq (j) had not changed since 1948 [24], and gained mass since 1985 [7]. Consistent with the hypothesis of perturbations at the terminus propagating thinning up-glacier because of changes in the stress regime [1,7,38], the three glaciers we measured as having undergone the largest terminus retreat, Umiammakku Sermiat (a) and Sermeq Silarleq (f) in Uummannaq Bay, and Kangilerngata Sermia (m) in Disko Bay, also accounted for~84% of the regional mass loss between 1985 and 2015 [7]. Changes in glacier flow velocity across the study region also demonstrated spatial and temporal variability between 1999 and 2012 [39].
Within individual years, we observed regular summer retreat and winter advance behavior at all of the termini, with the minimum and maximum positions generally occurring in August/September and April/May, respectively (Figure 6a;b). The spring onset of the retreat season we observed is in agreement with previous work that observed coincident increases in spring velocity [1,2], potentially related to runoff or ice mélange breakup [10]. Using the DBEST seasonal component, our average annual amplitude in terminus position of 156 ± 17 m since 2000 for the Disko and Uummannaq Bay region of West Greenland is smaller than previous observations of seasonal change in Northwest Greenland [9,40].  [1,2,14,41]. Previous studies indicated that the velocity of Umiammakku Sermiat (a) had yet to stabilize in 2009 [1,2]. Our extended record through 2015 indicated the end of an active retreat period in December 2009, when we observed the terminus of Umiammakku Sermiat (a) return to seasonal and cumulative behavior similar to pre-change point conditions (e.g., pre-2003). This suggests a likely concurrent decrease in surface speed from the end date of change in December 2009 to the end of our 2015 record. In general, the behavior of Umiammakku Sermiat (a) has been well documented and has been described as providing a good example of stepwise changes in terminus positions, generally associated with movement between underlying fjord and bedrock topographic pinning points [9,40], where the rate of movement between these pinning points is controlled by the individual fjord forcings (e.g., ice mélange, ocean temperatures at depth). More specifically, Catania et al. [4] correlate the magnitude and duration of retreat at Umiammakku Sermiat (a) with overdeepening, or the distance behind the terminus located on a retrograde bed slope. Although Sermeq Silarleq (f) remained stable between August 2003 and December 2009, we then observed the start of a second change duration and increase in annual amplitude that persisted until April 2012. The start date of this second change duration followed the largest observed August SST anomaly in Uummannaq Bay in 2009 [15]. While the trigger for the onset of active retreat seems to correspond with regional oceanic anomalies, the asynchronous response between Umiammakku Sermiat (a) and Sermeq Silarleq (f) at the end of the record further emphasizes the geometric controls on terminus behavior which remain unique to each glacier's fjord, including bed slope [4], runoff at the grounding line and ice mélange [10], ice front shape, and the relationship between calving speed and subaqueous melt rate [11].
Kangerluarssup Sermia (d) behaved differently from the other glaciers with change points in Uummannaq Bay regarding that it retreated more gradually over the 2000-2015 period, which has previously been noted as another dominant type of terminus retreat behavior [9,40]. The start date of the change duration at Kangerluarssup Sermia (d) also occurred in April 2002, but showed little variability in annual amplitude during the period of active retreat. However, this start date of change detected in 2002 reinforces the potential connection to the regional winter oceanographic influences seen at the start of other active terminus retreat periods of glaciers in Uummannaq Bay. The steepest trend in cumulative retreat at Kangerluarssup Sermia (d) occurs during the second change duration, with a start date October 2009. This trend persists until immediately entering the third change duration, starting in in April 2011 through the end of the study record. A lack of bathymetric pinning points which can promote calving [42] and narrow overdeepenings behind the terminus [4] are typically observed at glaciers with more gradual trends in terminus behavior, what we observed at Kangerluarssup Sermia (d). As such, it is predicted that Kangerluarssup Sermia (d) may continue to gradually retreat until reaching an underlying topographic pinning point such as what has been found at glaciers in West [4] and Northwest Greenland [9].
In contrast, the start date of active terminus change at Perlerfiup Sermia (e) did not begin until August 2008. Furthermore, it immediately entered a second, accelerated retreat change duration beginning in January 2013. Unlike Kangerluarssup Sermia (d) which showed little variability in annual amplitude during the period of active retreat, this second change at Perlerfiup Sermia (e) occurred in combination with an increase in annual amplitude in 2013 that continued through the end of the 2015 record. Here, accelerated terminus retreat at Perlerfiup Sermia (e) followed multiple years of anomalously low sea ice conditions in the region [28]. Perlerfiup Sermia (e) and Kangerluarssup Sermia (d) have deep fjords (>800 m) with the terminus of Perlerfiup Sermia (e) grounded at the end of the deep fjord in shallow, cold, fresh water [11]. Rignot et al. [11] observed these deeper fjords to be generally ice free, suggesting a less constrictive fjord geometry easier for icebergs to exit. Ice-free fjords also point to the glacier calving rate, the size of icebergs relative to that of the fjord, and the strength of the fjord circulation to move them out [43]. Such a lack of icebergs in combination with low sea ice formation years observed throughout the 2007-2012 change durations identified here, may have led to a weaker mélange formation in the winter and the onset of active terminus retreat at Perlerfiup Sermia (e) and Kangerluarssup Sermia (d). While both Perlerfiup Sermia (e) and Kangerluarssup Sermia (d) remained in a period of active retreat at the end of our 2015 study record, proximal prograde bed slopes behind their 2016 terminus position predicts they are likely to stabilize soon [4].

Active Glaciers-Disko Bay
In Disko Bay, there were similar start dates (2001-2003) of the first change duration at the actively retreating (>300 m cumulative retreat, 2000-2015) glaciers (m, n, o, p, and q) compared to those in Uummannaq Bay. The termini of glaciers (m, n, o, p, and q) also underwent a second change duration with a start date between 2007 and 2012 (Figure 6b), similar in timing to that observed at Perlerfiup Sermia (e) in Uummannaq Bay (Figure 6a).
Similar to the actively retreating glaciers in Uummannaq Bay, at Eqip Sermia (n), Sermeq Avannarleq (o) and Saqqarliup Sermia (q) in Disko Bay, the start date of the first change duration occurred in spring of 2002 and 2003. However, the start date of the first change durations at Kangilerngata Sermia (m) and Alanngorliup Sermia (p) differs from the rest of the actively changing glaciers in Disko and Uummannaq Bays regarding that they occurred in the prior year (2001) and began at the end of summer (September). This suggests that the topographic conditions of the bed/fjord at Kangilerngata Sermia (m) and Alanngorliup Sermia (p) may have led to terminus positions more vulnerable to environmental triggers such as the increased influence of the warm Irminger current in Disko Bay since the 1990s [44,45], and ultimately allowing for the onset of active retreat. Previous work has suggested Kangilerngata Sermia (m) and Eqip Sermia (n), in particular, have been dominated by submarine melting as a primary driver of mass loss [46].
The start date of the second change durations was between September 2007 and October 2009 at glaciers (m, n, o, and q), with the active retreat periods at Eqip Sermia (n) and Saqqarliup Sermia (q) persisting through the end of the record in 2015. Eqip Sermia (n) has a central bump on which it remained pinned, but also multiple overdeepened channels within its fjord that has led to unequal retreat along the terminus [45]. Once pushed past this bump, Eqip Sermia (n) is likely to continue to retreat through a deep upstream trough [4]. Given Eqip Sermia (n) was still undergoing a period of active retreat at the end of our 2015 study period and its known upstream bed topography, it may be of particular interest in more recent years as oceanic triggers have continued to create conditions for maintaining an active retreat trend. At Alanngorliup Sermia (p), the second duration of change started in 2012, similar in timing and behavior to the observations of Perlerfiup Sermia (e) in Uummannaq Bay, both of which did not retreat between 1985 and 1999, but almost entirely within the 2000-2015 period. This suggests similar topographic vulnerabilities of Alanngorliup Sermia (p) and Perlerfiup Sermia (e) to the lows in regional sea ice concentrations that occurred between 2007 and 2012. In general, we suggest all of the glaciers that underwent active terminus change in Disko Bay follow the common stepwise behavior [9,40] associated with the movement of termini between periods of decreased retreat while being pinned on underlying topographic highs, and increased retreat while progressing towards another pinning point.

Stable Glaciers-Disko and Uummannaq Bays
The termini of seven glaciers (b, c, g, h, j, k, and l) followed no discernable trend over the 2000-2015 period, retreating less than 300 m over the 2000-2015 period (Figure 7). These stable glaciers (b, c, g, h, j, k, and l) in both bays are found in a wide range of fjord depths (300-900 m); however, several of the termini are grounded in warmer Atlantic waters (>2.5 • C) [11]. Two glaciers, Kangerlussuup Sermia (c) and Sermeq Kujalleq (j) cumulatively advanced over the entire 1985-2015 record ( Table 2). The stability of Kangerlussuup Sermia (c) and Sermeq Kujalleq (j) was found to be correlated with the presence of an extensive prograde bed slope behind their terminus [4]. Since 2000, three additional glaciers Kangilleq (g), Sermilik (h), and Sermeq Avannarleq (k) advanced, although this recent advance is approximately one-third of their average annual amplitude in terminus position. The behavior of these five advancing glaciers (c, g, h, j, and k) is in agreement with previous work which found little terminus change since 1964 [11]. Furthermore, Kangerlussuup Sermia (c), Sermeq Kujalleq (j), and Sermeq Avannarleq (k) thickened since 1985 [7], and Sermeq Kujalleq (j), Sermeq Avannarleq (k), Sermeq Kujalleq (l) decreased in flow speed over the 1999-2012 period [39].
In general, we found a wide range in average annual amplitude at all of the stable glaciers. Kangilliup Sermia (b) may be of particular interest because it shows the overall largest average annual amplitude of 480 m compared to all study glaciers, both stable and active. The annual amplitude of Kangilliup Sermia (b) noted here is in agreement with the finding of a strong seasonal signal when it was compared to the other large (>3 km width) glaciers in Uummannaq Bay [1]. However, its seasonal variability was calculated to be~1 km using a higher temporal resolution dataset [15]. This seasonal cycle observed at the terminus of Kangilliup Sermia (b) was found to be most impacted by ice mélange when compared with influences, such as runoff at the grounding line and ocean temperatures [10]. Kangilliup Sermia (b) cumulatively retreated 1119 m over the 1985-2015 period with the majority (~75%) of this retreat occurring prior to 2000, yet Kangilliup Sermia (b) is considered stable here because it retreated <300 m (289 m) over the 2000-2015 period, likely owing to the more recent position of the terminus being grounded upstream of topographic highs [4]. Fjord bathymetry and ocean temperatures at depth indicated that the rate of ice melt relative to the rate of upstream ice advection might not be large enough to initiate terminus retreat at Kangilliup Sermia (b) [11].
For both the stable and actively retreating glaciers identified in this study, individual glacier geometry seems to both create and hinder the opportunity for environmental triggers to onset active terminus retreat. The front positions digitized in this study [47] could be used for further analysis to quantitatively compare terminus changes with, for example, recent bathymetric and bed data [48] or velocity measurements [39]. More research needs to be done in order to determine what is considered a significant amount of change in glacier terminus position and over what time period it needs to occur in order to be considered an abrupt event. Future accuracy assessment of the application of DBEST to detection of change points in glacier termini could also include the use of simulated data and testing for DBEST's ability to detect previously identified or known tipping points in glacier retreat (advance) behavior.

Conclusions
We found a significant positive relationship between average annual amplitude in terminus position and total cumulative retreat among the 10 glaciers with a cumulative retreat >300 m since 2000 in Disko and Uummannaq Bays, West Greenland. We found that the majority (nine) of the 10 glaciers that retreated >300 m also experienced periods of "active retreat" as defined by change points detected within the DBEST program. In general, cumulative retreat was greater during the periods of active change compared to inactive periods at these nine glaciers, while annual amplitude did not vary as much whether within or outside an identified duration of change. This suggests a good performance of the DBEST algorithm in identifying significant transitions in trends and points to a connection between seasonal terminus variability and long-term retreat. We found qualitatively coincident change point start times with both peaks in sea surface temperatures and lows in regional sea ice concentrations. In agreement with previous studies, marine conditions may be of particular importance as triggers of future outlet glacier terminus retreat (advance) trends, with subsequent duration and magnitude of individual glacier response to these triggers likely controlled by unique fjord/bedrock topography. This suggests only if fjord/bedrock conditions are ideal will the glacier's terminus position be vulnerable to an active change onset from anomalous conditions in typical climate triggers for retreat (e.g., SST and sea ice cover). This supports our observation of a strong spatial variability and unique timing/magnitude of terminus change, despite subjection to similar regional environmental changes in atmospheric temperature, SST, and sea ice cover across the study region. Glacier termini found in this study to have been remaining within a period of active retreat at the end of the 2015 record are particularly important because these environmental factors have only continued to create potentially more vulnerable conditions for active terminus retreat. The seasonal, decadal, and change point analysis of outlet glacier termini we present here will be particularly useful in future projections of non-linearity in glacier terminus retreat as climate warming across Greenland continues.