1. Introduction
Observations of the Greenland Ice Sheet (GrIS) mass balance over the past four decades have revealed accelerating ice loss, contributing over 10 mm to global sea level rise (Mouginot and others, Reference Mouginot2019). This trend is projected to continue in the 21st century, with high-emission scenarios likely to induce a global sea level rise of 90 ± 50 mm beyond the present committed mass loss (Goelzer and others, Reference Goelzer2020). Mass loss is primarily driven by decreases in surface mass balance and increases in ice discharge, but precise partitioning between these two mechanisms is subject to large uncertainty in climate forcings (Fox-Kemper and others, Reference Fox-Kemper2023) and thus remains a target of active research. Lately, mass loss through discharge or glacier dynamics has been proposed as an important driver of mass loss in both historical observations and future projections (Mouginot and others, Reference Mouginot2019; Choi and others, Reference Choi, Morlighem, Rignot and Wood2021). Thus, understanding the mass loss caused by the ice dynamic response to climatic forcing is critical to predicting the future evolution of the GrIS.
Dynamic mass change tracked via ice thickness change is primarily driven by glacier motion, via ice deformation and basal sliding in response to stress disequilibrium, particularly due to interannual to decadal-scale changes in ice frontal geometry from calving events (Nick and others, Reference Nick, Vieli, Howat and Joughin2009; Christian and others, Reference Christian2020). Over the past two decades, observations have revealed widespread retreat of outlet glaciers (Moon and others, Reference Moon, Gardner, Csatho, Parmuzin and Fahnestock2020; Goliber and others, Reference Goliber2022) primarily caused by the intrusion of comparatively warming North Atlantic water into fjords and submarine melting at the termini (Slater and others, Reference Slater2020; Wood and others, Reference Wood2021). These retreats trigger ice flow accelerations and along-flow divergence, leading to thinning caused by ice dynamics that propagates upstream, in some cases penetrating dozens of kilometers inland (Pritchard and others, Reference Pritchard, Arthern, Vaughan and Edwards2009; Wang and others, Reference Wang, Li and Zwally2012; Khan and others, Reference Khan2013; Felikson and others, Reference Felikson, Catania, Bartholomaus, Morlighem and Noël2021).
Despite its widespread occurrence, the thinning caused by ice dynamics (hereafter referred to as dynamic thickness change) exhibits complex temporal and spatial patterns even among neighboring glaciers subject to similar oceanic forcing (McFadden and others, Reference McFadden, Howat, Joughin, Smith and Ahn2011; Csatho and others, Reference Csatho2014; Khan and others, Reference Khan2013, Reference Khan2014). This implies the influence of local factors, such as fjord geometries and boundary conditions. Recent studies have highlighted the role of fjord width and depth on glacier stability (Bassis and Jacobs, Reference Bassis and Jacobs2013; Enderlin and others, Reference Enderlin, Howat and Vieli2013; Carr and others, Reference Carr, Stokes and Vieli2014; Haseloff and Sergienko, Reference Haseloff and Sergienko2018; Steiger and others, Reference Steiger, Nisancioglu, Åkesson, De Fleurian and Nick2018; Frank and others, Reference Frank, Åkesson, De Fleurian, Morlighem and Nisancioglu2022), which collectively govern the force balance structure and thus the terminus response to perturbations (Carnahan and others, Reference Carnahan, Catania and Bartholomaus2022). Although the terminus exerts critical control over inland flow dynamics, other hydro-mechanical processes are also important, including basal hydrologic processes that regulate ice flow dynamics. Basal lubrication caused by surface meltwater drainage has been extensively documented across the GrIS, resulting in seasonal acceleration and deceleration of ice flow (van de Wal and others, Reference van de Wal2008; Bartholomew and others, Reference Bartholomew2010; Chandler and others, Reference Chandler2013; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). While most studies focus on flow velocity, dynamic thickness change caused by basal lubrication has also been observed (Bevan and others, Reference Bevan, Luckman, Khan and Murray2015), and yet the records are comparatively sparse. Moreover, how the dynamic thickness of glaciers at various dynamical states responds to these basal perturbations remains uncharacterized (Zheng, Reference Zheng2022). Aside from observational studies, numerical simulations generally represent basal processes via parameterization known as sliding laws. However, it remains unclear how individual terms in the sliding laws, such as the effective pressure dependence, affect the simulated dynamic thickness change and its rate of change in different geometric configurations (Joughin and others, Reference Joughin, Smith and Schoof2019; Barnes and Gudmundsson, Reference Barnes and Gudmundsson2022; Felikson and others, Reference Felikson, Nowicki, Nias, Morlighem and Seroussi2022). This limitation hinders our progress in better initializing ice-sheet models (Aschwanden and others, Reference Aschwanden, Aøalgeirsdóttir and Khroulev2013) and therefore short-term projections of future ice loss (Goelzer and others, Reference Goelzer2018).
In this study, we examine the interplay between basal processes and glacier geometries in controlling patterns of dynamic thickness change. Specifically, we investigate two distinct types of basal perturbations that produce differing spatio-temporal impacts on ice thickness change. The first type involves variations in basal drag due to changes in ice overburden pressure. Ice overburden pressure is directly determined by the ice thickness, yet its impact on dynamic elevation change is rarely explored systematically (Habermann and others, Reference Habermann, Truffer and Maxwell2013; Joughin and others, Reference Joughin, Smith and Schoof2019). Nonetheless, it has been identified as a critical component in the tidewater glacier cycle, where frontal retreat leads to ice thinning, reduced effective pressure and basal drag, flow acceleration, and further thinning of a glacier (Benn and others, Reference Benn, Hulton and Mottram2007; Pfeffer, Reference Pfeffer2007). The second type is a localized perturbation of basal drag at the inland portion of the glacier, most commonly due to a change in effective pressure through a change in basal pore pressure. Observational studies have shown occurrences of localized dynamic elevation change far from the terminus, possibly caused by supraglacial lake drainages or changes in basal hydrologic system (Bevan and others, Reference Bevan, Luckman, Khan and Murray2015; Stevens and others, Reference Stevens2022). At fast-flowing outlet glaciers where basal sliding dominates over vertical deformation, the localized basal variability can have non-local effects on flow velocity and dynamic elevation change where theoretical consideration may fall short (Gudmundsson, Reference Gudmundsson2003; Sergienko and Hulbe, Reference Sergienko and Hulbe2011; Sergienko, Reference Sergienko2013), and therefore a numerical-model-based systematic characterization of dynamic thickness change throughout the glacier domain is much needed.
Here we investigate these two processes using numerical experiments on various idealized Greenland-like outlet glaciers. Using idealized glacier geometries that are broadly representative of multitudes of real-world glaciers allows a generalizable study of how different forcings affect the evolution of ice-surface elevation. It minimizes the tailoring of simulations to highly specific glacier characteristics, e.g. fjord size and shape, bed topography or basal drag. Recent studies have used idealized glacier simulation to examine glacier mass loss bias from terminus forcing temporal frequency (Felikson and others, Reference Felikson, Nowicki, Nias, Morlighem and Seroussi2022), terminus response to topographic features (Frank and others, Reference Frank, Åkesson, De Fleurian, Morlighem and Nisancioglu2022) and the impact of meltwater inputs on downstream ice velocity (Poinar and others, Reference Poinar, Dow and Andrews2019). In this study, we similarly construct a suite of idealized synthetic glaciers with variations in glacier geometric parameters and basal boundary conditions, referring to each constructed glacier as a ‘synthetic glacier testbed’ or simply ‘testbed’. For each testbed, we test and characterize the impact of changes in ice overburden pressure and localized basal lubrication on dynamic thickness change.
2. Methodology
2.1 Model setup
We utilized the Ice-sheet and Sea-level System Model (ISSM) to conduct the numerical experiments. ISSM is a state-of-the-art finite element package that can simulate glacier and ice-sheet scale flow dynamics (Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012) and we refer readers to Larour and others (Reference Larour, Seroussi, Morlighem and Rignot2012) for details of the modeling package and governing equations. To simulate the outlet glacier flow, we employed the 2D Shallow Shelf Approximation (MacAyeal, Reference MacAyeal1989) of ice flow physics on both grounded and floating ice. A uniform triangular meshing with a spatial resolution of 200 m was adopted throughout the model domain (12 km × 60 km). To account for the evolution of the grounding line position, we implemented a sub-element migration scheme where the sliding law coefficient at partially grounded elements scaled with the fraction of the grounded area (Gladstone and others, Reference Gladstone, Lee, Vieli and Payne2010). While the grounding line migrates dynamically according to the hydrostatic criterion, we prescribed the calving front migration enabled by the level set method in ISSM (Bondzio and others, Reference Bondzio2016).
We used a time-independent surface mass balance (SMB) across all the experiments and testbeds. This is because the impact of SMB variability on ice dynamic thickness occurs at timescales longer than our decadal-scale model runs (Christian and others, Reference Christian2020), precluding an ability to test SMB effects. We used Glen's flow law with n = 3 for all simulations. We assumed a uniform ice temperature of $-3^\circ$
C. Below we will provide a summary of forcings, model geometry and experimental designs. For mathematical details, please refer to the Appendix B.2. Variables defined throughout the main text and appendix can be found in Table 1.
Table 1. Parameters in synthetic testbeds and experiments
‘Variable parameters’ refers to values of a variable that differs across synthetic testbeds. Readers can refer to Table A1 in the supplementary material for the parameters grouped by each testbed.
2.2 Synthetic glacier testbeds
We adapted and modified the idealized Greenland outlet glacier geometry from Felikson and others (Reference Felikson, Nowicki, Nias, Morlighem and Seroussi2022), which itself was based on the Marine Ice Sheet Model Intercomparison Project geometry (Asay-Davis and others, Reference Asay-Davis2016, MISMIP). The calving front was initially located at 56.5 km from the influx boundary. We prescribed an across-flow bed topography similar to Felikson and others (Reference Felikson, Nowicki, Nias, Morlighem and Seroussi2022), but the differences are that in our model, the bed was flat in the along-flow direction and the width of the trough w c(x) narrowed quadratically along flow in the upper reaches of the model domain (Eqn (B.7)). Nonetheless, as an extended inquiry to findings we will discuss later, we also briefly investigated the influence of bed roughness on dynamic thickness change patterns (Fig. 2d), where we performed additional simulations using a bed with fractal roughness.
For model initialization, we adopted a Weertman sliding law (Weertman, Reference Weertman1957) describing sliding over a hard bed:
Here τb is basal shear stress, m is a prescribed constant assuming certain sliding mechanics, C w is the prescribed Weertman law coefficient field defined in Eqn (B.8), and vb is the sliding velocity. We used the sliding law and assumed m = 1 for three primary reasons: first, its simplicity makes it the most commonly used sliding law and exponent in ice-sheet modeling, and hence our findings will be relevant for modelers; second, the Weertman sliding law does not incorporate dependence on effective pressure and so when run against simulations with Budd sliding law, it can help isolate the impact of overburden pressure on dynamic thinning; third, the Weertman sliding law is valid at the high effective pressure limit, as both the Schoof and Tsai sliding law formulations (Schoof, Reference Schoof2005; Tsai and others, Reference Tsai, Stewart and Thompson2015) asymptotically approach the Weertman formulation at higher effective pressure.
To construct a suite of testbeds, we varied the width W of the fjord at the narrower end, the grounding line depth B gl (zero at sea level) and the sliding law coefficient C w, producing in total 18 testbeds as illustrated in Figure 1. To the first order, the prescribed sliding law coefficient magnitudes control mean basal drag levels near the termini (Table A2).
Figure 1. Synthetic testbeds and examples. The top panel shows three variables of interest. (1) Sliding law coefficient. (2) Grounding line depth and frontal geometry. (3) Fjord width. With the flow domain length fixed, the grounding line depth is adjusted via changing bedrock slope β, where testbeds with deep grounding line and floating termini (‘Deep’) have greater bed slope (β + = −0.012), and the ones with shallow grounding lines and fully grounded termini (‘Shallow’) have lesser bed slope (β − = −0.005). Four examples of testbeds are shown in the bottom panel, with the steady-state ice speed colored and superimposed on the surface.
We allowed the testbed glaciers, over a maximum of 500 simulation years, to reach their steady-state defined as dh/dt < 0.01m a−1 everywhere in the flow domain. At steady state, testbed glaciers with shallower grounding line depths were grounded across the whole domain, whereas testbeds with deeper grounding line depth developed floating sections up to 12 km long (Fig. 1 and Table A1). This is broadly consistent with Northern Greenland outlet glaciers (Hill and others, Reference Hill, Rachel Carr, Stokes and Hilmar Gudmundsson2018). For simplicity, we refer to glaciers with deep grounding lines and floating termini as ‘deep testbeds’, and their fully grounded shallow counterparts with shallow grounding lines as ‘shallow testbeds’. The 18 testbeds differ significantly in their average and maximum flow velocity near the terminus (Fig. 1 and Table A2).
2.3 Experiment design
For each testbed glacier, one control run and two perturbation experiments were conducted, and all simulations started at the same initial state, the steady state after model relaxation.
2.3.1 Control run
Previously studies have shown strong correlation between the evolution of terminus position and flow dynamics in certain glaciers (Nick and others, Reference Nick, Vieli, Howat and Joughin2009; Cheng and others, Reference Cheng, Morlighem, Mouginot and Cheng2022), but simulating terminus motion is known to be a challenging task due to a variety of under-constrained processes involved (Benn and others, Reference Benn, Hulton and Mottram2007; Bassis and Jacobs, Reference Bassis and Jacobs2013; Robel, Reference Robel2017; Slater and others, Reference Slater, Nienow, Goldberg, Cowton and Sole2017; Choi and others, Reference Choi, Morlighem, Wood and Bondzio2018; Slater and others, Reference Slater2019; An and others, Reference An2021). Therefore in this study, we did not aim to reproduce a sequence of terminus positions comparable to observational records. Instead, we forced the terminus in all testbeds to retreat identically throughout all the experiments.
After a testbed glacier is initialized to its steady state, we forced the calving front to retreat at a time-variable rate described by a triangular function that spans 16 years (grey box in Fig. 2A). The calving front experiences an accelerating retreat for eight years, decelerates for 8 years, and stabilizes. We designed this pattern to represent a smoothed-step decadal retreat of a calving front, broadly similar to the observed terminus retreats of many outlet glaciers around GrIS in the past 20 years, where the early 2000s marked the onset of widespread retreat, followed by a period of relative stability in the late 2000s through early 2010s (Khazendar and others, Reference Khazendar2019). Details regarding the control run can be found in Appendix A.3.3.
Figure 2. Testbeds and experiment designs. (a) Control run. The terminus is forced to retreat at a time-variable rate according to the triangular function (orange). (b) Overburden pressure experiment. The basal drag τ b decreases as a result of diffusive thinning from the retreating terminus. (c) Localized basal perturbation experiment. In addition to changes in overburden pressure due to thinning, a Gaussian-shaped region of lower sliding law coefficient is applied transiently 24.5 km upstream of the terminus. The magnitudes ϕ of the two types of temporal variability (‘Transient Pulse’ and ‘Diffused Pulse’) are shown in brown. The perturbation locally induces upstream thinning (blue) and downstream thickening (red). (d) Experiment with a rough bed. Zero in the elevation offset means no change concerning the original constant bed slope. Both the overburden pressure and localized basal perturbation experiment are repeated on a testbed glacier with a rough bed.
2.3.2 Overburden pressure experiment
The basal drag of a glacier depends on the contact area between the ice and the bedrock. It is regulated by a competition between the opening of cavities from sliding over bumps or melting and creep closure of ice (Cuffey and Paterson, Reference Cuffey and Paterson2010; Schoof, Reference Schoof2010), which manifests as varying effective pressure. To account for the dependence on the pressure, a sliding law alternative to Weertman's law, commonly known as Budd's law (Budd and others, Reference Budd, Keage and Blundy1979), is used:
where C b is the coefficient for the Budd sliding law and N is the effective pressure defined as the difference between ice overburden pressure ρ ig H and pore water pressure p w, i.e. N = ρ ig H − p w; m and q are sliding law exponents where we assume m = q = 1. In Budd's formulation, initial thinning near the glacier terminus will reduce the ice overburden pressure and hence the effective pressure N, reducing the basal drag and causing acceleration. The acceleration can lead to flux divergence that further reduces the ice overburden pressure, potentially precipitating positive feedback.
We investigated the impact of the varying overburden pressure on dynamic thinning and hence we refer to this experiment as the ‘overburden pressure experiment’. This is effectively the same simulation as the control run (section 2.3.1) but with Budd sliding law. To mimic the Budd sliding law with models initialized with the Weertman sliding law, we forced the terminus to retreat in the same fashion as in the control run, meanwhile adjusting the basal drag coefficient C w to compensate for changes in ice overburden pressure (for derivation details, see Appendix B.3.2):
where ρ i is the ice density, H(0) represents ice thickness values at the start of the experiment or the end of the model relaxation, and $\hat {C_{\rm b}}$
is the equivalent Weertman sliding law coefficient in Budd's formulation at steady state, i.e. $\hat {C_{\rm b}} = C_{{\rm w}0}/( \rho _{{\rm i}} g H( 0) ) ^{1/2m}$
. This amounts to representing Eqn (2) by modifying Eqn (1). As discussed above, in all experiments outlined in Figure 2 we assumed m = q = 1, but we also explored more plastic bed rheology (i.e. m = 5, Fig. A.3) and compared results to the linear viscous case in the discussion.
2.3.3 Localized basal perturbation experiment
In addition to the overburden pressure change discussed above, we considered the impact due to local drainage of meltwater to the bed. It was represented ideally by a localized basal drag reduction as a Gaussian-shaped patch of lower sliding law coefficient, centered 24.5 km behind the initial calving front. We used this location because it was immediately upstream of the most retreated grounding line in our control runs so that the localized perturbation remained engaged throughout the simulations.
We considered two types of temporal variability, Transient Pulse and Diffused Pulse, to represent the temporal variation of perturbation magnitude (Fig. 2c). Transient Pulse is a short-lived perturbation lasting for 0.1 years, which we designed to loosely represent the response of an efficient subglacial drainage system to supraglacial lake drainage or a rain event. The Diffused Pulse spanned 2 years with a lower peak value and integrated to the same total slipperiness perturbation as the Transient Pulse (Eqn (B.19)). We chose 2 years as a bounding case to provide a substantial contrast with the Transient Pulse signal. It was not designed based on observations of any specific glaciers, although we would discuss certain observations and model inferences that suggest a similarly prolonged period of reduced basal drag. There are a total of eight perturbation cycles and hence 16 years of perturbation. Details regarding the localized basal perturbation experiment can be found in Appendix B.3.3.
2.4 Bed constructed with fractal roughness
Glacier beds around GrIS are wavy at a range of length scales. This waviness is well characterized by fractal roughness (Jordan and others, Reference Jordan2017), meaning the asperity height at various wavelengths can be described by a Hurst exponent in a power law. To investigate the impact of bed roughness on dynamic thickness change, we generated a randomly rough surface superimposed onto a sloped flat bed (Mona Mahboob Kanafi, 2023), with a Hurst exponent of 0.8 and a root-mean-square roughness of 70 m (Fig. 2d). Similar values were used by Christian and others (Reference Christian, Robel and Catania2022) for the GrIS and are within the range of roughness estimates from radar observation (Jordan and others, Reference Jordan2017). The specified mean roughness stipulates the average height of bed bumps; in our glacier domain, the bumps that the grounding line retreats over are less than 100 m in height. The results are discussed in section 3.3.
2.5 Estimating frontal resistive stress loss
The diverse geometries and mean basal drag levels considered produce various stress balance regimes and changes in stress balance in response to the calving front and grounding line retreat. To quantitatively assess the changes, we follow the calculation outlined in van der Veen and Whillans (Reference van der Veen and Whillans1989) and Carnahan and others (Reference Carnahan, Catania and Bartholomaus2022) to estimate the stress components. The stress balance states that the gravitational driving stress of a glacier is approximately in balance with the sum of the basal shear stress, and the longitudinal, and lateral resistive stress gradients.
We define frontal resistive stress as the sum of the lateral, longitudinal and basal resistive stress from the current grounding line to the ice front. Hence, we define the frontal loss of resistive stress as the total change in the resistive stress throughout the model runs. Mathematical details are presented in Appendix B.4. The results are presented in section 3.4 and discussed in section 4.2.
3. Results
3.1 Overburden pressure experiment
As the terminus retreats, in all testbeds, dynamic thinning originated near the terminus and diffused upstream, and the largest degree of thinning was found behind the grounding line. If we isolate the thinning induced by overburden pressure feedback, for fully grounded testbed glaciers with shallower grounding lines, the sliding law correction for ice overburden pressure added a maximum of 97 m over 16 years, or 6 m a−1 (Fig. 3) and all grounding lines remained grounded throughout (e.g. Fig. 3a). Model testbeds with deep grounding lines (Figs 3b–d) showed a substantially larger degree of thinning accompanied by continued grounding line retreat. The deep narrow testbed with high basal drag (Fig. 3d) showed the most thinning, 250 m over the 16-year model run or an average thinning rate of 16 m a−1.
Figure 3. Dynamic thickness change due to changes in ice overburden pressure. All 18 testbeds are represented as colored circles in a 3 × 6 grid separated by the grounding line depths. The circular marker represents both the maximum dh/dt observed along the center flow line (marker size) and the attenuation distance of diffusive thinning (color). A shorter attenuation distance suggests stronger thinning attenuation. All values can be found in Tables A4 and A5. Four selected testbed glaciers are shown in greater detail. The lateral profiles show the evolution of ice thickness from the overburden pressure experiment, whereas the line plot at the top of each subplot shows the thickness change isolated (ΔH) from the effect of ice overburden pressure (i.e. ΔH = H(overburden pressure exp.) − H(control) as in Fig. 2). Black lines show the lateral profiles at the new steady states.
The colored circles in Fig. 3 illustrate how the maximum dh/dt and attenuation distance varies across fjord widths, mean basal drag levels and frontal geometries. Attenuation distance is defined as the distance from the ice front where the cumulative thickness change has dropped to 36.8% (e-folding length 1/e) of the total thickness change. At all testbed glaciers, attenuation distance was primarily controlled by the mean basal drag: high basal drag corresponded to larger thickness change attenuation, and vice versa. Maximum thinning rate, however, exhibited a more nuanced relationship with geometry and basal condition. At testbed glaciers with high mean basal drag (e.g. mean basal drag near the terminus >60 kPa in Table A2), the effect of fjord width was more pronounced, with narrow testbed experiencing greater maximum thinning rate up to 16 m a−1 despite less grounding line retreat, and wide testbed experiencing <10 m a−1 thinning. Conversely, at testbeds with lower mean basal drag (e.g. mean basal drag <30 kPa in Table A2), differences in fjord width did not result in variances in max thinning rate (10.4 − 10.5 m a−1).
3.2 Localized basal perturbation experiment
We present the results of the localized basal perturbation experiment as their difference in dynamic thickness change from the ice overburden pressure experiment. Since the localized basal perturbation experiment accounts for overburden pressure change by design (Fig. 2c), we are merely isolating the thinning caused by the localized basal perturbation alone. Immediately after it is introduced, the perturbation caused transient thickening on the downstream glacier and transient thinning on the upstream portion, regardless of the magnitude or duration of the forcing (Figs 4, 5). This dipole pattern is consistent with the results of previous theoretical studies (Gudmundsson, Reference Gudmundsson2003; Sergienko and Hulbe, Reference Sergienko and Hulbe2011; Sergienko, Reference Sergienko2013).
Figure 4. Spatio-temporal patterns of dynamic thickness change at deep and narrow testbed glaciers in response to the two types of localized basal perturbation pulses. The space-time plots (essentially a Hovmöller diagram) are created by plotting the thickness change (colors) along the center flow line (y-axis) over time (x-axis). All the results presented here account for the changes in ice overburden pressure on the basal drag. The relative grounding line position on the top plots (labeled ‘Δ GL(m)’) is the difference in grounding line position between the control run and the experiment run; the solid line ‘Grounding line’ only shows the grounding line from the experiment run for visual simplicity. The y-axis label ‘Distance to front’ refers to the ice front location at t = 0. The thin vertical dotted line marks the end of frontal retreat and local perturbation. The cyan dotted line marks the perturbation location. The two types of pulse forcings are shown at the top of each panel. The amplitudes of the pulses are illustrative and thus not to scale. (a) A testbed glacier with low mean basal drag (τ b) forced with Transient Pulse. (b) A testbed glacier with high τ b forced with Transient Pulse. (c) A testbed glacier with low τ b forced with Diffused Pulse. (d) A testbed glacier with high τ b forced with the Diffused Pulse.
Figure 5. Spatio-temporal patterns of dynamic thickness change at deep and wide testbed glaciers in response to the two types of localized basal perturbation pulses. Graphic features are identical to Figure 4. (a) A testbed glacier with low mean basal drag (τ b) forced with Transient Pulse. (b) A testbed glacier with high τ b forced with Transient Pulse. (c) A testbed glacier with low τ b forced with Diffused Pulse. (d) A testbed glacier with high τ b forced with Diffused Pulse.
Over multiple perturbation cycles, the amplitude of the transient response increased as ice flow sped up and the glacier thinned. The maximum observed thinning or thickening did not exceed 20 m concerning the state before the perturbation was engaged. Within each perturbation cycle, thickening and thinning at the site relaxed more quickly in testbed glaciers with lower mean basal drag and, consequently, higher flow speeds. The relaxation is particularly visible when the model is perturbed by the Transient Pulse (e.g. Fig. 4). Between testbeds, the dipole amplitudes showed amplitude differences of less than 12 m near the perturbation site (Table A3). At both deep and shallow testbed glaciers, we observed generally similar patterns in the dipole amplitude and its temporal variation. Therefore, for simplicity of presentation, we show the results of the localized basal perturbation experiment for only the deep testbeds, and all the ensuing qualitative discussions apply to shallow testbed glaciers as well unless indicated otherwise. Results from selected shallow testbeds can be found in the Appendix (Figs A5, A6).
Over time, trends in dynamic thickness change emerged both near and far from the perturbation site. Widespread thinning occurred 5–15 km upstream of the perturbation, while downstream, variable patterns of thickening and thinning occurred at different testbeds. At testbeds with lower mean basal drag (A and C in both Figs 4, 5), thinning propagated farther outward from the perturbation site, whereas at testbeds with higher mean basal drag (B and D in both Figs 4, 5), these attenuated closer. The total degree of far-field thinning over the long term depends on the type of perturbation pulse used, with the Diffused Pulse resulting in generally twice as much thinning or thickening as the Transient Pulse.
More substantial differences in spatio-temporal patterns can be observed in the downstream trunk, particularly after several perturbation cycles. We present a few examples here. For the narrow testbed with a low mean basal drag level (Fig. 4a), the basal perturbation incited initial thickening in the downstream trunk that was, within ~10 years, overridden by the diffusive thinning from the trunk upstream. Similarly, in the first 5 years of the experiment, the grounding line advanced slightly before retreating by about 40 m, relative to the control run. A qualitatively similar pattern can be observed in the narrow testbed with a high mean basal drag level (Fig. 4b), but in this case, net thinning (relative to the control run) emerged near the grounding line after the third perturbation cycle. This thinning reached ~3 m and diffused upstream; unlike in the low-basal-drag testbed, the thinning continued after the perturbations ceased, spreading throughout the domain.
When forced with the Diffused Pulse, these two testbeds exhibited similar spatial and temporal patterns (Figs 4c,d). However, there was more thickening and less thinning and the grounding lines advanced farther.
Figure 5 shows results on wide testbeds. Here, the spatiotemporal patterns were generally similar to those observed in narrow testbeds, except that the upstream and downstream thickness changes were more polarized, with the upstream dominantly thinning and the downstream dominantly thickening throughout the perturbation cycles (with the minor exception of the low-basal-drag testbed in Fig. 5a). An extreme example is the testbed glacier with a high mean basal drag level forced with the Diffused Pulse (Fig. 5d), where the downstream thickening was not overtaken by upstream thinning years after the perturbation had stopped (in contrast to Fig. 5c, e.g.). It is noteworthy that the grounding lines in testbed glaciers with a low basal drag level (Figs 5a,c) moved much more rapidly and extensively, with advance and retreat ranging from ~200 to 400 m – an order of magnitude greater than in high-basal-drag testbeds. In all experiments, regardless of patterns, the maximum thickness change caused by the localized basal perturbation did not exceed 12 m over the 26 years of the simulation run (see Table A3).
3.3 Influence of bed roughness
Due to the asymmetry of grounding line flux dynamics at prograde and retrograde sections of the bed (Schoof, Reference Schoof2007), we hypothesize that an idealistic smooth terminus retreat can translate into episodes of fast and slow grounding line movement as it retreats over the bed asperities, potentially giving rise to a different timescale of variability in dynamic thickness change time series observed across GrIS (Csatho and others, Reference Csatho2014). We explored this possibility with two additional simulations of the overburden pressure experiment and localized basal perturbation experiment, using a testbed with high mean basal drag in a narrow fjord with fractal roughness throughout the bed (Fig. 2d). The resulting grounding line movement is characterized by step-wise retreats, corresponding to faster and slower periods of thickness change (Figs 6a, 7). We also observe that grounding line retreat stabilizes on the lee side of the bed bumps (Figs 6a,b) that stops further thinning after calving front perturbation ceases, in contrast to the original flat bed simulation (Fig. 7b).
Figure 6. Dynamic thickness change over an undulating bed. (a) Ice thickness, grounding line and calving front change over time. Smooth multi-year front retreat causes step changes in the grounding line, temporally matching the periods of faster and slower dynamic thinning. Time series are extracted at the location marked as a red circle in B and C. Colored dots over the grounding line are the same as those dots in panel B but are plotted here to better visualize the retreat distance. (b) Lateral profiles of basal topography and ice surface elevation along the glacier centerline (the horizontal dotted line in panel c). (c) Dynamic thickness change rate (contours) at the last time step (year 16) superimposed onto the basal topography (colors) near the ice front and grounding line. Ice at the central topographic low becomes ungrounded and experiences a low thinning rate; ice at the topographic high nearby undergoes a much higher thinning rate.
Figure 7. Comparing dynamic thickness change over a flat and an undulating bed forced by localized basal perturbation. The dotted line box outlines the time and space where thinning diverges after perturbation stops. (a) Isolated thickness change due to the localized basal perturbation at a rough bed. (b) Same but at a flat bed (Fig. 4b repeated).
For the rough bed, dynamic thickness change rates also exhibit spatial heterogeneity. Here we observe the topographic low behind grounding line attains flotation near the end of simulation (Fig. 6c) and the thinning rate dwindles, at 0−4 m a−1, while its neighboring topographic high experiences 8−12 m a−1 of thinning.
3.4 Stress loss and correlation with thinning
We examine the correlation between the frontal resistive stress loss, the magnitude of dynamic thinning and the grounding line retreat distance of deep testbed glaciers in the overburden pressure experiments, as they exhibit the most dynamic changes. Figure 8a shows that when integrated over the central flowline, the total glacier thinning magnitude positively correlates with frontal resistive stress loss (r 2 = 0.97). Testbed glaciers with lower basal drag experience larger grounding line retreats and total thinning, while no clear clustering pattern exists for fjord width.
Figure 8. Relationships between total thinning, maximum thinning, grounding line retreat and frontal resistive stress loss at the end of perturbations (simulation year = 16) for deep testbeds in the overburden pressure experiment. Each marker represents a distinct testbed. R 2 values report the goodness of fit of selected data by a linear regression model. (a) Relationship between total thinning versus grounding line retreat distance (triangles), and total thinning versus frontal resistive stress loss (circles). (b) Relationship between the spatial maximum thinning and grounding line retreat distance (triangles) and frontal resistive stress loss (circles). (c) Detail of (b) with only the three testbeds with narrow fjords. The dashed lines with arrows point to testbeds of increasing mean basal drag. Sizes of markers are enlarged concerning (b) for better presentation.
When examining maximum thinning rather than total thinning, frontal resistive stress loss is a stronger predictor. Figure 8b shows that the R 2 value for the positive correlation between the maximum thinning and frontal resistive stress is 0.99, significantly larger than that of the correlation concerning grounding line retreat. Testbed glaciers with narrow fjord widths generally experience higher frontal resistive stress loss, while no clear clustering pattern exists for basal drag.
Figure 8c shows that, specifically at testbeds with narrow fjords, lower mean basal drag results in greater grounding line retreat yet lower spatial maxima in thinning. In fact, at narrow fjords, grounding line retreat anti-correlates with the spatial maxima in thinning; this is not the case in moderate-width and wide testbeds, as shown in the trends across sets of the larger-sized triangles in Figure 8b, as these testbeds do not exhibit either a monotonically positive or negative trend. Figure 8c, therefore, highlights the strong geometric impact on the maximum dynamic thinning.
4. Discussion
4.1 Grounding line position correlates with dynamic thinning
Our experiments show that the grounding line positions correlate better with dynamic thinning rates than the ice front position does (Fig. A2), a commonly used observable in both modeling and observational studies (Bondzio and others, Reference Bondzio2017; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). We ran all testbed simulations with the same ice front position forcing but obtained a wide range of thinning degrees and variability (Figs 3, 4, 5), suggesting the limited predictive power of ice front position alone. Most thinning is observed behind the grounding line, as model results for Pine Island Glacier also showed (Joughin and others, Reference Joughin, Smith and Schoof2019) despite the significant difference in Antarctic glacier geometry from the Greenlandic counterpart. Similar dynamics were observed at Kangerlussuag Glacier (Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017) where the termini stabilized but the glacier continued to thin dynamically as the grounding line retreated, even as the glacier rested on a prograde bed. At Sermeq Kujalleq (Jakobshavn), migration of the unknown grounding zone and ungrounding was argued to partly explain the abnormally high thinning rates (Hurkmans and others, Reference Hurkmans2012).
The simulated movement of the grounding line is highly dependent on the choice of sliding law (Brondex and others, Reference Brondex, Gagliardini, Gillet-Chaulet and Durand2017). Therefore, knowledge of the specific bed rheology and sliding mechanics is crucial to accurately reproduce grounding line movements from observations. Our experiments with the Weertman and Budd sliding laws are two bounding cases for the magnitude of grounding line retreat (Brondex and others, Reference Brondex, Gagliardini, Gillet-Chaulet and Durand2017). In that study, greater retreat distance of the grounding line was found to correlate with greater thinning; our results reproduce this finding for multiple glacier geometries and mean basal drag levels.
The crucial role of grounding lines in dynamic thickness change is also highlighted in our localized basal perturbation experiments. We found that, across testbed glaciers of varying widths and sliding laws, downstream elevation change patterns strongly correlate with relative grounding line movement. One striking example is the pronounced thinning near the grounding line as the grounding line retreats relative to its initial position (e.g. Fig. 4b). This thinning nearly overtakes the local thickening signal immediately downstream of the perturbation near the end of the experiment. Similarly, continued relative grounding line advance causes downstream thickening (e.g. Fig. 5d). Despite repeated forcing, the diversity of grounding line movements and dynamic thickness change patterns suggests that one must consider both grounding line movement and glacier geometry when interpreting thickness change records, with all else assumed equal.
Despite the critical role of grounding line movement, its sensitivity to basal topographic undulation (Figs 6, 7, and Enderlin and others (Reference Enderlin2016)) implies that more dramatic or subdued dynamic thinning near the grounding line is possible depending on the bed roughness (Thomas and others, Reference Thomas, Frederick, Krabill, Manizade and Martin2009). Dynamic thinning can also happen when the grounding line is fairly stable due to bed asperities while the ice front retreats (Fig. 6a, year 10–12, e.g.) as the glacier geometry continues to adapt to the new ice front position. At a minimum, we stress the role of the grounding line either in initiating or expressing dynamic thickness change, even if the perturbation is localized tens of kilometers upstream of the terminus.
4.2 Resistive stress change influences the spatial variation of dynamic thinning
Our results show that while the grounding line position is strongly correlated with centerline-integrated total thinning and average thinning rate (Fig. 8a), it gives far less insight into the spatial pattern of thinning, here represented by the spatial maximum in thinning (Fig. 8b). Resistive stress change is the more important variable for spatial variations in thinning.
Despite the same frontal retreat forcing, the force balance response differs across different frontal and grounding line retreat outcomes. Specifically, calving of fully grounded testbed glaciers removes basal resistive stress, whereas at a floating terminus, the loss of the longitudinal stress gradient associated with calving is typically orders of magnitude less. Therefore, for the same prescribed terminus retreat, fully grounded testbed glaciers should experience more thinning. This explains the pronounced difference in the maximum thinning rate at glaciers with high basal shear stress but different fjord width (Fig. 3), as the differences in the loss of resistive stress are significant, from 500 to 1300 MPa m (Fig. 8). Indeed, observations of grounded outlet glaciers in West Greenland suggest that fully grounded glaciers undergo higher-magnitude dynamical changes than those with floating termini (McFadden and others, Reference McFadden, Howat, Joughin, Smith and Ahn2011). Furthermore, most GrIS outlet glacier fjord widths observed by Wood and others (Reference Wood2021) are similar to our 4 km narrow testbed (Fig. A7). Thus, knowledge of the glacier stress state is likely necessary to explain locally observed high-magnitude thinning.
Further evidence of the sensitivity of basally supported glaciers to grounding line retreat can be observed in the localized basal perturbation experiment. At testbed glaciers with high mean basal drag, pervasive thinning originating near the grounding line (as seen near year 10 in Fig. 4b) highlights this sensitivity. In contrast, testbeds with low basal stress (e.g. Fig. 4a) undergo the same magnitude of grounding line retreat yet lack this diffusive thinning. The potential for higher-stressed glaciers to undergo dramatic thinning echoes the modeled high sensitivity of the ice loss at the East Antarctic Ice Sheet to a basal thermal state transition, where inversions identify large basal areas with high basal drag (Dawson and others, Reference Dawson, Schroeder, Chu, Mantelli and Seroussi2022).
4.3 Longer-duration basal perturbations incite greater thickness changes
The localized basal perturbation experiment emulates two types of drainage efficiency (Moon and others, Reference Moon2014), which produce contrasting examples of dynamical thickness changes both near and far downstream of the perturbation. The Diffused Pulse, which is a basal drag reduction whose peak value is 10 times less than its transient counterpart, actually induces a larger magnitude of thickening/thinning immediately downstream/upstream of the perturbation. Furthermore, it prolongs the initial grounding line advance period, resulting in continued downstream thickening, which is particularly visible in wide testbeds (Fig. 5). These results emphasize the disproportionately larger impact of extended basal drag reduction on the glacier state.
The reasons for a long-lasting lower basal drag can be diverse. For instance, modeling of Helheim hydrology shows elevated pore pressure and low effective pressure during winter from frictional dissipation from high sliding speed (Sommers and others, Reference Sommers2023). A subglacial drainage system may fail to channelize due to insufficient meltwater discharge or lack of meltwater forcing variability (Schoof, Reference Schoof2010), or high ice-overburden pressure limits sizes of cavity (Doyle and others, Reference Doyle2014; de Fleurian and others, Reference de Fleurian2016), although the latter is more likely to occur in the accumulation zone where ice thickness is over 1 km. Additionally, multi-year inversions on surge glaciers experiencing thermal state switches triggered by surface meltwater have inferred basal drag changes on inter-annual timescales (Dunse and others, Reference Dunse2015; Gong and others, Reference Gong2018). The synthetic pulses spanning 0.1 and 2 years used in this study can also be interpreted as lower and upper bounds of timescale, and efficient drainage can develop over a variety of timescales (Vijay and others, Reference Vijay2021). Generally, the disproportionately larger impact from a long-lasting perturbation should not be overlooked. Additionally, previous investigations into the drainage system efficiency on flow dynamics have focused primarily on ice velocity patterns. We complement this knowledge by suggesting that, when interpreting the dynamic elevation change records, future studies should also consider the possible impact of prolonged basal lubrication even if the total magnitude of basal lubrication is relatively small.
4.4 Implications for diffusive thinning propagation
In our testbeds, mean basal drag level primarily and grounding line depth, to a lesser extent, control ice velocity (Table A2). For example, the narrow testbed with a high mean basal drag has a maximum flow speed of less than 1 km per year, which is only 30% of the speed of its low-mean-basal-drag counterpart. The speed at which the diffusive thinning propagates from the terminus roughly scales with how quickly diffusive thinning can propagate, which is typically 5–8 times the ice flow velocity (van de Wal and Oerlemans, Reference van de Wal and Oerlemans1995; van der Veen, Reference van der Veen2001). With high ice velocity due to low mean basal drag, longitudinal stretching rapidly transmits upstream and leads to widespread thinning. A similar mechanism has been proposed to explain far-reaching inland acceleration at Sermeq Kujalleq due to low basal drag (Bondzio and others, Reference Bondzio2017).
Previous studies (Felikson and others, Reference Felikson2017, Reference Felikson, Catania, Bartholomaus, Morlighem and Noël2021) have used Peclet numbers to identify large undulations in basal topography, known as ‘knickpoints’ as limits to upstream thinning propagation. Provided a simplified flux-geometry assumption, the derived Peclet numbers measure the relative importance between diffusion – which can migrate upstream – and downstream advection. While this offers a valuable static map view of where diffusive thinning diminishes, our simulations show that glacier dynamics conditioned by geometry and basal conditions determine the spatial extent of thinning on a decadal timescale, which may occur far downstream of major knickpoints in real-world glaciers (e.g. near the grounding line). Our results complement previous studies by suggesting that the glacier's dynamic state and its evolution can also play a considerable role in mapping upstream thinning extent. Furthermore, our simulations show that while glaciers with low mean basal drag can propagate diffusive thinning far inland, similar to gentle bed topography discussed in Felikson and others (Reference Felikson, Catania, Bartholomaus, Morlighem and Noël2021), glaciers with narrow fjords and higher mean basal drag levels can lose almost the same amount of mass during the same period (the smallest magenta dot in Fig. 8a), despite its strong thinning attenuation which concentrates behind the grounding line. The more delayed recovery of grounding line retreat after the front stops retreating suggests that these glaciers may have even higher mass loss potential (e.g. the black profile of Fig. 3d testbed at its new steady state).
4.5 Implications for ice-sheet modeling
Our work has useful implications for future modeling studies. We have shown in Fig. 3 that thinning magnitude depends sensitively on the sliding law, where the addition of ice overburden pressure feedback causes large variability in thinning. The choice of exponent in the sliding law may also add uncertainty to projected ice loss. To explore the effect of the exponent, we perform one additional overburden pressure experiment where we set m = 5, corresponding to a more plastic bed where an increase in sliding velocity has a more limited impact on the basal drag strengthening. Simulation results (Fig. A3) show that the thinning pattern and magnitude resemble more the Weertman case (without overburden pressure dependence), and the difference in grounding line migration from the control run in Figure 3 is negligible. This can also be seen from Eqn (3) where in the limit of perfect plasticity, i.e. m → ∞, the sliding law coefficient C remains constant and thus is effectively Weertman sliding law. This suggests substantial differences in ice mass loss projection due to the choice of the exponent alone in the same sliding law. Since Weertman and Budd's sliding law remain the most commonly employed sliding laws in glacier and ice-sheet scale modeling (e.g. Bondzio and others, Reference Bondzio2017; Goelzer and others, Reference Goelzer2020; Dawson and others, Reference Dawson, Schroeder, Chu, Mantelli and Seroussi2022) our results echo previous findings that sliding laws can critically influence ice mass loss projections (Brondex and others, Reference Brondex, Gagliardini, Gillet-Chaulet and Durand2017). Our work contributes to the knowledge by showing that in a wide range of glacier geometries and basal boundary conditions, grounding area change is a decent proxy for total dynamic thinning (Fig. 8a), and therefore grounding area movement can potentially be used as a constraint to calibrate the choices of sliding law when initializing large-scale ice-sheet models.
Additionally, it is important for studies using idealized glacier setups to be cautious when initializing glaciers with steady-state frontal geometries, such as fully grounded or floating termini. Our simulations reveal substantial thinning differences between glaciers with deep or shallow grounding lines (Fig. A4), which can bias the identification of primary controls suggested in Felikson and others (Reference Felikson, Nowicki, Nias, Morlighem and Seroussi2022), for instance. We advocate for future modeling studies to consider various dimensions of glacier geometries when constructing idealized models.
5. Conclusion
Our study explores the effect of ice overburden pressure and local basal slipperiness perturbations on dynamic thickness change of Greenland-like testbed glaciers, in an effort to constrain potential factors that may be driving dynamic thickness changes across Greenland glaciers.
We find that changes in both overburden pressure and basal slipperiness can induce dynamic thickness change which correlates well with grounding line migration. We find relationships between grounding line position and domain-wide thinning, and between front-to-grounding-line resistive stress loss and maximum thinning rate, but we find great variability from testbed to testbed in dynamic thinning rates despite consistent ice-front position histories. Thus, although ice-front position is readily observable, it should be used with caution for prediction or diagnosis of glacier dynamic thinning patterns.
We find changes in ice overburden pressure alone can be responsible for over 100 m of dynamic thinning as terminus continuously retreats over a decade, particularly at glaciers with narrow fjords and high basal drag levels. Basal lubrication perturbations have a diagnostic dipole shape that could be identified in maps of dh/dt. The time duration of a basal forcing has greater efficacy on surface elevation than its magnitude.
Finally, we find that on wavy-bedded glaciers, a uniform retreat of a calving front can produce episodic grounding line retreats, which manifest as short-duration undulations in dynamic elevation. In light of all these findings, we stress the importance of incorporating knowledge of bed topography, grounding line locations and stress estimates in any interpretation of observed dynamic thickness changes.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2024.50
Data
The scripts to run ISSM simulations and recreate the figures can be found on GitHub (https://github.com/alastairyang/ThinningTestbedPublic.git). The simulation output data are available on Zenodo (https://doi.org/10.5281/zenodo.10564805). ISSM is publicly available at https://issm.jpl.nasa.gov/.
Acknowledgements
This work was supported by the NASA Cryospheric Sciences grant ‘Integration of ICESat-2 Observations into Ice Sheet Elevation Change Records to Investigate Ice Sheet Processes’ (80NSSC21K0915). We thank Climate Data Toolbox (Greene and others, Reference Greene2019) and Scientific Color Maps (Crameri, Reference Crameri2018) for their data visualization tools.
Open access
