Glacier geometry and ﬂow speed determine how Arctic marine-terminating glaciers respond to lubricated beds

. Basal conditions directly control the glacier sliding rate and the dynamic discharge of ice. Recent glacier destabilization events indicate that some marine-terminating glaciers quickly respond to lubricated beds with increased ﬂow speed, but the underlying physics, especially how this vulnerability relates to glacier geometry and ﬂow characteristics, remains unclear. This paper presents a 1D physical framework for glacier dynamic vulnerability assuming sudden basal lubrication as an initial perturbation. In this new model, two quantities determine the scale and the areal ex-tent of the subsequent thinning and acceleration after the bed is lubricated: Péclet number ( Pe ) and the product of glacier speed and thickness gradient (dubbed J 0 in this study). To validate the model, this paper calculates Pe and J 0 using multi-sourced data from 1996 to 1998 for outlet glaciers in the Greenland ice sheet and Austfonna ice cap, Svalbard, and compares the results with the glacier speed change during 1996/1998–2018. Glaciers with lower Pe and J 0 are more likely to accelerate during this 20-year span than those with higher Pe and J 0 , which matches the model prediction. A combined factor of ice thickness, surface slope, and initial ﬂow speed physically determines how much and how fast glaciers respond to lubricated beds in terms of speed, elevation, and terminus change. marine-terminating glaciers in the GrIS and the Austfonna ice cap, Svalbard respectively. The results show that Pe /(cid:96) and J 0 correlate with the ﬂow speed change during 1996/1998 and 2018, especially for nonsurge glaciers and glaciers without signiﬁcant terminus retreat, matching the model predic-tion. Glaciers with thick ice and a fast ﬂow result in low Pe and negative J 0 , and reduced basal friction leads to initial speedup and thinning, which can propagate further inland via diffusion. For glaciers in the GrIS subject to ocean–ice inter-actions, this new model indicates multiple feedback cycles that make glaciers more sensitive to changing basal conditions. Finally, this study highlights glaciers classiﬁed as vulnerable to lubricated beds (low Pe and high | J 0 | ) but with no signiﬁcant speed change during the past two decades. Fre-quent monitoring is suggested for these because they might be more prone to future instabilities and affect the pro-jected sea level


Introduction
Marine-terminating glaciers worldwide have undergone significant acceleration, retreat, and mass loss in past decades (e.g., Vaughan et al., 2013;Cook et al., 2016;Carr et al., 2017;Catania et al., 2020;Williams et al., 2021). At the Greenland ice sheet (GrIS), the dynamic discharge of marine-terminating glaciers accounts for 66 % of the region's total mass loss (Mouginot et al., 2019). For the other Arctic regions outside of the GrIS, surface mass balance contributes more mass loss than dynamic discharge (Catania et al., 2020), but several rapid acceleration events also dominate the local land ice budget (e.g., McMillan et al., 2014;Strozzi et al., 2017a;Willis et al., 2018;Haga et al., 2020).
The acceleration and dynamic thinning of marineterminating glaciers have been attributed to at least two sources: basal lubrication driven by surface melt accessing the bed, and terminus perturbation driven by ice-ocean interactions (Carr et al., 2013). Multiple observations suggest the warming of subsurface ocean waters as the primary and widespread driver across the outlet glaciers in the GrIS (e.g., Nick et al., 2009;Walsh et al., 2012;Tedstone et al., 2013;Catania et al., 2020;Wood et al., 2021;Williams et al., 2021). As a result, models of glacier dynamic loss for estimating sea-level rise usually focus on the glacier terminus and overlook the changing basal conditions (e.g., Nick et al., 2013). Outside of the GrIS, the primary drivers of the dynamic ice loss remain largely uncertain (Carr et al., 2017;Strozzi et al., 2017b), although significant melt-induced lubrication and speedup events have been identified around the Arctic (e.g., Sundal et al., 2011;Dunse et al., 2015;Zheng et al., 2019;Seddik et al., 2019). To date, the response to basal lubrication is mostly studied on a seasonal scale, which links to changing subglacial hydrology within a year (e.g., Zwally et al., 2002;Bartholomew et al., 2010;Sundal et al., 2013;Hewitt, 2013;Rathmann et al., 2017;King et al., 2018;Williams et al., 2020). However, there have been concerns and observations that the bed conditions can evolve and affect glacier Published by Copernicus Publications on behalf of the European Geosciences Union. 1432 W. Zheng: How marine-terminating glaciers respond to lubricated beds dynamics over multiple years. For example, the recent extensive formation of supraglacial lakes in the GrIS can open new moulins that last for years and contribute to long-lived speedups (e.g., Hoffman et al., 2018). During several glacier surge events in Svalbard and the Russian Arctic, initial flow speedups have also created highly crevassed glacier surface and led to further speedups by additional meltwater routing (e.g., Dunse et al., 2015;Strozzi et al., 2017a;Zheng et al., 2019;Sánchez-Gámez et al., 2019). These events potentially alter the basal conditions by allowing meltwater to reach the bed in all seasons, but their interannual impact is still less constrained (Kehrl et al., 2017). In addition, glacier geometry plays a vital role in how a glacier responds to an external perturbation (Carr et al., 2017;Kehrl et al., 2017), but only terminus disruption has been physically well documented and explained (McFadden et al., 2011;Felikson et al., 2017Felikson et al., , 2021. Whether some glaciers are more sensitive to basal lubrication than others owing to their geometry is not well known. To better understand how much and how fast a glacier responds to basal lubrication and its relationship to glacier geometry, this paper presents a physical model with a 1D framework along flowlines formulating the subsequent change (in terms of both glacier speed and surface elevation) after the glacier bed is suddenly lubricated. We use an existing glacier perturbation model (Zheng et al., 2019) and replace the initial thinning condition with a step reduction of basal friction along the glacier channel. This new model identifies key parameters that dominate elevation change rate and ice flow acceleration. Then, using data from the GrIS and Austfonna ice cap, Svalbard, we derive these parameters for each glacier basin and compare them with glacier speed change over 20 years. The entire processing workflows, including data fetching, all calculations, and figure scripts, are available on Github (https://github.com/whyjz/pejzero, Zenodo DOI: https://doi.org/10.5281/zenodo.5641953) and are compiled as a Jupyter Book ready to be cloud executed using the MyBinder service for full reproducibility (Project Jupyter et al., 2018;Executable Books Community, 2020).

Model development
We build the model on the perturbation theory developed by Nye (1963), Bindschadler (1997), Felikson et al. (2017), and Zheng et al. (2019). Our goal in this section is to formulate the change rate of ice elevation ( dH 1 dt ) and ice speed ( dU 1 dt ) after a glacier bed is lubricated permanently. In this new model, basal lubrication is considered as a sudden perturbation without initial elevation change. The variables defined in the model are listed in Table 1.

Perturbation due to a permanent change of basal conditions
We set up a glacier with the following initial values along its 1D flowline profile: speed (U 0 ), thickness (H 0 ), flux (q 0 ), surface slope (α 0 ), and bed friction term (K 0 ). These values vary along the flowline distance x (positive towards downstream) and do not vary with time (i.e., steady state). We assume that the ice is purely sliding on the bed and does not have any internal deformation, that is, The glacier speed can be further represented using the hardbed sliding law (Weertman, 1957): where m is the flow-law constant and is set to 3 in this study. At t = 0, the bed condition changes and the friction term becomes K 0 + K 1 . The second term K 1 denotes the amount of change and is positive for a lubrication scenario. There is no initial elevation change associated with this event. We also assume that this is a one-time change and is uniform over the glacier, so K 1 is a constant. The subsequent change of speed, elevation, flux, and slope are represented as U 1 , H 1 , q 1 , and α 1 respectively. Unlike K 1 , these quantities vary along the glacier channel and over time. Assuming zero local surface mass balance and zero local stress imbalance (e.g., Felikson et al., 2017;Zheng et al., 2019), the conservation of mass can be expressed as Taking the total derivative of q 1 with respect to t yields If we assume a much more gentle slope of the bedrock than that of the ice surface (Felikson et al., 2017;Zheng et al., 2019), the surface slope can also be expressed as the first derivative of ice thickness: Plugging Eqs. (1), (2), (4), and (5) into Eq. (3) yields where C 0 = ∂q 0 ∂H (7) Table 1. Variables used in the perturbation model defined in this study. In the dimension column, L is length and T is time. The prime notation used in the text denotes the partial derivative of a variable with respect to x; for example, H 1 ≡ ∂H 1 ∂x .

Variable Definition Dimension
x Distance along a 1D glacier flowline towards terminus L t Time after perturbation T m Flow law constant None U 0 (x) Glacier speed before perturbation LT −1 U 1 (x, t) Change of glacier speed after perturbation Basal friction term before perturbation L m−1 T −1 K 1 Change of basal friction after perturbation (assumed constant) L m−1 T −1 H 0 (x) Glacier thickness before perturbation L H 1 (x, t) Change of glacier thickness after perturbation L α 0 (x) Glacier slope before perturbation None α 1 (x, t) Change of glacier slope after perturbation None q 0 (x) Flux before perturbation Change of flux after perturbation Péclet number, see Eq. (14) None Characteristic length (length of perturbation) L and D 0 = ∂q 0 ∂α .
Since H 1 (t = 0, x) = 0, where 2.2 J 0 and Péclet number (Pe) Equation (9) indicates that the ratio of K 1 to K 0 and the value of J 0 both determine the initial elevation change rate. In a lubricating scenario, both K 1 and K 0 are positive, and the initial response of dH 1 dt is inversely proportional to J 0 . To relate J 0 to the glacier speed change, we start from the change of flux and assume U 1 >>H 1 . This can be justified by many observations of glacier destabilization as the amount of speed change is usually one to two orders of magnitude higher than the amount of elevation change (e.g., McMillan et al., 2014;Willis et al., 2018). Therefore, As K 1 is a constant, we can derive the glacier speed change using Eqs. (4), (5), (7), (8), and (11): At t = 0, Similar to the elevation change, J 0 is inversely proportional to the initial glacier acceleration. If two glacier beds are lubricated with the same amount of K 1 /K 0 , the glacier with a higher absolute value of J 0 (i.e., |J 0 |) will be more unstable as the initial speed and elevation change rates are higher. Also, from Eq. (9) we can predict a subsequent elevation change after t = 0. At this point, the last three terms in Eq. (6) begin to take part in the elevation change rate. The second term of Eq. (6) represents an exponential decay of the change rate, and the third and the fourth terms indicate advective and diffusive migration of elevation perturbation, respectively. The coefficients of the latter two terms determine the relative strength between advection and diffusion, with the ratio defined as the Péclet number, Pe: where is the length of a perturbation. If Pe is much higher than 0, forward advection will dominate, and any perturbation of ice thickness will only propagate downstream. This prohibits destabilization in the upper stream if thinning or glacier retreat initiates near the terminus. On the other hand, if Pe ∼ 0 or is negative, either diffusion or backward advection takes place, and a thickness perturbation at the terminus can propagate upstream, changing the dynamics of the entire glacier. Hence, we consider a glacier with a low Pe more vulnerable than one with a high Pe.
W. Zheng: How marine-terminating glaciers respond to lubricated beds Combining Eqs. (1), (2), (7), and (8) with Eq. (14), we can express Pe in terms of ice speed, elevation, and surface slope (see Sect. S4 of Zheng et al., 2019, Eqs. 11 to 16 for derivation details): The final assumption we adopt in the model is α 0 = ∂α 0 ∂x ≈ 0 as the estimated value from the data we use in this paper is essentially small. For example, α 0 of the glacier profile shown in Fig. 3 is only around 10 −6 -10 −7 m −1 for the first 100 km, and the last term in Eq. (15) is roughly an order less than the sum of the first three terms. In practice, this assumption might be necessary because α 0 is highly sensitive to local slope change and may not reflect the glacier's overall mechanism to dissipate the perturbation. With this assumption, Eq. (15) can be reduced to: Note we now express the Péclet number as the form of Pe as we do not focus on a particular perturbation length and instead plan to evaluate the general tendency for the ice flow to dissipate any length of perturbations. Compared with the past models, the expression of Pe in this model has an additional term U 0 U 0 , implying its relationship to spatially changing basal conditions. This extra dependency on glacier speed also suggests that Pe is a changing variable and needs to be re-calculated if ice flow speeds up or slows down (see Discussion for more details).
With the same assumption about α 0 , the expression of J 0 (Eq. 10) can be also reduced to: which is proportional to the product of ice speed and the gradient of ice thickness along the flowline. A typical glacier thins toward the terminus, corresponding to a negative H 0 and J 0 . According to Eq. (13), a negative J 0 suggests that when a lubricating scenario takes place, the glacier might speed up to accommodate the change. From Eq. (9) we can see that the glacier will also get thickened (except at the divide) as thicker ice is sliding and replaces thinner ice downstream.
To summarize, two parameters Pe and J 0 are derived from the 1D basal lubrication model. J 0 represents the strength of initial response to basal lubrication, and Pe gives insights into the mode of mass transport after elevation change occurs. Glaciers with a high |J 0 | and a low Pe (∼ 0 or negative) are more vulnerable to basal lubrication as reduced friction can lead to a high initial acceleration and elevation change rate, which will then propagate to the entire glacier via diffusion or negative advection.

Data and methods for validating the model
To test if the model is suitable for evaluating marineterminating glaciers, we derive observed Pe/ and J 0 for outlet glaciers in the GrIS and Austfonna ice cap, Svalbard. These two regions are selected primarily because surface elevation, bed elevation, and glacier speed data necessary for our calculation are publicly available. We compare the results with the NASA MEaSUREs ITS_LIVE glacier velocity  and frontal retreat records  spanning over 20 years and determine if both Pe/ and J 0 are indicative of the vulnerability to basal lubrication.

Greenland ice sheet
We use the data set published with Felikson et al. (2021), which provides well-constrained flowline data for Greenland's 141 marine-terminating glaciers and their branches (187 basins in total, Fig. 1). These glaciers scatter around the ice sheet and provide a diverse sampling over various climate and oceanic factors. The data set contains six primary flowline shapes for each glacier basin, with vertices sampled every 50 m along the flowlines. We use surface elevations at each vertex, sampled from the Greenland Ice Mapping Project (GIMP, Howat et al., 2014). The GIMP surface elevations come from multiple remote sensing sources and are coregistered with elevations acquired by the Ice, Cloud, and land Elevation Satellite (ICESat), thus best representing the ice sheet elevations during 2003-2009. The flowline vertices also contain the glacier bed elevations, sampled from the BedMachine v3 subglacial topography . Although BedMachine v3 uses the source data collected from 1993 to 2016, we assume that the bed elevations are stable over time and can represent any year in that period. To acquire U 0 and glacier speed change at each flowline vertex, we manually sample the annually mosaicked ITS_LIVE glacier speed data from 1998 and 2018, respectively. The ITS_LIVE data are derived from Landsat 4, 5, 7, and 8 images using the autoRIFT feature tracking software Lei et al., 2021). Finally, each vertex of a flowline has the following key parameters: surface elevation, bed elevation, glacier speed in 1998, and speed difference between 1998 and 2018.
We prepare and process the input data for each flowline using the following steps: 1. As the 1998 speed data do not cover the entire ice sheet, we remove flowlines with only 20 speed readings or less from the input list.
2. Locate vertices with NoData speed values along the flowlines and perform a linear interpolation to fill the missing values. We do not extrapolate the glacier speed; therefore, the NoData vertices at both ends of the flowline are still preserved after this step.
3. Remove flowlines with only 280 valid vertices or less from the input list. A valid vertex should contain all key parameters and no NoData Values.
4. To avoid the effect of small sloping change, we smooth the surface elevation, bed elevation, glacier speed data, and their derivatives using the Savitzky-Golay filter with a window size of 251 vertices (12.5 km) (Savitzky and Golay, 1964;Felikson et al., 2021). We do not apply the smoothing filter to data 0-3 km from the terminus owing to insufficient sampling points within the window size. These unfiltered data will not be used for the next step.
5. Derive Pe/ and J 0 along each flowline using Eqs. (17), (16), and parameters representative of the glacier geometry/speed from 1998. As we empirically derive Pe/ and J 0 for each basin and compare them on the same plot, the results will be insensitive to the selected value of m and the sliding law (Felikson et al., 2021).
6. We also need to include glacier retreat in our analysis to better distinguish the relationship between basal lubrication and speedups. Thus, we use the Greenland Marine-Terminating Glacier Retreat Data , which contain temporal evolution of terminus positions derived from Landsat 5, 7, and 8 images for 226 glaciers. We use QGIS to manually measure the terminus retreat between 1988 and 2018 for each glacier at its center flowline.
7. Compare Pe/ and J 0 with the frontal retreat and the speed change between 1998 and 2018.

Austfonna ice cap, Svalbard
We perform the same analysis for the marine-terminating glaciers of Austfonna, a polythermal ice cap located in NE Svalbard. Austfonna is only about 100 km wide and is considered to have a more uniform climate and oceanic factors than the GrIS, but its marine outlet glaciers have exhibited diverse speed histories in the past 20 years (Fig. 2). For instance, the glacier speed of basin-3 (Storisstraumen) increased 45-fold during the past two decades, likely triggered by feedback between summer melt, crevasse formation, and basal lubrication (McMillan et al., 2014;Dunse et al., 2015;Gong et al., 2018). The other surge-type glaciers include basin-1 (Bråsvellbreen) and basin-17 (Etonbreen), and the last surge periods of both glaciers are around 1938 (Schytt, 1969;Hagen et al., 1993;Dowdeswell et al., 2008). The other glaciers of Austfonna do not have a surge history, but many of them (e.g., basin-2, -5, -7, and -10) have also significantly increased the flow speed since 1996, as seen from Fig. 2.
To calculate Pe/ and J 0 , we use the Ice Thickness Models Intercomparison eXperiment (ITMIX) data set (Farinotti et al., 2017), hosted by the International Association of Cryospheric Sciences (IACS). We use the Austfonna DEM from 1996 (Moholdt and Kääb, 2012), velocity from 1995-1996 (Dowdeswell et al., 2008), and ice thickness from 1996 (Dowdeswell et al., 1986;Farinotti et al., 2017), all included in the ITMIX data set. Eight out of 11 marine-terminated glaciers (basin-1, -3, -4, -5, -6, -7, -10, and -17; Fig. 2) are selected for the analysis; the exceptions are basin-2, -8, and -9 owing to their small length roughly equal to the smoothing window. We construct six glacier flowlines based on the 2018 glacier velocity from the ITS_LIVE annual velocity mosaics (Fig. 2b) and sample ITMIX glacier elevation, ice thickness, and glacier speed data every 50 m along each flowline. Then we follow the same workflow for the outlet glaciers in GrIS (see the previous section) and finally compare Pe/ and J 0 with the ITS_LIVE 2018 glacier speeds. On the other hand, glacier 0277 (Alangordliup Sermia, 68.95 • N, 50.22 • W, ∼ 30 km south of Jakobshavn Isbrae; Figs. 1 and 4) is more stable than glacier 0001 as the amount of the speed change in the past two decades is only 0-40 m yr −1 and is constrained at the first 6 km from the terminus. Pe/ ranges from 2 × 10 −4 to 6 × 10 −4 m −1 within the first 10 km, which is roughly 10 times the values from Jakobshavn Isbrae. Also, the slow glacier speed in 1998 directly results in low |J 0 | (only ∼ −10 m yr −1 ) compared with Jakobshavn Isbrae. These results are supportive for glacier 0277 having a stable condition during the study period. Interestingly, the speed change pattern resembles Pe/ at the first 10 km. The glacier might have dealt with frontal or basal perturbation through advection, as indicated by a large Pe.

Pe & J 0 versus glacier speed change
This study focuses on the parameter variations close to the glacier terminus for two reasons. First, crevasses and moulins are more likely to form at the terminus region via hydrofracture than at higher elevations (Poinar et al., 2015). Therefore, additional meltwater routing can alter the basal conditions over multiple years. In addition, the J 0 tends to be small for the upper regions of all glaciers where the ice flow is slow, making this metric less distinctive from one basin to another. Considering the limit of spatial smoothing in the processing workflow (see Sect. 3.1, step 4), we select the data at 3 km from the terminus for the following intercomparison. Owing to the incomplete spatial coverage of ITS_LIVE data from 1998, only 104 out of 187 GrIS glacier basins have valid values of Pe/ and J 0 at this terminus distance (Fig. 1).
Although most of these glaciers have sped up during the 20 years, 26 glaciers slowed down by up to −522 m yr −1 (Fig. 5a). To illustrate Pe/ and J 0 at the glacier front and their varying direction along the flowline, we plot J 0 versus Pe/ using ball-head-pin-like curves. Each curve represents the average J 0 and Pe/ values of one glacier basin at 3-5 km with the head mark located at 3 km, color-coded based on the speed change at 3 km as well (Fig. 5b). Glaciers with low speed change (pale color curves) tend to cluster around the area where J 0 ≈ 0 m yr −1 and Pe/ > 0.0001 m −1 . Most of these curves are near horizontal on the plot, indicating a small change of J 0 and a large change of Pe/ at the glacier front. On the other hand, glaciers with high speed change (warm-or cold-color curves, including accelerated and decelerated glaciers) seem to cluster together in a different area where J 0 is much more negative and Pe/ ≈ 0 m −1 . These curves generally show a vertical orientation indicating changing J 0 and constant low Pe/ along the glacier flowline.
To further illustrate the clustering trend, we arbitrarily select a speed change threshold of ±300 m yr −1 and classify the glaciers based on their speed change at 3 km. The threshold value is determined in order to give each classification roughly the same number of samples. Note that glaciers with significant slowdown or speedup are classified into the same group because a glacier vulnerable to basal lubrication would also be sensitive to recovering basal conditions from a lubrication event. Fifty-four glaciers have an absolute value of speed change ≥ 300 m yr −1 , and the other 50 glaciers are considered more stable with an absolute value of speed change < 300 m yr −1 . We plot J 0 against Pe/ using the val- ues from 3 km as well as the Gaussian kernel density estimates of each classification for both J 0 and Pe/ . The results using all glaciers (Fig. 6a) indicate that two classifications have a slightly different distribution for J 0 and Pe/ . The unstable glaciers (red group) have J 0 and Pe/ distributions peaked at ∼ −200 m yr −1 and ∼ 0.00003 m −1 respectively, whereas the peaks of stable glaciers (blue group) shift to higher values to ∼ −50 m yr −1 for J 0 and ∼ 0.00013 m −1 for Pe/ . Despite the peak shift being small compared with the distributions themselves, the difference is statistically significant found by the Kolmogorov-Smirnov statistic (p = 0.003 and 0.006 for J 0 and Pe/ , respectively; see the Fig. 6 Jupyter Book page for details). As the major cause of glacier speedups has been attributed to terminus retreat (King et al., 2020), we further group the glaciers by the distance of terminus retreat. Figure 6b plots only the glaciers with a retreat larger than 0.5 km, and the separation between two glacier groups in terms of J 0 and Pe/ becomes much less significant (p = 0.366 and 0.050 respectively). On the other hand, Fig. 6c shows the J 0 and Pe/ distributions for glaciers with a retreat ≤ 0.5 km. Only seven glaciers with little or no retreat have accelerated over 300 m yr −1 , which might not be enough to determine the significance of separation for Pe/ (p = 0.255). However, it is clear to see two groups divided by the value of J 0 in the plot (p = 2 × 10 −5 ).
Among all the glacier outlets, glacier 0207 (65.17 • N, 41.16 • W; Fig. 1a) seems to be unusual as it is the only one with J 0 over 100 m yr −1 (Fig. 5b). A positive J 0 requires a decreasing ice thickness from terminus to upstream (Eq. 17), which seems to indicate a steeper bed than the glacier surface. In our analysis, glacier 0207 is considered an example with large |J 0 |, despite having a different sign than the other glaciers. However, additional data from more glaciers would be required to fully characterize the glaciers with positive J 0 and their sensitivity to basal lubrication.
We adopt the same method from Fig. 5 to plot the results from eight marine-terminating glaciers in Austfonna (Fig. 7). As all glaciers have accelerated at the terminus for the past two decades, we adjust the color code so that blue represents low change and other colors represent high change. For basin-3 and -5, there are no valid measurements at 3 km, and we only mark the first valid measurements from 7.3 and 6.7 km, respectively, as single points on the plot. Similar to the GrIS, glaciers with higher speed change (basin-3, -5, -7, and -10) roughly occupy the lower left side of the panel where Pe/ and J 0 are small or more negative, and glaciers with lower speed change (basin-1, -4, -6, and -17) fall on a different corner, where Pe/ and J 0 are larger. However, |J 0 | of all eight glaciers have values between 0 and 10 m yr −1 , much lower than those from the GrIS (see Fig. 5, where |J 0 | ranges from 0 to over 1500 m yr −1 ). This is because all eight glaciers are slowly moving in 1996, resulting in a low |J 0 | according to Eq. (17). It is also interesting that basin-1 (Bråsvellbreen) and basin-5 have similar J 0 versus Pe/ , but basin-5 has a speed change roughly ten times greater than basin-1. Since basin-1 had a surge record back in 1936-1938(Schytt, 1969 and is currently in the quiescent stage, its low J 0 and Pe/ values might indicate a future instability when a surge is triggered internally or externally.

Separation of glacier groups on the J 0 versus Pe/ plot
For the GrIS, the J 0 versus Pe/ plots (Figs. 5-6) seem to capture the characteristics of glaciers vulnerable to basal lubrication. GrIS glaciers with more negative J 0 and Pe/ in 1996-1998 are more likely to speed up in the next 20 years. This tendency is not obvious for glaciers with a significant retreat (Fig. 6b), possibly because instead of changing basal conditions, it is the retreat dominating the flow dynamics of these glaciers. Nevertheless, the separation of peaks at different Pe/ values suggests that the diffusive strength of ice flow might still play a major role in the decades-long speed changes. For glaciers with a more stable terminus position, the clear separation of J 0 (Fig. 6c) supports basal lubrication as a primary cause of glacier speedups and likely highlights the importance of the initial response to changing basal conditions. Even if a glacier has a low P 0 and allows diffusiondominating dynamics, a small |J 0 | prohibits much elevation and speed change under a lubrication scenario (Eq. 9). The other factor that might affect the tendency of Pe/ and |J 0 | we see in Fig. 6 is whether the basal lubrication takes place within the study period. Although melt-induced speedups are common for GrIS glaciers (e.g., van de Wal et al., 2008;Bartholomew et al., 2010;Kehrl et al., 2017;Rathmann et al., 2017;Seddik et al., 2019), not all 104 glacier basins analyzed here have been studied well enough to identify when and where a glacier responds to changing basal conditions. A glacier with low Pe/ and high |J 0 | can remain stable if no basal lubrication has taken place in the past decades, mixing up the distributions in Fig. 6. Also, the Pe/ and |J 0 | patterns at the glacier front may not represent the entire glacier length. If moulins can form at a higher elevation than previously thought owing to the warming climate and widespread supraglacial lakes (Hoffman et al., 2018), it might be necessary to reassess the Pe/ and |J 0 | patterns within a wider range of flowline distance.
The three surge-type glaciers of Austfonna (basin-1, -3, and -17) do not cluster in Fig. 7. Compared with the other two glaciers, basin-3 has a higher flow speed in 1995-1996, resulting in a slightly more negative J 0 . It has had an unusually long surge evolution over two decades as well: the ice flow speed has gradually increased since the mid-1990s (Dowdeswell et al., 2008;McMillan et al., 2014) and reached a peak velocity of ∼ 6500 m yr −1 in 2013 (Dunse et al., 2015). The sustaining high flow speed has been attributed to meltwater routing through crevasses formed during the speedup (Dunse et al., 2015;Gong et al., 2018). As the additional supply of surface melt can alter the behavior of a surge-type glacier by reaching a steady state balancing mass and enthalpy conservation through thinner and faster-moving ice (Benn et al., 2019), basin-3 may have entered an ice stream-like regime with higher sensitivity to changing basal conditions. This might explain why basin-3 is away from the other two surge-type but currently quiescent glaciers on Fig. 7. Nevertheless, additional analysis and tests will be required before inferring a general vulnerability for surge-type glaciers to basal lubrication.

Characteristics of glaciers vulnerable to basal lubrication
Both Pe and J 0 depend on the ice thickness and flow speed (Eqs. 16 and 17), but with a different relationship. Assuming a monotonous decrease in glacier thickness toward the ter- minus (i.e., no overdeepening area), glaciers with thicker ice and a faster flow yield lower Pe and higher |J 0 | and thus are more susceptible to basal lubrication. However, the thickness change along the flowline (H ) has a competing contribution to Pe and J 0 : a greater change (i.e., a more negative H ) increases both Pe and |J 0 |. Glaciers with a very negative H may not likely be activated through intense diffusion (which is also suggested to be true for a terminus perturbation scenario in Felikson et al., 2021), but a collapse-like destabilization is still possible at the terminus or a localized region along the glacier owing to its high |J 0 |. For a glacier with an overdeepening zone, increased ice thickness again lowers Pe and raises |J 0 |, making the glacier more vulnerable to basal lubrication at the overdeepened area. Despite having an additional associating factor J 0 in the model, inferences made to the Péclet number in this study is similar to the previous models based on the perturbation theory. A low or negative Pe allows an ocean-induced ter-minus perturbation to propagate to a certain inland distance where Pe becomes larger (Felikson et al., 2017(Felikson et al., , 2021. For an outlet glacier in the GrIS, it is probably common to have terminal perturbation and changing basal conditions in effect at the same time (as indicated by Jakobshavn Isbrae for example; Joughin et al., 2008;Khazendar et al., 2019;Riel et al., 2021). In this case, Pe reflects a general vulnerability to elevation perturbations and is indistinguishable from the source forcing. On the other hand, J 0 as a new term in the lubrication-induced perturbation model seems to be exclusively related to the basal sensitivity as indicated by Fig. 6b and c. Nevertheless, J 0 might still be a key factor for a glacier only subject to the ocean-ice interaction. As warm subsurface water-induced thinning debuttresses the glacier front and increases the longitudinal stretching and glacier speed (Holland et al., 2008), new crevasses can provide additional routes for surface melt accessing and lubricating the bed (Gagliardini and Werder, 2018;Gong et al., 2018). The The average J 0 versus Pe/ at 3-5 km from the terminus. The value at 3 km is marked with a big dot. Each line uses the same color code from (a) and is labeled by the glacier number (see Fig. 2). Note that for basin-3 and -5, there is no valid measurement at 3 km, and only the first valid measurements (at 7.3 and 6.7 km respectively) are plotted.
investigation for surface strain rates indicates that these new crevasses can propagate to up to 1600 m high, corresponding to ∼ 50 km away from the terminus for GrIS outlet glaciers (Poinar et al., 2015). Thus, J 0 can be used to evaluate the latter mechanism's impact, specifically for subsequent ice flow acceleration or the feedback to additional thinning. This oceanic forcing-induced basal lubrication seems to be important for marine-terminating glaciers to switch to and maintain a fast flow over the years.
Our model does not consider the melt production from strain heating or geothermal heating at the bed. If included, induced glacier speedup owing to basal lubrication can generate energy to melt basal ice (Strozzi et al., 2017b), further increasing K 1 and leading to a higher glacier speed and thinning rate than what Eqs. (6) and (12) indicate.

Feedback from basal lubrication
One implication from the lubrication-induced perturbation model is that both Pe and J 0 are temporally changing variables. As glacier speed increases due to basal lubrication, Pe will be closer to zero, and |J 0 | will become larger, making the glacier more sensitive to any following change in basal conditions. A potential example to illustrate this feedback is the Vavilov ice cap, Severnaya Zemlya, Russia. The western marine outlet of this ice cap was moving at less than 1 km yr −1 with no apparent summer speedups just before a surge-like collapse took place (Willis et al., 2018). The collapse initiated when the terminus advanced into weak marine sediments, bringing the glacier speed to a maximum of ∼ 9 km yr −1 in summer 2015, and then the ice flow started to slow down but with a significant seasonal variation (Willis et al., 2018;Zheng et al., 2019). During the collapse, the Péclet number at the thinning center reduced by at least 40 %, suggesting a shift to an ice stream-like dynamic regime more susceptible to basal lubrication (Zheng et al., 2019). However, this transition is not necessarily irreversible as the bed would also be sensitive to a refreezing or efficient draining event, increasing Pe and decreasing |J 0 | back to the preperturbation level. Such a cycle is probably happening on a yearly basis as many GrIS glaciers have seasonal speedups closely bonded to summer melt (e.g., Palmer et al., 2011;Sundal et al., 2013;Rathmann et al., 2017). Still, for a multiyear destabilization like Vavilov, it is uncertain whether such a shutdown can completely revert the glacier dynamics because the ice thickness, another critical parameter controlling Pe and J 0 in our model, has changed a great deal during the collapse as well.
As noted in the case of Vavilov, dynamic thinning caused by the lubricated bed can also consequently change the dynamic regime and create another feedback loop. Thinning would decrease the ice thickness and increase the surface slope, potentially raising Pe and |J 0 |. Unlike the acceleration feedback from the previous paragraph, this feedback circle seems to be milder as a glacier would gradually switch to advection and prevent further inland thinning. However, dynamic thinning may also contribute to glacier retreat (Thomas and Bentley, 1978;Wood et al., 2021) and the subsequent debuttressing and speedup events. The net effect for dynamic thinning to glacier vulnerability to basal conditions remains ambiguous based on this view and suggests a future research topic as dynamic discharge in GrIS will likely continue to contribute significant ice loss in the near future (Mouginot et al., 2019;Choi et al., 2021).

Conclusions
Based on the new lubrication-induced 1D perturbation model, we show that a lubricated bed can initiate a thinning perturbation and destabilize the entire glacier if a particular combination of glacier thickness, thickness gradient, and flow speed is met. The model identifies two controlling physical quantities Pe/ (the Péclet number divided by the characteristic length) and J 0 (essentially the product of glacier speed and thickness gradient). We use observational data from 1996-1998 and derive these numbers for 104 and 8 marine-terminating glaciers in the GrIS and the Austfonna ice cap, Svalbard respectively. The results show that Pe/ and J 0 correlate with the flow speed change during 1996/1998 and 2018, especially for nonsurge glaciers and glaciers without significant terminus retreat, matching the model prediction. Glaciers with thick ice and a fast flow result in low Pe and negative J 0 , and reduced basal friction leads to initial speedup and thinning, which can propagate further inland via diffusion. For glaciers in the GrIS subject to ocean-ice interactions, this new model indicates multiple feedback cycles that make glaciers more sensitive to changing basal conditions. Finally, this study highlights glaciers classified as vulnerable to lubricated beds (low Pe and high |J 0 |) but with no significant speed change during the past two decades. Frequent monitoring is suggested for these glaciers because they might be more prone to future instabilities and affect the projected sea level rise.
Executable research compendium (ERC). The pejzero Zenodo archive (https://doi.org/10.5281/zenodo.5641953; Zheng, 2022) links to a Binder-ready Github repository. To launch the Binder server, click the "launch binder" button in the repository readme (https://github.com/whyjz/pejzero, last access: 23 February 2022) or the rocket button in the corresponding Jupyter Book pages (https://whyjz.github.io/pejzero/workflows/Fig3-4.html or other Fig*.html, last access: 23 February 2022). All the Jupyter Notebooks (in the "workflows" folder) are executable on the Binder server without installing additional dependencies. Alternatively, users can download the entire Zenodo archive, install dependencies specified in environment.yml and the postBuild file, and execute the code and Notebooks on a local machine.
Competing interests. The author has declared that there are no competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. Kudos to the advice from the Jupyter Meets the Earth (JMTE) team for making the study's scientific workflow as sharable and reproducible as possible. Comments from Matthew Pritchard, Michael Willis, and Facundo Sapienza greatly improved the quality of the paper. The author thanks Denis Felikson for addressing questions about the use of the GrIS flowline data. The author also acknowledges the National Snow and Ice Data Center QGreenland package (Moon et al., 2021) as an extremely convenient tool for inter-comparing the GrIS data sets and preparing maps presented in this paper. The final appreciation goes to two anonymous reviewers and their constructive and insightful suggestions for the entire study.
Financial support. This study is partially supported by the Jupyter Meets the Earth (JMTE) program, an NSF EarthCube funded project (grant nos. 1928406 and 1928374).
Review statement. This paper was edited by Ben Marzeion and reviewed by two anonymous referees.