2.1 Model meshes

The MPAS mesh specification is general enough to describe unstructured meshes on most two-dimensional 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 high-quality cells because cells tend towards equi-dimensional 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 two-dimensional manifold spaces, they are convenient for describing geophysical locations, either on planar projections or directly on a sphere. Because they are unstructured, meshes can contain varying mesh resolution and can be culled to only retain regions of interest. Planar meshes can easily be made periodic by taking advantage of the 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 C-grid discretization (Arakawa and Lamb, 1977) for advection, meaning state variables (ice thickness and tracer values) are located at Voronoi cell centers, and flow variables (transport velocity, u_n) are located at cell edge midpoints (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 variable-resolution centroidal Voronoi meshes. Additionally, the JIGSAW(GEO) mesh-generation tool (Engwirda, 2017a, b) can be used to efficiently generate high-quality, variable-resolution meshes with data-based 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.

Figure 1MALI grids. (a) Horizontal grid with cell center (blue circles), edge midpoint (red triangles), and vertices (orange squares) identified for the center cell. Scalar fields (H, T) are located at cell centers. Advective velocities (u_n) and fluxes are located at cell edges. (b) Vertical grid with layer midpoints (blue circles) and layer interfaces (red triangles) identified. Scalar fields (H, T) are located at layer midpoints. Fluxes are located at layer interfaces.

Download

4.1 Conservation of momentum

Treating glacier ice as an incompressible fluid in a low-Reynolds-number flow, the conservation of momentum in a Cartesian reference frame is expressed by the Stokes flow equations, for which the gravitational driving stress is balanced by gradients 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².

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{ϵ}}_{ij}}$ is the strain rate tensor and μ_e is the “effective” non-Newtonian ice viscosity given by Nye's generalization of Glen's flow law (Glen, 1955):

In Eq. (5), A is a temperature-dependent rate factor, n is an exponent commonly taken as 3 for polycrystalline glacier ice, and γ is an ice “stiffness” factor (inverse enhancement factor 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{ϵ}}}_{\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 reduced-order approximation applied and are discussed further below.

4.2 Reduced-order 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 – reduced-order 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 first-order-accurate “Blatter–Pattyn” approximation and a zero-order “shallow ice approximation” as described in more detail in the following sections.

4.2.1 First-order velocity solver and coupling

Ice sheets typically have a small aspect ratio and small surface and bed slopes. These characteristics imply that reduced-order approximations of the Stokes momentum balance may apply over large areas of the ice sheets, potentially allowing for significant computational savings. Formal derivations involve 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, reduced-order models with a chosen degree of accuracy (relative to the original Stokes flow equations) can be derived by retaining terms of the appropriate order in δ. For example, the first-order-accurate Stokes approximation is arrived at by retaining terms of 𝒪(δ¹) and lower (the reader is referred to Schoof and Hindmarsh, 2010, and Dukowicz et al., 2010, for additional discussion³).

Using the notation of Perego et al. (2012) and Tezaur et al. (2015a) ⁴, the first-order-accurate Stokes approximation (also referred to as the Blatter–Pattyn approximation; see Blatter, 1995; Pattyn, 2003) is expressed through the following system of PDEs:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{rcl}-\mathrm{\nabla }\cdot \left(\mathrm{2}{\mathit{\mu }}_{\text{e}}{\dot{\mathit{ϵ}}}_{\mathrm{1}}\right)+\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial x}}& =& \mathrm{0},\\ -\mathrm{\nabla }\cdot \left(\mathrm{2}{\mathit{\mu }}_{\text{e}}{\dot{\mathit{ϵ}}}_{\mathrm{2}}\right)+\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial y}}& =& \mathrm{0},\end{array}\right.\end{array}$$

where ∇⋅ is the divergence operator, $s\equiv s(x,y)$ represents the ice sheet upper surface, and the vectors ${\dot{\mathit{ϵ}}}_{\mathrm{1}}$ and ${\dot{\mathit{ϵ}}}_{\mathrm{2}}$ are given by

$$\begin{array}{}\text{(10)}& {\dot{\mathit{ϵ}}}_{\mathrm{1}}={\left(\mathrm{2}{\dot{\mathit{ϵ}}}_{xx}+{\dot{\mathit{ϵ}}}_{yy},{\dot{\mathit{ϵ}}}_{xy},{\dot{\mathit{ϵ}}}_{xz}\right)}^T\end{array}$$

and

$$\begin{array}{}\text{(11)}& {\dot{\mathit{ϵ}}}_{\mathrm{2}}={\left({\dot{\mathit{ϵ}}}_{xy},{\dot{\mathit{ϵ}}}_{xx}+\mathrm{2}{\dot{\mathit{ϵ}}}_{yy},{\dot{\mathit{ϵ}}}_{yz}\right)}^T.\end{array}$$

Akin to Eqs. (5) and (6), μ_e in Eq. (9) represents the effective viscosity but for the case of the first-order stress balance with an effective strain rate given by

$$\begin{array}{}\text{(12)}& {\dot{\mathit{ϵ}}}_{\text{e}}\equiv {\left({\dot{\mathit{ϵ}}}_{xx}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{yy}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{xx}{\dot{\mathit{ϵ}}}_{yy}+{\dot{\mathit{ϵ}}}_{xy}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{xz}^{\mathrm{2}}+{\dot{\mathit{ϵ}}}_{yz}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

rather than by Eq. (6), and with individual strain rate terms given by

$$\begin{array}{ll}\text{(13)}& {\displaystyle }& {\displaystyle }{\dot{\mathit{ϵ}}}_{xx}={\displaystyle \frac{\partial u}{\partial x}},{\dot{\mathit{ϵ}}}_{yy}={\displaystyle \frac{\partial v}{\partial y}},{\dot{\mathit{ϵ}}}_{xy}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right),{\displaystyle }& {\displaystyle }{\dot{\mathit{ϵ}}}_{xz}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial u}{\partial z}},{\dot{\mathit{ϵ}}}_{yz}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial v}{\partial z}}.\end{array}$$

At the upper surface, a stress-free boundary condition is applied,

$$\begin{array}{}\text{(14)}& {\dot{\mathit{ϵ}}}_{\mathrm{1}}\cdot \mathit{n}={\dot{\mathit{ϵ}}}_{\mathrm{2}}\cdot \mathit{n}=\mathrm{0},\end{array}$$

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”):

$$\begin{array}{}\text{(15)}& \begin{array}{ll}u=v=\mathrm{0},& \text{no slip}\\ \mathrm{2}{\mathit{\mu }}_{\text{e}}{\dot{\mathit{ϵ}}}_{\mathrm{1}}\cdot \mathit{n}+\mathit{\beta }{u}^m=\mathrm{0},& \text{sliding,}\\ \mathrm{2}\mathit{\mu }{\dot{\mathit{ϵ}}}_{\mathrm{2}}\cdot \mathit{n}+\mathit{\beta }{v}^m=\mathrm{0},& \end{array}\end{array}$$

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,

$$\begin{array}{}\text{(16)}& \mathrm{2}{\mathit{\mu }}_{\text{e}}{\left({\dot{\mathit{ϵ}}}_{\mathrm{1}}\cdot \mathit{n},{\dot{\mathit{ϵ}}}_{\mathrm{2}}\cdot \mathit{n},\mathrm{0}\right)}^T-\mathit{\rho }g(s-z)\mathit{n}={\mathit{\rho }}_{\text{o}}gmax(z,\mathrm{0})\mathit{n},\end{array}$$

where ρ_o is the density of ocean water and n the outward normal vector to the lateral boundary (i.e., parallel to the (x,y) plane) so that lateral boundaries above sea level are effectively stress free and lateral boundaries submerged within the ocean experience hydrostatic pressure due to the overlying column of ocean water.

We solve these equations using the Albany-LI momentum balance solver, which is built using the Albany and Trilinos software libraries discussed above. The mathematical formulation, discretization, solution methods, verification, and scaling of Albany-LI are discussed in detail in Tezaur et al. (2015a). Albany-LI implements a classic finite-element discretization of the first-order approximation. At the grounding line, the basal friction coefficient β can abruptly drop to zero within an element of the mesh. This discontinuity is resolved by using a higher-order Gauss quadrature rule on elements containing the grounding line, which corresponds to the sub-element parameterization SEP3 proposed in Seroussi et al. (2014). Additional exploration of solver scalability and demonstrations of solver robustness on large-scale, high-resolution, realistic problems are discussed in Tezaur et al. (2015b). The efficiency and robustness of 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 ill-conditioned linear systems.

Table 1Correspondence between the MPAS Voronoi tessellation and its dual Delaunay triangulation used by Albany. Key MALI variables that are natively found at each location are listed. Note that variables are interpolated from one location to another as required for various calculations.

Download Print Version | Download XLSX

The Albany-LI first-order velocity solver written in C++ is coupled to MPAS written in Fortran using an interface layer. Albany uses a three-dimensional mesh extruded from a basal triangulation and composed of prisms or tetrahedra (see Tezaur et al., 2015a). When coupled to MPAS, the basal triangulation is 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 three-dimensional ice temperature, all at cell centers (Table 1), are passed from MPAS to Albany. Optionally, Dirichlet velocity boundary conditions can also be passed. After the velocity solve is complete, Albany 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 finite-element mesh that Albany will use. Lateral boundaries in Albany are applied at cell centers (triangle nodes) that do not contain dynamic ice on the MPAS mesh and that are adjacent to the last cell of the MPAS mesh that does contain dynamic ice. This one element extension is required to support 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 first-order 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.

Figure 2Correspondence between MPAS and Albany meshes and the application of boundary conditions for the first-order velocity solver. Solid black lines are cells on the Voronoi mesh and dashed gray lines are triangles on the Delaunay triangulation. Light blue Voronoi cells contain dynamic ice and gray cells do not. Dark blue circles are Albany triangle nodes that use variable values directly from the colocated MPAS cell centers. White circles are extended node locations that receive variable values as described in the text based on whether they are terrestrial, floating marine, or grounded marine locations. Red triangles indicate Voronoi cell edges on which velocities (u_n) are required for advection.

Download

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{‾}}{\mathit{u}}$ is depth-averaged 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 finite-volume method such that fluxes are calculated for each edge of each cell to calculate $\mathrm{\nabla }\cdot H\stackrel{\mathrm{‾}}{\mathit{u}}$. Specifically, we use a first-order upwind method that applies the normal velocity on each edge (u_n) and an upwind value of cell-centered ice thickness. Note that with the Blatter–Pattyn velocity solver, normal velocity is interpolated from cell centers to edges using the finite-element basis functions in Albany. In the shallow ice approximation velocity solver, normal velocity is calculated natively at edges. The MPAS Framework includes a higher-order flux-corrected 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 first-order 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).

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 mass-conserving 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.7 Verification and real-world application

To verify the implementation of the distributed drainage model, we use the nearly exact solution described by Bueler and van Pelt (2015). The problem configuration uses distributed drainage only on a two-dimensional, radially symmetric ice sheet of radius 22.5 km with parabolic ice sheet thickness and a nontrivial sliding profile. Bueler and van Pelt (2015) showed that this configuration allows for nearly exact reference values of W and P_w to be solved at steady state from an ordinary differential equation initial value problem with very high accuracy. We follow the test protocol of Bueler and van Pelt (2015) and initialize the model with the near-exact solution and then run the model forward for 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 first-order convergence (Fig. 3).

Figure 3Error in subglacial hydrology model for radial test case with the near-exact solution described by Bueler and van Pelt (2015) for different grid resolutions. (a) Error in water thickness; x symbols indicate maximum error, and squares indicate mean error. Average error in water thickness decays as 𝒪(Δx^0.97). (b) Error in water pressure, with same symbols. Average error in water pressure decays as 𝒪(Δx^1.02).

Download

To check the model implementation of channels, we use comparisons to other more mature drainage models through the Subglacial Hydrology Model Intercomparison Project (SHMIP)⁵. Steady-state solutions of the drainage system effective pressure, water fluxes, and channel development for an idealized ice sheet with varying magnitudes of meltwater input (SHMIP experiment suites A and B) compared between MALI and other models of similar complexity (GlaDS, Elmer) are very similar.

To demonstrate a real-world application of the subglacial hydrology model, we perform a stand-alone 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 first-order 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 down-glacier and is greatest in fast-flowing 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 subglacial-hydrology–ice-dynamics simulations are beyond the scope of this paper; we merely mean to demonstrate plausible behavior from the subglacial hydrology model for a realistic ice-sheet-scale problem.

Figure 4Subglacial hydrology model results for the 20 km resolution Antarctic ice sheet. Demonstration of subglacial hydrology model capability using a 20 km resolution simulation of Antarctica (too coarse resolution for scientific validity but sufficient for demonstrating model capabilities). (a) Grounded basal ice speed calculated by the first-order velocity solver optimized to surface velocity observations. This field and the calculated basal melt are the forcings applied to the stand-alone subglacial hydrology model. (b) Water flux in the distributed system calculated by the subglacial hydrology model at steady state. (Ice dynamics is prescribed.)

Download

5.2 Iceberg calving

MALI includes a few simple methods for removing ice from calving fronts during each model time step.

1. All floating ice is removed.

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

3. All floating ice thinner than a specified threshold is removed.

4. The calving front is maintained at its initial location by adding or removing ice after thickness evolution is complete. 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 spin-up).

5. Eigencalving scheme (Levermann et al., 2012). Calving front retreat rate, C_v, is proportional to the product of the principal strain rates ($\dot{{\mathit{ϵ}}_{\mathrm{1}}},\dot{{\mathit{ϵ}}_{\mathrm{2}}}$) if they are both extensional:

   $$\begin{array}{}\text{(57)}& C_{\text{v}}=K_{\mathrm{2}}\dot{{\mathit{ϵ}}_{\mathrm{1}}}\dot{{\mathit{ϵ}}_{\mathrm{2}}}\text{for}\dot{{\mathit{ϵ}}_{\mathrm{1}}}>\mathrm{0}\text{and}\dot{{\mathit{ϵ}}_{\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 sub-grid motion of the calving front by forcing floating cells adjacent to open ocean to remain dynamically inactive until ice thickness there reaches 95 % of the minimum thickness of all floating neighbors. This is an ad hoc alternative to methods tracking the calving front position at sub-grid 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 real-world application of the eigencalving parameterization, we perform a 1000-year spin-up 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 1000-year spin-up, we apply steady forcing of present-day 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₂ parameter tuned individually for large ice shelves and a minimum calving front thickness threshold of 100 m. This spin-up, 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).

Figure 5Demonstration of eigencalving capability using a 20 km resolution simulation of Antarctica (too coarse resolution for scientific validity but sufficient for demonstrating model capabilities). (a) Modeled ice extent and surface speed after optimization. The white contour line in each plot is the grounding line. Areas colored red exceed the maximum speed shown in the color bar. Gray areas are ice-free regions of the computational domain. (b) Modeled ice extent and surface speed after 1000 years with evolving velocity, geometry, and temperature and active eigencalving, plotted as in (a).

Download

6.1 Optimization

MALI includes an optimization capability through its coupling to the Albany-LI 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 right-hand side (RHS) is a cost function associated with the misfit between modeled and observed surface velocities, the second term on the RHS is a cost function associated with the ice stiffness factor, γ (see 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 trade-off between a smooth β field and one with higher-frequency oscillations (that may capture more spatial detail at the risk of over-fitting the observations). The optimal values of α_β and α_γ can be chosen through a standard L-curve analysis. The optimization problem is solved using the limited-memory BFGS method, as implemented in the Trilinos package ROL⁶, on the reduced-space problem. The functional gradient is computed using the adjoint method.

An example application of the optimization capability applied to a realistic, whole-ice-sheet 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 observational-based 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 quasi-equilibrium 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 run-time-generated simulation diagnostics and statistics output at user-specified time intervals. Support tools included with the code release allow for the definition of any number or combination of 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⁷. Python-based 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 ground-penetrating radar profile lines). In Sect. 8 below we demonstrate the analysis capability applied to an idealized simulation of the Antarctica ice sheet.

Table 2Standard model diagnostics available for an arbitrary number of predefined geographic regions.

Download Print Version | Download XLSX

7.1 Halfar analytic solution

In Halfar (1981, 1983), Halfar described an analytic solution for the time-evolving geometry of a radially symmetric, isothermal dome of ice on a flat bed with no accumulation flowing under the shallow ice approximation. This 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 zero-accumulation member of a family of analytic solutions, but we apply the original Halfar test here.

Figure 6(a) Root mean square error in ice thickness as a function of grid cell spacing for the Halfar dome after 200 years shown with black dots. The order of convergence is 0.78. The red square shows the RMS thickness error for the variable-resolution mesh shown in (b) with 1000 m spacing around the margin. (b) Mesh with resolution that varies linearly from 1000 m grid spacing beyond a radius of 20 km (thick white line) to 5000 m at a radius of 3 km (thin white line). The ice thickness initial condition for the Halfar problem is shown. This mesh requires 1265 cells for the 200-year duration Halfar test case, while a uniform 1000 m resolution mesh requires 2323 cells.

Download

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 first-order 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 variable-resolution 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) mesh-generation 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 variable-resolution meshes. The variable-resolution mesh has about half the cells of the 1000 m uniform-resolution 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 non-physical three-dimensional compensatory heat source necessary for its implementation. While we hope to use the verification of Bueler et al. (2007) in the future, for now we use the EISMINT intercomparison suites to test our implementation of thermal evolution and thermomechanical coupling.

The first phase (Huybrechts et al., 1996) (sometimes called EISMINT1) prescribes evolving ice geometry and temperature, 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 three-dimensional models reported by Huybrechts et al. (1996). MALI achieves similar good agreement for basal homologous temperature at the center of the dome with a value of −13.09 ^∘C compared with $-\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 steady-state solution from experiment A (the warmer air temperature case) as the initial condition to perturbations in the surface air temperature or surface mass balance forcings (experiments B, C, and D). Because these experiments are thermomechanically coupled, they test model ice dynamics and thickness and temperature evolution, as well as their coupling. There is no analytic solution to these experiments, but 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 long-studied feature of the EISMINT2 intercomparison is the cold “spokes” that appear in the basal temperature field of all models in experiment F and, for some models, experiment A (Payne et al., 2000; Saito et al., 2006; Bueler et al., 2007; Brinkerhoff and Johnson, 2013). MALI with shallow ice velocity exhibits cold spokes for experiment F but not experiment A (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 three-dimensional, advective–diffusive description of an enthalpy formulation for energy conservation; it uses the finite-element 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.

Table 3EISMINT2 results for MALI shallow ice model. For each experiment, the model name EISMINT2 refers to the mean and range of models reported in Payne et al. (2000), and we assume that the range reported by Payne et al. (2000) is symmetric about the mean. For experiments B, C, and D reported values are the change from experiment A results. MALI results that lie outside the range of values in Payne et al. (2000) are italicized.

Download Print Version | Download XLSX

Figure 7(a) Basal homologous temperature (K) for EISMINT2 experiment A. (b) Same for experiment F. Figures are plotted following Payne et al. (2000).

Download

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 parallel-sided ice slab with constant thickness and inclination and prescribed ice dynamics decoupled from thermodynamics, resulting in effectively one-dimensional vertical experiments.

Figure 8The result of (a) basal temperature (T_b), (b) basal melt rate (a_b), and (c) basal water thickness (H_w) for enthalpy benchmark experiment A. The green line from 150 to 170 ka in (b) indicates the analytical results (overlapped with model results).

Download

7.3.2 Experiment B

In experiment B, a 200 m thick, 4^∘ downward-inclined 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 steady-state temperate ice layer thickness is achieved (Fig. 9).

Figure 9The vertical distribution of (a) enthalpy (E), (b) temperature (T), and (c) water content (ω) for enthalpy benchmark experiment B. The green lines indicate the analytical results (overlapped with model results).

Download

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 ISMIP-HOM

The Ice Sheet Model Intercomparison Project-Higher Order Models (ISMIP-HOM) is a set of community benchmark experiments for testing higher-order approximations of ice dynamics (Pattyn et al., 2008). Tezaur et al. (2015a) describe results from the Albany-LI velocity solver for ISMIP-HOM experiments A (flow over a bumpy bed) and C (ice stream flow). For all configurations of both tests, Albany-LI results were within 1 standard deviation of the mean of first-order models presented in Pattyn et al. (2008) and showed excellent agreement with the similar first-order 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 stand-alone Albany-LI code it is using.

Figure 10Grid resolution convergence for the MISMIP3d Stnd experiment with (gray squares) and without (black circles) grounding line parameterization.

Download

Figure 11Results of the MISMIP3d P75R and P75S experiments from MALI at increasing grid resolution: (a) 2000 m, (b) 1000 m, (c) 500 m, (d) 250 m. Results for 250 m without grounding line parameterization (e) are also shown for reference. Plots follow conventions of Figs. 5 and 6 in Pattyn et al. (2013). Upper plots show steady-state grounding line positions for steady-state spin-up (black), P75S (red), and P75R (blue) experiments. Lower plots show grounding line position with time for P75R (red) and P75S (blue) at y=0 km (top curves) and y=50 km (bottom curves). 500 and 250 m results are nearly identical. 250 m results without grounding line parameterization are intermediate of those at 1000 and 2000 m resolution with grounding line parameterization.

Download

7.5 MISMIP3d

The Marine Ice Sheet Model Intercomparison Project-3d (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 steady-state ice sheet from a spatially uniform positive surface mass balance, with a prescribed flow rate factor A (no temperature calculation or coupling) and prescribed basal friction for a nonlinear basal friction law. A marine ice sheet forms with an unbuttressed floating ice shelf that terminates at a fixed ice front at the edge of 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 depth-integrated shallow shelf or L1L1/L2L2 approximations, hybrid shallow ice–shallow shelf approximation, or the complete Stokes equations; there were no three-dimensional first-order approximation models included. This relatively simple experiment revealed a number of key features necessary to accurately model even a simple marine ice sheet. Insufficient grid resolution prevented the reversibility of the steady-state grounding line position after experiments P75S and P75R. Reversibility required a grid resolution well below 1 km without a sub-grid 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 steady-state 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 Albany-LI 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 steady-state grounding line in the Stnd experiment for a range of resolutions against our highest-resolution 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 highest-resolution 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 three-dimensional first-order 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 sub-grid 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 highest-resolution 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 three-dimensional first-order 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 (Asay-Davis et al., 2016). These results are included in the MISMIP+ results paper in preparation and not shown here.