Numerical simulation of glacier terminus evolution using the dual action principle for momentum balance

Abstract

The momentum conservation equation for glacier flow can be described through minimization of an action functional. Several software packages for glacier flow modeling use this action principle in the design of numerical solution procedures. We derive here an equivalent dual action principle for the shallow stream approximation and implement this model using the finite element method. The key feature of the dual action is that the flow law and friction law are both inverted, which changes the character of the non-linearities. This altered character makes it possible to implement numerical solvers for the dual form that work even when the ice thickness or strain rate are exactly equal to zero. Solvers for the primal form typically fail on such input data and require regularization of the problem. This robustness makes it possible to implement iceberg calving in a simple way: the modeler sets the ice thickness to zero in the desired area. We provide several demonstrations and a reference implementation.

Introduction

On space and timescales >100 m and 1 d, glaciers flow like a viscous, incompressible fluid with a power-law rheology (Greve and Blatter, 2009). Ice flow is slow enough that the fluid inertia is negligible compared to viscous and gravitational forces, that is, the flow occurs at very low Reynolds and Froude numbers. There are multiple equivalent ways of expressing the momentum balance equations: a conservation law, a variational form, a partial differential equation (PDE). Each of these forms is best suited to a different type of numerical method. The momentum balance equation for low-Reynolds number viscous fluid flow can also be derived as the optimality conditions for the velocity to be the critical point of a certain action functional (Dukowicz and others, 2010). The action functional has units of energy per unit time and can be interpreted as the rate of dissipation of thermodynamic free energy (Edelen, 1972). For many problems – low-Reynolds number flow, heat conduction, saturated groundwater flow, steady elasticity – the action is a convex functional of the unknown field.

The existence of an action principle is a special property of a very restricted class of differential equations. Action principles are not just of theoretical interest – we can use them to design faster, more robust numerical solvers. First, a convex action principle implies that the second derivative is symmetric and positive-definite. These properties guarantee convergence for Newton-type algorithms. They also mean that we can use specialized methods, such as Cholesky factorization or the conjugate gradient method, to solve the linear systems of equations for the search direction in each step (Nocedal and Wright, 2006). These methods are not applicable to more general classes of linear systems. Second, the action principle offers a way to measure how well an approximate solution matches the true solution and it is distinct from, say, the square norm of the residual. The theory of convex optimization then provides us with a way to measure how close we are to convergence using only the current solution guess and search direction by evaluating the Newton decrement. In Shapero and others (2021), we showed how to use this theory to design physics-based convergence criteria.

This work follows in the footsteps of Dukowicz and others (2010) in studying action principles for glacier flow. Our main contribution is the derivation of an alternative dual action principle, distinct from that presented in Dukowicz and others (2010), from which the momentum conservation equations can be derived. The most important feature is that the dual action principle has favorable numerical properties for shear-thinning flows such as glacier dynamics. Solving the primal form of the problem requires regularization around zero strain rate, velocity, and thickness in order to smooth away infinite values. This regularization makes the momentum balance problem solvable, but it remains poorly conditioned and introduces other non-physical artifacts. The dual problem requires no regularization. We have implemented solvers for this dual form that still converge even when the thickness and strain rate are zero. As a consequence, we were able to simulate iceberg calving by setting the ice thickness to zero in part of the glacier. We illustrate these advantages in the final section with a numerical implementation and several demonstrations.

The main advantage of the approach we propose here is that it offers a new way to handle ice-free regions. Several strategies already exist in the literature on numerical ice flow modeling for handling ice-free regions. One can set a minimum ice thickness, which regularizes away the problem but introduces mass-balance errors. The BISICLES model uses a finite volume discretization, special handling of the terminus in assembling the stiffness matrix, a regularized Picard-type linearization and an artificial friction in ice-free areas (Cornford and others, 2013). The Ice Sheet and Sea-Level System Model uses the level-set method (Osher and Sethian, 1988; Bondzio and others, 2016). This approach introduces an additional scalar field which evolves according to a certain differential equation. The zero contour of this scalar field represents the glacier terminus. One can then ‘turn off’ the physics in the ice-free region where the momentum balance equation ceases to be well-posed. Finally, the Elmer/Ice model has used both (1) direct remeshing of the 3-D geometry (Todd and others, 2018), so ice-free areas are not included in the computational domain at all, and (2) coupling to a discrete-element model (Benn and others, 2017). These approaches are effective but come with their own drawbacks and implementation challenges. For example, using the level-set method requires solving a challenging, non-linear hyperbolic problem, the eikonal equation. The remeshing approach taken in Elmer/ice, on the other hand, requires projecting the solution between different computational meshes. The dual form that we describe here has its own challenges but we claim that these are easier to overcome than those of existing approaches.

In the following, we will assume familiarity with (1) the PDEs describing glacier flow, (2) variational calculus and the derivation of the Euler–Lagrange equations of a generic functional and (3) convex analysis and convex duality theory. For background reading, we refer the reader to Greve and Blatter (2009) for glacier dynamics, Weinstock (1974) for variational calculus and Boyd and Vandenberghe (2004) for convex optimization.

Theory

The shallow stream equations

Here we review the differential equations that are commonly used to describe glacier flow. In the next section, we will show how the momentum balance equation has a minimization principle. We will focus exclusively on the shallow stream approximation (SSA), which is commonly applied to model fast-flowing outlet glaciers and ice streams. The SSA model is derived by (1) expanding the Stokes equations in the aspect ratio and taking only the lowest-order terms, and (2) depth-averaging the equations, which assumes that the horizontal velocity varies much more in the longitudinal directions than with depth (Greve and Blatter, 2009). For a complete list of all symbols and units used in the following, see Table 1.

Table 1. Variable, symbol, physical units and tensor rank – 1 for vectors, 2 for matrices, etc.

The equations of motion are solved in a 2-D domain Ω. The main unknown to be solved for in the SSA is the depth-averaged ice velocity u. Some important intermediate quantities are the basal shear stress τ and the membrane stress tensor M, a rank-2 tensor with units of stress that results from applying the low-aspect ratio assumption to the full 3-D deviatoric stress tensor. The inputs to the problem include the thickness h, surface elevation s and fluidity factor A in Glen's flow law. The SSA momentum conservation equation is

(1)$$\nabla\cdot hM + \tau - \rho_{\rm I} gh\nabla s = 0,\; $$

where ρ _I is the ice density, and g the gravitational acceleration. In addition to Eqn (1), we need to know a constitutive relation between the membrane stress tensor and the depth-averaged strain rate tensor:

(2)$$\dot\varepsilon = {1\over 2}\left(\nabla u + \nabla u^\top\right).$$

In order to simplify the notation later, we introduce the dimensionless rank-4 tensor ${\scr C}$ defined by

(3)$${\scr C}\dot\varepsilon = {\dot\varepsilon + {\rm tr}( \dot\varepsilon) I\over 2}.$$

The tensor ${\scr C}$ plays a similar role to the elasticity tensor in linear elasticity. Moreover, we define the norm of a rank-2 tensor with respect to ${\scr C}$ as

(4)$$\vert \dot\varepsilon\vert _{{\scr C}}^2 = \dot\varepsilon\, \colon \, {\scr C}\dot\varepsilon.$$

Alternatively, in index notation, the square norm is $\vert \dot \varepsilon \vert _{\mathscr C}^2 = {\scr C}_{ijkl}\dot \varepsilon _{ij}\dot \varepsilon _{kl}$. With these notational conveniences in hand, the Glen flow law states that the membrane stress and strain rate are related by a power law:

(5)$$M = 2A^{- {1}/{n} }\vert \dot\varepsilon\vert _{\scr C}^{{1}/{n} - 1}{\scr C}\dot\varepsilon$$

where A is the depth-averaged fluidity coefficient and n ≈ 3 is the Glen flow law exponent.

Next, we need to provide some kind of sliding relation. We will assume a generalized power law with some exponent m, that is,

(6)$$\tau = -C\vert u\vert ^{{1}/{m} - 1}u.$$

Weertman sliding has m = n, while perfectly plastic sliding has m = ∞. Recent research suggests alternative forms that transition between Weertman-type sliding at low speeds and perfectly plastic sliding at higher speeds (Minchew and Joughin, 2020). For illustrative purposes Eqn (6) is sufficient, and we will describe how to incorporate alternatives in the discussion.

Finally, we need to supply a set of boundary conditions for the problem to be well-posed. Along the part of the boundary where ice is flowing into the domain, we assume that the ice velocity is known from observations; this is a Dirichlet boundary condition. At the glacier terminus, which we will denote by Γ, the membrane stresses at the cliff face are balanced by pressures from any proglacial water body. This is a Neumann boundary condition:

(7)$$ - \! hM\cdot\nu = {1\over 2}\left(\rho_{\rm I}gh^2 - \rho_{\rm W}gh_{\rm W}^2\right)\nu,\; $$

where h _W is the water depth and ν is the unit outward-pointing normal vector to the terminus.

We can then combine Eqns (1)–(6) and the boundary condition (7) to arrive at a second-order, non-linear elliptic system of differential equations for u.

For modeling a floating ice shelf, the friction term in Eqn (1) is zero. Moreover, assuming the ice is in hydrostatic equilibrium allows us to write the surface elevation in terms of the thickness:

(8)$$s = ( 1 - \rho_{\rm I} / \rho_{\rm W}) h.$$

As a result, the momentum conservation equation reduces to

(9)$$\nabla\cdot hM - {1\over 2}\varrho g\nabla h^2 = 0,\; $$

where ϱ = (1 − ρ _I/ρ _W)ρ _I is the reduced density of ice over sea water.

Marine ice sheets flow from the continent and into the ocean, where they go afloat. If we are to model a marine ice sheet with the SSA, we must distinguish between a grounded and a floating region. This results in a free boundary problem where ice goes afloat once condition (8) is satisfied. The boundary separating grounded from floating ice is known as the grounding line x _g. The possibility of a marine ice-sheet instability that could dramatically increase the rate of discharge of ice into the ocean has led to a substantial amount of research into grounding line dynamics (Schoof, 2007; Durand and others, 2009; Favier and others, 2012).

To complete our description of the dynamics, the thickness evolves in time according to the following depth-averaged mass conservation equation:

(10)$${\partial h\over \partial t} + \nabla \cdot hu = \dot a - \dot m,\; $$

where $\dot a$ is the rate of ice accumulation and $\dot m$ of melting or ablation. We assume that the ice thickness and velocity along the inflow boundary are known.

Primal action principles

The main idea behind action or minimization principles is that some PDEs really express the fact that their solutions are extrema of a given action functional. The momentum balance equation of glacier flow is one such PDE. Many publications in the glacier flow modeling literature have explored the advantages of using action principles to describe the momentum balance (Bassis, 2010; Dukowicz and others, 2010; Brinkerhoff and Johnson, 2013; Shapero and others, 2021). Here we briefly review these concepts as they pertain to the SSA momentum balance.

Given a particular action functional, the PDE that expresses the condition that a field is an extremum can be calculated mechanically in terms of the integrand. This PDE is called the Euler–Lagrange equation for the functional. We will not repeat it here but see Weinstock (1974). On the other hand, if we are given a PDE, it may or may not have an action functional at all. In other words, being the Euler–Lagrange equations for some action functional is a special property of only a restricted class of PDEs. There is no rote procedure to determine what the action functional might be, but in attempting to construct one, the main mathematical hurdle to overcome is computing the anti-derivatives of certain terms in the differential equation.

The constitutive relation (5) for M can be expressed as the derivative of a certain scalar quantity. In index notation,

(11)$$M_{ij} = {\partial\over \partial\dot\varepsilon_{ij}}\left({2n\over n + 1}A^{- {1}/{n} }\vert \dot\varepsilon\vert _{\scr C}^{{1}/{n} + 1}\right)$$

if we think of the strain rate tensor as an independent variable and briefly forget that it is the symmetrized velocity gradient. Likewise, for the sliding relation (6), we can observe that

(12)$$\tau_i = -{\partial\over \partial u_i}\left({m\over m + 1}C\vert u\vert ^{{1}/{m} + 1}\right).$$

Using Eqns (11) and (12), one can show that the SSA momentum balance are the Euler–Lagrange equations for the following action functional (Dukowicz and others, 2010):

(13)$$\eqalign{J( u) = & \int_\Omega\left\{{2n\over n + 1}hA^{-{1}/{n}}\vert \dot\varepsilon\vert _{{\scr C}}^{{1}/{n} + 1} + {m\over m + 1}C\vert u\vert ^{{1}/{m} + 1}\right. \cr & \left. \vphantom{{2n\over n + 1}} \qquad + \rho_{\rm I} gh\nabla s\cdot u\right\}\, {\rm d}x + {1\over 2}\int_\Gamma\left(\rho_{\rm I}gh^2 - \rho_{\rm W}gh_{\rm W}^2\right)u\cdot\nu\; {\rm d}\gamma}$$

Note that the action has units of energy per unit time, or power. The final summand of Eqn (13) enforces the Neumann boundary condition at the ice terminus. The Dirichlet boundary condition along the ice inflow, on the other hand, has to instead be enforced by eliminating coefficients from the resulting non-linear system.

A mechanical computation of the second derivative of J shows that this functional is convex. The theory of non-equilibrium thermodynamics tells us that J represents the rate of dissipation of thermodynamic free energy (Edelen, 1972). Once again, the essential point is that finding a minimizer u of the functional J defined above is equivalent to finding a solution of the SSA differential equation.

Dual action principles

Action principles have appeared in glaciology before, but dual forms have not. The dual form of a problem is a distinct but equivalent expression of the same underlying physics. We will show in the following that, for the particular case of the SSA momentum balance, the dual form has some better numerical properties that make it worth investigating.

The dual form can be understood at two levels. First, we can show what the dual action functional is without describing how we came up with it. A reader who knows some variational calculus can derive the Euler–Lagrange equations for this functional and verify that the resulting equation set is equivalent to the primal form of the SSA. We present an exposition at this level below. On the other hand, this approach can feel like pulling a rabbit out of a hat; it does not answer how we arrived at the dual form in the first place or how we might derive the dual forms of other problems. This level requires some knowledge of general convex duality theory. We assume that few glaciologists will be interested in this level of detail and refer instead to sections 9.7 and 9.8 of Attouch and others (2014).

Equations (1)–(7) can be combined into a second-order differential equation for the velocity u. Solving this differential equation is equivalent to finding a minimizer of the functional J defined in Eqn (13). In deriving the differential equation, we used the constitutive relation and the sliding law (Eqns (5) and (6)) to eliminate the membrane stress M and basal stress τ.

We could instead keep these additional equations and unknowns. Instead of solving a second-order equation for u, we would then have an equivalent first-order system of equations for u, M and τ. One way to motivate the dual problem is to ask: is there an optimization problem that is equivalent to this first-order system? As we will show below, there is, but it has some different characteristics from the primal problem. One of the key features of dual forms in general is that they invert constitutive relations. Conventionally, we write the membrane stress tensor as a function of the strain rate tensor, and the basal shear stress as a function of the sliding velocity. The dual form inverts these relations: the strain rate tensor becomes a function of the membrane stress tensor, and the sliding velocity a function of the basal shear stress.

In Eqn (3), we define the rank-4 tensor ${\scr C}$ as a convenience for writing down how the membrane stress is a function of the strain rate. We can show explicitly that the inverse ${\scr A}$ of this tensor ${\scr C}$ is

(14)$${\scr A}M = {M - \frac{1}{d + 1} {\rm tr}( M) I\over 2}$$

where d = 2 is the space dimension. The tensor ${\scr A}$ plays an analogous role to the compliance tensor in linear elasticity. With this definition in hand, we can invert Eqn (5) to get $\dot \varepsilon = 2A\vert M\vert _{{\scr A}}^{n - 1}{\scr A}M$. The sliding law is much easier to invert: u = −C ^−m|τ|^m−1τ. We can then define a slipperiness coefficient K = C ^−m.

With these preparations in place, the dual form is

(15)$$\eqalign{ L( u,\; M,\; \tau) = & \int_\Omega \Bigg\{{2\over n + 1}hA\vert M\vert _{\mathscr A}^{n + 1} + {1\over m + 1}K\vert \tau\vert ^{m + 1} \cr & - hM\, \colon \, \dot\varepsilon( u) + \tau\cdot u - \rho_{\rm I}gh\nabla s\cdot u\Bigg\}{\rm d}x \cr & - {1\over 2}\int_\Gamma\left(\rho_{\rm I}gh^2 - \rho_{\rm W}gh_{\rm W}^2\right)u\cdot\nu\; {\rm d}\gamma}$$

We can then evaluate the variational derivatives of L with respect to u, M and τ, require that these derivatives are all zero, and show that the resulting equations are equivalent to the primal form of the problem. First, if we require that the variational derivative of L with respect to M along any perturbation N is zero, we find that

(16)$$\left\langle{\partial L\over \partial M},\; N\right\rangle = \int_\Omega\left\{2hA\vert M\vert _{{\scr A}}^{n - 1}{\scr A}M - h\dot\varepsilon( u) \right\}\, \colon \, N\, {\rm d}x = 0.$$

This equation is the variational form of the inverse of the constitutive relation (Eqn (5)). Next, taking the variational derivative of L with respect to τ along some perturbation σ, we get

(17)$$\left\langle{\partial L\over \partial\tau},\; \sigma\right\rangle = \int_\Omega\left\{K\vert \tau\vert ^{m - 1}\tau + u\right\}\cdot\sigma\, {\rm d}x = 0$$

which is the inverse of the sliding law of Eqn (6). Finally, let v be any perturbation to the velocity field such that v = 0 along the inflow boundary. Taking the variational derivative of L with respect to u along v gives

(18)$$\eqalign{\left\langle{\partial L\over \partial u},\; v\right\rangle & = \int_\Omega\left\{-hM\, \colon \, \dot\varepsilon( v) + \tau\cdot v - \rho_{\rm I}gh\nabla s\cdot v\right\}{\rm d}x \cr & \quad -{1\over 2}\int_{\Gamma}\left(\rho_{\rm I}gh^2 - \rho_{\rm W}gh_{\rm W}^2\right)v\cdot\nu\; {\rm d}\gamma}$$

which is not exactly what we want. We can use integration by parts to push the symmetric gradient of v over onto hM in the first term however:

(19)$$\eqalign{\ldots & = \int_\Omega\left\{\nabla\cdot hM + \tau - \rho_{\rm I}gh\nabla s\right\}\cdot v\; {\rm d}x \cr & \quad + \int_\Gamma\left(hM\cdot\nu - {1\over 2}\left(\rho_{\rm I}gh^2 - \rho_{\rm W}gh_{\rm W}^2\right)\nu\right)\cdot v\; {\rm d}x.}$$

If we require that this is equal to 0 for all perturbation fields v, we recover both conservation of membrane stress (Eqn (1)) and the boundary condition (Eqn (7)). As with the primal form, we have to enforce the inflow boundary condition by eliminating degrees of freedom. In effect, finding a critical point of the dual action functional is a constrained optimization problem; the velocity field acts like a Lagrange multiplier that enforces the constraint of stress conservation.

The important feature of the dual formulation is that the nature of the non-linearity has changed. In the primal action shown in Eqn (13), the non-linearity consists of the strain rate raised to the power 1/n + 1. Since n > 1, the non-linearity in the primal form has an infinite singularity in its second derivative around any velocity field with zero strain rate. In the dual form, however, the non-linearity consists of the stress tensor raised to the power n + 1. Around zero stress, the second derivative of the action with respect to the stress is zero instead of infinity; see Figure 1.

Figure 1. Viscous part P of the action is shown in blue and its second derivative in orange, in (a) for the primal problem as a function of the strain rate $\dot \varepsilon$ and in (b) for the dual problem as a function of the stress τ. The second derivative of the viscous dissipation goes to infinity near zero strain rate for the primal problem, but to zero near zero stress for the dual problem.

Demonstrations

Here we will describe several computational experiments for evaluating the dual form of SSA and our implementation of it. First, we will conduct a verification exercise in order to make sure that we correctly implemented the dual form of SSA. This demonstration is to give some assurance that we are solving the equations right. Next, we will conduct two numerical exercises to compare how well the dual form works compared to the primal form on problems in simple geometries. Finally, we will conduct two more experiments to show off the use of the dual form on realistic glacier geometries.

Verification on solvable test cases

The verification exercises we use are taken from those used to test the implementation of the primal form of SSA in the icepack package (Shapero and others, 2021). We compare numerical results on a sequence of grids to exactly solvable instances of SSA and check that the results converge with the expected order of accuracy. Finite element theory predicts that the L ²-norm difference between the exact solution and the solutions obtained using CG(k) finite elements is ${\scr O}( \delta x^{k + 1})$ where δx is the mesh spacing. If the slope in a log–log fit of error against mesh spacing deviates significantly from k + 1, this would indicate some mis-specification of the problem or bug in the solver.

The first test is to use the exact solution for the velocity of a floating ice shelf with thickness

(20)$$h = h_0 - \delta h \cdot x / L_x$$

in a domain of length L _x = 20 km. With a constant value of the fluidity coefficient A, the velocity in the x direction is

(21)$$u_x = u_0 + L_x A \left({\varrho g h_0\over 4}\right)^n\left(1 - \left(1 - {\delta h \cdot x\over h_0\cdot L_x}\right)^{n + 1}\right)$$

where ϱ = ρ _I(1 − ρ _I/ρ _W). We use a 2-D domain in order to make sure that the numerical solution, like the exact solution, has no variation in the y direction.

The second test case uses the same geometry but adds basal friction and assumes the ice thickness is above flotation. Solving the resulting boundary value problem analytically in the presence of friction is now much more difficult. Instead, we used the method of manufactured solutions – we picked the ice velocity, thickness and surface elevation, and generated a friction coefficient that would make this velocity an exact solution. To generate this friction coefficient we used the computer algebra system sympy (Meurer and others, 2017).

Comparison with primal form on slab glacier

As a more challenging test in a flowline setting, we consider a slab of ice of constant thickness flowing down an inclined slope into the ocean, where it goes afloat at the grounding line (Fig. 2). We compute steady states under this configuration using the primal and dual forms described above. We set the ice thickness at x = 0 to h = 500 m and the bedrock angle to $\alpha = 1^\circ$. The bed is given by the expression

(22)$$b( x) = 1500\, {\rm m} - x\tan{\alpha}.$$

For the material parameters, we set n = 3, A = 10⁻²⁴ Pa⁻³ s⁻¹, C = 10⁻⁶ Pa m^−1/3 s^1/3. Instead of setting inflow boundary conditions at x = 0, we enforce a zero extensional stress condition M = 0. This condition allows the ice sheet to tend to the uniform thickness slab solution for upstream of the grounding line. For the uniform slab, the extensional stresses are equal to zero and the frictional stresses at the base of the ice sheet balance the gravitational forces. As a result, the horizontal velocity for the uniform thickness slab is given by

(23)$$u = \left({\rho_{\rm I} g h \tan \alpha \over C}\right)^n.$$

A consequence of enforcing the slab solution at the left boundary is that, as we move upstream away from the grounding line, the strain rate tends to zero and one must regularize the primal formulation of the momentum balance equation.

Figure 2. Setup for modeling a slab of ice on an inclined bed flowing into the ocean. At x = 0 we enforce a thickness h = 500 m in order to approach a parallel slab of ice far upstream of the grounding line. The dotted line is sea level.

For this problem, we not only solve for the velocity u or the velocity–stress pair (u,  M), but also for the thickness h and the grounding line position x _g. We therefore complement the momentum balance equations with the mass balance Eqn (10) and the flotation condition (8), effectively yielding a free boundary problem.

Comparison against primal form on gibbous ice shelf

The next comparison exercise uses the synthetic ‘gibbous’ ice shelf test case from §5.3 of Shapero and others (2021). The domain consists of the intersection of two circles of different radii chosen to roughly mimic the overall size of Larsen C.

We first run a spin-up of the system to steady state using the coupled mass and momentum balance equations with both the primal and dual forms. We check that the velocity obtained by solving the dual form is within discretization error of the velocity obtained with the primal form, which offers an additional degree of verification that we are solving the equations right. We additionally evaluate the total wall clock time to run this experiment using both the primal and dual forms. The primal problem using CG(1) elements for the velocity has two unknowns for each vertex of the mesh. The dual problem using CG(1) × DG(0) elements for the velocity and stress has an additional three unknowns per triangle. The Euler formula (#vertices − #edges + #triangles ≈2) implies that there are approximately twice as many triangles as there are vertices. Consequently there are ~4× as many degrees of freedom when solving the dual problem as there are for the primal problem. Assuming naively that the time to solution scales linearly with the number of unknowns, we would then expect that solving the dual problem is 4× as expensive as solving the primal problem.

As a third and final phase of this experiment, we run the same simulation, but every 24 years we set the ice thickness to 0 in a prescribed region near the terminus. This forcing mimics the effect of a large iceberg calving event. Our prescribed evolution of the terminus is not a realistic representation of how calving works. Instead, we aim only to stress test the solver in order to see if it can handle regions of zero thickness.

Demonstration on calving of the Larsen C Ice Shelf

To test the dual form of SSA on a realistic problem, we will simulate the evolution of the Larsen C Ice Shelf from a nominal start date of 2015 for 40 years, including the calving of Iceberg A-68 in 2017 (Larour and others, 2021). This experiment uses the observed calving front positions from satellite imagery to set the terminus positions at the start of the simulating and after the calving event. The goal is not to implement a calving law as such. In all, the experiment proceeds in several steps:

1. (1) Estimate the fluidity field A from remote-sensing measurements of the thickness and velocity. This step uses the primal form of the momentum balance equation from icepack.

2. (2) Extrapolate the ice thickness and velocity onto a larger spatial domain, making the ice thickness 0 in ice-free areas.

3. (3) Run the simulation using the mass and dual momentum balance from the start date of 2015 until the calving event in 2017.

4. (4) Digitize the terminus position immediately after the calving event by hand and use the digitized terminus position to define an ice mask.

5. (5) Using this mask, set the ice thickness to zero over the spatial extent of the calved area.

6. (6) Run the simulation for 40 years after the calving event to see how the terminus advances again.

Demonstration on Kangerlussuaq Glacier

Our final test case is simulating Kangerlussuaq Glacier, a grounded outlet glacier on the east coast of Greenland. Kangerlussuaq is one of the top three contributors to the total discharge from Greenland (Enderlin and others, 2014; Mouginot and others, 2019). The purpose of this exercise is to demonstrate that we can simulate the evolution of a marine-terminating glacier, including the seasonal advance and retreat of the terminus in response to ocean-induced frontal ablation in summer, using the dual form. We do not aim to reproduce the exact calving history.

The exercise proceeds in several steps, similar to our approach for Larsen C:

1. (1) Estimate the slipperiness (the coefficient K in the sliding law u|_z=b = −K|τ _b|ⁿ⁻¹τ) from remote-sensing measurements of the ice thickness, surface elevation and velocity. This step uses the primal form of the momentum balance equation from icepack.

2. (2) Extrapolate the thickness, surface elevation, velocity and friction coefficient onto a large spatial domain that extends further down Kangerlussuaq Fjord.

3. (3) Run the simulation using the mass and dual momentum balance equations for 1 year in order to propagate out any initial transients. This stage uses only surface mass balance (SMB) and thus permits the glacier to advance down the fjord.

4. (4) Turn on a time-periodic ablation field near the terminus in order to represent the effects of summer melt and calving and run the simulation for a further 4 years. This ablation field forces the terminus to advance and retreat.

To initialize the simulation, we use version 3 of the BedMachine Greenland dataset for ice thickness and surface elevation (Morlighem and others, 2017) and the MEaSUREs annual velocity mosaic from 2015 to 2016 (Joughin and others, 2010) to infer the basal friction. To force the mass conservation Eqn (10), we need to provide a SMB field $\dot a$ and a melt rate $\dot m$.

We use an SMB field that varies linearly with elevation:

(24)$$\dot a \approx a_0 + {\delta a\over \delta s}\cdot s$$

where a ₀ is the SMB at sea level and δa/δs is the SMB lapse rate. To fit the parameters a ₀ and δa/δs, we used output from 2006 to 2021 of version 3.12 of the Modèle Atmosphérique Régional (Fettweis and others, 2020). This regional climate model has been tested extensively for the polar regions and for Greenland. The fit had r ² = 0.91, so a substantial fraction of the variance is explainable by surface elevation alone.

To set the melt rate $\dot m$, we first create a smoothed ice mask μ. The smoothed mask is required to be equal to 1 on the inflow boundary 0 on the outflow boundary, and have 0 normal derivative along the side walls. We then compute μ as the minimizer of

(25)$$J( \mu) = {1\over 2}\int_\Omega\left(( \mu - {\bf 1}_{\{ h > 0\} }) ^2 + \alpha^2\vert \nabla\mu\vert ^2\right){\rm d}x$$

where α is some smoothing length, and 1_{h>0}(x) is equal to 1 if h(x) > 0 and 0 if h(x) = 0. Here we choose α to be 1 km, so the mask field rapidly approaches 1 within roughly one ice thickness of the terminus. The mask field μ is recalculated in every time step. Finally, we set the melt rate at time t as

(26)$$\dot m = m_0( 1 - \mu) \min\{ 0,\; \cos( 2\pi t) \} $$

where m ₀ is a maximum melt rate that we have to choose. Although we do not employ the level-set method here directly, the approach outlined above is similar to using a level-set method.

The purpose of this exercise is to demonstrate that our solver for the dual form can simulate advance and retreat of a grounded tidewater glacier in response to melt forcing at the terminus. Again, our goal is not to validate a particular calving law.

Results

We implemented a solver for the dual form of the SSA using the Firedrake package (Ham and others, 2023). For more information on discretizing the dual form using finite elements and for strategies to solve the resulting finite-dimensional optimization problem, see the Appendix.

Verification on solvable test cases

We tested meshes with between 16 and 256 cells to a side and we used both CG(1) × DG(0) and CG(2) × DG(1) finite element pairs for the velocity and membrane stress, where CG and DG denote respectively continuous and discontinuous Galerkin elements. The relative errors in the L ² norm have the expected asymptotic convergence rates of ${\scr O}( \delta x^2)$ for linear velocity elements and ${\scr O}( \delta x^3)$ for quadratic in both the ice shelf and ice stream test cases; see Figure 3.

Figure 3. Relative L ²-norm errors for approximate solutions to the analytical ice shelf (a) and ice stream (b) test cases using our newly-developed solver for the dual form of SSA. The points show the error values from each experiment, the lines show a log–log fit of the errors against mesh size. The convergence rates were obtained from this log–log fit.

While finite element theory can predict the asymptotic convergence rates, it does not immediately give estimates of what the constant prefactor should be except in the most trivial of linear problems. The constants can only be evaluated empirically. In particular, the theory predicts that quadratic elements converge faster asymptotically than linear elements, but it cannot tell us how many cells per side are necessary for each to achieve the same accuracy. Figure 3 shows that a numerical solution obtained with only 16 cells per side and quadratic elements is roughly as accurate as a solution with 256 cells per side and linear elements.

Comparison with primal form on slab glacier

We solved the free boundary problem with a primal method that seeks the velocity and thickness in CG(2) × CG(2), and with a dual method that computes the velocity, membrane stress and thickness in the space CG(2) × DG(1) × CG(2). For the primal method, we need to include a regularization parameter $\epsilon$ in order to prevent singularities in the constitutive relation. For this exercise we solve a 1-D form of the equation, so the relevant term in the variational form of the momentum balance equation is

(27)$${\langle F( u) ,\; v\rangle = \int_\Omega \left\{2hA^{-1/n}\vert \partial_xu^2 + \epsilon^2\vert ^{( n - 2) /2}\partial_xu\cdot\partial_xv + \cdots\right\}{\rm d}x}$$

We consider a sequence of regularization parameters $\epsilon$ between 1 and 10⁻¹² a⁻¹. The results for the grounding line position are displayed in Table 2. The discrete problem is solved with Newton's method, and the initial guess for the values of the ice velocity, the ice thickness and the extensional stress are set equal to the slab solution, such that h = 500 m, u is equal to Eqn (23), and M = 0. The initial guess for the grounding line position is set to the point where the flotation condition (8) holds for the constant thickness slab. We plot the values of the relative Newton residual in Figure 4. The solution obtained with the dual form is as accurate as the primal solution using the lowest value of regularization. Moreover, the rate of convergence of the Newton solver for the primal formulation quickly decreases for low values of $\epsilon$. For values of $\epsilon$ equal to or lower than 10⁻¹⁴ a⁻¹, the relative Newton residual no longer reaches the minimum tolerance of 10⁻⁸ that we set for this problem.

Table 2. Results for the slab of ice flowing into the ocean

Values of the steady state grounding line position x _g and thickness at the grounding line for computations with the primal formulation with varying regularization parameters $\epsilon$ and with the dual formulation. We also present the number of Newton iterations required to converge.

Figure 4. Results for the slab of ice flowing into the ocean. Norm of the relative Newton residual for computations with the primal formulation with varying regularization parameters $\epsilon$ and with the dual formulation.

Figure 4 shows that using a larger value of the regularization parameter reduces the number of iterations needed to achieve convergence. However, using more regularization also increases the misfit between the computed velocity and the true velocity. The dual form makes no such compromise in accuracy but the solver still retains a high degree of efficiency.

Gibbous ice shelf

For the spin-up phase of the experiment, we did an initial run for 400 years on a mesh with a 5 km resolution, at which point the system is close to steady state. We then projected these fields to a finer mesh with a resolution of 2 km and use them as the initial state for a further 400 years of spin-up. The results are shown in Figures 5a–c and are identical to those obtained from the primal form of the problem up to discretization error.

Figure 5. Thickness (a), velocity (b) and magnitude of the membrane stress tensor (c) in steady state, and the thickness (d), magnitude of the velocity change (e) and magnitude of the stress change (f) immediately after the calving event. We remove a semi-circular segment from the end of the shelf with a prescribed center and radius.

When we used the spin-up phase of the experiment as a benchmark to measure the performance of the dual and primal solvers, we found that the dual problem required between 2.5 × and 2.7 × as much time. These results were consistent across different mesh resolutions and when run several times on multiple machines. Since the dual problem has 4 × as many unknowns, the added cost that we found experimentally is less than what we would expect if we naively assumed that cost is proportional to the number of degrees of freedom.

In the calving phase of the experiment, our solver for the dual problem still worked in ice-free areas. This feature offers the possibility of implementing physically-based calving models in a simple way. Figure 6 shows the evolution of the volume of ice in the shelf over the two spin-up phases and the calving phase. The 24 year recurrence interval is not enough time for the ice to advance back to the original edge of the computational domain. Using a longer interval would allow the calving terminus to advance back to its original position. Figures 5d–f show the thickness of the ice shelf immediately after the calving event, the magnitude of the change in speed, and the magnitude of the change in stress. In particular, the stress field shows a discontinuity at the new calving front as expected.

Figure 6. Total volume of ice in the shelf over time. The different spin-up and experimental phases are labeled. Note how the finer spin-up equilibrates to a smaller ice volume than the coarser spin-up.

We repeated this phase of the experiment using a comparable solver for the primal problem. When we use no minimum thickness at all, the solver for the primal form diverges as soon as there are any ice-free areas. To remedy this problem, we clamped the thickness from below at 1 mm. Figure 7 shows the number of Newton iterations necessary to obtain the desired level of convergence through two calving events using both the primal and dual forms. In each case, the number of iterations goes up after a calving event. As the system relaxes back, the number of iterations decreases again. The number of iterations required for the dual form is in general slightly greater.

Figure 7. Number of Newton iterations to compute the ice velocity at each step of the calving phase of the experiment using the primal form with the thickness clamped from below and using the dual form. Calving occurs every 24 years.

Larsen C Ice Shelf

Our solver for the dual form of SSA was successfully able to compute velocity and stress fields on this realistic test case even in ice-free areas, enabling effective simulation of calving events. The simulated terminus positions of Larsen C at the start of the simulation, immediately after the calving event and at the end are shown in Figure 8. The model can effectively handle the shock of a calving event and the subsequent readvance of the terminus. We find that after 40 years, the terminus has readvanced beyond its original position at some points and only half-way at others.

Figure 8. Calving terminus locations for Larsen C Ice Shelf prognostic simulation. The contours shown are at the start of the run, immediately after the simulated calving event, and several decades later when the ice shelf has readvanced closer to its original position.

We used the CG(1)/DG(0) pair for the velocity and stress as in the previous examples. To simulate the evolution of the glacier thickness, we used DG(1) elements together with the upwind numerical flux. We found that using a DG discretization was necessary to get a reasonable-looking thickness. When using continuous elements for the thickness, we found that the thickness field would develop spurious oscillations generated at the calving terminus. This finding is to be expected because continuous elements usually fare poorly at advecting sharp features like an advancing ice cliff.

Kangerlussuaq Glacier

Our solver for the dual form was able to simulate the advance and retreat and of the terminus of a real grounded glacier. We ran several instances of the experiment outlined above with different values of the maximum melt rate m ₀. In general, the total volume of ice in the simulated domain oscillates from summer lows to winter highs over a wide range of m ₀ values. With too low or too high a maximum melt rate, there is an additional secular trend in the volume time series as the glacier advances or retreats down the fjord. We found that taking m ₀ on the order of 30 km a⁻¹ makes the yearly-averaged volume roughly constant; see Figure 9. Spread over an inland distance of ~1 km in a 5 km-wide fjord for only the summer season, this gives a total discharge roughly of the same order as the observed value of 24 km³ a⁻¹ (King and others, 2018).

Figure 9. Total volume in km³ of ice in the computational domain, exhibiting summer troughs and winter peaks. The summer maximum melt rate m ₀ is tuned to give a roughly constant yearly average volume, although this simulation shows a small secular trend.

Figure 10 shows the evolution of the calving terminus from a minimum to the following maximum extent. The simulated terminus position oscillates by ~2.5–4 km seasonally, which is close to the observed variation (Schild and Hamilton, 2013). The true calving terminus of the glacier is upstream of the simulated calving terminus, which is likely a consequence of our initialization or other under-parameterized quantities. Additionally, the center line of the true calving terminus is slightly more retreated than the margins. The center line of the simulated terminus, on the other hand, is more advanced than the margins. This discrepancy shows that the ad hoc rule we used to remove ice mass near the terminus is imperfect. Several processes govern the terminus dynamics of Greenland outlet glaciers, including frontal ablation from ocean melt, stress-induced crevassing and calving and back-pressure from sea ice or ice melange in the fjord. We did not attempt to include a real calving law in this exercise. We are nonetheless able to simulate ice-free areas and the advance and retreat of the glacier terminus using the dual form of the momentum balance equations. Closing the gap between the simple demonstrative parameterizations used here and reality is the subject of future work. For example, one could add calving by setting the ice thickness to zero in areas near the glacier terminus where surface crevasses would penetrate to the water line according to the Nye criterion.

Figure 10. Simulated terminus position of Kangerlussuaq Glacier over one half-period, from approximately August at its most retreated to April at its most advanced. The colors of the contours show the time.

Discussion

The momentum balance equation for glacier flow has an alternative, dual expression of the same underlying physics but with different properties and several advantages. The most significant advantage is that the dual form remains solvable in the limit of zero ice thickness. Existing strategies for handling ice-free areas include alteration of the equations or solvers, level-set methods and re-meshing. The dual form accomplishes the same goal and we claim that the challenges of implementing it, while not trivial, are favorable compared to other strategies.

Our comparison of the primal and dual forms shows that using the dual form is more computationally expensive than the primal form because of the greater number of unknowns. On an experiment including calving, the dual form required only slightly more Newton iterations than the primal form with the thickness clamped from below. Whether using the dual form is preferable in general depends on what the simulation aims to achieve. If speed is the main concern then the primal form with clamping is faster. But it introduces a mass-balance error, even more so if the velocity computed in the fictitious ice layer develops a non-zero divergence. If this mass-balance error is not acceptable then the additional cost of using the dual form may be worth it.

Another key feature is that the dual form reverses the behavior of all the non-linearities around the zero-disturbance state. The primal formulation of the problem has a singularity (i.e. terms in the momentum balance equation go to ∞) in the limit as the strain rate goes to zero. Infinite singularities can only be dealt with by fudging the problem itself. In the dual form, however, this singularity becomes instead a degeneracy (terms that go to zero where the usual theory requires positivity). These degeneracies are still a challenge. But the problem, with no modifications, is amenable to solution by approximate Newton methods, as described in the Appendix. Trust region methods (Nocedal and Wright, 2006) might work as well and this remains to be explored.

The story becomes more complicated when we consider the interaction between the dependence on thickness and membrane stress or strain rate. The dual form possesses only degeneracies in the limit as the thickness or membrane stress go to zero. The primal form has a singularity when the strain rate goes to zero, but a degeneracy in the limit as the thickness goes to zero. Moreover, when the thickness goes to zero, the strain rate tends to also go to zero. We hypothesize that this mixture of singularity and degeneracy makes the primal form of the problem impossible to solve in the limit as the thickness goes to zero. We additionally hypothesize that the dual form remains solvable as the thickness goes to zero because it contains only degeneracies. But we have no proof either way and at this stage these hypotheses are at best educated guesses.

The dual formulation does come with several disadvantages. The number of unknowns in the dual formulation is greater than in the primal form, thus putting more pressure on computer memory. The resulting linear systems are indefinite rather than positive-definite. Finally, since the dual form is a mixed problem, it is possible to make bad choices of finite element basis, whereas almost any basis will work for the primal form. We did find, however, that the increased cost of solving the dual problem was not as high as one might expect just based on counting the number of degrees of freedom. We used a fairly naive solution approach (direct factorization) for the linear system in each step of Newton's method in the benchmark for both the primal and dual forms. There may be significant room for improvement on these benchmarks through the use of more sophisticated techniques such as Schur complement preconditioners that use static condensation of the stress degrees of freedom (Boffi and others, 2013).

There are several promising avenues of future work on this problem. Including the stress tensor as an unknown and the constitutive relation as an equation to be solved opens up several possibilities for modifying the physics. Since we do not need to explicitly solve for the stress tensor in terms of the strain rate tensor, we can easily implement composite flow laws like the Goldsby–Kohlstedt law (Goldsby and Kohlstedt, 2001). We could also add a term containing the time derivative of the stress tensor to the constitutive relation to implement Maxwell viscoelasticity. Both of these extensions have historically been difficult to achieve with conventional approaches to glacier flow modeling. Second, the solvability of the dual problem in the limit of zero ice thickness can expand the scope of glaciological data assimilation. For example, it may become possible to assimilate the entire time series of altimetry measurements from ICESat-2 into flow models in a way that constraints not just the elevation of grounded ice, but what areas are free of ice. Finally, more work remains to be done from the applied math side on optimal solution algorithms for these types of problems.

Conclusion

In this paper, we derived the dual form of the glacier momentum balance equation, implemented a numerical solver for it, and demonstrated its use on synthetic and real problems. The key advantage of the dual form is that the problem does not need to be regularized when the strain rate or thickness are equal to 0. The disadvantages are that (1) the dual form has more unknowns and (2) solvers for the resulting non-linear optimization problem require special tuning. Despite these additional costs, we argue that the dual form is worth considering as an alternative to the conventional primal form because of how easy it is to simulate terminus advance and retreat. We did not aim to study directly the holy grail problem of calving laws here. But making it easier to simulate terminus evolution is a virtual requirement for testing these calving laws with computer models.