Compatible Finite Elements for Glacier Modeling

Described in this article is the first application of two mixed finite-element methods to the equations of glacier evolution under different simplifying assumptions, along with a framework for the implicit solution of the coupled velocity-thickness equations. The first method uses Raviart–Thomas elements for velocity and piecewise constants for thickness and is a reframing of a classic staggered-grid finite-difference method to the case of unstructured triangular meshes. The second method uses Mardal–Tai–Winther elements for velocity and exhibits several desirable properties: second-order convergence of velocity and near-exact mass conservation while resolving both membrane and shear stresses in steep topography with thin ice.

Described in this article is the first application of two mixed finite-element methods to the equations of glacier evolution under different simplifying assumptions, along with a framework for the implicit solution of the coupled velocity-thickness equations.The first method uses Raviart-Thomas elements for velocity and piecewise constants for thickness and is a reframing of a classic staggered-grid finite-difference method to the case of unstructured triangular meshes.The second method uses Mardal-Tai-Winther elements for velocity and exhibits several desirable properties: second-order convergence of velocity and near-exact mass conservation while resolving both membrane and shear stresses in steep topography with thin ice.

G
laciers and ice sheets respond to climate change through interdependent surface melt and ice dynamics.Predicting such responses is critical to understanding future sea level rise, ecological perturbations, and hazards over multiple time scales.Due to the feedback between elevation and snowfall, glaciers often flow through steep places; indeed, glaciers in mountainous regions are currently losing mass at an outsized rate relative to their volume. 1Simultaneously, mass loss driven by ice dynamics is often most apparent at marine termini.Predicting glacier evolution requires models that can resolve the physics associated with both mountain and marine environments and that remain stable and mass conservative in the presence of extreme topography while still resolving membrane stresses.In this work, we explore the application of a class of mixed finite elements to glacier flow models that accomplish these goals.
Mixed finite elements posit that velocity uðx, z, tÞ and ice thickness Hðx, tÞ be discretized separately and that their discrete approximations live in different function spaces with differing degrees of regularity.The combinations of these spaces formulated such that differential operators map without approximation from one space to another are termed compatible.Geophysical models have long taken advantage of such schemes, which are free from spurious modes and have desirable conservation properties. 2 Indeed, staggered-grid finite-difference methods (e.g., the C-grid) can be interpreted as a specialization of compatible elements to quadrilaterals.
Ice-sheet modelers have also used compatible elements for shallow ice approximation (SIA), which simplifies the Stokes' equations by assuming that membrane-resistive and vertical-resistive stresses are negligible and that the pressure field is hydrostatic in the vertical direction, which is appropriate for slow flow and shallow aspect ratios.However, the models of intermediate complexity that do not neglect membrane stresses, such as the Blatter-Pattyn 3 approximation (BPA), (which assumes hydrostatic pressure in the vertical direction, negligible vertical-resistive stresses, and relatively small bed slopes), and shallow shelf approximations (SSAs) (like the BPA but also neglecting vertical shear stresses), are often more appropriate for marine termini and mountainous environments.Finite-element implementations of these approximations (when coupled to the evolution of ice thickness), to date, have used either specialized grids or equalorder elements with stabilization.Dual grids (e.g., in Hoffman et al. 4 ), when used in conjunction with the finite-volume method, are conservative and stable but can be geometrically restrictive or require additional computational bookkeeping.Equal-order finite elements with artificial viscosity or a streamline upwind Petrov-Galerkin formulation (e.g., in Dos Santos et al. 5 ) are simple, efficient, and work well for regions of smooth surface and bed topography but may yield oscillations in the presence of steep gradients and do not necessarily guarantee a positive thickness, even in the absence of melt.
In an effort to extend the advantages of dual-grid formulations to a standard triangular mesh, first, a compatible finite-element pair, namely, the lowestorder Raviart-Thomas element coupled with a zerothorder discontinuous Galerkin element, is described.We use full nonlinear coupling for time integration, which allows large time steps.Because this discretization cannot resolve the membrane stresses present in the BPA and the SSA, our scheme is modified to use the Mardal-Tai-Winther (MTW) element, which has the same favorable properties as the Raviart-Thomas space with additional smoothness.These elements were evaluated on well-understood test cases and against a manufactured solution to confirm expected rates of convergence.Finally, to test the model under the motivating cases of both extreme topography and a marine terminus with real topography, the model was applied to the simulation of conditions reflective of the Pleistocene (including termination in a glacial lake) for a mountain basin in western Montana, where it exhibits the desired properties of numerical stability and mass conservation.

EQUATIONS OF ICE FLOW
Ice-sheet models solve for the ice thickness Hðx, tÞ and velocity uðx, z, tÞ ¼ ½uðx, z, tÞ,vðx, z, tÞ T : The evolution of the thickness is governed by the principle of mass conservation, which is expressed as a conservation law where an overbar indicates a vertical average, _ a is the specific mass balance (i.e., the rate of ice-equivalent snowfall or melting) (see Table 1), n x is a unit-normal vector in the map plane, and r x ½@=@x 1 , @=@x 2 T : The velocity is computed by solving the BPA to the Stokes' equations, which balances viscous and gravitational driving stresses as where q i and g are ice density and gravitational acceleration, respectively; Sðx, tÞ ¼ Bðx, tÞ þ Hðx, tÞ is the surface elevation; and r x,z ½@=@x 1 , @=@x 2 , @=@z T : The elevation of the ice base is Bðx, tÞ ¼ max z sl À q i q w Hðx, tÞ, z topo ðxÞ with z sl representing the sea level and z topo ðxÞ representing the topography.As such, (3) applies to both grounded and floating ice._ 1 is the strain rate tensor subject to the simplifications of the BPA The viscosity, which depends inversely on the effective rate of strain and thus describes a shear-thinning fluid, is given by Glen's flow law with A representing ice hardness, n the flow exponent, and _ II the second invariant of the strain rate tensor.
At the surface and basal boundaries, we impose stress-free and sliding conditions with b representing a friction coefficient, N as the effective pressure, m and p as the exponents that reflect frictional nonlinearity, and where n x,z is the 3-D normal vector.At lateral domain boundaries, we impose a normal stress that corresponds to an impenetrable sidewall.
with a representing a (usually large) coefficient that governs the degree of normal "friction" imposed at the boundary.

COMPUTATIONAL MODELING OF ICE SHEETS AND GLACIERS
HðxÞ À1 : The spatial derivatives in this coordinate system for an arbitrary function fðx, zÞ are with subsequent differential operators inheriting these definitions.We note that this coordinate system is undefined if the thickness goes to zero; we conceptualize margins as "ice cliffs" of nonzero-but perhaps very small-thickness.

Weak Form
We discretize (1) and ( 3) by dividing the model domain X into single-layered prismatic cells X T ¼ X T Â ½0, 1 with lateral boundaries @X T : We obtain a weak form for this cell by multiplying with test functions with the appearance of thickness in (13) resulting from the terrain-following coordinate, and where we have used boundary conditions ( 7) and ( 8) to eliminate integrals at the bottom and top of the prism.

COMPATIBLE ELEMENTS FOR THE SIA
We first consider the SIA, which can be derived from (13) by neglecting horizontal derivatives of the velocity in (5).
To proceed, we must specify function spaces V 1 and V 2 :

Thickness Space
We restrict our choice of V 2 & L 2 ðXÞ to function spaces that are constant over an element (and thus discontinuous across element boundaries), namely, the zeroth-order discontinuous Galerkin space (DG0).This space is closely associated with the finite-volume method and inherits its mass conservation properties as well as its tendency to preserve the monotonicity of solutions when paired with a carefully constructed transport scheme.It also inherits the necessity of specifying a numerical flux to couple adjacent elements.Summing ( 12) over all triangles, we have ð where @X are the set of internal triangular mesh boundaries; x is the jump of the test function; and fHg is the average across a cell boundary, respectively; and d uH is the Lax-Friedrichs flux To maintain consistency between geometric variables, we discretize the bed topography z B and the ice base B such that they also live in V 2 ðXÞ:

Velocity Space
Next we select a function space for the velocity that is compatible with our chosen thickness space.The compatible spaces are constructed so that projection and differential operators commute, which means that the same result is obtained whether applying the operator before or after the approximation of the function with finite elements.This notion can be represented with the commutative diagram (16) where p 1 and p 2 project from continuous to discrete spaces.In this case, the relevant differential operator is the divergence, which appears in both ( 12) and (13).Because the SIA system requires that velocity (and its associated test function) have a well-defined divergence, we seek & Hðdiv, XÞ where Compatible elements are useful because they ensure that when we take an inner product between the divergence of a test function in V 1,x ðXÞ and a function in V 2 ðXÞ, as is the case in the second term of (13), that nonphysical high frequencies in the thickness field are not inadvertently canceled while computing a velocity solution, which would allow such spurious thickness modes to remain in the solution.Indeed, the so-called "checkerboard" instability has long been

COMPUTATIONAL MODELING OF ICE SHEETS AND GLACIERS
observed when attempting to solve Stokes' and shallow-water equations, which are structurally similar to our ice flow model.More formally, compatible elements naturally satisfy a discrete inf À sup condition. 6he common finite-element space CG1, when used for each velocity component and in conjunction with a DG0 thickness space, does not satisfy this condition, and our experimentation confirms that it leads to evident numerical artifacts.
The discrete function spaces that are subsets of HðdivÞ have vector-valued functions with normal components that are continuous across cell edges but are tangentially discontinuous.Several choices exist for the HðdivÞ elements that are compatible with the DG0 element.We use the lowest-order Raviart-Thomas element (RT1), which has one degree of freedom per edge representing the average normal velocity.When combined with a DG0 thickness, this generalizes the C-grid (which also has edge-normal velocity components) to unstructured triangular meshes.We note here that the flux F ¼ uH should also be in V 1,x ðXÞ, which is nontrivial due to the discontinuity in H: This condition can be satisfied either by taking the flux as an ancillary degree of freedom defined by an explicit L 2 projection, or by the specification of a numerical flux as in (15), which renders F Á n x continuous across element boundaries. 2 As described previously, we take the latter approach.
Although HðdivÞ regularity is sufficient for uðxÞ, the shear-stress terms in 1 require that uðx, 1Þ have a continuous first derivative in the vertical dimension (ignoring in this work the fact that the SIA can be solved directly).To this end, we define the tensor-product element with V 1,z & H 1 ð½0, 1Þ, where We specify V 1,z ð½0, 1Þ by directly enumerating its basis functions and require that it be spanned by zeroth-and fourth-order monomials.If we choose bases that are orthogonal with respect to the L 2 inner product, then the velocity vector can be expanded as which has the capacity to represent the exact SIA solution for a suitable quadrature rule.This form is particularly convenient in that the depth-averaged velocity appears directly as a degree of freedom, while the basis function corresponding to deformation integrates to zero with respect to 1: Again summing over all elements, the global weak form is where we retain the full-strain rate tensor, but recall that all terms in _ 1 involving horizontal derivatives are zero under the SIA.The gradients in (21) should be understood in an elementwise sense, which may lead to jumps across element boundaries for some finiteelement spaces (including one that we use later), but we neglect these unless explicitly written.We note that we apply integration by parts to the driving stress.Because S is discontinuous, we must include another numerical pseudoflux, and we take the average of values in cells on each side of the boundary.Unlike in many treatments of analogous shallow-water equations, we do not split the driving stress into terms involving thickness gradients from those involving basal gradients but rather treat these terms identically so that the scheme is well balanced in that it preserves "lake-at-rest" solutions.

Discretization in Time and Numerical Solutions
We discretize in time with a hÀ method, and as such must solve a nonlinear system of equations involving both H and u at each time step, for which we use a Picard iteration.The linearized system is ð COMPUTATIONAL MODELING OF ICE SHEETS AND GLACIERS where an asterisk denotes a quantity computed at the previous iteration.There are options with respect to which variable to lag when u, H, or S appear nonlinearly; it was found that solving for (versus holding fixed) diffusive terms leads to better stability.For example, in (22), we solve for the velocity, which depends on the thickness gradient.To avoid nonpositive thickness resulting from melting, at each Picard iterate we set any thickness H < 10 À12 to H ¼ 10 À12 : Although the proposed scheme avoids nonpositive thickness from ice flow, this projection does induce a mass conservation error in the sense that the volumetric change of ice over a finite time interval is not necessarily equal to the integrated surface's mass balance over the same interval, a property that holds for all time-integration schemes with such an inequality condition. 7he model was implemented using Firedrake. 8For integrals arising in the assembly of the finite-element matrix, we allow Firedrake to automatically determine a suitable quadrature rule.A direct method is used to solve the resulting linear system.There was no particular effort made at optimization, leaving an evaluation of computational efficiency for later work.

Verification With the Halfar Solution
We first evaluate the RT1-DG0 element on the analytical "Halfar solution," 9 which represents a radially symmetric ice mass that evolves with neither mass balance nor sliding.Beginning with an initial dome height of 3 km and a radius of 500 km, the model was allowed to evolve from t 0 ¼ 292 a to T ¼ 10t 0 with Dt ¼ t 0 =10: Figure 2(a) shows the L 2 norm of the error relative to the exact solution (due to the discontinuous nature of the DG0 space, it is difficult to produce a meaningful comparison to a continuous function, so we project the solution into a CG1 finite-element space before computing the norm).The resulting error is linearly proportional to cell size, which is consistent with the theory for DG0 (although it is worth noting that even higher-order schemes would exhibit similar convergence due to the inability of the Halfar solution to be represented accurately in a piecewise polynomial basis).Under the chosen parameters, errors in time are substantially smaller than those in space.Figure 2(b) shows the spatially resolved differences between the exact and true solution along a diagonal transect.The emergent pattern indicates a transport scheme that is somewhat too diffusive, which is expected with a Lax-Friedrichs flux.The volume change between the beginning and end of the simulation is proportional to machine precision times the number of time steps.This is also true if we use a time step Dt ¼ 9t 0 , such that our time integration occurs in just a single time step.Although such a solution is not particularly accurate, as shown in Figure 2(b), it remains nonoscillatory and mass conservative.

Verification of Mass Conservation in Steep Topography
To assess the performance of the method in steep topography, we also compare it against the Jarosch-Schoof--Anslow solution, 10 in which ice flows over a sheer cliff with surface mass balance defined such that the cliff is higher than the ice surface below it.Although such a circumstance violates the assumption of small bed slopes required in deriving the BPA and the SIA, it is nevertheless common in mountainous domains.This scenario is numerically challenging due to the tendency for models to produce spurious ice over discontinuities, leading to expansion of the ice sheet past the location at which an along-flow integration of the surface mass balance function would place in the margin.Although the exact solution for this problem is 1-D, the problem was simulated in two dimensions with a periodic boundary condition transverse to flow.Following 50 ka of time integration with Dt ¼ 20 a and a nominal cell size of Dx ¼ 125 m; we find a total modeled volume of 97.2% of the analytical steady-state volume, which is a similar error to those of previous works. 10We used h ¼ 1, which allowed us to far exceed a quite-restrictive Courant-Friedrichs-Lewy criterion.The modeled solution is similar to the exact one (Figure 3) but slightly thinner, again indicating an over-diffuse transport scheme.

COMPATIBLE ELEMENTS FOR THE BPA
We next turn to the development of a compatible finite element that also resolves membrane stresses, again restricting our search to a DG0 space for thickness.

COMPUTATIONAL MODELING OF ICE SHEETS AND GLACIERS
A tempting path would be to continue with RT1 elements as before.However, this approach will not be successful because membrane stresses require that the velocity vector u lives in H 1 , where which is to say that both components of the velocity vector must have well-defined first derivatives, a more stringent smoothness condition than that of HðdivÞ: This resulting compatibility requirement corresponds to the modified commutative diagram (25) which requires unusual finite-element spaces.The one element that was constructed to satisfy these requirements is the MTW element. 11This element spans a vector-valued third-degree polynomial space, with additional constraints of order-one polynomials along element boundaries and a constant divergence, which requires two degrees of freedom per edge representing normal moments plus one degree of freedom per edge representing the first tangential moment (Figure 1).With respect to the latter, the element is nonconforming in the sense that tangential velocity components are only continuous across element boundaries in the average rather than in a pointwise sense.The MTW element was originally developed for an incompressible Darcy-Stokes flow, which has substantial similarities to the BPA in that solutions live in H 1 in the Stokes flow limit and HðdivÞ in the Darcy limit.Furthermore, it was designed to be compatible with DG0 pressure fields.
The solution procedure for the MTW-DG0 scheme is identical to that which was described previously for the RT1 element, with two exceptions.First, boundary conditions are imposed by adding a weak form of ( 9) to (21).Second, we require explicit gradients of the ice thickness and the base to account for the extra derivative terms arising from the 1-coordinate transformation.We accomplish this efficiently through a projection, solving at each Picard iteration the linear equation for rx

COMPUTATIONAL MODELING OF ICE SHEETS AND GLACIERS
(and similarly for the gradient of H) which is positive definite and converges in just a few iterations of conjugate gradient.

Validity of the MTW Element in the SIA Limit
Because H 1 & HðdivÞ and because both elements are compatible with DG0, the MTW element should perform similarly to the RT1 element under the SIA's assumptions.To show that this is the case, we performed both numerical experiments described in the previous section with the RT1 element replaced by the MTW element.In both cases (Figures 2 and 3), we find the solutions to be nearly identical.We do not advocate for the MTW element to be used for SIA models as it has triple the degrees of freedom, but such tests indicate that it inherits the simpler element's stability and accuracy properties.

Intercomparison With ISMIP-HOM
The principal advantage of the MTW element relative to RT1 is its capacity to correctly resolve membrane stresses.We assess this by applying the scheme to the ISMIP-HOM intercomparison experiments, 12 which test models on simplified geometries at various length scales.Figure 4(a) and (b) shows the respective velocity solutions for ISMIP-HOM A, which describes nonslip flow over a bed that varies sinusoidally in both dimensions, and ISMIP-HOM C, which describes flow over a flat bed with sinusoidal basal traction, both evaluated at the L ¼ 40-km-length scale (the results are similar for other length scales that were omitted for brevity).
In both cases, the solutions produced by the MTW element agree well with the distribution of the benchmark ensemble.ISMIP-HOM F, which tests the model's prognostic solution capacity through periodic flow

COMPUTATIONAL MODELING OF ICE SHEETS AND GLACIERS
over an isolated Gaussian bump, was also performed.Figure 4(c) shows that the steady-state surface elevation agrees well with the benchmark ensemble.

Verification With the Method of Manufactured Solutions
We next assessed the convergence of the MTW-DG0 element pair under both spatial and temporal refinement using the method of manufactured solutions.We note that these tests were performed under the assumptions of the SSA in which we neglect the second, vertically varying basis function in V 2,z : The reason for this is that because the vertical basis is not complete, there is no notion of refinement under which convergence to an exact solution could be achieved (despite it being a good ansatz for vertical variation in ice velocity, as indicated by the results of the previous section).
To construct a manufactured solution, we assumed a bed elevation given by with L ¼ 10 4 m; a constant basal traction of b 2 ¼ 10 2 Pa a m À1 with p ¼ 0 and m ¼ 1, and an analytical ice surface with T ¼ 25 a and C ¼ 200 m: The manufactured thickness solution is then H mms ¼ S mms À B: The velocity solutions are specified so that the resulting artificial forcing functions are relatively small, with and k mms ¼ 10 1 m a À1 .Such solutions are consistent with the periodic boundary conditions that we apply.Firedrake's symbolic manipulation capabilities were used to automatically derive the residual forcings associated with the specified solutions.Figure 4(a) shows the error-under-grid refinement of the velocity solution at t ¼ 0: We find that the element combination achieves expected second-order convergence.Figure 5(b) shows the CG1-projected thickness error at t ¼ T under varying choices of grid and time-step refinement for h ¼ 1: We find that as long as the error is not overwhelmed by the spatial discretization, then our method achieves expected linear convergence.Using h ¼ 1=2, we achieve convergence proportional to Dt 2 but observe less robust convergence of the nonlinear solver when applied to domains with realistic topography.

APPLICATION TO REALISTIC TOPOGRAPHY
As a realistic test of the MTW-DG0 element and our solver, a two-stage simulation, with topography from Bear and Sweathouse Creeks in the Bitterroot

COMPUTATIONAL MODELING OF ICE SHEETS AND GLACIERS
Mountains of western Montana, USA, was performed.Although both drainages are currently unglacierized, Pleistocene glaciation incised deep canyons with extremely steep walls and hanging valleys into the metamorphosed granite of the Idaho Batholith, creating challenging conditions for stability and mass conservation.Glacial Lake Missoula also filled the broad canyon mouths, which we include to test the model's capacity for simulating floating ice (although we neglect calving).
A mesh was generated by automatically delineating watershed boundaries in the high-elevation mountains and manually delineating arbitrary domain edges in the lower valley (see Figure 6).The nominal spatial resolution is 150 m, for a total of 16,000 cells.A maximum time step of Dt ¼ 10 a was used, but this was occasionally reduced automatically (and then increased again) when the Picard iteration failed to converge within 10 iterations.For reference, each time step took roughly 4 s running on four Intel i7-11800H 2.3-GHz cores.The maximum and minimum elevations are approximately

COMPUTATIONAL MODELING OF ICE SHEETS AND GLACIERS
2300 m and 1100 m, respectively, with a lake level of 1500 m.We set b 2 ¼ 100 m À1 a, p ¼ 1, and set the water pressure to the greater of three quarters of overburden or the depth below lake level.Importantly, zero-flux boundary conditions were imposed along all the external edges, and thus, no mass entered or left the system except through surface mass balance.
In the first stage of our test, glaciers were built by imposing the elevation-dependent specific mass balance _ a ¼ 5 a À1 ðS À 2100 mÞ (30) which yielded accumulation on only the highest peaks, as shown in Figure 6(a) and (b).The system evolved under this mass balance function for 1000 a, reaching an approximate steady state [Figure 6(c) and (d)].The resulting maximum thickness is 750 m with a maximum velocity of 350 ma À1 .We note the development of substantial floating ice.We next set the specific mass balance to _ a ¼ 0 and allowed the system to evolve for an additional 9000 a, as shown in Figure 6(e) and (f).Because of the no-flux boundary condition, ice in the (former) accumulation area flowed into the valley and stagnated in a large, flat pool.The remaining ice in steep areas became very thin (often millimeters or less).Although not a scenario that mirrors anything of particular physical relevance, this allows us to evaluate the model's capacity to conserve mass under such conditions.Figure 7 shows a time series of total mass for both stages.In this second stage, the mass changed by approximately À1 m 3 at the end of the simulation, which is consistent with an accumulated round-off error.As such, the element described here has satisfactory massconservation properties and can simulate glaciation and deglaciation without numerical issues.

CONCLUSION
In this work, the first application of two mixed finite-element methods to the equations of glacier motion, along with a fully coupled implicit timeintegration method was described.The resulting models are stable, exhibit the expected errors, are locally mass conservative and positivity preserving, and are free from spurious thickness modes.They are also relatively simple in the sense that they operate only on the cells and edges of a standard triangular grid.As such, this model is a practical choice for modeling many glaciological scenarios.Both elements are well-suited for simulations with steep topography, and the MTW-DG0 element can also accurately characterize membrane stresses.
The model presented here has a few important limitations that also suggest some future goals.First, because we use a DG0 element for thickness, spatial accuracy of the method is linear in mesh resolution.This could be improved with higher-order discontinuous Galerkin elements for thickness but would require the development of new compatible elements; the implementation of such new elements is made easier with modern finite-element libraries such as Firedrake.Although we did not explore the accuracy of grounding line and calving front resolution, higher-order elements could be particularly useful for their resolution through p refinement.Finally, much work is still needed to develop linear and nonlinear solvers optimized for these finite-element pairs.

FIGURE 1 .
FIGURE 1. Degrees of freedom associated with (a) the mixed lowest-order Raviart-Thomas and zeroth-order discontinuous Galerkin element and (b) the mixed MTW and zeroth-order discontinuous Galerkin element.The solid arrows represent the moments of the velocity vector evaluated at the corresponding edge that lead to continuity between elements, while the dashed arrows represent the moments that lead to a nonconforming approximation.The black dot represents elementwise constant geometric variables.

FIGURE 2 .
FIGURE 2. (a) Error in computed thickness at t ¼ 10t 0 for the radially symmetric, zero-accumulation Halfar solution as a function of cell size for both the RT1-DG0 and MTW-DG0 elements.As expected, due to the DG0 discretization of thickness, both finite-element pairs exhibit linear convergence in Dx: (b) Computed solution for the RT1-DG0 element as well as the spatially resolved difference between predicted and modeled thickness.

FIGURE 3 .
FIGURE 3. Exact and modeled steady-state surface elevation for flow over a bedrock step assuming the SIA using both the RT1-DG0 and MTW-DG0 elements (which are indistinguishable) as well as the difference between the exact and modeled solutions.Also shown is the same test with membrane stresses included using the MTW-DG0 element, which deviates slightly from the SIA solution immediately in the vicinity of the step.

FIGURE 4 .
FIGURE 4. Velocity at y ¼ L=4 for L ¼ 40 km for the ISMIP-HOM intercomparisons experiments A (a) and C (b), and (c) the steady-state surface-elevation anomaly for the 'slip' variant of experiment F (note that we have retained the discontinuous character of the solution, accounting for the stair step pattern.Shaded region indicates two standard deviations in the ensemble of higher-order solutions submitted for the original ISMIP-HOM.

FIGURE 5 .
FIGURE 5. (a) Error in computed velocity at t ¼ 0 against a manufactured solution as a function of cell size, which exhibits second-order convergence.(b) Error in computed thickness at t ¼ T against a manufactured solution as a function of time step for multiple cell sizes using the shelfy-stream approximation.Note the expected linear convergence, except when spatial discretization error creates a plateau.The dotted line shows second-order convergence in time step for the finest spatial resolution with h ¼ 1=2, i.e., Crank-Nicholson.

FIGURE 6 .
FIGURE 6. Evolution of the thickness and and velocity fields for our realistic steep-topography experiment during the glacier's growth (t ¼ 200 a; a and b), at a steady glacial maximum (t ¼ 1000 a; c and d) and following an additional 1000 years of evolution in the absence of any accumulation or ablation (e and f).

FIGURE 7 .
FIGURE 7. Ice volume through time for our realistic steeptopography experiment.During the first 1000 years of the simulation, we build a glacier using the mass balance function described in the article.At t ¼ 1000 a; the surface mass balance is set to uniformly zero, which, in tandem with a no-flux boundary condition, should yield exact maintenance of volume through time (note the change in scale for both axes).At t ¼ 10, 000 a; we see a volume change of roughly À1 m 3 , which is consistent with accumulating a numerical roundoff error.

TABLE 1 .
Selected symbols used in this article.