Hostname: page-component-89b8bd64d-9prln Total loading time: 0 Render date: 2026-05-09T06:09:32.313Z Has data issue: false hasContentIssue false
  • English
  • Français

Tidewater glacier response to individual calving events

Published online by Cambridge University Press:  08 April 2022

Jason M. Amundson Open the ORCID record for Jason M. Amundson [Opens in a new window] ,
Martin Truffer Open the ORCID record for Martin Truffer [Opens in a new window]  and
Thomas Zwinger Open the ORCID record for Thomas Zwinger [Opens in a new window]
Jason M. Amundson*
Affiliation:
Department of Natural Sciences, University of Alaska Southeast, Juneau, AK, USA Institute for Atmospheric and Earth System Research, University of Helsinki, Helsinki, Finland
Martin Truffer
Affiliation:
Geophysical Institute, University of Alaska Fairbanks, Fairbanks, AK, USA
Thomas Zwinger
Affiliation:
CSC-IT Center for Science, Espoo, Finland
*
Author for correspondence: Jason M. Amundson, E-mail: jmamundson@alaska.edu
Rights & Permissions [Opens in a new window]

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.

Keywords

Glacier modelingice dynamicsiceberg calving

Information

Type
Article
Information
Journal of Glaciology , Volume 68 , Issue 272 , December 2022 , pp. 1117 - 1126
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
Copyright © The Author(s), 2022. Published by Cambridge University Press

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, Reference Bassis and Walker2012; Mercenier and others, Reference Mercenier, Lüthi and Vieli2018; Schlemm and Levermann, Reference Schlemm and Levermann2019). Several recent studies have suggested that submarine melting, which evolves in response to changing ocean conditions and subglacial discharge (e.g. Jenkins, Reference Jenkins2011; Motyka and others, Reference Motyka, Dryer, Amundson, Truffer and Fahnestock2013), can affect calving rates by shaping glacier termini and continuously modifying near-terminus stresses (e.g. O'Leary and Christoffersen, Reference O'Leary and Christoffersen2013; Cook and others, Reference Cook2014; Todd and Christoffersen, Reference Todd and Christoffersen2014; Krug and others, Reference Krug, Durand, Gagliardini and Weiss2015; Cowton and others, Reference Cowton, Todd and Benn2019; Ma and Bassis, Reference Ma and Bassis2019). 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, Reference Amundson2008; Nettles and others, Reference Nettles2008; Rosenau and others, Reference Rosenau, Schwalbe, Maas, Baessler and Dietrich2013; Cassotto and others, Reference Cassotto2019; Kane and others, Reference Kane2020). 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, Reference de Juan2010). 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, Reference Hill, Gudmundsson, Carr and Stokes2018) 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, Reference Cassotto2019).

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, Reference Pfeffer2007; Felikson and others, Reference Felikson2017). 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, Reference Williams, Hindmarsh and Arthern2012) 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 (Reference Gladstone2017) and (2) the shallow shelf approximation (SSA) (code adapted from Enderlin and others, Reference Enderlin, Howat and Vieli2013). 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, Reference Cuffey and Paterson2010) 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, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010; van der Veen, Reference van der Veen2013). Basal shear stresses are assumed to depend on velocity and effective pressure (ice-overburden pressure minus pore-water pressure) (Budd and others, Reference Budd, Keage and Blundy1979; Gladstone and others, Reference Gladstone2017); 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−1 and, since we are only modeling the lower reaches of a glacier, the surface mass-balance rate is set to −2 m a−1 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, Reference Catania, Stearns, Moon, Enderlin and Jackson2020, 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−3 a−1 is the flow law parameter for temperate ice (Cuffey and Paterson, Reference Cuffey and Paterson2010), $\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, Reference Gagliardini, Durand, Zwinger, Hindmarsh and Le Meur2010), 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 a1/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, Reference Enderlin, O'Neel, Bartholomaus and Joughin2018). 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, Reference Favier, Gagliardini, Durand and Zwinger2012).

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, Reference MacAyeal1989). 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, Reference Hindmarsh2012). 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, Reference van der Veen2013)

(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−1 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−1 in the full-Stokes simulations but just 0.57 a−1 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 (Reference Williams, Hindmarsh and Arthern2012), 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 (Reference Williams, Hindmarsh and Arthern2012) 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−2 d−1 (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 (Reference Walters1989), 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 (Reference Walters1989) 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, Reference Cuffey and Paterson2010). 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 0), 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, Reference Sergienko2010), 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, Reference Krug, Weiss, Gagliardini and Durand2014; Mercenier and others, Reference Mercenier, Lüthi and Vieli2019). 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, Reference Riel, Minchew and Joughin2021) 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, Reference Vallot2018) 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, Reference Cowton, Todd and Benn2019). 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 (Reference Cassotto2019), 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.

Acknowledgements

T. Zwinger was supported by Academy of Finland grant number 322978. We are grateful for the thoughtful and constructive comments provided by two anonymous reviewers.

Appendix A.

Here we use a perturbation analysis to quantify the frequency response of tidewater glaciers, with the aim of assessing the temporal and spatial scales over which calving events and other high frequency forcings affect ice flow. The derivation closely follows that of Williams and others (Reference Williams, Hindmarsh and Arthern2012), who investigated the frequency response of ice streams. We consider a glacier that rests on a horizontal bed, such that ∂h s/∂x = ∂H/∂x, and has a constant width. We define the strain rate as $\dot \varepsilon \equiv \partial U/\partial x$. To further simplify the analysis, we also assume that ice flows in the positive x-direction, the terminus is located at x = 0 and the flow is extensional and consequently |U|−2/3U = U 1/3 and $\left \vert \dot \varepsilon \right \vert ^{-2/3}\dot \varepsilon = \dot \varepsilon ^{1/3}$. Thus, the SSA stress balance and mass continuity equations simplify to

(A1)$$\Omega_a{\partial\over \partial x}\left(H\dot\varepsilon^{1/3}\right)- \left[\beta\left(H-{\rho_{\rm w}\over \rho_i}D\right) + \Omega_{b}H \right]U^{1/3} = H{\partial H\over \partial x},\; $$

and

(A2)$${\partial H\over \partial t} = \dot{B}-{\partial ( HU) \over \partial x},\; $$

where we have defined the constants Ωa ≡ 2A −1/3(ρi g) −1 and Ωb ≡ Ωa4−1/6W −4/3.

We now consider a linear perturbation around a steady-state solution,

(A3)$$\matrix{H = H_0 + H_1 \cr U = U_0 + U_1\cr\; \quad\quad\quad\quad \dot\varepsilon = {\partial U_0\over \partial x} + {\partial U_1\over \partial x} = \dot\varepsilon_0 + \dot\varepsilon_1 .}$$

Subscripts 0 and 1 refer to the datum state and perturbation, respectively.

We start with the momentum balance equation. Plugging Eqn (A3) into Eqn (A1) yields

(A4)$$\matrix{& \Omega_a{\partial\over \partial x}\left[( H_0 + H_1) ( \dot\varepsilon_0 + \dot\varepsilon_1) ^{1/3}\right]\cr & \quad- \left[\beta\left(H_0 + H_1-{\rho_{\rm w}\over \rho_i}D\right) + \Omega_{b}( H_0 + H_1) \right]( U_0 + U_1) ^{1/3} \cr & \quad = ( H_0 + H_1) \left[{\partial H_0\over \partial x} + {\partial H_1\over \partial x}\right]. }$$

According to the binomial approximation,

(A5)$$( \dot\varepsilon_0 + \dot\varepsilon_1) ^{1/3} = \dot\varepsilon_0^{1/3}\left(1 + {\dot\varepsilon_1\over \dot\varepsilon_0}\right)^{1/3}\approx \dot\varepsilon_0^{1/3}\left(1 + {1\over 3}{\dot\varepsilon_1\over \dot\varepsilon_0}\right) = \dot\varepsilon_0^{1/3} + {1\over 3}\dot\varepsilon_0^{-2/3}\dot\varepsilon_1,\; $$

and similarly for (U 0 + U 1) 1/3. Thus, Eqn (A4) can be written as

(A6)$$\matrix{& \Omega_a{\partial\over \partial x}\left[( H_0 + H_1) \left(\dot\varepsilon_0^{1/3} + {1\over 3}\dot\varepsilon_0^{-2/3}\dot\varepsilon_1\right)\right]\;-\cr & \quad\left[\beta\left(H_0 + H_1-{\rho_{\rm w}\over \rho_i}D\right) + \Omega_{b}( H_0 + H_1) \right]\left(U_0^{1/3} + {1\over 3}U_0^{-2/3}U_1\right) \cr & \qquad = ( H_0 + H_1) \left[{\partial H_0\over \partial x} + {\partial H_1\over \partial x}\right]. }$$

Expanding, eliminating the steady-state solution and dropping higher order perturbation terms yields

(A7)$$\matrix{& \Omega_a\left[{1\over 3}\dot\varepsilon_0^{-2/3}\left({\partial H_0\over \partial x}-{2\over 3}{H_0\over \dot\varepsilon_0}{\partial\dot\varepsilon_0\over \partial x} \right)\dot\varepsilon_1 + {H_0\over 3}\dot\varepsilon_0^{-2/3}{\partial \dot\varepsilon_1\over \partial x} + \dot\varepsilon_0^{1/3}{\partial H_1\over \partial x} + {1\over 3}\dot\varepsilon_0^{-2/3}{\partial\dot\varepsilon_0\over \partial x}H_1\right]\cr & \quad-( \beta + \Omega_{b}) U_0^{1/3}H_1 - \left[\beta\left(1-{\rho_{\rm w}\over \rho_i}{D\over H_0}\right) + \Omega_{b}\right]{H_0\over 3}U_0^{-2/3}U_1 \cr & \quad = H_0{\partial H_1\over \partial x} + H_1{\partial H_0\over \partial x}. }$$

Next we rearrange and group the perturbation terms, yielding

(A8)$$\matrix{\left\{\left[\beta\left(1-{\rho_{\rm w}\over \rho_i}{D\over H_0}\right) + \Omega_{b}\right]{H_0\over 3}U_0^{-2/3} \right\}U_1 = \left\{{1\over 3}\Omega_a\dot\varepsilon_0^{-2/3}{\partial\dot\varepsilon_0\over \partial x}-( \beta + \Omega_{b}) U_0^{1/3}-{\partial H_0\over \partial x}\right\}H_1 & \cr \quad + \left\{\Omega_a\dot\varepsilon_0^{1/3}-H_0 \right\}{\partial H_1\over \partial x} + \left\{{1\over 3}\Omega_a\dot\varepsilon_0^{-2/3}\left({\partial H_0\over \partial x}-{2\over 3}{H_0\over \dot\varepsilon_0}{\partial \dot\varepsilon_0\over \partial x} \right)\right\}\dot\varepsilon_1 & \cr\quad + \left\{\Omega_a{H_0\over 3}\dot\varepsilon_0^{-2/3} \right\}{\partial \dot\varepsilon_1\over \partial x}. }$$

Here and hereafter we use curly brackets to indicate functions that depend only on the datum state and the forcing frequency. To simplify the algebra later on, we label these terms in Eqn (A8) as C 1–C 5:

(A9)$$C_1U_1 = C_2H_1 + C_3{\partial H_1\over \partial x} + C_4\dot\varepsilon_1 + C_5{\partial \dot\varepsilon_1\over \partial x}.$$

Similarly, perturbation of the mass continuity equation (Eqn (A2)) yields

(A10)$${\partial\over \partial t}( H_0 + H_1) = \dot{B} - {\partial\over \partial x}\left[( H_0 + H_1) ( U_0 + U_1) \right].$$

After expanding, eliminating the steady-state solution and dropping higher order terms, this becomes

(A11)$${\partial H_1\over \partial t} = -\left\{{\partial H_0\over \partial x}\right\}U_1 -\{ H_0\} \dot\varepsilon_1 -\left\{\dot\varepsilon_0\right\}H_1 -\{ U_0\} {\partial H_1\over \partial x}.$$

Equations (A8) and (A11) are the SSA equivalent of the well-known kinematic wave equation. The classic kinematic wave equation is derived by invoking the shallow ice approximation, which neglects membrane stresses (Nye, Reference Nye1960; Pfeffer, Reference Pfeffer2007; Felikson and others, Reference Felikson2017). This approach also differs from that of Walters (Reference Walters1989), who accounted for membrane stresses when analyzing the effect of tidal variations on tidewater glacier flow but did not allow for perturbations in viscosity or ice thickness and did not investigate the frequency response.

We next express the perturbations in terms of periodic functions, and in doing so we assume that the perturbations in thickness and velocity (and strain rate) have the same wavenumber and frequency. The perturbations are applied at the terminus (x = 0), and thus

(A12)$$\matrix{H_1 \hskip-8pt& = H^\star {\rm e}^{i( kx + \omega t) }\cr U_1 \hskip-9pt& = U^\star {\rm e}^{i( kx + \omega t) }\cr \dot\varepsilon_1 \hskip-10pt& \! \! \! \! \! \! \! \! \! \! \!\! \! = ikU_1,\; }$$

where $^\star$ refers to the amplitude of the perturbation, ω is the frequency of the forcing at the terminus and k is the complex spatial wavenumber.

By inserting Eqn (A12) into Eqn (A11) we find that the perturbations are related to each other by the wavenumber and frequency:

(A13)$$H^\star = \left\{{-H_0k + i( {\partial H_0}/{\partial x}) \over \omega + U_0k-i\dot\varepsilon_0}\right\}U^\star.$$

By expanding Eqn (A12), it is also apparent that the wavelength λ, decay or e-folding length D L (defined here as being positive in the upstream direction) and phase speed v p are given by

(A14)$$\eqalign{\lambda & = {2\pi\over {\rm Re}( k) }\cr D_L & = -{1\over {\rm Im}( k) }\cr v_{\rm p} & = {\omega\over {\rm Re}( k) }. }$$

Since ice flow is defined as being in the positive x-direction, Re(k) > 0 and Im(k) < 0 imply propagation and decay in the upstream direction, respectively.

Inserting Eqns (A12)–(A13) into Eqn (A9) results in a cubic equation for the wave number:

(A15)$$\matrix{& \left\{C_5U_0\right\}k^3 + \left\{i( C_3H_0-C_4U_0-C_5\dot\varepsilon_0) + C_5\omega\right\}k^2 \cr & \quad + \left\{C_1U_0 + C_2H_0 + C_3{\partial H_0\over \partial x}-C_4\dot\varepsilon_0-iC_4\omega\right\}k \cr & \quad+ \left\{C_1\omega-i\left(C_1\dot\varepsilon_0 + C_2{\partial H_0\over \partial x}\right)\right\} = 0. }$$

The coefficients in Eqn (A15) vary spatially, indicating that the wavenumber also varies spatially. Using the datum state, Eqn (A15) could be used to determine k = k(x) and completely describe a wave as it propagates through a glacier. However, our goal is to develop a basic understanding of the response of tidewater glaciers to calving events. We therefore select characteristic values for the thickness, velocity and strain rate by using their respective values at the terminus. This is analogous to invoking the quasi-uniform approximation (e.g. Hindmarsh, Reference Hindmarsh2004), which allows for gradients in velocity to affect the viscosity but not the mean flow. We additionally assume that the thickness and strain rate gradients are small, i.e. ∂H 0/∂x ≈ 0 and $\partial \dot \varepsilon _0/\partial x\approx 0$, such that the terms involving ∂H 0/∂x and $\partial \dot \varepsilon _0/\partial x$ have relatively little impact on wave propagation. This latter assumption is supported by tests that we conducted in which we selected characteristic values of ∂H 0/∂x and $\partial \dot \varepsilon _0/\partial x$ and found only modest changes to the calculated wavenumbers.

Under these assumptions, Eqn (A15) reduces to

(A16)$$\matrix{& \left\{\Omega_a \dot\varepsilon_0^{-2/3}U_0\right\}k^3 + \left\{\omega\Omega_a \dot\varepsilon_0^{-2/3} + i\left(2\Omega_a\dot\varepsilon_0^{1/3}-3H_0\right)\right\}k^2 \cr & \quad + \left\{\left(\beta\left(1-{\rho_{\rm w}\over \rho_i}{D\over H_0}\right) + \Omega_{b}\right)U_0^{1/3} - 3( \beta + \Omega_{b}) U_0^{1/3}\right\}k \cr & + \left\{\left(\beta\left(1-{\rho_{\rm w}\over \rho_i}{D\over H_0}\right) + \Omega_{b}\right)U_0^{-2/3}( \omega-i\dot\varepsilon_0) \right\} = 0. }$$

The wavenumber k is determined by finding the roots of Eqn (A16) for a given glacier geometry, glacier flow and forcing frequency. Two of the three roots produce similar decay lengths and wavelengths but with opposite signs, representing waves that propagate in the upstream and downstream directions. Note that the perturbation analysis does not specify where ice exists, only that the perturbation occurs at x = 0. Of these two, we select the solution for which Im(k) < 0 and Re(k) > 0 and therefore the wave propagates upstream and decays in the upstream direction. The third root always has short-decay lengths and wavelengths, ${\cal O}( 100\hbox { m})$, and we suspect it is a result of the quasi-uniform approximation (i.e. using a constant velocity and nonzero strain rate). The root that we select produces decay lengths, phase velocities and wavelengths that are consistent with the SSA.

The high-frequency limit can be found by analyzing Eqn (A16) for large ω, which leaves

(A17)$$\left\{\Omega_a\dot\varepsilon_0^{-2/3}\right\}k^2 + \left\{\left(\beta\left(1-{\rho_{\rm w}\over \rho_i}{D\over H_0}\right) + \Omega_{b}\right)U_0^{-2/3}\right\} = 0.$$

Solving for the wave number yields

(A18)$$k = - i\left({\dot\varepsilon_0\over U_0}\right)^{1/3} \sqrt{\left(\beta\left(1-{\rho_{\rm w}\over \rho_i}{D\over H_0}\right) + \Omega_{b}\right)\Omega_a^{-1}},\; $$

where we have selected the negative root to ensure that perturbations decay in the upstream direction. Consequently, the wavelength and phase velocity diverge in the high-frequency limit and the decay length converges:

(A19)$$D_L \rightarrow \left({U_0\over \dot\varepsilon_0}\right)^{1/3} \sqrt{{2\over A^{1/3}\beta\rho_ig( 1-( {\rho_{\rm w}}/{\rho_i}) ( {D}/{H_0}) ) + ( {1}/{W}) ( {4}/{W}) ^{1/3}}}.$$

References

Amundson, JM and 5 others (2008) Glacier, fjord, and seismic response to recent large calving events, Jakobshavn Isbræ, Greenland. Geophysical Research Letters 35(22), L22501. doi: 10.1029/2008GL035281.CrossRefGoogle Scholar
Bassis, JN and Walker, CC (2012) Upper and lower limits on the stability of calving glaciers from the yield strength envelope of ice. Proceedings of the Royal Society A: Mathematical, Physical 468(2140), 913–931. doi: 10.1098/rspa.2011.0422.CrossRefGoogle Scholar
Budd, W, Keage, P and Blundy, N (1979) Empirical studies of ice sliding. Journal of Glaciology 23(89), 157–170. doi: 10.3189/S0022143000029804.CrossRefGoogle Scholar
Cassotto, R, and 6 others (2019) Nonlinear glacier response to calving events, Jakobshavn Isbræ, Greenland. Journal of Glaciology 65(249), 39–54. doi: 10.1017/jog.2018.90.CrossRefGoogle Scholar
Catania, G, Stearns, L, Moon, T, Enderlin, E and Jackson, R (2020) Future evolution of Greenland's marine-terminating outlet glaciers. Journal of Geophysical Research. Earth Surface 125, e2018JF004873. doi: 10.1029/2018JF004873.CrossRefGoogle Scholar
Cook, S and 7 others (2014) Modelling environmental influences on calving at Helheim Glacier in eastern Greenland. The Cryosphere 8, 827–841. doi: 10.5194/tc-8-827-2014.CrossRefGoogle Scholar
Cowton, TR, Todd, JA and Benn, DI (2019) Sensitivity of tidewater glaciers to submarine melting governed by plume locations. Geophysical Research Letters 46(20), 11219–11227. doi: 10.1029/2019GL084215.CrossRefGoogle Scholar
Cuffey, KM and Paterson, WSB (2010) The physics of glaciers. 4th ed., Elsevier, Amsterdam.Google Scholar
de Juan, J and 12 others (2010) Sudden increase in tidal response linked to calving and acceleration at a large Greenland outlet glacier. Geophysical Research Letters 37(12), L12501. doi: 10.1029/2010GL043289.CrossRefGoogle Scholar
Enderlin, E, Howat, I and Vieli, A (2013) High sensitivity of tidewater outlet glacier dynamics to shape. The Cryosphere 7, 1007–1015. doi: 10.5194/tc-7-1007-2013.CrossRefGoogle Scholar
Enderlin, EM, O'Neel, S, Bartholomaus, TC and Joughin, I (2018) Evolving environmental and geometric controls on Columbia Glacier's continued retreat. Journal of Geophysical Research. Earth Surface 123, 1528–1545. doi: 10.1029/2017JF004541.CrossRefGoogle Scholar
Favier, L, Gagliardini, O, Durand, G and Zwinger, T (2012) A three-dimensional full Stokes model of the grounding line dynamics: effect of a pinning point beneath the ice shelf. The Cryosphere 6(1), 101–112. doi: 10.5194/tc-6-101-2012.CrossRefGoogle Scholar
Felikson, D and 11 others (2017) Inland thinning on the Greenland ice sheet controlled by outlet glacier geometry. Nature Geoscience 10, 366–369. doi: 10.1038/ngeo2934.CrossRefGoogle Scholar
Gagliardini, O, Durand, G, Zwinger, T, Hindmarsh, RCA and Le Meur, E (2010) Coupling of iceshelf melting and buttressing is a key process in icesheets dynamics. Geophysical Research Letters 37(14), L14501. doi: 10.1029/2010GL043334.CrossRefGoogle Scholar
Gladstone, RM and 5 others (2017) Marine ice sheet model performance depends on basal sliding physics and sub-shelf melting. The Cryosphere 11, 319–329. doi: 10.5194/tc-11-319-2017.CrossRefGoogle Scholar
Hill, EA, Gudmundsson, GH, Carr, JR and Stokes, CR (2018) Velocity response of Petermann Glacier, northwest Greenland, to past and future calving events. The Cryosphere 12(12), 3907–3921. doi: 10.5194/tc-12-3907-2018.CrossRefGoogle Scholar
Hindmarsh, RCA (2004) A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling. Journal of Geophysical Research 109, F01012. doi: 10.1029/2003JF000065.CrossRefGoogle Scholar
Hindmarsh, RCA (2012) An observationally validated theory of viscous flow dynamics at the ice-shelf calving front. Journal of Glaciology 58(208), 375–387. doi: 10.3189/2012JoG11J206.CrossRefGoogle Scholar
Jenkins, A (2011) Convection-driven melting near the grounding lines of ice shelves and tidewater glaciers. Journal of Physical Oceanography 41, 2279–2294. doi: 10.1175/JPO-D-11-03.1.CrossRefGoogle Scholar
Kane, E and 6 others (2020) Impact of calving dynamics on Kangilernata Sermia, Greenland. Geophysical Research Letters 47(20), e2020GL088524. doi: 10.1029/2020GL088524.CrossRefGoogle Scholar
Krug, J, Durand, G, Gagliardini, O and Weiss, J (2015) Modelling the impact of submarine frontal melting and ice mélange on glacier dynamics. The Cryosphere 9, 989–1003. doi: 10.5194/tc-9-989-2015.CrossRefGoogle Scholar
Krug, J, Weiss, J, Gagliardini, O and Durand, G (2014) Combining damage and fracture mechanics to model calving. The Cryosphere 8(6), 2101–2117. doi: 10.5194/tc-8-2101-2014.CrossRefGoogle Scholar
Ma, Y and Bassis, JN (2019) The effect of submarine melting on calving from marine terminating glaciers. Journal of Geophysical Research. Earth Surface 124, 334–346. doi: 10.1029/2018JF004820.CrossRefGoogle Scholar
MacAyeal, DR (1989) Large-scale ice flow over a viscous basal sediment: theory and application to Ice Stream B, Antarctica. Journal of Geophysical Research 94, 4071–4087. doi: 10.1029/JB094iB04p04071.CrossRefGoogle Scholar
Mercenier, R, Lüthi, MP and Vieli, A (2018) Calving relation for tidewater glaciers based on detailed stress field analysis. The Cryosphere 12, 721–739. doi: 10.5194/tc-12-721-2018.CrossRefGoogle Scholar
Mercenier, R, Lüthi, MP and Vieli, A (2019) A transient coupled ice flow-damage model to simulate iceberg calving from tidewater outlet glaciers. Journal of Advances in Modeling Earth Systems 11, 3057–3072. doi: 10.1029/2018MS001567.CrossRefGoogle Scholar
Motyka, RJ, Dryer, WP, Amundson, J, Truffer, M and Fahnestock, M (2013) Rapid submarine melting driven by subglacial discharge, LeConte Glacier, Alaska. Geophysical Research Letters 40(19), 5153–5158. doi: 10.1002/grl.51011.CrossRefGoogle Scholar
Nettles, M and 12 others (2008) Step-wise changes in glacier flow speed coincide with calving and glacial earthquakes at Helheim Glacier, Greenland. Geophysical Research Letters 35(24), L24503. doi: 10.1029/2008GL036127.CrossRefGoogle Scholar
Nye, J (1960) The response of glaciers and ice-sheets to seasonal and climatic changes. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 256, 559–584. doi: 10.1098/rspa.1960.0127.Google Scholar
O'Leary, M and Christoffersen, P (2013) Calving on tidewater glaciers amplified by submarine frontal melting. The Cryosphere 7(1), 119–128. doi: 10.5194/tc-7-119-2013.CrossRefGoogle Scholar
Pfeffer, WT (2007) A simple mechanism for irreversible tidewater glacier retreat. Journal of Geophysical Research 112, F03S25. doi: 10.1029/2006JF000590CrossRefGoogle Scholar
Riel, B, Minchew, B and Joughin, I (2021) Observing traveling waves in glaciers with remote sensing: new flexible time series methods and application to Sermeq Kujalleq (Jakobshavn Isbræ), Greenland. The Cryosphere 15(1), 407–429. doi: 10.5194/tc-15-407-2021.CrossRefGoogle Scholar
Rosenau, R, Schwalbe, E, Maas, HG, Baessler, M and Dietrich, R (2013) Grounding line migration and high-resolution calving dynamics of Jakobshavn Isbræ, West Greenland. Journal of Geophysical Research. Earth Surface 118, 382–395. doi: 10.1029/2012JF002515.CrossRefGoogle Scholar
Schlemm, T and Levermann, A (2019) A simple stress-based cliff-calving law. The Cryosphere 13, 2475–2488. doi: 10.5194/tc-13-2475-2019.CrossRefGoogle Scholar
Sergienko, O (2010) Elastic response of floating glacier ice to impact of long-period ocean waves. Journal of Geophysical Research 115, F04028. doi: 10.1029/2010JF001721.CrossRefGoogle Scholar
Todd, J and Christoffersen, P (2014) Are seasonal calving dynamics forced by buttressing from ice mélange or undercutting by melting? Outcomes from full-Stokes simulations of Store Glacier, West Greenland. The Cryosphere 8, 2353–2365. doi: 10.5194/tc-8-2353-2014.CrossRefGoogle Scholar
Vallot, D and 9 others (2018) Effects of undercutting and sliding on calving: a global approach applied to Kronebreen, Svalbard. The Cryosphere 12(2), 609–625. doi: 10.5194/tc-12-609-2018.CrossRefGoogle Scholar
van der Veen, CJ (2013) Fundamentals of glacier dynamics. 2 ed., CRC Press, Boca Raton, FL.CrossRefGoogle Scholar
Walters, R (1989) Small-amplitude, short-period variations in the speed of a tide-water glacier in South-Central Alaska, U.S.A. Annals of Glaciology 12, 187–191. doi: 10.3189/S0260305500007175.CrossRefGoogle Scholar
Williams, CR, Hindmarsh, RCA and Arthern, RJ (2012) Frequency response of ice streams. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 468, 3285–3310. doi: 10.1098/rspa.2012.0180.Google ScholarPubMed
Figure 0

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).

Figure 1

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.

Figure 2

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.

Figure 3

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).

Figure 4

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.

Figure 5

Fig. 6. Dependence of the high-frequency decay length on the proximity to flotation (D/H0), 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.