## 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

and

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:

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

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:

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:

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, 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:

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

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

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)

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

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:

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.

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:

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

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.

## 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*:

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

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

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.

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

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

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

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/3}*U* = *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

and

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

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

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

According to the binomial approximation,

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

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

Next we rearrange and group the perturbation terms, yielding

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}:

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

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

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

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:

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

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:

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

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

Solving for the wave number yields

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: