A stabilized finite element method for calculating balance velocities in ice sheets

We present a numerical method for calculating vertically averaged velocity fields using a mass conservation approach, commonly known as balance velocities. This allows for an unstructured grid, is not dependent on a heuristic flow routing algorithm, and is both parallelizable and efficient. We apply the method to calculate depth-averaged velocities of the Greenland Ice Sheet, and find that the method produces grid-independent velocity fields for a sufficient parameterization of horizontal plane stresses on flow directions. We show that balance velocity can be used as the forward model for a constrained optimization problem that can be used to fill gaps and smooth strong gradients in InSAR velocity fields.


Introduction
Balance velocities are useful in evaluating the dynamics of ice sheets, as a means to fill missing velocity data (e.g., Joughin et al., 2010), and as an additional point of comparison for data-derived and modeled velocities (Bamber et al., 2000).Stemming from a statement of mass conservation, balance velocity provides an intuitive means for understanding the distribution of flux within an ice sheet.It has often provided estimates of velocity with better fidelity to data than even advanced ice sheet models, while relying on fewer assumptions.It also gives us the means to assess the distance from equilibrium of an extant ice sheet.
Heretofore, balance velocity has been calculated by applying discrete routing algorithms to spatially distribute flux.These have traditionally been drawn from the hydrological literature (e.g., Tarboton, 1997;Budd and Warner, 1996).To leading order, hydrological routing and glaciological routing are similar; flow directions in both cases are governed by driving stresses, which are determined by surface slope.In overland routing of liquid water, this method is appropriate.However, in glacial ice, the flow direction is also determined by horizontal plane stresses (and, to a lesser extent, vertical resistive stresses), and neglecting these terms yields an overconvergent pattern.This emphasis on local slopes also tends to exacerbate grid dependence, causing the same routing algorithm to produce markedly different velocity fields for different grid resolutions (LeBrocq et al., 2006).Algorithms overcome this by using a spatially averaged slope rather than a purely local slope, with smoothing lengths and the shape of the averaging filter derived heuristically (Testut et al., 2003) or from theoretical results of parameterizing horizontal plane stresses (Kamb and Echelmeyer, 1986).
The aim of this paper is to show how balance velocity can be accomplished by solving a partial differential equation for the conservation of mass using finite elements rather than discrete flow routing algorithms.An unstructured grid also allows for enhanced resolution in regions of special interest, analogous to the mesh refinement used by contemporary next-generation ice sheet models (Larour et al., 2012;Seddik et al., 2012;Brinkerhoff and Johnson, 2013), or for simply scaling grid size by ice thickness.This approach also makes the incorporation of horizontal plane stress gradients straightforward by parameterizing horizontal plane stresses by solving an additional linear system.To these ends, we present the governing equations and the method of their numerical solution with finite elements.We apply this method to the Greenland Ice Sheet and show that this approach yields quality and grid-independent balance velocity fields.
In addition to the novel, but basic, method for computing balance velocities, we also present a method by which balance velocities can be used to fill gaps and smooth spurious gradients in InSAR-derived velocity data (e.g., Joughin et al., 2010).This is often advantageous, since further applications, such as inversion for basal traction or computing local stress balances, depend on having a smooth and complete velocity field.The method relies on minimizing a misfit functional over the velocity field with respect to error bounded thickness, apparent surface mass balance, and flow direction.

Continuum formulation
For an incompressible fluid, conservation of mass is stated as where u is the three-dimensional fluid velocity field, with kinematic boundary conditions on the surface S and bed and respectively.Vertically integrating Eq. ( 1), applying the Leibniz rule, and substitution of Eqs.
(2) and (3) yields a vertically averaged statement for conservation of mass, commonly called the continuity equation with surface mass balance ȧ, basal melt m b , and thickness H . ∇ • is the divergence operator in the two horizontal directions, and u = [u, v] is the vertically averaged horizontal velocity vector.We henceforth drop the parallel bars, and assume that all vectors and operators work on the horizontal plane.This equation is well known to ice sheet modelers as the prognostic equation for evolving the geometry of an ice sheet.In this case, we assume an estimate of ∂ t H , and group it with the other source terms, yielding where 5) is often used to calculate H (Morlighem et al., 2011;Johnson et al., 2012).Here, we assume that H is known, and instead use Eq. ( 5) to calculate u.As stated, the system is underdetermined, with only one equation for both velocity components.For closure, we restate the problem in terms of flow direction N and speed U = u 2 (where • 2 denotes the standard L 2 norm), such that This gives the scalar equation for unknown U Flow direction is specified as the solution to the problems with boundary condition and The solution to Eq. ( 8) is equivalent to the application of a Gaussian average of variable length scale lH to the driving stress τ d of the type suggested by Kamb and Echelmeyer (1986).Theoretical work typically expresses stress coupling length scales in terms of ice thicknesses, hence the notation lH ; l is the number of ice thicknesses over which horizontal plane stress coupling should act.Flow direction N is then proportional to the smoothed driving stress τ s with unit normalization.In the case where the boundary of the computational domain corresponds to the complete boundary of an ice mass (balance velocity for all of Greenland, say), no boundary condition need be specified, as the solution is implicitly defined to be zero at the ice divide due to the problem geometry.When considering a partial domain, a Dirichlet condition must be specified once per flow line.

Discretization and stabilization
Equations ( 5), (8), and (10) are closed, and can be used to calculate balance velocity.We use the finite element method in order to discretize the governing equations.Equation ( 8) can be discretized with standard Galerkin methods (e.g., Zienkiewicz and Taylor, 2000).Its weak form is where φ is a vector valued test function, and we have used Eq. ( 9) to eliminate the boundary integrals induced through integration by parts.Equation (10) can be calculated from Eq. (8) and does not require discretization.Equation ( 5) is hyperbolic and requires stabilization in order to suppress spurious oscillations.We use the streamline upwind Petrov-Galerkin (SUPG) method as a stabilization technique (Brooks and Hughes, 1982).SUPG have been used with success for the continuity equation in the ice sheet modeling context extensively (Morlighem et al., 2011;Larour et al., 2012).This case differs from previous work in that we are here attempting to solve for velocity rather than thickness.This means that velocity and thickness switch roles in the stabilization scheme; U is advected by the pseudo-velocity NH .The SUPG weak form is where λ is a test function that accommodates the influx or outflux Dirichlet boundary condition if so specified, V = {λ ∈ H 1 , λ| = 0}; τ is a mesh-dependent stabilization parameter given by and h is the element circumradius.We use linear Lagrange finite elements for discretization.The inclusion of this unusual stabilization term is key to achieving meaningful numerical solutions; without it, the solutions are plagued by non-physical oscillations.This instability is likely the reason that this approach has not been seen in the literature previously.

Application to the Greenland Ice Sheet
We apply this balance velocity approach to the Greenland Ice Sheet.We used the 1 km gridded GLAS/ICESat data set (DiMarzio et al., 2007) for surface elevations and a bed DEM from Bamber et al. (2001) for bed elevations.Annual average surface mass balance rates are derived from RACMO (Ettema et al., 2009).We assume that basal melt is small compared to surface mass balance, and neglect it.We also assume that the ∂ t H is negligible, or that the ice sheet is in balance.This is doubtless an incorrect assumption in some regions of the ice sheet, but although estimates for this field exist (e.g., Pritchard et al., 2009), it is not yet possible to determine what proportion of this signal is a result of ice dynamics, as opposed to other mechanisms such as firn densification that should not be included here.

Grid dependence
In order to assess the degree of grid dependence exhibited by this solution method, we start with a very coarse mesh, with an element circumradius of h = 32H and calculate balance velocity over progressively finer meshes, essentially halving the element size at each iteration, down to an element circumradius of h = H or 500 m, whichever is greater.We do this for smoothing lengths l ∈ {0, 4, 10, 15}.The difference between the coarse solution and progressively finer solutions is shown in Fig. 1.We see that for smoothing lengths of l ∈ {4, 10, 15} the norm of the difference between the refined and unrefined solutions stops changing with increasing refinement.When l = 0, the solution continues to change as the mesh becomes more refined.This indicates that incorporating a parameterization of horizontal plane stress in flow routing can overcome the tendency for the flow field to overconverge, even for very finely resolved meshes.

Flow direction smoothing radius
Theoretical results from Kamb and Echelmeyer (1986) suggest that the value of l for an ice sheet should fall between 4 and 10 ice thicknesses (although this range is based on temperate ice).Previous studies of horizontal plane stress coupling lengths for ice sheets typically indicate a value of l at the high end of this range (LeBrocq et al., 2006;Fricker et al., 2000), and often even higher (Testut et al., 2003;Joughin et al., 1997), in order to achieve heuristically good results.Identifying the optimal horizontal plane stress coupling length is also complicated by the fact that l should almost certainly be spatially variable.Nevertheless, we present balance velocities for l ∈ {4, 10, 15}, for a mesh size of h = H , which, based on results from the previous paragraph, should be a sufficiently small mesh size such that any smoothing of the flow is due to horizontal plane stress coupling rather than a lack of mesh detail.Figure 2 gives the balance velocity for the Greenland Ice Sheet at these length scales and mesh sizes, as well as the observed surface velocity.l = 4 produces an obviously overconvergent flow field, as evidenced by the abundance of discrete and overly narrow ice streams.l = 10 produces a better result, and we can see that most of the main flow features of the ice sheet are captured.velocity fields.The northeastern ice stream, while apparent, is less significant than indicated by observations.At l = 15, features begin to wash out, most notably the characteristic multi-pronged ice streams of Kangerdlugssuaq glacier.

Application: physics-based interpolation of the surface velocity
Here, we present an application of our new technique for determining the balance velocity.The application is one that relies on many thousands of evaluations of the continuity equation in order to numerically optimize model output.It is conceptually and mathematically similar to the technique described by Morlighem et al. (2011), but with balance thicknesses exchanged for balance velocities.For reasons of computational expense, our example could not be done without the advances presented earlier in this paper.
Geophysical data describing the cryosphere are in many cases incomplete or inconsistent with physical law.For example, take the surface velocity data of Joughin et al. (2010).They are characterized by large gaps in coverage and a highly variable structure in regions having low speed (less than ∼ 20 m a −1 ).Attributed to regions of high accumulation, high surface slopes, or incomplete satellite data, these problem regions frustrate many efforts that depend on complete coverage, or smoothness of the data.Applications affected might include inversion for basal traction (Morlighem et al., 2013;Brinkerhoff and Johnson, 2013) or calculations involving derivatives, such as resolving the stress balance (Van der Veen, 2013).
In order to use such data, practitioners are often required to smooth and/or interpolate the data.The fundamental procedure of interpolation is to generate a function that (1) is continuously valued over a given domain, (2) obeys some fundamental functional form between data points, and (3) adheres to observed values where data exist, with the understanding that such data are subject to error.Standard interpolation techniques often use polynomials as an interpolant.Physics-based interpolation differs by using solutions to the mass conservation partial differential equation (PDE) as the interpolating function.It is convenient to formulate this procedure as an optimization problem, which minimizes some measure of misfit between data values under the constraint of mass conservation.In particular, we are interested in minimizing the misfit between (possibly incomplete) velocity observations and balance velocities.This is expressed symbolically as where I is a misfit functional, F a functional that imposes continuity, and R a Tikhonov regularization used to impose a specified smoothness on the parameters.We depart from the previous notation by introducing balance velocity U m N , and observed velocity, u o , in order to keep the quantities being compared clear.We define the observed speed Finding the saddle point of Eq. ( 14) is known as PDE-constrained optimization.

Functional forms
I can take on a variety of forms.Here, we write a linear combination of least squares and log-least squares, or where e is the domain over which velocity observations exist.F is defined using a Lagrange multiplier λ to enforce where λ is a Lagrange multiplier.Note that this PDE constraint is still hyperbolic and requires the special numerical treatment defined previously in this paper.R is a Tikhonov regularization term that penalizes large gradients in the values of explanatory parameters; f ∈ P ≡ {F, H, N}.We adopt the following form: for i in the space of explanatory parameters.ξ i is a regularization parameter.

Solution method
Consider the following, simplified, form of the PDE constrained optimization problem: In practice we add a logarithm squared of the mismatch and regularization on each of the variables.However, this discussion neglects the terms to clarify the procedure that follows.Because each of the fields appearing in the continuity equation is measured in some way, we express the uncertainties in the measurements as follows: Thus, we state that the admissible spaces for the explanatory variables are defined by their assumed errors.Note that any choice within this range is assumed equally valid.
The mass conservation constraint, or forward model, is solved in two stages.First the directions of flow, N, are estimated from smoothed driving stress directions using the solution to Eq. ( 8).In regions where the direction of flow has been observed, N is replaced with the observed direction.The entire field is then smoothed to avoid large discontinuities on the boundaries between observed and estimated directions.The smoothing used takes the same form as Eq. ( 8).
Equation ( 12) is used to express the stabilized form of the forward model.The original problem, Eq. ( 18), can now be restated in terms of the stabilized PDE constraint as where the Lagrange multiplier plays the role of a test function.To simplify the mathematics to follow, identify λ = λ+τ ∇•N H λ and recover the original form stated in Eq. ( 18), the λ replacing λ.
We then take the first variation (formally a Gâteaux derivative) of I[U , H, F, N ; λ ] with respect to each of its parameters.For instance, the variation with respect to the thickness We note that a complete variation would have considered the error structure in observed speed, U o , as well, but given the large areas of missing data, we did not include this in the analysis.
After varying the functional with respect to all terms, the result is where we have ignored the dependence of λ on N and H .We also ignore variation with respect to U .Note that we can immediately identify individual terms specifying search directions (g i ) for each of the variables i ∈ {H, NF }, as well as the forward and adjoint models.
A few practical concerns arise, and are addressed as follows.
1. δN is ambiguous, because it is a vector.However, only one component of a normalized vector is independent; i.e., n 2 x + n 2 y = 1 can be solved for an unknown.In this example, the variation is always done on δn y .

Regularization is applied to each of the variables as
shown in Eq. ( 17).L-curve analysis suggests that values of ξ i between 10 7 and 10 8 are reasonable.In this example, all values were set to 10 7 .
3. In order to explain our approach, we present a simplified differentiation process.In practice, the complexity of the stabilization terms, the inclusion of the logarithmic mismatch function, and the introduction of regularization on the variables lead us to opt for automatic differentiation available through the FEniCS library that we use for finite element discretization (Logg et al., 2012).4. To make direct comparison of speeds, we need to estimate vertically averaged velocities from surface velocity (Joughin et al., 2010).To do so, we construct a function that approximates the role of deformation in the observed surface velocity.The function makes velocities above 120 m a −1 almost entirely due to sliding (surface velocity is the vertical average), and velocities below 25 m a −1 nearly entirely due to deformation (surface velocity is 80% of the vertical average).A smooth transition between the two end members is given by the logistic function 5. The weighting between logarithmic and linear terms in the misfit functional of Eq. ( 15) is set to be α = β = 0.5.Under this weighting choice, in fast flowing regions, the linear misfit is dominant, while in slow flowing regions, the logarithmic misfit is more important.

Errors and numerical details
For the ice thickness field, data are drawn from Bamber et al. (2013).These data represent the reduction and interpolation of hundreds of individual radar tracks into a map having complete coverage.Bamber et al. (2013) report errors along tracks of zero.Here, we use ±35 m along tracks, to reflect that there may be some error in the measurements.Off the tracks, we use the same values reported in Bamber et al. (2013).Ettema et al. (2009) provide surface mass balance, the only term used in our apparent mass balance, F .Because this is only part of the apparent mass balance, and because these The errors in the direction of the velocity reflect both differences from smoothed steepest descent where there are no velocity observations, as well as errors in the velocity observations.We assumed these to be in the range ±5 • .
All results were computed on an unstructured finite element mesh with an average spacing between nodes of 2 km.The optimization was done by using the gradients, g H , g F , g N , to drive the quasi-Newton bounded optimization technique, BFGS B (Nocedal and Wright, 2006).The optimization was terminated when the value of the objective function ceased to change appreciably, less than 0.5 % through searches along each of the gradients.

Results and discussion
We focus on results from the south of Greenland, where the velocity coverage is poor.Differences between observed and modeled speeds are shown in Figs. 3 and 4, respectively.The general structure of the observations is preserved, and the transitions between areas of no data and data are seamless.Much of the noisy signal that is apparent near the ice divide in the observed velocity is smoothed over in the interpolated data set.In the interpolated data there are numerous linear features that track the flow.These are not present in the original data and reflect the nature of the algorithm, which accumulates ice flux along flow lines.The interpolation scheme also diffuses the channelized nature of flow in the lower Jakobshavn area, and perhaps in other outlet glaciers as well.
Our approach also provides thickness and effective mass balance (F ) values that satisfy the continuity equation (Figs. 5 and 6).The changes made in order to uphold continuity are quite significant but still within the assumed error structure of the fields.In order to reproduce the observed speed in the outlet glaciers, thinner ice is required.This is due to the modeled velocity being too low; dividing the flux by a smaller thickness would increase the velocity.The bias toward slower ice could result from accumulation being too low, or velocity directions not being convergent enough.Apparent mass balance demonstrates that the search algorithm is utilizing this field to delimit streaming behavior by creating gradient in mass balance across the margins.
Changes in the direction of flow, N, were less significant due to the low errors assumed in this field.There was little systematic change in values and it is difficult to interpret how the optimization process impacts the values.
Moreover, the results demonstrate that it is difficult to uphold continuity and match the observed velocities.It is likely that the optimization is finding its way into a local minimum that is difficult to get out of.Once in this minimum, systematic changes in the surface mass balance and thickness fields are made in a manner that is not likely to be physically plausible, but is reasonable in terms of the stated error bounds.The technique presented here should improve in its utility as the coverage of fundamental data sets increases, and uncertainties decrease.Eventually, the minimum reached from the initial point will better correspond to a global rather than local one.One application of this approach will be to provide self-consistent initialization data for prognostic ice sheet modeling.Because the continuity is upheld by the data with a Lagrange multiplier, we are guaranteed that the combination of thickness, mass balance, and velocity produced by this method will not produce the strong gradients in model output produced by data in which flux divergence does not equal apparent mass balance (Perego et al., 2014).

Conclusions
We presented a novel numerical method for calculating the balance velocity of an ice sheet using the finite element method.This approach is an advance over classical routing techniques because it is not dependent on a heuristic routing algorithm and relies solely on a continuum conservation law and a theoretically motivated parameterization of flow directions.An unstructured grid easily allows for variable spatial resolutions.This method is made possible by two specific insights.First, flow directions that include horizontal plane stresses can be calculated by applying a spatially variable diffusion operator to the driving stress.Second, the balance ve-locity equations can be viewed as an advection equation with a pseudo-velocity field specified by thickness and flow direction, with velocity as the advected quantity.This problem is unstable.We use the streamline upwind Petrov-Galerkin method to make it tractable.
We applied this method to the Greenland Ice Sheet.Balance velocities were calculated over a number of different mesh resolutions, and we found that for given sufficient horizontal plane stress coupling distances, the solution shows grid independence.We also showed the balance velocity field calculated for theoretically justifiable smoothing lengths on detailed meshes.The resulting balance velocity compares favorably with a satellite-measured velocity field.
Additionally, we presented a numerical method that uses adjoint-based optimization to both fill data gaps and smooth spurious gradients present in an InSAR-derived velocity data set.This method is conceptually similar to Morlighem et al. (2011), but minimizes the misfit between balance velocities and observation, as opposed to thickness.We showed that we can find a balance velocity that matches InSAR data well, but does not possess gaps or strong gradients, while remaining within specified error bounds for input data fields.Despite this, we also find that upholding mass conservation requires surface mass balance and thickness fields that are distinctly less smooth than those reported.Regardless, this PDE-constrained interpolation technique promises to be a useful tool for providing smooth and continuous velocity data that conform well to observations.

Figure 2 .
Figure 2. Balance velocity solution for a mesh size of h = H and l ∈ {4, 10, 15} as well as InSAR surface velocities.

Figure 3 .
Figure 3. Surface speed of ice from observations reported in Joughin et al. (2010).

Figure 4 .
Figure 4. Final surface speeds, computed through the optimization of the speed constrained by the continuity equation described in this paper.

Figure 5 .
Figure 5. Differences between the ice thicknesses reported in Bamber et al. (2013) and the thicknesses found at the end of the optimization procedure.

Figure 6 .
Figure 6.Apparent surface mass balance determined at the end of the optimization procedure.