MPASSeaice (v1.0.0): seaice dynamics on unstructured Voronoi meshes
 ^{1}Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, USA
 ^{2}Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, CO, USA
 ^{3}Computer, Computational, and Statistical Sciences Division, Los Alamos National Laboratory, Los Alamos, NM, USA
 ^{4}Science Technology Policy Institute, Washington, DC 20006, USA
 ^{1}Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, USA
 ^{2}Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, CO, USA
 ^{3}Computer, Computational, and Statistical Sciences Division, Los Alamos National Laboratory, Los Alamos, NM, USA
 ^{4}Science Technology Policy Institute, Washington, DC 20006, USA
Correspondence: Adrian K. Turner (akt@lanl.gov)
Hide author detailsCorrespondence: Adrian K. Turner (akt@lanl.gov)
We present MPASSeaice, a seaice model which uses the Model for Prediction Across Scales (MPAS) framework and spherical centroidal Voronoi tessellation (SCVT) unstructured meshes. As well as SCVT meshes, MPASSeaice can run on the traditional quadrilateral grids used by seaice models such as CICE. The MPASSeaice velocity solver uses the elastic–viscous–plastic (EVP) rheology and the variational discretization of the internal stress divergence operator used by CICE, but adapted for the polygonal cells of MPAS meshes, or alternatively an integral (“finitevolume”) formulation of the stress divergence operator. An incremental remapping advection scheme is used for mass and tracer transport. We validate these formulations with idealized test cases, both planar and on the sphere. The variational scheme displays lower errors than the finitevolume formulation for the strain rate operator but higher errors for the stress divergence operator. The variational stress divergence operator displays increased errors around the pentagonal cells of a quasiuniform mesh, which is ameliorated with an alternate formulation for the operator. MPASSeaice shares the sophisticated column physics and biogeochemistry of CICE and when used with quadrilateral meshes can reproduce the results of CICE. We have used global simulations with realistic forcing to validate MPASSeaice against similar simulations with CICE and against observations. We find very similar results compared to CICE, with differences explained by minor differences in implementation such as with interpolation between the primary and dual meshes at coastlines. We have assessed the computational performance of the model, which, because it is unstructured, runs with 70 % of the throughput of CICE for a comparison quadrilateral simulation. The SCVT meshes used by MPASSeaice allow removal of equatorial model cells and flexibility in domain decomposition, improving model performance. MPASSeaice is the current seaice component of the Energy Exascale Earth System Model (E3SM).
Sea ice, the frozen surface of the sea at high latitudes, is an important component of the Earth climate system. Rejection of salt during seaice formation helps drive the thermohaline circulation (Killworth, 1983), and its high reflectivity increases planetary albedo and can help drive the polar amplification of climate change through an albedo feedback mechanism (Ingram et al., 1989). Numerical modeling of seaice dynamics and thermodynamics is an important tool in understanding global climate. One of the most popular seaice models currently in use is CICE (Hunke et al., 2015). CICE approximates the seaice cover as a continuous fluid and uses an elastic–viscous–plastic (EVP) rheology to describe the relationship between stress and strain in that fluid (Hunke and Dukowicz, 1997). This rheology adds a numerical elasticity to the viscous–plastic rheology of Hibler (1979) to allow explicit timestepping and improved parallelization scalability of the algorithm. CICE uses a quadrilateral structured mesh and has been used with both displacedpole (Smith et al., 1995) and tripolar (Murray, 1996) grids. The strain rate and internal ice stress divergence operators used by CICE are based on a variational principle (Hunke and Dukowicz, 1997, 2002).
Recently, unstructured mesh climate models have been gaining popularity (e.g., Ringler et al., 2013; Wang et al., 2014). While all the interior vertices of a structured mesh enjoy the same topological connectivity, the vertices of an unstructured mesh have arbitrary topological connectivity (Bern and Plassmann, 2000). This allows unstructured meshes to concentrate their model degrees of freedom in regions of interest, improving computational efficiency, while avoiding the difficulties of open boundary conditions in limiteddomain regional models. The cost of this flexibility in mesh adaptivity, however, is that data access is less straightforward and requires an explicit accounting of the mesh connectivity. Several unstructured seaice models have been developed. Hutchings et al. (2004) implemented and demonstrated a finitevolume, cellcentered discretization. The UnstructuredGrid CICE (UGCICE; Gao et al., 2011) model, which is built on the FVCOM framework (Chen et al., 2009), also uses a finitevolume formulation. In contrast, finiteelement discretizations have been implemented by Lietaer et al. (2008), by Danilov et al. (2015) in the FiniteElement Sea Ice Model (FESIM), and by Mehlmann and Korn (2021) in the seaice component of the Icosahedral Nonhydrostatic Weather and Climate Model (ICON). UGCICE, FESIM, and ICON all use triangular elements in their meshes.
The Model for Prediction Across Scales (MPAS) framework is another recently developed unstructured modeling framework, which uses a spherical centroidal Voronoi tessellation to form a mesh (Ringler et al., 2010). Several climate model components have been built with the MPAS modeling framework, including ocean (MPASO; Ringler et al., 2013), atmosphere (MPASA; Skamarock et al., 2012), and land ice (MALI; Hoffman et al., 2018) models. Here, we describe a new MPAS model for sea ice: MPASSeaice. This model uses a “B” Arakawatype grid (Arakawa and Lamb, 1977) with seaice velocity defined on cell vertices and tracers defined at cell centers. Grid cells in MPAS meshes are polygons with four or more sides, rather than triangles. This allows MPASSeaice to either use the same mesh as structured quadrilateral models such as CICE or match the mesh used by unstructured ocean models such as MPASOcean. Matching the ocean mesh is required for coupling seaice components within some global climate models such as CESM (Danabasoglu et al., 2020) and E3SM (Golaz et al., 2019), and simplifies the ocean–seaice coupling methodology.
We have implemented two discretizations of the momentum equation for MPASSeaice. The first is based on the variational scheme used by CICE (Hunke and Dukowicz, 2002). The second uses the integral (“finitevolume”) form of the relevant operators. We implement this second operator method as a comparison to help investigate sources of error in the standard variational scheme. Within the variational scheme we also investigate several enhancements with the aim of reducing errors associated with the variational scheme. MPASSeaice uses the same EVP rheology and column physics as CICE. Seaice tracer transport uses an incremental remapping scheme (Lipscomb and Hunke, 2004; Lipscomb and Ringler, 2005) modified for the MPAS mesh. In Sect. 2, we describe the modeling approach used in MPASSeaice, focusing on the solution of the seaice momentum equation and tracer transport. In Sect. 3 we validate this new model, with both idealized test cases and global simulations. In Sect. 4 we consider the computational performance of the model and conclude in Sect. 5. We focus in this paper on simulations with quasiuniform global meshes; variableresolution meshes will be considered in later publications.
2.1 The MPAS framework
The MPAS mesh uses a spherical centroidal Voronoi tessellation (SCVT), as described in Ringler et al. (2010), and consists of a primary mesh of the Voronoi cells tessellated on the sphere and a dual triangular Delaunay mesh, formed from joining the Voronoi cell centers into a triangulation. Figure 1 shows a schematic of part of a MPAS mesh with the primary mesh shown as a solid line and the dual mesh shown as a dashed line. The primary mesh cells have their Voronoi generating points coincident with the centroid of the cell. The mesh consists of three types of points arranged on the sphere: the Voronoi cell center points, vertex points of the Voronoi cell, and edge points at the midpoint of the Voronoi cell edges (see Fig. 1). The mesh is constructed so that a line joining neighboring cell centers is perpendicular to the edge that line passes through, with the edge equidistant from the two cell centers. On a typical quasiuniform MPAS mesh, the majority of the cells are hexagons, but at least 12 pentagons are needed to complete the tessellation when cells cover the entire sphere. Fewer pentagons would be required for an ocean–seaice mesh where only part of the sphere is covered in cells. In general, the cells are not regular polygons and may consist of polygons with edge numbers greater than or equal to four. As well as quasiuniform grids, meshes can be generated with regions of enhanced resolution, allowing computational effort to be focused in regions of interest. The MPAS mesh standard can also represent the quadrilateral meshes used by CICE, where four instead of three edges meet at each vertex. With these quadrilateral meshes the dual mesh is also quadrilateral.
2.2 Velocity solver
MPASSeaice uses a “B” Arakawatype grid (Arakawa and Lamb, 1977) with both components of velocity defined at cell vertices and with seaice concentration, volume, and other tracers defined at cell centers (see Fig. 1). When using CICElike quadrilateral meshes, the velocity solver algorithm of MPASSeaice reduces to that of CICE, allowing CICE and MPASSeaice to use identical test cases and supporting rapid testing and development of MPASSeaice.
In CICE the velocity components are aligned with the quadrilateral mesh. This is not possible, in general, with MPASSeaice since a SCVT MPAS mesh does not have edges with perpendicular directions as in a quadrilateral mesh. Instead, the velocity components at a given MPAS vertex are defined as eastwards (u) and northwards (v), irrespective of the orientation of edges joining that vertex. Such a definition, however, would result in a convergence of v components of velocity at the geographic poles and strong metric terms in the velocity solution. Consequently, we rotate the coordinate system so that the pole of u and v lies on the geographical Equator at 0^{∘} longitude. The polar regions then have the smallest metric effects on the globe.
To prognose seaice velocity we solve the same seaice momentum equation as CICE (Hibler, 1979; Hunke and Dukowicz, 1997):
Here m is the mass of snow and ice per unit area; u is the seaice velocity; σ is the ice internal stress tensor; τ_{a} and τ_{w} are the horizontal stresses due to atmospheric winds and ocean currents, respectively; $\widehat{\mathit{k}}$ is the unit vector normal to the Earth surface; f is the Coriolis parameter; g is the acceleration due to gravity; and H_{o} is the ocean surface height. The secondtolast term represents the Coriolis force, and the last term represents the force due to the ocean surface tilt. Only the internal stress divergence and ocean surface tilt terms depend on horizontal differential operators. During coupled simulations the ocean model provides the ocean surface tilt term, whereas in noncoupled simulations we assume that the ocean currents are in geostrophic balance so that
where u_{o} is the geostrophic component of the ocean surface velocity. Consequently, only the internal stress divergence depends on the properties of the horizontal grid, and only adaptations to this stress term are required to adapt the velocity solver of CICE to MPAS meshes. The other terms in the momentum equation are solved in an identical way to CICE.
Determination of the divergence of the internal stress can be broken down into three stages:

Determine the strain rate tensor from the velocity field.

Determine the stress tensor at a point, through a constitutive relation, from the strain rate tensor at that point.

Calculate the divergence of this stress tensor.
As in CICE we use an elastic–viscous–plastic (EVP) rheology (Hunke and Dukowicz, 1997) for the constitutive relation. This step does not depend on the details of the horizontal mesh, and we use the same formulation as CICE. We develop two schemes to calculate the strain rate tensor and the divergence of internal stress on MPAS meshes. A variational scheme is based on that used in CICE (Hunke and Dukowicz, 2002), whereas a finitevolume scheme uses the lineintegral forms of the symmetric gradient and divergence operators. These schemes are described in the following sections.
2.2.1 Variational scheme
We develop a variational scheme for calculating the stress divergence based on that of Hunke and Dukowicz (2002) but adapted for arbitrarily shaped and sided convex polygons. The principal change needed to adapt Hunke and Dukowicz (2002) to polygonal cells is a generalization of the basis functions from bilinear to a basis compatible with polygonal cells. The variational scheme is based on the fact that over the entire domain, Ω, the total work done by the internal stress is equal to the dissipation of mechanical energy:
Here $\dot{\mathit{\u03f5}}$ is the strain rate tensor, and the integrals are area integrals over the whole model domain. For simplicity we ignore boundary effects, which would add additional terms to the righthand side of this equation. The stress divergence operator is derived from the functional
where
In this functional the stress, σ, is taken as a separate parameter to the velocity, u. The derivation proceeds by determining a discretized version of the I[u,σ] functional and taking the variation of this functional with respect to one of the discretized velocities. The discretized functional is given by
where the functional has the discretized velocity components, u_{i} and v_{i}, and the discretized stresses, σ_{j}, as parameters; u is the eastwards component of u, while v is the northwards, and the two components of the stress divergence have been written as ${F}_{u}=(\mathbf{\nabla}\cdot \mathit{\sigma}{)}_{u}$ and ${F}_{v}=(\mathbf{\nabla}\cdot \mathit{\sigma}{)}_{v}$. The work done in the whole domain has been split into a sum over the contributions to the work done in each cell ($i=\mathrm{1},\mathrm{\dots},{n}_{d}$) on the dual Delaunay mesh. The contribution to the work done in each dual cell is an area integral over each dual cell. These dual cells consist of either triangles (for SCVT meshes), or quadrilaterals (for quadrilateral meshes) surrounding a single primary mesh vertex point where the discretized velocity is defined. The dissipation of mechanical energy, D, has been written as a function of the discretized velocity components and discretized stresses. The simplest assumption for the velocity and stress divergence components for the work part of the discretized functional is that these quantities are constant within the dual cell. This is assumed by Hunke and Dukowicz (2002), but we later improve this assumption. With this initial assumption the discretized functional becomes
Taking the variation in this functional with respect to a single discretized velocity component, u_{j}, we get the following Euler–Lagrange equation
where
and ${\dot{u}}_{j}=\frac{\partial {u}_{j}}{\partial x}$ within the Euler–Lagrange equation. In the following the strain rate in D will be given directly in terms of u_{i} and v_{i} rather than ${\dot{u}}_{i}$ and ${\dot{v}}_{i}$, so the second term of the Euler–Lagrange equation is zero. Then
and consequently the u component of the stress divergence for vertex j is given by
F_{v} is obtained in a similar way by taking the variation in I with respect to v_{j}. The dissipation of mechanical energy, D, can be split into three terms:
with
We calculate the contribution to F_{u} and F_{v} from D_{1}. Similar contributions come from D_{2} and D_{3}. Using the expression for ${\dot{\mathit{\u03f5}}}_{\mathrm{11}}$ in terms of the velocity components and latitude θ, D_{1} becomes
where x and y are locally Cartesian coordinates, with x in the rotated eastwards direction and y in the rotated northwards direction, and r is the radius of the Earth. The second term in $\dot{\mathit{\u03f5}}$ accounts for the metric effects of the curved domain (Batchelor, 1967). The integral can be broken up into a sum over the n_{p} cells in the primary mesh:
where the integral is over the interior area of the kth cell. To perform this integral we use a set of Lagrangian basis functions, 𝒲_{l}, to represent functions within a cell of the primary mesh. These basis functions are such that if a function, ψ, has a value of ψ_{l} at vertex l of a cell, then the value of the function at a position (x,y) within the cell can be approximated as
where the sum is over the n_{v} vertices of the cell in the primary mesh. Necessary properties for these basis functions are that
across the cell and that
We provide two options for the choice of basis functions, 𝒲_{l}: Wachspress basis functions and piecewise linear (PWL) basis functions. Both basis functions have a value of one on vertex l and zero on the other vertices of a cell and are linear on the cell boundaries. The Wachspress basis function for the ith vertex of a polygon with n sides is given by (Dasgupta, 2003)
where
and
l_{j}(x,y) is the line equation for the j polygon edge such that
where a_{j} and b_{j} are defined by the condition that ${l}_{j}(x,y)=\mathrm{0}$ along the jth edge of the polygon. When written out, 𝒲_{i} becomes a rational polynomial of the form
where 𝒫^{(m)}(x,y) is an mdegree polynomial in x,y. The integrals of the Wachspress basis function within a cell are performed using the eighthorder quadrature rules of Dunavant (1985). For quadrilateral meshes the Wachspress basis functions reduce to the bilinear basis functions used in CICE. For the four vertices of the quadrilateral cells, these are given by
where ξ_{1} and ξ_{2} are the transformed unit square coordinates of the cell (Hunke, 2001).
PWL basis functions divide the polygonal cell into subtriangles and use a linear basis within each subtriangle (Bailey et al., 2008). To divide the polygonal cell into subtriangles, a point is chosen within the cell and subtriangles formed using this point and two adjacent vertices. The central point in the cell, x_{c}, is chosen as
where the sum is over the n_{v} vertices of the cell each with position x_{i}. The simplest choice for the α_{i} is to set them all equal to the inverse of the number of cell vertices, $\mathrm{1}/{n}_{v}$. Example basis functions for the Wachspress and PWL basis functions are shown in Fig. 2.
Using those basis functions to expand σ_{11} (with basis index l), u, and v (with basis index m), Eq. (15) can be written as
where the derivative with respect to x has been taken inside the summation. Rearranging, this becomes
In moving the integral, we have assumed that θ, the latitude, is constant in the cell. The terms involving integrals are now only a function of the geometry of the mesh and can be calculated once for each cell during the initialization phase of the model run. Defining
and
we have
Taking the variation with respect to a discretized velocity component at a particular vertex point, j, as in Eq. (11), now gives us the contribution from D_{1} to the components of the stress divergence tensor at that velocity point:
Only cells that border the vertex point j contribute to the k sum over cells. The total stress divergence at the point j is then the sum of the contributions from D_{1}, D_{2}, and D_{3}:
All that remains is to determine the stress for each cell at its vertices. As in the formulation in CICE, each cell has its own stress value at its vertices, so each vertex has several values of the stress, each corresponding to a different surrounding cell. The stresses are calculated from the strain rate tensor at each vertex using the constitutive relation (see Sect. 2a of Hunke and Dukowicz, 2002). Including metric effects (Batchelor, 1967), the strain rate tensor is given by
The strain rate tensor at cell vertex l is then given by
The derivatives of the basis functions are taken at cell vertex l. For the variational method precomputed values for the variables ${\mathcal{S}}_{lm}^{x}$, ${\mathcal{S}}_{lm}^{y}$, 𝒯_{lm}, ${\left(\frac{\partial {\mathcal{W}}_{m}}{\partial x}\right}_{l}$, and ${\left(\right)}_{\frac{\partial {\mathcal{W}}_{m}}{\partial y}}l$ must be stored, resulting in a total of $\mathrm{5}{n}_{v}^{\mathrm{2}}$ values stored per cell.
One issue with the above derivation is the ambiguity of defining the correct dual cell area, A_{j}, for Eq. (11) since the cell center positions do not appear anywhere else in the formulation. In Sect. 3.1.2, in a unit sphere test case with analytical input stress fields, it is evident that using the default MPAS dual cell area for A_{j} results in large errors in the calculated stress divergence field around the 12 pentagonal cells found on a quasiuniform SCVT mesh. An alternative to using the MPAS dual cell area and assuming constant stress and velocity within a dual cell for the work equation is to assume that the velocity and divergence are given by the same basis functions as used for the dissipation of mechanical energy. Then the work done over the domain
Taking the variation in I with respect to the discretized velocity component u_{j}
where the sum over i cells is reduced to those surrounding the j vertex. This generates an undesirable implicit problem necessitating a global matrix inversion, which can be avoided if the stress divergence varies slowly spatially. This is the equivalent of lumping the mass matrix in the finite element method (Reddy, 1993). Then
where the sum over m is over the vertices in each cell surrounding the j vertex. This suggests an alternative to A_{j} given by
We find that approximating A_{j} with ${A}_{j}^{\prime}$ reduces the error in the stress divergence operator as compared to using the standard dual mesh areas and requires no additional calculations of integrals at initialization.
2.2.2 Finitevolume scheme
An alternative method of deriving operators is with an integral method where the divergence theorem is used to equate the integral form of the operator to a flux integrated around a closed loop. One potential advantage of the finitevolume scheme is that it solves the conservative form of the momentum equation and can handle nonsmooth solution features (such as sharp fronts) consistently.
For the finitevolume scheme we use line integrals around cells in the primary and dual meshes to calculate the strain rate tensor and the stress divergence, respectively. To determine the strain rate tensor we start from the generalized divergence theorem
where n is the normal to the closed surface S of domain Ω, and ⊗ is the tensor product. The symmetric version of this operator is then obtained as
The strain rate at a point is then obtained from the limit
where the integral is around a closed loop, S, with area A and normal vector n, and v is the seaice velocity. To determine the strain rate tensor at the centers of the primary mesh, we take this integral around the edges of the cells in the primary mesh. First the cell is projected onto a flat tangent plane perpendicular to the vector joining the center of the sphere to the cell center. We take the seaice velocity at a cell edge as the average of the values on the two vertices forming that edge projected onto the tangent plane:
Here, A is the area of the primary cell, the summation is over the n_{e} edges of the primary cell, n_{i} is the normal vector to the edge i that lies in the tangent plane, v_{i} is the edge velocity, and l_{i} is the length of edge i. We use the tangential projection of the velocity and account for metric terms separately. The full strain rate tensor including these metric terms is (Batchelor, 1967)
where the prime symbol signifies a strain rate without metric terms. The stress, which is determined from the strain rate tensor using the constitutive relation, is now defined on cell centers. To find its divergence we use the divergence theorem
or
for the stress divergence at a point. The divergence of internal stress is determined at primary cell vertices (where the velocity is defined and momentum equation solved) by performing a sum around the edges of the dual mesh on a tangent projected plane, tangential to the primary cell vertex. The vertices of the dual mesh are the cell centers of the primary mesh where the strain rate has been determined. The stress divergence at primary cell vertices is then given by
where A_{d} is the area of the dual mesh cell, the sum is over the n_{c} vertices of the dual mesh, l_{i} is the length of the i edge of the dual mesh, and n_{i} is a normal vector to the i edge on the projected plane. Metric terms in the divergence of a secondorder tensor, like the stress tensor, have two contributions: the first comes from the varying grid cell size and the second from the varying directions of the coordinate basis vectors (Malvern, 1969). Equation (45) accounts for the first of these effects, and in order to account for the second effect the stress divergence becomes
where the components of σ are approximated as the average of the values on the dual mesh vertices.
2.3 Transport
To transport seaice fractional area and various tracers, MPASSeaice uses an incremental remapping (IR) algorithm similar to that described by Dukowicz and Baumgardner (2000), Lipscomb and Hunke (2004), and Lipscomb and Ringler (2005). The Lipscomb and Hunke (2004) scheme was designed for structured quadrilateral meshes and is implemented in CICE (Hunke et al., 2015). The Lipscomb and Ringler (2005) scheme was implemented on a structured SCVT global mesh consisting of quasiregular hexagons and 12 pentagons.
For MPASSeaice the IR scheme was generalized to work on either the standard MPAS mesh (hexagons and other ngons of varying sizes, with a vertex degree of 3 as in Lipscomb and Ringler, 2005) or a quadrilateral mesh (with a vertex degree of 4 as in Lipscomb and Hunke, 2004, and CICE). Since the CICE IR transport code assumes a structured mesh, but MPAS meshes are unstructured, the IR scheme had to be rewritten from scratch. Most of the MPASSeaice IR code is now meshagnostic, but a small amount of code is specific to quad meshes as noted below.
Here we review the conceptual framework of incremental remapping as in Hunke et al. (2015) and describe features specific to the MPASSeaice implementation. IR is designed to solve equations of the form
where $\mathit{u}=(u,v)$ is the horizontal velocity; m is mass or a masslike field (such as density or fractional seaice concentration); and T_{1}, T_{2}, and T_{3} are tracers. These equations describe conservation of quantities such as mass and internal energy under horizontal transport. Sources and sinks of mass and tracers (e.g., ice growth and melting) are treated separately from transport.
In MPASSeaice, the fractional ice area in each thickness category is a masslike field whose transport is described by Eq. (47). (Henceforth, “area” refers to fractional ice area unless stated otherwise.) Ice and snow thickness, among other fields, are type 1 tracers obeying equations of the form of Eq. (48), and the ice and snow enthalpy in each vertical layer are type 2 tracers obeying equations like Eq. (49), with ice or snow thickness as their parent tracer. When run with advanced options (e.g., active melt ponds and biogeochemistry (BGC)), MPASSeaice advects tracers up to type 3. Thus, the masslike field is the “parent field” for type 1 tracers, type 1 tracers are parents of type 2, and type 2 tracers are parents of type 3.
Incremental remapping has several desirable properties for seaice modeling:

It is conservative to within machine roundoff.

It preserves tracer monotonicity. That is, transport produces no new local extrema in fields like ice thickness or internal energy.

The reconstructed mass and tracer fields vary linearly in x and y. This means that remapping is secondorder accurate in space, except where gradients are limited locally to preserve monotonicity.

There are economies of scale. Transporting a single field is fairly expensive, but additional tracers have a low marginal cost, especially when all tracers are transported with a single velocity field as in CICE and MPASSeaice.
The time step is limited by the requirement that trajectories projected backward from vertices are confined to the cells sharing the vertex (i.e., three cells for the standard MPAS mesh and four for the quad mesh). This is what is meant by incremental as opposed to general remapping. This requirement leads to a Courant–Friedrichs–Lewy (CFL)like condition,
where Δx is the grid spacing, and Δt is the time step. For highly divergent velocity fields, the maximum time step may have to be reduced by a factor of 2 to ensure that trajectories do not cross.
The IR algorithm consists of the following steps:

Given mean values of the ice area and tracer fields in each grid cell and thickness category, construct linear approximations of these fields. Limit the field gradients to preserve monotonicity.

Given ice velocities at grid cell vertices, identify departure regions for the transport across each cell edge. Divide these departure regions into triangles, and compute the coordinates of the triangle vertices.

Integrate the area and tracer fields over the departure triangles to obtain the area, volume, and other conserved quantities transported across each cell edge.

Given these transports, update the area and tracers.
Since all fields are transported by the same velocity field, the second step is done only once per time step. The other steps are repeated for each field.
With advanced physics and BGC options, MPASSeaice can be configured to include up to ∼ 40 tracer fields, each of which is advected in every thickness category and many of which are defined in each vertical ice or snow layer. In order to accommodate different tracer combinations and make it easy to add new tracers, the tracer fields are organized in a linked list that depends on which physics and BGC packages are active. The list is arranged with fractional ice area first, followed by the type 1 tracers, type 2 tracers, and finally type 3 tracers. In this way, values computed for parent tracers are always available when needed for computations involving child tracers.
The MPASSeaice version of the incremental remapping transport scheme has several advantages over the one in CICE. First, as already mentioned, the scheme has been generalized to work for SCVT as well as quadrilateral meshes. Second, the MPASSeaice scheme treats T_{3} tracers more accurately by using more accurate integration formulas (see Eq. 60).
We next describe the IR algorithm in detail, pointing out features that are new in MPASSeaice.
2.3.1 Reconstructing area and tracer fields
The fractional ice area and all tracers are reconstructed in each grid cell (quadrilaterals, hexagons, or other ngons) as functions of $\mathit{r}=(x,y)$ in a cellbased coordinate system. On spherical grids, r lies in a local plane that is tangent to the sphere at the Voronoi cell center. The state variable for ice area, denoted as $\stackrel{\mathrm{\u203e}}{a}$, should be recovered as the mean value when integrated over the cell:
where A_{C} is the grid cell area. Equation (52) is satisfied if a(r) has the form
where ∇a is a cellcentered gradient, α_{a} is a coefficient between 0 and 1 that enforces monotonicity, and $\stackrel{\mathrm{\u203e}}{\mathit{r}}$ is the cell centroid:
Similarly, tracer means should be recovered when integrated over a cell:
These equations are satisfied when the tracers are reconstructed as
where the tracer barycenter coordinates ${\stackrel{\mathrm{\u0303}}{\mathit{r}}}_{n}$ are given by
The integrals in Eq. (57) can be evaluated by applying quadrature rules for linear, quadratic, and cubic polynomials as described in Sect. 2.3.3.
Monotonicity is enforced by van Leer limiting (van Leer, 1979). The reconstructed area and tracers are evaluated at cell vertices, and the coefficients α are reduced as needed so that the reconstructed values lie within the range of the mean values in the cell and its neighbors. When α=1, the reconstruction is secondorder accurate in space. When α=0, the reconstruction reduces locally to firstorder.
2.3.2 Locating departure triangles
The next step is to identify the departure region associated with fluxes across each cell edge and to divide the departure region into triangles. Figure 4a illustrates the geometry for the standard MPAS mesh. The edge has vertices V_{1} and V_{2}. Each edge is oriented such that one adjacent cell (C_{1}) is defined to lie in the left halfplane and the other (C_{2}) in the right halfplane. The departure points D_{1} and D_{2} are found by projecting velocities backward from V_{1} and V_{2}. The shaded departure region is a quadrilateral containing all the ice transported across the edge in one time step. In addition to C_{1} and C_{2}, the departure region can include side cells C_{3} and C_{4}. The side cells share edges E_{1} to E_{4} and vertices V_{3} to V_{6} with the central cells C_{1} and C_{2}.
The edges and vertices in Fig. 4a are defined in a coordinate system lying in the local tangent plane at the midpoint of the main edge, halfway between V_{1} and V_{2}. These coordinates are precomputed at initialization. During each time step, departure triangles are found by locating D_{1} and D_{2} in this coordinate system and then looping through the edges to identify any intersections of line segment $\stackrel{\mathrm{\u203e}}{{D}_{\mathrm{12}}}$ (i.e., the segment joining D_{1} and D_{2}) with the various edges. If $\stackrel{\mathrm{\u203e}}{{D}_{\mathrm{12}}}$ intersects the main edge, then the departure region consists of two triangles (one each in C_{1} and C_{2}) rather than a quadrilateral (as shown in Fig. 4b). If $\stackrel{\mathrm{\u203e}}{{D}_{\mathrm{12}}}$ intersects any of edges E_{1} to E_{4}, the departure region includes triangles in side cells.
Each departure triangle lies in a single grid cell, and there are at most four such triangles. There are two triangles in the central cells (either one each in C_{1} and C_{2} or a quadrilateral that can be split into two triangles) and up to two triangles in side cells. The triangle vertices are a combination of cell vertices (V_{1} and V_{2}), departure points (D_{1} and D_{2}), and intersection points (points where $\stackrel{\mathrm{\u203e}}{{D}_{\mathrm{12}}}$ crosses an edge).
Figure 4b shows the geometry for a quadrilateral mesh. In this figure the departure region consists of two triangles, but it could also be a quadrilateral as in Fig. 4a. For the quad mesh there are two additional side cells (C_{5} and C_{6}), edges (E_{5} and E_{6}), and vertices (V_{7} and V_{8}). The search algorithm is designed such that the code used to find departure triangles for the standard mesh is also applied to the quad mesh. For quad meshes only, there is additional logic to find intersection points and triangles associated with the extra edges and cells. This is the only meshspecific code in the runtime IR code. For the quad mesh there are at most six departure triangles: two in the central cells and one in each of the four side cells. If the edges meet at right angles as shown in the figure, the maximum is five triangles, but this is not a mesh requirement.
Once triangle vertices have been found in edgebased coordinates, they are transformed to cellbased coordinates, i.e., coordinates in the local tangent plane of the cell containing each triangle. (Coefficients for these transformations are computed at initialization.) Triangle areas are computed as
where the (x_{i},y_{i}) are the triangle vertices.
2.3.3 Integrating the transport
Next, ice area and areatracer products are integrated in each triangle. The integrals have the form of Eq. (52) for area and Eq. (55) for tracers. Since each field is a linear function of (x,y) as in Eqs. (53) and (56), the areatracer products are quadratic, cubic, and quartic polynomials, respectively, for tracers of type 1, 2, and 3.
The integrals can be evaluated exactly by summing over values at quadrature points in each triangle. Polynomials of quadratic or lower order are integrated using the formula
The quadrature points are located at ${{\mathit{x}}^{\mathbf{\prime}}}_{i}=({\mathit{x}}_{\mathrm{0}}+{\mathit{x}}_{i})/\mathrm{2}$, where x_{0} is the triangle midpoint, and x_{i} are the three vertices. The products involving type 2 and type 3 tracers are cubic and quadratic polynomials, which can be evaluated using a similar formula with six quadrature points:
where x_{1i} and x_{2i} are two sets of three quadrature points, arranged symmetrically on the three medians of the triangle, and w_{1} and w_{2} are weighting factors. Coefficients and weighting factors for these and other symmetric quadrature rules for triangles were computed by Dunavant (1985). These integrals are computed for each triangle and summed over edges to give fluxes of ice area and areatracer products across each edge.
2.3.4 Updating area and tracer fields
The area transported across edge k for a given cell can be denoted as Δa_{k} and the areatracer products as Δ(aT_{1})_{k}, Δ(aT_{1}T_{2})_{k}, and Δ(aT_{1}T_{2}T_{3})_{k}. The new ice area at time n+1 is given by
where the sum is taken over cell edges k, with a positive sign denoting transport into a cell and a negative sign denoting outward transport. The new tracers are given by
Dukowicz and Baumgardner (2000) showed that Eq. (62) satisfies tracer monotonicity since the newtime tracer values are areaweighted averages of oldtime values.
2.4 Column physics
CICE has sophisticated vertical physics and biogeochemical schemes, which include vertical thermodynamics schemes (Bitz and Lipscomb, 1999; Turner et al., 2013; Turner and Hunke, 2015), several meltpond parameterizations (Flocco et al., 2010; Holland et al., 2012; Hunke et al., 2013), a deltaEddington radiation scheme (Briegleb and Light, 2007; Holland et al., 2012), schemes for transport in thickness space (Lipscomb, 2001), representations of mechanical redistribution (Lipscomb et al., 2007), and seaice BGC (Jeffery and Hunke, 2014; Jeffery et al., 2016, 2020). To include these developments in MPASSeaice, the column physics and BGC in CICE have been extracted into a separate library, which performs column calculations on individual grid cells with no knowledge of the details of the horizontal mesh used by the host model. This column package, now called Icepack, is used by the most recent version of CICE (Hunke et al., 2018). MPASSeaice uses a version of the column package that was forked from CICE version 5.1.2.
3.1 Velocity solver
To validate our implementation of the MPASSeaice velocity solver we investigate several idealized test cases.
3.1.1 Operators on planar meshes
In the first test case, we define analytical test fields on a unit square planar mesh with regularly shaped cells and determine the error generated by the strain and stress divergence operators when acting on these fields. Since the strain rate operators consist of a combination of spatial derivatives in the x and y directions, we examine the error generated by these spatial derivative operators instead of directly by the strain rate operator. We examine meshes with both square and hexagonal cells. The test analytical field for which spatial derivatives are calculated is given by
with A=2.56. This field gives sufficient variation to provide an adequate test of the operators, while the symmetry of f(x,y) between the x and y directions allows an accurate comparison between the $\frac{\partial}{\partial x}$ and $\frac{\partial}{\partial y}$ operators. Figure 5 shows the scaling of the L_{2} error norm for the spatial derivatives against grid resolution. For the spatial derivatives examined in this section and the strain rates examined in Sect. 3.1.2, we calculate the L_{2} norm with the following methods. For the variational operators we integrate the square of the error across the domain using the variational basis functions defined in Sect. 2.2.1 for the quantity of interest. The L_{2} norms for the variational derivatives and strains are given by
where the sum over i is a sum over cells in the region of interest, and the area integral is over cell i, performed by splitting the polygon into subtriangles and using the eighthorder integration rules of Dunavant (1985). The modeled quantity of interest is determined within the interior of the cell from the basis functions, 𝒲_{l}, and quantity of interest, f_{il}, on vertex l of cell i. For strains calculated with the Wachspress basis function we perform the integration with the Wachspress basis functions and likewise for the PWL basis functions; $\widehat{f}(x,y)$ represents the analytical values of the field of interest within the cell. For the finitevolume strain operators we perform line integrals of the square of the error around the dual mesh cell surrounding primary cell points. The L_{2} norms for the finitevolume derivatives and strains are given by
where the sum i is over the primary cells, and the sum j is over the edges of the dual cell surrounding cell i. Coordinate χ signifies the position along edge j, with f_{j}(χ) calculated from a linear interpolation from the strain values at each end of the edge; ${\widehat{f}}_{j}\left(\mathit{\chi}\right)$ are the analytical values of the field of interest along the edge. The line integrals are performed with seventhorder Gauss–Lobatto quadrature rules. These formulas emulate how the stress values, derived from the strain values, are used in their respective stress divergence operators. Each doubling of resolution should reduce the L_{2} error norm by a factor of 2 for firstorder accurate schemes and a factor of 4 for secondorder accurate schemes. Figure 5 also shows idealized first and secondorder scaling gradients as dotted lines. Evident from the figure is that the spatial derivative operators for the variational methods (Wachspress and PWL) are only firstorder accurate for both square and hexagonal meshes, while the finitevolume derivative operators are secondorder accurate with much lower errors than the variational methods. This is understandable since the variational derivative operators only use velocity values from vertices on the same cell as the vertex that the derivative is being calculated for. These velocity values can then only occupy one halfplane with respect to the derivative point, effectively creating a onesided stencil for the operator. Velocity values for the finitevolume operator, however, surround the derivative point since here the derivative point is at the center not the side of the cell. This effectively creates a stencil surrounding the derivative location. Small differences exist between the Wachspress and PWL variational basis functions as well. Differences between errors generated for the x and y spatial derivatives are evident for the hexagonal mesh. These differences are caused by a difference in spatial symmetry of the mesh in these two directions for the hexagonal mesh, whereas the square mesh has the same spatial symmetry in both directions. These results are confirmed by an error analysis with a Taylor series expansion of the methods. One possible way to improve the accuracy of the variational operators is to average for each surrounding cell the different onesided values of the derivatives calculated for a single vertex to create a multisided stenciled operator from the onesided variational ones. Figure 5 shows the error scaling for this averaging method for the Wachspress (“Wachs. Avg”) and PWL (“PWL Avg”) basis functions. Secondorder convergence is achieved in the y direction with this averaging, but interestingly, only averaging the PWL basis function results in secondorder convergence in the x direction, while the Wachspress basis function shows only firstorder convergence with much larger errors in the x direction. This difference between the two directions is explained by the different geometric symmetry found in the two directions, with, for this test case, the x direction parallel to a line from a cell center to a vertex and the y direction parallel to a line from a cell center to edge center. The derivatives at cell centers calculated by the finitevolume method can also be averaged to the cell vertices to create another operator (“FV Avg” in Fig. 5), creating a secondorder accurate method, although this averaging increases the error relative to the original finitevolume scheme.
Next, with the same meshes we examine the error scaling of the stress divergence operator. The same function as above, f(x,y), is used for the u and v components of velocity, and the analytical strain rates are calculated from these velocity components. The analytical internal ice stresses and stress divergence are calculated from this strain rate assuming a linear constitutive relation of the form ${\mathit{\sigma}}_{ij}={\dot{\mathit{\u03f5}}}_{ij}$. The analytical internal stresses are used as input to the stress divergence operators, and the output is compared to the analytical stress divergence. Figure 6 shows the scaling of the L_{2} error norm of the stress divergence operators with grid resolution. For the L_{2} error norm calculated for the stress divergence operators in this section and in Sect. 3.1.2, we use
where the sum is over either grid cells or vertices in the region of interest; A_{i} is the area of either the primary cell or the dual cell surrounding the vertex; and f_{i} and ${\widehat{f}}_{i}$ are the model and analytical values of the field of interest, respectively. Both the variational and finitevolume methods show secondorder convergence of errors for the square cell mesh, with the finitevolume scheme showing significantly lower errors. With the PWL basis functions the variational method shows secondorder convergence, while with Wachspress basis functions the variational method shows a varying order of convergence with nearsecondorder convergence at low resolution but with the order decreasing to first with higher resolutions. This suggests there is a small source of firstorder error when using the Wachspress basis functions. For the hexagonal cell mesh, the finitevolume method shows only firstorder convergence with significantly higher error than the variational method. This is because the line integral around vertices for the finitevolume method involves an integral around a triangle. Unlike the integral around a square or hexagon, there are no opposite sides to a triangle. For the finitevolume integrals around squares or hexagons, opposite edges of the polygon lead to cancellation of firstorder error, which does not occur for integrals around triangles. Smaller differences between the x and y directions are visible than with the spatial derivative operators examined above. Again, these results are confirmed with a Taylor series expansion of the methods.
In summary, for regular planar meshes, for the gradient operators, the finitevolume and averaged PWL and averaged finitevolume schemes show a higher order of error convergence and lower absolute errors than the variational schemes, while the averaged Wachspress scheme shows a higher order of convergence than the unaveraged Wachspress scheme, except for hexagonal cell meshes when the x derivative is being calculated. Conversely, for the stress divergence operator and hexagonal cell meshes, the variational schemes show lower absolute errors and better error convergence than the finitevolume scheme, while for square cell meshes the order of convergence is the same between the variational and finitevolume schemes, with the finitevolume scheme producing lower absolute errors. Within the variational scheme, use of the PWL basis functions displays more consistent error convergence characteristics than the Wachspress scheme.
3.1.2 Operators on a unit sphere
Having examined the spatial operators on planar meshes, we now examine the effect of metric terms introduced by using the operators on a sphere. To do this we assume an analytical velocity field on a unit sphere and derive analytical strain rate and stress divergence fields from those velocity fields (again assuming a linear constitutive relation of the form ${\mathit{\sigma}}_{ij}={\dot{\mathit{\u03f5}}}_{ij}$). We use spherical harmonic functions, Y, for the analytical velocity fields with $u(\mathit{\theta},\mathit{\varphi})={Y}_{l=\mathrm{5}}^{m=\mathrm{3}}(\frac{\mathit{\pi}}{\mathrm{2}}\mathit{\theta},\mathit{\varphi})$ and $v(\mathit{\theta},\mathit{\varphi})={Y}_{l=\mathrm{4}}^{m=\mathrm{2}}(\frac{\mathit{\pi}}{\mathrm{2}}\mathit{\theta},\mathit{\varphi})$, where u, v, the latitude (θ), and the longitude (ϕ) are on the rotated geographical mesh (see Sect. 2.1). This choice of spherical harmonic, with different values of m and l for u and v, produces a sufficiently varying velocity field to test the strain rate and stress divergence operators. Figure 7a and b show these analytical velocity fields, while Fig. 7c and d and e, f, and g show the derived analytical fields for stress divergence and strain rate, respectively.
Errors generated for the strain rate component ${\dot{\mathit{\u03f5}}}_{\mathrm{11}}$ for this test case are displayed in Fig. 7h, l, and p. These show that the error is significantly lower for the finitevolume scheme compared to the two variational schemes. They also show that the error in the variational scheme consists of alternating signed error values within a cell, while the error for the finitevolume scheme is enhanced around the 12 pentagonal cells found in the quasiuniform SCVT mesh. These features are clearer in Fig. 7i, m, and q, which show the detail around one of the pentagonal cells. A histogram of error values for the various strain rate operators is presented in Fig. 8a. Evident are the much larger error values for the variational schemes compared to the finitevolume scheme. The scaling of L_{2} error norm with grid resolution is shown in Fig. 9a. L_{2} error norm is calculated for regions of the grid with latitude $\left\mathit{\theta}\right>\mathrm{20}{}^{\circ}$ (so that the poles of the rotated mesh are excluded) for four different grid resolutions. Evident from the scaling figure is that the variational schemes display firstorder convergence, and the finitevolume scheme shows secondorder convergence and much lower error levels. Also shown in this figure are results for the operator averaging methods described in Sect. 3.1.1. Here the PWL and finitevolume operator averaging show similar error scaling characteristics with closetosecondorder convergence. Also apparent is that a lower improvement in error for the Wachspress variational scheme is achieved with averaging, where convergence is closer to first order. Figure 9c shows the scaling of the L_{∞} norm, defined as ${max}_{i}{f}_{i}\widehat{{f}_{i}}$, which shows firstorder convergence for all strain rate operator methods. As for the hexagonal cell planar mesh, for both the L_{2} and L_{∞} norms, averaging the finitevolume scheme increases the error.
Next we examine the stress divergence operators on the unit sphere. Figure 7j, n, and r show the error for the two variational schemes and the finitevolume scheme. Away from the pentagonal cells, the variational methods display smaller errors than the finitevolume scheme but show a significant enhancement of error at the pentagonal cells, something not found with the finitevolume scheme. As for the strain rate operator, the two variational stress divergence schemes show very similar error patterns on the unit sphere. These results are more easily seen in Fig. 7k, o, and s, which again show the detail around one of the pentagonal cells. The error enhancement around the pentagonal cell is almost entirely due to the choice of area denominator in Eq. (11). For the area denominator Fig. 7 uses the triangle area of the dual cell. Figure 10 again shows the error for the stress divergence operator on the unit sphere for the two variational schemes, but using the alternate formulation for the area denominator defined in Eq. (37). As can be seen with this alternate formulation there is no enhancement in error around the pentagonal cells. Figure 8b shows a histogram of the errors generated by the stress divergence operators. The enhanced error around the pentagonal cells for the variational scheme with the original area denominator formulation is present here as the high error tails in the distribution. The improvement for the alternate area denominator formulation is also apparent, with no such tails in the distribution for this formulation. Scaling for the stress divergence L_{2} error norm against grid resolution is shown in Fig. 9b, where all methods display firstorder accuracy, and the improvement in error from the alternate area denominator formulation is evident. The L_{∞} error norm scaling for the stress divergence operators is shown in Fig. 9d, where the effect of the enhanced error for the variational schemes around pentagonal cells is evident as the failure in L_{∞} error convergence. The alternate area denominator methods show much better convergence at low resolutions, but convergence slows for higher resolutions. The finitevolume method shows linear convergence for all grid resolutions.
Overall, the results for the unit sphere test case show similar error characteristics to the planar test case of Sect. 3.1.1. The largest effect of moving from a regular planar mesh to a unit sphere is that the variational scheme shows significantly enhanced errors for the stress divergence operator for the irregular cells surrounding the 12 pentagonal cells present on the unit sphere. This effect is ameliorated by using the alternate area denominator scheme (Eq. 37).
3.1.3 Velocity solver in a square domain
Since the MPAS framework and MPASSeaice support quadrilateral grids, direct comparisons can be made between MPASSeaice and CICE. For idealized planar test cases it is possible to set up MPASSeaice to have a virtually identical velocity solver algorithm to CICE. This is achieved by using the variational scheme with Wachspress basis functions and defining the u and v velocity component directions in the same sense as CICE. To compare MPASSeaice to CICE, we use a simple test case, similar to that used in Hunke (2001). This test case has a square planar domain of size 80 km. Ice thickness is fixed at 2 m, while ice concentration increases linearly in the eastwards direction from zero at the western boundary to one at the eastern boundary, and no snow is present. Only the velocity solver is active, with no advection or column physics. The sea ice is forced by atmospheric winds and ocean currents. The atmospheric wind forcing has the form
while the ocean currents have the form
where x and y are the horizontal position, and L_{x} and L_{y} are the domain size in the u and v directions, respectively. These forcing velocity fields are shown in Fig. 11. Seaice velocity was simulated for four time steps (each of length 1 h), which was sufficient time for the ice state to relax to the elliptical yield curve. Figure 12 shows a comparison of the modeled eastwards velocity and stress divergence component between CICE and MPASSeaice. In this comparison MPASSeaice uses an identical quadrilateral mesh to CICE. The eastwards component of wind stress pushes the sea ice against the east model boundary, and it is here that significant internal seaice stresses are present (see Fig. 12e). The figure also shows that the three MPASSeaice schemes for calculating stress divergence are capable of reproducing the results of CICE. As expected, the variational scheme with the Wachspress basis function best reproduces the results of CICE since this algorithm is most similar to CICE. Differences with CICE appear as noise, a function of incompletely damped elastic waves from the EVP rheology (Hunke, 2001). Figure 13 shows similar results for the same test case but with MPASSeaice using a regular hexagonal mesh. Here differences between the finitevolume and variational scheme with the PWL basis and the variational scheme with the Wachspress basis are larger than for quadrilateral meshes, but still small. The majority of the differences between the methods, such as the blue linear feature in Fig. 13b and c, are caused by differences in calculated strain rate. Compared to the quadrilateral case, for hexagonal cells there are larger differences between the derivatives of the Wachspress and PWL basis functions at cell vertices (used in Eq. 33). As can be seen from Fig. 14, all the schemes have stress states that lie within or on the elliptical yield curve for both quadrilateral and hexagonal meshes. A banding structure to the principal stresses can be seen for the quadrilateral meshes. Each band corresponds to grid cells in a vertical column in the top righthand corner of the domain.
3.2 Transport
To verify that the incremental remapping transport scheme works as expected, we ran two test cases on a global spherical grid, following Lipscomb and Ringler (2005). In each case there is a steady eastward velocity field given by $\mathit{u}=({u}_{\mathrm{0}}\mathrm{cos}\mathit{\theta},\mathrm{0})$, where ${u}_{\mathrm{0}}=\left(\mathrm{2}\mathit{\pi}R\right)/\left(\mathrm{12}\phantom{\rule{0.125em}{0ex}}\mathrm{d}\right)$, and R is the Earth's radius. We first advect a circular region of ice that has initial concentration given by a cosine bell within a distance $R/\mathrm{3}$ of a central point on the Equator, and a=0 elsewhere. The initial ice thickness is h=1. We compare results from the IR scheme to a simple firstorder upwind scheme to demonstrate the improvements in numerical diffusion gained by increasing the order of transport scheme. For both the IR scheme and the firstorder upwind scheme, the model was run at several grid resolutions for 12 d, at which time a perfect advection scheme would give a solution equal to the initial condition. For a grid resolution of 120 km, Fig. 15a shows equatorial cross sections of a for the initial condition, the upwind solution, and the IR solution. Figure 16a–c show the spatial distribution of ice concentration before and after the experiment for the same resolution. As expected, the upwind solution is very diffuse, while the IR scheme does a much better job of preserving the initial shape.
Next, we advect a slotted cylinder with initial concentration a=1, initial thickness h=1, and radius $R/\mathrm{2}$, also centered on the Equator. We set $a=\mathrm{0},h=\mathrm{0}$ for $r>R/\mathrm{2}$ and also in a slot of width $R/\mathrm{6}$ and length $\mathrm{5}R/\mathrm{6}$, with the long axis perpendicular to the flow. The model was run for 12 d at several resolutions. Figure 15b shows the initial condition and the upwind and IR solutions along the Equator at a resolution of 120 km, while Fig. 16d–f show the spatial distribution of the ice concentration before and after the experiment for the same resolution. Again, the upwind scheme is very diffusive; all traces of the slot vanish. The IR scheme does well at maintaining the initial plateaus and a distinct slot, although diffusion into the slot raises the minimum concentration from 0 to ∼ 0.2. Early in the IR simulation there is truncation at the leading and trailing edges of the cylinder, where the gradient is limited, but advection continues thereafter with little change in shape. The maximum concentration is just above 1 because the discretized velocity field is slightly convergent, and diffusion is small. On a plane (not shown) with steady $\mathit{u}=({u}_{\mathrm{0}},\mathrm{0})$, the discretized velocity field is nonconvergent, and a remains bounded by [0,1]. Ice thickness, having been initialized to h=1 everywhere, remains h=1 everywhere (within roundoff), showing that both schemes preserve tracer monotonicity as expected.
Figure 17 shows the L_{2} error norm of the 12 d solution for four grid resolutions ranging from 60 to 480 km (where resolution is taken as the mean distance between neighboring cell centers). Figure 17 shows that the IR solution converges with closetosecondorder accuracy (indicated by the dotted diagonal line) for the cosine bell and converges slightly below firstorder accuracy for the slotted cylinder. This slow convergence for the slotted cylinder is the result of the ice concentration discontinuity at the cylinder edge, which becomes sharper as the distance between neighboring grid cells decreases when resolution increases. The upwind scheme converges more slowly than the IR scheme, with larger errors at all resolutions. Some alongmotion asymmetry is visible in the IR solutions. This is also visible in IR solutions in Lipscomb and Ringler (2005) and is expected since IR is an upstreambased method.
3.3 Column physics
To validate the column physics in MPASSeaice we make use of the fact that CICE and MPASSeaice share identical code. CICE and MPASSeaice were run with identical forcing and with dynamics disabled. Results from the two models were bitforbit identical, indicating a correct implementation of the column physics in MPASSeaice.
3.4 Global simulations
To validate the full MPASSeaice model in a global setting we perform standalone global simulations and compare simulation results to observational datasets and to the results of simulations conducted with the CICE model (Hunke et al., 2015; release version 5.1.2). We choose version 5.1.2 of CICE to compare against to keep the comparison as clean as possible since this is the CICE version where the CICE and MPASSeaice column physics codes diverged. To aid the comparison to CICE we run both models with a 1^{∘} displaced pole quadrilateral mesh. We perform the simulation for 50 years from 1958 to 2007. Settings for the column physics are the standard ones for CICE (Hunke et al., 2015). For atmospheric and oceanic forcing we repeat the methods used by Hunke et al. (2013) and Hunke and Holland (2007). Air temperature, air specific humidity, and air velocity at 10 m height and 6hourly frequency are taken from the Coordinated Oceanice Reference Experiments (CORE) Corrected InterAnnual Forcing Version 2.0 (Large and Yeager, 2009; Griffies et al., 2009). Monthly climatologies of precipitation (Griffies et al., 2009) and cloudiness (Röske, 2001) are also used. Downwelling shortwave radiation is calculated from the monthly climatology of cloudiness using the Arctic Ocean Model Intercomparison Project (AOMIP) shortwave forcing formula (Hunke et al., 2015). Downwelling longwave radiation is calculated according to Rosati and Miyakoda (1988). Oceanic inputs, consisting of sea surface salinity, initial sea surface temperature, currents, seasurface slope, and deep ocean heat flux, come from monthly mean output of 20 years of a Community Climate System Model (CCSM) climate run (b30.009; Collins et al., 2006). The sea surface temperature is determined by a thermodynamic ocean mixed layer parameterization as used in Hunke et al. (2013). All input forcing fields are interpolated linearly in time, although the MPAS forcing functionality can be easily extended to allow interpolation in time with arbitrary order. To get good agreement between CICE and MPASSeaice it was necessary to fix several implementation errors in the CICE 5.1.2 forcing scheme. First, CICE 5.1.2 incorrectly repeats the rotation from geographical to coordinate directions for ocean current climatology data. Second, the ocean current forcing routine, rather than reading the surface layer ocean current data for all 12 months of the climatology, instead reads in the first 12 vertical layers for January. These issues have been fixed in CICE 6+.
Figure 18a–f compare total seaice extent between the MPASSeaice and CICE models and observational values for the years 1988 to 2007 inclusive. The observational seaice extent values for the Northern Hemisphere (Cavalieri and Parkinson, 2012; Parkinson et al., 1999) and Southern Hemisphere (Parkinson and Cavalieri, 2012; Zwally et al., 2002) show excellent agreement with both models, with the seasonal cycle of seaice extent well represented in both hemispheres. The largest discrepancy occurs in the Southern Hemisphere, where austral summertime seaice extent is too low in both models.
Figure 19 shows a similar agreement, comparing seaice concentration from Special Sensor Microwave/Imager (SSMI) observations using the NASATeam method (Cavalieri et al., 1996) to model results for summer and winter periods in both hemispheres. Minor differences are present in both models at the Arctic ice edge during winter and in the pack interior in summer. In general the seaice extent is well reproduced. In the Southern Hemisphere the seaice extent is reasonably reproduced in summer by both models, with more significant differences in the pack interior. As expected from Fig. 18, larger differences are found between the model results and observations in the Southern Hemisphere summer, where ice concentration is particularly underrepresented in the models in the Weddell Sea. In general, agreement is much closer between the two models than between the models and observations. This is expected given the similarity of the models and model forcing. Differences between MPASSeaice and CICE are explained by a number of differences in implementation between these models. Firstly, since MPASSeaice removes land cells, interpolation between tracer (T) and velocity (U) points (cell centers and vertices in MPAS parlance) does not include zero values for land cells, unlike in CICE. Weights in this interpolation also do not sum exactly to one for CICE since the CICE interpolation scheme mixes T and U cell areas. Secondly, CICE determines the grid angle with respect to geographical coordinates for its T points from averaging over the angle values for surrounding U points. This generates errors in wind and ocean current forcing directions around the North Pole. Finally, for ocean forcing, MPASSeaice and CICE have different orders of operation for interpolation in time, space, and rotation from geographical to model coordinate directions, generating small differences in forcing values.
Total seaice volume for the Northern Hemisphere is compared between the models and the PanArctic Ice Ocean Modeling and Assimilation System (PIOMAS) assimilated data product (Schweiger et al., 2011) in Fig. 18g–i. Both models and PIOMAS have the expected seasonal cycle with a similar variation between summer and winter and a decrease in total volume over time. Only small differences exist between the two models and the PIOMAS product in terms of absolute ice volume. Figure 20 shows the spatial patterns of seaice thickness for the Northern Hemisphere in summer, autumn, and winter compared to observations of seaice thickness from ICESat (see Fig. 20 in this paper; Yi and Zwally, 2009). ICESat data are available from set periods from 2003 to 2008 during these seasons, and the model climatological maps are generated for the same periods. ICESat observations exclude sea ice with concentration less than 20 %, so seaice thicknesses were excluded from the model results in the comparison where model ice concentration was less than 20 %. Similar spatial patterns of sea ice are found in the results of both models and the ICESat observations, with thicker sea ice along the Canadian archipelago coast and thinner sea ice everywhere at the end of the summer melt season. Both MPASSeaice and CICE have excess seaice thickness compared to ICESat observations in the Beaufort Sea and western Arctic basin and a deficit of seaice thickness in the Eurasian basin.
Finally, we perform simulations with MPASSeaice on a quasiuniform SCVT mesh with 30 km cell separation. The mesh is prepared using the Jigsaw tool (Engwirda, 2017), and the equatorial region is removed for computational efficiency. Figure 21 shows differences in northern and southernhemispheric total seaice extent and volume between this 30 km SCVT and the 1^{∘} quadrilateral mesh used above. Results are shown as a percent difference between the meshes with both simulations using the Wachspress variational scheme. Agreement is generally very good, with a difference of only a few percent in extent and volume between the meshes, with the quadrilateral mesh having smaller extent than the SCVT mesh in the Northern Hemisphere and larger extent in the Southern Hemisphere, while volume is lower in both hemispheres for the quadrilateral mesh. The differences also have a strong seasonal variation. We compare several of the other operator methods with the Wachspress variational scheme, all using the quasiuniform SCVT mesh. While using the alternate area denominator significantly reduced the errors surrounding pentagonal cells in the unit sphere test case in Sect. 3.1.2, it has almost no effect on total ice extent or volume in basinscale simulations. Compared to the Wachspress variational scheme, the PWL variational and finitevolume schemes have less effect on ice extent and a similar effect on ice volume as compared to the difference between the quasiuniform SCVT and 1^{∘} quadrilateral meshes. Differences are also strongly seasonal for these changes in operator, especially for total volume.
Through the MPAS framework, MPASSeaice incorporates code parallelization through domain decomposition and message passing with the Message Passing Interface (MPI) library. To assess the computational performance of MPASSeaice when run on multiple processors, we perform a strongscaling performance analysis (Fig. 22a). Fixing the grid resolution and run duration, we vary the number of processors and compare the time taken to perform the simulation. We measure model performance in simulated years per day of wall clock time (SYPD). The SYPD metric excludes time spent initializing or finalizing the model simulation but includes time spent reading timevarying forcing data. No output data files are written during these simulations. Simulations were performed as in Sect. 3.4 with a duration of 10 d. Two MPAS meshes are compared. The first is a global quasiuniform SCVT mesh with celltocell distance of 60 km (QU60km). This mesh has 114 539 cell centers and 234 609 vertices. The MPAS framework allows arbitrary regions of the domain to be removed. We use this capability with MPASSeaice by removing equatorial cells where sea ice does not form. This significantly decreases the size of the computational domain and increases computational performance. We use this capability in the second mesh compared, which is the same as the QU60km mesh but with equatorial cells removed, resulting in a mesh with 33 070 cell centers and 69 482 vertices, ∼ 29 % and ∼ 30 % of original cell centers and vertices in the global QU60km mesh. To determine which cells to keep and which to remove, simulations with the full mesh are used to determine a mask of cells where sea ice had ever existed during those simulations. These cells, as well as a buffer region of surrounding cells extending 1000 km further, are kept in the reduced mesh. Removal of equatorial cells produces extra domain boundaries midocean. These are treated as regular land–ocean boundaries, but because of the buffer region mentioned above, during physically reasonable simulations sea ice will not encounter them. Simulations are performed on the Anvil machine at the Laboratory Computing Resource Center at Argonne National Laboratory, which consists of 240 nodes with 36 cores per node.
Load balancing is an issue with seaice modeling since the presence of sea ice at high latitudes requires more computation for cells located in these regions than in equatorial regions. If the computational domain is partitioned without taking this into account, processors computing highlatitude cells will take longer to compute a time step than processors containing equatorial cells, which will have to wait periodically for highlatitude processors to catch up. This waiting time is wasted and contributes to poor performance. We have implemented three methods to deal with this issue. The first is removing equatorial cells as discussed above. In addition, two improved methods for partitioning the domain across processors have been developed. The standard domain partitioning method for MPAS uses the metis tool (Karypis and Kumar, 1999) to evenly divide the domain amongst processors without taking into account load balancing (labeled “simple” in Fig. 22). The two improved domain partitioning methods aim to give fewer cells and work to processors computing highlatitude regions. The first improved partitioning method (labeled “weight” in Fig. 22) adds a weight to cells in the partition; metis then aims, during the partitioning, to set the sum of weights in each partition equal. Giving highlatitude cells a larger weight than equatorial cells means fewer cells are included in highlatitude partitions, improving load balancing; metis requires this weight to be an integer, which we set to the nearest integer to (1+fn), where f is the fraction of the time the cell in question contains sea ice from a previous simulation. We find that a value of 4 for n maximizes performance. The second improved partitioning method (labeled “region” in Fig. 22) divides the globe into a polar and equatorial region based on the ice presence mask derived from previous simulations. These two regions are individually partitioned, and each processor's domain consists of one partition from the equatorial region and one from the polar region. A processor's computational domain then consists of two discontiguous regions, one polar and one equatorial.
Figure 22a shows slightly sublinear strongscaling performance for MPASSeaice for both meshes and for the “region” and “weight” partitioning methods below around 400 cores. Above about 400 cores the computational cost of exchanging halo information between processors begins to dominate, and linear scaling is no longer expected. The “simple” partition method also mostly shows nearlinear scaling, except at around 10 processors, when load balancing issues begin to affect performance relative to the other partition methods. As expected, removing the equatorial cells reduces this effect. Both the “region” and “weight” partitioning methods improve load balancing by around the same amount. Choice of equatorial mesh removal and partition method can affect performance by up to a factor of around 4. The “region” partition method seems to underperform at high processor number once we reach the limit of strong scaling. This was found to be caused by a deficiency in the generated partitions where, at high processor number, boundaries between partitions would become tortuous. Figure 22b shows the same performance results plotted against the average number of cells per core. Sublinear scaling is found when the number of cells per core is greater than ∼ 300, with a degradation in performance for simulations with fewer cells per core than this value. A similar result was found for FiniteElement/volumE Sea iceOcean Model (FESOM) (Koldunov et al., 2019).
Comparison of computational performance between MPASSeaice and CICE is nontrivial. While MPASSeaice can be used with quadrilateral meshes, its primary use is expected to be with SCVT meshes consisting primarily of hexagonal cells. The ratio of velocity points to cell centers for these meshes is approximately 2, whereas the quadrilateral grids used by CICE have approximately the same number of velocity points as cell centers. Comparing performance of MPASSeaice using SCVT meshes to CICE using quadrilateral meshes then depends on whether the performance is compared based on the number of cell centers or velocity points. To give an approximate idea of the relative performance of MPASSeaice and CICE we compare the SYPD achieved on Anvil for 10 d of the simulation described in Sect. 3.4 with both MPASSeaice and CICE using a fully global 1^{∘} quadrilateral grid. Simulation throughput for MPASSeaice and CICE is listed in Table 1 as simulated years per day (SYPD) for the whole model and for the velocity solver, transport, and column schemes. For the total model MPASSeaice achieved approximately 70 % of the CICE throughput. The better computational performance of CICE is expected since the unstructured mesh in MPASSeaice necessitates less efficient memory access patterns. As a percentage of CICE model performance, the MPASSeaice velocity solver displayed the poorest throughput and the column physics the most competitive throughput.
We have described a new seaice model, MPASSeaice, and successfully validated the velocity solver and transport schemes in idealized test cases, on both planar and spherical grids. These schemes are closely based on those implemented on the quadrilateral grid used in the CICE seaice model but adapted for the polygonal cells of MPAS meshes. When using the variational scheme with Wachspress basis functions and a quadrilateral MPAS mesh, the velocity solver of MPASSeaice replicates the velocity solver algorithm of CICE, allowing rapid testing and validation.
We developed several other schemes for the strain rate and stress divergence spatial operators to compare with the variational Wachspress scheme. We find that, while the variational scheme, with the alternate area denominator formulation, has excellent error characteristics for the stress divergence operator, the onesided stencil of the variational strain rate operators results in poor error characteristics. The finitevolume scheme shows the opposite, with good error characteristics for the strain rate operators, but asymmetric integrals around the dual triangles of the SCVT mesh result in larger errors for the finitevolume stress divergence operator. This suggests that the variational scheme could be improved by modifying its strain rate operator to have a twosided stencil. We investigated several averaging techniques to implement a twosided stencil for the variational scheme, which resulted in improved error characteristics for the variational strain rate operator. For basinscale seaice simulations, however, these alternate operators had only a small effect on simulation results.
MPASSeaice and CICE share the sophisticated suite of column physics and BGC originally developed in CICE, again allowing the rapid development of MPASSeaice. Global simulations with realistic forcing have validated MPASSeaice against similar simulations with CICE and against observations for seaice concentration, extent, and volume. MPASSeaice has been coupled into the Energy Exascale Earth System Model (Golaz et al., 2019; Burrows et al., 2020), and the validation experiments described here give confidence in the seaice results from E3SM simulations. MPASSeaice shows powerlaw strongscaling performance with a nearly linear exponent, and the flexibility in mesh partitioning afforded by its use of an unstructured mesh allows efficient load balancing.
Future work will assess the fidelity and performance of MPASSeaice on variableresolution meshes and examine more recent metrics for evaluating seaice dynamics, such as new statistical metrics of linear kinematic features (e.g., Hutter and Losch, 2020). Potential challenges with variableresolution meshes include assessing the resolution invariance of seaice rheologies and developing resolutioninvariant versions, efficiently using future heterogeneous computing architectures, as well as generating efficient domain partitions of highly variableresolution meshes.
MPASSeaice v1.0.0 is released as part of the Energy Exascale Earth System Model (E3SM) version 2 and is available at https://doi.org/10.11578/E3SM/dc.20210927.1 (E3SM Project, DOE, 2021). Model forcing data are generated by scripts included in the E3SM repository from data available at https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html (NOAA, 2021). Other datasets used in this paper are available at https://doi.org/10.5281/zenodo.6230907 (Turner, 2022).
AKT developed and implemented the velocity solver and performed all simulations and analysis. WHL developed and implemented the transport scheme. ECH, AKT, and NJ packaged the thermodynamics into the Icepack library used by MPASSeaice. DE provided support for the analysis of the horizontal operators. TDR contributed to the development of the velocity solver. JDW provided software support. DWJ provided support for use of the MPAS framework. AKT prepared the manuscript with contributions from all coauthors.
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This research was supported as part of the Energy Exascale Earth System Model (E3SM) project, funded by the US Department of Energy, Office of Science, Office of Biological and Environmental Research. The National Center for Atmospheric Research is a major facility sponsored by the National Science Foundation under cooperative agreement no. 1852977. This research used a highperformance computing cluster provided by the BER Earth System Modeling program and operated by the Laboratory Computing Resource Center at Argonne National Laboratory. The authors would like to thank Michael Duda of the National Center for Atmospheric Research for his contributions to the MPAS software infrastructure and Sara Calandrini (LANL) for help with mesh generation. Adrian K. Turner would like to thank the Linda Hall Library for access to their collection.
This research has been supported by the Office of Science (grant no. 89233218CNA000001).
This paper was edited by Riccardo Farneti and reviewed by Sergey Danilov and two anonymous referees.
Arakawa, A. and Lamb, V. R.: Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, Methods in Computational Physics: Advances in Research and Applications, 17, 173–265, https://doi.org/10.1016/B9780124608177.500094, 1977. a, b
Bailey, T. S., Adams, M. L., Yang, B., and Zika, M. R.: A piecewise linear finite element discretization of the diffusion equation for arbitrary polyhedral grids, J. Comput. Phys., 227, 3738–3757, https://doi.org/10.1016/j.jcp.2007.11.026, 2008. a
Batchelor, G. K.: An introduction to fluid dynamics, 1st edn., Cambridge University Press, ISBN 9780521663960, https://doi.org/10.1017/CBO9780511800955, 1967. a, b, c
Bern, M. and Plassmann, P.: Mesh Generation, in: Handbook of Computational Geometry, 1st edn., edited by: Sack, J. R. and Urrutia, J., chap. 10, pp. 291–332, NorthHolland, Amsterdam, ISBN 9780444825377, https://doi.org/10.1016/B9780444825377.X50001, 2000. a
Bitz, C. M. and Lipscomb, W. H.: An energyconserving thermodynamic model of sea ice, J. Geophys. Res.Oceans, 104, 15669–15677, https://doi.org/10.1029/1999JC900100, 1999. a
Briegleb, B. P. and Light, B.: A DeltaEddington multiple scattering parameterization for solar radiation in the sea ice component of the Community Climate System Model, Tech. Rep. NCAR/TN472+STR, National Center for Atmospheric Research, Boulder, Colorado USA, https://doi.org/10.5065/D6B27S71, 2007. a
Burrows, S. M., Maltrud, M., Yang, X., Zhu, Q., Jeffery, N., Shi, X., Ricciuto, D., Wang, S., Bisht, G., Tang, J., Wolfe, J., Harrop, B. E., Singh, B., Brent, L., Baldwin, S., Zhou, T., CameronSmith, P., Keen, N., Collier, N., Xu, M., Hunke, E. C., Elliott, S. M., Turner, A. K., Li, H., Wang, H., Golaz, J.C., BondLamberty, B., Hoffman, F. M., Riley, W. J., Thornton, P. E., Calvin, K., and Leung, L. R.: The DOE E3SM v1.1 Biogeochemistry Configuration: Description and Simulated EcosystemClimate Responses to Historical Changes in Forcing, J. Adv. Model. Earth Sy., 12, e2019MS001766, https://doi.org/10.1029/2019MS001766, 2020. a
Cavalieri, D. J. and Parkinson, C. L.: Arctic sea ice variability and trends, 1979–2010, The Cryosphere, 6, 881–889, https://doi.org/10.5194/tc68812012, 2012. a
Cavalieri, D. J., Parkinson, C. L., Gloersen, P., and Zwally, H. J.: Sea ice concentrations from Nimbus7 SMMR and DMSP SSM/ISSMIS passive microwave data, version 1, Tech. Rep., NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado USA, https://doi.org/10.5067/8GQ8LZQVL0VL, 1996. a
Chen, C., Gao, G., Qi, J., Proshutinsky, A., Beardsley, R. C., Kowalik, Z., Lin, H., and Cowles, G.: A new highresolution unstructured grid finite volume Arctic Ocean model (AOFVCOM): An application for tidal studies, J. Geophys. Res.Oceans, 114, C08017, https://doi.org/10.1029/2008JC004941, 2009. a
Collins, W. D., Bitz, C. M., Blackmon, M. L., Bonan, G. B., Bretherton, C. S., Carton, J. A., Chang, P., Doney, S. C., Hack, J. J., Henderson, T. B., Kiehl, J. T., Large, W. G., McKenna, D. S., Santer, B. D., and Smith, R. D.: The Community Climate System Model Version 3 (CCSM3), J. Climate, 19, 2122–2143, https://doi.org/10.1175/JCLI3761.1, 2006. a
Danabasoglu, G., Lamarque, J.F., Bacmeister, J., Bailey, D. A., DuVivier, A. K., Edwards, J., Emmons, L. K., Fasullo, J., Garcia, R., Gettelman, A., Hannay, C., Holland, M. M., Large, W. G., Lauritzen, P. H., Lawrence, D. M., Lenaerts, J. T. M., Lindsay, K., Lipscomb, W. H., Mills, M. J., Neale, R., Oleson, K. W., OttoBliesner, B., Phillips, A. S., Sacks, W., Tilmes, S., van Kampenhout, L., Vertenstein, M., Bertini, A., Dennis, J., Deser, C., Fischer, C., FoxKemper, B., Kay, J. E., Kinnison, D., Kushner, P. J., Larson, V. E., Long, M. C., Mickelson, S., Moore, J. K., Nienhouse, E., Polvani, L., Rasch, P. J., and Strand, W. G.: The Community Earth System Model Version 2 (CESM2), J. Adv. Model. Earth Sy., 12, e2019MS001916, https://doi.org/10.1029/2019MS001916, 2020. a
Danilov, S., Wang, Q., Timmermann, R., Iakovlev, N., Sidorenko, D., Kimmritz, M., Jung, T., and Schröter, J.: FiniteElement Sea Ice Model (FESIM), version 2, Geosci. Model Dev., 8, 1747–1761, https://doi.org/10.5194/gmd817472015, 2015. a
Dasgupta, G.: Interpolants within convex polygons: Wachpress' shape functions, J. Aerospac. Eng., 16, 1–8, https://doi.org/10.1061/(ASCE)08931321(2003)16:1(1), 2003. a
Dukowicz, J. K. and Baumgardner, J. R.: Incremental remapping as a transport/advection algorithm, J. Comput. Phys., 160, 318–335, https://doi.org/10.1006/jcph.2000.6465, 2000. a, b
Dunavant, D. A.: High degree efficient symmetrical Gaussian quadrature rules for the triangle, Int. J. Numer. Meth. Eng., 21, 1129–1148, https://doi.org/10.1002/nme.1620210612, 1985. a, b, c
E3SM Project, DOE: Energy Exascale Earth System Model v2.0, DOE Code [code], https://doi.org/10.11578/E3SM/dc.20210927.1, 2021. 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, 2017. a
Flocco, D., Feltham, D. L., and Turner, A. K.: Incorporation of a physically based melt pond scheme into the sea ice component of a climate model, J. Geophys. Res.Oceans, 115, C08012, https://doi.org/10.1029/2009JC005568, 2010. a
Gao, G., Chen, C., Qi, J., and Beardsley, R. C.: An unstructuredgrid, finitevolume sea ice model: Development, validation, and application, J. Geophys. Res.Oceans, 116, C00D04, https://doi.org/10.1029/2010JC006688, 2011. a
Golaz, J.C., Caldwell, P. M., Van Roekel, L. P., Petersen, M. R., Tang, Q., Wolfe, J. D., Abeshu, G., Anantharaj, V., AsayDavis, X. S., Bader, D. C., Baldwin, S. A., Bisht, G., Bogenschutz, P. A., Branstetter, M., Brunke, M. A., Brus, S. R., Burrows, S. M., CameronSmith, P. J., Donahue, A. S., Deakin, M., Easter, R. C., Evans, K. J., Feng, Y., Flanner, M., Foucar, J. G., Fyke, J. G., Griffin, B. M., Hannay, C., Harrop, B. E., Hoffman, M. J., Hunke, E. C., Jacob, R. L., Jacobsen, D. W., Jeffery, N., Jones, P. W., Keen, N. D., Klein, S. A., Larson, V. E., Leung, L. R., Li, H.Y., Lin, W., Lipscomb, W. H., Ma, P.L., Mahajan, S., Maltrud, M. E., Mametjanov, A., McClean, J. L., McCoy, R. B., Neale, R. B., Price, S. F., Qian, Y., Rasch, P. J., Reeves Eyre, J. E. J., Riley, W. J., Ringler, T. D., Roberts, A. F., Roesler, E. L., Salinger, A. G., Shaheen, Z., Shi, X., Singh, B., Tang, J., Taylor, M. A., Thornton, P. E., Turner, A. K., Veneziani, M., Wan, H., Wang, H., Wang, S., Williams, D. N., Wolfram, P. J., Worley, P. H., Xie, S., Yang, Y., Yoon, J.H., Zelinka, M. D., Zender, C. S., Zeng, X., Zhang, C., Zhang, K., Zhang, Y., Zheng, X., Zhou, T., and Zhu, Q.: The DOE E3SM Coupled Model Version 1: Overview and Evaluation at Standard Resolution, J. Adv. Model. Earth Sy., 11, 2089–2129, https://doi.org/10.1029/2018MS001603, 2019. a, b
Griffies, S. M., Biastoch, A., Böning, C., Bryan, F., Danabasoglu, G., Chassignet, E. P., England, M. H., Gerdes, R., Haak, H., Hallberg, R. W., Hazeleger, W., Jungclaus, J., Large, W. G., Madec, G., Pirani, A., Samuels, B. L., Scheinert, M., Gupta, A. S., Severijns, C. A., Simmons, H. L., Treguier, A. M., Winton, M., Yeager, S., and Yin, J.: Coordinated Oceanice Reference Experiments (COREs), Ocean Model., 26, 1–46, https://doi.org/10.1016/j.ocemod.2008.08.007, 2009. a, b
Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, https://doi.org/10.1175/15200485(1979)009<0815:ADTSIM>2.0.CO;2, 1979. a, b
Hoffman, M. J., Perego, M., Price, S. F., Lipscomb, W. H., Zhang, T., Jacobsen, D., Tezaur, I., Salinger, A. G., Tuminaro, R., and Bertagna, L.: MPASAlbany Land Ice (MALI): a variableresolution ice sheet model for Earth system modeling using Voronoi grids, Geosci. Model Dev., 11, 3747–3780, https://doi.org/10.5194/gmd1137472018, 2018. a
Holland, M. M., Bailey, D. A., Briegleb, B. P., Light, B., and Hunke, E.: Improved sea ice shortwave radiation physics in CCSM4: The impact of melt ponds and aerosols on arctic sea ice, J. Climate, 25, 1413–1430, https://doi.org/10.1175/JCLID1100078.1, 2012. a, b
Hunke, E., Allard, R., Bailey, D., Blain, P., Craig, T., Damsgaard, A., Dupont, F., DuVivier, A., Grumbine, R., Holland, M., Jeffery, N., Lemieux, J.F., Roberts, A., Turner, M., and Winton, M.: CICE Consortium/Icepack version 1.1.0 (Version Icepack1.1.0), Zenodo [data set], https://doi.org/10.5281/zenodo.1891650, 2018. a
Hunke, E. C.: Viscousplastic sea ice dynamics with the EVP model: Linearization issues, J. Comput. Phys., 170, 18–38, https://doi.org/10.1006/jcph.2001.6710, 2001. a, b, c
Hunke, E. C. and Dukowicz, J. K.: An elasticviscousplastic model for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867, https://doi.org/10.1175/15200485(1997)027<1849:AEVPMF>2.0.CO;2, 1997. a, b, c, d
Hunke, E. C. and Dukowicz, J. K.: The elasticviscousplastic sea ice dynamics model in general orthogonal curvilinear coordinates on a sphere – Incorporation of metric terms, Mon. Weather Rev., 130, 1848–1865, https://doi.org/10.1175/15200493(2002)130<1848:TEVPSI>2.0.CO;2, 2002. a, b, c, d, e, f, g
Hunke, E. C. and Holland, M. M.: Global atmospheric forcing data for Arctic iceocean modeling, J. Geophys. Res.Oceans, 112, C04S14, https://doi.org/10.1029/2006JC003640, 2007. a
Hunke, E. C., Hebert, D. A., and Lecomte, O.: Levelice melt ponds in the Los Alamos sea ice model, CICE, Ocean Model., 71, 26–42, https://doi.org/10.1016/j.ocemod.2012.11.008, 2013. a, b, c
Hunke, E. C., Lipscomb, W. H., Turner, A. K., Jeffery, N., and Elliott, S.: CICE: the Los Alamos sea ice model documentation and software user's manual version 5.1, Tech. Rep., Los Alamos National Laboratory, LACC06012, https://github.com/CICEConsortium/CICEsvntrunk/blob/main/cicedoc/cicedoc.pdf (last access: 4 April 2022), 2015. a, b, c, d, e, f
Hutchings, J. K., Jasak, H., and Laxon, S. W.: A strength implicit correction scheme for the viscousplastic sea ice model, Ocean Model., 7, 111–133, https://doi.org/10.1016/S14635003(03)000404, 2004. a
Hutter, N. and Losch, M.: Featurebased comparison of sea ice deformation in leadpermitting sea ice simulations, The Cryosphere, 14, 93–113, https://doi.org/10.5194/tc14932020, 2020. a
Ingram, W. J., Wilson, C. A., and Mitchell, J. F. B.: Modeling climate change: An assessment of sea ice and surface albedo feedbacks, J. Geophys. Res.Atmos., 94, 8609–8622, https://doi.org/10.1029/JD094iD06p08609, 1989. a
Jeffery, N. and Hunke, E. C.: Modeling the winterspring transition of firstyear ice in the western Weddell Sea, J. Geophys. Res.Oceans, 119, 5891–5920, https://doi.org/10.1002/2013JC009634, 2014. a
Jeffery, N., Elliott, S., Hunke, E. C., Lipscomb, W. H., and Turner, A.: Biogeochemistry of CICE: The Los Alamos sea ice model documentation and software user's manual, zbgc_colpkg modifications to version 5.0., Tech. Rep., LAUR1627780, Los Alamos National Laboratory, Los Alamos, NM, USA, https://doi.org/10.2172/1329842, 2016. a
Jeffery, N., Maltrud, M. E., Hunke, E. C., Wang, S., Wolfe, J., Turner, A. K., Burrows, S. M., Shi, X., Lipscomb, W. H., Maslowski, W., and Calvin, K. V.: Investigating controls on sea ice algal production using E3SMv1.1BGC, Ann. Glaciol., 61, 51–72, https://doi.org/10.1017/aog.2020.7, 2020. a
Karypis, G. and Kumar, V.: A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs, SIAM J. Sci. Comput., 20, 359–392, https://doi.org/10.1137/S1064827595287997, 1999. a
Killworth, P. D.: Deep convection in the World Ocean, Rev. Geophys., 21, 1–26, https://doi.org/10.1029/RG021i001p00001, 1983. a
Koldunov, N. V., Aizinger, V., Rakowsky, N., Scholz, P., Sidorenko, D., Danilov, S., and Jung, T.: Scalability and some optimization of the FinitevolumE Sea ice–Ocean Model, Version 2.0 (FESOM2), Geosci. Model Dev., 12, 3991–4012, https://doi.org/10.5194/gmd1239912019, 2019. a
Large, W. G. and Yeager, S. G.: The global climatology of an interannually varying air–sea flux data set, Clim. Dynam., 33, 341–364, https://doi.org/10.1007/s0038200804413, 2009. a
Lietaer, O., Fichefet, T., and Legat, V.: The effects of resolving the Canadian Arctic Archipelago in a finite element sea ice model, Ocean Model., 24, 140–152, https://doi.org/10.1016/j.ocemod.2008.06.002, 2008. a
Lipscomb, W. H.: Remapping the thickness distribution in sea ice models, J. Geophys. Res.Oceans, 106, 13989–14000, https://doi.org/10.1029/2000JC000518, 2001. a
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, b, c, d
Lipscomb, W. H. and Ringler, T. D.: An incremental remapping transport scheme on a spherical geodesic grid, Mon. Weather Rev., 133, 2335–2350, https://doi.org/10.1175/MWR2983.1, 2005. a, b, c, d, e, f
Lipscomb, W. H., Hunke, E. C., Maslowski, W., and Jakacki, J.: Ridging, strength, and stability in highresolution sea ice models, J. Geophys. Res.Oceans, 112, C03S91, https://doi.org/10.1029/2005JC003355, 2007. a
Malvern, L. E.: Introduction to the Mechanics of Continuous Medium, PrenticeHall, Inc., Englewood Cliffs, ISBN 9780134876030, 1969. a
Mehlmann, C. and Korn, P.: Seaice dynamics on triangular grids, J. Comput. Phys., 428, 110086, https://doi.org/10.1016/j.jcp.2020.110086, 2021. a
Murray, R. J.: Explicit generation of orthogonal grids for ocean models, J. Comput. Phys., 126, 251–273, https://doi.org/10.1006/jcph.1996.0136, 1996. a
NOAA: Version 2 forcing for Coordinated Oceanice Reference Experiments (CORE), NOAA [data set], https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html, last access: 28 November 2021. a
Parkinson, C. L. and Cavalieri, D. J.: Antarctic sea ice variability and trends, 1979–2010, The Cryosphere, 6, 871–880, https://doi.org/10.5194/tc68712012, 2012. a
Parkinson, C. L., Cavalieri, D. J., Gloersen, P., Zwally, H. J., and Comiso, J. C.: Arctic sea ice extents, areas, and trends, 1978–1996, J. Geophys. Res.Oceans, 104, 20837–20856, https://doi.org/10.1029/1999JC900082, 1999. a
Reddy, J. N.: An introduction to the finite element method, 2nd edn., McGrawHill, New York, ISBN: 9780070513556, 1993. a
Ringler, T., Thuburn, J., Klemp, J., and Skamarock, W.: A unified approach to energy conservation and potential vorticity dynamics for arbitrarilystructured Cgrids, J. Comput. Phys., 229, 3065–3090, https://doi.org/10.1016/j.jcp.2009.12.007, 2010. 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
Rosati, A. and Miyakoda, K.: A general circulation model for upper ocean simulation, J. Phys. Oceanogr., 18, 1601–1626, https://doi.org/10.1175/15200485(1988)018<1601:AGCMFU>2.0.CO;2, 1988. a
Röske, F.: An atlas of surface fluxes based on the ECMWF reanalysis. A climatological data set to force global ocean general circulation models, Tech. Rep., 323, MaxPlanckInst. Für Meteorol., Hamburg, Germany, http://hdl.handle.net/21.11116/000000032E854 (last access: 4 April 2022), 2001. a
Schweiger, A., Lindsay, R., Zhang, J., Steele, M., Stern, H., and Kwok, R.: Uncertainty in modeled Arctic sea ice volume, J. Geophys. Res.Oceans, 116, C00D06, https://doi.org/10.1029/2011JC007084, 2011. 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
Smith, R., Kortas, S., and Meltz, B.: Curvilinear coordinates for global ocean models, Tech. Rep. LAUR951146, Los Alamos National Laboratory, https://lanlprimo.hosted.exlibrisgroup.com/permalink/f/1hmhmmc/01LANL_ALMA5199299490003761 (last access: 4 April 2022), 1995. a
Turner, A.: MPASSeaice v1.0.0 data set, Zenodo [data set], https://doi.org/10.5281/zenodo.6230907, 2022. a
Turner, A. K. and Hunke, E. C.: Impacts of a mushylayer thermodynamic approach in global seaice simulations using the CICE seaice model, J. Geophys. Res.Oceans, 120, 1253–1275, https://doi.org/10.1002/2014JC010358, 2015. a
Turner, A. K., Hunke, E. C., and Bitz, C. M.: Two modes of seaice gravity drainage: A parameterization for largescale modeling, J. Geophys. Res.Oceans, 118, 2279–2294, https://doi.org/10.1002/jgrc.20171, 2013. a
van Leer, B.: Towards the ultimate conservative difference scheme. V. A secondorder sequel to Godunov's method, J. Comput. Phys., 32, 101–136, https://doi.org/10.1016/00219991(79)901451, 1979. a
Wang, Q., Danilov, S., Sidorenko, D., Timmermann, R., Wekerle, C., Wang, X., Jung, T., and Schröter, J.: The Finite Element Sea IceOcean Model (FESOM) v.1.4: formulation of an ocean general circulation model, Geosci. Model Dev., 7, 663–693, https://doi.org/10.5194/gmd76632014, 2014. a
Yi, D. and Zwally, H. J.: Arctic sea ice freeboard and thickness, version 1., Tech. Rep., NSIDC: National Snow and Ice Data Center, Boulder, Colorado USA [data set], https://doi.org/10.5067/SXJVJ3A2XIZT, 2009. a
Zwally, H. J., Comiso, J. C., Parkinson, C. L., Cavalieri, D. J., and Gloersen, P.: Variability of Antarctic sea ice 1979–1998, J. Geophys. Res., 107, 3041, https://doi.org/10.1029/2000JC000733, 2002. a