MPAS-Albany Land Ice (MALI): A variable resolution ice sheet model for Earth system modeling using Voronoi grids

. We introduce MPAS-Albany Land Ice (MALI), a new, variable resolution land ice model that uses unstructured Voronoi grids on a plane or sphere. MALI is built using the Model for Prediction Across Scales (MPAS) framework for developing variable resolution Earth System Model components and the Albany multi-physics code base for solution of coupled systems of partial-differential equations, which itself makes use of Trilinos solver libraries. MALI includes a three-dimensional, ﬁrst-order momentum balance solver (“Blatter-Pattyn") by linking to the Albany-LI ice sheet velocity solver, as well as an 5 explicit shallow ice velocity solver. Evolution of ice geometry and tracers is handled through an explicit ﬁrst-order horizontal advection scheme with vertical remapping. Evolution of ice temperature is treated using operator splitting of vertical diffusion and horizontal advection and can be conﬁgured to use either a temperature or enthalpy formulation. MALI includes a mass-conserving subglacial hydrology model that supports distributed and/or channelized drainage and can optionally be coupled to ice dynamics. Options for calving include “eigencalving”, which assumes calving rate is proportional to extensional strain rates. 10 MALI is evaluated against commonly used exact solutions and community benchmark experiments and shows the expected accuracy. We report ﬁrst results for the MISMIP3d benchmark experiments for a Blatter-Pattyn type model and show that results fall in-between those of models using Stokes ﬂow and L1L2 approximations. We use the model to simulate a semi-realistic Antarctic Ice Sheet problem for 1100 years at 20 km resolution. MALI is the glacier component of the Energy Exascale Earth System Model (E3SM) version 1, and we describe current and planned coupling to other components.

While these efforts represent significant steps forward, next-generation ISMs continue to confront new challenges. These come about as a result of (1) applying ISMs to larger (whole-ice sheet), higher-resolution (regionally O(1 km) or less), and more realistic problems, (2) adding new or improved sub-models of critical physical processes to ISMs, and (3) applying ISMs as partially or fully coupled components of ESMs. The first two challenges relate to maintaining adequate performance and 15 robustness, as increased resolution and/or complexity have the potential to increase forward model cost and/or degrade solver reliability. The latter challenge relates to the added complexity and cost associated with optimization workflows, which are necessary for obtaining model initial conditions that are realistic and compatible with forcing from ESMs. These challenges argue for ISM development that specifically targets the following model features and capabilities: 1. parallel, scalable, and robust, linear and nonlinear solvers 20 2. variable and / or adaptive mesh resolution 3. computational kernels based on flexible programming models, to allow for implementation on a range of High-Performance Computing (HPC) architectures 1 4. adjoint capabilities for use in high-dimensional parameter field optimization and uncertainty quantification Based on these considerations, we have developed a new land ice model, the MALI model, which is composed of three 25 major components: 1) model framework, 2) dynamical cores for solving equations of conservation of momentum, mass, and energy, and 3) modules for additional model physics. The model leverages existing and mature frameworks and libraries, namely the Model for Prediction Across Scales (MPAS) framework and the Albany and Trilinos solver libraries. These have allowed us to take into consideration and address, from the start, many of the challenges discussed above. We discuss each of these components in more detail in the following sections.
The MPAS Framework provides the foundation for a generalized geophysical fluid dynamics model on unstructured spherical and planar meshes. On top of the framework, implementations specific to the modeling of a particular physical system (e.g., land ice, ocean) are created as MPAS cores. To date, MPAS cores for atmosphere (Skamarock et al., 2012), ocean Petersen et al., 2015Petersen et al., , 2018, shallow water (Ringler et al., 2011), sea ice , and land ice have 5 been implemented. At the moment the land ice model is limited to planar meshes due to the planar formulation of the flow models; however, we have an experimental implementation of the flow model for spherical coordinates that enables runs on spherical meshes. The MPAS design philosophy is to leverage the efforts of developers from the various MPAS cores to provide common framework functionality with minimal effort, allowing MPAS core developers to focus on development of the physics and features relevant to their application. 10 The framework code includes shared modules for fundamental model operation. Significant capabilities include: -Description of model data types. MPAS uses a handful of fundamental Fortran derived types for basic model functionality. Core-specific model variables are handled through custom groupings of model fields called pools, for which custom accessor routines exist. Core-specific variables are easily defined in XML syntax in a Registry, and the framework parses the Registry, defines variables, and allocates memory as needed. 15 -Description of the mesh specification. MPAS requires 36 fields to fully describe the mesh used in a simulation. These include the position, area, orientation, and connectivity of all cells, edges, and vertices in the mesh. The mesh specification can flexibly describe both spherical and planar meshes. More details are provided in the next section.
-Distributed memory parallelization and domain decomposition. The MPAS Framework provides needed routines for exchanging information between processors in a parallel environment using Message Passing Interface (MPI). This 20 includes halo updates, global reductions, and global broadcasts. MPAS also supports decomposing multiple domain blocks on each processor to, for example, optimize model performance by minimizing transfer of data from disk to memory. Shared memory parallelization through OpenMP is also supported, but the implementation is left up to each core.
-Parallel input and output capabilities. MPAS performs parallel input and output of data from and to disk through the 25 commonly used libraries of NetCDF, Parallel NetCDF (pnetcdf), and Parallel Input/Output (PIO) (Dennis et al., 2012).
The Registry definitions control which fields can be input and/or output, and a framework streams functionality pro- This allows explicit definition of model epochs in terms of years, months, days, hours, minutes, seconds, and fractional seconds and can be set to three different calendar types: Gregorian, Gregorian no leap, and 360 day. This flexibility helps enable multi-scale physics and simplifies coupling to ESMs. To manage the complex date/time types that ensue, MPAS framework provides routines for arithmetic of time intervals and the definition of alarm objects for handling events (e.g., when to write output, when the simulation should end).

5
-Run-time configurable control of model options. Model options are configured through namelist files that use standard Fortran namelist file format, and input/output are configured through streams files that use XML format. Both are completely adjustable at run time.
-Online, run-time analysis framework. See section Section 6.2 for examples.
Additionally, a number of shared operators exist to perform common operations on model data. These include geometric 10 operations (e.g., length, area, and angle operations on the sphere or the plane), interpolation (linear, barycentric, Wachspress, radial basis functions, spline), vector and tensor operations (e.g., cross products, divergence), and vector reconstruction (e.g., interpolating from cell edges to cell centers). Most operators work on both spherical and planar meshes.

Model Meshes
The MPAS mesh specification is general enough to describe unstructured meshes on most two-dimensional manifold spaces, 15 however most applications use centroidal Voronoi tesselations (Du and Gunzburger, 2002) on a sphere or plane. This manuscript focuses on applications with planar centroidal Voronoi meshes, with some additional consideration of spherical centroidal Voronoi meshes. Voronoi meshes are constructed by specifying a set of generating points (cell centers), and then partitioning the domain into cells that contain all points closer to each generating point than any other. Edges of Voronoi cells are equidistant between neighboring cell centers and perpendicular to the line connecting those cell centers. A planar Voronoi tesselation 20 is the dual graph of a Delaunay triangulation, which is a triangulation of points in which the circumcircle of every triangle contains no points in the point set. Voronoi meshes that are centroidal (the Voronoi generator is also the center of mass of the cell) have favorable properties for some geophysical fluid dynamic applications (Ringler et al., 2010) and maintain high quality cells because cells tend towards equidimensional aspect ratios, and mesh resolution (where non-uniform) changes smoothly.
On both planes and spheres, Voronoi tesselations tend toward perfect hexagons as resolution is increased. Note that, while the 25 MPAS mesh specification supports quadrilateral grids, such as traditional rectangular grids, they are described as unstructured, which introduces significant overhead in memory and calculation over regular rectangular grid approaches.
Because MPAS meshes are two-dimensional manifold spaces, they are convenient for describing geophysical locations, either on planar projections or directly on a sphere. Because they are unstructured, meshes can contain varying mesh resolution and can be culled to only retain regions of interest. Planar meshes can easily be made periodic by taking advantage of the 30 unstructured mesh specification. MPAS meshes are static in time. The vertical coordinate, if needed by a core, is extruded from the base horizontal mesh. Each core chooses its own vertical coordinate system. A comprehensive suite of tools for the generation of centroidal Voronoi tesselations on a plane or sphere has been developed, as well as tools for modifying existing a b x y z Figure 1. MALI grids. a) Horizontal grid with cell center (blue circles), edge midpoint (red triangles), and vertices (orange squares) identified for the center cell. Scalar fields (H, T ) are located at cell centers. Advective velocities (un) and fluxes are located at cell edges. b) Vertical grid with layer midpoints (blue circles) and layer interfaces (red triangles) identified. Scalar fields (H, T ) are located at layer midpoints.
The basic unit of the MPAS mesh specification is the cell. A cell has area and is formed by three or more sides, which are referred to as edges. The endpoints of edges are defined by vertices. Figure 1 illustrates the relationships between these mesh primitives. The MPAS mesh specification utilizes 36 fields that describe the position, orientation, area, and connectivity of the 5 various primitives. Only four of these fields (x, y, z cell positions and connectivity between cells) are necessary to describe any mesh, but the larger set of fields in the mesh specification provide information that is commonly used for routine operations.
This avoids the need for the model to calculate these fields internally, speeding up the process of model initialization and integration.
MALI typically uses centroidal Voronoi meshes on a plane. Spherical Voronoi meshes can also be used, but little work has 10 been done with such meshes to date. MALI employs a C-grid discretization (Arakawa and Lamb, 1977) for advection, meaning state variables (ice thickness and tracer values) are located at Voronoi cell centers, and flow variables (transport velocity, u n ) are located at cell edge midpoints (Figure 1). MALI uses a sigma vertical coordinate (specified number of layers, each with a spatially uniform layer thickness fraction, see (Petersen et al., 2015) for more information): where s is surface elevation, H is ice thickness, and z is the vertical coordinate.
A set of tools supporting the MPAS Framework includes tools for generating uniform and variable resolution centroidal Voronoi meshes. Additionally, the JIGSAW(GEO) mesh generation tool (Engwirda, 2017a, b) can be used to generate high quality variable resolution meshes with data-based density functions very efficiently.
Albany is an open source, C++ multi-physics code base for the solution and analysis of coupled systems of partial-differential equations (PDEs) . It is a finite element code that can (in three spatial dimensions) employ unstructured meshed comprised of hexahedral, tetrahedral, or prismatic elements. Albany is designed to take advantage of the computational mathematics tools available within the Trilinos suite of software libraries (Heroux et al., 2005) and it uses template-based 5 generic programming methods to provide extensibility and flexibility (Pawlowski et al., 2012). Together, Albany and Trilinos provide parallel data structures and I/O, discretization and integration algorithms, linear solvers and preconditioners, nonlinear solvers, continuation algorithms, and tools for automatic differentiation (AD) and optimization. By formulating a system of equations in the residual form, Albany employs AD to automatically compute the Jacobian matrix, as well as forward and adjoint sensitivities. Albany can solve large-scale PDE-constrained optimization problems, using Trilinos optimization pack-10 age ROL, and it provides uncertainty quantification capabilities through the Dakota framework (Adams et al., 2013). It is a massively parallel code by design and recently it has been adopting the Kokkos (Edwards et al., 2014a) programming model to provide manycore performance portability (Demeshko et al., 2018) on major HPC platforms. Albany provides several applications including LCM (Laboratory for Computational Mechanics) for solid mechanics problems, QCAD (Quantum Computer Aided Design) for quantum device modeling, and LI (Land Ice) for modeling ice sheet flow. We refer to the code that discretizes 15 these diagnostic momentum balance equations as Albany-LI. Albany-LI was formerly known as Albany/FELIX (Finite Elements for Land Ice eXperiments), and described by Tezaur et al. (2015a, b) and Tuminaro et al. (2016) under that name. Here, these tools are brought to bear on the most complex, expensive, and fragile portion of the ice sheet model, the solution of the momentum balance equations (discussed further below).

20
The "dynamical core" of the the MALI ice sheet model solves the governing equations expressing the conservation of momentum, mass, and energy.

Conservation of Momentum
Treating glacier ice as an incompressible fluid in a low-Reynolds number flow, the conservation of momentum in a Cartesian reference frame is expressed by the Stokes-flow equations, for which the gravitational-driving stress is balanced by gradients 25 in the viscous stress tensor, σ ij : where x i is the coordinate vector, ρ is the density of ice, and g = is acceleration due to gravity 2 .
2 In Equation 2 and elsewhere we use indicial notation, with summation over repeat indices.
Deformation results from the deviatoric stress, τ ij , which relates to the full stress tensor as for which − 1 3 σ kk is the mean compressive stress and δ ij is the Kroneker delta (or the identity tensor). Stress and strain rate are related through the constitutive relation,

5
where˙ ij is the strain-rate tensor and η e is the "effective", non-Newtonian ice viscosity given by Nye's generalization of Glen's flow law (Glen, 1955), In Equation 5, A is a temperature dependent rate factor, n is an exponent commonly taken as 3 for polycrystalline glacier ice, and γ is an ice "stiffness" factor (inverse enhancement factor) commonly used to account for other impacts on ice rheology, 10 such as impurities or crystal anisotropy (see also Section 6.1). The effective strain rate˙ e is given by the second invariant of the strain-rate tensor, The strain rate tensor is defined by gradients in the components of the ice velocity vector u i : 15 Finally, the rate-factor A follows an Arrhenius relationship in which A o is a constant, (T * ) is the absolute temperature (i.e., corrected for the dependence of melt temperature on ice pressure), Q a is the activation energy for crystal creep, and R is the gas constant.
Boundary conditions required for the solution of Eq. 2 depend on the form of reduced-order approximation applied and are 20 discussed further below.

Reduced-order Equations
Ice sheet models solve Eq. 2-8 with varying degrees of complexity in terms of the tensor components in Eq. 2-7 that are accounted for or omitted, based on geometric scaling arguments. Because ice sheets are inherently "thin" -their widths are several orders of magnitude larger than their thickness -reduced-order approximations of the full momentum balance are 25 often appropriate (see, e.g., Dukowicz et al., 2010;Schoof and Hewitt, 2013) and, importantly, can often result in considerable computational cost savings. Here, we employ two such approximations, a first-order-accurate "Blatter-Pattyn" approximation and a zero-order, "shallow-ice approximation" as described in more detail in the following sections.

First-Order Velocity Solver and Coupling
Ice sheets typically have a small aspect ratio, small surface and bed slopes, and vertical pressure distributions that are very nearly hydrostatic. These characteristics imply that reduced-order approximations of the Stokes momentum balance may apply over large areas of the ice sheets, potentially allowing for significant computational savings. Formal derivations involve non-dimensionalizing the Stokes momentum balance and introducing a geometric scaling factor, δ = H/L, where H and L rep-5 resent characteristic vertical and horizontal length scales (often taken as the ice thickness and the ice sheet span), respectively.
Upon conducting an asymptotic expansion, reduced-order models with a chosen degree of accuracy (relative to the original Stokes flow equations) can be derived by retaining terms of the appropriate order in δ. For example, the first-order accurate Stokes approximation is arrived at by retaining terms of O(δ 1 ) and lower (The reader is referred to Schoof and Hindmarsh (2010) and Dukowicz et al. (2010) for additional discussion 3 ). 10 Using the notation of Perego et al. (2012) and Tezaur et al. (2015a) 4 , the first-order accurate Stokes approximation (also referred to as the "Blatter-Pattyn" approximation, see Blatter, 1995;Pattyn, 2003) is expressed through the following system of PDEs, where ∇· is the divergence operator, s ≡ s(x, y) represents the ice sheet upper surface, and the vectors˙ 1 and˙ 2 are given by Akin to Equations 5 and 6, η f in Equation 9 represents the effective viscosity but for the case of the first-order stress balance with an effective-strain rate is given by At the upper surface, a stress-free boundary condition is applied, 25 3 In practice, additional scaling parameters describing the ratio of deformation to sliding velocity may also be introduced. 4 Vectors and tensors are given in bold rather than using indices. Note that, in a slight abuse of notation, we have switched from using x 1 , x 2 , x 3 to denote the three coordinate directions to x, y, z.
with n the outward normal vector at the ice sheet surface, z = s(x, y). At the bed, z = b(x, y), we apply no slip or continuity of basal tractions ("sliding"), where β is a linear-friction parameter and m ≥ 1. In most applications we set m = 1 (see also Section 5.1.6).
On lateral boundaries, a stress boundary condition is applied, where ρ o is the density of ocean water and n the outward normal vector to the lateral boundary (i.e., parallel to the (x, y) plane), so that lateral boundaries above sea level are effectively stress free and lateral boundaries submerged within the ocean experience hydrostatic pressure due to the overlying column of ocean water.
We solve these equations using the Albany-LI momentum balance solver, which is built using the Albany and Trilinos 10 software libraries discussed above. The mathematical formulation, discretization, solution methods, verification, and scaling of Albany-LI are discussed in detail in Tezaur et al. (2015a). Albany-LI implements a classic finite element discretization of the first-order approximation. At the grounding line, the basal friction coefficient β can abruptly drop to zero within an element of the mesh. This discontinuity is resolved by using an higher-order Gauss quadrature rule on elements containing the grounding line, which corresponds to the sub-element parametrization SEP3 proposed in Seroussi et al. (2014). Additional exploration of solver scalability and demonstrations of solver robustness on large scale, high-resolution, realistic problems are discussed in Tezaur et al. (2015b). The efficiency and robustness of nonlinear solvers are achieved using a combination of the Newton method (damped with a line search strategy when needed) and of a parameter continuation algorithm for the numerical regularization of the viscosity. The scalability of linear solvers is obtained using a multilevel preconditioner (see Tuminaro et al. (2016)) specifically designed to target shallow problems characterized by meshes extruded in the vertical dimension, like 20 those found in ice sheet modeling. The preconditioner has been demonstrated to be particularly effective and robust even in the presence of ice shelves that typically lead to highly ill-conditioned linear systems. Because the momentum balance solver is ≥ 95% of the cost of a typical forward model time step, the model performance reported on in Tezaur et al. (2015a, b) and Tuminaro et al. (2016) is generally representative of overall MALI performance.
The Albany-LI first-order velocity solver written in C++ is coupled to MPAS written in Fortran using an interface layer. 25 Albany uses a three-dimensional mesh extruded from a basal triangulation and composed of prisms or tetrahedra (see Tezaur et al. (2015a)). When coupled to MPAS, the basal triangulation is the part of the Delaunay triangulation, dual to an MPAS Voronoi mesh, that contains active ice and it is generated by the interface. Bed topography, ice lower surface, ice thickness, basal friction coefficient (β), and three-dimensional ice temperature, all at cell centers (Table 1), are passed from MPAS to Albany. Optionally, Dirichlet velocity boundary conditions can also be passed. After the velocity solve is complete, Albany   30 returns the x and y components of velocity at each cell center and blue layer interfaces, the normal component of velocity at each cell edge and blue layer interfaces, and viscous dissipation at blue each cell vertex and layer midpoints. The interface code defines the lateral boundary conditions on the finite element mesh that Albany will use. Lateral boundaries in Albany are applied at cell centers (triangle nodes) that do not contain dynamic ice on the MPAS mesh and that are adjacent to the last cell of the MPAS mesh that does contain dynamic ice. This one element extension is required to support calculation of normal velocity on edges (u n ) required for advection of ice out of the final cell containing dynamic ice ( Figure 2). The interface identifies three types of lateral boundaries for the first-order velocity solve: terrestrial, floating marine, and grounded 5 marine. Terrestrial margins are defined by bed topography above sea level. At these boundary nodes, ice thickness is set to a small ice minimum thickness value ( = 1 m). Floating marine margin triangle nodes are defined as neighboring one or more triangle edges that satisfy the hydrostatic floatation criterion. At these boundary nodes, we need to ensure the existence of a realistic calving front geometry, so we set ice thickness to the minimum of thickness at neighboring cells with ice. Grounded marine margins are defined as locations where the bed topography is below sea level, but no adjacent triangle edges satisfy 10 the floatation criterion. At these boundary nodes, we apply a small floating extension with thickness . For all three boundary types, ice temperature is averaged from the neighboring locations containing ice.

Shallow-Ice Approximation Velocity Solver
A similar procedure to that described above for the first-order accurate Stokes approximation can be used to derive the socalled "shallow-ice approximation" (SIA) (Hutter, 1983;Fowler and Larson, 1978;Morland and Johnson, 1980;Payne et al., 15 2000), in this case by retaining only terms of O(δ 0 ). In the case of the SIA, the local gravitational driving stress is everywhere balanced by the local basal traction and the horizontal velocity as a function of depth is simply the superposition of the local basal sliding velocity and the integral of the vertical shear from the ice base to that depth: where b is the bed elevation and u b is the sliding velocity.

20
SIA ice sheet models typically combine the momentum and mass balance equations to evolve the ice geometry directly in the form of a depth-integrated, two-dimensional diffusion problem (Hindmarsh and Payne, 1996;Payne et al., 2000). However we implement the SIA as an explicit velocity solver that can be enabled in place of the more accurate first order solver, while keeping the rest of the model identical. The purpose of the SIA velocity solver is primarily for rapid testing, so the less efficient explicit implementation of Eq. 17 is not a concern.
We implement Eq. 17 in sigma coordinates on cell edges, where we only require the normal component of velocity, u n : where x n is the normal direction to a given edge and u bn is sliding velocity in the normal direction to the edge. We average 5 A and H from cell centers to cell edges. ds dxn is calculated as the difference in surface elevation between the two cells that neighbor a given edge divided by the distance between the cell centers; on a Voronoi grid, cells edges are midway between cell centers by definition. The surface slope component tangent to an edge (required to complete the calculation of ∇s) is calculated by first interpolating surface elevation from cell centers to vertices.

10
Conservation of mass is used to conduct ice sheet mass transport and evolution. Assuming constant density to write conservation of mass in volume form, the equation relates ice thickness change to the divergence of mass and sources and sinks: where H is ice thickness, t is time,ū is depth-averaged velocity,ȧ is surface mass balance, andḃ is basal mass balance. Botḣ a andḃ are positive for ablation and negative for accumulation.
Eq. 19 is used to update thickness in each grid cell on each time step using a forward Euler, fully explicit time evolution scheme. Eq. 19 is implemented using a finite volume method, such that fluxes are calculated for each edge of each cell to calculate ∇ · Hū. Specifically, we use a first-order upwind method that applies the normal velocity on each edge (u n ) and an 5 upwind value of cell centered ice thickness. Note that with the First Order velocity solver, normal velocity is interpolated from cell centers to edges using the finite element basis functions in Albany. In the shallow ice approximation velocity solver, normal velocity is calculated natively at edges. MPAS Framework includes a higher-order flux-corrected transport scheme  for which we have performed some initial testing, but is not routinely used in the Land Ice core at this time.
Tracers are advected layer by layer with a similar equation: where Q t is a tracer quantity (e.g., temperature -see below), l is layer thickness, andṠ represents any tracer sources or sinks. While any number of tracers can be included in the model, the only one to be considered here is temperature, due to its important effect on ice rheology through Eq. 8 and will be discussed further in the following section.
Because we employ a sigma vertical coordinate system with fixed layer fractions, after Eqs. 19 and 20 are applied, a vertical 15 remapping operation is required. Overlaps between the newly calculated layers and the target sigma layers are calculated for each grid cell. Assuming uniform values within each layer, mass, energy, and other tracers are transferred between layers based on these overlaps to restore the prescribed sigma layers while conserving mass and energy.

Conservation of Energy
Conservation of energy is expressed through the three-dimensional, advective-diffusive heat equation: with thermal conductivity k and heat capacity c. In Equation 21, the rate of temperature change (left-hand side) is balanced by diffusive, advective, and internal (viscous dissipation -see Equation 28 for Φ) source terms (first, second, and third terms on the right-hand side, respectively). In MALI we solve an approximation of Equation 21, 25 in which horizontal diffusion is assumed negligible (van der Veen, 2013, p. 280) and k is assumed constant and uniform. The viscous dissipation term Φ is discussed further below (Section 4.4.2).
Temperatures are staggered in the vertical relative to velocities and are located at the centers of nz − 1 vertical layers, which are bounded by nz vertical levels (grid point locations). This convention allows for conservative temperature advection, since the total internal energy in a column (the sum of ρcT ∆z over nz − 1 layers) is conserved under transport. The upper surface temperature T s and the lower surface temperature T b , coincident with the surface and bed grid points, give a total of nz + 1 temperature values within each column.
Equation 22 is solved using an operator splitting technique. At each time step, we first perform an implicit vertical solve accounting for the diffusion and dissipation terms, and then we explicitly advect horizontally the resulting temperature.

5
Using a "sigma" vertical coordinate, the vertical diffusion portion of Equation 22 can be discretized as: In σ-coordinates, the central difference formulas for first partial derivatives at the upper and lower interfaces of layer k are whereσ k is the value of σ at the midpoint of layer k, halfway between σ k and σ k+1 . The second partial derivative, defined at 10 the midpoint of layer k, is then given by By inserting Equation (24) into Equation (25), we obtain the discrete form of the vertical diffusion term in Equation 22: . (26) To simplify some expressions below, we define the following coefficients associated with the vertical temperature diffusion,

Viscous Dissipation
The source term from viscous dissipation in Equation 22 is given by the product of the stress and strain rate tensors: The change to deviatoric stress on the right-hand side of Equation 28 follows from terms related to the mean compressive stress (or pressure) dropping out due to incompressibility. Analogous to the effective strain rate given in Equation 6, the effectivedeviatoric stress is given by 25 which can be combined with Equations 28 and 6 to derive an expression for the viscous dissipation in terms of effective deviatoric stress and strain, Finally, an analog to Equation 4 gives 5 which can be used to eliminate˙ e in Equation 30 and arrive at an alternate expression for the dissipation based on only two scalar quantities The viscous dissipation source term is computed within Albany-LI at MPAS cell vertices and then reconstructed at cell centers in MPAS. 10 For the SIA model, dissipation can be calculated in sigma coordinates as which can be combined with Eq. 17 to make: We calculate Φ on cell edges following the procedure described for Eq. 18, and then interpolate Φ back to cell centers to solve 15 Eq. 22.

Vertical Temperature Solution
The vertical diffusion portion of Equation 22 is discretized according to where a k and b k are defined in (27), n is the current time level, and n + 1 is the new time level. Because the vertical diffusion 20 terms are evaluated at the new time level, the discretization is backward-Euler (fully implicit) in time.
The temperature T 0 at the upper boundary is set to min(T air , 0), where the mean-annual surface air temperature T air is a two-dimensional field specified from observations or climate model output.
At the lower boundary, for grounded ice there are three potential heat sources and sinks: (1) the diffusive flux from the bottom surface to the ice interior (positive up), (2) the geothermal flux F g , prescribed from a spatially variable input file (based on observations), and (3) the frictional heat flux associated with basal sliding, where τ b and u b are 2D vectors of basal shear stress and basal velocity, respectively, and the friction law from Equation 15 becomes If the basal temperature T nz < T pmp (where T pmp is the pressure melting point temperature), then the fluxes at the lower boundary must balance, so that the energy supplied by geothermal heating and sliding friction is equal to the energy removed by vertical diffusion. If, 10 on the other hand, T nz = T pmp , then the net flux is nonzero and is used to melt or freeze ice at the boundary: where M b is the melt rate and L is the latent heat of melting. Melting generates basal water, which may either be stored at the bed locally, serve as a source for the basal hydrology model (See Section 5.1), or may simply be ignored. If basal water is present locally, T nz is held at T pmp .

15
For floating ice the basal boundary condition is simpler: T nz is simply set to the freezing temperature T f of seawater.
Optionally, a melt rate can be prescribed at the lower surface.
Rarely, the solution for T may exceed T pmp for a given internal layer. In this case, T is set to T pmp , excess energy goes towards melting of ice internally, and the resulting melt is assumed to drain to the bed immediately.
If (40) applies, we compute M b and adjust the basal water depth. When the basal water goes to zero, T nz is set to the 20 temperature of the lowest layer (less than T pmp at the bed) and flux boundary conditions apply during the next time step.

Horizontal Advection
Temperature advection in any individual layer k is treated using tracer advection, as in Equation 20 above, where the ice temperature T k is substituted for the generic tracer Q. After horizontal transport, the surface and basal mass balance is applied to the top and bottom ice surfaces, respectively. Because layer transport and the application of mass balance terms results in 25 an altered vertical-layer spacing with respect to σ coordinates, a vertical remapping scheme is applied. This conservatively transfers ice volume and internal energy between adjacent layers while restoring σ layers to their initial distribution. Internal energy divided by mass gives the new layer temperatures.
Physical processes currently implemented in MALI are a mass-conserving subglacial hydrology model and a small number of basic schemes for iceberg calving. These are described in more detail below.

Subglacial Hydrology
Sliding of glaciers and ice sheets over their bed can increase ice velocity by orders of magnitude and is the primary control on 5 ice flux to the oceans. The state of the subglacial hydrologic system is the primary control on sliding (Clarke, 2005;Cuffey and Paterson, 2010;Flowers, 2015), and ice sheet modelers have therefore emphasized subglacial hydrology and its effects on basal sliding as a critical missing piece of current ice sheet models (Little et al., 2007;Price et al., 2011).
MALI includes a mass-conserving model of subglacial hydrology that includes representations of any or all of water storage in till, distributed drainage, and channelized drainage and is coupled to ice dynamics. The model is based on the model of 10 Bueler and van Pelt (2015) with an additional component for channelized drainage and modified for MALI's unstructured horizontal grid. While the implementation follows closely that of Bueler and van Pelt (2015), the model and equations are summarized here along with a description of the features unique to the application in MALI.

Till
The simple till component represents local storage of water in subglacial till without horizontal transport within the till. Evolution of the effective water depth in till, W till is therefore a balance of delivery of meltwater, m b , to the till, drainage of water out of the till at rate C d (mass leaving the subglacial hydrologic system, for example, to deep groundwater storage), and overflow to the distributed drainage system, γ: In the model, meltwater (from either the bed or drained from the surface), is first delivered to the till component. Water in 20 excess of the the maximum storage capacity of the till, W max till , is instantaneously transferred as a source term to the distributed drainage system through the γ t term.

Distributed drainage
The distributed drainage component is implemented as a "macroporous sheet" that represents bulk flow through linked cavities that form in the lee of bedrock bumps as the glacier slides over the bed (Flowers and Clarke, 2002;Hewitt, 2011;Flowers, 25 2015). Water flow in the system is driven by the gradient of the hydropotential, φ, defined as where P w is the water pressure in the distributed drainage system. A related variable, the ice effective pressure, N , is the difference between ice overburden pressure and water pressure in the distributed drainage system, P w : The evolution of the area-averaged cavity space is a balance of opening of cavity space by the glacier sliding over bedrock bumps and closing through creep of the ice above. The model uses the commonly used assumption (e.g. Schoof, 2010;Hewitt, 5 2011; Werder et al., 2013;Hoffman and Price, 2014) that cavities always remain water filled (c.f. Schoof et al., 2012), so cavity space can be represented by the effective water depth in the macorporous sheet, W : where c s is bed roughness parameter, W r is the maximum bed bump height, c cd is creep scaling parameter representing geometric and possibly other effects, and A b is the ice flow parameter of the basal ice. 10 Water flow in the distributed drainage system, q, is driven the hydropotential gradient and is described by a general powerlaw: where k q is a conductivity coefficient. The α 1 and α 2 exponents can be adjusted so that Eq. 45 reduces to commonly used water flow relations, such as Darcy flow, the Darcy-Weisbach relation, and the Manning equation.

Channelized drainage
The inclusion of channelized drainage in MALI is an extension to the model of Bueler and van Pelt (2015). The distributed drainage model ignores dissipative heating within the water, which in the real world leads to melting of the ice roof, and the formation of discrete, efficient channels melted into the ice above when the distributed discharge reaches a critical threshold (Schoof, 2010;Hewitt, 2011;Werder et al., 2013;Flowers, 2015). These channels can rapidly evacuate water from the dis-20 tributed drainage system and lower water pressure, even under sustained meltwater input (Schoof, 2010;Hewitt, 2011;Werder et al., 2013;Hoffman and Price, 2014;Flowers, 2015).
The implementation of channels follows the channel network models of Werder et al. (2013) and Hewitt (2013). The evolution of channel area, S, is a balance of opening and closing processes as in the distributed system, but in channels the opening mechanism is melting caused by dissipative heating of the ice above: where c cc is the creep scaling parameter for channels.
The channel opening rate, the first term in Eq. 46, is itself a balance of dissipation of potential energy, Ξ, and sensible heat change of water, Π, due to changes in the pressure-dependent melt temperature. Dissipation of potential energy includes energy produced by flow in both the channel itself and a small region of the distributed system along the channel: where s is the spatial coordinate along a channel segment, Q is the flow rate in the channel, and q c is the flow in the distributed drainage system parallel to the channel within a distance l c of the channel. The term adding the contribution of dissipative melting within the distributed drainage system near the channel is included to represent some of the energy that has been 5 ignored from that process in the description of the distributed drainage system and allows channels to form even when channel area is initially zero if discharge in the distributed drainage system is sufficient (Werder et al., 2013). The term representing sensible heat change of the water, Π, is necessitated by the assumption that the water always remains at the pressure-dependent melt temperature of the water. Changes in water pressure must therefore result in melting or freezing: where c t is the Clapeyron slope and c w is the specific heat capacity of water. The pressure-dependent melt term can be disabled in the model.
Water flow in channels, Q, mirrors Eq. 45: where k Q is a conductivity coefficient for channels.

Drainage component coupling
Eqs. 41-49 are coupled together by describing the drainage system with two equations, mass conservation and pressure evolution. Mass conservation of the subglacial drainage system is described by where V d is water velocity in the distributed flow, D d is the diffusivity of the distributed flow, and δ(x c ) is the Dirac delta 20 function applied along the locations of the linear channels.
Combining Eq. 50 and Eq. 44 and making the simplification that cavities remain full at all times yields an equation for water pressure within the distributed drainage system, P w : where φ 0 is an englacial porosity used to regularize the pressure equation. Following Bueler and van Pelt (2015), the porosity 25 is only included in the pressure equation and is excluded from the mass conservation equation.
Any of the three drainage components (till, distributed drainage, channelized drainage) can be deactivated at runtime. The most common configuration currently used is to run with distributed drainage only.

Numerical implementation
The drainage system model is implemented using Finite Volume Methods on the unstructured grid used by MALI. State variables (W, W till , S, P w ) are located at cell centers and velocities and fluxes (q, V d , Q) are calculated at edge midpoints.
Channel segments exist along the lines joining neighboring cell centers. Eq. 50 is evaluated by summing tendencies from discrete fluxes into or out of each cell. First-order upwinding is used for advection. At land-terminating ice sheet boundaries, 5 P w = 0 is applied as the boundary condition. At marine-terminating ice sheet boundaries, the boundary condition is P w = −ρ w gz b , where ρ w is ocean water density. The drainage model uses explicit forward Euler time-stepping using Eqs. 41, 50, 46, and 51. This requires obeying advective and diffusive Courant-Friedrichs-Lewy (CFL) conditions for distributed drainage as described by Bueler and van Pelt (2015), as well as an additional advective CFL condition for channelized drainage, if it is active. 10 We acknowledge that the non-continuum implementation of channels can make the solution grid-dependent, and grid convergence may therefore not exist for many problems (Bueler and van Pelt, 2015). However, for realistic problems with irregular bed topography, we have found dominant channel location is controlled by topography, mitigating this issue.

Coupling to ice sheet model
The subglacial drainage model is coupled to the ice dynamics model through a basal friction law. Currently, the only option is 15 a modified Weertman-style power law (Bindschadler, 1983;Hewitt, 2013) that adds a term for effective pressure to Eq. 15: where C 0 is a friction parameter. Implementation of a Coulomb friction law (Schoof, 2005;Gagliardini et al., 2007) and a plastic till law (Tulaczyk et al., 2000;Bueler and van Pelt, 2015) are in development. When the drainage and ice dynamics components are run together, coupling of the systems allows the negative feedback described by (Hoffman and Price, 2014) 20 where elevated water pressure increases ice sliding and increased sliding opens additional cavity space, lowering water pressure.
The meltwater source term, m, is calculated by the thermal solver in MALI. Either or both of the ice dynamics and thermal solvers can be disabled, in which case the relevant coupling fields can be prescribed to the drainage model.

Verification and Real-world Application
To verify the implementation of the distributed drainage model, we use the nearly exact solution described by Bueler and van 25 Pelt (2015). The problem configuration uses distributed drainage only on a two-dimensional, radially-symmetric ice sheet of radius 22.5 km with parabolic ice sheet thickness and a nontrivial sliding profile. Bueler and van Pelt (2015) showed that this configuration allows for nearly-exact reference values of W and P w to be solved at steady state from an ordinary differential equation initial value problem with very high accuracy. We follow the test protocol of Bueler and van Pelt (2015) and initialize the model with the near-exact solution and then run the model forward for one month, after which we evaluate model error due to drift away from the expected solution. Performing this test with the MALI drainage model, we find error comparable to that found by Bueler and van Pelt (2015) and approximately first order convergence (Figure 3).
To check the model implementation of channels, we use comparisons to other, more mature drainage models through the Subglacial Hydrology Model Intercomparison Project (SHMIP) 5 . Steady state solutions of the drainage system effective pressure, water fluxes, and channel development for an idealized ice sheet with varying magnitudes of meltwater input (SHMIP 5 experiment suites A and B) compared between MALI and other models of similar complexity (GlaDS, Elmer) are very similar.
To demonstrate a real-world application of the subglacial hydrology model, we perform a standalone subglacial hydrology simulation of the entire Antarctica Ice Sheet on a uniform 20 km resolution mesh (Figure 4). We force this simulation with basal sliding and basal melt rate after optimizing the first-order velocity solver optimized to surface velocity observations ( Figure   4a). We then run the subglacial hydrology model to steady state with only distributed drainage active and using standard 10 parameter values. Results show a subglacial hydrologic state that is reasonable. For example, the subglacial water flux shows similar spatial patterns to the basal sliding speed (Figure 4), suggesting a basal friction law based on the subglacial hydrologic state could be configured to yield realistic ice velocity. Calibrating parameters for the subglacial hydrology model and a basal friction law and performing coupled subglacial-hydrology/ice-dynamics simulations are beyond the scope of this paper; we merely mean to demonstrate plausible behavior from the subglacial hydrology model for a realistic ice-sheet-scale problem.

5
MALI includes a few simple methods for removing ice from calving fronts during each model time step: 1. All floating ice is removed.
2. All floating ice in cells with an ocean bathymetry deeper than a specified threshold is removed.
3. All floating ice thinner than a specified threshold is removed. 4. The calving front is maintained at its initial location by adding or removing ice after thickness evolution is complete. 10 This option does not conserve mass or energy but provides a simple way to maintain a realistic ice shelf extent (e.g., for model spinup). 5. "Eigencalving" scheme (Levermann et al., 2012). Calving front retreat rate, C v , is proportional to the product of the principal strain rates (˙ 1 ,˙ 2 ) if they both are extensional: The eigencalving scheme can optionally also remove floating ice at the calving front with thickness below a specified thickness threshold (Feldmann and Levermann, 2015). In practice we find this is necessary to prevent formation of 5 tortuous ice tongues and continuous, gradual extension of some ice shelves along the coast.
Ice that is eligible for calving can be removed immediately or fractionally each time step based on a calving timescale. To allow ice shelves to advance as well as retreat, we implement a simple parameterization for sub-grid motion of the calving front by forcing floating cells adjacent to open ocean to remain dynamically inactive until ice thickness there reaches 95% of the minimum thickness of all floating neighbors. This is an ad hoc alternative to methods tracking the calving front position 10 at sub-grid scales (Albrecht et al., 2011;Bondzio et al., 2016). In Section 8 below, we demonstrate the eigencalving scheme applied to a realistic Antarctic ice sheet simulation. More sophisticated calving schemes are currently under development.
6 Additional Capabilities

Optimization
MALI includes an optimization capability through its coupling to the Albany-LI momentum-balance solver described in Sec-15 tion 4.2.1. We provide a brief overview of this capability here, while referring to Perego et al. (2014) for a complete description of the governing equations, solution methods, and example applications. In general, our approaches are similar to those reported on for other advanced ice sheet modeling frameworks already described in the literature (e.g., Goldberg and Sergienko, 2011;Larour et al., 2012;Gagliardini et al., 2013;Brinkerhoff and Johnson, 2013;Cornford et al., 2013) and focus (here) primarily on optimizing the model velocity field relative to observed surface velocities. Briefly, we consider the optimization where the first term on the right-hand side (RHS) is a cost function associated with the misfit between modeled and observed surface velocities, the second term on the RHS is a cost function associated with the ice stiffness factor, γ (see Equation 5), and the third and fourth terms on the RHS are Tikhonov regularization terms given by σ u is an estimate for the the standard deviation of the uncertainty in the observed ice surface velocities and the parameter c γ controls how far the ice stiffness factor is allowed to stray from unity in order to improve the match to observed surface velocities. The regularization parameters α β > 0 and α γ > 0 control the tradeoff between a smooth β field and one with higherfrequency oscillations (that may capture more spatial detail at the risk of over-fitting the observations). The optimal values of α β and α γ can be chosen through a standard L-curve analysis 6 . The optimization problem is solved using the Limited Memory BFGS method, as implemented in the Trilinos package ROL 7 , on the reduced-space problem. The functional gradient is computed using the adjoint method.

5
An example application of the optimization capability applied to a realistic, whole-ice-sheet problem is given below in Section 8. Hoffman et al. (2018) present another application to the assimilation of surface velocity time series in Western Greenland.
We note that our optimization framework has been designed to be significantly more general than implied by Equation 54.
While not applied here, we are able to introduce additional observational-based constraints (e.g., mass balance terms) and 10 optimize additional model variable (e.g., the ice thickness). These are necessary, for example, when targeting model initial conditions that are in quasi-equlibrium with some applied climate forcing. These capabilities are discussed in more detail in Perego et al. (2014).

Simulation Analysis
As with other climate model components built using the MPAS framework, MALI supports the development and application 15 of "analysis members", which allow for a wide range of run-time-generated simulation diagnostics and statistics output at user specified time intervals. Support tools included with the code release allow for the definition of any number or combination of pre-defined "geographic features" -points, lines ("transects"), or areas ("regions") -of interest within an MPAS mesh.
Features are defined using the standard GeoJSON format (Butler et al., 2016) and a large existing database of globally defined features is currently supported 8 . Python-based scripts are available for editing GeoJSON feature files, combining or splitting 20 them, and using them to define their coverage within MPAS mesh files. Currently, MALI includes support for standard ice sheet model diagnostics (see Table 2) defined over the global domain (by default) and / or over specific ice sheet drainage basins and ice shelves (or their combination). Support for generating model output at points and along transects will be added in the future (e.g., vertical samples at ice core locations or along ground-penetrating radar profile lines). In Section 8 below we demonstrate the analysis capability applied to an idealized simulation of the Antarctica ice sheet.

Model Verification and Benchmarks
MALI has been verified by a series of configurations that test different components of the code. In some cases analytic solutions are used, but other tests rely on intercomparison with community benchmarks that have been run previously by many different ice sheet models. 6 We note that, while not done here, additional scalar parameters could be placed in front of R(β) and R(γ) in Equation 54 to differentially weight their contributions to the cost function 7 https://trilinos.org/packages/rol/ 8 https://github.com/MPAS-Dev/geometric_features

Halfar analytic solution
In (Halfar, 1981(Halfar, , 1983, Halfar described an analytic solution for the time-evolving geometry of a radially-symmetric, isothermal dome of ice on a flat bed with no accumulation flowing under the shallow ice approximation. This test provides an obvious test of the implementation of the shallow-ice velocity calculation and thickness evolution schemes in numerical ice sheet mod-10 els and a way to assess model order of convergence Egholm and Nielsen, 2010).  showed the Halfar test is the zero accumulation member of a family of analytic solutions, but we apply the original Halfar test here.
In our application we use a dome following the analytic profile prescribed by Halfar (1983) with an initial radius of 21213.2 m and an initial height of 707.1 m. We run MALI with the shallow ice velocity solver and isothermal ice for 200 years and then compare the modeled ice thickness to the analytic solution at 200 years. We find the root mean square error in model thickness decreases as model grid spacing is decreased (Figure 5a). The order of convergence of 0.78, consistent with the first-order approximation used for advection.
We also use this test to assess the accuracy of simulations with variable resolution. We perform an addition run of the Halfar test using a variable resolution mesh that has 1000 m cell spacing beyond a radius of 20 km that transitions to 5000 m 5 cell spacing at a radius of 3 km (Figure 5b), generated with the JIGSAW(GEO) mesh generation tool (Engwirda, 2017a, b).
Root mean square error in thickness for this simulation is similar to that for the uniform 1000 m resolution case (Figure 5a), providing confidence in the advection scheme applied to variable resolution meshes. The variable resolution mesh has about half the cells of the 1000 m uniform resolution mesh.

EISMINT
European Ice Sheet Modeling Initiative (EISMINT) model intercomparison consisted of two phases designed to provide community benchmarks for shallow-ice models. Both phases included experiments that grow a radially symmetric ice sheet on a flat bed to steady state with a prescribed surface mass balance. The EISMINT intercomparisons test ice geometry evolution and ice temperature evolution with a variety of forcings. Bueler et al. (2007) describe an alternative tool for testing thermo-5 mechanical shallow-ice models with artificially constructed exact solutions. While their approach has the notable advantage of providing exact solutions, we have not implemented the non-physical three-dimensional compensatory heat source necessary for its implementation. While we hope to use the verification of Bueler et al. (2007) in the future, for now we use the EISMINT intercomparison suites to test our implementation of thermal evolution and thermomechanical coupling.
The first phase (Huybrechts et al., 1996) (sometimes called EISMINT1) prescribes evolving ice geometry and temperature, 10 but the flow rate parameter A is set to a prescribed value so there is no thermomechanical coupling. We have conducted the Moving Margin experiment with steady surface mass balance and surface temperature forcing. Following the specifications described by Huybrechts et al. (1996), we run the ice sheet to steady state over 200 ka. We use the grid spacing prescribed by Huybrechts et al. (1996) (50 km), but due to the uniform Voronoi grid of hexagons we employ, we have a slightly larger number of grid cells in our mesh (1080 vs 961). At the end of the simulation, the modeled ice thickness at the center of the dome by 15 MALI is 2976.7 m, compared with a mean of 2978.0 ± 19.3 m for the ten three-dimensional models reported by Huybrechts et al. (1996). MALI achieves similar good agreement for basal homologous temperature at the center of the dome with a value of −13.09 • C, compared with −13.34 ± 0.56 • C for the six models that reported temperature in Huybrechts et al. (1996).
The second phase of EISMINT (Payne et al., 2000) (sometimes called EISMINT2), uses the basic configuration of the EISMINT1 Moving Margin experiment but activates thermomechanical coupling through Eq. 8. Two experiments (A and F) 20 grow an ice sheet to steady state over 200 ka from an initial condition of no ice, but with different air temperature boundary conditions. Additional experiments use the steady state solution from experiment A (the warmer air temperature case) as the initial condition to perturbations in the surface air temperature or surface mass balance forcings (experiments B, C, and D).
Because these experiments are thermomechanically coupled, they test model ice dynamics and thickness and temperature evolution, as well as their coupling. There is no analytic solution to these experiments, but ten different models contributed 25 results, yielding a range of behavior against which to compare additional models. Here we present MALI results for the five such experiments that prescribe no basal sliding (experiments A, B, C, D, F). Our tests use the same grid spacing as prescribed by Payne et al. (2000) (25 km), again with a larger number of grid cells in our mesh (4464 vs 3721). Payne et al. (2000) report results for five basic glaciological quantities calculated by ten different models, which we have summarized here with the corresponding values calculated by MALI (Table 3). All MALI results fall within the range of 30 previously reported values, except for volume change and divide thickness change in experiment C and melt fraction change in experiment D. However, these discrepancies are close to the range of results reported by Payne et al. (2000), and we consider temperature evolution and thermomechanical coupling within MALI to be consistent with community models, particularly given the difference in model grid and thickness evolution scheme. Table 3. EISMINT2 results for MALI shallow-ice model. For each experiment, model name "EISMINT2" refers to mean and range of models reported in Payne et al. (2000), where we assume the range reported in by Payne et al. (2000)  A long-studied feature of the EISMINT2 intercomparison is the cold "spokes" that appear in the basal temperature field of all models in Experiment F and, for some models, experiment A (Payne et al., 2000;Saito et al., 2006;Bueler et al., 2007;Brinkerhoff and Johnson, 2013). MALI with shallow-ice velocity exhibits cold spokes for experiment F but not experiment A ( Figure 6). Bueler et al. (2007) argue these spokes are a numerical instability that develops when the derivative of the strain heating term is large. Brinkerhoff and Johnson (2013) demonstrate that the model VarGlaS avoids the formation of these 5 cold spokes. However, that model differs from previously analyzed models in several ways: it solves a three-dimensional, advective-diffusive description of an enthalpy formulation for energy conservation; it uses the Finite Element Method on unstructured meshes; conservation of momentum and energy are iterated on until they are consistent (rather than lagging energy and momentum solutions as in most other models). At present, it is unclear which combination of those features is responsible for preventing the formation of the cold spokes. These tests only require a single diagnostic solve of velocity, and thus results through MALI match those of the standalone Albany-LI code it is using.

MISMIP3d
The Marine Ice Sheet Model Intercomparison Project-3d (MISMIP3d) is a community benchmark experiment testing grounding line migration of marine ice sheets and includes nontrivial effects in all three dimensions . The exper-5 iments use a rectangular domain that is 800 km long in the longitudinal direction and 50 km wide in the transverse direction, with the transverse direction making up half of a symmetric glacier. The bedrock forms a sloping plane below sea level. The first phase of the experiment (Stnd) is to build a steady state ice sheet from a spatially uniform positive surface mass balance, with a prescribed flow rate factor A (no temperature calculation or coupling) and prescribed basal friction for a nonlinear basal friction law. A marine ice sheet forms with an unbuttressed floating ice shelf that terminates at a fixed ice front at the edge of  km. Participating models used depth-integrated shallow-shelf or L1L1/L2L2 approximations, hybrid shallow-ice/shallow-shelf approximation, or the complete Stokes equations; there were no three-dimensional first-order approximation models included.
This relatively simple experiment revealed a number of key features necessary to accurately model even a simple marine ice sheet. Insufficient grid resolution prevented reversibility of the steady state grounding line position after experiments P75S and P75R. Reversibility required grid resolution well below 1 km without a subgrid parameterization of grounding line position, and grids a couple times coarser with a grounding line parameterization Gladstone et al., 2010). The steady state grounding line position in the Stnd experiment was dependent on the stress approximation employed, with Stokes model calculating grounding lines the farthest upstream and models that simplify or eliminate vertical shearing (e.g., shallow shelf) 5 having grounding lines farther downstream, by up to 100 km. With these features resolved, numerical error due to grounding line motion is smaller than errors due to parameter uncertainty .
We find MALI using the Albany-LI first-order velocity solver is able to resolve the MISMIP3d experiments satisfactorily compared to the Pattyn et al. (2013) benchmark results when using a grid resolution of 500 m with grounding line parameterization. Results at 1 km resolution with grounding line parameterization are close to fully resolved. We first assess grid Thus for marine ice sheets with similar configuration to the MISMIP3d test, we recommend using MALI with the grounding line parameterization and a resolution of 1 km or less.
The transient results using the MALI three-dimensional first-order stress balance look most similar to those of the "SCO6" 25 L1L2 model presented by Pattyn et al. (2013) in that it takes about 50 years for the grounding line to reach its most advanced position during P75S. In contrast, the Stokes models took notably longer and the models with reduced or missing representation of membrane stresses reached their furthest advance within the first couple decades .
In addition to MISMIP3d, we have used MALI to perform the MISMIP+ experiments (Asay-Davis et al., 2016). These results are included in the MISMIP+ results paper in preparation and not shown here.  To demonstrate a large scale, semi-realistic ice sheet simulation that exercises many of the model capabilities discussed above, we apply MALI to the InitMIP-Antarctica experiments 9 . The overall goal of InitMIP is to explore the impact of different ice sheet model initialization approaches on simulated ice sheet evolution. Following model initialization -via spin-up, optimization approaches, or both -three 100 year forward model experiments are conducted in order to examine model response to (1) 5 an unforced control run, (2) an idealized surface mass balance perturbation, and (3) an idealized sub-ice shelf melt perturbation.
Additional details of the experiments are described in Seroussi and others (2018) 10 and additional details on the broader Ice Sheet Model Intercomparison Project (ISMIP6, part of CMIP6) that InitMIP is a part of are described in Nowicki et al. (2016) and Goelzer and others (2018).
The Antarctica model configuration we demonstrate here uses a 20 km uniform resolution mesh, an initial internal ice 10 temperature field from Van Liefferinge and Pattyn (2013), and the optimization capability described above in Section 6.1 (note that we optimize both β and γ in this case), along with observed surface velocities from Rignot et al. (2011) to obtain a realistic model initial state (Figure 9a,b). Then, we perform a 1000 yr spin-up with evolving velocity, geometry, and temperature, applying steady forcing of present-day estimates for surface mass balance and submarine melting (Lenaerts et al. (2012) and Rignot et al. (2013), respectively) (Figure 9c,d). For temperature boundary conditions, we apply the steady geothermal flux 15 field from Shapiro and Ritzwoller (2004) and the surface (2 meter) air temperature field from Lenaerts et al. (2012). We use the "eigencalving" calving option (discussed in Section 5.2) with the K 2 parameter tuned by ice shelf and include a minimum calving front thickness threshold of 100 m. This standard spin-up, albeit much too short to come to full to equilibrium, removes a substantial portion of the largest model transients (Figure 9d) and retains an ice sheet extent and ice speed field qualitatively very similar to the observations (Figure 9a,c). The model state at this point serves as our final initial condition from which we 20 run the control and sub-ice shelf melt perturbation experiments mentioned above. Simulation outputs for these two experiments are shown in Figure 10.
In Figure 10 we use the online-analysis capability discussed above in Section 6.2 to present some of the broad features of the initMIP sub-ice shelf melt perturbation experiment. The left column of Figure  the majority of floating ice for the PT system. As expected due to a thinning ice shelf leading to reduced buttressing, the PT system initially shows an increase in flux across the grounding line. A seemingly counterintuitive result for the PT system is the subsequent rapid decrease in grounding-line flux after year 31, which closer examination reveals is due to grounding line 9 Additional realistic applications of MALI to Greenland simulations are discussed in Shannon et al. (2013) and Edwards et al. (2014b). 10 http://www.climate-cryosphere.org/wiki/index.php?title=InitMIP-Antarctica retreat into an area of relatively less slippery bed, whereby velocities at the grounding line slow down by ∼ 50 − 100%. For the Ross and Filchner-Ronne (FR) systems, we see a more muted decrease in grounded area and VAF as a result of the melt perturbation, which is anticipated as the prescribed melt rate perturbations applied to these systems is more than four times smaller. An additional counterintuitive result observed for the Amery system is that, despite a substantial loss of floating ice area (and to a much lesser extent, grounded ice area), the overall VAF for this system increases following the melt perturbation, 5 relative to the control run. A detailed examination indicates that the overall floating and grounded ice area loss for this system occurs in smaller, isolated systems of this region rather than within the main Amery ice shelf. The increase in volume above floatation in the run with increased basal melting relative to the control run can be explained for the same reason as in PTgrounding line retreat into a stickier bed leads to decreased flux across the grounding line and retention of mass within the ice sheet interior. Lastly, we note the strongly stochastic nature of the (annually averaged) calving flux, which varies over time by 10 up to an order of magnitude for some systems. Overall, relative to the control run, calving fluxes are eventually lower for the basal-melt perturbation experiment, which can be explained by the overall decrease in ice shelf thickness under conditions of increased basal melting.
We note that, based on this coarse-resolution model configuration, our non-equilibrated initial conditions, and the idealized climate forcing applied we do not attempt to interpret results from these experiments as predictive of real-world conditions 15 (e.g., grounding line migration will be severely restricted and unrealistic due to our coarse resolution mesh). Rather, our aim here is simply to demonstrate that when applied to large scale, whole-ice-sheet simulations on realistic geometries, MALI is robust and evolves reasonably during multi-century-length, free-running simulations. Future work will include running the complete set of initMIP experiments, taking advantage of the variable-resolution mesh capability to place adequate resolution where needed, like at and near to grounding lines. (E3SM). E3SM is an Earth System Model with atmosphere, land, ocean, and sea ice components, linked through a coupler that passes the necessary fields (e.g., model state, mass, momentum, and energy fluxes) between the components. E3SM, which branched from the Community Earth System Model (CESM, version 1.2 beta10) in 2014, targets high resolution global 25 simulations, and all components have a variable resolution mesh capability. The ocean Petersen et al., 2015Petersen et al., , 2018 and sea ice  components are also built on the MPAS Framework. Because the coupling between E3SM and MALI is currently still fairly rudimentary, we include only a few additional details below and leave a more detailed description to future work. Having all three of these E3SM components in the MPAS framework has simplified adding and maintaining them within E3SM, because developments in the component driver code and build and configuration scripts made  Ongoing and future work on MALI and E3SM coupling includes: passing subglacial discharge at terrestrial ice margins to the land runoff model in E3SM; passing surface runoff calculated in E3SM to the land ice model (for use as a source term in the subglacial hydrology model); two-way coupling between the ocean and a dynamic MALI model 11 ; discharge of icebergs 5 (solid ice flux from MALI) to the coupler and from there to the ocean and sea ice models.

Conclusions and Future Work
We have described MPAS-Albany Land Ice (MALI), a higher-order, thermomechanically coupled ice sheet model using unstructured Voronoi meshes. MALI takes advantage of the MPAS framework for development of unstructured grid Earth system model components and the Albany and Trilinos frameworks for parallel, performance-portable, implicit solution of the chal-lenging higher-order ice sheet momentum balance through the Albany-LI velocity solver. Together, these tools provide an accurate, efficient, scalable, and portable ice sheet model targeted for high resolution ice sheet simulations within a larger Earth system modeling framework run on tens of thousands of computing cores, and MALI makes up the ice sheet component of the Energy Exascale Earth System Model version 1.
MALI includes three-dimensional Blatter-Pattyn and shallow-ice velocity solvers, a standard explicit mass evolution scheme, 5 a thermal solver that can use either a temperature or enthalpy formulation, and an adaptive time stepper. Physical processes represented in the model include subglacial hydrology and calving. The model includes a mass-conserving subglacial hydrology model that can represent combinations of water drainage through till, distributed systems, and channelized systems, and can be coupled to ice dynamics. A handful of basic calving schemes are currently implemented, including the physically based "eigencalving" method. 10 We have demonstrated accuracy of the various model components through commonly used exact solutions and community benchmarks. Of note, we presented the first results for the MISMIP3d benchmark experiments using a Blatter-Pattyn model, and the results are intermediate to those of Stokes and L1L2 models, as might be expected. We also showed simulation results for a semi-realistic Antarctic Ice Sheet configuration at coarse resolution, and this capability was facilitated by the optimization tools within Albany-LI. 15 A number of model improvements are planned over the next five years, focused heavily on improved representation of ice sheet physical processes and Earth system coupling. An implicit subglacial hydrology model based on existing such models (Werder et al., 2013;Hoffman and Price, 2014) is under development using the Albany framework. It will include optimization capabilities, a technique which has yet to be applied to subglacial hydrology beyond a spatially-average, 0-dimension application (Brinkerhoff et al., 2016). The difficulty in subglacial drainage parameter estimation remains one of the primary reasons 20 drainage models have yet to be widely applied in ice sheet models. Improved calving schemes are also under development using a continuum damage mechanics approach (Borstad et al., 2012;Duddu et al., 2013;Bassis and Ma, 2015). Additionally, solid earth processes affecting ice sheets are planned future development, including gravitational, elastic and viscous effects.
Higher-order advection, through the flux-corrected transport  and/or incremental remapping (Dukowicz et al., 2010;Lipscomb and Hunke, 2004;Turner et al., 2018) schemes that are already implemented in MPAS, and semi-25 implicit time-stepping are planned. Finally, a high priority is completing coupling between ice sheet, ocean, and sea ice models in E3SM.

Code availability
MPAS releases are available at https://mpas-dev.github.io/ and model code is maintained at https://github.com/MPAS-Dev/MPAS-Release/releases. MPAS-Albany Land Ice is included in MPAS version 6. The Digital Object Identifier for MPAS v6 is (Zenodo 30 link to be added when release is complete). The Albany library is developed openly at https://github.com/gahansen/Albany, and the Trilinos library is developed openly at https://github.com/trilinos/Trilinos. Region definitions for analysis are openly maintained at https://github.com/MPAS-Dev/geometric_features.