the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
MPASAlbany Land Ice (MALI): a variableresolution ice sheet model for Earth system modeling using Voronoi grids
Mauro Perego
Stephen F. Price
William H. Lipscomb
Tong Zhang
Douglas Jacobsen
Irina Tezaur
Andrew G. Salinger
Raymond Tuminaro
Luca Bertagna
We introduce MPASAlbany Land Ice (MALI) v6.0, a new variableresolution 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 variableresolution Earth system model components and the Albany multiphysics code base for the solution of coupled systems of partial differential equations, which itself makes use of Trilinos solver libraries. MALI includes a threedimensional firstorder momentum balance solver (Blatter–Pattyn) by linking to the AlbanyLI ice sheet velocity solver and an explicit shallow ice velocity solver. The evolution of ice geometry and tracers is handled through an explicit firstorder horizontal advection scheme with vertical remapping. The evolution of ice temperature is treated using operator splitting of vertical diffusion and horizontal advection and can be configured to use either a temperature or enthalpy formulation. MALI includes a massconserving 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 that the calving rate is proportional to extensional strain rates. MALI is evaluated against commonly used exact solutions and community benchmark experiments and shows the expected accuracy. Results for the MISMIP3d benchmark experiments with MALI's Blatter–Pattyn solver fall between published results from Stokes and L1L2 models as expected. We use the model to simulate a semirealistic Antarctic ice sheet problem following the initMIP protocol and using 2 km resolution in marine ice sheet regions. MALI is the glacier component of the Energy Exascale Earth System Model (E3SM) version 1, and we describe current and planned coupling to other E3SM components.
During the past decade, numerical ice sheet models (ISMs) have undergone a renaissance relative to their predecessors. This period of intense model development was initiated following the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2007), which pointed to deficiencies in ISMs of the time as being the single largest shortcoming with respect to the scientific community's ability to project future sea level rise stemming from ice sheets. Model maturation during this period, which continued through the IPCC's Fifth Assessment Report (IPCC, 2013) and to the present day, has focused on improvements to ISM “dynamical cores” (including the fidelity, discretization, and solution methods for the governing conservation equations; e.g., Bueler and Brown, 2009; Schoof and Hindmarsh, 2010; Goldberg, 2011; Perego et al., 2012; Leng et al., 2012; Larour et al., 2012; Aschwanden et al., 2012; Cornford et al., 2013; Gagliardini et al., 2013; Brinkerhoff and Johnson, 2013), ISM model “physics” (for example, the addition of improved models of basal sliding coupled to explicit subglacial hydrology, e.g., Schoof, 2005; Werder et al., 2013; Hewitt, 2013; Hoffman and Price, 2014; Bueler and van Pelt, 2015; and ice damage, fracture, and calving, e.g., Åström et al., 2014; Bassis and Ma, 2015; Borstad et al., 2016; Jiménez et al., 2017), and the coupling between ISMs and Earth system models (ESMs) (e.g., Ridley et al., 2005; Vizcaíno et al., 2008, 2009; Fyke et al., 2011; Lipscomb et al., 2013). These “nextgeneration” ISMs have been applied to communitywide experiments focused on assessing (i) the sensitivity of ISMs to idealized and realistic boundary conditions and environmental forcing and (ii) the potential future contributions of ice sheets to sea level rise (see, e.g., Pattyn et al., 2013; Nowicki et al., 2013a, b; Bindschadler et al., 2013; Shannon et al., 2013; Edwards et al., 2014b).
While these efforts represent significant steps forward, nextgeneration ISMs continue to confront new challenges. These come about as a result of (1) applying ISMs to larger (wholeicesheet), higherresolution (regionally 𝒪(1 km) or less), and more realistic problems, (2) adding new or improved submodels 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 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:

parallel, scalable, and robust linear and nonlinear solvers;

variable and/or adaptive mesh resolution;

computational kernels based on flexible programming models to allow for implementation on a range of highperformance computing (HPC) architectures; and ^{1}

automatic differentiation capability for the computation of adjoint sensitivities to be used in highdimensional parameter field optimization and uncertainty quantification.
Based on these considerations, we have developed a new land ice model, MALI, which is composed of three 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 (Ringler et al., 2013; Petersen et al., 2015, 2018), shallow water (Ringler et al., 2011), sea ice (Turner et al., 2018), and land ice have 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 the development of the physics and features relevant to their application.
The framework code includes shared modules for fundamental model operation. Significant capabilities include the following.

Description of model data types. MPAS uses a handful of fundamental Fortranderived types for basic model functionality. Model variables specific to an MPAS core are handled through custom groupings of model fields called pools, for which custom accessor routines exist. Corespecific variables are easily defined in XML syntax in a registry, and the framework parses the registry, defines variables, and allocates memory as needed.

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 a MessagePassing Interface (MPI). This 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 the transfer of data from disk to memory. Shared memory parallelization through OpenMP is also supported, but the implementation is left up to each MPAS core.

Parallel input and output capabilities. MPAS performs parallel input and output of data from and to disk through the 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 provides easy runtime configuration of what fields are to be written to what file name and at what frequency through an XML streams file. The MPAS Framework includes additional functionality specific to providing a flexible model restart capability.

Advanced timekeeping. MPAS uses a customized version of the timekeeping functionality of the Earth System Modeling Framework (ESMF), which includes a robust set of time and calendar tools used by many Earth system models (ESMs). This allows for the 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 multiscale physics and simplifies coupling to ESMs. To manage the complex date–time types that ensue, the MPAS Framework provides routines for the arithmetic of time intervals and the definition of alarm objects for handling events (e.g., when to write output, when the simulation should end).

Runtimeconfigurable control of model options. Model options are configured through namelist files that use standard Fortran namelist file format, and input/output is configured through streams files that use XML format. Both are completely adjustable at run time.

Online runtime analysis framework. See Sect. 6.2 for examples.
Additionally, a number of shared operators exist to perform common operations on model data. These include geometric 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.
2.1 Model meshes
The MPAS mesh specification is general enough to describe unstructured meshes on most twodimensional manifold spaces; however, most applications use centroidal Voronoi tessellations (Du and Gunzburger, 2002) on a sphere or plane. This paper 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 tessellation 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 highquality cells because cells tend towards equidimensional aspect ratios, and mesh resolution (where nonuniform) changes smoothly. On both planes and spheres, Voronoi tessellations tend toward perfect hexagons as resolution is increased. Note that while the 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 twodimensional 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 unstructured mesh specification and, for most operations, periodic cell relationships are handled the same as for neighboring cell relationships. MPAS meshes are static in time. The vertical coordinate, if needed by an MPAS core, is extruded from the base horizontal mesh. Each MPAS core chooses its own vertical coordinate system. A comprehensive suite of tools for the generation of centroidal Voronoi tessellations on a plane or sphere has been developed, as have tools for modifying existing meshes (e.g., removing unneeded cells, coordinate transformations, etc.) and converting some common unstructured mesh formats (e.g., Triangle; Shewchuk, 1996) to the MPAS specification.
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 end points 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 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 provides 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 been done with such meshes to date. MALI employs a Cgrid 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 (Fig. 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 variableresolution centroidal Voronoi meshes. Additionally, the JIGSAW(GEO) meshgeneration tool (Engwirda, 2017a, b) can be used to efficiently generate highquality, variableresolution meshes with databased density functions. Density functions that are a function of observed ice velocity or its spatial derivatives and/or distance to the existing or potential future grounding line position have been used.
Albany is an open source, C++ multiphysics code base for the solution and analysis of coupled systems of partial differential equations (PDEs) (Salinger et al., 2016). It is a finiteelement 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 templatebased 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 of the discrete PDE residual, as well as forward and adjoint sensitivities. Albany can solve largescale PDEconstrained optimization problems using the Trilinos optimization package 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 these diagnostic momentum balance equations as AlbanyLI. AlbanyLI was formerly known as Albany/FELIX (Finite Elements for Land Ice eXperiments) and is 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).
The “dynamical core” of the MALI ice sheet model solves the governing equations expressing the conservation of momentum, mass, and energy.
4.1 Conservation of momentum
Treating glacier ice as an incompressible fluid in a lowReynoldsnumber 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 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}.
Deformation results from the deviatoric stress, τ_{ij}, which relates to the full stress tensor as
for which $\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\sigma}}_{kk}$ is the mean compressive stress and δ_{ij} is the Kronecker delta (or the identity tensor). Stress and strain rate are related through the constitutive relation,
where $\dot{{\mathit{\u03f5}}_{ij}}$ is the strain rate tensor and μ_{e} is the “effective” nonNewtonian ice viscosity given by Nye's generalization of Glen's flow law (Glen, 1955):
In Eq. (5), A is a temperaturedependent rate factor, n is an exponent commonly taken as 3 for polycrystalline glacier ice, and γ is an ice “stiffness” factor (inverse enhancement factor related to the commonly used enhancement factor E_{f} by $\mathit{\gamma}={E}_{\text{f}}^{\frac{\mathrm{1}}{n}}$) used to account for other impacts on ice rheology, such as impurities or crystal anisotropy (see also Sect. 6.1). The effective strain rate ${\dot{\mathit{\u03f5}}}_{\text{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}:
Finally, the rate factor A follows an Arrhenius relationship,
in which A_{o} is a constant, T^{*} is the temperature (relative to the pressure melting point), 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 reducedorder approximation applied and are discussed further below.
4.2 Reducedorder equations
Ice sheet models solve Eqs. (2)–(8) with varying degrees of complexity in terms of the tensor components in Eqs. (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 – reducedorder approximations of the full momentum balance are 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 firstorderaccurate “Blatter–Pattyn” approximation and a zeroorder “shallow ice approximation” as described in more detail in the following sections.
4.2.1 Firstorder velocity solver and coupling
Ice sheets typically have a small aspect ratio and small surface and bed slopes. These characteristics imply that reducedorder approximations of the Stokes momentum balance may apply over large areas of the ice sheets, potentially allowing for significant computational savings. Formal derivations involve nondimensionalizing the Stokes momentum balance and introducing a geometric scaling factor, $\mathit{\delta}=H/L$, where H and L represent characteristic vertical and horizontal length scales (often taken as the ice thickness and the ice sheet span), respectively. Upon conducting an asymptotic expansion, reducedorder 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 firstorderaccurate Stokes approximation is arrived at by retaining terms of 𝒪(δ^{1}) and lower (the reader is referred to Schoof and Hindmarsh, 2010, and Dukowicz et al., 2010, for additional discussion^{3}).
Using the notation of Perego et al. (2012) and Tezaur et al. (2015a) ^{4}, the firstorderaccurate 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\equiv s(x,y)$ represents the ice sheet upper surface, and the vectors ${\dot{\mathit{\u03f5}}}_{\mathrm{1}}$ and ${\dot{\mathit{\u03f5}}}_{\mathrm{2}}$ are given by
and
Akin to Eqs. (5) and (6), μ_{e} in Eq. (9) represents the effective viscosity but for the case of the firstorder stress balance with an effective strain rate given by
rather than by Eq. (6), and with individual strain rate terms given by
At the upper surface, a stressfree boundary condition is applied,
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 Sect. 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 AlbanyLI momentum balance solver, which is built using the Albany and Trilinos software libraries discussed above. The mathematical formulation, discretization, solution methods, verification, and scaling of AlbanyLI are discussed in detail in Tezaur et al. (2015a). AlbanyLI implements a classic finiteelement discretization of the firstorder 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 a higherorder Gauss quadrature rule on elements containing the grounding line, which corresponds to the subelement parameterization SEP3 proposed in Seroussi et al. (2014). Additional exploration of solver scalability and demonstrations of solver robustness on largescale, highresolution, realistic problems are discussed in Tezaur et al. (2015b). The efficiency and robustness of the nonlinear solvers are achieved using a combination of the Newton method (damped with a line search strategy when needed) and a parameter continuation algorithm for the numerical regularization of the viscosity. The scalability of the 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 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 illconditioned linear systems.
The AlbanyLI firstorder velocity solver written in C++ is coupled to MPAS written in Fortran using an interface layer. Albany uses a threedimensional 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 part of the Delaunay triangulation, dual to an MPAS Voronoi mesh, that contains active ice and is generated by the interface. The bed topography, ice lower surface, ice thickness, basal friction coefficient (β), and threedimensional 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 returns the x and y components of velocity at each cell center and layer interface, the normal component of velocity at each cell edge and layer interface, and viscous dissipation at each cell vertex and layer midpoint.
The interface code defines the lateral boundary conditions on the finiteelement 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 the calculation of normal velocity on edges (u_{n}) required for the advection of ice out of the final cell containing dynamic ice (Fig. 2). The interface identifies three types of lateral boundaries for the firstorder velocity solve: terrestrial, floating marine, and grounded 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 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.
4.2.2 Shallow ice approximation velocity solver
A similar procedure to that described above for the firstorderaccurate 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., 2000), in this case by retaining only terms of 𝒪(δ^{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.
SIA ice sheet models typically combine the momentum and mass balance equations to evolve the ice geometry directly in the form of a depthintegrated, twodimensional 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 firstorder 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. (16) is not a concern.
We implement Eq. (16) in sigma coordinates on cell edges for which we only require the normal component of velocity, u_{n}:
where x_{n} is the normal direction to a given edge and ${\mathit{u}}_{{b}_{\text{n}}}$ is sliding velocity in the normal direction to the edge. We average A and H from cell centers to cell edges. $\frac{\text{d}s}{\text{d}{x}_{\text{n}}}$ 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.
4.3 Conservation of mass
Ice sheet mass transport and evolution is conducted using the principle of conservation of mass. Assuming constant density to write the 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, $\stackrel{\mathrm{\u203e}}{\mathit{u}}$ is depthaveraged velocity, $\dot{a}$ is surface mass balance, and $\dot{b}$ is basal mass balance. Both $\dot{a}$ and $\dot{b}$ are positive for ablation and negative for accumulation.
Equation (18) is used to update thickness in each grid cell on each time step using a forward Euler, fully explicit time evolution scheme. Eq. (18) is implemented using a finitevolume method such that fluxes are calculated for each edge of each cell to calculate $\mathrm{\nabla}\cdot H\stackrel{\mathrm{\u203e}}{\mathit{u}}$. Specifically, we use a firstorder upwind method that applies the normal velocity on each edge (u_{n}) and an upwind value of cellcentered ice thickness. Note that with the Blatter–Pattyn velocity solver, normal velocity is interpolated from cell centers to edges using the finiteelement basis functions in Albany. In the shallow ice approximation velocity solver, normal velocity is calculated natively at edges. The MPAS Framework includes a higherorder fluxcorrected transport scheme (Ringler et al., 2013) for which we have performed some initial testing, but is not routinely used in MALI at this time.
Tracers are advected horizontally layer by layer with a similar equation:
where Q_{t} is a tracer quantity (e.g., temperature; see below), l is layer thickness, and $\dot{S}$ 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.
Vertical advection of tracers is included through a vertical remapping operation. On the upper and lower domain boundaries, the grid moves to follow the material, and in the interior we maintain fixed layer fractions that need to be updated on each time step after Eqs. (18) and (19) are applied. The model does not explicitly calculate vertical velocity, but the appropriate vertical transport of tracers occurs during this vertical remapping operation. We employ a firstorder vertical remapping method. 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.
4.4 Conservation of energy
Conservation of energy within glaciers can be formulated in terms of temperature or enthalpy (internal energy) (Aschwanden et al., 2012; Kleiner et al., 2015). The enthalpy formulation has the advantage of eliminating the need for tracking the cold–temperate transition surface, as both cold (below the pressure melting point) and temperate (at the pressure melting point) ice regions are handled with the same equations. MALI includes both temperature and enthalpy formulations. In both cases, an operator splitting technique is used. At each time step, an implicit vertical solve accounting for the diffusion and dissipation terms (described below) is performed, followed by explicit advection of the resulting temperature or enthalpy field (described above in Sect. 4.3). We describe the temperature formulation in detail, followed by a briefer description of the enthalpy formulation that uses a similar procedure. Note that the thermal model described here shares a common lineage with that of the Community Ice Sheet Model, and parts of the description below are therefore similar to the documentation of the thermal solver in the Community Ice Sheet Model (Price et al., 2015; Lipscomb et al., 2018).
4.4.1 Temperature formulation
Conservation of energy can be expressed in terms of temperature through the threedimensional, advective–diffusive heat equation:
with thermal conductivity k and heat capacity c. In Eq. (20), the rate of temperature change (lefthand side) is balanced by diffusive, advective, and internal (viscous dissipation; see Eq. 25 for Φ) source terms (first, second, and third terms on the righthand side, respectively). In MALI we solve an approximation of Eq. (20),
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.
Temperatures are staggered in the vertical relative to velocities and are located at the centers of n_{z−1} vertical layers, which are bounded by n_{z} 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 n_{z−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 n_{z+1} temperature values within each column.
As mentioned above, Eq. (21) is solved by first performing an implicit vertical solve accounting for the diffusion and dissipation terms (described below), followed by explicit advection of the resulting temperature field. The method for evolving ice temperature and default parameter value choices are adapted from the implementation in the Community Ice Sheet Model (Price et al., 2015; Lipscomb et al., 2018), which is in turn based on the Glimmer model (Rutt et al., 2009). The choice of constant k with a temperate ice value (Table A1) will lead to underestimation of conduction in cold ice. Relaxation of this assumption is planned for future releases of MALI.
Vertical diffusion
Using a “sigma” vertical coordinate, the vertical diffusion portion of Eq. (21) 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 ${\stackrel{\mathrm{\u0303}}{\mathit{\sigma}}}_{k}$ is the value of σ at the midpoint of layer k, halfway between σ_{k} and σ_{k+1}. The second partial derivative, defined at the midpoint of layer k, is then given by
By inserting Eq. (23) into Eq. (24), we obtain the discrete form of the vertical diffusion term in Eq. (21):
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 Eq. (21) is given by the product of the stress and strain rate tensors:
The change to deviatoric stress on the righthand side of Eq. (25) follows from terms related to the mean compressive stress (or pressure) dropping out due to incompressibility. Analogous to the effective strain rate given in Eq. (6), the effective deviatoric stress is given by
which can be combined with Eqs. (25) and (6) to derive an expression for the viscous dissipation in terms of effective deviatoric stress and strain:
Finally, an analog to Eq. (4) gives
which can be used to eliminate ${\dot{\mathit{\u03f5}}}_{\text{e}}$ in Eq. (27) and arrive at an alternate expression for the dissipation based on only two scalar quantities:
The viscous dissipation source term is computed within AlbanyLI at MPAS cell vertices and then reconstructed at cell centers in MPAS.
For the SIA model, dissipation can be calculated in sigma coordinates as
which can be combined with Eq. (16) to make
We calculate Φ on cell edges following the procedure described for Eq. (17) and then interpolate Φ back to cell centers to solve Eq. (21).
Vertical temperature solution
The vertical diffusion portion of Eq. (21) is discretized according to
where a_{k} and b_{k} are defined in Eq. (25), n is the current time level, and n+1 is the new time level. Because the vertical diffusion 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}_{\mathrm{air}},\mathrm{0})$, where the mean annual surface air temperature T_{air} is a twodimensional 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 bedparallel vectors of basal shear stress and basal velocity, respectively, and the friction law from Eq. (14) becomes
If the basal temperature ${T}_{{n}_{z}}<{T}_{\mathrm{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, on the other hand, ${T}_{{n}_{z}}={T}_{\mathrm{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 Sect. 5.1), or may simply be ignored. If basal water is present locally, ${T}_{{n}_{z}}$ is held at T_{pmp}.
For floating ice the basal boundary condition is simpler: ${T}_{{n}_{z}}$ 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 the melting of ice internally, and the resulting melt is assumed to drain to the bed immediately.
If Eq. (36) applies, we compute M_{b} and adjust the basal water depth. When the basal water goes to zero, ${T}_{{n}_{z}}$ is set to the temperature of the lowest layer (less than T_{pmp} at the bed) and flux boundary conditions apply during the next time step.
Temperature advection
Temperature advection in any individual layer k is treated using tracer advection, as in Eq. (19) 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 an altered vertical layer spacing with respect to σ coordinates, a vertical remapping scheme is applied to provide the necessary vertical advection of temperature. 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.
Enthalpy formulation
The specific enthalpy (internal energy), E, in ice sheets and glaciers can be expressed as a combination of ice temperature (T) and liquid water fraction (water content; ω) (Aschwanden et al., 2012):
where T_{ref} is the reference temperature, c is the heat capacity of ice, L is the latent heat of fusion, and E_{pmp} is the specific enthalpy at the pressure melting point for different vertical locations (T_{pmp}(z)) defined as
The balance equation for enthalpy reads
where K is the diffusivity of ice defined differently in cold and temperate ice:
where ν is the water diffusivity in temperate ice, which is generally taken as an empirical small number due to a lack of knowledge (Greve and Blatter, 2016).
The implementation of the enthalpy model follows that of the temperature model. Vertical diffusion is as described above but replacing T with E. Viscous dissipation remains unchanged. Boundary conditions in the vertical enthalpy solution follow those applied for the temperature formulation above, but cast in terms of E. Advection is as described above for temperature but replacing T with E. Verification of the enthalpy model is described below in Sect. 7.3.
Additional physical processes currently implemented in MALI are a massconserving subglacial hydrology model and a small number of basic schemes for iceberg calving. These are described in more detail below.
5.1 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 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 massconserving model of subglacial hydrology that includes representations of water storage in till, distributed drainage, and channelized drainage and is coupled to ice dynamics. The model is based on the model of Bueler and van Pelt (2015) but modified for MALI's unstructured horizontal grid and with an additional component for channelized drainage. 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.
5.1.1 Till
The simple till component represents local storage of water in subglacial till without horizontal transport within the till. The evolution of the effective water depth in till, W_{till}, is therefore a balance of the delivery of meltwater, m_{b}, to the till, the 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 excess of the maximum storage capacity of the till, ${W}_{\text{till}}^{max}$, is instantaneously transferred as a source term to the distributed drainage system through the γ_{t} term.
5.1.2 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, 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 areaaveraged cavity space is a balance of the opening of cavity space by the glacier sliding over bedrock bumps and closing through creep of the ice above. The model uses the common assumption (e.g., Schoof, 2010; Hewitt, 2011; Werder et al., 2013; Hoffman and Price, 2014) that cavities always remain water filled (cf. Schoof et al., 2012), so cavity space can be represented by the effective water depth in the macroporous sheet, W:
where c_{s} is bed roughness parameter, W_{r} is the maximum bed bump height, c_{cd} is a creep scaling parameter representing geometric and possibly other effects, and A_{b} is the ice flow parameter of the basal ice.
Water flow in the distributed drainage system, q, is driven by the hydropotential gradient and is described by a general power law:
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.
5.1.3 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 the 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 distributed 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 the dissipation of potential energy, Ξ, and sensible heat change of water, Π, due to changes in the pressuredependent 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 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 the sensible heat change of the water, Π, is necessitated by the assumption that the water always remains at the pressuredependent 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 pressuredependent 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.
5.1.4 Drainage component coupling
Equations (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 function applied along the locations of the linear channels.
Combining Eqs. (50) and (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 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 run time. The most common configuration currently used is to run with distributed drainage only.
5.1.5 Numerical implementation
The drainage system model is implemented using finitevolume 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. Equation (50) is evaluated by summing tendencies from discrete fluxes into or out of each cell. Firstorder upwinding is used for advection. At landterminating ice sheet boundaries, P_{w}=0 is applied as the boundary condition. At marineterminating ice sheet boundaries, the boundary condition is ${P}_{\text{w}}={\mathit{\rho}}_{\text{w}}g{z}_{\text{b}}$, where ρ_{w} is ocean water density. The drainage model uses explicit forward Euler time stepping using Eqs. (41), (50), (46), and (50). 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.
We acknowledge that the noncontinuum 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 that the dominant channel location is controlled by topography, mitigating this issue.
5.1.6 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 a modified Weertmanstyle power law (Bindschadler, 1983; Hewitt, 2013) that adds a term for effective pressure to Eq. (14):
where C_{0} is a friction parameter. Implementations 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 for the negative feedback described by Hoffman and Price (2014) in which 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.
5.1.7 Verification and realworld application
To verify the implementation of the distributed drainage model, we use the nearly exact solution described by Bueler and van Pelt (2015). The problem configuration uses distributed drainage only on a twodimensional, 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 nearexact solution and then run the model forward for 1 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 firstorder convergence (Fig. 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}. Steadystate solutions of the drainage system effective pressure, water fluxes, and channel development for an idealized ice sheet with varying magnitudes of meltwater input (SHMIP experiment suites A and B) compared between MALI and other models of similar complexity (GlaDS, Elmer) are very similar.
To demonstrate a realworld application of the subglacial hydrology model, we perform a standalone subglacial hydrology simulation of the entire Antarctic ice sheet on a uniform 20 km resolution mesh (Fig. 4). We force this simulation with basal sliding and basal melt rate after optimizing the firstorder velocity solver to surface velocity observations (Fig. 4a). We then run the subglacial hydrology model to steady state with only distributed drainage active and using standard parameter values and prescribed ice dynamic forcing. Though this mesh is too coarse to provide scientifically valid results, the modeled subglacial hydrologic state is reasonable. For example, the subglacial water flux increases downglacier and is greatest in fastflowing outlet glaciers and ice streams (Fig. 4b), as expected from theory and seen in other subglacial hydrology models (e.g., Le Brocq et al., 2009). Calibrating parameters for the subglacial hydrology model and a basal friction law and performing coupled subglacialhydrology–icedynamics simulations are beyond the scope of this paper; we merely mean to demonstrate plausible behavior from the subglacial hydrology model for a realistic icesheetscale problem.
5.2 Iceberg calving
MALI includes a few simple methods for removing ice from calving fronts during each model time step.

All floating ice is removed.

All floating ice in cells with an ocean bathymetry deeper than a specified threshold is removed.

All floating ice thinner than a specified threshold is removed.

The calving front is maintained at its initial location by adding or removing ice after thickness evolution is complete. When ice is completely lost in a grid cell through evolution, it is replaced with a thin layer of ice (default value of 1 m). This does not conserve mass or energy but provides a simple way to maintain a realistic ice shelf extent (e.g., for model spinup).

Eigencalving scheme (Levermann et al., 2012). Calving front retreat rate, C_{v}, is proportional to the product of the principal strain rates ($\dot{{\mathit{\u03f5}}_{\mathrm{1}}},\dot{{\mathit{\u03f5}}_{\mathrm{2}}}$) if they are both extensional:
$$\begin{array}{}\text{(57)}& {C}_{\text{v}}={K}_{\mathrm{2}}\dot{{\mathit{\u03f5}}_{\mathrm{1}}}\dot{{\mathit{\u03f5}}_{\mathrm{2}}}\phantom{\rule{0.25em}{0ex}}\text{for}\phantom{\rule{0.25em}{0ex}}\dot{{\mathit{\u03f5}}_{\mathrm{1}}}>\mathrm{0}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}\dot{{\mathit{\u03f5}}_{\mathrm{2}}}>\mathrm{0}.\end{array}$$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 the formation of 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 subgrid 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 at subgrid scales (Albrecht et al., 2011; Bondzio et al., 2016). In Sect. 8 below, we demonstrate the eigencalving scheme applied to a realistic Antarctic ice sheet simulation. More sophisticated calving schemes are currently under development.
To demonstrate a realworld application of the eigencalving parameterization, we perform a 1000year spinup of Antarctica with evolving velocity, geometry, and temperature and active eigencalving (Fig. 5). For the purposes of this demonstration, we use the same uniform, 20 km resolution mesh used in Fig. 4, which is too coarse to accurately resolve grounding line dynamics (see Sect. 7.5) and therefore should not be interpreted as a scientifically realistic simulation. We use an initial internal ice temperature field from Van Liefferinge and Pattyn (2013) and the optimization capability described below in Sect. 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 (Fig. 5a). For the 1000year spinup, we apply steady forcing of presentday estimates for surface mass balance and submarine melting (Lenaerts et al., 2012, and Rignot et al., 2013, respectively). For temperature boundary conditions, we apply the steady geothermal flux field from Shapiro and Ritzwoller (2004) and the surface (2 m) air temperature field from Lenaerts et al. (2012). We apply eigencalving calving with the K_{2} parameter tuned individually for large ice shelves and a minimum calving front thickness threshold of 100 m. This spinup, albeit much too short to come to full equilibrium, allows ample time for migration of the calving front and grounding line and removes a substantial portion of the largest model transients. It demonstrates that a tuned eigencalving parameterization is capable of maintaining stable and realistic calving front positions in MALI during ice sheet evolution (Fig. 5a, b), consistent with its implementation in other models (Levermann et al., 2012; Feldmann and Levermann, 2015).
6.1 Optimization
MALI includes an optimization capability through its coupling to the AlbanyLI momentum balance solver described in Sect. 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 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., 2015) and we focus here primarily on optimizing the model velocity field relative to observed surface velocities. Briefly, we consider the optimization functional
where the first term on the righthand 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 Eq. 5), and the third and fourth terms on the RHS are Tikhonov regularization terms given by
σ_{u} is an estimate for 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 overfitting the observations). The optimal values of α_{β} and α_{γ} can be chosen through a standard Lcurve analysis. The optimization problem is solved using the limitedmemory BFGS method, as implemented in the Trilinos package ROL^{6}, on the reducedspace problem. The functional gradient is computed using the adjoint method.
An example application of the optimization capability applied to a realistic, wholeicesheet problem is given below in Sect. 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 Eq. (52). While not applied here, we are able to introduce additional observationalbased constraints (e.g., mass balance terms) and optimize additional model variables (e.g., the ice thickness). These are necessary, for example, when targeting model initial conditions that are in quasiequilibrium with some applied climate forcing. These capabilities are discussed in more detail in Perego et al. (2014).
6.2 Simulation analysis
As with other climate model components built using the MPAS Framework, MALI supports the development and application of “analysis members”, which allow for a wide range of runtimegenerated simulation diagnostics and statistics output at userspecified time intervals. Support tools included with the code release allow for the definition of any number or combination of predefined “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^{7}. Pythonbased scripts are available for editing GeoJSON feature files, combining or splitting 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 groundpenetrating radar profile lines). In Sect. 8 below we demonstrate the analysis capability applied to an idealized simulation of the Antarctica ice sheet.
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.
MALI currently includes 86 automated system regression tests that run the model for various problems with analytic solutions or community benchmarks. In addition to checking the accuracy of model answers, some of the tests check that model restarts give bitforbit exact answers with longer runs without restart. Some others check that the model gives bitforbit exact answers on different numbers of processors. All but 20 of the longerrunning tests are run every time new features are added to the code, and these tests each also include a check for answer changes. The verification and benchmark descriptions below are the most important examples from the larger test suite.
7.1 Halfar analytic solution
In Halfar (1981, 1983), Halfar described an analytic solution for the timeevolving geometry of a radially symmetric, isothermal dome of ice on a flat bed with no accumulation flowing under the shallow ice approximation. This provides an obvious test of the implementation of the shallow ice velocity calculation and thickness evolution schemes in numerical ice sheet models and a way to assess model order of convergence (Bueler et al., 2005; Egholm and Nielsen, 2010). Bueler et al. (2005) showed that the Halfar test is the zeroaccumulation 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 21 213.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 that the root mean square error in model thickness decreases as model grid spacing is decreased (Fig. 6a). The order of convergence, 0.78, is somewhat lower than expected from the firstorder methods used for advection and time evolution.
We also use this test to assess the accuracy of simulations with variable resolution. We perform an additional run of the Halfar test using a variableresolution mesh, which has 1000 m cell spacing beyond a radius of 20 km that transitions to 5000 m cell spacing at a radius of 3 km (Fig. 6b), generated with the JIGSAW(GEO) meshgeneration tool (Engwirda, 2017a, b). The root mean square error in thickness for this simulation is similar to that for the uniform 1000 m resolution case (Fig. 6a), providing confidence in the advection scheme applied to variableresolution meshes. The variableresolution mesh has about half the cells of the 1000 m uniformresolution mesh.
7.2 EISMINT
The 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 thermomechanical shallow ice models with artificially constructed exact solutions. While their approach has the notable advantage of providing exact solutions, we have not implemented the nonphysical threedimensional 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, 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 kyr. 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 MALI is 2976.7 m compared with a mean of 2978.0±19.3 m for the 10 threedimensional 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 $\mathrm{13.34}\pm \mathrm{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) grow an ice sheet to steady state over 200 kyr from an initial condition of no ice, but with different air temperature boundary conditions. Additional experiments use the steadystate 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 10 different models contributed 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 10 different models, which we have summarized here with the corresponding values calculated by MALI (Table 3). All MALI results fall within the range of 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.
A longstudied 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 (Fig. 7). Bueler et al. (2007) argue that 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 cold spokes. However, that model differs from previously analyzed models in several ways: it solves a threedimensional, advective–diffusive description of an enthalpy formulation for energy conservation; it uses the finiteelement method on unstructured meshes; and 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.
7.3 Enthalpy benchmarks
Kleiner et al. (2015) present a set of benchmark experiments to test numerical models of ice sheet enthalpy evolution. The experiments use designs that allow for comparison to analytic solutions. Here we only give very brief descriptions of the enthalpy benchmark experiments. Details can be found in Kleiner et al. (2015). The benchmark includes two different experiments for testing the capability of transient behavior and horizontal and vertical advection of the enthalpy model. Both experiments use a parallelsided ice slab with constant thickness and inclination and prescribed ice dynamics decoupled from thermodynamics, resulting in effectively onedimensional vertical experiments.
7.3.1 Experiment A
Both heat advection and frictional heating are neglected in this experiment. Heat diffusion is the only controlling process in the redistribution of enthalpy; i.e., the enthalpy balance equation simplifies to
To test the transient ability of the enthalpy model, this experiment was designed to run for a time period of 300 kyr. During the run, the geothermal heat flux (G) is constant over time, but the surface temperature (upper Dirichlet boundary condition) changes in three different time intervals.
(${T}_{\text{s}}=\mathrm{5}$ ^{∘}C, from personal communication with Thomas Kleiner, compared to ${T}_{\text{s}}=\mathrm{10}$ ^{∘}C, incorrectly stated in Kleiner et al., 2015.)
During the time period of 100–150 ka, when the surface temperature rises to −5 ^{∘}C, the glacier base becomes temperate and the basal boundary condition changes from Neumann type ($\partial T/\partial z=G$) to Dirichlet type (T=T_{pmp}). Then the surface temperature switches back to its initial value, −30 ^{∘}C, for testing the reversibility of the enthalpy model. In this experiment the basal water content produced by basal melting is allowed to freely accumulate in order to test the basal melt rate calculation.
From Fig. 8a, we can clearly see that the basal temperature becomes steady (−10 ^{∘}C) at around 50 ka, and then starts to rise to the pressure melting point after the prescribed change in surface temperature to −5 ^{∘}C at 100 ka. The model then keeps the glacier base temperate until 225 ka, with the prescribed change in the surface temperature boundary condition back to −30 ^{∘}C at 150 ka. The subglacial water layer starts to accumulate from around 110 ka when the basal ice starts to melt (Fig. 8b) and reaches a maximum layer thickness of around 133.9 m at around 155 ka (Fig. 8c). After the surface gets colder again, it gradually decreases and disappears completely at 225 ka. From the comparison with the analytical basal melt rate result during 150–170 ka (green line), we can clearly see that the enthalpy model of MALI captures the features of basal melting (water content production).
7.3.2 Experiment B
In experiment B, a 200 m thick, 4^{∘} downwardinclined slab is used as the model domain. A particular objective of this experiment is to test the model ability to find the correct position of the cold–temperate ice transition interface. The horizontal velocity is given as an analytical (shallow ice approximation type) expression, and the vertical velocity is set to be constant (i.e., thermomechanically decoupled). In addition, the geothermal heat flux is set to be zero during the model run so that the englacial strain heating is the only energy source in the enthalpy balance equation. Initialized from an isothermal field (−1.5 ^{∘}C), our model spins up for several thousand years with a constant surface temperature of −3 ^{∘}C until a steadystate temperate ice layer thickness is achieved (Fig. 9).
From Fig. 9 we can see that our enthalpy model can predict very close enthalpy, temperature, and basal water content results (blue lines) compared to analytical solutions (green lines; almost overlapped). Using a uniform vertical resolution of 1 m, MALI simulates a temperate ice layer thickness of 19 m, nearly identical to the analytical output. This experiment also shows that MALI can accurately compute the englacial enthalpy distribution in the presence of nonzero horizontal ice advection.
7.4 ISMIPHOM
The Ice Sheet Model Intercomparison ProjectHigher Order Models (ISMIPHOM) is a set of community benchmark experiments for testing higherorder approximations of ice dynamics (Pattyn et al., 2008). Tezaur et al. (2015a) describe results from the AlbanyLI velocity solver for ISMIPHOM experiments A (flow over a bumpy bed) and C (ice stream flow). For all configurations of both tests, AlbanyLI results were within 1 standard deviation of the mean of firstorder models presented in Pattyn et al. (2008) and showed excellent agreement with the similar firstorder model formulation of Perego et al. (2012). These tests only require a single diagnostic solve of velocity, and thus results through MALI match those of the standalone AlbanyLI code it is using.
7.5 MISMIP3d
The Marine Ice Sheet Model Intercomparison Project3d (MISMIP3d) is a community benchmark experiment testing the grounding line migration of marine ice sheets and includes nontrivial effects in all three dimensions (Pattyn et al., 2013). The experiments 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 steadystate 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 the domain. From this steady state, the P75R perturbation experiment reduces basal friction by a maximum of 75 % across a Gaussian ellipse centered where the Stnd grounding line position crosses the symmetry axis. The perturbation is applied for 100 years, resulting in a curved grounding line that is advanced along the symmetry axis. After the completion of P75S, a reversibility experiment named P75R removes the basal friction perturbation and allows the ice sheet to relax back towards the Stnd state.
Pattyn et al. (2013) report results from 33 models of varying complexity applied at resolutions ranging from 0.1 to 20 km. Participating models used depthintegrated shallow shelf or L1L1/L2L2 approximations, hybrid shallow ice–shallow shelf approximation, or the complete Stokes equations; there were no threedimensional firstorder 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 the reversibility of the steadystate grounding line position after experiments P75S and P75R. Reversibility required a grid resolution well below 1 km without a subgrid parameterization of grounding line position and grids a couple of times coarser with a grounding line parameterization (Pattyn et al., 2013; Gladstone et al., 2010). The steadystate grounding line position in the Stnd experiment was dependent on the stress approximation employed, with the Stokes model calculating grounding lines the farthest upstream and models that simplify or eliminate vertical shearing (e.g., shallow shelf) 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 (Pattyn et al., 2013).
We find that MALI using the AlbanyLI Blatter–Pattyn 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 convergence by comparing the position of the steadystate grounding line in the Stnd experiment for a range of resolutions against our highestresolution configuration of 250 m (Fig. 10). With the grounding line parameterization, the grounding line positions at 500 and 250 m resolution are very similar (differing by less than the grid resolution), whereas without the grounding line parameterization the grounding line positions in our two highestresolution simulations still differ by 6 km. The converged grounding line position for the Stnd experiment with MALI is 533 km. The grounding line position from our threedimensional firstorder stress approximation model falls between that of the L1L2 and Stokes models reported by Pattyn et al. (2013), consistent with the intermediate level of approximation of our model. The dissertation work by Leguy (2015) reported similar results for the Blatter–Pattyn velocity solver in the Community Ice Sheet Model (Lipscomb et al., 2018) when using a subgrid grounding line parameterization.
Reversibility of the P75S and P75R experiments shows the same grid resolution requirement of 500 m, while the 1 km simulation with grounding line parameterization is close to reversible (Fig. 11). Our highestresolution 250 m simulation without grounding line parameterization does show reversibility at the end of P75R (not shown), but the results differ somewhat from the runs with grounding line parameterization due to the differing starting position determined from the Stnd experiment. 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 threedimensional firstorder stress balance look most similar to those of the “SCO6” 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 of decades (Pattyn et al., 2013).
In addition to MISMIP3d, we have used MALI to perform the MISMIP+ experiments (AsayDavis et al., 2016). These results are included in the MISMIP+ results paper in preparation and not shown here.
To demonstrate a largescale, semirealistic ice sheet simulation that exercises many of the model capabilities discussed above, we describe MALI results from the initMIPAntarctica experiments. 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 spinup, optimization approaches, or both – three 100year forward model experiments are conducted in order to examine model response to (1) an unforced control run, (2) an idealized surface mass balance perturbation, and (3) an idealized subiceshelf melt perturbation. Here we show results from simulations 1 and 3 for brevity. Additional details of the experiments are described in http://www.climatecryosphere.org/wiki/index.php?title=InitMIPAntarctica (last access: 10 September 2018), 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 et al. (2018). Additional realistic applications of earlier versions of MALI to Greenland simulations are discussed in Shannon et al. (2013) and Edwards et al. (2014b).
The Antarctica model configuration we use here has 2 km resolution near grounding lines and in regions of marine (below sea level) bedrock in West Antarctica and regions of East Antarctica where presentday ice thickness is less than 2500 m to ensure that the grounding line remains in the fineresolution region even under full retreat of West Antarctica and large parts of East Antarctica. The resolution then slowly coarsens to 20 km elsewhere, but maintaining no greater than 6 km resolution on ice shelves (Fig. 12a). The 2 km resolution in regions of potential grounding line migration was chosen as a balance between acceptable accuracy (Sect. 7.5) and computational cost. The mesh has 1 642 490 horizontal grid cells and uses 10 vertical layers, which are finest near the bed (4 % of total thickness) and coarsen towards the surface (23 % of total thickness).
Basal friction (β) and ice stiffness factor (γ) fields are optimized as described above in Sect. 6.1 to best allow the model to match observed surface velocities from Rignot et al. (2011). Calving position is fixed as described in option 2 in Sect. 5.2. To avoid energy conservation issues related to maintaining a fixed calving front position, we use a steadyinternal ice temperature field from Van Liefferinge and Pattyn (2013). Over the short timescales investigated here (200 years), we expect the assumption of fixed ice temperature to be a minor uncertainty (Seroussi et al., 2013). From this initial state, we perform a 99year relaxation with evolving velocity and geometry (Fig. 12), applying steady forcing of presentday estimates for surface mass balance and submarine melting (Lenaerts et al., 2012, and Rignot et al., 2013, respectively). This relaxation removes a substantial portion of the largest model transients (Fig. 12c, d). The model state at this point serves as our initial condition from which we run the control and subiceshelf melt perturbation experiments mentioned above. In this initial state, the volume above floatation mass loss for the entire Antarctic ice sheet is 602 Gt yr^{−1}. This is substantially larger than the current best estimate for Antarctic ice sheet mass loss of 109 Gt yr^{−1} for the 1992–2017 period, as well as the larger value of 219 Gt yr^{−1} for the 2012–2017 period (IMBIE team, 2018), and results largely from retreat and thinning in the Thwaites Glacier basin during the 99 years of relaxation (Fig. 12d).
The control simulation lost mass above floatation equivalent to a 167 mm sea level rise after 100 years. The subiceshelf melt perturbation experiment yielded the equivalent of a 250 mm sea level rise, an 83 mm sea level rise addition beyond the control run. During the control run, the Ross and Filchner–Ronne ice shelves experienced a modest slowdown with some acceleration near the grounding line (Fig. 13a). The Amery and Thwaites ice shelves exhibited significant speedup, with the acceleration propagating inland from Thwaites. Thickness changes in the control run are modest with pronounced thinning occurring only on Thwaites Glacier and inland of Cook Ice Shelf (Fig. 13b). The only noticeable grounding line changes in the control run are a slight retreat at Thwaites Glacier and a slight advance at Princess Ragnhild Coast in Queen Maud Land. In the subiceshelf melt perturbation experiment there is much more pronounced speedup at Thwaites and Pine Island glaciers, propagating far inland (Fig. 13c), with corresponding ice thinning up to 1000 m (Fig. 13d). All other ice shelves show increased thinning or reduced thickening relative to the control run, consistent with the application of increased ice shelf basal melt rates. The regions of the Ross and Filchner–Ronne ice shelves near the grounding line experience greater ice speedup than in the control run and some associated grounding line retreat, but little additional ice thinning. The Totten Glacier region exhibits significant ice speedup and thinning in the subiceshelf melt perturbation experiment.
To summarize the disparate regional behavior seen in the experiments, simulation statistics for selected drainage basins for the control and subiceshelf melt perturbation experiments are shown in Fig. 14 using the analysis capability described in Sect. 6.2. The applied basal melt forcing for selected basins can be seen in Fig. 14a. The prescribed subiceshelf melt perturbation according to the initMIP protocol increases for 40 years and then remains constant. The continued growth in total ice shelf basal melt after year 40 in Fig. 14a reflects increasing ice shelf area as the grounding line retreats, while we have forced the ice shelf calving front to remain fixed. Ice shelf thinning from increased basal melt results in increased flux across the grounding line, with the largest increase in the Thwaites–Pine Island catchment where the ice shelf basal melt perturbation was largest (Fig. 14b). In turn, increased flux across the grounding line leads to ice sheet mass loss (Fig. 14c) and grounding line retreat (Fig. 14d). The largest changes occur in the Thwaites–Pine Island catchment. Mass loss and grounding line retreat occur in all basins for the control run but are greater in the perturbation experiment, as expected. Note that the initMIP experiments use schematic forcing, and results should not be interpreted as realistic ice sheet or sea level projections. Our aim here is to demonstrate that when applied to largescale, wholeicesheet simulations on realistic geometries, MALI is robust and evolves reasonably during multicenturylength, freerunning simulations.
MALI is the current land ice model component of the US Department of Energy's Energy Exascale Earth System Model (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 highresolution global simulations, and all components have a variableresolution mesh capability. The ocean (Ringler et al., 2013; Petersen et al., 2015, 2018) and sea ice (Turner et al., 2018) 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 by one MPAS component can easily be leveraged by the others. Note that each component of E3SM can be run on differing numbers of processors within the coupled model, including the individual MPAS cores.
Physics at the ice sheet–atmosphere interface are handled by the snow model within the E3SM land model (ELM; Zhu et al., 2018; Ricciuto et al., 2018). ELM's snow model calculates ice sheet surface mass balance using a surface energy balance model and, at each coupling interval, MALI passes the current ice sheet extent and surface elevation through the coupler to ELM. The coupler then returns the surface mass balance and surface temperature calculated by ELM to MALI. These fields are used within MALI as boundary conditions to the mass and thermal evolution equations (Sect. 4.3 and 4.4). Currently, runoff from surface melting is calculated within ELM and routed directly through E3SM's runoff model, rather than being passed to and used by MALI. The subglacial discharge model discussed above in Sect. 5.1 is not currently coupled to the rest of E3SM.
Ongoing and future work on MALI and E3SM coupling includes the following: 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); twoway coupling between the ocean and a dynamic MALI model^{8}; and discharge of icebergs (solid ice flux from MALI) to the coupler and from there to the ocean and sea ice models.
Detailed analysis of the performance and scalability of the AlbanyLI velocity solver for idealized test cases and realistic highresolution applications to Greenland and Antarctica has been reported by Tezaur et al. (2015a, b) and Tuminaro et al. (2016). Because the momentum balance solver is ≥95 % of the cost of a typical forward model time step outside of I/O, the previously reported model performance is generally representative of overall MALI performance. To provide some additional context we summarize MALI performance for the highresolution Antarctica application described in the previous section. That mesh contained 1 642 490 horizontal grid cells and 11 vertical interfaces (10 vertical levels) at which the two horizontal components of velocity are solved. The simulations described above were run on the Edison Cray XC30 supercomputer^{9} at the National Energy Research Scientific Computing Center (NERSC). Computational nodes on Edison each contain two 12core Intel “Ivy Bridge” processors operating at 2.4 GHz and 64 GB DDR3 1866 MHz memory. Simulations were done using 6400 processors. The control simulation averaged 5.26 simulated years per wallclock hour (SYPH) over 2031 time steps, and the subiceshelf melt perturbation experiment averaged 4.61 SYPH over 2181 time steps. The differing performance is partially due to the higher number of time steps required by the perturbation experiment due to faster maximum ice velocity forcing the adaptive time stepper to take smaller time steps, but may also be a symptom of varying machine performance; performance during different segments of the simulations varied from 3.25–7.38 SYPH, presumably due to the usage of different node layouts on the machine and varying I/O performance. On average the velocity solve took 91.9 % of the computational time, writing output took 7.5 %, and all other operations the remaining 0.6 %. Total simulation cost for the 100year simulations was 122 000 core hours for the control run and 139 000 core hours for the perturbation simulation. For reference, the highresolution E3SM configuration (25 km resolution in atmosphere and land components, varying 18 to 6 km resolution in ocean and sea ice components) runs at 0.12 SYPH using 52 000 processors on Edison. While MALI could be run with substantially fewer processors to match the slower throughput of E3SM, the current optimal processor layout for highresolution E3SM could run MALI at the processor count we have done here without incurring any additional expense due to latent processors during model time stepping. At the resolution described here, MALI's computational cost of 1400 processor hours per simulated year would be insignificant compared to the cost of highresolution E3SM at 448 000 processor hours per simulated year. Of course running MALI within E3SM would restrict simulation lengths to those used for the coupled model (decades to centuries), which are too short for the investigation of many ice sheet science questions. At E3SM low resolution (100 km resolution in atmosphere and land components, varying 60 to 30 km resolution in ocean and sea ice components), the computational cost of MALI within E3SM remains a minor cost.
We have described MPASAlbany Land Ice (MALI), a higherorder, thermomechanically coupled ice sheet model using unstructured Voronoi meshes. MALI takes advantage of the MPAS Framework for the development of unstructured grid Earth system model components and the Albany and Trilinos frameworks for a parallel, performanceportable, implicit solution of the challenging higherorder ice sheet momentum balance through the AlbanyLI velocity solver. Together, these tools provide an accurate, efficient, scalable, and portable ice sheet model targeted for highresolution 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 threedimensional Blatter–Pattyn and shallow ice velocity solvers, a standard explicit mass evolution scheme, 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 massconserving subglacial hydrology model that can represent combinations of water drainage through till, distributed systems, and channelized systems, and it can be coupled to ice dynamics. A handful of basic calving schemes are currently implemented, including the physically based eigencalving method.
We have demonstrated the 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 semirealistic Antarctic ice sheet configuration at coarse resolution, and this capability was facilitated by the optimization tools within AlbanyLI.
A number of model improvements are planned over the next 5 years, focused heavily on improved representation of ice sheet physical processes and Earth system coupling. An implicit subglacial hydrology model based on such existing models (Werder et al., 2013; Hoffman and Price, 2014) is under development using the Albany framework. It will include optimization capabilities, a technique that has yet to be applied to subglacial hydrology beyond a spatially average, zerodimension application (Brinkerhoff et al., 2016). The difficulty in subglacial drainage parameter estimation remains one of the primary reasons 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 for future development, including gravitational, elastic, and viscous effects. Higherorder advection, through the fluxcorrected transport (Ringler et al., 2013) and/or incremental remapping (Dukowicz et al., 2010; Lipscomb and Hunke, 2004; Turner et al., 2018) schemes that are already implemented in MPAS, and semiimplicit time stepping are planned. Finally, a high priority is completing the coupling between ice sheet, ocean, and sea ice models in E3SM.
MPAS releases are available at https://mpasdev.github.io/ (last access: 10 September 2018) and model code is maintained at https://github.com/MPASDev/MPASRelease/releases (last access: 10 September 2018). MPASAlbany Land Ice is included in MPAS version 6.0. The digital object identifier for MPAS v6.0 is https://doi.org/10.5281/zenodo.1219886. MPAS is openly developed at https://github.com/MPASDev/MPASRelease (last access: 10 September 2018). The Albany library is developed openly at https://github.com/gahansen/Albany (last access: 10 September 2018), and the Trilinos library is developed openly at https://github.com/trilinos/Trilinos (last access: 10 September 2018). Region definitions for analysis are openly maintained at https://github.com/MPASDev/geometric_features (last access: 10 September 2018).
The authors declare that they have no conflict of interest.
This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the US Department of Energy or the United States Government.
Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the US Department of Energy's National Nuclear Security Administration under contract DENA0003525.
We thank additional contributors to the MALI code, Michael Duda,
Dominikus Heinzeller, Benjamin Hills, and Adrian Turner, as well as
the developers of the MPAS, Albany, and Trilinos libraries. Support for this work
was provided through the Scientific Discovery through Advanced Computing
(SciDAC) program and the Energy Exascale Earth System Model (E3SM) project
funded by the US Department of Energy (DOE), Office of Science, Biological
and Environmental Research, and Advanced Scientific Computing Research
programs. Development of the subglacial hydrology model was supported by a
grant to Matthew J. Hoffman from the Laboratory Directed Research and
Development Early Career Research Program at Los Alamos National Laboratory
(20160608ECR). This research used resources of the National Energy Research
Scientific Computing Center, a DOE Office of Science user facility supported
by the Office of Science of the US Department of Energy under contract
no. DEAC0205CH11231, and resources provided by the Los Alamos National
Laboratory Institutional Computing Program, which is supported by the US
Department of Energy National Nuclear Security Administration under contract
no. DEAC5206NA25396.
Edited by: Jeremy
Fyke
Reviewed by: Stephen Cornford and Thomas Zwinger
Adams, B., Bauman, L., Bohnhoff, W., Dalby, K., Ebeida, M., Eddy, J., Eldred, M., Hough, P., Hu, K., Jakeman, J., Swiler, L., and Vigil, D.: DAKOTA, A Multilevel Parallel ObjectOriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis: Version 5.4 User's Manual, Sandia Technical Report SAND20102183, 2013. a
Albrecht, T., Martin, M., Haseloff, M., Winkelmann, R., and Levermann, A.: Parameterization for subgridscale motion of iceshelf calving fronts, The Cryosphere, 5, 35–44, https://doi.org/10.5194/tc5352011, 2011. a
Arakawa, A. and Lamb, V. R.: Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, Academic Press, Inc., New York, 1977. a
AsayDavis, X. S., Cornford, S. L., Durand, G., GaltonFenzi, B. K., Gladstone, R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.: Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP+), ISOMIP v. 2 (ISOMIP+) and MISOMIP v. 1 (MISOMIP1), Geosci. Model Dev., 9, 2471–2497, https://doi.org/10.5194/gmd924712016, 2016. a
Aschwanden, A., Bueler, E., Khroulev, C., and Blatter, H.: An enthalpy formulation for glaciers and ice sheets, J. Glaciol., 58, 441–457, 2012. a, b, c
Åström, J. A., Vallot, D., Schäfer, M., Welty, E. Z., O'neel, S., Bartholomaus, T. C., Liu, Y., Riikilä, T. I., zwinger, T., Timonen, J., and Moore, J. C.: Termini of calving glaciers as selforganized critical systems, Nat. Geosci., 7, 874–878, 2014. a
Bassis, J. N. and Ma, Y.: Earth and Planetary Science Letters, Earth Planet. Sc. Lett., 409, 203–211, 2015. a, b
Bindschadler, R.: The importance of pressurized subglacial water in separation and sliding at the glacier bed, J. Glaciol., 29, 3–19, 1983. a
Bindschadler, R. A., Nowicki, S., AbeOuchi, A., Aschwanden, A., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Khroulev, C., Levermann, A., Lipscomb, W. H., Martin, M. A., Morlighem, M., Parizek, B. R., Pollard, D., Price, S. F., Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W. L.: Icesheet model sensitivities to environmental forcing and their use in projecting future sea level (the SeaRISE project), J. Glaciol., 59, 195–224, 2013. a
Blatter, H.: Velocity and StressFields in Grounded Glaciers – a Simple Algorithm for Including Deviatoric Stress Gradients, J. Glaciol., 41, 333–344, 1995. a
Bondzio, J. H., Seroussi, H., Morlighem, M., Kleiner, T., Rückamp, M., Humbert, A., and Larour, E. Y.: Modelling calving front dynamics using a levelset method: application to Jakobshavn Isbræ, West Greenland, The Cryosphere, 10, 497–510, https://doi.org/10.5194/tc104972016, 2016. a
Borstad, C., Khazendar, A., Scheuchl, B., Morlighem, M., Larour, E., and Rignot, E.: A constitutive framework for predicting weakening and reduced buttressing of ice shelves based on observations of the progressive deterioration of the remnant Larsen B Ice Shelf, Geophys. Res. Lett., 43, 2027–2035, https://doi.org/10.1002/2015GL067365, 2016. a
Borstad, C. P., Khazendar, a., Larour, E., Morlighem, M., Rignot, E., Schodlok, M. P., and Seroussi, H.: A damage mechanics assessment of the Larsen B ice shelf prior to collapse: Toward a physicallybased calving law, Geophys. Res. Lett., 39, L18502, https://doi.org/10.1029/2012GL053317, 2012. a
Brinkerhoff, D. J. and Johnson, J. V.: Data assimilation and prognostic whole ice sheet modelling with the variationally derived, higher order, open source, and fully parallel ice sheet model VarGlaS, The Cryosphere, 7, 1161–1184, https://doi.org/10.5194/tc711612013, 2013. a, b, c, d
Brinkerhoff, D. J., Meyer, C. R., Bueler, E., Truffer, M., and Bartholomaus, T. C.: Inversion of a glacier hydrology model, Ann. Glaciol., 57, 84–95, https://doi.org/10.1017/aog.2016.3, 2016. a
Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, J. Geophys. Res., 114, 1–21, 2009. a
Bueler, E. and van Pelt, W.: Massconserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635, https://doi.org/10.5194/gmd816132015, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m
Bueler, E., Lingle, C. S., KallenBrown, J. A., Covey, D. N., and Bowman, L. N.: Exact solutions and verification of numerical models for isothermal ice sheets, J. Glaciol., 51, 291–306, https://doi.org/10.3189/172756505781829449, 2005. a, b
Bueler, E., Brown, J., and Lingle, C.: Exact solutions to the thermomechanically coupled shallowice approximation: effective tools for verification, J. Glaciol., 53, 499–516, https://doi.org/10.3189/002214307783258396, 2007. a, b, c, d
Butler, H., Daly, M., Doyle, A., Gillies, S., Hagen, S., and Schaub, T.: The GeoJSON Format, Tech. rep., Hobu Inc., available at: https://www.rfceditor.org/rfc/rfc7946.txt (last access: 10 September 2018), 2016. a
Clarke, G. K.: Subglacial Processes, Annu. Rev. Earth Pl. Sc., 33, 247–276, https://doi.org/10.1146/annurev.earth.33.092203.122621, 2005. a
Cornford, S. L., Martin, D. F., Graves, D. T., Ranken, D. F., Le Brocq, A. M., Gladstone, R. M., Payne, A. J., Ng, E. G., and Lipscomb, W. H.: Adaptive mesh, finite volume modeling of marine ice sheets, J. Comput. Phys., 232, 529–549, 2013. a
Cornford, S. L., Martin, D. F., Payne, A. J., Ng, E. G., Le Brocq, A. M., Gladstone, R. M., Edwards, T. L., Shannon, S. R., Agosta, C., van den Broeke, M. R., Hellmer, H. H., Krinner, G., Ligtenberg, S. R. M., Timmermann, R., and Vaughan, D. G.: Centuryscale simulations of the response of the West Antarctic Ice Sheet to a warming climate, The Cryosphere, 9, 1579–1600, https://doi.org/10.5194/tc915792015, 2015. a
Cuffey, K. and Paterson: The Physics of Glaciers, ButterworthHeinneman, Amsterdam, 4th edn., 2010. a
Demeshko, I., Watkins, J., Tezaur, I. K., Guba, O., Spotz, W. F., Salinger, A. G., Pawlowski, R. P., and Heroux, M. A.: Toward performance portability of the Albany finite element analysis code using the Kokkos library, Int. J. High Perform C., https://doi.org/10.1177/1094342017749957, online first, 2018. a
Dennis, J. M., Edwards, J., Loy, R., Jacob, R., Mirin, A. A., Craig, A. P., and Vertenstein, M.: An applicationlevel parallel I/O library for Earth system models, Int. J. High Perform C., 26, 43–53, https://doi.org/10.1177/1094342011428143, 2012. a
Du, Q. and Gunzburger, M.: Grid generation and optimization based on centroidal Voronoi tessellations, Appl. Math. Comput., 133, 591–607, 2002. a
Duddu, R., Bassis, J. N., and Waisman, H.: A numerical investigation of surface crevasse propagation in glaciers using nonlocal continuum damage mechanics, Geophys. Res. Lett., 40, 3064–3068, https://doi.org/10.1002/grl.50602, 2013. a
Dukowicz, J. K., Price, S. F., and Lipscomb, W. H.: Consistent approximations and boundary conditions for icesheet dynamics from a principle of least action, J. Glaciol., 56, 480–496, https://doi.org/10.3189/002214310792447851, 2010. a, b, c
Edwards, H. C., Trott, C. R., and Sunderland, D.: Kokkos: Enabling manycore performance portability through polymorphic memory access patterns, J. Parallel Distr. Com., 74, 3202–3216, https://doi.org/10.1016/j.jpdc.2014.07.003, 2014a. a
Edwards, T. L., Fettweis, X., Gagliardini, O., GilletChaulet, F., Goelzer, H., Gregory, J. M., Hoffman, M., Huybrechts, P., Payne, A. J., Perego, M., Price, S., Quiquet, A., and Ritz, C.: Effect of uncertainty in surface mass balance–elevation feedback on projections of the future sea level contribution of the Greenland ice sheet, The Cryosphere, 8, 195–208, https://doi.org/10.5194/tc81952014, 2014b. a, b
Egholm, D. L. and Nielsen, S. B.: An adaptive finite volume solver for ice sheets and glaciers, J. Geophys. Res., 115, F01006, https://doi.org/10.1029/2009JF001394, 2010. a
Engwirda, D.: JIGSAWGEO (1.0): locally orthogonal staggered unstructured grid generation for general circulation modelling on the sphere, Geosci. Model Dev., 10, 2117–2140, https://doi.org/10.5194/gmd1021172017, 2017a. a, b
Engwirda, D.: JIGSAW(GEO): Unstructured grid generation for geophysical modelling, available at: https://github.com/dengwirda/jigsawgeomatlab (last access: 10 September 2018), 2017b. a, b
Feldmann, J. and Levermann, A.: Collapse of the West Antarctic Ice Sheet after local destabilization of the Amundsen Basin, P. Natl. Acad. Sci. USA, 112, 14191–14196, https://doi.org/10.1073/pnas.1512482112, 2015. a, b
Flowers, G. E.: Modelling water flow under glaciers and ice sheets, P. Roy. Soc. AMath. Phy., 471, 20140907, https://doi.org/10.1098/rspa.2014.0907, 2015. a, b, c, d
Flowers, G. E. and Clarke, G. K.: A multicomponent coupled model of glacier hydrology 1. Theory and synthetic examples, J. Geophys. Res., 107, 2287, https://doi.org/10.1029/2001JB001122, 2002. a
Fowler, A. C. and Larson, D. A.: On the flow of polythermal glaciers I. Model and preliminary analysis, P. Roy. Soc. AMath. Phy., 363, 217–242, https://doi.org/10.1098/rspa.1983.0054, 1978. a
Fyke, J. G., Weaver, A. J., Pollard, D., Eby, M., Carter, L., and Mackintosh, A.: A new coupled ice sheet/climate model: description and sensitivity to model physics under Eemian, Last Glacial Maximum, late Holocene and modern climate conditions, Geosci. Model Dev., 4, 117–136, https://doi.org/10.5194/gmd41172011, 2011. a
Gagliardini, O., Cohen, D., Råback, P., and Zwinger, T.: Finiteelement modeling of subglacial cavities and related friction law, J. Geophys. Res., 112, F02027, https://doi.org/10.1029/2006JF000576, 2007. a
Gagliardini, O., Zwinger, T., GilletChaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a newgeneration ice sheet model, Geosci. Model Dev., 6, 1299–1318, https://doi.org/10.5194/gmd612992013, 2013. a, b
Gladstone, R. M., Lee, V., Vieli, A., and Payne, A. J.: Grounding line migration in an adaptive mesh ice sheet model, J. Geophys. Res., 115, F04014, https://doi.org/10.1029/2009JF001615, 2010. a
Glen, J. W.: The Creep of Polycrystalline Ice, P. Roy. Soc. AMath. Phy., 228, 519–538, https://doi.org/10.1098/rspa.1955.0066, 1955. a
Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., AbeOuchi, A., Aschwanden, A., Calov, R., Gagliardini, O., GilletChaulet, F., Golledge, N. R., Gregory, J., Greve, R., Humbert, A., Huybrechts, P., Kennedy, J. H., Larour, E., Lipscomb, W. H., Le clec'h, S., Lee, V., Morlighem, M., Pattyn, F., Payne, A. J., Rodehacke, C., Rückamp, M., Saito, F., Schlegel, N., Seroussi, H., Shepherd, A., Sun, S., van de Wal, R., and Ziemen, F. A.: Design and results of the ice sheet model initialisation experiments initMIPGreenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460, https://doi.org/10.5194/tc1214332018, 2018. a
Goldberg, D. N.: A variationally derived, depthintegrated approximation to a higherorder glaciological flow model, J. Glaciol., 57, 157–170, 2011. a
Goldberg, D. N. and Sergienko, O. V.: Data assimilation using a hybrid ice flow model, The Cryosphere, 5, 315–327, https://doi.org/10.5194/tc53152011, 2011. a
Greve, R. and Blatter, H.: Comparison of thermodynamics solvers in the polythermal ice sheet model SICOPOLIS, Polar Science, 10, 11–23, 2016. a
Halfar, P.: On the dynamics of the ice sheets, J. Geophys. Res., 86, 11065–11072, https://doi.org/10.1029/JC088iC10p06043, 1981. a
Halfar, P.: On the Dynamics of the Ice Sheets 2, J. Geophys. Res., 88, 6043–6051, 1983. a, b
Heroux, M., Bartlett, R., Howle, V., Vicki, E., Hoekstra, R., Hu, J., Kolda, T., Lehoucq, R., Long, K., Pawlowski, R., Phipps, E., Salinger, A., Thornquist, H., Tuminaro, R., Willenbring, J., Williams, A., and Stanley, K.: An overview of the Trilinos project, ACM Trans. Math. Softw., 31, 397–423, 2005. a
Hewitt, I. J.: Modelling distributed and channelized subglacial drainage: the spacing of channels, J. Glaciol., 57, 302–314, https://doi.org/10.3189/002214311796405951, 2011. a, b, c, d
Hewitt, I. J.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sc. Lett., 371–372, 16–25, 2013. a, b, c
Hindmarsh, R. C. and Payne, A. J.: Timestep limits for stable solutions of the icesheet equation, Ann. Glaciol., 23, 74–85, 1996. a
Hoffman, M. and Price, S.: Feedbacks between coupled subglacial hydrology and glacier dynamics, J. Geophys. Res.–Earth, 119, 414–436, 2014. a, b, c, d, e
Hoffman, M. J., Perego, M., Andrews, L. C., Price, S. F., Neumann, T. A., Johnson, J. V., Catania, G., and Lüthi, M. P.: Widespread moulin formation during supraglacial lake drainages in Greenland, Geophys. Res. Lett., 45, 778–788, https://doi.org/10.1002/2017GL075659, 2018. a
Hutter, K.: Theoretical glaciology; material science of ice and the mechanics of glaciers and ice sheets, Reidel Publishing Co., Terra Scientific Publishing Co., Tokyo, 1983. a
Huybrechts, P., Payne, T., and The EISMINT Intercomparison Group: The EISMINT benchmarks for testing icesheet models, Ann. Glaciol., 23, 1–12, 1996. a, b, c, d, e
IMBIE team: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 558, 219–222, https://doi.org/10.1038/s415860180179y, 2018. a
IPCC: Climate Change 2007 – The Physical Science Basis: Working Group I Contribution to the Fourth Assessment Report of the IPCC, Cambridge University Press, Cambridge, UK, New York, NY, USA, 2007. a
IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK, New York, NY, USA, https://doi.org/10.1017/CBO9781107415324, 2013. a
Jiménez, S., Duddu, R., and Bassis, J.: An updatedLagrangian damage mechanics formulation for modeling the creeping flow and fracture of ice sheets, Comput. Methods Appl. Mech. Eng., 313, 406–432, 2017. a
Kleiner, T., Rückamp, M., Bondzio, J. H., and Humbert, A.: Enthalpy benchmark experiments for numerical ice sheet models, The Cryosphere, 9, 217–228, https://doi.org/10.5194/tc92172015, 2015. a, b, c, d
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res., 117, F01022, https://doi.org/10.1029/2011JF002140, 2012. a, b
Le Brocq, A., Payne, A., Siegert, M., and Alley, R.: A subglacial waterflow model for West Antarctica, J. Glaciol., 55, 879–888, https://doi.org/10.3189/002214309790152564, 2009. a
Leguy, G. R.: The Effect of a Basalfriction Parameterization on Groundingline Dynamics in Icesheet Models, PhD thesis, New Mexico Institute of Mining and Technology, 2015. a
Lenaerts, J. T. M., Van Den Broeke, M. R., Van De Berg, W. J., Van Meijgaard, E., and Kuipers Munneke, P.: A new, highresolution surface mass balance map of Antarctica (1979–2010) based on regional atmospheric climate modeling, Geophys. Res. Lett., 39, L04501, https://doi.org/10.1029/2011GL050713, 2012. a, b, c
Leng, W., Ju, L., Gunzburger, M., Price, S., and Ringler, T.: A parallel highorder accurate finite element nonlinear Stokes ice sheet model and benchmark experiments, J. Geophys. Res., 117, F01001, https://doi.org/10.1029/2011JF001962, 2012. a
Levermann, A., Albrecht, T., Winkelmann, R., Martin, M. A., Haseloff, M., and Joughin, I.: Kinematic firstorder calving law implies potential for abrupt iceshelf retreat, The Cryosphere, 6, 273–286, https://doi.org/10.5194/tc62732012, 2012. a, b
Lipscomb, W. H. and Hunke, E. C.: Modeling Sea Ice Transport Using Incremental Remapping, Mon. Weather Rev., 132, 1341–1354, https://doi.org/10.1175/15200493(2004)132<1341:MSITUI>2.0.CO;2, 2004. a
Lipscomb, W. H., Fyke, J. G., Vizcaíno, M., Sacks, W. J., Wolfe, J., Vertenstein, M., Craig, A., Kluzek, E., and Lawrence, D. M.: Implementation and Initial Evaluation of the Glimmer Community Ice Sheet Model in the Community Earth System Model, J. Climate, 26, 7352–7371, 2013. a
Lipscomb, W. H., Price, S. F., Hoffman, M. J., Leguy, G. R., Bennett, A. R., Bradley, S. L., Evans, K. J., Fyke, J. G., Kennedy, J. H., Perego, M., Ranken, D. M., Sacks, W. J., Salinger, A. G., Vargo, L. J., and Worley, P. H.: Description and Evaluation of the Community Ice Sheet Model (CISM) v2.1, Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd2018151, in review, 2018. a, b, c
Little, C. M., Oppenheimer, M., Alley, R. B., Balaji, V., Clarke, G. K. C., Delworth, T. L., Hallberg, R., Holland, D. M., Hulbe, C. L., Jacobs, S., Johnson, J. V., Levy, H., Lipscomb, W. H., Marshall, S. J., Parizek, B. R., Payne, A. J., Schmidt, G. A., Stouffer, R. J., Vaughan, D. G., and Winton, M.: Toward a new generation of ice sheet models, Eos T. Am. Geophys. Un., 88, 578–579, https://doi.org/10.1029/2007EO520002, 2007. a
Morland, L. W. and Johnson, I. R.: Steady motion of ice sheets., J. Glaciol., 25, 229–246, 1980. a
Nowicki, S., Bindschadler, R. A., AbeOuchi, A., et al.: Insights into spatial sensitivities of ice mass response to environmental change from the SeaRISE ice sheet modeling project II: Greenland, J. Geophys. Res., 118, 1025–1044, 2013a. a
Nowicki, S., Bindschadler, R. A., AbeOuchi, A., et al.: Insights into spatial sensitivities of ice mass response to environmental change from the SeaRISE ice sheet modeling project I: Antarctica, J. Geophys. Res., 118, 1002–1024, 2013b. a
Nowicki, S. M. J., Payne, A., Larour, E., Seroussi, H., Goelzer, H., Lipscomb, W., Gregory, J., AbeOuchi, A., and Shepherd, A.: Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6, Geosci. Model Dev., 9, 4521–4545, https://doi.org/10.5194/gmd945212016, 2016. a
Pattyn, F.: A new threedimensional higherorder thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes, J. Geophys. Res., 108, 1–15, https://doi.org/10.1029/2002JB002329, 2003. a
Pattyn, F., Perichon, L., Aschwanden, A., Breuer, B., de Smedt, B., Gagliardini, O., Gudmundsson, G. H., Hindmarsh, R. C. A., Hubbard, A., Johnson, J. V., Kleiner, T., Konovalov, Y., Martin, C., Payne, A. J., Pollard, D., Price, S., Rückamp, M., Saito, F., Soucek, O., Sugiyama, S., and Zwinger, T.: Benchmark experiments for higherorder and fullStokes ice sheet models (ISMI–HOM), The Cryosphere, 2, 95–108, https://doi.org/10.5194/tc2952008, 2008. a, b
Pattyn, F., Perichon, L., Favier, L., et al.: Groundingline migration in planview marine icesheet models: results of the ice2sea MISMIP3d intercomparison, J. Glaciol., 59, 410–422, 2013. a, b, c, d, e, f, g, h, i, j
Pawlowski, R. P., Phipps, E. T., Salinger, A. G., Owen, S. J., Siefert, C. M., and Staten, M. L.: Automating embedded analysis capabilities and managing software complexity in multiphysics simulation, Part II: Application to partial differential equations, Sci. ProgrammingNeth., 20, 327–345, https://doi.org/10.3233/SPR20120351, 2012. a
Payne, A. J., Huybrechts, P., Calov, R., Fastook, J. L., Greve, R., Marshall, S. J., Marsiat, I., Ritz, C., Tarasov, L., and Thomassen, M. P. A.: Results from the EISMINT model intercomparison: The effects of thermomechanical coupling, J. Glaciol., 46, 227–238, 2000. a, b, c, d, e, f, g, h, i, j, k
Perego, M., Gunzburger, M., and Burkardt, J.: Parallel finiteelement implementation for higherorder icesheet models, J. Glaciol., 58, 76–88, https://doi.org/10.3189/2012JoG11J063, 2012. a, b, c
Perego, M., Price, S., and Stadler, G.: Optimal initial conditions for coupling ice sheet models to Earth system models, J. Geophys. Res.Earth, 119, 1–24, https://doi.org/10.1002/2014JF003181.Received, 2014. a, b
Petersen, M. R., Jacobsen, D. W., Ringler, T. D., Hecht, M. W., and Maltrud, M. E.: Evaluation of the arbitrary LagrangianEulerian vertical coordinate method in the MPASOcean model, Ocean Model., 86, 93–113, https://doi.org/10.1016/j.ocemod.2014.12.004, 2015. a, b, c
Petersen, M. R., AsayDavis, X., Jacobsen, D., Jones, P., Maltrud, M., Ringler, T. D., Turner, A. K., Van Roekel, L., Veneziani, M., Wolfe, J., and Wolfram, P. J.: An evaluation of the ocean and sea ice climate of E3SM using MPAS and interannual COREII forcing, J. Adv. Model. Earth Sy., in review, 2018. a, b
Price, S., Flowers, G., and Schoof, C.: Improving Hydrology in Land Ice Models, Eos T. Am. Geophys. Un., 92, p. 164, 2011. a
Price, S., Lipscomb, W., Hoffman, M., Hagdorn, M., Rutt, I., Payne, T., Hebeler, F., and Kennedy, J. H.: Community Ice Sheet Model (CISM) v2.0.5 Documentation, Tech. rep., Los Alamos National Laboratory, Los Alamos, NM, available at: https://cism.github.io/data/cism_documentation_v2_1.pdf (last access: 10 September 2018), 2015. a, b
Ricciuto, D., Sargsyan, K., and Thornton, P.: The Impact of Parametric Uncertainties on Biogeochemistry in the E3SM Land Model, JAMES, 10, 297–319, https://doi.org/10.1002/2017MS000962, 2018. a
Ridley, J. K., Huybrechts, P., and Gregory, J. M.: Elimination of the Greenland ice sheet in a high CO_{2} climate, J. Climate, 18, 3409–3427, 2005. a
Rignot, E., Mouginot, J., and Scheuchl, B.: Ice Flow of the Antarctic Ice Sheet, Science, 333, 1427–1430, 2011. a, b
Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: IceShelf Melting Around Antarctica, Science, 341, 266–270, 2013. a, b
Ringler, T., Petersen, M., Higdon, R. L., Jacobsen, D., Jones, P. W., and Maltrud, M.: A multiresolution approach to global ocean modeling, Ocean Model., 69, 211–232, https://doi.org/10.1016/j.ocemod.2013.04.010, 2013. a, b, c, d
Ringler, T. D., Thuburn, J., Klemp, J. B., and Skamarock, W. C.: A unified approach to energy conservation and potential vorticity dynamics for arbitrarilystructured Cgrids, J. Computat. Phys., 229, 3065–3090, https://doi.org/10.1016/j.jcp.2009.12.007, 2010. a
Ringler, T. D., Jacobsen, D., Gunzburger, M., Ju, L., Duda, M., and Skamarock, W.: Exploring a Multiresolution Modeling Approach within the ShallowWater Equations, Mon. Weather Rev., 139, 3348–3368, https://doi.org/10.1175/MWRD1005049.1, 2011. a
Rutt, I. C., Hagdorn, M., Hulton, N. R. J., and Payne, A. J.: The Glimmer community ice sheet model, J. Geophys. Res., 114, F02004, https://doi.org/10.1029/2008JF001015, 2009. a
Saito, F., Abeouchi, A., and Blatter, H.: European Ice Sheet Modelling Initiative (EISMINT) model intercomparison experiments with firstorder mechanics, J. Geophys. Res., 111, 1–9, https://doi.org/10.1029/2004JF000273, 2006. a
Salinger, A. G., Bartlett, R. A., Bradley, A. M., Chen, Q., Demeshko, I. P., Gao, X., Hansen, G. A., Mota, A., Muller, R. P., Nielsen, E., Ostien, J. T., Pawlowski, R. P., Perego, M., Phipps, E. T., Sun, W., and Tezaur, I. K.: Albany: Using componentbased design to develop a flexible, generic multiphysics analysis core, Int. J. Multiscale Com., 14, 415–438, 2016. a
Schoof, C.: The effect of cavitation on glacier sliding, P. Roy. Soc. AMath. Phy., 461, 609–627, 2005. a, b
Schoof, C.: Icesheet acceleration driven by melt supply variability, Nature, 468, 803–806, https://doi.org/10.1038/nature09618, 2010. a, b, c
Schoof, C. and Hewitt, I.: Icesheet dynamics, Annu. Rev. Fluid Mech., 45, 217–239, 2013. a
Schoof, C. and Hindmarsh, R. C. A.: ThinFilm Flows with Wall Slip: An Asymptotic Analysis of Higher Order Glacier Flow Models, Q. J. Mech. Appl. Math., 63, 73–114, 2010. a, b
Schoof, C., Hewitt, I. J., and Werder, M. A.: Flotation and free surface flow in a model for subglacial drainage, Part 1. Distributed drainage, J. Fluid Mech., 702, 126–156, https://doi.org/10.1017/jfm.2012.165, 2012. a
Seroussi, H., Morlighem, M., Rignot, E., Khazendar, A., Larour, E., and Mouginot, J.: Dependence of centuryscale projections of the Greenland ice sheet on its thermal regime, J. Glaciol., 59, 1024–1034, https://doi.org/10.3189/2013JoG13J054, 2013. a
Seroussi, H., Morlighem, M., Larour, E., Rignot, E., and Khazendar, A.: Hydrostatic grounding line parameterization in ice sheet models, The Cryosphere, 8, 2075–2087, https://doi.org/10.5194/tc820752014, 2014. a
Shannon, S. R., Payne, A. J., Bartholomew, I. D., Van Den Broeke, M. R., Edwards, T. L., Fettweis, X., Gagliardini, O., GilletChaulet, F., Goelzer, H., Hoffman, M. J., Huybrechts, P., Mair, D. W. F., Nienow, P. W., Perego, M., Price, S. F., Smeets, C. J. P. P., Sole, A. J., van de Wal, R. S. W., and Zwinger, T.: Enhanced basal lubrication and the contribution of the Greenland ice sheet to future sealevel rise, P. Natl. Acad. Sci. USA, 110, 14156–14161, 2013. a, b
Shapiro, N. and Ritzwoller, M.: Inferring surface heat flux distributions guided by a global seismic model: particular application to Antarctica, Earth Planet. Sc. Lett., 223, 213–224, 2004. a
Shepherd, A., Ivins, E. R., Geruo, A., and Barletta, V. R.: A reconciled estimate of icesheet mass balance, Science, 338, 1183–1189, 2012. a
Shewchuk, J.: Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, Workshop on Applied Computational Geometry, 1148, 203–222, 1996. a
Skamarock, W. C., Klemp, J. B., Duda, M. G., Fowler, L. D., Park, S.H., and Ringler, T. D.: A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and CGrid Staggering, Mon. Weather Rev., 140, 3090–3105, https://doi.org/10.1175/MWRD1100215.1, 2012. a
Tezaur, I. K., Perego, M., Salinger, A. G., Tuminaro, R. S., and Price, S. F.: Albany/FELIX: a parallel, scalable and robust, finite element, firstorder Stokes approximation ice sheet solver built for advanced analysis, Geosci. Model Dev., 8, 1197–1220, https://doi.org/10.5194/gmd811972015, 2015a. a, b, c, d, e, f
Tezaur, I. K., Tuminaro, R. S., Perego, M., Salinger, A. G., and Price, S. F.: On the Scalability of the Albany/FELIX firstorder Stokes Approximation ice Sheet Solver for LargeScale Simulations of the Greenland and Antarctic ice Sheets, Procedia Comput. Sci., 51, 2026–2035, https://doi.org/10.1016/j.procs.2015.05.467, 2015b. a, b, c
Tulaczyk, S., Barclay, W., and Engelhardt, F.: Basal mechanics of Ice Stream B, West Antarctica: 1. Till mechanics, J. Geophys. Res., 105, 463–481, 2000. a
Tuminaro, R., perego, M., Tezaur, I., Salinger, A., and Price, S. F.: A matrix dependent/algebraic multigrid approach for extruded meshes with applications to ice sheet modeling, SIAM J. Sci. Comput., 38, 504–532, 2016. a, b, c
Turner, A. K., Lipscomb, W. H., Hunke, E. C., Jacobsen, D. W., Jeffery, N., Ringler, T. D., and Wolfe, J. D.: MPASSeaice: a new variable resolution seaice model, J. Adv. Model Earth Sy., in preparation, 2018. a, b, c
van der Veen, C. J.: Fundamentals of Glacier Dynamics, CRC Press, Boca Raton, FL, 2nd edn., 2013. a
Van Liefferinge, B. and Pattyn, F.: Using iceflow models to evaluate potential sites of million yearold ice in Antarctica, Clim. Past, 9, 2335–2345, https://doi.org/10.5194/cp923352013, 2013. a, b
Vizcaíno, M., Mikolajewicz, U., Gröger, M., MaierReimer, E., Schurgers, G., and Winguth, A. M. E.: Longterm ice sheet–climate interactions under anthropogenic greenhouse forcing simulated with a complex Earth System Model, Clim. Dynam., 31, 665–690, 2008. a
Vizcaíno, M., Mikolajewicz, U., Jungclaus, J., and Schurgers, G.: Climate modification by future ice sheet changes and consequences for ice sheet mass balance, Clim. Dynam., 34, 301–324, 2009. a
Werder, M. A., Hewitt, I. J., Schoof, C. G., and Flowers, G. E.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res.Earth, 118, 2140–2158, https://doi.org/10.1002/jgrf.20146, 2013. a, b, c, d, e, f, g
Zhu, Q., Riley, W., Tang, J., Collier, N., Hoffman, F., Randerson, J., Yang, X., and Bisht, G.: Representing carbon, nitrogen, and phosphorus interaction in the E3SM Land Model v1: Model development and global benchmarking, JAMES, in review, 2018. a
For example, traditional CPUonly architectures and MPI programming models versus CPU + GPU, hybrid architectures using MPI for nodal communication, and OpenMP or CUDA for onnode parallelism.
In practice, additional scaling parameters describing the ratio of deformation to sliding velocity may also be introduced.
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.
https://shmip.bitbucket.io/ (last access: 10 September 2018).
https://trilinos.org/packages/rol/ (last access: 10 September 2018).
https://github.com/MPASDev/geometric_features (last access: 10 Septemebr 2018).
Coupling to a static Antarctic ice sheet with ocean circulation in subiceshelf cavities is supported in E3SM version 1.0.0.
More information about Edison can be found at http://www.nersc.gov/users/computationalsystems/edison/configuration/ (last access: 10 September 2018).
 Abstract
 Introduction
 MPAS Framework
 The Albany software library
 Conservation equations
 Additional model physics
 Additional capabilities
 Model verification and benchmarks
 Realistic application: Antarctic ice sheet perturbation experiment
 Coupling to Energy Exascale Earth System Model
 Model performance
 Conclusions and future work
 Code availability
 Appendix A
 Competing interests
 Disclaimer
 Acknowledgements
 References
 Abstract
 Introduction
 MPAS Framework
 The Albany software library
 Conservation equations
 Additional model physics
 Additional capabilities
 Model verification and benchmarks
 Realistic application: Antarctic ice sheet perturbation experiment
 Coupling to Energy Exascale Earth System Model
 Model performance
 Conclusions and future work
 Code availability
 Appendix A
 Competing interests
 Disclaimer
 Acknowledgements
 References