Data Assimilation and Prognostic Whole Ice Sheet Modelling with the Variationally Derived, Higher Order, Open Source, and Fully Parallel Ice Sheet Model VarGlaS

We introduce a novel, higher order, finite element ice sheet model called VarGlaS (Variational Glacier Simulator), which is built on the finite element framework FEniCS. Contrary to standard procedure in ice sheet modelling, VarGlaS formulates ice sheet motion as the minimization of an energy functional, conferring advantages such as a consistent platform for making numerical approximations, a coherent relationship between motion and heat generation, and implicit boundary treatment. VarGlaS also solves the equations of enthalpy rather than temperature, avoiding the solution of a contact problem. Rather than include a lengthy model spinup procedure, VarGlaS possesses an automated framework for model inversion. These capabilities are brought to bear on several benchmark problems in ice sheet modelling, as well as a 500 yr simulation of the Greenland ice sheet at high resolution. VarGlaS performs well in benchmarking experiments and, given a constant climate and a 100 yr relaxation period, predicts a mass evolution of the Greenland ice sheet that matches present-day observations of mass loss. VarGlaS predicts a thinning in the interior and thickening of the margins of the ice sheet.


Introduction
Models have become an important tool in the study of glacier and ice sheet physics, with applications to the prediction of cryosphere/climate interactions, sea level rise, and fundamental questions of ice dynamics.In recent years, while conceptual and theoretical advances in the development of ice sheet models have been made, computational constraints limited practical ice sheet models to low-order asymptotic ap-proximations of ice physics.Early models (and many modern ones) relied upon the shallow ice approximation, sacrificing accuracy in regions of complex flow for efficiency at large scales (Hutter, 1983;Schäfer et al., 2008).Models including higher order physics were also used, but the increased computational demand made operating at high resolutions infeasible.Additionally, many models were based on finite difference schemes, which made variable resolutions over different regions in the computational domain difficult.Recent increases in computing power and the availability of parallel libraries, in tandem with the finite element method and unstructured meshes, have made possible the use of higher order physical approximations coupled with the high resolution necessary to resolve fine-scale features of glacier flow.The ice sheet modelling community has been quick to take advantage of these advances, for example, in Larour et al. (2012), Seddik et al. (2012), Gillet-Chaulet et al. (2012), Leng et al. (2012), and Bueler and Brown (2009), among others.In this work, we present a new thermomechanically coupled, prognostic ice sheet model called VarGlaS (Variational Glacier Simulator).As stated above, there are already a few examples of so-called "next-generation" ice sheet models in existence, but this model differs in implementation strategy in several critical regards.Namely, these are the use a of a variational principle in the model formulation, the solution of an enthalpy equation, the use of a kinematic boundary condition, and the use of automatic differentiation for data assimilation.The implementation of these features is streamlined by building the model upon the framework offered by the open-source FEniCS project (Logg et al., 2012).FEniCS offers access to a variety of finite element solvers and tools, as well as a very deep integration of automatic differentiation, which makes many of the more unique capabilities of VarGlaS possible.
VarGlaS treats the solution to the momentum balance (Stokes' equations) as the minimization of an energy functional.The existence of a variational principle for nonlinear Stokes' flow was shown by Bird (1960).When applied to ice sheets, the method consists in minimizing the dissipation of gravitational potential energy by viscous and frictional heat generation.This approach has been suggested by Schoof (2006), Dukowicz (2012), andBassis (2010).This type of treatment stands in contrast to the heretofore standard treatment of the velocity field, which is to explicitly account for the balance of viscous stresses and forcing by gravity, as is done in most other ice sheet models (e.g.Larour et al., 2012;Seddik et al., 2012;Leng et al., 2012;Bueler and Brown, 2009;Rutt et al., 2009).Viewing the problem as a variational minimization problem confers a number of advantages.The most important advantage is that the momentum balance is uniformly derived from a single scalar conservation statement.When approximations to the physics are made, they are made to the scalar quantity, and these changes are automatically propagated through the rest of the model.This is particularly useful when automatic differentiation is available.As noted above, the FEniCS software offers strong automatic differentiation capabilities, with particularly strong support for the calculation of Gâteaux derivatives, which is the necessary operation for deriving equations of motion from a variational principle (Dukowicz, 2012).The procedure of generating the code for a new approximation to the Stokes' equations is as simple as making a change to the variational principle.This makes an extension of the model to different asymptotic approximations straightforward.Other potential advantages are that the variational principle is coordinateindependent, streamlining the transition to a curvilinear or geographic coordinate system.Boundary conditions are also implicitly defined within the variational principle, meaning that the often complex process of imposing boundary conditions is simplified.The use of a variational principle also confers several computational advantages.For example, the operators derived from the variational principle are guaranteed to be at least symmetric and semi-definite.Also, these operators are already in the appropriate form for use with the finite element method, requiring no further manipulation.
Treatment of the energy balance by VarGlaS also differs from standard methods.Typically, temperature is the variable of interest in the energy balance (Larour et al., 2012;Seddik et al., 2012;Greve and Hutter, 1995;Rutt et al., 2009;Pattyn, 2003).Computing the temperature field is a contact problem where the temperature must be constrained to remain below the phase boundary.Different methods have been employed to enforce this constraint, such as treating temperate and cold ice as two separate fluids (Greve and Hutter, 1995), manipulating heat sources and sinks such that heat sources are applied to the temperature equation when below the melting point and to calculate a melt rate when the ice is at the pres-sure melting point (Rutt et al., 2009), or solving a contact problem (Zwinger et al., 2007).To avoid inconsistency and heuristics, we eschew the temperature formulation in favour of an enthalpy treatment that tracks total internal energy density rather than sensible heat.This eliminates the need for special numerical treatment of the cold-temperate transition surface, at the expense of introducing a nonlinearity in energy diffusion (Aschwanden et al., 2012).While enthalpy is more straightforward computationally, the temperature field is still necessary for the computation of ice rheology and for interpretation of model results.Temperature can be recovered from enthalpy in a straightforward way through a bijection between enthalpy, temperature, and water content.
VarGlaS is equipped with a kinematic boundary condition that allows for the evolution of ice geometry.Many models calculate change in surface elevation as the flux divergence of the vertically averaged horizontal velocity field (Larour et al., 2012;Rutt et al., 2009;Bueler and Brown, 2009).We treat it as the advection of the ice surface by the surface velocity field, as is done in Seddik et al. (2012) and Leng et al. (2012).The two forms are equivalent, and both are numerically unstable.VarGlaS, similar to other modern finite element models, uses streamline upwind Petrov-Galerkin finite elements to stabilize the free surface problem (Larour et al., 2012;Seddik et al., 2012).Additionally, transitions between glaciated and ice-free regions of the model domain produce numerical instabilities which need to be addressed with methods beyond upwinding in order to maintain higher order accuracy.To this end, VarGlaS also introduces a discontinuity capturing scheme to maintain stability in the presence of large gradients.VarGlaS uses a unique and fully explicit total variation diminishing Runge-Kutta scheme for time discretization.This method guarantees that no new spurious extrema are generated by the time-stepping scheme, and provides second order in time accuracy without necessitating a complex and expensive coupling between the momentum balance and time evolution.
For both diagnostic and prognostic modelling, it is important to constrain the model to match observed values of state variables as closely as possible.To this end, VarGlaS possesses tools for data assimilation, again based on the automatic differentiation capabilities of FEniCS.With this functionality in place, calculating the adjoint state of a model and the gradient of a given objective function with respect to arbitrary parameters under the constraint that a forward model be satisfied is simple and automated.We use this capability in two complementary ways.First, we have been able to invert for sliding velocity (or basal traction) to best match surface velocities, a classic problem in glacier modelling (MacAyeal, 1993;Goldberg and Sergienko, 2011;Larour et al., 2005;Gudmundsson and Raymond, 2008;Morlighem et al., 2010;Brinkerhoff et al., 2011).In contrast to many implementations of data assimilation in ice sheet modelling (MacAyeal, 1993;Morlighem et al., 2010), VarGlaS uses an unsimplified adjoint (Goldberg and Sergienko, 2011;Brinkerhoff et al., The Cryosphere, 7, 1161-1184, 2013 www.the-cryosphere.net/7/1161/2013/2011).We can also minimize the total imbalance in mass continuity with respect to basal topography (e.g.Morlighem et al., 2011), although this procedure is presented in a separate paper (Johnson et al., 2013).These methods are important for long-term simulations, as velocity assimilation reproduces a plausible velocity that is of leading order importance in calculating surface rates of change, while mass conservation assimilation helps to eliminate strong transients resulting from estimates of basal topography that are incoherent relative to the ice physics and observed data.
The paper is structured as follows.Section 2.1 discusses the continuum mechanical formulation of the model physics.Section 2.2 deals with the numerical implementation of the model physics, and the difficulties arising from their discretization.Section 3 involves the application of the model to a few numerical experiments including well-known benchmarks involving idealized geometry, as well as the entire Greenland ice sheet.Finally, in Sect. 4 we discuss how the combination of advances has allowed VarGlaS to successfully simulate both numerical benchmarks and continentalscale ice dynamics, as well as the fundamental limitations of the model and what directions should be explored in order to overcome these.

Model
VarGlaS solves for the three-dimensional ice sheet velocity, enthalpy, and geometry through time.All three of these variables are strongly coupled.We first present the continuum formulation of ice sheet physics, followed the numerical treatment for each.

Momentum balance
Our development of a variational principle for the momentum balance largely follows Dukowicz (2012).The variational principle for a power law rheology with linear basal sliding under the constraints of incompressibility and bed impenetrability is where u is the ice velocity and ˙ the rate of strain tensor, P the pressure, η(˙ 2 ) the strain-rate-dependent ice viscosity, g the gravitational vector, β 2 the basal sliding coefficient, h the ice thickness, r a factor determining the relationship between basal traction and thickness, m a factor determining the nonlinearity of the frictional term, n the outward normal vector, and P e , the environmental pressure, either atmospheric or water (see Table 1 for a complete listing of symbols used in this manuscript).The expression is minimized over the ice domain with boundaries , where B is the grounded portion of the ice sheet.E is the non-grounded portion of the ice sheet, and is defined by Each of the additive terms in Eq. ( 1) has a specific meaning.
Terms integrated from left to right over are viscous dissipation, gravitational potential energy, and the incompressibility constraint, respectively.Terms under the first boundary integral are frictional heat dissipation and the impenetrability constraint.Note that basal traction is scaled by thickness compared to the standard form of the sliding law.For r = 1, this eliminates the covariance between basal traction and pressure evident in previously computed basal traction fields (e.g.Larour et al., 2012).For r = 0, it is obvious that this form reduces to the standard sliding law.The term under the second boundary integral is the pressure at the ice-air or ice-water interface.Note that, at this point, VarGlaS does not allow these boundary domains to change, implying a static grounding line.The constitutive relationship for ice given by Glen (1955) gives a viscosity of where ˙ 2 is defined to be the square of the second invariant of the strain rate tensor, and b(T , ω) is a temperature-and water-content-dependent rate factor: Here, E is an enhancement factor, with a(T , ω), Q * (T ) and R parameters, and T * is temperature-corrected for pressure melting point dependence.The traditional momentum balance form of the Stokes' equations can be recovered (in weak form) by taking the variation of Eq. ( 1).This is the functional that VarGlaS minimizes in order to solve the Stokes' problem.
The Stokes' functional is a relatively complete statement of ice physics (the only assumptions being negligible inertial terms), but it includes four degrees of freedom per computational node (three velocity components and pressure) and is a saddle point problem due to the presence of the Lagrange multiplier pressure terms.A considerable simplification can be made to the Stokes' functional by expressing vertical velocities in terms of horizontal ones through the incompressibility and bed impenetrability constraints, that is www.the-cryosphere.net/7/1161/2013/The Cryosphere, 7, 1161-1184, 2013 with boundary condition Substitution of these expressions into A yields an unconstrained and positive definite integro-differential functional which is equivalent to Eq. (1).However, the integral terms that result from the vertical integration of the mass conservation equation are undesirable.Standard methods for the numerical solution of partial differential equations (PDEs) are not equipped to handle integral terms of this type, so we seek a simplification that eliminates them.In order to derive the functional associated with the so-called "first-order" equations of ice sheet motion (Blatter, 1995;Pattyn, 2003), two assumptions must be made.First, bed slopes are small, which is also equivalent to assuming cryostatic pressure.Second, horizontal gradients of vertical velocity are small compared to other components of the strain rate tensor.This eliminates vertical velocity terms in the interior of the ice, as well as at the surface and grounded basal boundaries.After these assumptions and some manipulation, the first-order functional is where u is the velocity vector in the horizontal directions, S the elevation of the ice surface, and ˙ 2 1 the first-order strain rate tensor given by Pattyn (2003) and Dukowicz (2012).The first-order assumption does not allow a decoupling of the vertical velocity where the ice is not grounded, which is a direct consequence of assuming that vertical resistive stresses are negligible.Thus we assume, for the first-order model, that any ice beyond the grounding line instantly calves and is no longer considered.Thus we drop the last term in Eq. ( 6).Since the first-order equations are only associated with horizontal velocity components, this formulation yields significant computational savings, as well as desirable numerical properties such as guaranteed positive definiteness.Vertical velocity for the grounded ice sheet is recovered from Eqs. (4) and (5).

Enthalpy
VarGlaS uses an enthalpy formulation of the energy balance (Aschwanden et al., 2012).Enthalpy methods track total internal energy, rather than sensible heat, which corresponds bijectively to temperature for ice below the pressure melting point, and to water content for ice at the pressure melting point.The enthalpy equation is a typical advection-diffusion equation with a nonlinear diffusivity: where H is enthalpy, ρ ice density, and Q strain heat generated by viscous dissipation, given by the dissipative term in the momentum balance functional.κ is an enthalpydependent diffusivity given by where k is the thermal conductivity of cold ice and C p is heat capacity.ν is the diffusivity of enthalpy in temperate ice, and can also be thought of as a parameterization of the subgridscale intraglacial flow of liquid water.It is not clear what the value of ν should be.Both Hutter (1982) and Aschwanden et al. (2012) have suggested that it be a function of both water content and gravity, but intra-glacial liquid modelling is beyond the scope of this work, and we set this value to ν k C p .This implies that heat does not move diffusively within temperate ice, and that any heat generation immediately goes towards melting.The definitions for cold and temperate ice are as follows: where H m is the pressure melting point expressed in enthalpy, and C w is the heat capacity of liquid water, γ the dependence of the melting point on pressure, T 0 the triple point temperature of water, and L the latent heat of fusion for water.
At the ice surface, we specify a Dirichlet boundary condition corresponding to surface temperature.At the basal boundary, we apply the Neumann boundary condition: where q g is geothermal heat flux, assumed known, q f frictional heat generated by basal sliding, and M b the basal melt rate.Frictional heat is given by the frictional term in the momentum balance functional.Note that in temperate ice, where κ(H ) is nearly zero (no diffusion), this relation defines the basal melt rate.In cold ice, a value must be specified for the basal melt rate (which can be negative to account for basal freeze-on).We usually take this to be zero.Enthalpy is uniquely related to temperature and liquid water in the following way: where ω is fractional water content, while H m and T m are the pressure melting points expressed in enthalpy and temperature, respectively.

Dynamic boundaries
The ice sheet geometry evolves over time according to the kinematic boundary condition where ȧ is the accumulation rate.

Marine outlet treatment
For the Stokes' model, VarGlaS currently treats the grounding line in the simplest way possible, which is to keep its location fixed.At this point, ice can not become ungrounded.For transient runs, the geometry of calving fronts is fixed so that mass loss due to calving is effectively proportional to the velocity at the calving front.At the scale of outlet glaciers, this is a limitation and is a major priority in ongoing development, but the physics of the shelf is treated.As mentioned, the first-order approximation is not conducive to the treatment of a coupled sheet-shelf system, and so ice flowing past a pre-specified grounding line is assumed to calve immediately.

Numerical methods
In the following sections, we discuss how the continuum mechanical equations discussed above are discretized in order to be made computationally tractable.

Finite element discretization using FEniCS
VarGlaS is built upon the finite element package FEniCS (Logg et al., 2012).FEniCS is a powerful development environment for performing finite element modelling, including strong support for symbolic automatic differentiation, native parallel support and parallel interface with linear algebra solvers such as PETSc (Balay et al., 2012) and Trillinos (Heroux et al., 2005), and automatic code generation and compilation for compiled performance from an interpreted language interface.The Python scripting environment makes the generation and linking of new code straightforward.We find that this interface provides a level of extensibility that makes VarGlaS promising for distributed development and rapid prototyping of models for additional components of the cryosphere.
FEniCS has a large library of finite elements available.We use only one, the continuous, linear Lagrange finite element, defined over an unstructured triangular mesh.This choice of element is unstable for advection-dominated equations, such as the kinematic boundary condition and the enthalpy equation (in most cases), as well as for Stokes' equations due to the pressure term.We cover stabilization procedures in the following sections.
The velocity field and enthalpy equations are both nonlinear.These are each solved by using a relaxed Newton's method (e.g.Deuflhard, 2004) with a Jacobian calculated by automatic differentiation: where U i is the solution vector at the i-th iteration, U a solution update, J the Jacobian matrix, and F the system of nonlinear equations.R is a relaxation parameter that arbitrarily shortens the step size in order to improve numerical stability.The amount of damping required is specific to the problem, but we find that a relaxation parameter between 0.7 and 1.0 is typically sufficient to achieve convergence.We specify both a relative and absolute tolerance as convergence criteria for Newton's method.The solution is considered converged if the L ∞ norm of U is less than 10 −6 m a −1 or the L ∞ norm of U U n is less than 10 −3 .In order to resolve the coupling between enthalpy and velocity, we use a fixed point iteration.Each of these nonlinear equations is solved independently, and the result is iteratively used as input in calculating the other variable.Convergence is assumed when both the velocity and temperature updates are less than 10 −6 m a −1 and 10 −6 K, respectively.

Mesh refinement
The model domain is discretized using a tetrahedral mesh which is unstructured in the horizontal dimensions, and structured in the vertical.In order to equidistribute discretization error, we use the anisotropic error metric according to Habashi et al. (2000): where e(c) is a cell-wise error estimate, E a given mesh cell, x i an edge in E and M a metric tensor, in this case defined by V and are the respective eigenvectors and eigenvalues of the Hessian matrix of the field over which error is to be equidistributed (Habashi et al., 2000).For all the meshes presented forthwith, we use the Hessian of an observed velocity norm (either observed or modelled) in calculating error metrics.A discrete approximation for each component of the Hessian matrix is obtained iteratively for each level of mesh refinement by solving the variational problem where Hess ij denotes the components of the Hessian and U is the surface speed.With error estimates in hand, we isotropically refine all cells that are above a specified proportion of the average error.In order to account for the directional nature of the velocity field, we incorporate anisotropy by using Gauss-Seidel iterations (Habashi et al., 2000) to solve approximately an elasticity problem, with computed edge errors as "spring constants".This mixed isotropic-anisotropic technique yields high quality and efficient meshes with both the structural simplicity of isotropic refinement, as well as the better error to mesh size ratio of anisotropic techniques.An example of a mesh created with this method is shown in Fig. 1.Note that this algorithm does not allow mesh coarsening, so the refinement procedure is initialized with a coarse mesh.

Data assimilation and regularization
Many physical quantities of leading order relevance to glacier and ice sheet flow are either practically impossible to collect, or are point measurements which cannot generally be extrapolated to a broader spatial context.Examples of the former include historic variables such as a detailed record of surface temperature or ice impurity content at deposition.Examples of the latter include basal water pressure, basal temperature, enhancement factors, and geothermal heat flux.A particularly important parameter which must usually be estimated is the coefficient of basal traction, which relates basal shear stress to sliding velocity.In many cases, sliding makes up nearly all of a glacier's surface velocity (e.g.Weis et al., 1999).Thus, any model that wishes to reproduce plausible velocity and thermal structures must parameterize traction.The availability of widespread surface velocity data, and the conceptually simple relationship between surface and bed velocities have made the inversion of surface velocities for basal traction a popular choice for performing this parameterization (MacAyeal, 1993;Goldberg and Sergienko, 2011;Larour et al., 2005;Gudmundsson and Raymond, 2008;Morlighem et al., 2010;Brinkerhoff et al., 2011).
We have implemented basal traction inversion in VarGlaS using a partial differential equation-constrained optimization procedure.In the following, we illustrate the method using surface velocity in the cost functional, and basal traction as the control variable, but the procedure is analogous to any choice of objective function of control variable.The fundamental concept behind this method is to define a scalar objective function, to calculate its gradient, and to use standard minimization techniques to find the minimum.We use a general form for the definition of the cost functional I .Examples include a linear cost functional: or a logarithmic one We require that the velocity field obtained by this functional satisfies the equations of motion by imposing the momentum equations via a Lagrange multiplier: where δ implies the first variation operator, and A is one of the energy functionals defined in Sect.2.1.1.λ is a Lagrange multiplier used to enforce the forward model as a constraint.
Taking the variation of I with respect to u, β 2 , and λ yields, respectively, a forward model, an adjoint model, and an expression for the gradient of the objective function with respect to β 2 , which is expressed in terms of u and λ.Note that no simplifying assumptions about the nature of the forward model are made.In particular, when possible we use the full adjoint, calculated via automatic differentiation, rather than making the assumption that the viscosity does not depend on u, as is done in many inversion procedures (e.g.Goldberg and Sergienko, 2011;Larour et al., 2012).In the case where strong mismatches between the modelled and surface velocity exist, stability of the inversion numerics necessitates fixing the viscosity, and using an incomplete adjoint, as in Goldberg and Sergienko (2011), but only for the first few iterations.
In order to impose a minimum bound on the smoothness of the solution, we add a Tikhonov regularization term, which penalizes wiggles in the control variable.This regularization is of the form where α is a positive weighting tensor.This value is different for different objective functions and different model domains, and the means to determine its value is not obvious.In the glaciological literature, L-curve analysis is a popular choice that involves solving the inverse problem with increasing regularization until a notable change in the rate of objective function increase is found (e.g.Gillet-Chaulet et al., 2012).It is assumed that this inflexion point represents the optimal balance between regularization and minimization of the objective function.We prefer to use heuristics, assuming (and imposing) a maximum rate of spatial variability in the basal traction field based on factors such as ice thickness.Note that applying Tikhonov regularization on the gradient in this way is equivalent to applying an anisotropic diffusion operator to the control variable.
With a means of efficiently computing the objective function and its gradient with respect to the control variable in hand, we can use any number of optimization algorithms to minimize I.We use the quasi-Newton algorithm L BFGS B (Nocedal and Wright, 2000).We used a parallel implementation of the L BFGS B algorithm derived from that appearing in the dolfin-adjoint library (Farrell et al., 2012).Termination criterion for the optimization routine is essentially heuristic, with the optimization procedure terminating upon the objective function reaching a valley.The definition of a reliable convergence criterion is a subject of ongoing research, with the methods of Habermann et al. (2012) particularly promising.

Time evolution
We use two different algorithms for the discretization of time.
For the enthalpy equation, we use a semi-implicit Crank-Nicholson time-stepping scheme, with an arbitrary Lagrangian-Eulerian (ALE) treatment of the convective velocities in order to compensate for the moving mesh (Donea et al., 2005) (ALE operates by subtracting the mesh velocity from the fluid velocity).This semi-implicit method is acceptable because of the relatively minor nonlinear coupling between enthalpy and the velocity field; linearization is achieved by calculating the nonlinear dependence with the values of the previous time step.Crank-Nicholson provides second-order accuracy in time, provided that the Courant-Friedrich-Lewy criterion is satisfied, where D is the vector of cell dimensions in each direction and C is a constant, usually taken to be 1 2 .Note that this condition can be restrictive in some outlet glaciers, where the combination of 1 km-scale spatial resolution and > 1 km a −1 -scale velocities requires time steps on the order of a month.
For the free surface equation, the nonlinear coupling between velocity and surface elevation is stronger, so linearization is not a desirable option, and solving the full nonlinear problem implicitly is inefficient.We instead choose to use a fully explicit total variation diminishing Runge-Kutta (TVD-RK)-type scheme of second or third order (Gottlieb and Shu, 1998).The total variation diminishing property implies that no spurious oscillations should be created as a result of the time discretization.This must be coupled with a non-oscillatory spatial discretization (such as streamline upwinding) in order to maintain a stable solution.The TVD property is important in suppressing spurious oscillations near sharp margins (such as those characterizing glacial termini).Being a Runge-Kutta method, we must solve the momentum balance once for each order of accuracy, but we find that the added stability and accuracy of using a higher order explicit method makes this increase in computational overhead worthwhile.

Stabilization
Both the enthalpy and free surface equations are hyperbolic, and the standard centred Galerkin finite element method gives rise to spurious oscillations.In order to provide stabilization, we apply streamline upwind Petrov-Galerkin (SUPG) methods (Brooks and Hughes, 1982).For the enthalpy equation, this consists of adding an additional diffusion term of the form where K is a tensor valued diffusivity defined by where D is a cell size metric.Alternatively, we can view this stabilization as using skewed finite element test functions to weight the advective portion of the governing equation.Since the time derivative is implicitly defined, there is no need to apply upwind weighting to the time derivative or source terms, and because linear elements are used, applying this weighting to the diffusive component would necessitate second derivatives of test functions, which are always zero for linear elements.
For the free surface equation, a similar procedure is used, where modified Galerkin test functions are used to discretize the equations.Note that in this case, due to the explicit time-stepping scheme, the augmented weighting function is applied consistently to the entire residual, including the time derivative.In addition to streamline upwinding, we apply a shock-capturing artificial viscosity in order to smooth the sharp discontinuities that occur at the ice boundaries, where the model domain switches from ice to ice-free regimes (Donea and Huerta, 2003).This additional term is given by where K shock is the nonlinear residual-dependent scalar Here, R is the residual of the original free surface equation.
For the Stokes' equations to remain stable, it is necessary either to satisfy or to circumvent the Ladyzhenskaya-Babuska-Brezzi (LBB) condition.The typical way of doing this is to use a mixed second order in velocity, and first order in pressure finite element (the Taylor-Hood element).While VarGlaS has the capacity to use this formulation, we find that the additional degrees of freedom introduced by the higher order elements lead to an unacceptable loss of computational performance.Instead, we circumvent this condition by using a Galerkin least squares (GLS) formulation of the Stokes' functional where τ gls is a stabilization parameter (Baiocchi et al., 1993).For a linear Stokes' problem, the usual value for τ gls is Since τ gls is a function of the ice viscosity, τ gls should rightly be nonlinear.However, we have found through experimentation that ignoring the strain rate dependence of the viscous term yields acceptable results and much better numerical stability.Thus we use where η is some linear estimate of η.We have found η = 10 3 × b(T ) to yield an appropriate blend of fidelity to the governing equations and stabilization.Note that this has the effect of adding a diffusive term over pressure to the conservation of mass equation.

Parallelism
VarGlaS has been developed to take full advantage of the innate parallel capabilities of PETSc (Balay et al., 2012) and FEniCS (Logg et al., 2012), on which it is built.All computationally intensive components of the model are compatible with parallel usage, such as the nonlinear solvers, time stepping, and optimization.VarGlaS exhibits good scaling between 1 and 16 cores, the largest cluster to which the authors have access.Parallel efficiency for a nonlinear solution of the first-order equations over all of Greenland for over a million degrees of freedom is shown in Fig. 2.

Verification with the method of manufactured solutions
We use the method of manufactured solutions (MMS, Salari and Knupp, 2000) to construct analytical solutions to the momentum balance and surface evolution equations in order to determine the extent to which our numerical solutions reproduce the exact results.This procedure is known as verification, which ensures that the given equations are being solved correctly and consistently (this is held in contrast to validation, which is the procedure of ensuring that the appropriate equations for a given physical scenario are being solved).The use of MMS in the context of glaciological models has been previously developed by Bueler et al. (2005), with more recent results by Sargent and Fastook (2010) and Leng et al. (2013).We adopt a similar approach to these works, but like Leng et al. (2013) and unlike Bueler et al. (2005), we explicitly forego verifying thermomechanical coupling.The main idea behind MMS is to select an arbitrary solution, and then construct a source function that produces that solution.To that end, we select the following functions for the components of the ice velocity and free surface: Also, we have that B(x) = S(x, t = ∞) − 10 3 , L = 10 4 , and Amp = 500, which corresponds to the geometric set-up of We begin by verifying the computation of the velocity field, which is a pseudo-steady-state calculation; the velocity field's only dependence on time is through the problem geometry.This implies that we can verify the velocity field independent of time, and for this procedure set t = 0.
We seek vector-valued functions F(x, t) and G(x, t) such that the functions û(x, t) that we have selected minimize the functional where diag(G) is a diagonal matrix with the components of g on the diagonal.Note that we have eliminated the Lagrange multiplier P from the functional, since û is constructed to be divergence-free.Taking the variation with respect to û yields the Stokes' equations with arbitrary forcings within the model domain and at the boundary.Since û is analytical, we can invert for F and G.The procedure is the same for finding the forcing functions for the first-order equations, but we use the first-order functional Eq. ( 6).
While the functions used in the MMS are arbitrary, we note a few features of our selection.We have eliminated P from the formulation by a specific choice of velocity field.
One may view the pressure as itself being an arbitrary forcing function, and in this case it becomes part of F and G. Also, we have chosen v = 0 because it greatly simplifies the analytical source terms, and also because it allows us to carry out an important test of rotational symmetry by swapping the definitions of û and v (and their spatial arguments).
With the analytic forcing functions F and G in hand, as well as their corresponding analytical solutions, we can solve the discrete model over progressively finer meshes and evaluate whether the error in the numerical solution tends towards zero as the discrete domain better approximates the continuum.In practice, computational intensity limits the maximum level of refinement, and it is sufficient to show that a discretization converges towards the analytical solution at close to the theoretical order of accuracy, which for linear finite elements is O( x 2 ) (Salari and Knupp, 2000).We solved both of these problems over meshes ranging from x = 5×10 3 m to x = 40 m, two orders of magnitude.The resulting error profiles are shown in Fig. 3.The slope of the relative L 2 norm for both approximations is close to m = 2 in logarithmic space, which is the theoretically correct result for linear elements.This confirms that our discretization procedure shows theoretical order of accuracy convergence and thus correctly solves the governing equations for the velocity field.
Verification for our free surface evolution scheme proceeds similarly, except that refinements to the computational domain are carried out both in time and space.We seek a compensatory accumulation field ċ, such that the kinematic boundary condition holds for the pre-defined, analytical functions Ŝ and û.Again, it is a simple matter to determine a closed form expression for ċ.We use the same analytical functions as defined above, but allow t to vary.To test the convergence of our time-stepping scheme, we progressively refine t and run the model to t = τ using the compensatory accumulation ċ in place of ȧ in Eq. ( 14) and evaluate the error in surface elevation at that time.Note that because we are using an explicit time-stepping algorithm, the Courant-Friedrich-Lewy (CFL) condition applies, limiting the minimum cell size that can be used for a given time step.Simultaneously, we require a certain spatial accuracy such that errors in the temporal discretization are not overwhelmed by errors resulting from the spatial discretization.Thus, as we refine the time step, we also refine the spatial resolution such that the Courant number C = 0.1 for all simulations.This ensures that instability due to a violation of the CFL condition does not occur, and that an appropriately fine spatial resolution is used for a given time step.We refine between t = 50 a and t = 0.1 a, which correspond to spatial resolutions between x = 10 4 m and x = 10 -1 10 0 10 1 10 2 ∆t (a) 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1  3. The theoretical behaviour of our RK2 algorithm implies a single-step error of O( t 3 ) and a cumulative error over time of O( t 2 ) (Gottlieb and Shu, 1998).These errors correspond to a slope of m = 3 and m = 2 respectively, and we expect the slope of the error in surface height at time τ with respect to t to be between these two values.Indeed, we find that m = 2.8 for both the L 2 and L ∞ norms.This result shows that our timestepping scheme possesses theoretical convergence properties, and can be considered verified.
We noted earlier that our choice of velocity fields could be used to test symmetry.We performed both of the above convergence tests for a swapped û and v (with associated swapped coordinates x and y), and produced convergence results that were identical to machine precision.Thus we conclude that VarGlaS possesses the appropriate properties of rotational symmetry.

Numerical experiments
In order to assess the correctness and efficiency of model, we apply it to a number of well-known ice sheet modelling benchmark experiments, before turning it towards a largescale data assimilation and prognostic time-stepping simulation for the whole Greenland ice sheet.

ISMIP-HOM
The Ice Sheet Model Intercomparison Project-Higher Order Model (ISMIP-HOM) benchmarks are a widely used test of higher order model capabilities (Pattyn et al., 2008).In order to verify our model performance, we run ISMIP-HOM tests A, B, C, D, and F using both the first-order and Stokes' equations for the momentum balance.All meshes for the ISMIP-HOM experiments used structured grids with 20 elements in each horizontal dimension, and 10 elements in the vertical dimension, for 4851 degrees of freedom.
ISMIP-HOM A and ISMIP-HOM B simulate steady ice flow with no basal slip over a sinusoidally varying bed with periodic boundary conditions.B varies from A in that the basal topography does not vary in the transverse direction, whereas in A it does.Note that despite the fact that B is uniform in one spatial dimension, we still treat it as a threedimensional problem and use the same discretization and boundary conditions (periodic, namely) as for A. Figures 4  and 5 show the simulated surface velocity for all length scales outlined in the benchmarks.
ISMIP-HOM C and ISMIP-HOM D simulate steady ice flow with sinusoidally varying basal traction over a flat bed with periodic boundary conditions.Again, D varies from C in that only C is transversely varying.Figures 6 and 7 show the simulated surface velocity for all length scales outlined in the benchmark.This experiment specifies that r = 0 in the sliding law.
Experiments A-D were run in serial with a 2.2 GHz laptop computer.Runs using the first-order approximation took approximately 10 s, while runs using the Stokes' approximation took around 60 s.
After running the ISMIP-HOM C experiment forward, we have in hand the velocity field predicted by the model for a given basal traction.We used this as an opportunity to test the inverse capabilities of the model, and to invert for a known basal traction.We performed the inversion using the first-order approximation.Starting from an initial guess of a uniform basal traction of 10 3 Pa a m −1 , we allow the inverse model to predict the basal traction field that produces the velocity field of the forward model (which we know to be a sinusoid).In assessing misfits, we used the linear cost functional Eq. ( 20), because the velocity scale varies relatively little over the domain.This problem was unregularized, so α = 0. Figure 8 shows both the rate of convergence as well as the "observed" and modelled basal tractions and surface velocities along the L 4 transect of the ISMIP-HOM model domain using the results for the first-order approximation at the L = 10 km length scale.Results for other length scales are similar.
ISMIP-HOM F simulates unsteady ice flow over a Gaussian bump with periodic boundary conditions under slip and non-slip conditions, and evaluates the surface geometry and velocity as they relax to steady state.Experiment F also mandates a linear rheology, so n = 1. Figure 9 shows the simulated surface velocities and elevations for the slip and nonslip cases, using a time step of t = 1 a.These experiments were run on a 4-core machine with 8 GB of memory and 2.2 GHz per processor.Using the first-order approximation, the 500 a run took approximately 50 min.Using the Stokes' approximation, the run took 5.5 h.

EISMINT II
The European Ice Sheet Model INTercomparison II (EIS-MINT II) benchmarks are an older set of benchmarks than ISMIP-HOM, and were designed for use with models employing the shallow ice approximation (Payne et al., 2000).While ISMIP-HOM generally tests the accuracy of the momentum balance scheme over varying length scales, the EISMINT experiments are more focused towards assessing the time-dependent mass and energy balances.EISMINT II stresses the dynamic evolution of the system and provides a test for the long-term stability of our time-stepping scheme.We know of no other higher order model operating on an unstructured grid that has demonstrated a capacity to run forward on such timescales.
In EISMINT II A, a radially symmetric surface mass balance and temperature field are imposed on an initially icefree, flat bed.The model geometry and temperature field are allowed to evolve for 200 ka.Our chosen mesh had a cell size of 25 km in the horizontal dimension, and 10 vertical layers, for a total of approximately 32 000 degrees of freedom.We used a time step of 10 a, and on a 16-core machine with 400 GB of memory and 2.9 GHz per processor, this run (and EISMINT II F) took approximately 65 h.At the end of this period, the total energy and mass were changing by less than 10 4 % a −1 , implying near-steady-state conditions.the ice sheet were 3647 m and 255.4K respectively.These are important metrics for intercomparison performance, and lie near the benchmark means of 3688 m and 255.605K.
EISMINT II F is identical to experiment A, except that imposed surface temperatures are 15 • C colder throughout the model domain.The results of this experiment are similar to those of experiment A, albeit with a thicker ice sheet, and a slightly different basal temperature profile.The resulting thickness and basal temperature fields are given in Fig. 10 An interesting difference between the resulting temperature field here and that documented in the EISMINT II paper is that VarGlaS does not predict the breakdown in radial symmetry that occurred in all of the finite difference models of the intercomparison.We suspect that this is due to VarGlaS using an unstructured grid, which alleviates some of symptoms of grid dependency seen in the original experiment.This hypothesis is supported by the work of Saito et al. (2006), who provides evidence that suggests that the "spoking" phenomenon seen in Payne et al. (2000) is numerical in origin.It could be the case that using an unstructured grid removes the grid-induced pathways that allow these spokes to form.

Greenland
We applied VarGlaS to a large-scale problem in glaciology, namely the transient simulation of the Greenland ice sheet.The strategy in so doing was to initialize the model using measured present-day geometry, apply data assimilation tools to obtain an initial estimate of the basal traction field, and then allow Greenland's geometry, velocity, and temperature to evolve over 500 yr.We performed all simulations of Greenland using the first-order approximation for the momentum balance (see Eq. 6).

Data
We relied on the SeaRISE model set-up for input data (SeaRISE, 2012).Bedrock and surface geometry were from Bamber et al. (2001) 2010) for a surface velocity target.The Joughin data set is incomplete.The gaps were filled with balance velocities, with gradients between the two reduced by systematically exploring the uncertainties in the accumulation rate.Specifically, we used a steepest descent algorithm to minimize the misfit between the edges of the InSAR velocities and the corresponding balance velocities by varying the surface mass balance subject to the constraint that it remains within its reported point-wise error bounds.

A mesh for Greenland
The boundary of Greenland was digitized using the 1 m contour of the Bamber et al. (2001) thickness data.We created an initial two-dimensional (map plane) mesh by imposing a 2 km element size at the margins, grading to a variable but much coarser resolution at the centre of the ice sheet.This ensured that the mesh captured the complexity of the boundary while maintaining appropriate coarseness in the interior.We refined the mesh using the techniques of Sect.2.2.2.We extruded this footprint over 10 vertical layers.Greenland, when more highly resolved in the vertical dimension, demonstrates convergence problems during Newton's method solution process.Similar issues have been reported by other highresolution, higher order models, such as Elmer/Ice (Seddik et al., 2012) and ISSM (Larour et al., 2012).Larour et al. (2012) in particular demonstrated that this convergence issue is a result of very low aspect ratio elements producing poor matrix conditioning.This is a significant limitation, and attempts to overcome it are ongoing.Preliminary experimentation suggests that using the shallow ice approximation as a physically motivated preconditioner may help ameliorate this problem, but we have thus far been unsuccessful in applying this scheme at a continental scale.The generated mesh which is used for the data assimilation experiment has around 4.4 × 10 5 nodes, corresponding to 8.8 × 10 5 degrees of freedom for the first-order model.The horizontal resolution for  this mesh is variable, but grades from 500 m at the edges of the domain to approximately 100 km over the interior of the ice sheet.We used a coarser mesh with 4 × 10 4 nodes for our transient run, corresponding to 8 × 10 4 degrees of freedom for the first-order model.This corresponds to a minimum horizontal resolution of 5000 m up to approximately 100 km in the interior.

Data assimilation
We calculated a basal traction field using the techniques of Sect.2.2.3.We begin by calculating steady-state velocity and enthalpy fields for an arbitrary basal traction (with an initial guess of 4 Pa a m−2), with r equal to unity (which implies that basal traction is linearly scaled by thickness; this effectively eliminates the dependence of sliding speed on normal force, and eliminates the covariance between β 2 and h).With this initial state in hand, we ran the BFGS algorithm using a fixed viscosity (and an incomplete adjoint) for the first 10 iterations, before switching to a full adjoint.The selection of the regularization parameter α is motivated by the results of Balise and Raymond (1985), which indicate that below the one ice thickness length scale, variations in basal traction do not propagate to the surface.As such, we select α = h 2 , which corresponds to applying smoothing with a Gaussian kernel with standard deviation h to the basal traction field.
The temperature field was also recomputed every 50 evaluations of the objective function in order to maintain thermal equilibrium.After the first 10 iterations, the velocity field was visually indistinguishable from that of the data, and the convergence between temperature and velocity fields became a fixed point iteration on the enthalpy field.The BFGS algorithm was allowed to run for 200 evaluations of the objective function, which took approximately 3.5 h on a 16-core machine with 2.9 GHz per processor and 400 GB of memory.
Convergence of the algorithm is shown in Fig. 11.The observed and modelled velocities, along with basal traction and temperature fields, are shown in Fig. 12.To illustrate some of the fine-scale detail of both the mesh and the data assimilation result, Fig. 13 shows a close-up of Kangerdlugssuaq  glacier in eastern Greenland.The modelled velocity field matches the observed field closely, and the RMS error for the whole ice sheet is approximately 36 m a−1.For outlet glaciers like Kangerdlugssuaq, we see that surface velocity can be explained by spatially variable basal traction field composed of both low traction streaming features, and sticky pinning points that slow flow.Basal temperature is also related to basal traction, where fast sliding is associated with a melted bed.

Prognostic run
After performing the data assimilation procedure outlined above, we allowed the ice sheet to evolve through time for 500 yr under the C1 (constant climate) SeaRISE experiment, with a time step of 0.1 a.On a 16-core machine with 400 GB of memory and 2.9 GHz per processor, this run took 27 h.The present-day measured surface elevation of Greenland is not in a steady state with respect to the basal topography and the inverted velocity field, which is unsurprising given the errors associated with the basal topography, the surface velocity field, and the approximations made by the model physics.Because of this mismatch, large transient signals propagate through the system at the beginning of the run.We monitored the size of these transients as the L ∞ norm of the ∂ t S field, given by Fig. 14.Note the exponential decay rate; this gives some estimate of how much relaxation a model requires in order to eliminate transients.Results from Pritchard et al. (2009) show that the average ∂ t S of the Greenland ice sheet (GrIS) is −0.84 ma.An ice sheet model should have relaxed to around this level before any conclusions should be drawn from additional forcing being applied to it.For VarGlaS, with the specified data inputs, this process takes about 100 a.We also monitored the mass of the entire ice sheet throughout time.After the initial transient period of ice increase, we found an annual average total ice decrease of 10 −3 % a −1 (2.63 × 10 3 Gt a−1), which is in orderof-magnitude agreement with the GRACE-derived ice loss of approximately 6 × 10 −3 % a −1 (17.7 × 10 3 Gt a−1) (Baur et al., 2009).As seen in Fig. 15, the qualitative distribution of the basal temperature field changes relatively little over 500 yr.Although relative to the magnitude of their initial values, thickness and velocity also change relatively little; they show an interesting qualitative pattern that explains the mass loss of Fig. 14.We see that over time the combination of surface mass balance and basal traction tends towards thinning the interior of the ice sheet, while thickening the edges.This change is most notable where ice stream development is apparent, such as Jakobshavn and Kangerdlugssuaq, as well as the outlet glaciers of the southern GrIS.This increase in thickness generates a steeper ice surface at outlets, elevating the flux through the lateral boundaries.Since ice loss due to surface mass balance is constant (because surface mass balance itself is constant), the change in total mass for the GrIS is driven wholly by this increased boundary flux.This pattern of thickness and velocity evolution is also seen in comparable simulations of the GrIS, such as in Seddik et al. (2012) and Larour et al. (2012).This pattern is problematic, because it suggests an opposite redistribution of mass from that measured through remote sensing (Pritchard et al., 2009;Baur et al., 2009).More work is needed to determine if this mass loss rate is believable.

Discussion
VarGlaS performs well in a number of standard benchmarking experiments and can also simulate the evolution of the GrIS using higher order physics and thermomechanical coupling, as well as an advanced treatment of time evolution.
For the diagnostic and isothermal ISMIP-HOM benchmark experiments, both VarGlaS first-order and Stokes' solvers perform well with respect to existing benchmark results, with our first-order solver generally predicting values close the the Stokes' mean, and our Stokes' solver generally predicting velocities slightly slower than those reported by the benchmark.Our data assimilation procedure is able to reproduce effectively these simple imposed basal traction fields through model inversion.Both first-order and Stokes' solvers performed well on the prognostic ISMIP-HOM F, yielding velocity and surface elevation fields that are in good agreement with other models.Results for the EISMINT II experiments are similar, with model results comparing favourably to mean values from the original publication.Also, thickness and temperature fields compare well with the results of Pattyn (2003) and Saito et al. (2003), who also applied higher order models to these experiments.Note that this is the first time that the EISMINT II experiments have been performed with a higher order finite element model using an unstructured mesh.
VarGlaS' behaviour over the entire Greenland ice sheet echoes what has been determined by various investigators in the past, which is that after the relaxation of a strong transient signal derived from the incompatibility of flawed basal topography, surface velocities, and surface mass balance data, the ice sheet seems to be losing mass on the order of 10 −3 % a −1 , which is in agreement with contemporary estimates, although the qualitative distribution of this mass loss is not entirely believable.Performing simulations of this temporal length and at this spatial resolution has only been made possible by a combination of advances in ice sheet modelling technology, namely variable spatial resolution, data assimilation, and parallelism.
We used an adjoint model derived from automatic differentiation to invert the model at continental scale, yielding plausible velocity fields.Inversion was a critical step for a few reasons; insofar as the measured velocity and surface elevation fields are correct and the geometry of the ice sheet is self-consistent, the inversion procedure minimizes transient signals, and allows the model to reach a self-consistent state more quickly than would be otherwise possible with a less sophisticated starting procedure.This reasoning extends to  the calculation of a starting enthalpy field, which would be of a much lower quality without incorporating the very significant heat source due to friction at the bed.
We employed an anisotropically refined, variable resolution mesh in order to minimize superfluous degrees of freedom in slowly varying regions of the ice sheet while maintaining detailed solutions in regions of large velocity gradients.Although variable resolution modelling is not impossible with finite differences (Colella et al., 2000;Cornford et al., 2013), it is applicable to finite elements in a straightforward way.
Parallelism was another critical component in modelling the whole Greenland ice sheet.The number of degrees of freedom is simply too large for one processor to handle in a reasonable amount of time.We found that we retained better than 50 % parallel efficiency for nearly 1 million degrees of freedom and 16 processors, using an iterative solver.This corresponds to a speed-up of around a factor of 10.With the increasing availability of large computers with many processors, the benefit of incorporating parallelism into model de-sign is clear.Diagnostic modelling of very large ice sheets at resolutions similar to those of contemporary data products is possible with even modest computers.For example, in this paper, we perform data assimilation over the entire Greenland ice sheet at the horizontal scale of a few ice thicknesses using 16 cores in 3.5 h, which we believe is reasonable.For prognostic simulations with this same resolution, particularly past the century timescale, larger clusters are necessary.For example, for prognostic simulations with a 0.1 a time step and the same spatial resolution as the data assimilation exercise, using the same machine, we can simulate one year of model time in one day of real time.This is not practical.Fortunately, there now exist many clusters with many thousands of processors, and this type of computing power is becoming more available, such that, with sufficient parallel efficiency, solutions of problems at long timescales and at high spatial and temporal resolution are plausible.
VarGlaS currently does not possess a detailed treatment of marine terminal processes, namely grounding line migration dynamics and a prognostic calving law.Each of these Objective Function Value Fig. 11.rate of L BFGS B algorithm for surface velocity assimilation over the GrIS, using basal traction as a control variable.Vertical lines are where the enthalpy equations were recalculated presents its own computational and theoretical challenges.It is well known that grounding line dynamics operates at a scale that is typically sub-grid relative to the whole ice sheet field equations (Nowicki, 2007;Favier et al., 2012).In many cases, relatively small changes in the position of the grounding line probably do not greatly affect the evolution of an ice sheet at a continental scale, particularly in cases like Greenland where outlets are spatially localized, and topography is the primary control on outlet configuration.Nevertheless, scenarios such as the potential collapse of the West Antarctic ice sheet due to a fundamentally unstable bed geometry provide a strong motivation for getting the physics right, and higher order physics is mandatory in so doing.Var-GlaS is in a good position to incorporate a detailed treatment of grounding line dynamics.The detailed spatial resolution necessary for capturing the physics can be managed by the existing mesh refinement code.Additionally, VarGlaS also possesses robust free surface stabilization that will certainly be necessary for performing simulations of grounding line migration on complex real world topography.A prognostic calving would also be necessary to simulate accurately the grounding line dynamics.VarGlaS currently does not have the capacity to move the lateral bounds of its computational domain.Leaving the thickness at the lateral boundaries fixed is equivalent to making the assumption that any flux through the boundary in excess of the present rate calves immediately, or is otherwise removed from the ice sheet, and that the height of that boundary is fixed.This may be valid for continental scale model runs where the geometry of the ice sheet is not expected to adjust dramatically.This treatment is certainly insufficient for regional-scale experiments such as modelling the response of inland glaciers to ice shelf collapse.
Data assimilation is an essential part of correctly modelling ice surface velocities.Without using inverse methods to estimate the value of the basal traction field, the observed pattern of surface velocities is not well reproduced, and the present-day ice configuration is not (and should not be) tured.Additionally, without relying on inverse methods, long and computationally expensive spin-up procedures are required, which are not feasible from a processing standpoint, even with the efficient structure of VarGlaS and other modern ice sheet models.Simultaneously, we must recognize the limitations of inverting for the basal traction field.The inversion is ill-posed and often lacks a physically motivated stopping point for optimization algorithms.Even with the inclusion of a regularization term, there generally exist multiple solutions for the basal traction field that produce plausible surface velocity results (although many qualitative features must exist in all solutions).This lack of uniqueness makes drawing conclusions about basal conditions at specific points from inverted models tenuous.Also, the inversion procedure does not allow for time dependency of the basal traction field.Basal traction is believed to be fundamentally linked With the major increases in efficiency gained from parallelism and anisotropic mesh refinement, it is tempting to model on increasingly detailed meshes, simply due the philosophy that this will give us more detailed and thus more meaningful results.In an ideal scenario, the data from which we draw our surface elevations, bed elevations, surface velocities, and surface mass balances are effectively error-free and at spatial resolutions much finer than the grids on which we model.This may have been the case in the past, when the errors derived from using the shallow ice approximation were larger than those inherent in data, and low spatial resolution grids could average from a number of data points falling within a given pixel.It is certainly not the case now.It is simple to generate meshes with a horizontal resolution of near an ice thickness, but it is not simple to determine how to interpolate accurately from 5 km thickness data down to this resolution.Nevertheless, in some heavily studied re- gions of the GrIS, data at thickness-scale horizontal resolution exist (e.g.Jakobshavn, Russell-Isunnguata Sermia), and these should be incorporated when possible.For the future, we intend to incorporate error metrics over all of VarGlaS's input data into our mesh generation procedure, in order to avoid obtaining spurious results and using unnecessary computational resources as a result of over-resolving in regions where the data do not have a commensurate level of detail.

Conclusions
In this paper, we have introduced a novel next-generation ice sheet model called VarGlaS.VarGlaS is built upon the finite element package FEniCS, and borrows heavily on the innate capabilities of FEniCS, with such features as automatic differentiation of the adjoint state, an interface to a variety of www.the-cryosphere.net/7/1161/2013/The Cryosphere, 7, 1161-1184, 2013 efficient solvers, and demonstrably scalable parallelism.Var-GlaS eschews the common stress balance formulation of ice flow physics in favour of formulation as the problem of minimizing a scalar variational principle representing the conversion of gravitational potential energy into heat under the constraint of incompressibility.We use an enthalpy formulation for the energy balance equations, exchanging the temperature equation's contact problem for additional nonlinearity.VarGlaS treats conservation of mass using the kinematic boundary condition.We applied VarGlaS to the ISMIP-HOM benchmarks for higher order models, as well as several of the EISMINT II experiments.VarGlaS performed well in all of these benchmarks, proving that it correctly solves the field equations.Additionally, we applied VarGlaS's data assimilation to one of the benchmark experiments, showing that it can recover a simple basal traction field from an observed velocity field.
We then turned model to diagnostic and prognostic simulation of the Greenland ice sheet.We began by solving for basal traction using interferometrically derived surface velocities.Using this basal traction field, and other data from the SeaRISE data set, we solved for the geometry, temperature, and velocity of Greenland 500 yr into the future.We found that, after a brief period of relaxing transient signals, the model predicts a 10 −3 % a −1 decrease in total mass over the 500 yr period.

Fig. 1 .
Fig. 1.Observed surface velocity projected onto an anisotropically refined mesh of Greenland's northeast ice stream.

Fig. 2 .
Fig. 2.Parallel efficiency computed for one complete Newton solve of the first-order equations for meshes with different degrees of freedom.All linear solves were performed using parallel GMRES (generalized minimal residual method) preconditioned with Hypre-AMG.

Fig. 3 .
Fig. 3. Error between numerical and analytical solutions versus mesh cell size and time step.The slopes of these lines represent the convergence rates of the spatial and temporal discretization methods used by VarGlaS.
Figure  10shows the resulting thickness and basal temperature fields.Thickness and basal temperature at the centre point of
with updated basal topography in the Jakobshavn region from CReSIS, surface temperatures from www.the-cryosphere.net/7

Fig. 8 .
Fig.8.Convergence profile and modelled basal traction and velocity from inverting ISMIP-HOM C using the first-order approximation.RMS between observed and modelled velocities is 5 × 10 −3 m a −1 .

Fig. 10 .
Fig. 10.Thickness and basal temperature fields for EISMINT II A and F at 200 ka.

Fig. 12 .
Fig. 12. Modelled and measured surface velocities, basal traction, and basal temperature of the GIS after assimilation of surface velocity.RMS velocity mismatch is 36 m a−1.

Fig. 13 .
Fig. 13.Modelled and measured surface velocities, basal traction, and basal temperature of Kangerdlugssuaq glacier in eastern Greenland after assimilation of surface velocity.

)Fig. 14 .
Fig. 14.The L ∞ norm of the surface elevation change field and the total ice mass over a 500-year model run.

Fig. 15 .
Fig. 15.Velocity, surface elevation, and basal temperature change through a 500-year simulation of the Greenland Ice Sheet.

Table 1 .
Table of symbols.