Tidewater glacier response to individual calving events

Abstract

Tidewater glaciers have been observed to experience instantaneous, stepwise increases in velocity during iceberg-calving events due to a loss of resistive stresses. These changes in stress can potentially impact tidewater glacier stability by promoting additional calving and affecting the viscous delivery of ice to the terminus. Using flow models and perturbation theory, we demonstrate that calving events and subsequent terminus readvance produce quasi-periodic, sawtooth oscillations in stress that originate at the terminus and propagate upstream. The stress perturbations travel at speeds much greater than the glacier velocities and, for laterally resisted glaciers, rapidly decay within a few ice thickness of the terminus. Consequently, because terminus fluctuations due to individual calving events tend to be much higher frequency than climate variations, individual calving events have little direct impact on the viscous delivery of ice to the terminus. This suggests that the primary mechanism by which calving events can trigger instability is by causing fluctuations in stress that weaken the ice and lead to additional calving and sustained terminus retreat. Our results further demonstrate a stronger response to calving events in simulations that include the full stress tensor, highlighting the importance of accounting for higher order stresses when developing calving parameterizations.

1. Introduction

As tidewater glaciers flow into the ocean they lose mass via submarine melting and iceberg calving. These processes, collectively referred to as frontal ablation, influence tidewater glacier stability by modifying the magnitude and distribution of stresses within a glacier. Moreover, iceberg calving is itself a function of near-terminus stresses, and therefore near-terminus stresses must be carefully accounted for in calving parameterizations used in long-timescale glacier simulations (e.g. Bassis and Walker, 2012; Mercenier and others, 2018; Schlemm and Levermann, 2019). Several recent studies have suggested that submarine melting, which evolves in response to changing ocean conditions and subglacial discharge (e.g. Jenkins, 2011; Motyka and others, 2013), can affect calving rates by shaping glacier termini and continuously modifying near-terminus stresses (e.g. O'Leary and Christoffersen, 2013; Cook and others, 2014; Todd and Christoffersen, 2014; Krug and others, 2015; Cowton and others, 2019; Ma and Bassis, 2019). Calving events also shape glacier termini and modify near-terminus stresses. However, unlike submarine melting, iceberg calving is a discrete process that can potentially produce large, transient fluctuations in stress that propagate upstream.

Full-glacier thickness calving events at some of Greenland's largest glaciers have been observed to cause nearly instantaneous, stepwise increases in glacier velocity of up to 30% followed by a gradual decay back to pre-calving velocities during terminus readvance (Amundson and others, 2008; Nettles and others, 2008; Rosenau and others, 2013; Cassotto and others, 2019; Kane and others, 2020). In addition, the sensitivity of these glaciers to tidal forcings along their termini is impacted by changes in the near-terminus stress state following calving events (de Juan and others, 2010). The magnitude of the velocity response is not directly related to iceberg size, as large, multiple-kilometer scale calving events from ice shelves may produce little change in velocity (Hill and others, 2018) while relatively small calving events from grounded termini can produce large responses if the calving events occur from regions of the glacier that are particularly prone to rapid stress redistribution (Cassotto and others, 2019).

Thinning disturbances that originate near tidewater glacier termini can propagate upstream and lead to flow instability and retreat, particularly for glaciers that are close to flotation and that have overdeepened beds (e.g. Pfeffer, 2007; Felikson and others, 2017). Whether individual calving events can trigger tidewater glacier instability through this mechanism remains unclear. By changing a glacier's stress state, calving events increase both ice velocities and strain rates. Increases in velocity cause advective thickening due to rapid delivery of ice toward the terminus, which should promote stability by increasing traction and reducing strain rates, whereas increases in strain rates lead to dynamic thinning and destabilization. How these competing processes compare, and evolve in time, will dictate the longer term impact of calving events on glacier stability. In this study, we use glacier flowline models and a perturbation analysis (based on the work of Williams and others, 2012) to begin developing a mechanistic understanding of the factors controlling how tidewater glaciers respond to calving events while also providing some insights into the ways in which transient, calving-induced stress fluctuations might impact long-timescale glacier behavior.

2. Model description

In order to investigate glacier response to calving events, we run simulations using both (1) the full-Stokes capabilities of Elmer/Ice in a set-up similar to Gladstone and others (2017) and (2) the shallow shelf approximation (SSA) (code adapted from Enderlin and others, 2013). Comparison of these models allows us to address the impacts of higher order stresses on glacier response. In addition, the SSA lends itself to a perturbation analysis that we use to help interpret the model results.

We use a similar set-up for both models. The model domain consists of the lowermost 10 km of a 5-km wide glacier that has a water depth of 600 m at the terminus and a constant bed slope that can be prograde, flat or retrograde. We model ice flow along the central flowline using the Stokes equations for a viscous fluid in which the rheology is described by Glen's flow law (Cuffey and Paterson, 2010) and for simplicity we assume that the ice is temperate everywhere. Later we will discuss how this choice of temperature affects the model results. Lateral stresses are accounted for by integrating from the glacier margins to the centerline (Gagliardini and others, 2010; van der Veen, 2013). Basal shear stresses are assumed to depend on velocity and effective pressure (ice-overburden pressure minus pore-water pressure) (Budd and others, 1979; Gladstone and others, 2017); we assume an efficient subglacial hydrological system and therefore the phreatic surface is horizontal and at sea level. The velocity at the upstream end of the domain is set to 4000 m a⁻¹ and, since we are only modeling the lower reaches of a glacier, the surface mass-balance rate is set to −2 m a⁻¹ everywhere. The model domain and parameter values were chosen in order to produce glacier geometries and flow speeds that roughly mimic Greenland outlet glaciers (e.g. Catania and others, 2020, and references therein).

2.1. Full Stokes model

We define a coordinate system in which x is the horizontal coordinate in the down-glacier direction, with x = 0 corresponding to the initial terminus position, and y is the elevation relative to sea level. The velocity vector is u = 〈u _x,  u _y〉. We use the Cauchy stress tensor $\sigma _{ij} = \sigma ^{\prime }_{ij}-P\delta _{ij}$, where $\sigma ^\prime _{ij}$ is the deviatoric stress tensor, P is the isotropic pressure and δ_ij is the Kronecker delta. Conservation of mass and momentum dictate that

(1)$${\partial u_i\over \partial x_i} = 0$$

and

(2)$${\partial \sigma_{ij}\over \partial x_j}-f_i-\rho g\delta_{iy} = 0,\; $$

where f is a body force that is used to parameterize lateral drag, ρ is the density of ice, g is the (scalar) gravitational acceleration and Einstein notation is used. The deviatoric stress is related to the strain rate by Glen's flow law:

(3)$$\sigma^\prime_{ij} = A^{-1/3}\dot\varepsilon_{\rm e}^{-2/3}\dot\varepsilon_{ij},$$

where A = 75 MPa⁻³ a⁻¹ is the flow law parameter for temperate ice (Cuffey and Paterson, 2010), $\dot \varepsilon _{ij} = ( \partial u_i/\partial x_j + \partial u_j/\partial x_i) /2$ is the strain rate, $\dot \varepsilon _{\rm e} = ( \dot \varepsilon _{ij}\dot \varepsilon _{ij}/2) ^{1/2}$ is the effective strain rate and we have assumed a flow law exponent of n = 3. We use a heuristic parameterization for the lateral resistance; the parameterization is derived by considering the balance of forces in a long, weak-bedded ice stream or ice shelf and neglecting longitudinal stresses (Gagliardini and others, 2010), which results in

(4)$$f_i = {1\over W}\left({4\over AW}\right)^{1/3}\left\vert u \right\vert ^{-2/3}u_i,\; $$

where W is the glacier half-width.

Along the glacier surface and the subaerial portion of the terminus we invoke a traction-free boundary condition:

(5)$$\sigma_{ij}\hat{n}_j = 0,$$

where $\hat {n}$ is the outward pointing surface unit normal vector. For grounded portions of the glacier, (1) the bed-parallel shear stress σ_b depends on the effective pressure N and opposes the bed-parallel velocity u _b and (2) the component of the velocity that is perpendicular to the bed must equal zero:

(6)$$\left.\eqalign{ \sigma_{\rm b} = \hat{t}_i( \sigma_{ij}\hat{n}_j) = -\beta N \vert u_{\rm b}\vert ^{-2/3}u_{\rm b} \cr u_i\hat{n}_i = 0 }\right\} x\leq x_{\rm g},\; $$

where $\hat {t}$ is the tangential unit vector that points in the downstream direction, β is the basal roughness factor, N is the effective pressure and x _g is the location of the grounding line. We use a constant basal roughness factor of β = 0.0022 m^−1/3 a^1/3, which yields basal shear stresses in the range of 10–50 kPa (consistent with inversions for basal shear stresses beneath tidewater glaciers; e.g. Enderlin and others, 2018). For portions of the glacier that are in contact with the ocean, the shear stress is zero and the normal stress equals the hydrostatic pressure of the seawater:

(7)$$\left.\eqalign{ \hat{t}_i( \sigma_{ij}\hat{n}_j) = 0 \cr \hat{n}_i( \sigma_{ij}\hat{n}_j) = \rho_{\rm w}gy }\right\} x > x_{\rm g},\; $$

where ρ_w is the density of seawater. Finally, the velocity at x = −10 km (10 km upstream from the terminus) is set to $u = \langle {4000\, \hbox {m\ a}^{-1}, \, 0\rangle }$, independent of depth.

The glacier surface evolves according to

(8)$${\partial h_{\rm s}\over \partial t} + u_x{\partial h_{\rm s}\over \partial x} = u_y + \dot{b},$$

and, for portions of the glacier that are floating, the bottom of the glacier evolves according to

(9)$${\partial h_{\rm b}\over \partial t} + u_x{\partial h_{\rm b}\over \partial x} = u_y,\; $$

where h _s and h _b are the surface and bed elevations, $\dot {b}$ refers to the width-averaged surface mass-balance rate (we set the basal mass-balance rate to 0), and 〈u _x,  u _y〉 corresponds to the velocity either at the surface or at the bottom of the glacier. After model spin-up and between calving events (see Section 2.3) the terminus is free to advance at a rate dictated by the velocity, and the grounding line is tracked by solving a contact problem (Favier and others, 2012).

2.2. Shallow shelf approximation

Assuming that velocities and strain rates are independent of depth and vertically integrating the momentum and mass conservation equations yields the commonly used shallow shelf/stream approximation (SSA) (MacAyeal, 1989). Additional transverse integration produces 1-D forms of these equations that are suitable for investigating tidewater glacier behavior but tend to underestimate tidewater glacier velocities because they do not fully capture the complex stress fields observed near their termini (Hindmarsh, 2012). Despite this concern, the SSA remains the dominant model used for long-timescale simulations (years to decades) due to its computational efficiency. Our goal here is twofold: (1) by comparing full-Stokes and SSA simulations we can quantify how well the SSA model performs when simulating transient changes in terminus position and (2) the simplified equations of the SSA lend themselves to a perturbation analysis that provides insights into how glaciers respond to calving events and other stress perturbations. We later compare results from the perturbation analysis to the SSA simulations.

The 1-D momentum equation along the centerline is given by (van der Veen, 2013)

(10)$$\eqalign{2{\partial\over \partial x}\left(HA^{-1/3}\left\vert {\partial U\over \partial x}\right\vert ^{-2/3}{\partial{U}\over \partial x}\right)- \beta N\left\vert U\right\vert ^{-2/3}U \cr- {H\over W}\left({4\over AW}\right)^{1/3}\vert U\vert ^{-2/3}U = \rho gH{\partial h_{\rm s}\over \partial x},\;} $$

where H is the glacier thickness and U is the depth-averaged centerline velocity in the x-direction. Although here the flow rate factor A is depth-averaged, we are assuming that the ice is temperate everywhere and therefore the value of A in the SSA simulations is the same as in the full-Stokes simulations. A Dirichlet boundary condition is used to prescribe a velocity of U = 4000 m a⁻¹ at x = −10 km (as also done in the full-Stokes simulations) and a Neumann boundary condition is prescribed at the terminus (x = L) by subtracting the depth-integrated hydrostatic pressure from the depth-integrated glaciostatic pressure, converting the resulting resistive stress into a deviatoric stress, and inserting the result into Glen's flow law to yield

(11)$$\left.{\partial U\over \partial x}\right\vert _{x = L} = A\left[{\rho gH_L\over 4}\left(1-{\rho_{\rm w}\over \rho}{D^2\over H_L^2}\right)\right]^3,\; $$

where H _L is the terminus thickness and D is the submerged depth of the terminus.

The glacier thickness evolves according to the associated 1-D mass continuity equation for a glacier with constant width:

(12)$${\partial H\over \partial t} = \dot{b} - {\partial( UH) \over \partial x}.$$

A moving grid is used to track the grounding line, and as with the full-Stokes simulations the terminus is free to advance at its flow rate after model spin-up and between calving events.

2.3. Simulations

Both the full-Stokes and SSA models were run to steady-state with a fixed length of L = 10 km for three different bed slopes (prograde, flat and retrograde). After spinning up the models, we simulated calving by removing the lowermost 200 m of the glacier and then allowing the terminus to re-advance to its pre-calving position. We tested the effect of calving event size on glacier response and found that, although this changed the overall magnitude of the glacier response, it did not affect the key results. In order to relate our results to what would be expected from stress-perturbations at the terminus, we also ran simulations with the SSA model in which (1) the terminus fluctuated around a mean length of L = 10 km in response to calving events and glacier flow and (2) we kept the terminus position fixed but applied a periodic forcing to the terminus (Section 5).

3. Results

In the full-Stokes simulations, calving events produced instantaneous changes in velocity and strain rate that were greatest near the glacier termini and rapidly decayed over a distance of a few kilometers (Fig. 1), with the strain rate perturbations decaying more rapidly than the velocity perturbations. All three glacier geometries experienced increases in velocity of close to 20% in the near terminus region, whereas the strain rate (plotted here as a gradient in the depth-averaged horizontal velocity) increased by about a factor of two. The highest velocities and strain rates occurred immediately following the calving events. As the termini re-advanced the velocities and strain rates gradually evolved back toward their pre-calving values. The largest response occurred for the simulation with the retrograde bed, but the perturbation decay length was greatest for the prograde bed.

Fig. 1. Results from the full-Stokes simulations. (a)–(c) Steady-state velocity fields for the retrograde, flat and prograde beds and corresponding fractional changes in (d)–(f) depth-averaged velocity and (g)–(i) gradient in depth-averaged velocity following a 200 m long calving event. The open triangles in (d)–(i) indicate the observed e-folding lengths of the velocity and strain rate perturbations, and the filled triangles indicate the e-folding length predicted by the perturbation analysis (Section 5).

Glacier response to calving events is dominated by changes in strain rate, which can be clearly seen by splitting the glacier thinning rates (Eqn (12)) into advective thickening and dynamic thinning terms:

(13)$${\partial H\over \partial t} \approx - U{\partial H\over \partial x} - H{\partial U\over \partial x},\; $$

recalling that U is the depth-averaged horizontal velocity and noting that we omitted the surface mass-balance rate because it is much smaller than the other terms during these transient events. Changes in ice velocity have very little impact on the advective thickening rate, whereas changes in strain rate cause the dynamic thinning rate to increase by about a factor of two (Fig. 2).

Fig. 2. Evolution of the (a)–(c) advective thickening rate (d)–(f) dynamic thinning rate, and (g)–(i) rate of thickness change for the full-Stokes simulations following a 200 m calving event. The top, middle and bottom rows correspond to the retrograde, flat and prograde bed geometries.

Simulations run using the SSA model exhibited similar, but less pronounced, behavior than those run with the full-Stokes model (Fig. 3). The two models had similar initial velocity and thickness profiles, but the full-Stokes model yielded greater strain rates in the near-terminus region. The gradient in the depth-averaged velocity, ∂U/∂x, reached 0.87 a⁻¹ in the full-Stokes simulations but just 0.57 a⁻¹ in the SSA simulations. As a result, the terminus was also thinner in the full-Stokes simulations (even being superbuoyant) and had a steeper surface slope. Due to these differences in geometry and strain rate, the velocity, strain rate and deviatoric stress perturbations from the calving events were $\sim \!50\percnt$ larger in the full-Stokes simulations than in the SSA simulations but also had shorter decay lengths. For both, these changes in stress are on the order of 10% of the pre-calving stress.

Fig. 3. Comparison of the full-Stokes and SSA simulations for a flat bed and 200 m long calving event: (a), (b) depth-averaged velocity and thickness profiles prior to the calving event and (c), (d) changes in velocity and depth-averaged longitudinal deviatoric stress relative to the pre-calving values (profile spacing is 0.5 d and the profiles span 7 d). The open triangles in (c) indicate the observed e-folding lengths of the velocity perturbations, and the filled triangles indicate the e-folding length predicted by the perturbation analysis (Section 5). We observed similar discrepancies between the full-Stokes and SSA simulations when using sloping beds.

4. Fluctuations in near-terminus stresses

Our model results highlight the impact of near terminus stresses on variations in tidewater glacier flow. As a glacier terminus calves and re-advances, it experiences periodic, sawtooth variations in longitudinal stresses due to changes in lateral and basal shear stresses, driving stress and the terminus stress balance. To further demonstrate this, we ran two sets of additional simulations with the SSA model. In the first, the glacier experienced calving events of 200 m in length and the terminus oscillated around a mean position of x = 10 km (‘calving simulation’). In the second, we kept the terminus position fixed but imposed a sinusoidally varying stress at the terminus (‘fixed-terminus simulation’).

In order to determine the amplitude of the equivalent stress oscillations to impose in the fixed-terminus simulation, we first rearrange Eqn (10) and integrate to find the deviatoric longitudinal stress at x = 10 km in the calving simulation right as the terminus is about to calve from an advanced position of x = c:

(14)$$\tilde\sigma_{xx}^{\prime} = {1\over \tilde H}\left(H_L\sigma_{xx, L}^\prime - \int_{\tilde x}^c \tau\, {\rm d}x\right),\; $$

where $\sigma _{xx}^\prime$ is the deviatoric longitudinal stress, tildes are used to indicate values at x = 10 km, and as before subscript L refers to terminus values. For convenience we have defined

(15)$$\eqalign{\tau \equiv {1\over 2} \left(\beta N\left\vert U\right\vert ^{-2/3}U + {H\over W}\left({4\over AW}\right)^{1/3}\vert U\vert ^{-2/3}U + \rho gH{\partial h_{\rm s}\over \partial x}\right).}$$

The deviatoric longitudinal stress at the terminus is the term in brackets in Eqn (11), i.e.

(16)$$\sigma_{xx, L}^\prime = {1\over 4}\rho g H_L\left(1-{\rho_{\rm w}\over \rho}{D^2\over H_L^2}\right).$$

Thus, we determine the amplitude of the imposed stress oscillations by subtracting the deviatoric longitudinal stress at x = 10 km when the terminus has advanced to x = c from the deviatoric longitudinal stress at the same location at the end of the model spin-up (i.e. when the glacier is in a steady-state). The period of the oscillations is determined by the periodicity of calving events in the calving simulation, which in this case is 13 d.

The variations in deviatoric stress at the terminus for the two simulations are shown in Figure 4a. The magnitude of the stress variations imposed on a fixed terminus are somewhat larger than those for the terminus that advances and retreats because the imposed stress variations account for drag, whereas in the calving simulation the deviatoric stress at the terminus is simply due to changes in ice thickness and water depth. The fact that these two curves have similar amplitude indicates that the stress balance is primarily being affected by changes in the balance of glaciostatic and hydrostatic stresses along the terminus. The changes in velocity and deviatoric stress that occur in the fixed-terminus simulation are lower than those in the calving simulation (Figs 4b, c). However, due to the nature of the forcing (sinusoidal vs sawtooth) the velocities and strain rates remain high for a longer period of time in the fixed-terminus simulation. The agreement between these two simulations, especially with respect to the decay length, provides further support for the notion that glacier response to calving events can be modeled as periodic oscillations in near-terminus stresses. This may provide an avenue for incorporating the time-varying effects of iceberg calving into calving parameterizations without needing to model individual calving events.

Fig. 4. Variations in flow in the SSA model due to periodic forcing at the terminus and calving events: (a) deviatoric stress at the terminus (t = 0 is the first time step after a calving event), (b) changes in depth-averaged velocity and (c) changes in depth-averaged deviatoric stress. In (b) and (c), the differences are taken relative to the most advanced position in the calving scenario, the calving profiles represent the days following a calving event, the profile spacing is 0.5 d, and the profiles span 7 d. In (b), the open triangles indicate the observed e-folding lengths of the velocity differences, and the filled triangle indicates the value predicted from the perturbation analysis (Section 5).

5. Tidewater glacier response to periodic terminus forcings

Removal of ice via iceberg calving, and subsequent terminus readvance prior to the next calving event, represents a quasi-periodic perturbation to the glacier stress balance. Given that the flow response to calving can be approximately characterized with a harmonic stress perturbation at the terminus, we apply perturbation theory to the SSA model (see Appendix A) in order to better understand what controls the spatial and temporal scales of calving-induced flow perturbations.

Our work closely follows that of Williams and others (2012), who modeled the frequency response of ice streams, except that we include effective pressure in our basal shear stress parameterization and we simultaneously include the effects of both basal and lateral shear stresses (Williams and others (2012) treated basally and laterally resisted ice streams separately). The perturbation analysis involves (1) assuming a horizontal bed, (2) applying a linear perturbation to Eqns (10) and (12), (3) expressing the perturbed equations in terms of periodic forcings in thickness, velocity and strain rate, all of which are assumed to have the same frequency and wavenumber and (4) re-casting the perturbed momentum and mass continuity equations in terms of an algebraic equation that relates the wavenumber to the forcing frequency. The wavenumber varies spatially as a result of spatial variations in glacier geometry and flow. We select characteristic values of geometry and flow and compute the wavenumber for various forcing frequencies, which we then use to determine the decay length, wavelength and phase velocity. We compare the results to perturbations that are applied to the SSA model.

In Figure 5, we show how the decay length, phase velocity and wavelength depend on frequency as predicted by the perturbation analysis for the spun-up SSA model. We then perturb the SSA model by applying a sinusoidally varying deviatoric stress at the terminus and compute the resulting decay length, phase velocity and wavelength. We find good agreement between the perturbation analysis and the SSA simulations. For frequencies associated with calving events from fast flowing tidewater glaciers (i.e. periods of days to weeks), the perturbation analysis and SSA simulations both suggest that the decay length is on the order of a few kilometers and is nearly constant for frequencies $\gtrsim$10⁻² d⁻¹ (i.e. those associated with calving events), the phase velocity is 1–2 orders of magnitude faster than the glacier flow speeds and the wavelength approaches typical glacier lengths. The perturbation analysis systematically overpredicts the decay length by ~30% (see also Figs 3 and 4), which we attribute to (1) selecting characteristic values of the velocity, strain rate and thickness and therefore not accounting for their spatial variations and (2) selecting these characteristic values at the terminus (where velocities and strain rates are high and the ice is at, or is approaching, flotation).

Fig. 5. Decay length, phase velocity and wavelength for periodic forcings applied to the glacier terminus. PA and SSA represent the results from the perturbation analysis and from applying a periodic forcing to the SSA model, respectively. The black-dashed line in (a) indicates the high-frequency limit of the decay length.

We analyze the perturbation analysis in the high-frequency limit and find that as the angular frequency ω → ∞, the decay length D _L converges to

(17)$$D_L \rightarrow \left({U_0\over \dot\varepsilon_0}\right)^{1/3} \sqrt{{2\over A^{1/3}\beta\rho g( 1-( {\rho_{\rm w}}/{\rho}) ( {D}/{H_0}) ) + ( {1}/{W}) ( {4}/{W}) ^{1/3}}},\; $$

where subscript 0 refers to characteristic values selected from the datum state. The wavelength and phase velocity diverge (become infinitely large) in the high-frequency limit. To evaluate Eqn (17), we use the terminus velocity and thickness and compute the characteristic strain rate by using the terminus boundary condition (Eqn (11)). The theoretical decay length is essentially identical to that of Walters (1989), who applied a perturbation analysis to quantify tidewater glacier response to tides. That analysis differs from ours in that it (1) used a different parameterization for basal sliding, (2) neglected lateral drag and (3) treated the perturbations as step changes in stress and did not address the dependency on frequency. Our analysis extends that of Walters (1989) by suggesting that the decay length equation is applicable for forcings with periods $\lesssim$100 d. Thus the theoretical decay limit provides a useful metric for assessing how a fast-flowing glacier will respond to calving events and other high-frequency perturbations.

The high-frequency decay limit depends on the ratio of ice velocity to strain rate, basal resistance (basal roughness factor and proximity to flotation), glacier width (Fig. 6) and ice stiffness. For narrow glaciers, lateral shear stresses cause the decay length to be short (i.e. a few kilometers) and insensitive to the proximity to flotation. For these glaciers the decay length is primarily a function of the ratio of the ice velocity to the strain rate. This ratio will typically be large for flat, weak-bedded glaciers (high velocities and low strain rates) and small for steep, rapidly deforming glaciers that have high strain rates. As glacier width increases the decay length also increases (up to tens of kilometers) and becomes more sensitive to both the velocity–strain rate ratio and the proximity to flotation. Finally, Eqn (17) indicates that the decay length depends inversely on the flow rate factor, which depends strongly on temperature (Cuffey and Paterson, 2010). The decay length is greatest for cold, stiff ice. All else being equal, a decrease in ice temperature from 0 to −10^°C results in an increase in the decay length of ~40%.

Fig. 6. Dependence of the high-frequency decay length on the proximity to flotation (D/H ₀), the velocity–strain rate ratio at the terminus ($U_0/\dot {\varepsilon }_0$) and the glacier width. In (a)–(c), the glacier width is 5, 50 and 500 km, respectively.

The large-phase velocities and wavelengths that we have calculated suggest that changes in stress that originate at the terminus—both during calving events and subsequent readvance—are transmitted upstream almost instantaneously. Moreover, the short-decay lengths of fast-flowing glaciers mean that calving-induced stress fluctuations do not propagate far upstream and that the glacier quickly readjusts back toward its pre-calving state as it readvances and gains traction. Thus, we suggest that the primary mechanism by which individual calving events can affect long-term glacier stability is by fatiguing the near-terminus ice, similar to the fatiguing of ice shelves by ocean waves (e.g. Sergienko, 2010), leading to increased rates of calving and sustained terminus retreat. Our model predicts that near-terminus stresses may fluctuate by ~10% around the mean background stress state. The effect of these stress fluctuations on ice stiffness, microcracks and fracture propagation and calving remains uncertain.

6. Conclusions

Tidewater glaciers experience quasi-periodic fluctuations in stress resulting from rapid geometric changes associated with iceberg-calving events. Using idealized glacier geometries designed to mimic narrow (≤5 km wide) tidewater glaciers, we simulated the response to calving events in both full-Stokes and SSA models by removing the lowermost portion of the glaciers and observing the resulting changes in velocity and strain rate. In the full-Stokes simulations the near-terminus velocities temporarily increased by ~20%, but the strain rates increased by more than 100%. As the termini re-advanced the velocity and strain rates decayed back toward their initial, pre-calving values. The decay length of the velocity and strain rate perturbations was always within a few ice thicknesses of the termini, with the strain rate perturbations having a shorter decay length than the velocity perturbations. Glaciers that calved down a retrograde bed exhibited shorter decay lengths than those that calved on flat or seaward-sloping beds, likely because larger thickness and velocity gradients of glaciers on retrograde beds limits the ability of rapid, near-terminus stress fluctuations to reach far upstream. Glaciers on retrograde beds also experienced larger changes in flow than those that calved on flat or seaward sloping beds, although differences in flow response become small when normalized by the pre-calving flow rates. The SSA model produced similar patterns of speed up and subsequent slowdown, but the velocity and strain rate perturbations were smaller than those observed in the full-Stokes simulations. We suggest that the enhanced response of the full-Stokes model is a result of the higher order terms contributing to weaker ice and causing the glacier to be closer to floatation.

Glacier response to calving is a result of changes in the glacier force balance. Removal of ice causes the deviatoric stress at a fixed Eulerian point to increase due to a loss of resistive stresses and an increase in the deviatoric stress at the terminus (terminus retreat produces a thicker terminus). As the glacier readvances the deviatoric stress at the terminus decreases and the glacier regains traction. These stress changes, which can have an amplitude on the order of 10% of the background stress, are transmitted upstream nearly instantaneously.

To demonstrate that iceberg calving and subsequent terminus readvance can be viewed as quasi-periodic stress fluctuations at the terminus, we compared flow variations of a glacier that was allowed to advance and retreat to one that had a fixed terminus position that was subjected to sinusoidally varying stresses. These simulations exhibited very similar behavior. We further showed that at high frequencies, such as those associated with calving events, glacier response to calving events can be described to zeroth-order by a perturbation analysis of the SSA approximation. We find that the decay length of narrow tidewater glaciers is strongly controlled by lateral shear stresses and is relatively insensitive to the proximity to flotation. For wider glaciers the range of possible decay lengths increases greatly and the decay length becomes highly sensitive to a glacier's proximity to flotation. Furthermore, the decay length depends on ice stiffness, with perturbations being felt farther upstream for glaciers composed of cold, stiff ice.

Our analysis provides an important step toward assessing the impact of individual calving events on tidewater glacier stability. Stress fluctuations associated with calving events have short-decay lengths and are transmitted at speeds much greater than the ice velocity. These fluctuations appear to have little direct impact on the delivery of ice to the terminus because the glacier quickly readjusts its stress field as the terminus readvances following a calving event. Thus the primary mechanism by which calving-generated stress fluctuations can affect glacier stability is by changing the near terminus stresses, leading to the weakening of ice through the development of microcracks and fractures (e.g. Krug and others, 2014; Mercenier and others, 2019). If this weakening leads to subsequent calving events prior to the terminus readvancing to its pre-calving position, then the stress fluctuations associated with individual calving events can contribute to lower frequency flow variability (e.g. Riel and others, 2021) that may be out of sync with climate variations and, if so, help to trigger glacier instability by bringing the glacier surface into a warmer climate.

We note that we conducted flowline simulations, used a heuristic parameterization of lateral shear stresses, and did not account for small-scale topography. Simulations of changes in stress due to submarine melting (e.g. Vallot and others, 2018) indicate that the location of submarine melt plumes is important. Melting along the glacier centerline appears to have less impact on stress, and therefore calving, then melting along the margins (Cowton and others, 2019). We suspect that the location of calving events will have similar effects on stress distribution and transmission and may partly explain why our simulations are unable to explain the nonlinear glacier response to calving events observed by Cassotto and others (2019), which likely also requires detailed knowledge of the bedrock topography. Moreover, a full 3-D model would likely yield higher lateral shear stresses than predicted by the heuristic parameterization, resulting in greater rates of deformation and thinner ice. We expect that this would produce a stronger response to calving events than in our simulations, akin to the full-Stokes model producing a stronger response than the SSA model in our flowline simulations. Future work should investigate the impact of the spatial distribution of calving events (both horizontal as well as vertical) on glacier response and to account for high-frequency stress fluctuations when developing calving parameterizations that are suitable for long-timescale simulations.