Assessing the effects of fjord geometry on Greenland tidewater glacier stability

Abstract

Tidewater glaciers frequently advance and retreat in ways uncoupled from climate forcing. This complicates the task of forecasting the evolution of individual glaciers and the overall Greenland ice sheet, much of which is drained by tidewater glaciers. Past observational research has identified a set of processes collectively known as the tidewater glacier cycle (TGC) to describe tidewater glacier evolution in four stages: the advancing stage, the extended stage, the retreating stage and the retreated stage. Once glacier retreat is initiated, the TGC is thought to depend largely on the glacier's calving rate, which is controlled by fjord geometry. However, there has been little modeling or systematic observational work on the topic. Measuring calving rates directly is challenging and thus we developed an averaged von Mises stress state at the glacier terminus as a calving rate proxy that can be estimated from surface velocities, ice thickness, a terminus position and subglacial topography. We then analyzed 44 tidewater glaciers in Greenland and assessed the current state in the TGC for them. Of the 44 glaciers, we find that fjord geometry is causing instability in ten cases, vs stability in seven, with 11 in rapid retreat and 16 have been historically stable.

1. Introduction

It has long been known that tidewater glaciers advance and retreat out of sync with land-terminating glaciers and external ocean and climate forcing (Post, 1975; Meier and Post, 1987; Pfeffer, 2003). This observation has led to the formulation of the tidewater glacier cycle (TGC) (Pfeffer, 2007; Pollard and DeConto, 2009; Post and others, 2011; Brinkerhoff and others, 2017), a combination of processes that proceed in four archetypal phases, as described in Brinkerhoff and others (2017). In the advancing stage, development and advection of a shoal at the front reduces calving and submarine melting, causing glacier thickening and advance. Eventually the glacier enters an extended phase, in which accumulation and ablation are in balance and further advance is halted. A glacier enters the retreating phase when it can no longer maintain sufficient thickness to remain grounded on the shoal and the glacier retreats into progressively deeper water; at which point dramatic unstable rapid retreat takes place. Retreat ends when the terminus approaches a position that reduces the calving in the absence of sedimentation; possibly at a pinning point (temporary narrowing of the fjord), or else the terminus effectively re-grounds on bedrock. Multi-decadal modeling studies with sedimentation are able to reproduce this cycle (Nick and others, 2007; Amundson, 2016; Brinkerhoff and others, 2017) even in the absence of variations in climate (Brinkerhoff and others, 2017).

Typical timescales and rates of advance/retreat for tidewater glaciers is an area of active research. Catania and others (2018) report rates of retreat for Greenland tidewater glaciers of up to 500 m a⁻¹. Brinkerhoff and others (2017) show a simulated retreat phase lasting ~100 years and an advance phase ~1000 years, with the terminus advancing or retreating ~5 km in that time (50 m a⁻¹ retreat and 5 m a⁻¹ advance). Carlson and others (2017) report the Columbia Glacier in Alaska has retreated ~20 km in 30 years (667 m a⁻¹). Pearce and others (2022) report an advance rate of Kangiata Nunaata Sermia of 115 m a⁻¹ during the Little Ice Age (12th and 13th centuries CE). In general, the retreat phase of the TGC is thought to happen on a decadal or centennial timescale, and the advance phase is about an order of magnitude slower.

In this observational study, we focus on detecting glaciers beginning, sustaining or finishing the retreat phase using the ITS_LIVE surface velocity dataset from 1985 through 2018, or 33 years (Gardner and others, 2019). We use annual terminus positions, thereby sidestepping the complex issue of seasonal variability. Our time series are not long enough to investigate the advance phase. Sedimentation is an important part of the TGC at the century timescale (Brinkerhoff and others, 2017), but one we can safely ignore in this study covering rapid retreat of numerous glaciers over a few decades.

Greenland's outlet glaciers are currently at diverse stages in the TGC. For example, after advancing 800 m from 1973 to 2000, Sermeq Silarleq in central west Greenland retreated 5 km from 2000 through 2019; but just 47 km to the south, Store Glacier has remained remarkably stable during the same time frame (Cheng and others, 2021a). Interestingly, the TGC suggests that tidewater glacier termini can only remain in stable equilibrium, neither advancing or retreating, at places where further advance would cause a negative feedback: at fjord mouths, at pinning points (temporary narrowing of the fjord) and at other places involving change in the fjord width (Mercer, 1961).

Many glaciers around Greenland have been retreating since 2000 (Murray and others, 2015), suggesting that conditions required for TGC advance are currently rare. While the TGC has been used as a post-hoc explanation for advancing (McNeil and others, 2021) and retreating glaciers, it has not been used to systematically investigate how a glacier's current behavior reflects the TGC phases. Quantifying where a glacier currently falls in the TGC could help to understand how it might respond to future environmental change.

1.1. Basics of advance and retreat

Wood and others (2021) present a model for Greenland in which the position L of a tidewater glacier's grounded terminus is a result of four competing processes causing advance or retreat: advection of ice downstream (q _f) leads to terminus advance, whereas frontal melt (q _m), calving (q _c) and thinning-induced retreat (q _s) (Felikson and others, 2017) lead to terminus retreat. Adopting the convention that positive sign means advancing terminus for all advance/retreat rates (unit: m s⁻¹), mass balance at the ice front requires:

(1)$$\Delta L = L - L_0 = \int_{t_0}^{t} ( q_{\rm f} + q_{\rm m} + q_{\rm c} + q_{\rm s}) \, {\rm d}t,\; $$

where L is the current terminus position, L ₀ is the terminus position at a reference time t ₀ and t is the current time. The values of q _s computed by Wood and others (2021) and shared in that paper's supplement are at least an order of magnitude smaller than q _c for the glaciers in this study, and we can therefore ignore q _s. Wood and others (2021) observe or model all the terms of Eqn (1) except for q _c.

Wood and others (2021) use a parameterization for frontal melt, similar to Slater and others (2019), in which subglacial discharge and thermal forcing are both derived from an ensemble of MITgcm runs (Xu and others, 2013; Rignot and others, 2016): because ocean gridcells are too large to resolve fjords, model-based ocean temperatures at the mouth of each fjord are translated into temperatures inside the fjord at the calving front. (See section titled Thermal forcing in Wood and others (2021), page 8 of 10.) Thinning-induced retreat q _s is calculated using a simple geometrically derived relationship for grounding line migration rate as a function of surface elevation change (Thomas and Bentley, 1978).

Disregarding frontal melt for a moment, stability may be investigated systematically by evaluating the relationship between terminus position and calving rate q _c for each glacier. If a currently stable glacier is about to enter the retreat phase of the TGC, then q _c would be expected to increase as the glacier begins to retreat, potentially leading to runaway retreat; whereas if the glacier is stabilizing, q _c would be expected to decrease as the glacier retreats, causing retreat to slow.

If one could measure or model all four components of Eqn (1), glaciers about to enter or exit the rapid retreat phase of the TGC could be identified by investigating the correlation between observed changes in q _c vs observed changes in terminus position. However, such observations are challenging because of limitations in observational capabilities and uncertainties in process models. While calving rates can be measured with localized systems (Walter and others, 2020; Taylor and others, 2022), earth observing satellites do not provide high enough resolution in time to apply these techniques over a wide region. For this reason, a proxy $\bar {\sigma _{\scriptscriptstyle T}}$, representing relative levels of expected q _c and computable from remote-sensing data, is used to evaluate glacier stability instead of the calving rate q _c. This proxy relies on the von Mises calving law (von Mises, 1913; Morlighem and others, 2016; Choi and others, 2018) as a simple but reasonable model for tidewater glaciers (Section 6.1).

Even with a reliable proxy for q _c, efforts to find a relationship between q _c and observed retreat ΔL will fail: due to warming in the ocean around Greenland, frontal melt has recently become the dominant process driving retreat (Slater and others, 2019; Wood and others, 2021), with calving playing a secondary role. Most glaciers are retreating, and we cannot immediately conclude that observed retreat is due to tidewater glacier instability. Some glaciers continue to retreat even as they move into shallower water, for example Lille Glacier (Fig. 13). To demonstrate a correlation between glacier retreat and calving rate, it is first necessary to estimate and remove the amount of retreat caused by changes in frontal melt rate q _m. We use the empirical model of Slater and others (2019) for that task. Note that calving rates are affected by frontal melt; and because the model of Slater and others (2019) is empirical, q _m will include changes in frontal melt and calving due to changes in ocean conditions. This leaves the proxy for q _c representing only changes in calving due to fjord geometry, the main driver of the TGC.

2. Methodology

We develop an averaged proxy $\bar {\sigma _{\scriptscriptstyle T}}$ for the calving rate q _c, which can be computed from readily available observations of surface velocities, ice thickness, a terminus line and subglacial topography. The proxy is derived from the von Mises calving law (von Mises, 1913; Morlighem and others, 2016; Choi and others, 2018), which this study shows in Section 7 can be tuned to be consistent with observations.

Values of $\bar {\sigma _{\scriptscriptstyle T}}$ and L in the recent history of each glacier are regressed against each other, allowing diagnosis of the glacier's stability and current stage in the TGC. In this study we use annual averages of L, thereby sidestepping issues of seasonal melting driven by submarine discharge; and we also remove the dominant effects of ocean warming. If $\bar {\sigma _{\scriptscriptstyle T}}$ is found to decrease as the glacier retreats, then the fjord geometry at that point in space is destabilizing, providing evidence that the glacier may be entering the retreat phase of the TGC; whereas if $\bar {\sigma _{\scriptscriptstyle T}}$ increases the opposite is true and the fjord geometry is now stabilizing the glacier, suggesting the glacier may be finishing the retreat phase of the TGC and moving up onto land. If $\bar {\sigma _{\scriptscriptstyle T}}$ remains the same as the glacier retreats, the glacier may be retreating through a section of the fjord with nearly constant cross-sectional width and depth, typical for glaciers in the middle of rapid retreat. In summary:

(2)$$\matrix{ \displaystyle{{\rm d}\bar{\sigma_{\scriptscriptstyle T}}\over {\rm d}L} < 0 & \hbox{destabilizing fjord geometry},\; \cr \displaystyle{{\rm d}\bar{\sigma_{\scriptscriptstyle T}}\over {\rm d}L} = 0 & \hbox{in rapid retreat},\hfill\; \cr \displaystyle{{\rm d}\bar{\sigma_{\scriptscriptstyle T}}\over {\rm d}L} > 0 & \hbox{stabilizing fjord geometry}.\hfill }$$

Finally, some glaciers have remained stable in recent years and it is not possible to tell from observations whether $\bar {\sigma _{\scriptscriptstyle T}}$ would go up or down if the terminus were to retreat.

3. Paper organization

The main hypothesis of this paper is that a glacier's stability can be assessed by observing changes in calving rate vs advance and retreat of the terminus. This hypothesis is tested by evaluating both quantities based on observations and models, and then evaluating the degree of linear correlation between them. Section 5 evaluates advance and retreat of the terminus, Section 6 evaluates changes in calving rate and Section 6.7 presents the main regression between those two quantities. Evaluating these quantities involves a combination of observations, models and statistical procedures with complex interdependencies, as illustrated by the organizational chart in Figure 2 and supported by datasets in Figure 1. Table 1 lists all symbols used in this paper. Here is a systematic summary of the following sections:

• • Section 4 describes the datasets used in this study. Some readers might wish to begin reading at Section 5 while referring to Section 4 as needed for reference.

• • Section 5 develops the terminus residual $l_{\epsilon }$, which represents the amount of calving due to fjord geometry effects as opposed to ocean warming. Supporting concepts are developed in Sections 5.1 and 5.2, and the terminus residual is presented in Section 5.3.

• • Section 6 develops the calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$, beginning in Section 6.1 with the von Mises calving law (von Mises, 1913; Morlighem and others, 2016; Choi and others, 2018). It then proceeds to regress that against the terminus residual $l_{\epsilon }$, which is the main method used to evaluate glacier stability in this paper.

• • Section 7 is a side note showing empirically that the von Mises calving law is, in fact, a reasonable calving model for the glaciers in this study. It is not required for later sections, and thus the reader might wish to skip directly to Section 8.

• • Section 8 shows how to use the calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ to classify glaciers as stable so far, in retreat, destabilizing or stabilizing (Fig. 11), and applies it to the 44 glaciers in this study. Examples and discussion of each category are provided. Excerpts of the glacier analyses are included in this paper with the full results available in Supplement S1. Table 2 summarizes the per-glacier results in a single table.

• • Sections 9 and 10 synthesize and discuss the results further.

• • Appendix A documents a novel numerical technique for evaluating line integrals on gridded data. Developed in conjunction with this study, it may be applicable for other projects.

Figure 1. Datasets used in this study. Each dataset is represented by a yellow tag, used in Figure 2. See Section 4 for further details.

Figure 2. Models and methods used in this paper: blue ovals are theoretical models, gray rectangles are methods and green rounded rectangles are methods that produce an end result of this study. Arrows represent dependencies, for example Up Area values (Section 5.2) are required to produce Terminus Residuals. Section 5 presents the frontal melt model by Slater and others (2019) driven by ocean warming, and uses it to remove effects of ocean warming from terminus data, resulting in terminus residuals. Section 6 introduces the von Mises calving law and derives $\bar {\sigma _{\scriptscriptstyle T}}$, a proxy for calving rate, which it regresses against terminus residuals to provide diagnostics on tidewater glaciers. Section 7 uses the data to show why the von Mises calving law is a reasonable model.

Table 1. Symbols used in this paper, organized by section where they are introduced

Table 2. Summary of results per glacier

Glaciers are grouped by their final categorization (Destabilizing, Stabilizing, Stable or In Retreat). Columns are ID: ID of glacier as found in Rignot and Mouginot (2012); Name: name of glacier as found in Wood and others (2021); Latitude: latitude (degrees north) of glacier combined with indication of its location on the east (E) or west (W) side of Greenland; Retreat: total amount of retreat (m) over the study period (negative for retreat, positive for advance); ν: relationship between terminus residuals and $\bar {\sigma _{\scriptscriptstyle T}}$ (Section 6.7); p-value: level of statistical significance of ν (Section 6.7); Mean $\bar {\sigma _{\scriptscriptstyle T}}$: mean value of $\bar {\sigma _{\scriptscriptstyle T}}$ for this glacier's terminus across all years.

4. Datasets

We used the following grids and datasets.

4.1. Local MEaSUREs grids

The MEaSUREs Greenland Ice Velocity dataset NSIDC-481 (Joughin and others, 2010, 2020) has already been constructed to cover many Greenland glaciers, with local high-resolution grids defined in areas with glacier activity (Fig. 3). Regridding the other datasets (below) to these local grids allows for detailed study of individual glaciers while omitting most of the interior of the ice sheet. They also allow for cross-referencing with other datasets and studies that also use the same grids. Each glacier in our analysis was identified as falling on a single local grid from NSIDC-481 (a MEaSUREs grid). For glaciers located on more than one MEaSUREs grid, the most appropriate grid for that glacier was determined by hand based on the distance from the center of the grid for each glacier. Glaciers that did not fall within a MEaSUREs grid were removed from the selection.

Figure 3. Local high-resolution grids (green rectangles) defined by the MEaSUREs dataset, NSIDC-0481.

4.2. ITS_LIVE surface velocities

Annual average surface velocities (advection rate) from 1985 through 2018, necessary to compute the von Mises stress, were obtained from the ITS_LIVE dataset (Gardner and others, 2019), and regridded to the local MEaSUREs grids. An annually averaged dataset was used to avoid complexities of seasonality, surges and short-term variability; however, use of a higher temporal resolution dataset for surface velocities might produce improved statistical precision. Mouginot and others (2017) also provide surface velocity datasets derived from satellite Landsat-8, Sentinel-1 and RADARSAT-2 data, which might be useful in similar future studies, for example in Antarctica.

4.3. BedMachine v3 subglacial topography

Subglacial topography required for this computation was provided by BedMachine v3 (Morlighem and others, 2017) and regridded to the local MEaSUREs grids. BedMachine was chosen as our best understanding of the bed underneath the Greenland ice sheet, given available data and models. Although BedMachine does not supply uncertainty estimates, it is widely believed that knowledge of the bed is least certain near each glacier terminus.

4.4. Terminus lines

Terminus positions may be computed from satellite images by tracing the terminus, either manually (Wood and others, 2021) or via machine learning (Cheng and others, 2021b; Goliber and others, 2022), with newer machine-learning approaches greatly expanding the quantity of available terminus traces. Wood and others (2021) provide one of the two main theoretical models for this study (Section 1.1), and we re-use terminus traces from it to maintain compatibility.

Slater and others (2019) provide the other main theoretical model for this study, but it represents terminus positions as a 1-D scalar distance up the fjord. This implicitly assumes a specific model of a fjord as a long narrow channel with length much greater than its approximately uniform width, which is not always reasonable. There is no simple automated way to delineate center lines, and some fjords have complex geometry not well described by a simple 1-D model. Therefore, spatial analysis in this study is conducted on a full 2-D map of the fjord. For compatibility, the scalar terminus positions of Slater and others (2019) are cross-referenced against information obtained from the 2-D terminus lines of Wood and others (2021).

4.5. Modeled frontal melt

We use two datasets from Slater and others (2019): annual scalar terminus positions l _s and annual frontal melt rate q _m, modeled as $q_{\rm m} = Q^{0.4} \hbox {TF}$, where Q is subglacial discharge due to surface meltwater runoff and basal melt and $\hbox {TF}$ is the thermal forcing in the fjord. These values are based on model ocean outputs of MITgcm (Adcroft and others, 2004). The model of Slater and others (2019) may significantly underestimate submarine melting (Sutherland and others, 2019; Catania and others, 2020; Jackson and others, 2022); but to first order we do not expect that to affect results because frontal melt is empirically calibrated to observations (Section 5.1).

4.6. Selection of glacier set

As described above, this study uses the modeling framework from Slater and others (2019), data from Wood and others (2021) and Gardner and others (2019) and grids from Joughin and others (2010). Therefore glaciers need to be present in all four datasets, resulting in a set of 44 glaciers available for the study as shown in the results (Fig. 4). Although this procedure reduces the number of glaciers for analysis, it maximizes the ability to compare and cross-reference results with previous studies. Geographic representation of glaciers, classifying by regions as defined by Wood and others (2021), is: central-west Greenland (11 of 14 total tidewater glaciers), northeast (1 of 14), northwest (15 of 64), southeast (16 of 56) and southwest (1 of 12). This study has no geographic representation in the central-east (35 total) or north (12 total) regions of Greenland. We note that the datasets we used do not include many glaciers in the southwest of Greenland.

Figure 4. Location and stability assessment of the 44 Greenland tidewater glaciers in this study. Of the 44 glaciers, 16 glaciers are stable, 7 are stabilizing, 10 are destabilizing and 11 are in retreat. Subglacial topography is from BedMachine v3 (Morlighem and others, 2014) and surface speeds from ITS_LIVE (Gardner and others, 2019).

5. Frontal melt and terminus residuals

This section develops the terminus residual $l_{\epsilon }$, which represents the amount of calving due to fjord geometry effects as opposed to ocean warming. Supporting concepts are developed in Sections 5.1 and 5.2, and the terminus residual is presented in Section 5.3.

5.1. Frontal melt model

Slater and others (2019) state that warming oceans are currently the primary driver of tidewater glacier retreat in Greenland. Based on data, they provide a glacier-by-glacier relationship between the change of the scalar terminus position l _p and frontal melt rate q _m, that is, they empirically derive κ and β, based on data averaged over 5 year intervals such that:

(3)$$l_{\rm p} = \kappa q_{\rm m} + \beta.$$

This relationship is based on the process of ice front undercutting/frontal melt only, modeled because it cannot be observed directly via remote sensing. The value $q_{\rm m} = Q^{0.4} \hbox {TF}$ from Slater and others (2019) represents the ocean heat available to drive melting. As a proxy for subglacial discharge Q, Slater and others (2019) used surface meltwater runoff estimated by the regional climate model RACMO2 (Noël and others, 2018); and for $\hbox {TF}$, they used the monthly EN4 dataset from the Hadley Center consisting of observed subsurface ocean temperature and salinity profiles (Good and others, 2013).

The linear model of Slater and others (2019) incorporates frontal melt from ocean warming but ignores the calving effects due to glacier geometry. Glaciers close to each other will experience similar changes in ocean temperature, but different fjord geometry could cause them to behave differently in spite of similar ocean forcing. Therefore, the model can predict advance or retreat of glaciers as a whole within a region due to ocean warming, but cannot predict the behavior of individual glaciers, which also depends on fjord geometry (Morlighem and others, 2017).

In recent years, ocean warming has become the dominant process causing glaciers in this study to retreat (Slater and others, 2019; Wood and others, 2021). In order to study the secondary effect of fjord geometry, the effects of the dominant process must first be removed from the data. We use the model of Slater and others (2019) to estimate the amount of retreat caused by ocean warming and subtract that out of the total retreat, leaving a terminus residual (Section 5.3) in which retreat due to calving is the dominant process.

5.2. Computing glacier retreat

Spatial analysis in this study is conducted on a full 2-D map of the fjord. In place of a scalar terminus position L, the scalar up area A _T is used, defined as the entire ice-covered area upstream of the glacier terminus T for which the basal topography is below sea level. This avoids assumptions about fjords, their linear geometry or center lines.

The up area is calculated as follows (Fig. 5). Using GIS software, a rough polygon is manually drawn around the fjord by hand, and a single point is identified in the upper reaches of the fjord (the up point). The fjord is determined in raster form by identifying gridcells within the polygon with bed below sea level. The terminus line is extended to the full width of the fjord and rasterized, to produce the set of gridcells on the terminus. A raster flood fill algorithm is then used, starting from the up point, to identify all the gridcells of the fjord that are upstream of the terminus. The up area is computed by summing the areas of these upstream grid cells. Because of the way these rough polygons are manually drawn, A _T does not include far-upstream areas of some fjords. Therefore, up area may be used for relative comparison between termini, but not as an absolute measure of how much fjord ‘remains’ ice covered before a glacier becomes land terminating.

Figure 5. Aerial map of AP Bernstorff Glacier in Southeast Greenland, with terminus as of 2005. Digitized terminus datasets typically come in vector format (black line on top of red gridcells), which is rasterized (red gridcells). To help the computer determine the extent of the fjord, we drew a rough polygon around the fjord by hand (red shaded area), and identified a point (red star) that is upstream of all expected termini used in this study. Based on these inputs and bathymetry from BedMachine, the computer was able to delineate the extent of the fjord (green) as those gridcells that are below sea level and reachable from the identified point via flood fill.

5.3. Terminus residuals

We examine the relationship between fjord geometry and glacier retreat due to calving rates, an effect that Slater and others (2019) determined to be secondary to ocean warming. In order to see this effect in the data, it is essential to remove the dominant effect of ocean warming. This takes place in two steps: calibration and computation.

5.3.1. Calibration and computation

This study describes observed terminus state based on data from Wood and others (2021) using up area A _T (Section 5.2); whereas Slater and others (2019) describe observed terminus state as a linear position l _s along a center line. Assuming a constant fjord width w, there is a linear relationship between the two:

(4)$$A_T = w l_{\rm s} + b.$$

We determine w and b empirically via linear regression, where b is an arbitrary constant that depends on zero points chosen for l _s and A _T. These coefficients are then used to convert observed up area A _T to observed l _w, an effective terminus position calibrated to the same scale and offset used in Slater and others (2019):

(5)$$l_{\rm w} = ( A_T-b) /w.$$

Using Eqns (3) and (5), the predicted terminus position l _p (Slater and others, 2019) based solely on observed advection vs increases in frontal melt due to ocean warming may be compared to the observed terminus position l _w, which is based on all processes affecting terminus position (Eqn (1)). We compute the terminus residual $l_{\epsilon }$ (Fig. 6) as:

(6)$$l_{\epsilon} = l_{\rm p} - l_{\rm w} = \kappa q_{\rm m} + \beta - {A_T-b\over w}.$$

$l_{\epsilon }$ will be affected by all processes except advection and frontal melt: calving (q _c) and thinning-induced retreat (q _s). The values of q _s computed by Wood and others (2021) and shared in that paper's supplement are at least an order of magnitude smaller than q _c for the glaciers in this study. Therefore to first order, $l_{\epsilon }$ is an estimate of advance or retreat due to decreases or increases in calving.

Figure 6. Computation of the terminus residual for AP Bernstorff glacier. Blue dots: terminus positions as predicted by a thermal forcing model from Slater and others (2019). Annual predictions are available because annual thermal forcing estimates are available; however, note that the Slater model coefficients are determined based on regressions involving 5 year averaged data. Orange plusses: terminus positions based on up area calculated from termini in Wood and others (2021) and calibrated to terminus positions from Slater and others (2019). Black lines: The terminus residual is the difference between the two predictions. The increasingly negative terminus residual means the glacier is retreating faster than Slater and others (2019) would predict based on thermal forcing alone, indicating a destabilizing influence of fjord geometry. The Fjord Map for this glacier (Fig. 15) confirms that runaway retreat is well underway.

Observations show that some glaciers have been stable since 2000, for example Rink Isbræ and Sermeq Avannarleq (Fig. 14). Even though they have been stable overall, their termini have still advanced or retreated by up to 600 m over the study period, where total retreat is computed based on a least squares fit through the annual terminus locations. Total retreat of <600 m over the study period of 1980–2020 is considered not significant because that is within the common range of natural variability for stable glaciers in this study. It is hypothesized that in the face of continued ocean warming, these glaciers might destabilize in the future. This study is unable to test that hypothesis because by design it contains no forward modeling component, and it only collects information when glacier termini move significant distances. Note that most glaciers in our study that are retreating today only began to do so ~2000 (Murray and others, 2015).

6. Classification by TGC stage

The current stage in the TGC for an individual glacier may be evaluated by computing the calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ based on the von Mises calving law, and regressing it against the terminus residuals (Section 5.3). This is developed as follows:

• • Sections 6.1–6.3: The von Mises tensile stress $\tilde {\sigma }$ (von Mises, 1913; Morlighem and others, 2016; Choi and others, 2018) is computed at every point on the glacier surface.

• • Section 6.4: $\tilde {\sigma }$ is integrated over the glacier terminus T to obtain $\sigma _{\scriptscriptstyle T}$, the von Mises stress at the terminus.

• • Sections 6.5 and 6.6: $\sigma _{\scriptscriptstyle T}$ is averaged across velocity fields of different years to obtain a single value $\bar {\sigma _{\scriptscriptstyle T}}$ for each year's terminus.

• • Section 6.7: $\bar {\sigma _{\scriptscriptstyle T}}$ is regressed against terminus residuals to determine whether each glacier stabilizes or destabilizes when it begins retreating.

6.1. von Mises calving law

The von Mises calving law (von Mises, 1913; Morlighem and others, 2016; Choi and others, 2018) predicts a glacier will calve when the tensile von Mises stress $\tilde {\sigma }$ at the terminus exceeds the ice's yield strength $\sigma _{\it \scriptsize \hbox {max}}$. The calving rate q _c is given by Morlighem and others (2016):

(7)$$q_{\rm c} = \displaystyle{\tilde{\sigma} \over \sigma_{\scriptsize \it \hbox{max}}} \vert \vert \vec{u} \vert \vert ,\; $$

where $\vec {u}$ is the vertically averaged horizontal velocity. We assume plug flow near the calving front (Greve and Blatter, 2009; Bassis and Ultee, 2019), making the vertically averaged velocity equal to surface velocity.

6.2. Tensile von Mises stress

The von Mises calving law requires computation of the tensile von Mises stress. In continuum mechanics, the strain rate tensor $\dot {\epsilon }$ may be computed from the gradient of the velocity $\vec {u}$ as:

(8)$$\dot{\epsilon} = {1\over 2} \left(\nabla \vec{u} + \nabla \vec{u}^T\right),\; $$

where $\nabla \vec {u}^T$ is the transpose of the rank 2 tensor $\nabla \vec {u}$. (See Gibbs and Wilson (1901), page 404, eqn 3, also Cajori (1928), volume II, page 135.) The scalar tensile strain rate $\bar {\dot {\epsilon }}$ (Morlighem and others, 2016) is defined as:

(9)$$\bar{\dot{\epsilon}}^2 = {1\over 2}\left(\max( 0,\; \dot{\epsilon}_1) ^2 + \max( 0,\; \dot{\epsilon}_2) ^2\right),\; $$

where $\dot {\epsilon }_1$ and $\dot {\epsilon }_2$ are the eigenvalues of $\dot {\epsilon }$. Glen's flow law, the constitutive relation used to model ice deformation and flow, relates the strain rate tensor $\dot {\epsilon }$ to the deviatoric or shear stress tensor σ:

(10)$$\dot{\epsilon} = \tilde{A} \sigma^n,\; $$

where $\tilde {A}$ is the temperature-dependent rate factor (s⁻¹ Pa⁻ⁿ), and n is typically assumed to be 3 (Behn and others, 2021). In this case (Morlighem and others, 2016), Glen's flow law is used with the scalar tensile strain rate $\bar {\dot {\epsilon }}$, and solved for the scalar tensile von Mises stress $\tilde {\sigma }$, to obtain:

(11)$$\tilde{\sigma} = \sqrt{3} \tilde{B} ( {\bar{\dot{\epsilon}}})^{1/n},\; $$

where $\tilde {B} = \tilde {A}^{-1/n}$ is the ice hardness (Greve and Blatter, 2009).

Figure 7 shows the von Mises stress computed on a grid for one velocity field. Disregarding processes other than calving for now, the von Mises calving law predicts that advancing glaciers will have $\tilde {\sigma } < \sigma _{\it \scriptsize \hbox {max}}$ and retreating glaciers will have $\tilde {\sigma } > \sigma _{\it \scriptsize \hbox {max}}$. As a catch-all parameter, $\sigma _{\it \scriptsize \hbox {max}}$ accounts not just for ice cliff properties and fjord geometry but all factors affecting calving, for example frozen melange in the fjord (Schlemm and Levermann, 2020).

Figure 7. (a) von Mises tensile stress $\tilde {\sigma }$ shown for Kangilleq and Sermeq Silarleq as computed by the PISM, based on a sample velocity field from 2018. (b) Ice velocity vectors and sample terminus (red line), used in conjunction with $\tilde {\sigma }$ to obtain calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$. Ice velocities downstream of the terminus do not reflect grounded ice, they could be an ice shelf or ice melange.

6.3. Computing von Mises stress

For each surface velocity map, we use the parallel ice sheet model (PISM; Khroulev and PISM Authors, 2022) to compute the tensile von Mises stress $\tilde {\sigma }$ for a given ITS_LIVE velocity field, using the PISM default constant ice hardness of $\tilde {B} = {68\, 082}\, {\hbox{Pa}\, \hbox{s}^{1/3}}$.

6.4. Integrating across the terminus

To obtain a single von Mises stress number for a glacier, the von Mises map computed in Section 6.3 is integrated across the glacier's terminus. The value $\sigma _{\scriptscriptstyle T}$ is defined as the average von Mises stress across the glacier terminus for an entire terminus line T of arbitrary shape:

(12)$$\sigma_{\scriptscriptstyle T} = { \displaystyle\int_T ( \tilde{\sigma} \vec{u} \cdot \vec{n}\,) {\rm d}s\over \displaystyle\int_T ( \vec{u} \cdot \vec{n}\,) {\rm d}s},\; $$

where $\vec {n}$ is the unit normal and ds is used for the line integral along T, using a rasterized terminus and a raster-based formulation for the line integral (Appendix A).

This definition of $\sigma _{\scriptscriptstyle T}$ is robust to missing velocity data near the edges of fast glacier flow and near the terminus, a common situation when using remote-sensing ice surface velocity data. If $\vec {u}$ is missing at some points along L, then the line integrals in the numerator and denominator will both be missing at the same points, and will therefore avoid biasing the result to first order. In this way, $\sigma _{\scriptscriptstyle T}$ is normalized by the amount of data that can be measured (Fig. 8). Because of missing data near the margins, the value of $\sigma _{\scriptscriptstyle T}$ depends more heavily on what is happening in the center of the fjord.

Figure 8. Aerial map of AP Bernstorff Glacier in Southeast Greenland showing incomplete data for ice velocities that happen in some cases. Annual ITS_LIVE velocity data within the fjord are overlaid on bedmap elevation and fjord bathymetry. Ice velocity data are not shown outside the fjord, where bedrock is above sea level. Terminus measurements within the year are shown in red, with three termini available in 1990, and just one each in 1996 and 2005. Velocity data coverage is sometimes incomplete, especially close to the terminus or near the margins of the glacial trough. Line integrals in this study disregard any portion of the terminus with missing data. Although the equation for $\sigma _{\scriptscriptstyle T}$ is robust to missing data at the terminus, it can still fail for lack of data, as in 1996.

Annual terminus lines from Wood and others (2021), manually digitized from LANDSAT 5, 7 and 8 imagery, were rasterized on the MEaSUREs grids and used for this analysis. By reusing data from Wood and others (2021), this study maximizes the ability to compare results with other recent work; however, it is also limited to glaciers included in that study.

In theory a 1-D calving rate q _c can be estimated directly by using $\sigma _{\scriptscriptstyle T}$ for $\tilde {\sigma }$ in Eqn (7). However, errors in estimating $\sigma _{\it \scriptsize \hbox {max}}$ lead to large uncertainties in the actual value of q _c, which is not needed anyway. Instead, the von Mises calving law suggests that $\sigma _{\scriptscriptstyle T}$ on average should be proportional to calving rate q _c. Even without knowing the coefficient of proportionality, this allows $\sigma _{\scriptscriptstyle T}$ to be used as a proxy for q _c without ever having to explicitly determine $\sigma _{\it \scriptsize \hbox {max}}$.

6.5. Stacking to obtain calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$

Change in surface velocity, not terminus position, is the dominant driver of annual variation in $\sigma _{\scriptscriptstyle T}$ (Fig. 9). To single out the effect of the position of the terminus in the fjord rather than surface velocity, $\sigma _{\scriptscriptstyle T}$ is computed using multiple velocity fields for each terminus, even if the terminus and velocity field are from different years. The result is then averaged to create $\bar {\sigma _{\scriptscriptstyle T}}$. For this procedure to work, there must be ice at the terminus so that $\sigma _{\scriptscriptstyle T}$ can be computed; which for retreating glaciers means the velocity field must pre-date the terminus position.

Figure 9. Calving proxy $\sigma _{\scriptscriptstyle T}$ value computed for one glacier (Hayes N); plotted by velocity year (year of the velocity field used) and terminus year (year of the terminus used), where the velocity year is always less than the terminus year. Although $\sigma _{\scriptscriptstyle T}$ varies due to the position of the terminus, the largest variation usually occurs due to changes in the overall ice velocity field: some years a glacier may be moving faster, whereas other years it may be moving more slowly. $\sigma _{\scriptscriptstyle T}$ is averaged across velocity fields of different years to obtain a single value $\bar {\sigma _{\scriptscriptstyle T}}$ for each year's terminus.

Most glaciers in this study were relatively stable until ~2000, after which they began to retreat en masse (Murray and others, 2015). Due to limited availability of data and the need for surface velocities to pre-date terminus positions, only the post-2000 period of retreat is studied. Therefore, only terminus/velocity pairs were used in which the terminus year was 2001 or later, and the velocity year was older than the terminus year.

6.6. Effect of surface velocity data quality on calving proxy

This study uses many approaches to be robust in spite of missing surface velocity data. It uses a robust integration technique along the terminus (Section 6.4), and then it uses an averaging technique analogous to stacking (Section 6.5), a well-established technique in seismology used to improve the signal-to-noise ratio of data. Each terminus line is used to integrate all the available velocity fields older than it, thereby decreasing the effect of a poor velocity field from any single year. Figure 9 shows each terminus behaving similarly no matter which velocity field it is applied to, adding confidence that poor quality velocity fields with missing data are not overwhelming the signal. Finally, termini are only applied to older velocity fields. Because most glaciers are retreating, this means that the newer terminus will typically be somewhat upstream of the end of an older velocity field and will likely be sampling an improved portion of that velocity field.

6.7. Analysis of $\bar {\sigma _{\scriptscriptstyle T}}$ and terminus residual

The terminus residual $l_{\epsilon }$ represents the amount of terminus advance/retreat that is not explained by thermal forcing alone (Section 5.3). With $l_{\epsilon }$ and $\bar {\sigma _{\scriptscriptstyle T}}$ it is now possible to evaluate whether the calving rate proxy $\bar {\sigma _{\scriptscriptstyle T}}$ increases or decreases as the glacier retreats. $l_{\epsilon }$ and $\bar {\sigma _{\scriptscriptstyle T}}$ are regressed against each other with a p-value significance threshold of 0.21 (see Section 8.5):

(13)$$\Delta l_{\epsilon} = \nu \Delta \bar{\sigma_{\scriptscriptstyle T}},\; $$

where ν is the regression coefficient indicating whether fjord geometry causes $\bar {\sigma _{\scriptscriptstyle T}}$ to increase or decrease as the glacier retreats.

If a glacier is susceptible to rapid retreat and has just begun to retreat, then the TGC predicts ν should be negative. That is, stress at the terminus $\bar {\sigma _{\scriptscriptstyle T}}$ increases as the glacier begins to retreat, causing a positive feedback that could lead to instability. Such a glacier could well continue to retreat, even if frontal melt rate were to stabilize or decrease. If on the other hand a glacier is in a stable configuration, then ν will be positive, meaning $\bar {\sigma _{\scriptscriptstyle T}}$ decreases as the glacier retreats. Such glaciers could be retreating in spite of their geometric stability due to the primary forcing from warming oceans; however, at this time the fjord geometry is helping stabilize the glacier and prevent runaway retreat.

If a glacier has already begun rapid retreat and is currently retreating through an area with little variation in fjord cross-sectional geometry, then $\bar {\sigma _{\scriptscriptstyle T}}$ will be about constant, even as $l_{\epsilon }$ changes. There is no relationship between $\bar {\sigma _{\scriptscriptstyle T}}$ and the terminus residual $l_{\epsilon }$, and hence the p-value value for ν will be high. Lack of statistical significance is correlated with glaciers already in rapid retreat, as was confirmed in our results. Note that in principle, lack of predictive power of the Slater regression must also be considered as a possible cause.

7. Validation of von Mises calving law

The von Mises calving law as a model may be validated by applying it to the data of Wood and others (2021) and evaluating the result for coherence and consistency. Wood and others (2021) measure or model all terms of the mass balance equation (Eqn (1)) except for the calving rate q _c – which is computed as a residual between observed terminus location L and the integrated effect of all other fluxes: ice advection ($q_{\rm f} = \Vert \vec {u}\Vert$), frontal melt (q _m) and thinning-induced retreat (q _s). Although $\sigma _{\it \scriptsize \hbox {max}}$ does not need to be computed for this study, it may still be determined from the Wood and others (2021) data and the definition of the von Mises calving law (Eqn (7)):

(14)$$\sigma _{\it \scriptsize \hbox {max}} = {\sigma_{\scriptscriptstyle T}\over q_{\rm c}} {{\rm d}L\over {\rm d}t},\; $$

where dL/dt is the rate of terminus advance or retreat.

Using this formula, $\sigma _{\it \scriptsize \hbox {max}}$ was estimated based on all available terminus lines of Wood and others (2021), using the velocity field from the same year as each terminus line. Figure 10a shows the result grouped by glacier. In this plot, $\sigma _{\it \scriptsize \hbox {max}}$ displays a two-tailed cumulative distribution function. This is to be expected for a value like $\sigma _{\it \scriptsize \hbox {max}}$ that is thought to be affected by a number of glacier-specific parameters such as ice shelves, melange characteristics, etc.; and would therefore be expected to converge on a normal distribution. Figure 10b shows that the average of $\sigma _{\it \scriptsize \hbox {max}}$ across all glaciers does not change much from year to year. These results are reasonable and coherent, even though q _c is a residual, and therefore incorporates all errors and biases from the various datasets and models used by Wood and others (2021). Overall estimated value for $\sigma _{\it \scriptsize \hbox {max}}$ is 300 ± 100 kPa.

Figure 10. Implied $\sigma _{\it \scriptsize \hbox {max}}$ parameter obtained by fitting $\sigma _{\scriptscriptstyle T}$ computed using same-year velocity and terminus measurements, to calving rate obtained by residuals of other quantities from Wood and others (2021) (Eqn (14)), and grouped by either glacier or year. The red line is the median, the box extends to the edge of the first and third quartiles, the whiskers extend to the furthest data point in the first and third quartiles and outliers are not shown. (a) $\sigma _{\it \scriptsize \hbox {max}}$ grouped by glacier. For most glaciers, $\sigma _{\it \scriptsize \hbox {max}}$ lies in the range 250–350 kPa, with some outliers. Occasional negative values of $\sigma _{\it \scriptsize \hbox {max}}$ are non-physical and caused by issues with Wood and other's data: $\sigma _{\scriptscriptstyle T}$ is always positive. Consistent value across most glaciers supports von Mises calving law as a reasonable model. (b) $\sigma _{\it \scriptsize \hbox {max}}$ across all glaciers grouped by year. Consistent year-to-year stability supports von Mises calving law as a reasonable model.

In some cases, $\sigma _{\it \scriptsize \hbox {max}}$ is estimated to be negative. That is a limitation of the Wood and others (2021) dataset and the nature of q _c as a residual, since $\sigma _{\scriptscriptstyle T}$ used in Eqn (14) by definition is always positive. The overall consistency of $\sigma _{\it \scriptsize \hbox {max}}$ suggests that the von Mises calving law is a useful model for Greenland tidewater glaciers; and that residual values of q _c (Wood and others, 2021) are not overwhelmed by noise, in spite of the multi-step process used to compute them. Overall, our results (Fig. 10) support the von Mises criterion as a reasonable calving model for Greenland tidewater glaciers.

8. Results and discussion

Of the 44 glaciers analyzed, 10 were determined to have a destabilizing fjord geometry (the glacier is calving more as it retreats), 7 a stabilizing fjord geometry (the glacier is calving less as it retreats), 16 were found to be stable so far (their termini have not moved much in the dataset) and 11 to have already entered the rapid retreat phase of the TGC (Figs 4 and 11). Each category is analyzed further below.

Figure 11. Glacier categorization flowchart. Glaciers that have moved <600 m over the study period are considered stable so far. Otherwise, a regression between the calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ and the terminus residual $l_{\epsilon }$ is performed. If that regression lacks significance at p-value of 0.21, then the glacier is considered to already be in rapid retreat. Otherwise, the sign of the regression coefficient ν distinguishes between destabilizing geometry (negative sign) vs stabilizing geometry (positive sign).

8.1. Destabilizing fjord geometry

Some glaciers have a regression coefficient ν < 0 (negative slope of line in column (b)), suggesting that they are entering the rapid retreat stage of the TGC. Their termini have retreated substantially (more than 600 m) since 2000; and they have retreated faster than thermal forcing would have predicted. Puisortoq N and Puisortoq S (Fig. 12) in Southeast Greenland are both canonical examples of retreat that has continued due to fjord geometry in spite of recent decreases in thermal forcing, suggesting that the retreat has become self-sustaining. Carlos Glacier on the west coast shows a similar pattern. Some glaciers show episodic retreat; for example, Eqip Sermia on the west coast. In this case, the episodic retreat is correlated to changes in thermal forcing, although it could also be due to pinning points.

Figure 12. Analysis of glaciers that destabilize upon retreat. (a) 5 year Slater relative terminus (blue) and melt (green crosses) used in Slater regressions; and annual Wood relative terminus (orange). Slater (blue) and Wood (orange) relative termini should be similar because they measure the same physical quantity. Predictions from the Slater thermal forcing model are not shown. (b) Regression of calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ vs relative terminus residuals as per Slater. (c) Reference map of fjord. Although thermal forcing has decreased since 2015, retreat has continued. Based on fjord geometry and recent decreases in retreat rate, Puisortoq N and Eqip Sermia might stabilize soon; however, that is speculation because the terminus has not yet had a chance to ‘see’ these potential pinning points, and thermal forcing could cause continued retreat in any case.

8.2. Stabilizing fjord geometry

Some glaciers have a stabilizing fjord geometry. This category is expected to be small because a glacier must have stabilizing fjord geometry but still be retreating anyway due to frontal melt, a condition that would happen near the end of the rapid retreat phase. Lille Glacier (Fig. 13) is a good example of this, as the terminus retreats into a narrow section at the head of the fjord. Ussing Braeer N may also fall into this category, although its geometry is more complex. This increasing stability and slowing down of retreat is the ultimate fate for many tidewater glaciers because fjords must become shallower close to their head, or narrower at a pinning point. In the past, glaciers may have come to rest long term at pinning points, but continually rising ocean temperatures make that less likely in the future. In this study we observe many glaciers slowing their retreat at pinning points; but we see no evidence of tidewater glaciers stabilizing anywhere but on land once they have begun rapid retreat. Some glaciers, for example Hayes North, have complex geometry and are a poor ‘fit’ for a linear regression.

Figure 13. Analysis of glaciers that stabilize upon retreat. (a) 5 year Slater relative terminus (blue) and melt (green crosses) used in Slater regressions; and annual Wood relative terminus (orange). Slater (blue) and Wood (orange) relative termini should be similar because they measure the same physical quantity. Predictions from the Slater thermal forcing model are not shown. (b) Regression of calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ vs relative terminus residuals as per Slater. (c) Reference map of fjord. Kujalleq's terminus has not moved enough to adequately sample changes in fjord geometry. And from the map, it Lille now terminates near the head of the fjord, where water becomes more shallow with further retreat.

8.3. Currently stable

Some glaciers have been stable during the study period, with termini that moved on average <600 m: the methods of this study revealed no new information about them, beyond their already known recent stability. Because the terminus stayed relatively stationary, no statistically significant relationship was found between terminus residual and $\bar {\sigma _{\scriptscriptstyle T}}$ (Fig. 14). The complete list of glaciers in this class is Anorituup Kangerlua N, Daugaard Jensen, Hayes M, Kangerlussup, Kangiata Nunaata, Kangilleq, Koge Bugt S, Rink Isbræ and Sermeq Avannarleq (Table 2).

Figure 14. Analysis of glaciers for which a least square fit of terminus position has retreated <600 m over the study period; and due to lack of sampling from terminus movement, were statistically insignificant. (a) 5 year Slater relative terminus (blue) and melt (green crosses) used in Slater regressions; and annual Wood relative terminus (orange). Slater (blue) and Wood (orange) relative termini should be similar because they measure the same physical quantity. Predictions from the Slater thermal forcing model are not shown. (b) Regression of calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ vs relative terminus residuals as per Slater. (c) Reference map of fjord.

To account for noisy data, the threshold for retreat was determined based on the slope of the least squares fit line through the terminus positions of Wood and others (2021) since 2000, implicitly assuming a constant retreat rate since that time. This procedure works in most cases, but it can yield erroneous results when retreat rate has not been constant. For example, Daugaard Jensen was historically considered stable, with concern it could soon destabilize (Bevan and others, 2012). Examination of the data suggests it advanced a modest 700 m before 2013, and since then has retreated almost 1 km. Daugaard Jensen is no longer, in fact a stable glacier – it is currently retreating. However this study erroneously classifies it currently stable because the overall retreat since 2000 has been modest. More sophisticated statistical techniques might be used to overcome this methodological deficiency.

8.4. In retreat

Finally, there are the glaciers for which no statistically significant relationship could be found between terminus residual and $\bar {\sigma _{\scriptscriptstyle T}}$, but have retreated at least 600 m over the study period. This happened for various reasons.

8.4.1. Change of behavior

As above, some glaciers changed behavior during the course of the study, confounding a single linear model. For example, AP Bernstorff (Fig. 15) retreated rapidly until 2005, after which it has remained stable – in spite of changes in thermal forcing both up and down. This is apparently caused by a shallowing of the fjord at the current terminus location. Herluf Trolle S and Mogens Heinesen C are other examples. Ummiammakku retreated rapidly until 2010, at which it stabilized on a pinning point. It is classified as in retreat by the systematic methods of this study because most of the data show it retreating: if it has truly stopped retreating, there have not yet been enough years of stability to statistically ‘overwhelm’ the previous years of retreat. Improvements to the methodology that weight recent behavior more strongly might be able to overcome this kind of limitation.

Figure 15. Glaciers that changed their behavior over the course of the study, confounding the linear model. All four of these glaciers retreated faster in the past but have since stabilized, or begun to stabilize. (a) 5 year Slater relative terminus (blue) and melt (green crosses) used in Slater regressions; and annual Wood relative terminus (orange). Slater (blue) and Wood (orange) relative termini should be similar because they measure the same physical quantity. Predictions from the Slater thermal forcing model are not shown. (b) Regression of calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ vs relative terminus residuals as per Slater. (c) Reference map of fjord.

8.4.2. Retreating steadily

Inngia, Kangerlussuaq, Mogens Heinesen N and Savissuaq WWW (Fig. 16) are retreating steadily through a uniform portion of the fjord, likely driven by thermal forcing and having already retreated off their stable terminal moraine before the start of this study. Fjord geometry does not affect retreat at this point in time because of this uniformity, which results in points without a clear linear relationship when plotted, and thus a lack of statistical significance when regressing for ν (the slope of the line in column (b)). However, they all show negative ν at less than statistical significance, suggesting that small variations in fjord geometry are affecting terminus position as the TGC hypothesis would suggest.

Figure 16. Glaciers that retreated steadily through a uniform portion of the fjord. (a) 5 year Slater relative terminus (blue) and melt (green crosses) used in Slater regressions; and annual Wood relative terminus (orange). Slater (blue) and Wood (orange) relative termini should be similar because they measure the same physical quantity. Predictions from the Slater thermal forcing model are not shown. (b) Regression of calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ vs relative terminus residuals as per Slater. (c) Reference map of fjord.

8.4.3. Complex fjord geometry

Some glacier termini inhabit broad regions of grounded ice without well-defined fjords, often fed by multiple glaciers upstream. In these cases, neither Slater's thermal forcing model nor the TGC seems to show much predictive power: Hayes NN and Uunartit, for example (Fig. 17). Some glaciers do exist in well-defined fjords but the terminus is close to a branching or merging point, for example Savissuaq WW (Fig. 17). Other glaciers simply lack data sufficient to build statistically meaningful results: for example, Wood and others (2021) provide only two terminus positions for Akullikassaap E and three for Anorituup Kangerlua N.

Figure 17. Glaciers with poorly defined or complex fjord geometry. (a) 5 year Slater relative terminus (blue) and melt (green crosses) used in Slater regressions; and annual Wood relative terminus (orange). Slater (blue) and Wood (orange) relative termini should be similar because they measure the same physical quantity. Predictions from the Slater thermal forcing model are not shown. (b) Regression of calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ vs relative terminus residuals as per Slater. (c) Reference map of fjord. Hayes NN and Uunartit exist in broad areas without clear fjord boundaries: the straight lines defining the ‘edges’ of these fjords are edges of the manually drawn polygons and do not represent any actual physical boundary. Savissuaq WW has a well-defined fjord, but complexity arises in this case as the terminus retreats through a branch point.

8.5. Edge cases and outliers

The choice of the threshold at 0.21 to separate glaciers in rapid retreat from ones that are stabilizing/destabilizing (Fig. 11) is somewhat arbitrary. Some glaciers show clear and consistent behavior and have small p-values, for example Puisortoq N (Fig. 12) and Lille (Fig. 13). Other glaciers show large p-values indicating no effect of changes in fjord geometry on continued retreat (Fig. 16). However, it is harder to classify the behavior of glaciers with p-values close 0.21.

Figure 18 shows the glaciers of highest p-value in each of the destabilizing and stabilizing categories, and the glacier of lowest p-value in the in retreat category (AP Bernstorff). All three of these glaciers are correctly classified, but are also edge cases for their categories, as evidenced by their marginal p-values.

• • Kangilernata is destabilizing; but it is retreating off an unusually broad shoal more than 2 km wide, creating a situation, similar to that of ongoing rapid retreat, in which the fjord cross section does not change much even as the glacier has retreated more than 2 km.

• • Uunartit is stabilizing as it reaches the end of the fjord; but the scenario is confounded because this fjord gets deeper even as it narrows, thereby reducing for now the amount that calving decreases as it retreats.

• • AP Bernstorff is in the mid of rapid retreat; however, retreat has recently slowed down as it has reached a shallower section of the fjord, resulting in overall more stabilizing behavior. This shows how changes in behavior over the study period can confound the methods of this study.

Figure 18. Edge case glaciers within each category: Kangilernata and Uunartit are destabilizing and stabilizing, respectively, and have the highest p-values in their classifications. AP Bernstorff is in retreat and has the lowest p-value in its classification. (a) 5 year Slater relative terminus (blue) and melt (green crosses) used in Slater regressions; and annual Wood relative terminus (orange). Slater (blue) and Wood (orange) relative termini should be similar because they measure the same physical quantity. Predictions from the Slater thermal forcing model are not shown. (b) Regression of calving proxy $\bar {\sigma _{T}}$ vs relative terminus residuals as per Slater. (c) Reference map of fjord.

Full results and classifications are provided in the Supplementary material, allowing the reader to compare other glaciers to these edge cases and to evaluate the potential effect of other p-value cutoffs. Although the cutoff value p = 0.21 provided accurate classifications in this case, we do not expect p = 0.21 to be fundamental to this method. The data sources in this study came with large and often unquantifiable uncertainties. If those were to be reduced in a future study, we would expect a smaller p-value cutoff to be appropriate.

In some cases, outliers can cause the regression to mis-classify. For example, Mogens Heinesen S (Fig. 19) is classified as stabilizing, but the regression data suggest it is actually destabilizing, except for two outlier data points from the 1980s.

Figure 19. Mogens Heinesen S, which is mis-classified due to two outlier points in the regression of $\sigma _{T}$ vs residuals (column b). (a) 5 year Slater relative terminus (blue) and melt (green crosses) used in Slater regressions; and annual Wood relative terminus (orange). Slater (blue) and Wood (orange) relative termini should be similar because they measure the same physical quantity. Predictions from the Slater thermal forcing model are not shown. (b) Regression of calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$ vs relative terminus residuals as per Slater. (c) Reference map of fjord.

9. Future work

This study offers encouraging preliminary results that could be improved in many ways: more glaciers in the study, more data per glaciers, more advanced machine-learning techniques, and more predictor variables. Lack of satellite data before 2000 is a persistent issues limiting the statistical techniques available.

Although we examined glaciers systematically in this study, only 44 of the over 200 Greenland tidewater outlet glaciers (Fahrner and others, 2021) were included, a consequence of relying on multiple previous studies for data. The limiting factor was the requirement that glaciers appear in both the Wood and Slater datasets, and also on a MEaSUREs grid. Although Wood and others (2021) provide data on the different factors driving terminus retreat, ultimately the only portion of that dataset used was the terminus lines. Recent efforts have produced abundant terminus traces through machine-learning techniques (Cheng and others, 2021b), which could in principle allow these methods to be run on more data and a larger set of glaciers. It might be possible to reduce noise by ‘stacking’ results of high-frequency terminus measurements from different seasons within 1 year. Similarly, it might be possible to use more than one velocity field per year.

Improving the data analysis of this study is another avenue for future research. The current study uses two sequential linear regressions: first the regression of Slater and others (2019), and then a regression of $\bar {\sigma _{\scriptscriptstyle T}}$ on terminus residuals. More typically, multiple linear regression would be used here. The use of sub-annual termini could add data for more robust statistics, but would also introduce more natural seasonal variability in terminus position, which would have to be accounted for; there is no obvious way to take an ‘average’ of multiple terminus lines. Recent efforts have produced abundant terminus traces through machine-learning techniques such as automated deep learning (Cheng and others, 2021b), which would in principle enable a larger number of glaciers for a study like this.

If there are enough data to support them, advanced machine-learning techniques could be applied to predict terminus position based on a range of predictor variables: subsurface runoff (Q), ocean/fjord thermal forcing (TF), $\bar {\sigma _{\scriptscriptstyle T}}$, air temperature and other climate drivers (Fahrner and others, 2021). These methods might be a reasonable way to use high-frequency (sub-annual) terminus lines and velocity fields.

The observational methods in this study rely on large amounts of data and are only applicable for the satellite era. For the study of tidewater glacier behavior in the past or future, modeling based on dynamic ice models such as PISM ( Khroulev and PISM Authors, 2022) bounded by Bed Machine (Morlighem and others, 2017) would be more appropriate. Although observations from before the satellite era are too sparse to use for the methods in this study, they would be invaluable in calibrating and validating physics based models, opening a window into the past.

Although this study focuses on Greenland only, it does not rely on any properties specific to Greenland; and given appropriate datasets, we believe its methods can be generalized to tidewater glaciers worldwide. Given appropriate data, these methods could help provide a stability assessment for tidewater glaciers in other regions such as Alaska and Antarctica.

10. Conclusions

Using the calving proxy $\bar {\sigma _{\scriptscriptstyle T}}$, we quantitatively identify Greenland tidewater glaciers for which fjord geometry is either adding to or detracting from terminus stability, and qualitatively match to expectations based on a visualization of fjord geometry. By showing general agreement to predictions, we provide quantitative support of the TGC as a useful model for understanding the behavior of tidewater glaciers in Greenland. Based on data from Wood and others (2021) and estimates of $\sigma _{\it \scriptsize \hbox {max}}$ from that data, we support the von Mises criterion as a reasonable calving model for Greenland tidewater glaciers (Fig. 10), with $\sigma _{\it \scriptsize \hbox {max}} \approx 300 \pm 100\, \hbox{kPa}$.

We confirm the general assertions of Wood and others (2021) and Slater and others (2019), that increased frontal melt due to ocean warming since 2000 is currently the dominant process driving tidewater glacier retreat in Greenland today. This dominant effect must be removed from the data in order to study calving dynamics and rapid retreat controlled by fjord geometry. Because frontal melt has only recently become dominant over calving for tidewater glacier retreat in Greenland (due to ocean warming), early work does not address ocean warming and instead focuses on calving as the primary mechanism of retreat, and does not address ocean warming (Post, 1975; Meier and Post, 1987).

Although statistically significant in many cases, the linear relationship between ocean thermal forcing and tidewater glacier retreat as developed by Slater and others (2019) should be used with caution because it does not account for the calving effects of fjord geometry inherent in the TGC. The linear relationship would suggest tidewater glaciers behave like land-terminating glaciers, advancing and retreating in lockstep with climate, which runs contrary to our understanding of the TGC (Pfeffer, 2007). For this reason, we suggest caution in using Slater and others (2019) to generate future extrapolated boundary conditions for a general circulation model, as was proposed in that study.

In spite of increasing frontal melt, not all Greenland glaciers are retreating. We hypothesize this is due to exceptionally stabilizing fjord geometry, which the methods of this study are unable to confirm or deny. Speculation on the future of currently stable glaciers might best be accomplished through modeling studies based on the measured bed geometry and idealized thermal forcings and frontal melt.

A number of glaciers confound the methods presented here. Some lack statistical significance for glaciers with complex bed geometries or ill-defined fjords. Some transition between regimes over time – either increasing or decreasing retreat rate quickly as in a surge-type glacier. These issues are problems in the current analysis, which is based on simple linear regressions with an assumption of stationarity. However, the satellite era for glaciers is short and overall lack of data could render useless efforts to use more powerful machine-learning techniques, which would require large datasets.

In spite of the complexity, the TGC is consistently supported by the evidence in this systematic study of glaciers. Glaciers that retreat faster than thermal forcing models would predict have increasing $\bar {\sigma _{\scriptscriptstyle T}}$ with retreat; and in these cases, the terminus is observed to be retreating through a section of the fjord that is widening or deepening, thereby generally confirming the TGC. However glaciers with less retreat than thermal forcing would show decreasing $\bar {\sigma _{\scriptscriptstyle T}}$, which can often be verified by observing pinning points, confirming the TGC as well.