Beyond the Stokes approximation: shallow visco-elastic ice-sheet models

Abstract

The hypothesis that ice-sheet evolution is only controlled by the long-term non-Newtonian viscous behavior of ice has been challenged by observations indicating that effects like brittle failure, stick-slip sliding, tides and wave action may affect ice-sheet evolution on sub-daily timescales. Over these timescales, the quasi-static-creep approximation is no longer appropriate and elastic effects become important. Simulating elastic effects in ice-sheet models over relevant timescales, however, remains challenging. Here, we show that by including a visco-elastic rheology and reintroducing the oft neglected acceleration term back into the ice-sheet stress balance, we can create a visco-elastic system where the velocity is locally determined and information propagates at the elastic wave speed. Crucially, the elastic wave speed can be treated like an adjustable parameter and set to any value to reproduce a range of phenomena, provided the wave speed is large compared to the viscous velocity. We illustrate the system using three examples. The first two examples demonstrate that the system converges to the steady-state viscous and elastic limits. The third example examines ice-shelf rifting and iceberg calving. This final example hints at the utility of the visco-elastic formulation in treating both long-term evolution and short-term environmental effects.

1. Introduction

Ice sheets have traditionally been viewed as thin films of non-Newtonian fluids, deforming slowly over centuries to millennia in response to gravitationally induced pressure gradients (Cuffey and Paterson, 2010). It is only recently that this view has been replaced with the understanding that ice sheets can undergo rapid changes over daily and shorter time scales. This has been pointedly illustrated by the unexpectedly rapid collapse of the Larsen B ice shelf in Antarctica and the subsequent acceleration of its tributary glaciers (Rignot and others, 2004; Scambos and others, 2004). Additional examples include the twice daily tidally driven ‘stick-slip’ surges of Whillans Ice Stream (e.g. Bindschadler and others, 2003; Winberry and others, 2011), the rapid step-like increase in ice velocity following calving events from the ice front of Helheim glacier (Nettles and others, 2008; Roeoesli and others, 2016) and the sudden acceleration of portions of ice sheets following rapid supra-glacial lake drainage (e.g. Das and others, 2008; Tedesco and others, 2013).

The ice-sheet process most closely associated with short-term elastic effects, however, is iceberg calving. Iceberg calving involves processes that range from rapid fracture propagation over time scales of minutes (or shorter) all the way to rifting and the decadal (or longer) detachment of tabular icebergs from ice shelves (e.g. Fricker and others, 2002; Benn and others, 2007; Bassis and others, 2008; Walker and others, 2013; Jeong and others, 2016). The challenge in accurately simulating calving and the response of ice sheets to environmental forcing over a range of time scales is all the more daunting given increased evidence that short-term atmospheric and oceanic drivers, like atmospheric rivers and ocean swell, play a role in triggering iceberg detachment (e.g. Massom and others, 2018; Francis and others, 2021; Lipovsky, 2022; Wille and others, 2022). These results hint at the need to simultaneously resolve both the glaciological effect of short-term processes on ice-sheet dynamics and the mechanical coupling between ice sheets and the atmosphere/ocean over time scales as short as minutes to hours.

The dilemma ice-sheet modelers face is that to inspire confidence in their ability to forecast and project future ice-sheet changes, numerical models must be able to reproduce past ice-sheet variability over time scales of centuries or longer. At the same time, models must be able to explain observations that show ice sheets and glaciers exhibit variations on diurnal and shorter time scales where elastic effects can be important. This points to a simultaneous need for higher temporal resolution and sufficient efficiency to perform long time scale simulations. At present, these two goals are accomplished using very different modeling approaches. For example, discrete element models simulate the short-term brittle failure of ice on time scales of seconds to minutes (Bassis and Jacobs, 2013; Benn and others, 2017; Åström and Benn, 2019; Benn and others, 2022). These shorter-term effects then enter as semi-heuristic parameterizations in today's larger-scale viscous ice-sheet models (e.g. DeConto and Pollard, 2016; Crawford and others, 2021). On tidal timescales, there is a long history of applying visco-elastic models, especially in a flowline context, to understand grounding line flexure and tidally modulated flow of ice streams (e.g. Reeh and others, 2003; Walker and others, 2012; Rosier and others, 2014, 2017; Wild and others, 2017). These models provide key insight into shorter timescale dynamics, but they are often too computationally expensive to integrate into larger-scale ice-sheet–ice-shelf studies.

Here we present a method which unifies the two seemingly disparate goals of accurate short time scale and long time scale behavior within a single modeling framework. To this end, we introduce an elastic component to the rheology of ice while also re-introducing the often neglected acceleration term into the stress balance. These subtle changes alter the system of equations from a set of elliptic equations, where information propagates instantaneously, to a hyperbolic system. This new system is analogous to the hyperbolic heat equations (Cattaneo, 1958) and, like the hyperbolic heat equation, a consequence of this change is that the velocity is locally determined by the wave speed. However, so long as we are not interested in wave propagation, by judiciously choosing the elastic wave speed we can use the same model to simulate phenomena across a range of time scales where the march to steady state is analogous to the iterations to convergence associated with solving large systems of non-linear equations in traditional ice-sheet models. Intriguingly, the hyperbolic visco-elastic system that we propose has analogies with the regularized elasto-visco-plastic rheology conventionally used for sea ice (Hunke and Dukowicz, 1997), although inertial terms play a larger role in sea-ice dynamics and sea-ice models use a very different constitutive relation.

To illustrate the approach, we start by considering a simple, prototypical problem: that of a one-dimensional ice shelf or ice stream. We first review the equations of conservation of momentum and the rheology governing the viscous deformation of a freely floating ice shelf in 1D (the two-dimensional case is treated in Appendix A). We then show that a change in the constitutive relationship to include elasticity transforms the system into a fully hyperbolic system of equations for deviatoric stress and velocity. We then analyze the various regimes of flow, showing that we can recover elastic, visco-elastic and viscous limits. Finally, we present a series of numerical examples illustrating the proposed method's ability to resolve both short-term elastic displacements and longer-term viscous flow. Our final example illustrates the flexibility of the approach by simulating rifting of ice shelves and iceberg detachment. This example shows that the hyperbolic system can be used to simulate both ice-shelf dynamics and iceberg drift, hinting that it may be possible to couple ice sheets into Earth System Models by considering ice sheets as part of the oceans and co-evolving the ice sheet and ocean system together. Readers uninterested in the mathematical or historical exposition may wish to skip directly to the Numerical Examples section.

2. Governing equations for an ice shelf

2.1. Conservation of momentum

We begin by considering the governing equations for a one-dimensional ice shelf. We focus on the one-dimensional case here to simplify the exposition, but include the fully two-dimensional case in Appendix A. Denoting the depth averaged deviatoric stress by τ _xx and the thickness of the ice shelf by h(x,  t), conservation of momentum can be written as follows (e.g. MacAyeal, 1989; Weis and others, 1999):

(1)$${\rm {{\it D}\over {\it Dt}}}\left(\rho h u\right) + {\partial\over \partial x}\left(- 2\, h \tau_{xx} + {1\over 2} \rho g' h^2\right) = F$$

In the above equations D/Dt denotes the material derivative, u(x,  t) is the horizontal velocity, ρ is the (assumed constant) density of ice and F(u,  x,  t) is an additional forcing (e.g. friction from the ice bed, lateral drag from embayment walls, ocean drag, Coriolis force, etc.). We have also defined the ‘reduced’ acceleration due to gravity,

(2)$$g' = \left(1 - {\rho\over \rho_{\rm w}}\right)g,\; $$

with ρ _w the density of water and g = 9.81 m/s² the acceleration due to gravity. Unlike previous studies, we have explicitly retained the material derivative. The above approximation relies on two assumptions: (1) the ice is assumed to be shallow, i.e. the ratio of a characteristic thickness to relevant horizontal length scales is small compared to unity and; (2) vertical shear along the ice shelf base is assumed to be sufficiently small that horizontal velocity is nearly independent of depth. Where the ice shelf terminates in a cliff in the ocean, we require continuity with ocean pressure, where for a freely floating ice shelf τ _xx = ρg′h/4 at the calving front.

Because ice is a very viscous fluid, the material time derivative in Eqn (1) is typically neglected leading to the ‘shallow shelf approximation’ (SSA) (MacAyeal, 1989; Weis and others, 1999).

(3)$${\partial\over \partial x}\left(- 2\, h \tau_{xx} + {1\over 2} \rho g' h^2\right) = F,\; $$

where F could represent drag from the walls or bed if the ice is grounded. By contrast, when τ _xx is small and F is the Coriolis force, Eqn (1) is analogous to the shallow water equations, suggesting a connection between ocean and ice-shelf dynamics that we will eventually return to in the Discussion (Pedlosky, 1987).

2.2. Viscous constitutive relationship

Over long time scales, deviatoric stresses are often related to velocities using a power-law quasi-viscous rheology of the form:

(4)$$\tau_{xx} = 2\eta {\partial u\over \partial x}$$

where η is the effective viscosity, given by

(5)$$\eta = {B\over 2} \left\vert {\partial u\over \partial x}\right\vert ^{{1}/{n}-1},\; $$

with flow law exponent n ~ 3 (Cuffey and Paterson, 2010). The rate parameter B is a function of the depth-averaged temperature. For simplicity and because it does not impact the treatment that follows, we assume B is constant and neglect the possible thermo-mechanical variations of B in our analysis. It is also possible to write the rheology in terms of the effective deviatoric stress:

(6)$$\eta = {B^n\over 2} \left\vert \tau_{xx}\right\vert ^{1-n}.$$

This later form of the rheology will be more convenient when we consider the strain and strain rate associated with both viscous and elastic deformation.

2.3. The shallow shelf approximation

Substituting Eqn (4) into Eqn (3) we arrive at the SSA:

(7)$${\partial\over \partial x}\left(- 2\, h \eta {\partial u\over \partial x} + {1\over 2} \rho g' h^2\right) = F.$$

The SSA is a non-linear elliptic equation for the velocity. Although the ice-shelf model is simpler than higher-order or full Stokes models, the non-local nature of the velocity is shared with its more complicated brethren. The key feature of Eqn (7) is that, as a consequence of dropping the material derivative, information propagates throughout the system infinitely fast and we now need to solve systems of non-linear equations to determine the velocity field. Numerically, discretizing the elliptic system results in systems of non-linear equations that need to be solved using iterative methods. Few modelers have escaped quietly weeping in despair when applying an existing model to a new situation and witnessing the solver fail to converge. Moreover, the elliptic system requires a combination of deviatoric stress and velocity boundary conditions along some portion of the domain to render the problem well posed. For a freely floating iceberg with stress boundaries and no inflow velocity specified, the problem is actually ill-posed unless additional conditions are imposed to eliminate arbitrary rigid body translation and rotation (Zarrinderakht and others, 2022).

2.4. Momentum diffusion, the parabolic limit

To avoid the problem of a non-local velocity field, it is tempting to simply re-introduce the time derivative of velocity back into Eqn (7) and solve the time dependent parabolic problem:

(8)$$\dot u + {1\over \rho h}{\partial\over \partial x} \left(-2\, h \eta {\partial u\over \partial x} + {1\over 2} \rho g' h^2\right) = 0.$$

Here we have retained the time derivative, denoted with the dot decoration, but have dropped the advective part of the material derivative. This is essentially the strategy used by Berg and Bassis (2022), although applied to a full Stokes system. It is, however, illustrative to estimate the time scale of momentum diffusion implied by Eqn (8). Dimensional analysis shows that the diffusion time scale is approximately given by:

(9)$$t_{\rm D} = {L^2\rho\over \eta_0},\; $$

where L is a characteristic horizontal dimension and η ₀ a reference characteristic viscosity. For L=100 km, ρ=910 kg m⁻³ and η ₀ = 10¹² − 10¹⁵ Pa s we find t _D < 10 s. This should be contrasted with the 50 s it takes a seismic p-wave with a velocity 2000 m s⁻¹ to traverse the domain. Compressional waves are amongst the fastest of several classes of seismic waves with the slower modes providing an even larger contrast in characteristic time scales between elastic wave propagation and momentum diffusion.

Moreover, information in the parabolic limit still propagates at infinite speed (Cattaneo, 1958). Of course, if we are only interested in the steady state, we can scale the time derivatives in various ways and simply march Eqn (8) forward in time until it reaches steady state. Variations of this trick of re-introducing acceleration terms are frequently used as numerical methods that apply various mass weighting and other parameterization to accelerate convergence of the integration in ‘pseudo time’ (Logan and others, 2017; Räss and others, 2020). However, as we show next, we can transform the (parabolic) diffusion equation into a hyperbolic system of equations that limits information to a characteristic wave speed of the system by including a visco-elastic rheology, allowing us to resolve processes at a greater range of time scales.

3. A visco-elastic ice-shelf model

3.1. Maxwell visco-elasticity

One model for a material that behaves elastically on short time scales and viscously on long times scales is the Maxwell visco-elastic rheology, which can be viewed as a viscous dashpot connected in series with an elastic spring. The Maxwell visco-elastic rheology can be written

(10)$${\rm {{\it D}\over {\it Dt}}} \left(\tau_{ik} - p\delta_{ik}\right) + {\mu\over \eta} \tau_{ik} = 2\mu \dot \varepsilon_{ik} + \lambda \dot \varepsilon_{\,jj}\delta_{ik},\; $$

where λ and μ are the moduli of elasticity and $\varepsilon _{ik}$ represent the components of the strain tensor and $\dot \varepsilon _{ik}$ the strain rate tensor. Einstein summation convention is assumed and indices (i,  j,  k) run over coordinates (x,  y,  z). The effective viscosity η is most conveniently defined in terms of the effective stress, as in Eqn (6), because the strain rate includes both viscous and elastic deformation. The material derivative D/Dt of the deviatoric stress takes the form:

(11)$${\rm {{\it D}\over {\it Dt}}} \tau_{ik} = \dot \tau_{ik} + u_j{\partial\tau_{ik}\over \partial_{x_j}} - W_{ij}\tau_{\,jk} + \tau_{ij}W_{\,jk}.$$

with W _ik the vorticity or spin tensor (Malvern, 1969) and the dot decoration represents time derivatives. The rather mysterious three extra terms that appear on the right-hand side of the above equation are necessary to ensure the time derivative of the stress deviator remain frame independent. In the interest of simplicity, we shall make the a priori assumption that rotational terms and advective terms are small so that we can replace all material derivatives with the time derivative and set:

$${\rm {{\it D}\tau_{ik}\over {\it Dt}}} \approx \dot \tau_{ik},\; \quad {Dp\over Dt}\approx \dot p,\; \quad {Du\over Dt}\approx \dot u.$$

This approximation is adequate for the small strain visco-elastic limit and for the slow, viscous limit, but will become inadequate in the large strain, visco-elastic limit. The approximation is not necessary and is made here primarily to simplify the exposition that follows.

Our task is to derive an expression for the Maxwell rheology of ice appropriate for use in our shallow shelf model. To this end, we neglect flexural stresses and assume for the one-dimensional flow geometry that τ _yy = τ _xy = 0. From Eqn (10) we conclude that $\dot \varepsilon _{xy}$ vanishes. Setting i = k in Eqn (10) and summing the three resulting equations yields an expression for the time evolution of pressure

(12)$$\dot p = - \kappa \left(\dot \varepsilon_{xx} + \dot \varepsilon_{yy} + \dot \varepsilon_{zz}\right),\; $$

where, again, the dot decoration denotes time derivatives and we have defined the bulk modulus of compression

(13)$$\kappa = \left(\lambda + {2\over 3}\mu\right).$$

Substituting Eqn (12) into Eqn (10), the longitudinal component of the Maxwell visco-elastic rheology can be expressed in the form

(14)$${\dot \tau_{xx}} + \left({\lambda\over \kappa}-1\right)\dot p + {\mu\over \eta} \tau_{xx} = 2\mu \dot \varepsilon_{xx}.$$

To make further progress, we need to determine an additional evolution equation for pressure. The simplest assumption corresponds to the incompressible limit, whence from Eqn (12) $\dot p\approx 0$. We consider a more general model that includes compressibility in Appendix A. In either case, we can eliminate $\dot p$ from Eqn (14) and express the Maxwell visco-elastic rheology solely in terms of deviatoric stress and horizontal strain rate

(15)$$\dot \tau_{xx} = {\left(2\eta \dot \varepsilon_{xx} - \tau_{xx}\right)\over T_{\rm r}}.$$

Here we have defined the Maxwell relaxation time

$$T_{\rm r} = {\eta\over G.}$$

and the effective shear modulus (see Appendix B):

(16)$$G = \left\{\matrix{ \mu \hfill & {\rm Incompressible}\hfill \cr \mu {3\lambda + 2\mu \over 3\lambda + 4\mu} \hfill & {\rm Compressible} \hfill }\right.$$

Upon depth-integrating Eqn (15), and then replacing all quantities with their depth-averaged value, including ρ, B and G, we find the system of equations (along with Eqn (1) with F=0):

(17)$$\dot u = {1\over \rho h}{\partial\over \partial x}\left(2\sigma_{xx} - \rho g' {h^2\over 2}\right)$$

(18)$$\dot \sigma_{xx} = {\left(2\eta h {\partial u \over \partial x} - \sigma_{xx} \right)\over T_{\rm r}},\; $$

where we have introduced the new variable σ _xx = hτ _xx, which corresponds to the depth-integrated deviatoric stress.

At first glance, adding a time-dependent term to the rheology appears to be catastrophic. To solve the visco-elastic system, we now need to time-step both Eqns (17) and (18) to steady state, and only once in steady state does it converge to the conventional SSA Eqns (3) and (4). However, as we show next, this change alters the structure of the differential equation in a fundamental way: it is now fully hyperbolic. The advantage of the hyperbolic system is that, not only does information propagate at finite speeds, but both σ _xx and u are entirely locally determined. Crucially, we shall see that the wave speed can act as a free parameter that is independent of the quasi-steady viscous velocity solution allowing us to make the wave speed as large (or small) as we wish to investigate a range of phenomena. If we are solely concerned with quasi-static ice-sheet velocities, the only requirement is that the ratio of the ice-sheet velocity to elastic wave speed – the Mach number – remain small.

3.2. Waves in the visco-elastic system

In the previous section, we asserted that the relaxation systems correspond to a hyperbolic relaxation system with elastic wave speed. We now show that this is true. Our approach here follows the usual approach in linear wave theory. Defining the vectors

(19)$${\bf u} = \left(\matrix{u \cr \sigma_{xx} \cr }\right),\; \quad {\bf q} = \left(\matrix{g' {\partial h\over \partial x} \cr {\sigma_{xx}\over T_{\rm r}} \cr }\right)$$

and the matrix

(20)$${\bf A} = \left(\matrix{0 & -{2\over \rho h} \cr\ - {2\eta h\over T_{\rm r}} & 0 }\right),\; $$

we can express Eqns (17) and (18) in matrix notation:

(21)$$\dot {{\bf u}} + {\bf A}{\partial {\bf u}\over \partial x} = -{\bf q}.$$

The eigenvalues of A are $c = \pm \sqrt {{4G\over \rho }}$, where c is the elastic wave speed. To see this, defining the matrix of right eigenvectors, along with its inverse

(22)$${\bf R} = \left(\matrix{1 & 1 \cr {\rho h c\over 2} & -{\rho h c\over 2} }\right),\; \quad {\bf R}^{-1} = {1\over 2} \left(\matrix{1 & {2\over \rho hc} \cr 1 & -{2\over \rho h c} }\right),\; $$

we can diagonalize the matrix A

(23)$${\bf \Lambda} = {\bf R}^{-1} {\bf A} {\bf R} = \left(\matrix{-c & 0 \cr 0 & c }\right).$$

Consistent with our previous assumptions, we hold ice thickness constant in time and pre-multiply Eqn (21) with R⁻¹:

(24)$$\left({\bf R}^{-1}{\dot{{\bf u}}}\right) + {\bf \Lambda}{\partial\over \partial x}\left({\bf R}^{-1}{\bf u}\right) = -{\bf R}^{-1}{\bf q} + {\bf \Lambda}{\partial {\bf R}^{-1}\over \partial x}{\bf u},\; $$

and then define

(25)$${\bf w} = {\bf R}^{-1}{\bf u},\; $$

and

(26)$${\bf p} = -{\bf R}^{-1}{\bf q} + {\bf \Lambda}{\partial {\bf R}^{-1}\over \partial x} {\bf u},\; $$

so that Eqn (24) can be written in form

(27)$$\dot{{\bf w}} + {\bf \Lambda} {\partial\over \partial x}{\bf w} = {\bf p}.$$

The left-hand side of Eqn (27) now represents a pair of advection equations for the quantities

$$w^{( 1) } = {1\over 2}\left(u + {\sigma_{xx}\over \rho h c}\right),\;$$

and

$$w^{( 2) } = {1\over 2}\left(u - {\sigma_{xx}\over \rho h c}\right).$$

The two quantities w ⁽¹⁾ and w ⁽²⁾ represent waves propagating leftward (in the negative x-direction) and rightward (in the positive x-direction) with velocities ±c superposed on the longer-term viscous state. These two advection equations are coupled together through the source term on the right-hand side. We could consider more elaborate situations that, e.g. do away with the hydrostatic approximation, permit the ice thickness to vary with waves and/or include the advective and rotational parts of the material derivatives to provide a more complete treatment of waves and the large deformation limit. However, even with these limitations, we show next that we can consider different limits of flow within the same model framework simply by considering different limits of dimensionless variables.

3.3. Non-dimensional groups, limiting behavior and the march to steady state

We define characteristic scales for ice-sheet thickness H, horizontal length L, time T and horizontal speed U such that

$$x = L \tilde x,\; \quad h = H \tilde h,\; \quad u = U \tilde u,\; \quad t = {L\over U} \tilde t.$$

We scale the deviatoric stress τ _xx with unit ρg′ H/4. With these definitions, we can define a characteristic viscosity scale

(28)$$\eta_0 = {B^n\over 2}\left({\rho g H\over 4}\right)^{1-n}.$$

Defining the non-dimensional Froude number,

$${\rm Fr} = \sqrt{{2U^2\over {{g' H}}}},\;$$

the non-dimensional version of Eqns (17) and (18) thus become:

(29)$$\dot{ {\tilde u} } - {1\over \tilde h {\rm Fr}^2}{\partial\over \partial \tilde x}\left(\tilde \sigma_{xx} - \tilde h^2\right) = 0$$

(30)$$\dot{ \tilde{\sigma}}_{xx} = \left({{\rm Fr}\over {\rm M}}\right)^2\left(\tilde h {\partial \tilde u\over \partial \tilde x} - {\rm Ar} {\tilde \sigma_{xx}\over \tilde \eta} \right),\; $$

where we have introduced the Mach number ${\rm M} = U/c = U\sqrt {\rho /4G}$, and the Argand number,

$${\rm Ar} = \left({\rho g'H \over 4 B}\right)^n{L\over U}.$$

Equations (29) and (30) appear to be controlled by the set of three dimensionless parameters (Fr, M, Ar), which in turn depend on the larger set of material parameters and scales (ρ,  g,  η ₀,  G,  H,  L,  U). However, linearizing Eqns (29) and (30) about a steady-state condition:

(31)$$\tilde u( x,\; t) = U( x) + \delta u( x,\; t) $$

(32)$$\tilde{\sigma}_{xx}( x,\; t) = \Sigma( x) + S( x,\; t) ,\; $$

taking $\tilde h$ and $\tilde \eta$ constant for the linearized system and defining $Q = \tilde h\delta u$, ${\rm Ar}' = {\rm Ar}/\tilde \eta$ and introducing the Deborah number:

$${\rm De} = {\rm M}^2/( {\rm Ar}'{\rm Fr}^2) ,\;$$

the linearized system can be written:

(33)$${\rm M}^2\ddot Q = {\partial^2 Q\over \partial \tilde x^2} -{ {{\rm M}^2\over {\rm De}}}\dot Q$$

(34)$${\rm M}^2\ddot S = {\partial^2 S\over \partial \tilde x^2} -{{{\rm M}^2\over {\rm De}}}\dot S,\; .$$

This makes it clear that the behavior of the system is actually only controlled by two dimensionless parameters: the Mach number M and the Deborah number De, which characterizes the relative importance of the material's viscosity to elasticity. The Deborah number can also be expressed as the ratio of the Maxwell relaxation time to the characteristic time scale De = Tr/T, so that large De corresponds to elastic materials and small De corresponds to more viscous materials. We now see that the parabolic limit corresponds to the case when M² ≪ 1 while Ar′ ⋅ Fr² remains order unity. The hyperbolic limit corresponds to the elastic limit where De ≫ 1, or, equivalently, Ar′ ⋅ Fr² ≪ 1. The elliptic limit corresponds to M² ≪ 1 and De ≫ M² (or Ar′ ⋅ Fr² ≪ 1). We also see that the last terms on the right-hand side of Eqns (33)–(34) are proportional to De⁻¹ and look like dissipation terms, a fact that we demonstrate next. Values of the non-dimensional numbers depend on the time scale and velocity scale of interest, a fact that we exploit in the numerical examples to use the model to simulate both elastic and viscous behavior.

3.4. Dissipation and relaxation to steady state

To investigate how waves attenuate over time, we consider a slightly more general relaxation system. This generalized system includes the possibility of frictional dissipation if there is friction along the bed or sides of the ice. We write the more general equation for conservation of momentum (compare with Eqn (1)) as

(35)$$\tilde h \dot {\tilde u} - {1\over {\rm Fr}^2} {\partial\over \partial \tilde x}\left(\tilde{\sigma}_{xx} - \tilde h^2\right)- {2\over {\rm Fr}^2}\tilde h {\partial \tilde b\over \partial \tilde x} - \tilde \beta^2 \tilde u = 0$$

where $\tilde \beta ^2$ is a non-dimensional Newtonian friction coefficient and $\tilde b$ represents the (non-dimensional) position of the ice-sheet base relative to sea level. The $\tilde \beta ^2$ coefficient is assumed zero when the ice is floating. To understand why the term h∂b/∂x appears, note that if we define the (non-dimensional) position of the ice-sheet surface as $\tilde {s} = \tilde {b} + \tilde {h}$, then

(36)$$\tilde h {\partial \tilde h\over \partial \tilde x} + \tilde h {\partial \tilde b\over \partial \tilde x} = \tilde h{\partial \tilde s\over \partial \tilde x}.$$

If the ice is freely floating,

(37)$$\tilde s = \left(1 - {\rho_i\over \rho_w}\right)\tilde h,\; $$

from which we see that for an ice shelf

(38)$$\tilde h {\partial \tilde h\over \partial \tilde x} + \tilde h{\partial \tilde b\over \partial \tilde x} = {1\over 2}\left(1 - {\rho_{\rm i}\over \rho_{\rm w}}\right){\partial ( \tilde h^2) \over \partial \tilde x}.$$

This is the form we used in the original momentum balance for an ice shelf. In what follows we use the more general momentum balance valid for both grounded and floating ice, with the understanding that when the ice is freely floating $\tilde \beta ^2$ is zero and b is related to s via Archimedes principle.

Linearizing Eqns (29) and (30) about constant $\tilde u$ and η′ and defining $\beta = \tilde \beta /h$ yields the system of equations:

(39)$$\ddot{Q} = {1\over {\rm M}^2} {\partial^2 Q\over \partial \tilde x^2} - {1\over {\rm De}} \dot Q - {{\beta}^2\over {\rm De}}Q - \beta^2\dot Q$$

(40)$$\ddot{S} = {1\over {\rm M}^2} {\partial^2 S\over \partial \tilde x^2} - {1\over {\rm De}} \dot S - {\beta^2\over {\rm De}}S - {\beta}^2\dot S$$

Resolving Q and S into Fourier modes, the solutions

(41)$$Q = \hat Q \exp{\left(-\ell \tilde t + ik\tilde x\right)},\; \quad S = \hat S \exp{\left(-\ell \tilde t + ik\tilde x\right)},\; $$

satisfy Eqns (39) and (40) provided

(42)$$\ell = {1\over 2}\left({1\over {\rm De}} + { \beta}^2 \right)\pm {1\over 2}\sqrt{\left({1\over {\rm De}} + {\beta}^2\right)^2 - 4\left({k^2\over {\rm M}^2} + {{\beta}^2\over {\rm De}}\right)},\; $$

showing that the real part of ℓ is a damping factor with amplitude primarily controlled by (De⁻¹ + β ²). This makes it clear that both β and the inverse of De control damping of the system.

In the fully elastic limit De ≫ 1 and

(43)$$\ell = { \beta^2\over 2} \pm {1\over 2}\sqrt{ \beta^4 - 4k^2/{\rm M}^2},\; $$

which shows that attenuation is solely due to friction and ℓ is purely imaginary when β vanishes. We obtain the same behavior with large friction (β ≫ 1) and De of order unity. In the long wavelength limit where k ≪ βM/2, dissipation is independent of k and ℓ = β ².

In the viscous limit (De ≪ 1), with β order unity or smaller,

(44)$$\ell = {1\over 2{\rm De}}\left(1 \pm \sqrt{1 - 4{k^2\over {\rm M}^2}{\rm De}^2}\right),\; $$

which shows that attenuation in the viscous limit is controlled by the inverse of De and that in the long wavelength limit where k ≪ M/(2De), ℓ = 1/De and dissipation is again independent of wavelength. This analysis shows that in the viscous limit (De ≪ 1) or when bottom friction is large (β ² ≫ 1), waves will rapidly attenuate and the system will converge quickly to the steady-state viscous solution. If on the other hand friction is small and we are in the elastic limit (large De), waves will propagate across the entire domain (perhaps multiple times) before attenuating.

For a fixed geometry and fixed material parameters, the chosen velocity scale of interest determines the three dimensionless parameters. For simulations that target elastic and visco-elastic effects, we can set all dimensionless numbers to an appropriate value for ice and step the system forward in time with short time steps to get both the short time-scale transient visco-elastic effects, and eventually, steady-state creep. Alternatively, because the steady state of the relaxation system in Eqns (17) and (18) corresponds to the viscous Stokes flow problem of Eqns (3) and (4), we may set the dimensionless numbers to convenient values (e.g. by slowing down elastic waves and thereby increasing M) and then step the relaxation system forward with moderate time steps as a means of iteratively solving for the steady-state creep velocity. In this case, the elastic wave propagation is less accurate, but serves to allow the velocity to converge to the same steady flow as more conventional methods. Alternatively, if all we are interested in is the elastic behavior, we can use an appropriate time scale for elastic processes and then adjust the Argand number to ‘speed up’ viscous processes to be more comparable to elastic processes.

4. Numerical examples

To make the previous discussion more concrete, we demonstrate that the relaxation system not only converges to the expected Stokes flow solution for an ice shelf, but that without modification, the same model can be used to study elastic effects. In the first two examples, we solve Eqns (29) and (30) subject to a velocity boundary condition at x = 0. We apply a deviatoric stress boundary condition at x = L that corresponds to continuity of normal traction between ice and water such that σ _xx = h ². This boundary condition reflects elastic waves, but other boundary conditions (e.g. absorbing) could also be explored. In the third example, we consider a fully visco-elastic system in two-dimensions with an elaborate rheology designed to mimic failure of ice associated with crevassing and rifting.

Many numerical discretization schemes of varying sophistication have been proposed in the literature as a means of solving the relaxation system (LeVeque, 2002). For illustrative purposes, we use the Lax–Wendroff approach (Appendix B). This fully explicit approach is stable for time steps less than the CFL criterion without the need to ‘upwind’ (LeVeque, 2002). A consequence of this discretization is additional (and unphysical) numerical dissipation: the numerical scheme provides its own pseudo friction to the problem. For most applications, this artificial dissipation spuriously removes energy and is undesirable. However, if all we are interested in is steady-state velocity profiles, the additional dissipation may actually be helpful. For the third example, we also implemented advection of tracers for ice thickness and plastic strain. To accommodate sharp gradients in these quantities, we implemented the superbee flux limiter (LeVeque, 1992). Flux limiters avoid spurious oscillations in hyperbolic problems with ‘shocks’ or discontinuities making them ideal tools to aid in resolving sharp ice-shelf features, like rifts and calving fronts.

The first two examples are illustrated using the same prototypical ice-shelf configuration that corresponds to a steady-state ice shelf with no accumulation/ablation, inflow velocity U _in of 1000 m⋅a⁻¹, inflow ice thickness H ₀ of 1400 m and ice shelf length 80 km. Under these simplified conditions, the steady-state ice thickness h _s and velocity u _s can be found analytically (e.g. van der Veen, 1999, p. 162–163). Material parameters used in all examples are shown in Table 1 and we used grid size Δx = 500 m in all examples.

Table 1. List of the parameter values used for the numerical experiments

4.1. Example 1: convergence of relaxation system to viscous steady-state solution

The purpose of this example is to demonstrate that the relaxation system converges to the appropriate viscous Stokes flow limit. A secondary goal is to examine how many time steps/iterations are required to do so. As we are interested in steady-state viscous velocity profiles, we set (Fr,  Ar) = 1 and apply a perturbation of the form

(45)$$\delta u = A\exp\left[-{\left({x-L\over \Delta x}\right)^2} \right],\; $$

to the steady-state velocity field. In this experiment, we could view model time steps as analogous to the iterations of more conventional iterative schemes, albeit potentially less efficient as we model physical transient behaviors.

Figures (1) and (2) show the ratio of the velocity perturbations to the steady-state velocities (left panels) and the ratio of the deviatoric stress perturbations to the steady-state deviatoric stress (right panels) at three instances of time for two choices of Mach number; M=1 in Figure (1) and M=0.1 in Figure (2). These two cases illustrate the decay of the perturbation to steady state for two different relaxation times. The units of time in both simulations are normalized to the wave velocity so that t = 1 is the time it takes a wave to propagate across the entire domain (at t = 0.5, the wave has propagated halfway). Distance has also been non-dimensionalized to the ice-shelf length. Thus, the CFL criterion to explicitly propagate a wave across the entire domain requires N time steps, where N is the number of grid points.

Figure 1. Decay of a perturbation in the velocity and stress for a freely floating ice shelf for Mach number M = 1. Left panels show the percent difference of the velocity at three times from the steady-state profile. Right panels show the percent difference of the stress at three times from the steady-state profile. Time is non-dimensionalized such that t = 1 corresponds to the time it takes for the pulse to propagate once across the ice-shelf length.

Figure 2. Decay of a perturbation in the velocity and stress for a freely floating ice shelf for Mach number M = 0.1. Left panels show the percent difference of the velocity at three times from the steady-state profile. Right panels show the percent difference of the stress at three times from the steady-state profile. Time is non-dimensionalized such that t = 1 corresponds to the time it takes for the pulse to propagate once across the ice-shelf length.

Figure (1) (with M = 1 and De ~ 1), the velocity and stress perturbations decrease in amplitude as they propagate across the domain. Once the perturbation reaches the leftmost boundary at t = 1, it reflects and begins propagating rightward. After propagating across the domain once, there is still significant energy in the pulse. By contrast, Figure (2) shows the same perturbation for M = 0.1 (De ~ 0.1). With the faster relaxation time, both the velocity and stress perturbations decay faster so that after the pulse has propagated across the domain once, the amplitude of the perturbation has decreased to a relative amplitude of less than 1%.

The above example provides some crude verification of the numerical method and shows that with a judiciously chosen relaxation time, the relaxation method does converge to the appropriate velocity and propagates waves at the expected velocities. In the next example, we illustrate the versatility of the relaxation approach with an example that exploits the ease with which we can also treat elastic problems.

4.2. Example 2: effect of a pulsating ice stream

In this example, we illustrate the versatility of our method using the same discretization as the previous example, but we consider the elastic limit. We impose a smooth acceleration in the velocity into the ice shelf at the grounding line as might be realized, for instance, in a stick-slip event, of the form:

(46)$$u_{{\rm in}} = U_{{\rm in}} + C \exp\left[-\left({t-t_0\over \Delta t}\right)^2\right],\; $$

where C is the amplitude of the perturbation and Δt is the duration. This corresponds to a smooth acceleration in flux into the ice shelf, starting at zero and increasing to a maximum of C before smoothly and symmetrically decelerating to the initial base velocity. For all simulations, we used a perturbation velocity C = 25  km⋅ a⁻¹ and Δt = 90 s. We then choose a velocity scale U ₀ = 30 m s⁻¹ and compute non-dimensional parameters using the values in Table 1, which results in the choice Fr = 1.1 and M = 0.014 and Ar = 3.7 × 10⁻⁶, corresponding to an elastic ice shelf with a small amount of viscous damping (De ~ 25). Choosing different velocity scales affects the numerical values of non-dimensional parameters, but so long as the ratio Fr/M and product Ar ⋅ M remain fixed the results are independent of the velocity scale.

Figures 3(a,b) show the displacement and velocity perturbation as a function of time at the grounding line whereas Figures 3(c,d) show displacement and velocity perturbation as a function of time at the calving front. The displacement and velocity at the calving front smoothly increases. This is followed by a small amplitude ‘ringing’ motion caused by reflection of the pulse off of the calving front. This shows that the same model used to predict steady-state viscous velocities can also be used to simulate ‘seismograms’ for short-time motion. This provides the potential to examine the response of ice sheets to a suite of forcings with overlapping time scales. This should be contrasted with the current approach that uses different models with different assumptions to explore the transient response of ice sheets over different time scales

Figure 3. Evolution of a perturbation to the displacement and velocity at the grounding line (left panels) and margin/calving front (right panels) as a function of time. Panels (a) and (c) show the displacement and velocity perturbations at the grounding line as a function of time and represent the input signature of the perturbation. Panels (b) and (d), by contrast, show the displacement and velocity perturbation as a function of time at the margin/calving front after the signal has propagated across the ice shelf.

4.3. Example 3: ice-shelf failure, rifting and iceberg detachment

For our final example, we return to the challenges of simulating failure and fracture of ice over glaciologically relevant time scales. This time we consider a two-dimensional version of the previous experiments with W = 40 km. Let (x,  y) denote the longitudinal and transverse coordinates and (u,  v) denote the longitudinal and transverse components of the velocity. We assume no inflow along the lateral margins and set the shear stress at the margins assuming a plastic sliding law, whence |σ _xy| = τ _m and v = 0 at y = ±W/2. This is accommodated by applying a Dirichlet boundary condition directly to τ _xy. We further assume an inflow velocity at x = 0 of the form:

(47)$$u( x = 0,\; y) = 4 U_{{\rm in}} ( 1-f) \left({y\over W}\right)^2 \, + \, fU_{{\rm in}}.$$

Here U _in = 1000 m a⁻¹ represents the peak inflow velocity at the grounding line and f = 0.7 parameterizes the magnitude of velocity decrease along the margins.

Following Bassis and others (2022), we adopt a visco-plastic rheology with a failure strength of the form η = min (η _p,  η _v) where η _v denotes the two-dimensional generalization of the effective viscous rheology defined by Eqn 6:

(48)$$\eta_{\rm v} = {B^n\over 2}\tau_{\rm e}^{1-n},\; $$

with $\tau _{\rm e} = h^{-1}\sqrt {\sigma _{xx}^2 + \sigma _{yy}^2 + \sigma _{xy}^2 + \sigma _{xx}\sigma _{yy}}$. Writing the effective viscous strain rate $\dot \epsilon _{\rm e} = ( \tau _{\rm e}/B) ^{-n}$, we define a failure strength τ _y and write the plastic viscosity $\eta _{\rm p} = {\tau _y\over \dot \epsilon _{\rm e}}$. Again, following Bassis and others (2022), we promote failure localization by tracking plastic strain $\epsilon _{\rm p}$ accumulated in faults and fractures and assuming the failure strength decreases with plastic strain:

(49)$$\tau_y = \max\left\{\tau_{{\rm c}}-\left(\tau_{{\rm c}} -\tau_{{\rm min}}\right){\epsilon_{\rm p}\over \epsilon_{{\rm crit}}},\; \tau_{{\rm min}}\right\}.$$

Here τ _c = 261 kPa (0.75 in dimensionless units) is the intact strength of ice, $\epsilon _{{\rm crit}} = 0.0125$ denotes the limit where ice is completely broken and τ _min = 17.4 kPa (0.05 in dimensionless units) represents any residual strength of failed ice and accounts for mélange filling rifts. We evolve $\epsilon _{\rm p}$ assuming $\dot \epsilon _{\rm p} = \dot \epsilon _{\rm e}$ where η _p < η _v. We further set $\epsilon _{\rm p} = 0$ within 10 km of the grounding line, assuming intact ice is flowing into the ice shelf and to avoid edge effects and artifacts associated specifying the velocity boundary condition.

We evolve the ice thickness, velocity and a scalar variable $\epsilon _{\rm p}$ that tracks accumulated strain in faults and fractures. Simulations are conducted with dimensionless parameters M = 0.1, Ar = 1 and Fr = 1.

Figure 4(a) shows the initial steady-state ice thickness and velocity obtained by time stepping the system of equations to steady state without failure. The subsequent panels show snapshots of the ice thickness and plastic strain $\epsilon _{\rm p}$ after failure is permitted, illustrating the initiation and propagation of rifts. In this example, rifts episodically initiate from the margins and propagate across the domain, eventually isolating an iceberg where a low viscosity and ice thickness effectively uncouple the iceberg from the shelf (see Supplementary Animation M1). However, the calving front remains quasi-stable throughout these simulations. Numerically, however, what is intriguing about these simulations is that we have used a purely explicit method that completely avoids solving systems of non-linear equations. Implementing an analogous elliptical solver for the velocity converged slowly, typically taking dozens of iterations to converge. Moreover, unlike the simulations in Bassis and others (2022), we did not need to introduce an artificial ‘water drag’ to avoid ill-conditioning when icebergs detach in the aftermath of rift propagation.

Figure 4. Evolution of rifts in an idealized ice shelf. Panel a shows a snapshot of the steady state ice thickness and speed at the beginning of the simulation. Panels b–d show snapshots at different points in time. Contours show fully failed regions. Supplementary animation M1 shows an animation over several centuries illustrating the sporadic detachment of bergs.

5. Discussion

One of the long-standing tenets of ice dynamics is that ice sheets deform slowly, like highly viscous fluids. The assumption of slow and steady flow results in the hierarchy of approximations to Stokes flow currently implemented in numerical ice-sheet models (e.g. Greve and Blatter, 2009). Here, we have broken with both the assumption that ice behaves purely like a viscous fluid and the assumption that flow is slow and steady. Adding visco-elasticity and re-introducing the acceleration terms into the force balance transform the equations of motion into a hyperbolic system where velocity is locally determined and information propagates at the finite elastic wave speed. The quasi-static viscous limit treated by traditional ice-sheet models corresponds to the flow in the visco-elastic system once the transient behaviors have attenuated. Because of this connection, if we are not interested in the dynamic elastic behavior at short time scales, we can treat the elastic wave speed as a numerical parameter with the time steps analogous to iterations in the non-linear solvers typically used in more traditional methods to find the velocity field.

Although including acceleration into the force balance may seem radical and, perhaps, ill-advised, recent studies have similarly re-introduced the acceleration term back into Stokes models as a numerical regularization needed to avoid numerical issues associated with ill-posedness from applying the flotation condition to 3D models of ice shelves (Berg and Bassis, 2022) and as a numerical trick to efficiently solve the non-linear Stokes equations that can be efficiently implemented on graphical processing units (GPUs). For example, Räss and others (2020) introduced the pseudo-inertial terms as a means of creating a matrix-free full Stokes solver that can be implemented on GPUs. The introduction of pseudo-time creates a purely viscous parabolic system that can be marched to steady state as a Stokes solver. This should be contrasted with our purely hyperbolic visco-elastic system where the time steps can correspond to transient elastic and visco-elastic behavior. The steady-state behavior, however, is identical in both methods.

The number of operations needed for a single time step in our explicit visco-elastic system scales like the number of degrees of freedom. This scaling is comparable to state-of-the-art algebraic multigrid methods. However, multiple time steps – iterations – are typically needed to converge. Detailed comparisons are often highly problem – and implementation – dependent. Given the ability to craft bespoke elliptic solvers, it seems unlikely that solving an explicit hyperbolic system will be more efficient than specialized non-linear Stokes solvers that only seek the long-time viscous flow of the ice. However, when there are fast time scales in the system or discontinuities, like when fractures are present, the explicit hyperbolic system may be advantageous and able to represent a wider range of physical processes. To be clear, although most ice-sheet modelers eventually confront situations where the elliptical velocity solvers maddeningly fail to converge, hyperbolic systems have their own intricacies, including ringing and spurious oscillations that also need to be managed.

In our treatment, we assumed that the ice thickness was constant over the time scales of wave propagation. This effectively filtered out gravity waves. We could build in more sophisticated treatments of elasticity in conjunction with more sophisticated force balances to simulate the full bestiary of waves possible in the system. Some of these waves, such as flexural gravity waves, have been implicated in ice-shelf demise (Banwell and MacAyeal, 2015; Lipovsky, 2018; Massom and others, 2018; Aster and others, 2021) and could be incorporated by relaxing the assumption that the ice is locally hydrostatically compensated. Modeling these effects may provide an avenue to simulate the effect of short-term environmental drivers on rifting and calving and explosive ice-shelf collapse (Bassis and others, 2008; Walker and others, 2013, 2015; Banwell and MacAyeal, 2015; Wille and others, 2022).

Both the promise and obstacles in simulating ice dynamics using hyperbolic visco-elastic systems are illustrated in our simulations of ice failure and rifting. Numerical accuracy in these simulations demanded time steps on the order of hours to days. The time step size of these simulations already pushes into the elastic regime. By contrast, in more traditional (elliptical) approaches (e.g. Bassis and others, 2022), each time step requires solving large systems of non-linear equations that converge slowly (if at all). The approach we used here sidestepped the convergence issue by treating the (highly non-linear) rheology explicitly, with the added benefit that the elastic component of the rheology serves as an added regularization (Hunke and Dukowicz, 1997).

Using elasticity to regularize the system and avoid numerical headaches has analogies in many fields. In fact, Cattaneo (1958) recognized that information travels at infinite speed in the heat equation and applied an analogous transformation to Newton's law of cooling. But the analogy that is closest in spirit to that of the one we propose is that of the elastic-visco-plastic sea-ice model. Here, Hunke and Dukowicz (1997) recognized that introducing an elastic element to the visco-plastic rheology used for sea-ice rheology was a numerically efficient trick that avoided numerical headaches associated with nearly rigid motion when the effective viscosity becomes large. However, in sea ice, inertial terms associated with sea-ice drift are part of the dominant balance, whereas inertial terms play a negligible role in the Stokes flow that traditional ice-sheet models invoke.

Although speculative, it might be possible to build on the connection between sea-ice and ocean models and simply set the elastic wave speed within the ice to be comparable to that of ocean or sea-ice dynamics and consider ice sheets and ice shelves as simply another ‘layer’ or component within an ocean or sea-ice model and co-evolve the ice and ocean together. Co-evolving the ice and ocean could potentially allow grounding lines (and the ocean cavity) to dynamically evolve with both tidal forcing and basal melting. In this type of framework, we could naturally simulate the interaction between the ice and shorter-term environmental drivers and even the evolution of icebergs as they detach from the ice shelf, drift and melt, one of the largest current challenges in coupling ice sheets into climate system models.

6. Conclusions

The challenge that ice-sheet modelers face is that ice-sheet evolution is controlled by internal and external processes with time scales that range from seconds to millennia. For example, short-term atmospheric and oceanic drivers, like atmospheric rivers and ocean swell, have been hypothesized to play a role driving ice-shelf collapse (e.g. Massom and others, 2018; Francis and others, 2021; Wille and others, 2022). As computational resources continue to expand, numerical ice-sheet models have progressed toward higher spatial and temporal resolution. As the spatial and temporal resolution of ice-sheet models continues to improve, the disparity between the fundamental time scales in an ice-sheet model and environmental forcing will continue to shrink. However, even in the most optimistic scenario, the disparity in time scales associated with brittle, elastic failure and longer-term ice-sheet evolution is likely to remain stark. Here, we argued that we can bridge this divide by ‘slowing down’ the fast time scale processes, associated with elastic waves so that these processes are more comparable to the slower viscous processes that ice-sheet models traditionally represent. This change allowed us to simulate short-term elastic and longer-term viscous effects within a single modeling framework.

Our approach could also be valuable in more regional, process-based models. For example, discrete element models resolve the fast elastic brittle failure (Bassis and Jacobs, 2013; Astrom and others, 2014). Time step restrictions, however, limit these simulations to durations of minutes to hours. So long as the longer term dynamics is not controlled by elastic wave propagation, by slowing down the fastest propagating waves it may be possible to extend simulations to weeks or even years. Similarly, given the vast disparity in time scales between ice-sheet evolution and ocean forcing, leaning into the faster time scales of ice-sheet evolution may provide a bridge where ice-sheet models are co-evolved with ocean models to allow finer scale coupling of dynamic ocean cavities and grounding line migration. This latter option, although speculative, could provide a more direct opportunity to directly simulate the effect of both short-term and longer-term atmospheric and oceanic forcing on ice-sheet evolution. Although we have focused on ice shelves and ice sheets, the system we propose is general and could also be applied to other systems. For example, we could slow down acoustic waves in the ocean as proposed by Salmon (2009) for studies of ocean circulation or use a similar formulation for geodynamics simulations that seek to simulate both brittle failure and longer term deformation.