dual-problems.tex

\NeedsTeXFormat{LaTeX2e}

\documentclass[review,oneside]{igs}
%\documentclass[twocolumn,letterpaper]{igs}

\usepackage{igsnatbib}

% check if we are compiling under latex or pdflatex
\ifx\pdftexversion\undefined
    \usepackage[dvips]{graphicx}
\else
    \usepackage[pdftex]{graphicx}
    \usepackage{epstopdf}
    \epstopdfsetup{suffix=}
\fi


\usepackage{amsmath}
%\usepackage{amsfonts}
%\usepackage{amsthm}
%\usepackage{amssymb}
\usepackage{mathrsfs}
%\usepackage{fullpage}
%\usepackage{mathptmx}
%\usepackage[varg]{txfonts}
\usepackage{color}
%\usepackage[charter]{mathdesign}
%\usepackage{float}
%\usepackage{hyperref}
%\usepackage[modulo, displaymath, mathlines]{lineno}
%\usepackage{setspace}
%\usepackage[titletoc,toc,title]{appendix}
%\usepackage{natbib}
\usepackage{booktabs}
\usepackage{multirow}

%\linenumbers
%\doublespacing

\begin{document}

\title[Terminus evolution via duality]{Numerical simulation of glacier terminus evolution using the dual action principle for momentum balance}

\author[Shapero and de Diego]{Daniel Shapero$^1$, Gonzalo Gonzalez de Diego$^2$}

\affiliation{
    $^1$Polar Science Center, Applied Physics Laboratory, University of Washington, Seattle, WA, USA \\
    $^2$Courant Institute of Mathematical Sciences, New York University, New York, NY, USA \\
    Correspondence: Daniel Shapero $<$shapero@uw.edu$>$
}

\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 \emph{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 nonlinearities.
    This altered character makes it possible to implement numerical solvers for the dual form that work \emph{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.
}

\maketitle

% --------------------
\section{Introduction}

On space and time scales greater than 100m and 1 day, glaciers flow like a viscous, incompressible fluid with a power-law rheology \citep{greve2009dynamics}.
Ice flow is slow enough that the fluid inertia is negligible compared to viscous and gravitational forces, i.e. the flow occurs at very low Reynolds and Froude number.
There are multiple equivalent ways of expressing the momentum balance equations: a conservation law, a variational form, a partial differential equation.
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 \emph{action functional} \citep{dukowicz2010consistent}.
The action functional has units of energy per unit time and can be interpreted as the rate of dissipation of thermodynamic free energy \citep{edelen1972nonlinear}.
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 \citep{nocedal2006numerical}.
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 \emph{Newton decrement}.
In \citet{shapero2021icepack}, we showed how to use this theory to design physics-based convergence criteria.

This work follows in the footsteps of \citet{dukowicz2010consistent} in studying action principles for glacier flow.
Our main contribution is the derivation of an alternative \emph{dual} action principle, distinct from that presented in \citet{dukowicz2010consistent}, from which the momentum conservation equations can be derived.
The most important feature is that \textbf{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.
\textbf{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, \textbf{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 \citep{cornford2013adaptive}.
The Ice Sheet and Sea-Level System Model (ISSM) uses the level set method \citep{bondzio2016modelling, osher1988fronts}.
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 3D geometry \citep{todd2018full}, so ice-free areas are not included in the computational domain at all, and (2) coupling to a discrete element model \citep{benn2017melt}.
These approaches are effective but come with their own drawbacks and implementation challenges.
For example, using the level-set method requires solving a challenging, nonlinear 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 partial differential equations 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 \citet{greve2009dynamics} for glacier dynamics, \citet{weinstock1974calculus} for variational calculus, and \citet{boyd2004convex} for convex optimization.



% --------------
\section{Theory}

\subsection{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 \emph{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 \citep{greve2009dynamics}.

\begin{table}[h]
    \begin{tabular}{cccc}
        \cline{1-4}
        name & symbol & units & rank \\
        \cline{1-4}
        velocity & $u$ & meters years${}^{-1}$ & 1 \\
        strain rate & $\dot\varepsilon$ & years${}^{-1}$ & 2 \\
        viscosity tensor & $\mathscr{C}$ & & 4 \\
        compliance tensor & $\mathscr{A}$ & & 4 \\
        membrane stress & $M$ & megapascals & 2 \\
        basal stress & $\tau$ & megapascals & 1 \\
        ice thickness & $h$ & meters &  \\
        surface elevation & $s$ & meters &  \\
        flow law exponent & $n$ & &  \\
        sliding law exponent & $m$ & &  \\
        fluidity coefficient & $A$ & megapascals${}^{-n}$ years${}^{-1}$ &  \\
        slipperiness coefficient & $K$ & megapascals${}^{-m}$ meters years${}^{-1}$ &
    \end{tabular}
    \caption{Variable, symbol, physical units, and tensor rank -- 1 for vectors, 2 for matrices, etc.}
    \label{tab:symbols}
\end{table}

The equations of motion are solved in a two-dimensional domain $\Omega$.
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 $\tau$ and the \emph{membrane stress} tensor $M$, a rank-2 tensor with units of stress that results from applying the low-aspect ratio assumption to the full 3D 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
\begin{equation}
    \nabla\cdot hM + \tau - \rho_I gh\nabla s = 0,
    \label{eq:membrane-stress-conservation}
\end{equation}
where $\rho_I$ is the ice density, and $g$ the gravitational acceleration.
In addition to \eqref{eq:membrane-stress-conservation}, we need to know a constitutive relation between the membrane stress tensor and the depth-averaged strain rate tensor
\begin{equation}
    \dot\varepsilon = \frac{1}{2}\left(\nabla u + \nabla u^\top\right).
    \label{eq:strain-rate}
\end{equation}
In order to simplify the notation later, we introduce the dimensionless rank-4 tensor $\mathscr{C}$ defined by
\begin{equation}
    \mathscr{C}\dot\varepsilon = \frac{\dot\varepsilon + \text{tr}(\dot\varepsilon)I}{2}.
    \label{eq:elasticity-tensor}
\end{equation}
The tensor $\mathscr{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 $\mathscr{C}$ as
\begin{equation}
    |\dot\varepsilon|_{\mathscr{C}}^2 = \dot\varepsilon : \mathscr{C}\dot\varepsilon.
\end{equation}
Alternatively, in index notation, the square norm is $|\dot\varepsilon|_{\mathscr C}^2 = \mathscr{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:
\begin{equation}
    M = 2A^{-\frac{1}{n}}|\dot\varepsilon|_{\mathscr C}^{\frac{1}{n} - 1}\mathscr{C}\dot\varepsilon
    \label{eq:constitutive-relation}
\end{equation}
where $A$ is the depth-averaged fluidity coefficient and $n \approx 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$, i.e.
\begin{equation}
    \tau = -C|u|^{\frac{1}{m} - 1}u.
    \label{eq:sliding-law}
\end{equation}
Weertman sliding has $m = n$, while perfectly plastic sliding has $m = \infty$.
Recent research suggests alternative forms that transition between Weertman-type sliding at low speeds and perfectly plastic sliding at higher speeds \citep{minchew2020toward}.
For illustrative purposes equation \eqref{eq:sliding-law} 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 $\Gamma$, the membrane stresses at the cliff face are balanced by pressures from any proglacial water body.
This is a Neumann boundary condition:
\begin{equation}
    -hM\cdot\nu = \frac{1}{2}\left(\rho_Igh^2 - \rho_Wgh_W^2\right)\nu,
    \label{eq:terminus-bc}
\end{equation}
where $h_W$ is the water depth and $\nu$ is the unit outward-pointing normal vector to the terminus.

We can then combine equations \eqref{eq:membrane-stress-conservation}-\eqref{eq:sliding-law} and the boundary condition \eqref{eq:terminus-bc} to arrive at a second-order, nonlinear elliptic system of differential equations for $u$.

For modelling a floating ice shelf, the friction term in \eqref{eq:membrane-stress-conservation} is zero.
Moreover, assuming the ice is in hydrostatic equilibrium allows us to write the surface elevation in terms of the thickness:
\begin{equation}
    s = (1 - \rho_I / \rho_W) h.
    \label{eq:flotation-condition}
\end{equation}
As a result, the momentum conservation equation reduces to
\begin{equation}
    \nabla\cdot hM - \frac{1}{2}\varrho g\nabla h^2 = 0.
\end{equation}
where $\varrho = (1 - \rho_I / \rho_W)\rho_I$ is the reduced density of ice over seawater.

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 the condition \eqref{eq:flotation-condition} 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 \citep{schoof2007mis, durand2009mis, favier2012mis}.

To complete our description of the dynamics, the thickness evolves in time according to the following depth-averaged mass conservation equation:
\begin{equation}
    \frac{\partial h}{\partial t} + \nabla \cdot hu = \dot a - \dot m,
    \label{eq:conservation-mass}
\end{equation}
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.


% --------------------------------
\subsection{Primal action principles}
\label{sec:primal-action-principles}

The main idea behind \emph{action} or \emph{minimization} principles is that some partial differential equations (PDE) really express the fact that their solutions are extrema of a given \emph{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 \citep{bassis2010hamilton, dukowicz2010consistent, brinkerhoff2013data, shapero2021icepack}.
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 \emph{Euler-Lagrange} equation for the functional.
We will not repeat it here but see \citet{weinstock1974calculus}.
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 \eqref{eq:constitutive-relation} for $M$ can be expressed as the derivative of a certain scalar quantity.
In index notation,
\begin{equation}
    M_{ij} = \frac{\partial}{\partial\dot\varepsilon_{ij}}\left(\frac{2n}{n + 1}A^{-\frac{1}{n}}|\dot\varepsilon|_{\mathscr C}^{\frac{1}{n} + 1}\right)
    \label{eq:M-anti-derivative}
\end{equation}
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 \eqref{eq:sliding-law}, we can observe that
\begin{equation}
    \tau_i = -\frac{\partial}{\partial u_i}\left(\frac{m}{m + 1}C|u|^{\frac{1}{m} + 1}\right).
    \label{eq:tau-anti-derivative}
\end{equation}
Using equations \eqref{eq:M-anti-derivative} and \eqref{eq:tau-anti-derivative}, one can show that the SSA momentum balance are the Euler-Lagrange equations for the following action functional \citep{dukowicz2010consistent}:
\begin{align}
    & J(u) = \nonumber\\
    & \qquad\int_\Omega\Bigg\{\frac{2n}{n + 1}hA^{-\frac{1}{n}}|\dot\varepsilon|_{\mathscr{C}}^{\frac{1}{n} + 1} + \frac{m}{m + 1}C|u|^{\frac{1}{m} + 1} \nonumber\\
    & \qquad\qquad\qquad\qquad\qquad\qquad\qquad + \rho_I gh\nabla s\cdot u\Bigg\}dx  \nonumber \\
    & \qquad\qquad + \frac{1}{2}\int_\Gamma\left(\rho_Igh^2 - \rho_Wgh_W^2\right)u\cdot\nu\; d\gamma
    \label{eq:ssa-primal-action}
\end{align}
Note that the action has units of energy per unit time, or power.
The final summand of equation \eqref{eq:ssa-primal-action} 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 nonlinear 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 \citep{edelen1972nonlinear}.
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.


% ------------------------------
\subsection{Dual action principles}
\label{subsec: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 \citet{attouch2014variational}.

Equations \eqref{eq:membrane-stress-conservation}-\eqref{eq:terminus-bc} 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 equation \eqref{eq:ssa-primal-action}.
In deriving the differential equation, we used the constitutive relation and the sliding law (equations \eqref{eq:constitutive-relation} and \eqref{eq:sliding-law}) to eliminate the membrane stress $M$ and basal stress $\tau$.

We could instead \textbf{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 $\tau$.
One way to motivate the dual problem is to ask: \textbf{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 equation \eqref{eq:elasticity-tensor}, we defined the rank-4 tensor $\mathscr{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 $\mathscr{A}$ of this tensor $\mathscr{C}$ is:
\begin{equation}
    \mathscr{A}M = \frac{M - \frac{1}{d + 1}\text{tr}(M)I}{2}
\end{equation}
where $d = 2$ is the space dimension.
The tensor $\mathscr{A}$ plays an analogous role to the compliance tensor in linear elasticity.
With this definition in hand, we can invert equation \eqref{eq:constitutive-relation} to get $\dot\varepsilon = 2A|M|_{\mathscr{A}}^{n - 1}\mathscr{A}M$.
The sliding law is much easier to invert: $u = -C^{-m}|\tau|^{m - 1}\tau$.
We can then define a slipperiness coefficient $K = C^{-m}$.

With these preparations in place, the dual form is:
\begin{align}
    & L(u, M, \tau) = \nonumber\\
    & \quad \int_\Omega\Bigg\{\frac{2}{n + 1}hA|M|_{\mathscr A}^{n + 1} + \frac{1}{m + 1}K|\tau|^{m + 1} \nonumber\\
    & \qquad\qquad - hM:\dot\varepsilon(u) + \tau\cdot u - \rho_Igh\nabla s\cdot u\Bigg\}dx \nonumber \\
    & \qquad - \frac{1}{2}\int_\Gamma\left(\rho_Igh^2 - \rho_Wgh_W^2\right)u\cdot\nu\;d\gamma
    \label{eq:ssa-dual-form}
\end{align}
We can then evaluate the variational derivatives of $L$ with respect to $u$, $M$, and $\tau$, 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
\begin{equation}
    \left\langle\frac{\partial L}{\partial M}, N\right\rangle = \int_\Omega\left\{2hA|M|_{\mathscr{A}}^{n - 1}\mathscr{A}M - h\dot\varepsilon(u)\right\} : N\; dx = 0.
\end{equation}
This equation is the variational form of the inverse of the constitutive relation (equation \eqref{eq:constitutive-relation}).
Next, taking the variational derivative of $L$ with respect to $\tau$ along some perturbation $\sigma$, we get
\begin{equation}
    \left\langle\frac{\partial L}{\partial\tau}, \sigma\right\rangle = \int_\Omega\left\{K|\tau|^{m - 1}\tau + u\right\}\cdot\sigma\;dx = 0
\end{equation}
which is the inverse of the sliding law of equation \eqref{eq:sliding-law}.
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
\begin{align}
    \left\langle\frac{\partial L}{\partial u}, v\right\rangle & = \int_\Omega\left\{-hM : \dot\varepsilon(v) + \tau\cdot v - \rho_Igh\nabla s\cdot v\right\}dx \nonumber\\
    & \qquad -\frac{1}{2}\int_{\Gamma}\left(\rho_Igh^2 - \rho_Wgh_W^2\right)v\cdot\nu\;d\gamma
\end{align}
which is not exactly what we want.
We can use integration by parts to push the symmetric gradient of $v$ over onto $h\,M$ in the first term however:
\begin{align}
    \ldots & = \int_\Omega\left\{\nabla\cdot hM + \tau - \rho_Igh\nabla s\right\}\cdot v\;dx \nonumber\\
    & \quad + \int_\Gamma\left(hM\cdot\nu - \frac{1}{2}\left(\rho_Igh^2 - \rho_Wgh_W^2\right)\nu\right)\cdot v\;dx.
\end{align}
If we require that this is equal to 0 for all perturbation fields $v$, we recover both conservation of membrane stress (equation \eqref{eq:membrane-stress-conservation}) and the boundary condition (equation \eqref{eq:terminus-bc}).
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.

\begin{figure}[h]
    \includegraphics[width=0.48\linewidth]{demos/singularity/primal.pdf}\includegraphics[width=0.48\linewidth]{demos/singularity/dual.pdf}
    \caption{The viscous part of the action for the primal problem as a function of the strain rate (left) and for the dual problem as a function of the stress (right).
    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.}
    \label{fig:primal-vs-dual}
\end{figure}

The important feature of the dual formulation is that \textbf{the nature of the nonlinearity has changed}.
In the primal action shown in equation \eqref{eq:ssa-primal-action}, the nonlinearity consists of the strain rate raised to the power $\frac{1}{n} + 1$.
Since $n > 1$, the nonlinearity 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 nonlinearity 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 \ref{fig:primal-vs-dual}.


% ----------------------
\section{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.


\subsection{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\citep{shapero2021icepack}.
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^2$-norm difference between the exact solution and the solutions obtained using $CG(k)$ finite elements is $\mathscr{O}(\delta x^{k + 1})$ where $\delta 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
\begin{equation}
    h = h_0 - \delta h \cdot x / L_x
\end{equation}
in a domain of length $L_x$ = 20km.
With a constant value of the fluidity coefficient $A$, the velocity in the $x$ direction is
\begin{equation}
    u_x = u_0 + L_x A \left(\frac{\varrho g h_0}{4}\right)^n\left(1 - \left(1 - \frac{\delta h \cdot x}{h_0\cdot L_x}\right)^{n + 1}\right)
\end{equation}
where $\varrho = \rho_I(1 - \rho_I/\rho_W)$.
We use a 2D 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 \citep{sympy}.


\subsection{Comparison with primal form on slab glacier}

\begin{figure}[h]
\centering
    \includegraphics[width = 0.8\linewidth]{demos/slab/figures/marine_ice_sheet.pdf}
    \caption{Setup for modelling a slab of ice on an inclined bed flowing into the ocean.
    At $x = 0$ we enforce a thickness $h = 500\,\mathrm{m}$ in order to approach a parallel slab of ice far upstream of the grounding line. The dotted line is sea level.}
    \label{fig:domain_parallel_slab}
\end{figure}

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 (figure \ref{fig:domain_parallel_slab}).
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\, \mathrm{m}$ and the bedrock angle to $\alpha = 1^\circ$. The bed is given by the expression
\begin{align}
    b(x) = 1500 - x\tan{\alpha}\,\mathrm{m}.
\end{align}
%
For the material parameters, we set $n = 3$, $A = 10^{-24}\,\mathrm{Pa}^{-3}\,\mathrm{s}^{-1}$, $C = 10^{-6}\,\mathrm{Pa}\,\mathrm{m}^{-1/3}\mathrm{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 far 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
\begin{align}\label{eq:u_slab}
    u = \left(\frac{\rho_I g h \tan(\alpha)}{C}\right)^n.
\end{align}
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 regularise the primal formulation of the momentum balance equation.

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 equation \eqref{eq:conservation-mass} and the flotation condition \eqref{eq:flotation-condition}, effectively yielding a free boundary problem.


\subsection{Comparison against primal form on gibbous ice shelf}

The next comparison exercise uses the synthetic ``gibbous'' ice shelf test case from \S5.3 of \citet{shapero2021icepack}
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 form.
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) \times DG(0)$ elements for the velocity and stress has an additional three unknowns per triangle.
The Euler formula (\#vertices - \#edges + \#triangles $\approx$ 2) implies that there are approximately twice as many triangles as there are vertices.
Consequently there are about 4$\times$ 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$\times$ 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.


\subsection{Demonstration on calving of 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 \citep{larour2021physical}.
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:
\begin{enumerate}
    \item 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.
    \item Extrapolate the ice thickness and velocity onto a larger spatial domain, making the ice thickness 0 in ice-free areas.
    \item Run the simulation using the mass and dual momentum balance from the start date of 2015 until the calving event in 2017.
    \item Digitize the terminus position immediately after the calving event by hand and use the digitized terminus position to define an ice mask.
    \item Using this mask, set the ice thickness to zero over the spatial extent of the calved area.
    \item Run the simulation for 40 years after the calving event to see how the terminus advances again.
\end{enumerate}


\subsection{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 \citep{enderlin2014improved, mouginot2019forty}.
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:
\begin{enumerate}
    \item Estimate the slipperiness (the coefficient $K$ in the sliding law $u|_{z = b} = -K|\tau_b|^{n - 1}\tau$) 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.
    \item Extrapolate the thickness, surface elevation, velocity, and friction coefficient onto a large spatial domain that extends further down Kangerlussuaq Fjord.
    \item Run the simulation using the mass and dual momentum balance equations for one year in order to propagate out any initial transients.
        This stage uses only surface mass balance and thus permits the glacier to advance down the fjord.
    \item Turn on a time-periodic ablation field near the terminus in order to represent the effects of summer melt and calving and ran the simulation for a further four years.
        This ablation field forces the terminus to advance and retreat.
\end{enumerate}

To initialize the simulation, we use version 3 of the BedMachine Greenland data set for ice thickness and surface elevation \citep{morlighem2017bedmachine} and the MEaSUREs annual velocity mosaic from 2015-2016 \citep{joughin2010greenland} to infer the basal friction.
To force the mass conservation equation \eqref{eq:conservation-mass}, we need to provide a surface mass balance (SMB) field $\dot a$ and a melt rate $\dot m$.

We use a surface mass balance field that varies linearly with elevation:
\begin{equation}
    \dot a \approx a_0 + \frac{\delta a}{\delta s}\cdot s
\end{equation}
where $a_0$ is the SMB at sea level and $\delta a/\delta s$ is the SMB lapse rate.
To fit the parameters $a_0$ and $\delta a/\delta s$, we used output from 2006-2021 of version 3.12 of the Mod\`ele Atmosph\'erique R\'egional (MAR) \citep{fettweis2020grsmbmip}.
This regional climate model has been tested extensively for the polar regions and for Greenland.
The fit had $r^2 = 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 $\mu$.
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 $\mu$ as the minimizer of
\begin{equation}
    J(\mu) = \frac{1}{2}\int_\Omega\left((\mu - \mathbf{1}_{\{h > 0\}})^2 + \alpha^2|\nabla\mu|^2\right)dx
\end{equation}
where $\alpha$ is some smoothing length, and $\mathbf{1}_{\{h > 0\}}(x)$ is equal to 1 if $h(x) > 0$ and 0 if $h(x) = 0$.
Here we choose $\alpha$ to be 1km, so the mask field rapidly approaches 1 within roughly one ice thickness of the terminus.
The mask field $\mu$ is recalculated in every timestep.
Finally, we set the melt rate at time $t$ as
\begin{equation}
    \dot m = m_0(1 - \mu)\min\{0, \cos(2\pi t)\}
\end{equation}
where $m_0$ 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.


% ---------------
\section{Results}

We implemented a solver for the dual form of the SSA using the Firedrake package \citep{FiredrakeUserManual}.
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.


\subsection{Verification on solvable test cases} \label{sec:linear-ice-shelf}

\begin{figure}[h]
    \begin{center}
        \includegraphics[width=0.75\linewidth]{demos/convergence-tests/results.pdf}
    \end{center}
    \caption{Relative $L^2$-norm errors for approximate solutions to the analytical ice shelf (left) and ice stream (right) 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.}
    \label{fig:linear-glacier-convergence-rate}
\end{figure}

We tested meshes with between 16 and 256 cells to a side and we used both $CG(1) \times DG(0)$ and $CG(2) \times 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^2$ norm have the expected asymptotic convergence rates of $\mathscr{O}(\delta x^2)$ for linear velocity elements and $\mathscr{O}(\delta x^3)$ for quadratic in both the ice shelf and ice stream test cases; see figure \ref{fig:linear-glacier-convergence-rate}.

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 \ref{fig:linear-glacier-convergence-rate} 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.


\subsection{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)\times CG(2)$, and with a dual method that computes the velocity, membrane stress and thickness in the space $CG(2)\times DG(1)\times 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 1D form of the equation, so the relevant term in the variational form of the momentum balance equation is
\begin{equation}
    \langle F(u), v\rangle = \int_\Omega \left\{2hA^{-1/n}|\partial_xu^2 + \epsilon^2|^{(n - 2)/2}\partial_xu\cdot\partial_xv + \ldots\right\}dx
\end{equation}
We consider a sequence of regularization parameters $\epsilon$ between $1\, \mathrm{yr}^{-1}$ and $10^{-12}\,\mathrm{yr}^{-1}$.
The results for the grounding line position are displayed in Table \ref{tab:slab}.
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\,\mathrm{m}$, $u$ is equal to \eqref{eq:u_slab}, and $M = 0$.
The initial guess for the grounding line position is set to the point where the flotation condition \eqref{eq:flotation-condition} holds for the constant thickness slab.
We plot the values of the relative Newton residual in figure \ref{fig:newton-its}.
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^{-14}\,\mathrm{yr}^{-1}$, the relative Newton residual no longer reaches the minimum tolerance of $10^{-8}$ that we set for this problem.

Figure \ref{fig:newton-its} 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.

%\renewcommand{\arraystretch}{1.25}
\begin{table}[h]
\centering
\caption{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.}
\label{tab:slab}
\begin{tabular}{ccccc}
\toprule
solver & $\epsilon\,\mathrm{yr}^{-1}$ & $x_g\,\mathrm{km}$ & $h(x_g)\,\mathrm{m}$ & iterations \\
\midrule
\input{demos/slab/tables/table_alpha1.00.txt}
\end{tabular}
\end{table}

\begin{figure}[h]
    \centering
    \includegraphics[width=\linewidth]{demos/slab/figures/newton_its_alpha1.00.pdf}
    \caption{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.}
    \label{fig:newton-its}
\end{figure}


\subsection{Gibbous ice shelf} \label{sec:gibbous-ice-shelf}

\begin{figure}[h]
    \begin{center}
        \includegraphics[width=0.85\linewidth]{demos/gibbous-ice-shelf/gibbous.pdf}
    \end{center}
    \caption{The 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.}
    \label{fig:gibbous}
\end{figure}

\begin{figure}[h]
    \begin{center}
        \includegraphics[width=0.85\linewidth]{demos/gibbous-ice-shelf/volumes.pdf}
    \end{center}
    \caption{Total volume of ice in the shelf over time.
    The different spin-up and experimental phases are labelled.
    Note how the finer spin-up equilibrates to a smaller ice volume than the coarser spin-up.}
    \label{fig:gibbous-calving-volumes}
\end{figure}

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

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$\times$ and 2.7$\times$ 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$\times$ 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 \ref{fig:gibbous-calving-volumes} 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.
Figure \ref{fig:gibbous}(d) shows the thickness of the ice shelf immediately after the calving event (e) and (f) show the magnitude of the change in speed and stress respectively.
In particular, the stress field shows a discontinuity at the new calving front as expected.

\begin{figure}[h]
    \begin{center}
        \includegraphics[width=0.75\linewidth]{demos/gibbous-ice-shelf/counts.pdf}
    \end{center}
    \caption{The 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.}
    \label{fig:gibbous-residuals}
\end{figure}

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 1mm.
Figure \ref{fig:gibbous-residuals} 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.


\subsection{Larsen C Ice Shelf}

\begin{figure}[h]
    \begin{center}
        \includegraphics[width=0.75\linewidth]{demos/larsen/contours.pdf}
    \end{center}
    \caption{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.}
    \label{fig:larsen-terminus-position}
\end{figure}

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 \ref{fig:larsen-terminus-position}.
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.

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.


\subsection{Kangerlussuaq Glacier}

\begin{figure}[h]
    \begin{center}
        \includegraphics[width=0.75\linewidth]{demos/kangerd/volumes.pdf}
    \end{center}
    \caption{Total volume in km${}^3$ of ice in the computational domain, exhibiting summer troughs and winter peaks.
    The summer maximum melt rate $m_0$ is tuned to give a roughly constant yearly average volume, although this simulation shows a small secular trend.}
    \label{fig:kangerd-volumes}
\end{figure}

\begin{figure}[h]
    \begin{center}
        \includegraphics[width=0.75\linewidth]{demos/kangerd/contours.pdf}
    \end{center}
    \caption{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.}
    \label{fig:kangerd-contours}
\end{figure}

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_0$.
In general, the total volume of ice in the simulated domain oscillates from summer lows to winter highs over a wide range of $m_0$ 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_0$ on the order of 30 km/yr makes the yearly-averaged volume roughly constant; see figure \ref{fig:kangerd-volumes}.
Spread over an inland distance of roughly 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${}^3$/yr \citep{king2018seasonal}.

Figure \ref{fig:kangerd-contours} shows the evolution of the calving terminus from a minimum to the following maximum extent.
The simulated terminus position oscillates by roughly 2.5-4km seasonally, which is close to the observed variation \citep{schild2013seasonal}.
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 centerline of the true calving terminus is slightly more retreated than the margins.
The centerline 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.



% ------------------
\section{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 implementating 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 moreso 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 nonlinearities around the zero-disturbance state.
The primal formulation of the problem has a singularity (i.e. terms in the momentum balance equation go to $\infty$) 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 \citep{nocedal2006numerical} 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 possess 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 \citep{boffi2013mixed}.

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 \citep{goldsby2001superplastic}.
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.


% ------------------
\section{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 nonlinear 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.



% ----------------------------------
\section{Code and data availability}

The complete source code used for the simulations described in this paper is available at: https://github.com/icepack/dual-problems.git.



% ------------------------
\section{Acknowledgements}

We thank Robert Baraldi for many helpful discussions.
DRS was supported by the US National Science Foundation (grant OAC-1835321) and National Aeronautics and Space Administration (grants 80NSSC20K0954 and 80NSSC21K1003).
GGdD was supported by the University of Oxford Mathematical Institute Graduate Scholarship.

% ---------------------
\bibliographystyle{igs}
\bibliography{dual-problems.bib}

\pagebreak



% -------------------------
\appendix\section{Appendix}

We have largely focused on the dual form of the SSA momentum balance as an alternative to the primal form with certain favorable numerical properties, such as solvability at zero thickness.
We do not have conclusive answers about the best way to discretize and solve the dual form of SSA.
In the following, we detail some of the techniques that we used.
A further publication will explore these issues in greater detail.

\subsection{Discretization by finite elements}
\label{app:discretization}

Roughly any conforming finite element basis is stable for the primal form of symmetric, positive-definite elliptic equations, such as the diffusion and elasticity equations, as long as the mesh is regular.
The most common choice is to use piecewise-continuous polynomials of a given degree $k$ on triangles, or the tensor product of polynomials on quads.
We will refer to this basis as $CG(k)$.
While dual formulations have many advantages, the main challenge to overcome is that most choices of basis are unstable -- the resulting linear systems are either singular or their inverses have unbounded norm in the limit as the mesh is refined.
For example, using $CG(k)$ elements for the temperature and the product $CG(k)^2$ for the flux is an unstable discretization of the dual form of the diffusion equation.
Making matters even harder, the SSA and other problems for a pair of vector and tensor fields have an additional invariant to enforce -- the symmetry of the stress tensor -- which can be difficult to achieve in practice.

The question of how to choose basis functions that give a stable discretization of dual problems is the subject of \emph{mixed} finite element methods.
This subject is covered in great detail in \citet{boffi2013mixed}.
There is, however, a wide chasm between the motivation for using dual formulations in most of the finite element literature and our reasons for applying them to glacier momentum balance.
The big motivating problem for dual formulations in the finite element literature is linear elasticity.
In that setting, the goal is to compute the stress tensor with high accuracy in order to make sure that it does not exceed some failure threshold for the material.
Using the dual form of the elasticity equations offers the promise of approximating the stress tensor with a higher order of accuracy than the primal form.
Finding stable finite element bases for the dual form of the elasticity equations is a holy grail problem because of its potential impact on engineering practice.

It might seem at first blush as if the heavy focus on finding stable discretizations of the dual form of the elasticity equations is beneficial to us because the SSA is formally similar to 2D elasticity, even though these equations have different provenance.
Our purpose for using the dual form, however, is not to obtain a more accurate resolution of the membrane stress tensor -- we are only interested in the dual form because of how it changes the character of the nonlinearities in the SSA.
With this goal in mind, there are several choices that we make differently from how they are done in the finite element literature.
These are of a technical nature and not of special interest to most glaciologists, but we include them here for the sake of completeness.
A typical dual formulation of elasticity would assume that:
\begin{enumerate}
    \item the displacements live in the function space $L^2(\Omega, \mathbb{R}^d)$, i.e. the space of square-integrable vector fields, and
    \item the stresses live in the space $H^{\text{div}}(\Omega, \mathbb{R}_{\text{sym}}^{d \times d})$ of square-integrable symmetric tensor fields whose divergences are also square-integrable.
\end{enumerate}
This $L^2 \times H^{\text{div}}$ formulation offers the best possible asymptotic accuracy for the stress tensor.
The dual form of the problem with these assumptions is different from what we wrote down in equation \eqref{eq:ssa-dual-form} -- the gradient of $u$ is instead pushed over as a stress divergence.
Moreover, with the $L^2 \times H^{\text{div}}$ form, Dirichlet boundary conditions become natural and Neumann conditions become essential.
Finding stable bases for the $L^2\times H^{\text{div}}$ form requires very sophisticated finite element bases.
At the simplest end of the spectrum, one can enrich the stress space by cubic bubbles \citep{brezzi1993mixed}.
A host of more complex approaches are possible \citep{arnold1984peers, arnold2002mixed}.

Although it is almost completely unheard of in the literature on mixed finite elements, we make a different but equally valid set of assumptions.
We instead assume that
\begin{enumerate}
    \item the velocities live in the function space $H^1(\Omega, \mathbb{R}^d)$ of vector fields that are square-integrable and have square-integrable derivatives, and
    \item the membrane stress tensor lives in the space $L^2(\Omega, \mathbb{R}^{d \times d}_{\text{sym}})$ of square-integrable symmetric tensor fields.
\end{enumerate}
With this $H^1 \times L^2$ dual form, Dirichlet conditions remain essential and Neumann conditions natural.
Finding a stable finite element basis is much more straightforward for the $H^1\times L^2$ form of the problem.
We use the space $CG(k)^d$ of continuous piecewise-polynomial vector fields for the velocities, and $DG(k)^{d\times d}_{\text{sym}}$ of discontinuous piecewise-polynomial symmetric tensor fields for the membrane stress.



% ---------------------------------------
\subsection{Solution by Newton-type methods}
\label{app:solution}

The finite element method reduces the infinite-dimensional optimization problems that we have described into finite-dimensional ones.
All that remains is to decide how to solve the resulting finite-dimensional optimization problems.

We can approximate a minimizer for the primal form of the action functional using standard Newton line search algorithms \citep{shapero2021icepack}.
But the primal form of the momentum balance equation has singularities in the limit as the strain rate tensor approaches 0.
When we calculate the derivative of the action, these singularities are multiplied by 0 in such a way that they become removable, i.e. they have a finite limit.
In floating-point arithmetic, however, evaluating an expression with a removable singularity does not always produce the right limit.
Moreover, the second derivative of the action does have genuine infinite singularities, and we need to be able to calculate the second derivative or some approximation to it in order to use Newton-type methods.
The usual remedy is to introduce a smoothing factor $\delta$ into the action that rounds off the behavior around $\dot\varepsilon = 0$.
Regularizing the action functional makes the minimization problem solvable but very ill-conditioned.
Additionally, for some simulations the ice thickness can go to zero, which makes the minimization problem difficult or impossible to solve numerically.
The usual remedy for this is to clamp the thickness from below at some fixed value, say 1m or 10m.
Where the ice thickness approaches zero, usually the strain rate does as well.
In these scenarios, we are certain to encounter the worst behavior possible associated with the singularity at zero strain rate.

The dual form, on the other hand, does not have infinite singularities around zero strain rate.
Instead, the action functional has \emph{degeneracies} -- terms that go to zero where, in a nicer problem, they would stay strictly positive.
(See again figure \ref{fig:primal-vs-dual}.)
Degeneracies are not good news either.
In order to use a Newton-type algorithm to find a critical point of the dual action $L$, we compute a search direction by solving the linear system
\begin{equation}
    d^2L\cdot\left(\begin{matrix} v \\ N \\ \sigma\end{matrix}\right) = -dL.
\end{equation}
We know that the second derivative of $L$ has the structure of a saddle-point matrix.
Usually one assumes that certain blocks of this matrix are symmetric and strictly positive-definite in order to guarantee the existence of a solution \citep{boffi2013mixed}.
When the problem is degenerate, we no longer have these guarantees.
We still know that $L$ has a unique saddle point because it is \emph{strictly} convex with respect to $M$ and $\tau$, the problem is that it fails to be \emph{strongly} or \emph{uniformly} convex.
There are workable remedies for this issue that do not degrade the conditioning of the problem to the same extent as regularization does for the primal problem.

Newton's method with line search guarantees second-order convergence for nice problems.
In the event that the second derivative has degeneracies, we can instead try to compute a search direction b solving the perturbed system
\begin{equation}
    \left(d^2L + \lambda\cdot d^2G\right)\left(\begin{matrix} v \\ N \\ \sigma\end{matrix}\right) = -dL
    \label{eq:regularized-newton-step}
\end{equation}
where $G$ is some strongly convex function of $M$ and $\tau$ and $\lambda$ is a small parameter.
For example, one reasonable choice is to take
\begin{equation}
    G = \frac{1}{2}\int_\Omega\left(\max\{h, h_{\text{min}}\}A'|M|_{\mathscr{A}}^2 + K'|\tau|^2\right)dx
\end{equation}
for some constants $A'$, $K'$ having the right units and for some minimum thickness $h_{\text{min}}$ on the order of 1-10m.
The addition of $d^2G$ regularizes the search directions.
It does \emph{not} regularize or perturb what solution we are looking for, only how we look for it.

Regularizing the search directions sacrifices the second-order convergence rate of Newton's method.
It does, however, achieve faster convergence than typical first-order quasi-Newton methods like BFGS \citep{nocedal2006numerical}.


\end{document}