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

This paper presents a simple method to overcome time step dependence of the solution arising when solving for an ice-shelf which departs significantly for hydrostatic equilibrium. This could be the case for instantaneous non-vertical icebergs calving or supraglacial lake drainage. This is quite a technical paper but as the problem might be encountered by other groups using different Stokes solvers, this brief communication certainly deserves to be published. The overall writing of the paper is quite good even if I think that there is some room for improvement.

This is a good point and we agree that the original figures were confusing. We have added text to the Figure 2 caption clarifying that the time step variability for the sea-spring with Navier Stokes solution is connected to the time evolution of the system. In the classic Stokes system, the velocity depends solely on the geometry of the ice shelf/glacier and internal properties (e.g., temperature). Hence, the dependence of the velocity field on time step is unphysical. However, when solving the full Navier-Stokes system, the velocity becomes time dependent and, like all numerical ODE integrations, we must take sufficiently small time steps to ensure numerical convergence when integrating with respect to time to find the numerical approximation to the solution. For the sea-spring + NS solution, the velocity tends towards zero as time step size decreases. As a consequence there is no deformation and, for very small times, the ice shelf behaves as a nearly rigid body as it approaches hydrostatic equilibrium. In the sea-spring + NS method, taking small time steps allows us to resolve the quasi-rigid body uplift of the ice shelf as it "bobs" in the water. We have added text in Section 4 of the manuscript clarifying the rigid body behavior at short time steps (lines 132-135). This is illustrated below, where we show a plot of the L2 norm of the greatest principal stress for the sea-spring + NS method. The points represent different time step sizes (as in the manuscript). The dashed line is created by choosing the smallest time step and numerically integrating the system in time. For small times, the solution obtained from taking a single step with different time step size and the solution obtained from numerical integration (with a very small time step size) are similar, but begin to differ as time step size increases because taking large time steps results in less accurate solutions. Thus, the variation from the sea-spring + NS solution at small time steps is consistent with the actual time evolution of the system. This is connected to our discussion in section 4 of the manuscript.
The viscosity has no timestep dependence for the sea-spring solution and it is the sea-spring+NS solution that has no time step dependence for effective strain-rate.
We have chosen to omit plotting the viscosity form the final manuscript, instead replacing it with a plot of vertical velocity because effective strain rate and viscosity display somewhat redundant information and to ease the exposition of this short manuscript. Showing the vertical velocity directly provides a more direct illustration of the problem with the vertical velocity becoming unphysically large as time step size decreases.
However, we include a plot of viscosity below that shows the minimum and maximum value of the viscosity at different time steps for the two methods. Plotting the maximum and minimum values better highlights the time step dependence of the viscosity at small time steps. Note that for small time step sizes (and times) with the sea-spring + NS solution, the motion is quasi-rigid (i.e., deformation is small), strain rates are small and the viscosity increases. The viscosity is ultimately limited by the regularization we use for the rheology at small strain rates. By contrast, for the sea-spring solution, as time step size decreases the increasingly large vertical velocity causes a large strain rate increase. Because the viscosity and strain rate are inversely related, this results in decreasing viscosities. This is even less clear for stress where both solutions are diverging but presents both a timestep dependence. I would expect more comments on the text on this and how from the figure one can decide which is the working solution.
We have clarified our language and removed most uses of the word "divergent" from the manuscript. The vertical velocity, which we now show in Figure 2 panel (a), becomes unphysically large for small time step sizes for the sea-spring method. Effective strain rate becomes approximately zero for small time step sizes for the sea-spring + NS method because the motion at small times is nearly rigid body. Effective shear stress tends to zero for the seaspring + NS method for the same reason as the effective strain rate. We also point out in section 4 that while Figure 3 shows higher maximum stresses for the sea-spring method, the L2 norm shows higher stresses for the sea-spring + NS method because of high negative compressive stresses (lines 136-137).

Smaller Points:
page 2, line 41: ", where u1 is a constant" Comma has been added. It is not clear from the text and the captions if what is plotted on these figures is the solution after the timestep following the calving event.
Thanks for pointing this out. We now state in the figure captions that the plotted solution is immediately after the (emulated) calving event.
What are the differences of setup between Fig. 2 c and d and Fig. 3?
The difference between Fig. 2 c and d and Fig. 3 is that we plot the L2 norm of the solution in Figure 2 and the L1 norm (maximum) of the solution in Figure 3. This is emphasized in the text and on figure y-axes. We plot the L1 norm (maximum) because this is the criterion that is often used in stress based calving laws, like the Nye zero stress. We show the L2 norm because it is (often) a more robust diagnostic of the behavior of the numerical solution.
I would suggest to modify Pa to MPa or kPa. For the x axis, the caption should tell that time step are varying from xx seconds to xx years?
We have modified the axes for stress to be kPa rather than Pa. Figure caption now states that the time step ranges from 1 second to 30 years.
page 5, line 95: not sure the second sentence of part 3.2 should start with "Furthermore"?
A sentence has been added after the equation defining the variable. page 6, line 113: "where the damping coefficient is" Correction has been made.
Responses to Referee comments by Christian Schoof on "Brief communication: Time step dependence (and fixes) in Stokes simulations of calving ice shelves" by Brandon Berg and Jeremy Bassis We thank Christian Schoof for his feedback on this manuscript. Our responses to comments are given below, with original recommendations/points in black and responses in red. When referenced, line numbers refer to the revised manuscript.

Major Recommendations:
We thank the reviewer for his comments and suggestions. We have incorporated most of the reviewers' excellent suggestions in the manuscript. We have vacillated slightly in our preferred terminology between ill-posed and stiff before settling on "unphysical" for reasons that are described in more detail in response to specific reviewer comments.
I would make clear the distinction between the poorly condition Stokes flow problem in Durand et al and the two flavours of actual ill-posedness seen when your boundary conditions permit force and/or torque balance to fail. It won't hurt to allude to the latter, even though I don't imagine many people are trying to solve Stokes flow problems for icebergsyou never know. It would also be reasonable to say that the sea spring mechanism (probably) works well for the second version of the actually ill-posed case, where big departures from equilibrium need never occur. This is a good point. Our emphasis here was really on pointing out that when the geometry evolves rapidly, as is the case for an iceberg calving event, the sea-spring method becomes problematic and can lead to numerical problems. And when these numerical inaccuracies are combined with, for example, stress based calving criteria, there is the possibility of introducing purely numerical calving instabilities. However, under most circumstances, the sea-spring method remains satisfactory. This is better emphasized in lines 69-76.
We have added additional text to clarify the fact that, for our choice of boundary conditions, global force and torque balance are not necessarily satisfied leading to an ill-posed problem (lines 54-56, lines 110-114). However, if we consider fixed velocity (Dirichlet) boundary conditions over a portion of the domain, there is no rigid body motion (translation or rotation) that can be added to the ice shelf. In this case we can still obtain large velocities that are time step dependent using the sea-spring method when the geometry departs significantly from hydrostatic equilibrium over a portion of the domain.
Starting with a geometry that exactly satisfies global and local force/torque balance and then introducing small changes to that geometry can result in large changes to the velocity. This is what we think the reviewer calls "stiff" or poorly conditioned. Small changes in the initial conditions (i.e. position of the ice water interface) lead to large changes in the velocity solution. This can be partly cured by adding, say, a quadratic drag force due to the water, as the reviewer notes. However, for realistic drag coefficients, this still results in exceptionally large velocities. In fact, for configurations that we tested, the velocities can exceed the speed of sound! Because of this and because of the fact that we wish to avoid any confusion between "ill-posed", "ill-conditioned", and "stiff", we have decided to rephrase and call this behavior "unphysical". We believe this captures the numerical issue accurately and avoids introducing additional jargon that glaciologists might not be as familiar with.
The paper says 'However, as we shall show, this assumption is problematic for applications where the ice departs from hydrostatic equilibrium.' I'd circle back to this at some point and point out that things may not be quite so dramatic as to say that the Stokes equations have nothing to say about what happens during these 'inertial' events; they do, but in modified form. I think this will also tie in to the discussion of how to formulate the inertial term in discrete form, see again under 'minor points' below.
As stated above, if we consider fixed velocity (Dirichlet) boundary conditions over a portion of the domain, there is no rigid body motion (translation or rotation) that can be added to the ice shelf and we still obtain large velocities that are time step dependent using the sea-spring method. The large velocities are tied to the bending of the ice that occurs in response to removal of ice at the calving front. As the reviewer notes, this is because the problem is "stiff".
But to simplify our discussion, we have removed equation (8) and the accompanying text regarding the separation of the velocity into viscous and uplift components. Instead, we focus on highlighting the "stiffness" of the problem and how small changes to the ice-ocean boundary location can cause large changes in the solution. In this way, we emphasize the importance of carefully treating the hydrostatic uplift without commenting on the exact nature of the decomposition of viscous and rigid body motion. However, we do add text in the manuscript stating that such a decomposition may be possible (lines 110-114).

Minor Points:
I would probably make a bit clearer how boundary conditions are important in determining whether an actual ill-posedness can occur in the sense of there being no solution to the Stokes flow problem. As the extended discussion above indicates, the partial Dirichlet conditions in the present paper ensure there is no issue of torque or horizontal force imbalance, but it others may run into these issues in their own research, and look to apply the method developed here. Also, as indicated in the second paragraph of this review, there may be other tricks to ensuring force balance.
We have added text to clarify this. In particular, we have noted that adding global constraints on force and torque balance is possible (lines 110-114). We do, however, note that because ice breaks, global force and torque balance would have to be considered on each intact segment and this becomes increasingly challenging to efficiently identify and manage. Hence, our solution of simply including the acceleration directly into the Stokes equations becomes a more practical solution. As noted in our previous response, we have also shifted our terminology to "unphysical" because unphysically large velocities are still possible when drag is included.
The decomposition in equation (8): for a sea spring model, I don't think you can argue that the sea spring term u(x, z)δt simply causes an additive term ∆zuplift/δt… We have streamlined and simplified this section and have eliminated this equation and explanation. We now focus on the "stiff" nature of the problem rather than a specific decomposition into viscous and uplift components.
Notation: there is quite a bit of randomness about which quantities are in boldface and which are not, especially when it comes to tensors (σ versus ε and I?). Make it consistent to please the eye. . .

Notation has been changed so that both vectors and tensors are all in boldface.
The '/2' should probably be inside the square root in the definition of the invariant εe just after equation (4) Error has been fixed.
Writing u(∆t) on the left-hand side of (8) is confusing as u has a well-defined meaning as the continuum solution of the Navier-Stokes problem as a function of (x, z, t), so changing the arguments of that function haphazardly to ∆t is bad form (and actually confused me quite a bit). For starters, the quantity on the left isn't u but a numerical approximation to it, solving a modified problem, so give it a different symbol, and be clear why you are using the ∆t argument on the left (your numerical algorithm thus constructed leads to a solution that turns out to depend on ∆t, whereas you would want it not to be dependent on ∆t. We have eliminated Equation (8) in response to other comments by the reviewer.
The numerical form of the acceleration term in equation (9): this is defensible when you have no rotational degrees of freedom in the rigid body motion, because the advection term for momentum u · ∇u is dominated by r · ∇r (see point 2 above), and in the absence of a rotational degree of freedom, ∇r = 0 so the advection term goes away. As soon as there is rotation, this is no longer true, and you are wel advised to retain the full inertial term ∂u/∂t + u · ∇u. I realize that the present paper does not allow for that possibility, but I think it is worth mentioning.
In our model, we are using an Arbitrary Langrangian Eulerian (ALE) formulation, which updates all mesh coordinates at every time step based on the velocity field. This led us to use the material derivative in Equation (9). But we have added text to clarify the fact that even in a Eulerian reference frame we could neglect the u · ∇u term because the velocity field does not contain a rigid body rotation (lines 117-118). We have also added text explicitly stating we are using an ALE formulation (lines 90-92).
Again, equation (9): I am actually not clear how you imagine you are computing this in a 'Lagrangian frame' to begin with, since you are solving, in discrete terms, an elliptic equation (or a parabolic equation with a backward Euler step, which is the same thing); if you genuinely are using a Lagrangian transformation here, please be explicit and speciifc. In terms of implementation, the reduced acceleration term in equation (9) is my biggest concern, even if I believe it to be leading-order correct (in the Reynolds number, see above) for the verticalmotion-only case discussed in the paper.
We are not quite sure that we understand the reviewer's question. In our implementation, we are solving a parabolic equation with a backward Euler step. The update to the ice geometry is done using a fully Langrangian formulation in which we update the mesh coordinates at every time step. Of course, we do need an initial condition for particle velocities. For the initial condition on velocity to use in the backward Euler step, we choose a uniform velocity field 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 motivated by the experimental design, in which we imagine an ice shelf that is initially perfectly at hydrostatic equilibrium, and thus should have nearly zero vertical velocity before calving. We have added text to section 4 of the manuscript to clarify the precise initial condition on velocity (lines 118-121). Abstract. The buoyancy boundary condition applied to floating portions of ice sheets and glaciers in Stokes models is numerically ill-posed ::::::: requires :::::: special ::::::::::: consideration when the glacier rapidly departs from hydrostatic equilibrium. This manifests in velocity solutions that diverge with decreasing ::::::: boundary :::::::: condition ::: can :::::::: manifest :: in ::::::: velocity ::::: fields ::: that ::: are ::::::::::: unphysically :::: (and ::::::: strongly) ::::::::: dependent ::: on time step size, contaminating diagnostic strain rate and :::::: thereby :::::::::::: contaminating :::::::: diagnostic : stress fields.
Copyright statement. TEXT

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

20
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.   1 1) , where the upper portion has calved away, emulating the "footloose" mechanism proposed by Wagner et al. (2014) where a :: an aerial portion of the calving front first detaches.

35
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 ::::::::::::: 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.

Boundary Conditions
To illustrate an example where the Stokes flow problem becomes ill-posed :::::: departs ::::: from ::::::::: hydrostatic :::::::::: equilibrium, we consider a two-dimensional floating ice shelf (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 ) :: and :: n :: is ::: the ::::::: normal ::::: vector ::::: along ::: Γ 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 50 ice-ocean interface cannot be balanced by internal stresses. In this case, : there is no :::::: unique : solution and vertical velocities are singular. In reality, of course, the ice will quickly re-adjust to hydrostatic equilibrium through rapid :: as :: a :::::::::: consequence ::: of buoyant uplift through the (nearly) inviscid ocean.
Inspecting Eq. (??) shows that in the limit of large ∆t, the uplift term becomes small compared to the viscous velocity. Thus, 85 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 90 to be in exact hydrostatic equilibriumwith . ::: We ::: set ::: the : 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(Rignot et al., , 2011Mouginot 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 95 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 (Alnaes et al., 2015), an open source finite element solver with a Python interface 100 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 a unstructured mesh with uniform grid spacing of 10 m near the calved portion of the domain and grid spacing of 40 m elsewhere.
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 outside of the Stokes range.
In addition to the divergence ::::::: behavior : of the L 2 norm, we also examine the maximum effective :::: value ::: of ::: the ::::::::: maximum 115 shear stress and greatest principle ::::::: principal : 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.

135
:::::::: Restoring ::: the :::::: inertial ::::: term effectively introduces a Newtonian damping term on the entire body of the glacier where the damping coefficient :: is 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 unstable) for long time steps. Therefore, we propose to use both damping terms so that the system of equations is numerically accurate for small :: all time step sizesand the velocity converges to the viscous limit for large time step 140 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. : .

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 diverges : is :::::::::: unphysical for small time steps when the domain departs from hydrostatic equilibrium. For model applications where changes in the domain are only due to viscous flow, the 165 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 170 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 results in a physically consistent solution ::::::: provides :::::::: consistent :::::::: solutions ::: for :: all :::: time ::::: steps.
Author contributions. BB identified the numerical issue with guidance from JNB. BB and JNB developed the proposed solution to the numerical issue. BB prepared the manuscript with contributions from JNB.