Brief communication: Time step dependence (and fixes) in Stokes simulations of calving ice shelves

The buoyancy boundary condition applied to floating portions of ice sheets and glaciers in Stokes models requires special consideration when the glacier rapidly departs from hydrostatic equilibrium. This boundary condition can manifest in velocity fields that are unphysically (and strongly) dependent on time step size, thereby contaminating diagnostic stress fields. This can be especially problematic for models of calving glaciers, where rapid changes in geometry cause configurations that suddenly depart from hydrostatic equilibrium and lead to inaccurate estimates of the stress field. Here we show that the unphysical behavior can be cured with minimal computational cost by reintroducing a regularization that corresponds to the acceleration term in the stress balance. This regularization provides consistent velocity solutions for all time step sizes.

1 Introduction

Stokes simulations are used in glaciology as a tool to determine the time evolution of glaciers (e.g., Gagliardini et al., 2013). Increasingly, these models are also used to examine the stress field within glaciers to better understand factors that control crevasse formation and the onset of calving events (Ma et al., 2017; Benn et al., 2017; Nick et al., 2010; Todd and Christoffersen, 2014; Ma and Bassis, 2019). This type of model can provide insight into the relationship between calving, climate forcing, and boundary conditions (e.g., Todd and Christoffersen, 2014; Ma et al., 2017; Ma and Bassis, 2019).

Here we show that a common method used to implement the ice–ocean boundary condition in Stokes models can result in solutions that are unphysically sensitive to the choice of simulation time step size. This behavior manifests in applications that allow for rapid changes in the model domain – a type of change associated with models that allow for instantaneous calving events or crevasses (Todd and Christoffersen, 2014; Todd et al., 2018; Yu et al., 2017). The time step dependence arises because for glaciers outside of hydrostatic equilibrium, the acceleration is not small as assumed in Stokes flow. We illustrate both the issue and the solution using an idealized ice shelf geometry (Fig. 1) where the upper portion has calved away, emulating the “footloose” mechanism proposed by Wagner et al. (2014) where an aerial portion of the calving front first detaches.

2 Problem description

2.1 Glacier stress balance

Denoting the velocity field in two dimensions by $\mathit{u}(x,z,t)=\left(u_x\right(x,z,t),u_z(x,z,t\left)\right)$ and pressure by $P(x,z,t)$, conservation of linear momentum can be written in the form

$$\begin{array}{}\text{(1)}& \mathbf{\nabla }\cdot \mathit{\sigma }+{\mathit{\rho }}_i\mathit{g}={\mathit{\rho }}_i{\displaystyle \frac{D\mathit{u}}{Dt}},\end{array}$$

where D∕Dt denotes the material derivative. The Cauchy stress is defined in terms of strain rate, effective viscosity, pressure, and the identity matrix I:

$$\begin{array}{}\text{(2)}& \mathit{\sigma }=\mathrm{2}\mathit{\eta }\mathit{ϵ}-P\mathbf{I},\end{array}$$

with strain rate tensor ϵ given by

$$\begin{array}{}\text{(3)}& {\mathit{ϵ}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial r_j}}+{\displaystyle \frac{\partial u_j}{\partial r_i}}\right).\end{array}$$

Here ρ_i is the density of ice, g is the acceleration due to gravity, and η is the effective viscosity of ice:

$$\begin{array}{}\text{(4)}& \mathit{\eta }={\displaystyle \frac{B}{\mathrm{2}}}{\mathit{ϵ}}_e^{\frac{\mathrm{1}}{n}-\mathrm{1}}.\end{array}$$

The effective viscosity is a function of the effective strain rate ${\mathit{ϵ}}_e=\sqrt{{\mathit{ϵ}}_{ij}{\mathit{ϵ}}_{ij}/\mathrm{2}}$, a temperature-dependent constant B, and the flow-law exponent n=3.

In the Stokes approximation, we drop the acceleration term on the right-hand side of Eq. (1), which is justified for most glaciological applications (e.g., Greve and Blatter, 2009). However, as we shall show, this assumption is problematic for applications where the ice rapidly departs from hydrostatic equilibrium.

Figure 1 A diagram showing the boundary conditions of an idealized floating ice shelf. The ice–ocean interface is subject to two normal stresses – the depth-dependent water pressure and the numerical damping force for stabilization to hydrostatic equilibrium. The dashed red line illustrates an iceberg that breaks off from the top of the calving front (exaggerated), reducing the freeboard and instantaneously perturbing the ice shelf from hydrostatic equilibrium.

2.2 Boundary conditions

To illustrate an example where the Stokes flow problem rapidly departs from hydrostatic equilibrium, we consider a two-dimensional floating ice shelf (Fig. 1). Along the inflow portion of the domain (Γ₁), shear traction vanishes and we specify the normal component of the velocity $\mathit{u}\cdot \widehat{\mathit{n}}=u_{\mathrm{1}}$, where u₁ is a constant and $\widehat{\mathit{n}}$ is the normal vector along Γ₁. At the ice–atmosphere boundary (Γ₂) we assume the surface is traction free. At the boundary between ice and ocean (Γ₃) the shear traction along the ice interface vanishes, and continuity of normal traction along the ice–ocean interface can be written as ${\mathit{\sigma }}_{nn}(x,t)=-{\mathit{\rho }}_{w}gb(x,t)$, where b(x,t) is the position of the ice–ocean interface.

Problems arise with this form if the glacier is not exactly in hydrostatic equilibrium because buoyancy forces along the ice–ocean interface cannot be balanced by internal stresses. In this case, there is no unique solution. In reality, of course, inertial effects are not negligible and the ice will quickly readjust to hydrostatic equilibrium as a consequence of buoyant uplift through the (nearly) inviscid ocean.

We can more appropriately specify the boundary condition for Stokes flow at the ice–ocean interface by writing it in the form

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}_{nn}(x,t)=-{\mathit{\rho }}_{w}g\left[b(x,t)+\mathrm{\Delta }z(x,t)\right]\text{on}{\mathrm{\Gamma }}_{\mathrm{3}},\end{array}$$

where Δz(x,t) is an a priori unknown isostatic adjustment that could potentially include a rigid body translation and, for a freely floating iceberg, rigid body rotation. Crucially, the rigid body motion must be determined as part of the solution to enable the local and global force balance to close.

The additional isostatic adjustment term Δz has a simple physical explanation: if normal stress was exactly hydrostatic, σ_nn=ρ_igH, where H is the ice shelf thickness. Equation (5) can then be solved for Δz to determine the position of the bottom interface needed for the forces to balance, which is exactly what is done in the shallow shelf approximation. The full Stokes approximation is more complex as internal stresses also contribute to the normal stress at the ice–ocean interface, but the location of the ice–ocean interface needs to be solved for as part of the solution to the problem, which we examine next.

2.3 Numerical stabilization of buoyant uplift

Different numerical methods use different techniques to estimate Δz in Eq. (5). For example, in Elmer/Ice, a popular package for modeling Stokes glacier flow, Durand et al. (2009) proposed an ingenious solution in which Δz is estimated based on a Taylor series of vertical position of the ice–ocean interface:

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }z=u_z(x,t)\mathrm{\Delta }t+O\left(\mathrm{\Delta }{t}^{\mathrm{2}}\right).\end{array}$$

This Taylor series transforms the isostatic adjustment into a time-step-dependent Newtonian velocity damping:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{nn}(x,t)=-{\mathit{\rho }}_{w}g\left[b(x,t)+u_z(x,t)\mathrm{\Delta }t\right].\end{array}$$

Here, the velocity u_z can include rigid body translation. The coefficient of the damping force in this approximation is proportional to the time step size Δt and vanishes in the limit of small Δt. We refer to the damping method given in Eq. (7) as the “sea-spring” method based on the nomenclature used in Elmer/Ice documentation.

We illustrate the time step dependence using a floating ice shelf as an example. In this case, global force balance is not guaranteed, leading to a formally ill-posed problem. However, the lack of a solution (which results in an exceptionally large velocity using the sea-spring stabilization) persists even when grounded ice is included in the domain because the velocity (and strain rate) solution can become unphysically sensitive to the position of the ice–ocean interface.

3 Calving-based convergence test

Figure 2 L₂ norm of (a) vertical velocity, (b) effective strain rate, (c) maximum shear stress, and (d) greatest principal stress immediately after the emulated calving event. The norm is calculated based on cell-averaged values. Solutions after a single time step are shown for sea-spring damping and sea-spring with Navier–Stokes (NS) for time step size ranging from 1 s to 30 years. The sea-spring with Navier–Stokes term gives physical solutions for all time steps but still shows variability with a time step consistent with the evolution of the system.

For our test, we implement an idealized rectangular ice shelf of thickness 400 m and length 10 km. This ice shelf is initialized to be in exact hydrostatic equilibrium. We set the inflow velocity for the upstream boundary condition to 4 km a⁻¹. These thickness and velocity parameters are broadly consistent with observations for the last 10 km of Pine Island ice shelf (Rignot et al., 2017, 2011; Mouginot et al., 2012; Paden et al., 2010). The temperature-dependent constant in Glen's flow law is chosen to be 1.4×10⁸ Pa s^1∕3, the value given by Cuffey and Paterson (2010) for −10 ^∘C.

To emulate the occurrence of a calving event that would perturb the ice shelf from hydrostatic equilibrium, a rectangular section of length 50 m and thickness 20 m is removed from the upper calving front of the glacier (Fig. 1). This type of calving behavior has been proposed as the trigger of a larger calving mechanism related to buoyant stresses on the ice shelf (Wagner et al., 2014). The numerical effects we document are not unique to this style of calving, and this mechanism is only meant to illustrate the numerical issues.

The problem is implemented in FEniCS (Alnæs et al., 2015), an open-source finite element solver with a Python interface that has been previously used for Stokes glacier modeling (Ma et al., 2017; Ma and Bassis, 2019). The problem is solved using an arbitrary Lagrangian–Eulerian formulation using mixed Taylor–Hood elements with quadratic elements for velocity and linear elements for pressure. The open-source finite element mesh generator Gmsh is used to generate an unstructured mesh with uniform grid spacing of 10 m near the calved portion of the domain and grid spacing of 40 m elsewhere.

In Fig. 2, we show the L₂ norm of vertical velocity, effective strain rate, maximum shear stress, and greatest principal stress immediately after the calving event, but computed using different values for Δt in Eq. (7). This shows the sensitivity of the vertical velocity, effective strain rate, maximum shear stress, and greatest principal stress to time step size when using the sea-spring boundary condition. In addition to the behavior of the L₂ norm, we also examine the maximum value of the maximum shear stress and greatest principal stress (Fig. 3). Maximum values may be a better predictor of the effect of time step dependence on the results of Stokes calving models. Because calving models often assume that calving is likely if a stress threshold is exceeded (Ma et al., 2017), outliers in stress are more important than a stress averaged over the entire domain.

Figure 3 Maximum (a) maximum shear stress and (b) greatest principal stress immediately after the emulated calving event. The maximum is calculated based on cell-averaged values. Solutions after a single time step are shown for sea-spring damping and sea-spring with Navier–Stokes (NS) for time step size ranging from 1 s to 30 years. Without corrections, maximum stress is overestimated in the Stokes model, which could lead to overestimation of glacier retreat due to calving.

4 Proposed solution – reintroduce acceleration term into stress balance

The velocity solution is unphysical because the neglected acceleration term is not actually small relative to the other terms in Eq. (1). This is because large velocities associated with hydrostatic adjustment rapidly change on timescales that are short compared to the internal deformation of the ice. It may be possible to separate a rigid body translation and rotation that satisfies global force and torque balance from the internal deformation, but this quickly becomes cumbersome and impractical when we include the potential for ice to break: global force and torque balance would have to be maintained on each intact portion of ice. Instead, we reintroduce the acceleration term directly to the Stokes equation and show that this regularizes the solution for small time steps. We use a simple first-order backwards differentiation scheme in a Lagrangian reference frame:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{D{\mathit{u}}_i}{Dt}}={\displaystyle \frac{{\mathit{u}}_i-{\mathit{u}}_{i-\mathrm{1}}}{\mathrm{\Delta }t}},\end{array}$$

where u_i and u_i−1 denote the velocity of Lagrangian fluid parcels at the current time step and the velocity at the previous time step, respectively. In the example considered here, the fixed horizontal velocity at the upstream boundary does not permit rigid body rotation and Eq. (8) remains a valid approximation even in an Eulerian coordinate system. As an initial condition on velocity, we assume a uniform velocity field of 4 km a⁻¹ (equal to the inflow velocity) in the horizontal direction and zero in the vertical direction. The choice of zero initial velocity in the vertical direction is consistent with the idea that the ice shelf has instantaneously been perturbed from hydrostatic equilibrium.

Restoring the inertial term effectively introduces a Newtonian damping term on the entire body of the glacier where the damping coefficient is $C=\mathrm{1}/\mathrm{\Delta }t$. Computational difficulty is not impacted by reintroducing the acceleration term in this way because the term is linear with respect to velocity. However, unless a fully implicit scheme was implemented, the solution becomes inaccurate (and unstable) for long time steps. Therefore, we propose to use both damping terms so that the system of equations is numerically accurate for all time step sizes.

When we include both damping terms, vertical velocity, effective strain rate, maximum shear stress, and greatest principal stress maintain physical values for both small and large time steps (Fig. 2). At small time steps the acceleration term dominates and the sea-spring with Navier–Stokes solution departs from the sea-spring solution. In this limit, the (nearly) rigid body isostatic adjustment of the ice shelf dominates the solution. By contrast, for large time steps, the sea-spring damping dominates and the two methods overlap. At intermediate time steps, both damping terms contribute as the solution transitions from the regime dominated by inertial effects to one where inertial effects are small. It is crucial to note that although the sea-spring with Navier–Stokes solution retains time step dependence for small time steps, the time step dependence now results from the physical evolution of the system: the solution resolves the acceleration and deceleration of the glacier as it bobs up and down in the ocean and approaches a steady state.

Although the sea-spring solution shows a smaller L₂ norm of greatest principal stress than the sea-spring with Navier–Stokes solution, this is due to large negative compressive stresses associated with being outside of hydrostatic equilibrium. If we instead examine the maximum of the stress fields, the sea-spring solution shows larger values for both maximum shear stress and greatest principal stress (Fig. 3). This is particularly evident for the maximum shear stress, which is overestimated by an order of magnitude at the shortest time step tested. In the footloose calving mechanism, when a portion of the upper calving front is removed, the front of the ice shelf becomes buoyant and produces increased shear stress upstream on the ice shelf (Wagner et al., 2014). This overprediction of stresses could cause a calving model to predict unphysical calving events due to numerical inaccuracies.

5 Conclusions

Our study shows that using a common numerical stabilization method of the ice–ocean boundary in Stokes glacier modeling, there is an explicit time step dependence of the solution that is unphysical for small time steps when the domain rapidly departs from hydrostatic equilibrium. For model applications where changes in the domain are only due to viscous flow, the time step dependence is not problematic as long as domains are (nearly) in hydrostatic equilibrium at the start of simulation. However, for applications where rapid changes to the model domain occur, such as when calving rules are implemented, sudden departure from hydrostatic equilibrium is not only possible but expected. In these cases, time step dependence of the solution will appear. This can contaminate solutions of the stress after calving, potentially leading to a cascade of calving events and an overestimate of calving flux if numerical artifacts are not addressed. However, the time step dependence can be easily cured with little computational cost by reintroducing the acceleration term to the Stokes flow approximation. The acceleration term regularizes the solution for small time step sizes and provides consistent solutions for all time steps.