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.


Introduction
Stokes simulations are used in glaciology as a tool to determine the time evolution of glaciers (e.g., Gagliardini et al., 2013). 10 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 15 solutions that are 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 (illustrated in Fig.   20 1), where the upper portion has calved away, emulating the "footloose" mechanism proposed by Wagner et al. (2014) where a aerial portion of the calving front first detaches.

Glacier Stress Balance
Denoting the velocity field by u(x, z, t) = (u x (x, z, t), u z (x, z, t)) and pressure by P (x, z, t), conservation of linear-momentum 25 can be written in the form: The Cauchy stress is defined in terms of strain rate, effective viscosity, pressure, and the identity matrix I: with strain rate tensor ε given by: Here ρ i is the density of ice, g is the acceleration due to gravity, and η is the effective viscosity of ice: The effective viscosity is a function of the effective strain rate ε e = √ ε ij ε ij /2, a temperature dependent constant B, and the flow-law exponent n = 3; the acceleration term on the right hand side of Eq. (1) denotes the material derivative.

35
In the Stokes limit, we drop the acceleration term from Eq. (1), an approximation which is justified for most glaciological applications (Greve and Blatter, 2009). However, as we shall show, this assumption is problematic for applications where the ice departs from hydrostatic equilibrium.

Boundary Conditions
To illustrate an example where the Stokes flow problem becomes ill-posed, we consider a two-dimensional floating ice shelf 40 ( Fig. 1). We specify the normal component of the velocity u ·n = u 1 where u 1 is a constant along the inflow portion of the domain (Γ 1 ). At the ice-atmosphere boundary (Γ 2 ) we assume the surface is traction free. At the boundary between ice and ocean (Γ 3 ) the shear traction along the ice-interface vanishes and continuity of normal traction along the ice-ocean interface can be written as σ nn (x, t) = −ρ 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 45 ice-ocean interface cannot be balanced by internal stresses. In this case there is no solution and vertical velocities are singular.
In reality, of course, the ice will quickly re-adjust to hydrostatic equilibrium through rapid buoyant uplift through the (nearly) inviscid ocean.
We can more accurately specify the boundary condition for Stokes flow at the ice-ocean interface by writing it in the form: where ∆z(x, t) is an a priori unknown uplift that must be determined as part of the solution to enable the full force balance to close.
The additional uplift term ∆z has a simple physical explanation: if normal stress was exactly hydrostatic, σ nn = ρ i gH 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. The Stokes limit is more complex as internal stresses also contribute to the normal stress at 55 the ice-ocean interface, but the location of the ice-ocean interface needs to be solved for as part of solution to the problem, which we examine next.

Numerical Stabilization of Buoyant Uplift
Different numerical methods apply different techniques to solve for ∆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 60 Taylor series of vertical position of the ice-ocean interface: This Taylor series transforms the uplift into a time step dependent Newtonian velocity damping: However, the coefficient of the damping force is proportional to ∆t and vanishes in the limit of small ∆t. In this limit, vertical 65 velocities are singular. We refer to the damping method given in Eq. (7) as the "sea-spring" method based on the nomenclature used in Elmer/Ice documentation. With this method, we can decompose the velocity u into a "viscous" component u visc and a hydrostatic uplift component ∆z uplif t , which we write in the form:

70
where ∆z uplif t is the vector form of the same displacement written in Eq. (5). Due to the dependence of the ice-ocean boundary stress on time step size ∆t , the total velocity becomes time step dependent.
Inspecting Eq. (8) shows that in the limit of large ∆t, the uplift term becomes small compared to the viscous velocity. Thus, the sea-spring damping method can provide a good approximation for the viscous velocity as long as (∆z uplif t )/(u visc ∆t) << 1, which is true if the glacier is close to hydrostatic equilibrium (∆z uplif t small) or a sufficiently long time step is used.

Test Design
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 with inflow velocity for the upstream boundary condition set to 4 km a −1 . These thickness and velocity parameters are broadly consistent with observations for the last 10 km of Pine Island ice shelf (Rignot et al., 2017, 80 2011; Mouginot et al., 2012;Paden et al., 2010Paden et al., , updated 2018. The temperature dependent constant in Glen's flow law is chosen to be 1.4 × 10 8 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 85 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 (Alnaes et al., 2015), an open source finite element solver with a Python interface that has been previously used for Stokes glacier modeling ( Figure 2 shows the sensitivity of the velocity field, effective strain rate, effective shear stress, and greatest principle stress to time step size when using the sea-spring boundary condition from Eq. (7). Furthermore, because of the coupling between 95 effective strain rate and viscosity, the viscosity for the majority of the domain becomes unphysically small as the time step decreases. This positive feedback between effective strain rate and viscosity is especially problematic as higher strain rate causes lower viscosity, and vice versa, leading to unphysical results. This problem can be alleviated by using a viscoelastic rheology when examination of short time scale behavior is desired. However, even for a purely viscous model, short time steps may be necessary to satisfy numerical stability criteria during hydrostatic adjustment that momentarily forces the model 100 outside of the Stokes range.

Divergent Behavior
In addition to the divergence of the L 2 norm, we also examine the maximum effective shear stress and greatest principle stress (Fig. 3). Maximum values may be a better predictor of the effect of time step dependence on the output 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.

4 Proposed Solution -Reintroduce Acceleration Term into Stress Balance
The velocity solution is ill-posed because the neglected acceleration term is not actually small relative to the other terms in Eq. (1): large velocities associated with hydrostatic adjustment rapidly change on time scales short compared to the internal deformation of the ice. We therefore reintroduce the acceleration term and show that this regularizes the solution for small time 110 steps. We use a simple first order backwards differentiation scheme in a Lagrangian reference frame where This effectively introduces a Newtonian damping term on the entire body of the glacier where the damping coefficient C = 1/∆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 115 unstable) for long time steps. Therefore, we propose to use both damping terms so that the system of equations is numerically accurate for small time step sizes and the velocity converges to the viscous limit for large time step sizes. Although this method rectifies the numerical inaccuracies present at short time steps with sea-spring damping, it does not address physical inaccuracies from using a rheology not suited for elastic effects. However, the numerical divergence exists independent of rheology and would need to be addressed even for a viscoelastic model.

120
When we include both damping terms, effective strain rate, effective principle stress, and greatest principle stress are consistent for both small and large time steps (Fig. 2). At small time steps the acceleration term dominates, so the sea-spring with physical evolution of the system: the solution resolves the acceleration and deceleration of the glacier as it evolves towards a steady-state as opposed to being a consequence of an ill-posed problem.
Notably, the sea-spring solution shows larger maximum stresses than the sea-spring with Navier Stokes solution (Fig. 3).
This is particularly evident for the effective shear stress, which is overestimated by an order of magnitude at the shortest time 130 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 over prediction of stresses could cause a calving model to predict unphysical calving events due to numerical inaccuracies.

Conclusions
Our study shows that using a common numerical stabilization method of the ice-ocean boundary in Stokes glacier modeling Competing interests. The authors declare that they have no conflicts of interest.