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, Reference MacAyeal1989; Weis and others, Reference Weis, Greve and Hutter1999):

(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, Reference MacAyeal1989; Weis and others, Reference Weis, Greve and Hutter1999).

(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, Reference Pedlosky1987).

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, Reference Cuffey and Paterson2010). 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, Reference Zarrinderakht, Schoof and Peirce2022).

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 (Reference Berg and Bassis2022), 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, Reference Cattaneo1958). 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, Reference Logan, Lavier, Choi, Tan and Catania2017; Räss and others, Reference Räss, Licul, Herman, Podladchikov and Suckale2020). 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.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, Reference Malvern1969) 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.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 (Reference Bassis, Berg, Crawford and Benn2022), 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 (Reference Bassis, Berg, Crawford and Benn2022), 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 (Reference Bassis, Berg, Crawford and Benn2022), 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.