MPAS-Seaice (v1.0.0): sea-ice dynamics on unstructured Voronoi meshes

We present MPAS-Seaice, a sea-ice model which uses the Model for Prediction Across Scales (MPAS) framework and spherical centroidal Voronoi tessellation (SCVT) unstructured meshes. As well as SCVT meshes, MPAS-Seaice can run on the traditional quadrilateral grids used by sea-ice models such as CICE. The MPAS-Seaice 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 (“finite-volume”) 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 finite-volume 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 quasi-uniform mesh, which is ameliorated with an alternate formulation for the operator. MPAS-Seaice 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 MPAS-Seaice 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 MPAS-Seaice allow removal of equatorial model cells and flexibility in domain decomposition, improving model performance. MPAS-Seaice is the current sea-ice component of the Energy Exascale Earth System Model (E3SM).

1 Introduction

Sea ice, the frozen surface of the sea at high latitudes, is an important component of the Earth climate system. Rejection of salt during sea-ice 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 sea-ice dynamics and thermodynamics is an important tool in understanding global climate. One of the most popular sea-ice models currently in use is CICE (Hunke et al., 2015). CICE approximates the sea-ice 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 time-stepping and improved parallelization scalability of the algorithm. CICE uses a quadrilateral structured mesh and has been used with both displaced-pole (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 limited-domain 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 sea-ice models have been developed. Hutchings et al. (2004) implemented and demonstrated a finite-volume, cell-centered discretization. The Unstructured-Grid CICE (UG-CICE; Gao et al., 2011) model, which is built on the FVCOM framework (Chen et al., 2009), also uses a finite-volume formulation. In contrast, finite-element discretizations have been implemented by Lietaer et al. (2008), by Danilov et al. (2015) in the Finite-Element Sea Ice Model (FESIM), and by Mehlmann and Korn (2021) in the sea-ice component of the Icosahedral Nonhydrostatic Weather and Climate Model (ICON). UG-CICE, 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 (MPAS-O; Ringler et al., 2013), atmosphere (MPAS-A; Skamarock et al., 2012), and land ice (MALI; Hoffman et al., 2018) models. Here, we describe a new MPAS model for sea ice: MPAS-Seaice. This model uses a “B” Arakawa-type grid (Arakawa and Lamb, 1977) with sea-ice 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 MPAS-Seaice to either use the same mesh as structured quadrilateral models such as CICE or match the mesh used by unstructured ocean models such as MPAS-Ocean. Matching the ocean mesh is required for coupling sea-ice components within some global climate models such as CESM (Danabasoglu et al., 2020) and E3SM (Golaz et al., 2019), and simplifies the ocean–sea-ice coupling methodology.

We have implemented two discretizations of the momentum equation for MPAS-Seaice. The first is based on the variational scheme used by CICE (Hunke and Dukowicz, 2002). The second uses the integral (“finite-volume”) 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. MPAS-Seaice uses the same EVP rheology and column physics as CICE. Sea-ice 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 MPAS-Seaice, focusing on the solution of the sea-ice 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 quasi-uniform global meshes; variable-resolution meshes will be considered in later publications.

2 Model description

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 quasi-uniform 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–sea-ice 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 quasi-uniform 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.

Figure 1 Schematic representation of three cells in an example MPAS mesh. The mesh is composed of cell centers (circles), cell edge points (triangles), and cell vertices (squares). The dual triangular Delaunay mesh is formed by joining cell centers (dashed line). MPAS-Seaice uses a “B” type grid, with both velocity components defined at cell vertices. Unlike traditional quadrilateral meshes, the directions of the velocity components in MPAS-Seaice are not in general aligned with the cell edges.

2.2 Velocity solver

MPAS-Seaice uses a “B” Arakawa-type grid (Arakawa and Lamb, 1977) with both components of velocity defined at cell vertices and with sea-ice concentration, volume, and other tracers defined at cell centers (see Fig. 1). When using CICE-like quadrilateral meshes, the velocity solver algorithm of MPAS-Seaice reduces to that of CICE, allowing CICE and MPAS-Seaice to use identical test cases and supporting rapid testing and development of MPAS-Seaice.

In CICE the velocity components are aligned with the quadrilateral mesh. This is not possible, in general, with MPAS-Seaice 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 sea-ice velocity we solve the same sea-ice momentum equation as CICE (Hibler, 1979; Hunke and Dukowicz, 1997):

$$\begin{array}{}\text{(1)}& m{\displaystyle \frac{\partial \mathit{u}}{\partial t}}=\mathbf{\nabla }\cdot \mathit{\sigma }+{\mathit{\tau }}_{\mathrm{a}}+{\mathit{\tau }}_{\mathrm{w}}-\widehat{\mathit{k}}\times mf\mathit{u}-mg\mathbf{\nabla }{H}_{\mathrm{o}}.\end{array}$$

Here m is the mass of snow and ice per unit area; u is the sea-ice 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 second-to-last 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 non-coupled simulations we assume that the ocean currents are in geostrophic balance so that

$$\begin{array}{}\text{(2)}& mg\mathbf{\nabla }{H}_{\mathrm{o}}=mf\widehat{\mathit{k}}\times {\mathit{u}}_{\mathrm{o}},\end{array}$$

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:

1. Determine the strain rate tensor from the velocity field.

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

3. 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 finite-volume scheme uses the line-integral 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:

$$\begin{array}{}\text{(3)}& \underset{\mathrm{\Omega }}{\int }\mathit{u}\cdot (\mathbf{\nabla }\cdot \mathit{\sigma })\mathrm{d}A=-\underset{\mathrm{\Omega }}{\int }({\mathit{\sigma }}_{\mathrm{11}}{\dot{\mathit{ϵ}}}_{\mathrm{11}}+\mathrm{2}{\mathit{\sigma }}_{\mathrm{12}}{\dot{\mathit{ϵ}}}_{\mathrm{12}}+{\mathit{\sigma }}_{\mathrm{22}}{\dot{\mathit{ϵ}}}_{\mathrm{22}})\mathrm{d}A.\end{array}$$

Here $\dot{\mathit{ϵ}}$ 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 right-hand side of this equation. The stress divergence operator is derived from the functional

$$\begin{array}{}\text{(4)}& I[\mathit{u},\mathit{\sigma }]=\underset{\mathrm{\Omega }}{\int }L(\mathit{u},\mathit{\sigma })=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(5)}& L(\mathit{u},\mathit{\sigma })=\mathit{u}\cdot (\mathbf{\nabla }\cdot \mathit{\sigma })+({\mathit{\sigma }}_{\mathrm{11}}{\dot{\mathit{ϵ}}}_{\mathrm{11}}+\mathrm{2}{\mathit{\sigma }}_{\mathrm{12}}{\dot{\mathit{ϵ}}}_{\mathrm{12}}+{\mathit{\sigma }}_{\mathrm{22}}{\dot{\mathit{ϵ}}}_{\mathrm{22}}).\end{array}$$

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

$$\begin{array}{}\text{(6)}& \begin{array}{rl}I& [u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}}]\\ & =\sum _i^{n_d}\underset{i}{\int }(u{F}_u+v{F}_v)\mathrm{d}A\\ & -D(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}})=\mathrm{0},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(7)}& \begin{array}{rl}I& [u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}}]\\ & =\sum _i^{n_d}(u_i{F}_{ui}+v_i{F}_{vi})A_i\\ & -D(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}})=\mathrm{0}.\end{array}\end{array}$$

Taking the variation in this functional with respect to a single discretized velocity component, u_j, we get the following Euler–Lagrange equation

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\displaystyle \frac{\partial L(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}})}{\partial u_j}}\\ & -{\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}\left({\displaystyle \frac{\partial L(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}})}{\partial {\dot{u}}_j}}\right)=\mathrm{0},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(9)}& \begin{array}{rl}L& (u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}})\\ & =\sum _i^{n_d}(u_i{F}_{ui}+v_i{F}_{vi})A_i\\ & -D(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}}),\end{array}\end{array}$$

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

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial u_j}}& \sum _i^{n_d}(u_i{F}_{ui}+v_i{F}_{vi})A_i\\ & ={\displaystyle \frac{\partial }{\partial u_j}}D(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}}),\end{array}\end{array}$$

and consequently the u component of the stress divergence for vertex j is given by

$$\begin{array}{}\text{(11)}& F_{uj}={\displaystyle \frac{\mathrm{1}}{A_j}}{\displaystyle \frac{\partial }{\partial u_j}}D(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}}).\end{array}$$

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:

$$\begin{array}{}\text{(12)}& D=D_{\mathrm{1}}+D_{\mathrm{2}}+D_{\mathrm{3}},\end{array}$$

with

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{D}_{\mathrm{1}}& =-\int {\mathit{\sigma }}_{\mathrm{11}}{\dot{\mathit{ϵ}}}_{\mathrm{11}}\mathrm{d}A,D_{\mathrm{2}}=-\int \mathrm{2}{\mathit{\sigma }}_{\mathrm{12}}{\dot{\mathit{ϵ}}}_{\mathrm{12}}\mathrm{d}A,\\ & D_{\mathrm{3}}=-\int {\mathit{\sigma }}_{\mathrm{22}}{\dot{\mathit{ϵ}}}_{\mathrm{22}}\mathrm{d}A.\end{array}\end{array}$$

We calculate the contribution to F_u and F_v from D₁. Similar contributions come from D₂ and D₃. Using the expression for ${\dot{\mathit{ϵ}}}_{\mathrm{11}}$ in terms of the velocity components and latitude θ, D₁ becomes

$$\begin{array}{}\text{(14)}& D_{\mathrm{1}}=-\int {\mathit{\sigma }}_{\mathrm{11}}\left[{\displaystyle \frac{\partial u}{\partial x}}-{\displaystyle \frac{v\tan \mathit{\theta }}{r}}\right]\mathrm{d}A,\end{array}$$

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{ϵ}}$ 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:

$$\begin{array}{}\text{(15)}& D_{\mathrm{1}}=-\sum _k^{n_p}\underset{k}{\int }{\mathit{\sigma }}_{\mathrm{11}}\left[{\displaystyle \frac{\partial u}{\partial x}}-{\displaystyle \frac{v\tan \mathit{\theta }}{r}}\right]\mathrm{d}A,\end{array}$$

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

$$\begin{array}{}\text{(16)}& \mathit{\psi }(x,y)=\sum _l^{n_v}{\mathit{\psi }}_l{\mathcal{W}}_l(x,y),\end{array}$$

where the sum is over the n_v vertices of the cell in the primary mesh. Necessary properties for these basis functions are that

$$\begin{array}{}\text{(17)}& \sum _l^{n_v}{\mathcal{W}}_l(x,y)=\mathrm{1}\end{array}$$

across the cell and that

$$\begin{array}{}\text{(18)}& {\mathcal{W}}_l(x,y)\left\{\begin{array}{ll}\mathrm{1}& \text{if}(x,y)\text{at vertex}l\\ \mathrm{0}& \text{if}(x,y)\text{at any other vertex.}\end{array}\right.\end{array}$$

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)

$$\begin{array}{}\text{(19)}& {\mathcal{W}}_i={\displaystyle \frac{{\mathcal{N}}_i}{{\sum }_j^n{\mathcal{N}}_j}},\end{array}$$

where

$$\begin{array}{}\text{(20)}& {\mathcal{N}}_i(x,y)={\mathit{\kappa }}_i\prod _{j\ne i,j\ne i+\mathrm{1}}^{j=n}{l}_j(x,y),\end{array}$$

and

$$\begin{array}{}\text{(21)}& {\mathit{\kappa }}_i={\mathit{\kappa }}_{i-\mathrm{1}}\left({\displaystyle \frac{a_{i+\mathrm{1}}(x_{i-\mathrm{1}}-x_i)+b_{i+\mathrm{1}}(y_{i-\mathrm{1}}-y_i)}{a_{i-\mathrm{1}}(x_i-x_{i-\mathrm{1}})+b_{i-\mathrm{1}}(y_i-y_{i-\mathrm{1}})}}\right);{\mathit{\kappa }}_{\mathrm{1}}=\mathrm{1};\end{array}$$

l_j(x,y) is the line equation for the j polygon edge such that

$$\begin{array}{}\text{(22)}& l_j(x,y)=\mathrm{1}-a_{j}x-b_{j}y,\end{array}$$

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

$$\begin{array}{}\text{(23)}& {\mathcal{W}}_i^{\left(n\right)}(x,y)={\displaystyle \frac{{\mathcal{P}}^{(n-\mathrm{2})}(x,y)}{{\mathcal{P}}^{(n-\mathrm{3})}(x,y)}},\end{array}$$

where 𝒫^(m)(x,y) is an m-degree polynomial in x,y. The integrals of the Wachspress basis function within a cell are performed using the eighth-order 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

$${\displaystyle }\begin{array}{rl}& w_{\mathrm{1}}=(\mathrm{1}-{\mathit{\xi }}_{\mathrm{1}})(\mathrm{1}-{\mathit{\xi }}_{\mathrm{2}})\\ & w_{\mathrm{2}}=(\mathrm{1}-{\mathit{\xi }}_{\mathrm{1}}){\mathit{\xi }}_{\mathrm{2}}\\ & w_{\mathrm{3}}={\mathit{\xi }}_{\mathrm{1}}(\mathrm{1}-{\mathit{\xi }}_{\mathrm{2}})\\ & w_{\mathrm{4}}={\mathit{\xi }}_{\mathrm{1}}{\mathit{\xi }}_{\mathrm{2}},\end{array}$$

where ξ₁ and ξ₂ are the transformed unit square coordinates of the cell (Hunke, 2001).

PWL basis functions divide the polygonal cell into sub-triangles and use a linear basis within each sub-triangle (Bailey et al., 2008). To divide the polygonal cell into sub-triangles, a point is chosen within the cell and sub-triangles formed using this point and two adjacent vertices. The central point in the cell, x_c, is chosen as

$$\begin{array}{}\text{(24)}& {\mathit{x}}_{\mathrm{c}}=\sum _i^{n_v}{\mathit{\alpha }}_i{\mathit{x}}_i,\end{array}$$

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.

Figure 2 Example basis functions for the (a–d) Wachspress and (e–h) PWL methods. The basis function is shown for the top vertex for a (a, e) square, (b, f) pentagonal, (c, g) hexagonal, and (d, h) heptagonal cell. Contour levels are drawn between 0.1 (dark blue) and 0.9 (red) at intervals of 0.1.

Using those basis functions to expand σ₁₁ (with basis index l), u, and v (with basis index m), Eq. (15) can be written as

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{D}_{\mathrm{1}}& =-\sum _k^{n_p}\underset{k}{\int }\\ & \left[\sum _l^{n_v}{\mathit{\sigma }}_{\mathrm{11}l}{\mathcal{W}}_l\cdot \sum _m^{n_v}\left(u_m{\displaystyle \frac{\partial {\mathcal{W}}_m}{\partial x}}-{\displaystyle \frac{\tan \mathit{\theta }}{r}}{v}_m{\mathcal{W}}_m\right)\right]\mathrm{d}A,\end{array}\end{array}$$

where the derivative with respect to x has been taken inside the summation. Rearranging, this becomes

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{D}_{\mathrm{1}}& =-\sum _k^{n_p}\sum _l^{n_v}\sum _m^{n_v}{\mathit{\sigma }}_{\mathrm{11}l}\\ & \times \left(u_m\underset{k}{\int }{\mathcal{W}}_l{\displaystyle \frac{\partial {\mathcal{W}}_m}{\partial x}}\mathrm{d}A-{\displaystyle \frac{\tan \mathit{\theta }}{r}}{v}_m\underset{k}{\int }{\mathcal{W}}_l{\mathcal{W}}_m\mathrm{d}A\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(27)}& {\mathcal{S}}_{lm}^x=\underset{k}{\int }{\mathcal{W}}_l{\displaystyle \frac{\partial {\mathcal{W}}_m}{\partial x}}\mathrm{d}A\end{array}$$

and

$$\begin{array}{}\text{(28)}& {\mathcal{T}}_{lm}=\underset{k}{\int }{\mathcal{W}}_l{\mathcal{W}}_m\mathrm{d}A,\end{array}$$

we have

$$\begin{array}{}\text{(29)}& D_{\mathrm{1}}=-\sum _k^{n_p}\sum _l^{n_v}\sum _m^{n_v}{\mathit{\sigma }}_{\mathrm{11}l}\left(u_m{\mathcal{S}}_{lm}^x-{\displaystyle \frac{\tan \mathit{\theta }}{r}}{v}_m{\mathcal{T}}_{lm}\right).\end{array}$$

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₁ to the components of the stress divergence tensor at that velocity point:

$$\begin{array}{}\text{(30)}& {\displaystyle }\begin{array}{rl}& (\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{u_j}^{D_{\mathrm{1}}}={\displaystyle \frac{\mathit{\delta }{D}_{\mathrm{1}}}{\mathit{\delta }{u}_j}}=-\sum _k^{n_p}\sum _l^{n_v}{\mathit{\sigma }}_{\mathrm{11}l}{\mathcal{S}}_{lj}^x\\ & (\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{v_j}^{D_{\mathrm{1}}}={\displaystyle \frac{\mathit{\delta }{D}_{\mathrm{1}}}{\mathit{\delta }{v}_j}}=\sum _k^{n_p}\sum _l^{n_v}{\mathit{\sigma }}_{\mathrm{11}l}{\displaystyle \frac{\tan \mathit{\theta }}{r}}{\mathcal{T}}_{lj}.\end{array}\end{array}$$

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₁, D₂, and D₃:

$$\begin{array}{}\text{(31)}& {\displaystyle }\begin{array}{rl}& (\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{u_j}=(\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{u_j}^{D_{\mathrm{1}}}+(\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{u_j}^{D_{\mathrm{2}}}+(\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{u_j}^{D_{\mathrm{3}}}\\ & (\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{v_j}=(\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{v_j}^{D_{\mathrm{1}}}+(\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{v_j}^{D_{\mathrm{2}}}+(\mathbf{\nabla }\cdot \mathit{\sigma }{)}_{v_j}^{D_{\mathrm{3}}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(32)}& {\displaystyle }\begin{array}{rl}& {\dot{\mathit{ϵ}}}_{\mathrm{11}}={\displaystyle \frac{\partial u}{\partial x}}-{\displaystyle \frac{v\tan \mathit{\theta }}{r}}\\ & {\dot{\mathit{ϵ}}}_{\mathrm{22}}={\displaystyle \frac{\partial v}{\partial y}}\\ & {\dot{\mathit{ϵ}}}_{\mathrm{12}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)+{\displaystyle \frac{u\tan \mathit{\theta }}{\mathrm{2}r}}.\end{array}\end{array}$$

The strain rate tensor at cell vertex l is then given by

$$\begin{array}{}\text{(33)}& {\displaystyle }\begin{array}{rl}& {\dot{\mathit{ϵ}}}_{\mathrm{11}l}=\sum _m^{n_v}{u}_m{\left.{\displaystyle \frac{\partial {\mathcal{W}}_m}{\partial x}}\right|}_l-{\displaystyle \frac{v_l\tan {\mathit{\theta }}_l}{r}}\\ & {\dot{\mathit{ϵ}}}_{\mathrm{22}l}=\sum _m^{n_v}{v}_m{\left.{\displaystyle \frac{\partial {\mathcal{W}}_m}{\partial y}}\right|}_l\\ & {\dot{\mathit{ϵ}}}_{\mathrm{12}l}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\sum _m^{n_v}{u}_m{\left.{\displaystyle \frac{\partial {\mathcal{W}}_m}{\partial y}}\right|}_l+\sum _m^{n_v}{v}_m{\left.{\displaystyle \frac{\partial {\mathcal{W}}_m}{\partial x}}\right|}_l\right)+{\displaystyle \frac{u_l\tan {\mathit{\theta }}_l}{\mathrm{2}r}}.\end{array}\end{array}$$

The derivatives of the basis functions are taken at cell vertex l. For the variational method pre-computed 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.\frac{\partial {\mathcal{W}}_m}{\partial y}\right|}_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 quasi-uniform 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

$$\begin{array}{}\text{(34)}& \begin{array}{rl}\sum _i^{n_d}\underset{i}{\int }(u{F}_u+v{F}_v)\mathrm{d}A& =\sum _i^{n_d}\left[\underset{i}{\int }\sum _l^{n_v}{u}_l{\mathcal{W}}_l\sum _m^{n_v}{F}_{um}{\mathcal{W}}_m\right.\\ & \left.+\underset{i}{\int }\sum _o^{n_v}{v}_o{\mathcal{W}}_o\sum _q^{n_v}{F}_{vq}{\mathcal{W}}_q\right]\mathrm{d}A.\end{array}\end{array}$$

Taking the variation in I with respect to the discretized velocity component u_j

$$\begin{array}{}\text{(35)}& \begin{array}{rl}\sum _i^{n_d}& \sum _m^{n_v}{F}_{um}\underset{i}{\int }{\mathcal{W}}_j{\mathcal{W}}_m\mathrm{d}A\\ & ={\displaystyle \frac{\partial }{\partial u_j}}D(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}}),\end{array}\end{array}$$

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

$$\begin{array}{}\text{(36)}& \begin{array}{rl}{F}_{uj}& \sum _i^{n_d}\underset{i}{\int }\sum _m^{n_v}{\mathcal{W}}_j{\mathcal{W}}_m\mathrm{d}A\\ & \simeq {\displaystyle \frac{\partial }{\partial u_j}}D(u_{\mathrm{1}},\mathrm{\dots },u_{n_d},v_{\mathrm{1}},\mathrm{\dots },v_{n_d},{\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_{n_{\mathit{\sigma }}}),\end{array}\end{array}$$

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

$$\begin{array}{}\text{(37)}& A_j^{\prime }=\underset{i}{\int }\sum _m^{n_v}{\mathcal{W}}_j{\mathcal{W}}_m\mathrm{d}A.\end{array}$$

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 Finite-volume 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 finite-volume scheme is that it solves the conservative form of the momentum equation and can handle non-smooth solution features (such as sharp fronts) consistently.

Figure 3 Contour integration lines used by the finite-volume stress divergence operator scheme. (a) Strain rate at cell centers (circle) are calculated from line integrals around primary mesh cells (solid line). (b) Stress divergence at cell vertices (square) are calculated from line integrals around the dual mesh cells (solid line). Directions of normal vectors used in the integrals are shown for both panels.

For the finite-volume 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

$$\begin{array}{}\text{(38)}& \underset{\mathrm{\Omega }}{\int }\left[\mathbf{\nabla }\mathit{v}\right]\partial \mathrm{\Omega }=\underset{S}{\oint }\left[\mathit{n}\otimes \mathit{v}\right]\partial S,\end{array}$$

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

$$\begin{array}{}\text{(39)}& \underset{\mathrm{\Omega }}{\int }\left[{\mathbf{\nabla }}_{\mathit{S}}\mathit{v}\right]\partial \mathrm{\Omega }=\underset{S}{\oint }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathit{n}\otimes \mathit{v}+\mathit{v}\otimes \mathit{n}\right]\partial S.\end{array}$$

The strain rate at a point is then obtained from the limit

$$\begin{array}{}\text{(40)}& \dot{\mathit{ϵ}}={\mathbf{\nabla }}_S\mathit{v}=\underset{A\to \mathrm{0}}{lim}{\displaystyle \frac{\mathrm{1}}{A}}\underset{S}{\oint }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathit{n}\otimes \mathit{v}+\mathit{v}\otimes \mathit{n}\right]\partial S,\end{array}$$

where the integral is around a closed loop, S, with area A and normal vector n, and v is the sea-ice 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 sea-ice velocity at a cell edge as the average of the values on the two vertices forming that edge projected onto the tangent plane:

$$\begin{array}{}\text{(41)}& {\dot{\mathit{ϵ}}}^{\prime }={\displaystyle \frac{\mathrm{1}}{A}}\sum _i^{n_{\mathrm{e}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\mathit{n}}_i\otimes {\mathit{v}}_i+{\mathit{v}}_i\otimes {\mathit{n}}_i\right]l_i.\end{array}$$

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)

$$\begin{array}{}\text{(42)}& {\displaystyle }\begin{array}{rl}& {\dot{\mathit{ϵ}}}_{\mathrm{11}}={\dot{\mathit{ϵ}}}_{\mathrm{11}}^{\prime }-{\displaystyle \frac{v\tan \mathit{\theta }}{r}}\\ & {\dot{\mathit{ϵ}}}_{\mathrm{22}}={\dot{\mathit{ϵ}}}_{\mathrm{22}}^{\prime }\\ & {\dot{\mathit{ϵ}}}_{\mathrm{12}}={\dot{\mathit{ϵ}}}_{\mathrm{12}}^{\prime }+{\displaystyle \frac{u\tan \mathit{\theta }}{\mathrm{2}r}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(43)}& \iint \mathbf{\nabla }\cdot \mathit{\sigma }\mathrm{d}A=\oint \left[\mathit{\sigma }\cdot \mathit{n}\right]\mathrm{d}l\end{array}$$

or

$$\begin{array}{}\text{(44)}& \mathbf{\nabla }\cdot \mathit{\sigma }=\underset{A\to \mathrm{0}}{lim}{\displaystyle \frac{\mathrm{1}}{A}}\oint \left[\mathit{\sigma }\cdot \mathit{n}\right]\mathrm{d}l\end{array}$$

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

$$\begin{array}{}\text{(45)}& (\mathbf{\nabla }\cdot \mathit{\sigma }{)}^{\prime }={\displaystyle \frac{\mathrm{1}}{A_{\mathrm{d}}}}\sum _i^{n_c}\left[{\mathit{\sigma }}_i\cdot {\mathit{n}}_i\right]l_i,\end{array}$$

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 second-order 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

$$\begin{array}{}\text{(46)}& {\displaystyle }\begin{array}{rl}& (\mathbf{\nabla }\cdot \mathit{\sigma }{)}_u=(\mathbf{\nabla }\cdot \mathit{\sigma }{)}_u^{\prime }-{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{12}}\tan \mathit{\theta }}{r}}\\ & (\mathbf{\nabla }\cdot \mathit{\sigma }{)}_v=(\mathbf{\nabla }\cdot \mathit{\sigma }{)}_v^{\prime }+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}\tan \mathit{\theta }}{r}}\end{array},\end{array}$$

where the components of σ are approximated as the average of the values on the dual mesh vertices.

2.3 Transport

To transport sea-ice fractional area and various tracers, MPAS-Seaice 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 quasi-regular hexagons and 12 pentagons.

For MPAS-Seaice the IR scheme was generalized to work on either the standard MPAS mesh (hexagons and other n-gons 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 MPAS-Seaice IR code is now mesh-agnostic, 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 MPAS-Seaice implementation. IR is designed to solve equations of the form

$$\begin{array}{}\text{(47)}& {\displaystyle \frac{\partial m}{\partial t}}=-\mathbf{\nabla }\mathbf{\cdot }\left(\mathit{u}m\right)\text{(48)}& {\displaystyle \frac{\partial \left(m{T}_{\mathrm{1}}\right)}{\partial t}}=-\mathbf{\nabla }\mathbf{\cdot }\left(\mathit{u}m{T}_{\mathrm{1}}\right),\text{(49)}& {\displaystyle \frac{\partial \left(m{T}_{\mathrm{1}}{T}_{\mathrm{2}}\right)}{\partial t}}=-\mathbf{\nabla }\mathbf{\cdot }\left(\mathit{u}m{T}_{\mathrm{1}}{T}_{\mathrm{2}}\right),\text{(50)}& {\displaystyle \frac{\partial \left(m{T}_{\mathrm{1}}{T}_{\mathrm{2}}{T}_{\mathrm{3}}\right)}{\partial t}}=-\mathbf{\nabla }\mathbf{\cdot }\left(\mathit{u}m{T}_{\mathrm{1}}{T}_{\mathrm{2}}{T}_{\mathrm{3}}\right),\end{array}$$

where $\mathit{u}=(u,v)$ is the horizontal velocity; m is mass or a mass-like field (such as density or fractional sea-ice concentration); and T₁, T₂, and T₃ 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 MPAS-Seaice, the fractional ice area in each thickness category is a mass-like 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)), MPAS-Seaice advects tracers up to type 3. Thus, the mass-like 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 sea-ice 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 second-order 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 MPAS-Seaice.

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,

$$\begin{array}{}\text{(51)}& {\displaystyle \frac{max\left(\right|\mathit{u}\left|\mathrm{\Delta }t\right)}{\mathrm{\Delta }x}}\le \mathrm{1},\end{array}$$

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:

1. 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.

2. 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.

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

4. 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, MPAS-Seaice 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 MPAS-Seaice 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 MPAS-Seaice scheme treats T₃ 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 MPAS-Seaice.

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 n-gons) as functions of $\mathit{r}=(x,y)$ in a cell-based 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{‾}}{a}$, should be recovered as the mean value when integrated over the cell:

$$\begin{array}{}\text{(52)}& \underset{A}{\int }a(x,y)\mathrm{d}A=\stackrel{\mathrm{‾}}{a}{A}_{\mathrm{C}},\end{array}$$

where A_C is the grid cell area. Equation (52) is satisfied if a(r) has the form

$$\begin{array}{}\text{(53)}& a\left(\mathit{r}\right)=\stackrel{\mathrm{‾}}{a}+{\mathit{\alpha }}_a\mathbf{\nabla }a\mathbf{\cdot }(\mathit{r}-\stackrel{\mathrm{‾}}{\mathit{r}}),\end{array}$$

where ∇a is a cell-centered gradient, α_a is a coefficient between 0 and 1 that enforces monotonicity, and $\stackrel{\mathrm{‾}}{\mathit{r}}$ is the cell centroid:

$$\begin{array}{}\text{(54)}& \stackrel{\mathrm{‾}}{\mathit{r}}={\displaystyle \frac{\mathrm{1}}{A_{\mathrm{C}}}}\underset{A}{\int }\mathit{r}\mathrm{d}A.\end{array}$$

Similarly, tracer means should be recovered when integrated over a cell:

$$\begin{array}{}\text{(55)}& {\displaystyle }\begin{array}{rl}\underset{A}{\int }a\left(\mathit{r}\right)T_{\mathrm{1}}\left(\mathit{r}\right)\mathrm{d}A& =\stackrel{\mathrm{‾}}{a}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}{A}_{\mathrm{C}},\\ \underset{A}{\int }a\left(\mathit{r}\right)T_{\mathrm{1}}\left(\mathit{r}\right)T_{\mathrm{2}}\left(\mathit{r}\right)\mathrm{d}A& =\stackrel{\mathrm{‾}}{a}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}{A}_{\mathrm{C}},\\ \underset{A}{\int }a\left(\mathit{r}\right)T_{\mathrm{1}}\left(\mathit{r}\right)T_{\mathrm{2}}\left(\mathit{r}\right)T_{\mathrm{3}}\left(\mathit{r}\right)\mathrm{d}A& =\stackrel{\mathrm{‾}}{a}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{3}}{A}_{\mathrm{C}}.\end{array}\end{array}$$

These equations are satisfied when the tracers are reconstructed as

$$\begin{array}{}\text{(56)}& {\displaystyle }\begin{array}{r}{T}_{\mathrm{1}}\left(\mathit{r}\right)={\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}+{\mathit{\alpha }}_{T\mathrm{1}}\mathbf{\nabla }{T}_{\mathrm{1}}\cdot (\mathit{r}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{1}}),\\ T_{\mathrm{2}}\left(\mathit{r}\right)={\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}+{\mathit{\alpha }}_{T\mathrm{2}}\mathbf{\nabla }{T}_{\mathrm{2}}\cdot (\mathit{r}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{2}}),\\ T_{\mathrm{3}}\left(\mathit{r}\right)={\stackrel{\mathrm{‾}}{T}}_{\mathrm{3}}+{\mathit{\alpha }}_{T\mathrm{3}}\mathbf{\nabla }{T}_{\mathrm{3}}\cdot (\mathit{r}-{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{3}}),\end{array}\end{array}$$

where the tracer barycenter coordinates ${\stackrel{\mathrm{̃}}{\mathit{r}}}_n$ are given by

$$\begin{array}{}\text{(57)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{1}}& ={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{a}{A}_{\mathrm{C}}}}\underset{A}{\int }\mathit{r}a\mathrm{d}A,\\ {\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{2}}& ={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{a}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}{A}_{\mathrm{C}}}}\underset{A}{\int }\mathit{r}a{T}_{\mathrm{1}}\mathrm{d}A,\\ {\stackrel{\mathrm{̃}}{\mathit{r}}}_{\mathrm{3}}& ={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{a}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{2}}{A}_{\mathrm{C}}}}\underset{A}{\int }\mathit{r}a{T}_{\mathrm{1}}{T}_{\mathrm{2}}\mathrm{d}A.\end{array}\end{array}$$

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 second-order accurate in space. When α=0, the reconstruction reduces locally to first-order.

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₁ and V₂. Each edge is oriented such that one adjacent cell (C₁) is defined to lie in the left half-plane and the other (C₂) in the right half-plane. The departure points D₁ and D₂ are found by projecting velocities backward from V₁ and V₂. The shaded departure region is a quadrilateral containing all the ice transported across the edge in one time step. In addition to C₁ and C₂, the departure region can include side cells C₃ and C₄. The side cells share edges E₁ to E₄ and vertices V₃ to V₆ with the central cells C₁ and C₂.

Figure 4 (a) Schematic showing transport across a cell edge on a standard MPAS mesh with three edges meeting at each vertex. The letters C, E, and V denote cell centers, edges, and vertices, respectively. Vertices are also shown with a square marker. Backward trajectories are shown as red arrows, ending in departure points D₁ and D₂. These backward trajectories define the departure region, which is shaded in blue. (b) Schematic showing transport across a cell edge on a quadrilateral MPAS mesh with four edges meeting at each vertex. Shaded regions are the departure regions.

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₁ and V₂. These coordinates are pre-computed at initialization. During each time step, departure triangles are found by locating D₁ and D₂ in this coordinate system and then looping through the edges to identify any intersections of line segment $\stackrel{\mathrm{‾}}{D_{\mathrm{12}}}$ (i.e., the segment joining D₁ and D₂) with the various edges. If $\stackrel{\mathrm{‾}}{D_{\mathrm{12}}}$ intersects the main edge, then the departure region consists of two triangles (one each in C₁ and C₂) rather than a quadrilateral (as shown in Fig. 4b). If $\stackrel{\mathrm{‾}}{D_{\mathrm{12}}}$ intersects any of edges E₁ to E₄, 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₁ and C₂ 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₁ and V₂), departure points (D₁ and D₂), and intersection points (points where $\stackrel{\mathrm{‾}}{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₅ and C₆), edges (E₅ and E₆), and vertices (V₇ and V₈). 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 mesh-specific 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 edge-based coordinates, they are transformed to cell-based 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

$$\begin{array}{}\text{(58)}& A_T={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left|(x_{\mathrm{2}}-x_{\mathrm{1}})(y_{\mathrm{3}}-y_{\mathrm{1}})-(y_{\mathrm{2}}-y_{\mathrm{1}})(x_{\mathrm{3}}-x_{\mathrm{1}})\right|,\end{array}$$

where the (x_i,y_i) are the triangle vertices.

2.3.3 Integrating the transport

Next, ice area and area-tracer 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 area-tracer 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

$$\begin{array}{}\text{(59)}& I={\displaystyle \frac{A_T}{\mathrm{3}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}f\left({{\mathit{x}}^{\mathbf{\prime }}}_i\right).\end{array}$$

The quadrature points are located at ${{\mathit{x}}^{\mathbf{\prime }}}_i=({\mathit{x}}_{\mathrm{0}}+{\mathit{x}}_i)/\mathrm{2}$, where x₀ 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:

$$\begin{array}{}\text{(60)}& I=A_T\left[w_{\mathrm{1}}\sum _{i=\mathrm{1}}^{\mathrm{3}}f\left({\mathit{x}}_{\mathrm{1}i}\right)+w_{\mathrm{2}}\sum _{i=\mathrm{1}}^{\mathrm{3}}f\left({\mathit{x}}_{\mathrm{2}i}\right)\right],\end{array}$$

where x_1i and x_2i are two sets of three quadrature points, arranged symmetrically on the three medians of the triangle, and w₁ and w₂ 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 area-tracer 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 area-tracer products as Δ(aT₁)_k, Δ(aT₁T₂)_k, and Δ(aT₁T₂T₃)_k. The new ice area at time n+1 is given by

$$\begin{array}{}\text{(61)}& a^{n+\mathrm{1}}=a^n+{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{C}}}}\sum _k\pm \mathrm{\Delta }{a}_k,\end{array}$$

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

$$\begin{array}{}\text{(62)}& {\displaystyle }\begin{array}{rl}{T}_{\mathrm{1}}^{n+\mathrm{1}}& ={\displaystyle \frac{a^n{T}_{\mathrm{1}}^n+\frac{\mathrm{1}}{A_{\mathrm{C}}}\sum _k\pm \mathrm{\Delta }(a{T}_{\mathrm{1}}{)}_k}{a^{n+\mathrm{1}}}},\\ T_{\mathrm{2}}^{n+\mathrm{1}}& ={\displaystyle \frac{a^n{T}_{\mathrm{1}}^n{T}_{\mathrm{2}}^n+\frac{\mathrm{1}}{A_{\mathrm{C}}}\sum _k\pm \mathrm{\Delta }(a{T}_{\mathrm{1}}{T}_{\mathrm{2}}{)}_k}{a^{n+\mathrm{1}}{T}_{\mathrm{1}}^{n+\mathrm{1}}}},\\ T_{\mathrm{3}}^{n+\mathrm{1}}& ={\displaystyle \frac{a^n{T}_{\mathrm{1}}^n{T}_{\mathrm{2}}^n{T}_{\mathrm{3}}^n+\frac{\mathrm{1}}{A_{\mathrm{C}}}\sum _k\pm \mathrm{\Delta }(a{T}_{\mathrm{1}}{T}_{\mathrm{2}}{T}_{\mathrm{3}}{)}_k}{a^{n+\mathrm{1}}{T}_{\mathrm{1}}^{n+\mathrm{1}}{T}_{\mathrm{2}}^{n+\mathrm{1}}}}.\end{array}\end{array}$$

Dukowicz and Baumgardner (2000) showed that Eq. (62) satisfies tracer monotonicity since the new-time tracer values are area-weighted averages of old-time 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 melt-pond parameterizations (Flocco et al., 2010; Holland et al., 2012; Hunke et al., 2013), a delta-Eddington 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 sea-ice BGC (Jeffery and Hunke, 2014; Jeffery et al., 2016, 2020). To include these developments in MPAS-Seaice, 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). MPAS-Seaice uses a version of the column package that was forked from CICE version 5.1.2.

3 Model validation

3.1 Velocity solver

To validate our implementation of the MPAS-Seaice 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

$$\begin{array}{}\text{(63)}& f(x,y)=\sin \left(\mathrm{2}A\mathit{\pi }x\right)\sin \left(\mathrm{2}A\mathit{\pi }y\right),\end{array}$$

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₂ 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₂ 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₂ norms for the variational derivatives and strains are given by

$$\begin{array}{}\text{(64)}& L_{\mathrm{2}}=\sqrt{{\displaystyle \frac{{\sum }_i{\int }_i({\sum }_l{\mathcal{W}}_l{f}_{il}-\widehat{f}(x,y){)}^{\mathrm{2}}\mathrm{d}A}{{\sum }_i{\int }_i\widehat{f}(x,y{)}^{\mathrm{2}}\mathrm{d}A}}},\end{array}$$

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 sub-triangles and using the eighth-order 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 finite-volume strain operators we perform line integrals of the square of the error around the dual mesh cell surrounding primary cell points. The L₂ norms for the finite-volume derivatives and strains are given by

$$\begin{array}{}\text{(65)}& L_{\mathrm{2}}=\sqrt{{\displaystyle \frac{{\sum }_i{\sum }_j{l}_j{\int }_{\mathrm{0}}^{\mathrm{1}}\left(f_j\right(\mathit{\chi })-{\widehat{f}}_j(\mathit{\chi }){)}^{\mathrm{2}}\mathrm{d}\mathit{\chi }}{{\sum }_i{\sum }_j{l}_j{\int }_{\mathrm{0}}^{\mathrm{1}}{\widehat{f}}_j(\mathit{\chi }{)}^{\mathrm{2}}\mathrm{d}\mathit{\chi }}}},\end{array}$$

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 seventh-order 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₂ error norm by a factor of 2 for first-order accurate schemes and a factor of 4 for second-order accurate schemes. Figure 5 also shows idealized first- and second-order scaling gradients as dotted lines. Evident from the figure is that the spatial derivative operators for the variational methods (Wachspress and PWL) are only first-order accurate for both square and hexagonal meshes, while the finite-volume derivative operators are second-order 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 half-plane with respect to the derivative point, effectively creating a one-sided stencil for the operator. Velocity values for the finite-volume 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 one-sided values of the derivatives calculated for a single vertex to create a multi-sided stenciled operator from the one-sided variational ones. Figure 5 shows the error scaling for this averaging method for the Wachspress (“Wachs. Avg”) and PWL (“PWL Avg”) basis functions. Second-order convergence is achieved in the y direction with this averaging, but interestingly, only averaging the PWL basis function results in second-order convergence in the x direction, while the Wachspress basis function shows only first-order 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 finite-volume method can also be averaged to the cell vertices to create another operator (“FV Avg” in Fig. 5), creating a second-order accurate method, although this averaging increases the error relative to the original finite-volume scheme.

Figure 5 Scaling of L₂ error norm against grid resolution for the model calculation of the gradient of an analytical function in the (a, c) x and (b, d) y directions for a planar test case with a regular (a, b) square cell and (c, d) hexagonal cell mesh. Grid resolution is the mean distance between cell centers. Ideal first- and second-order scaling gradients are also shown as dotted lines.

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{ϵ}}}_{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₂ error norm of the stress divergence operators with grid resolution. For the L₂ error norm calculated for the stress divergence operators in this section and in Sect. 3.1.2, we use

$$\begin{array}{}\text{(66)}& L_{\mathrm{2}}=\sqrt{{\displaystyle \frac{\sum A_i(f_i-\widehat{f_i}{)}^{\mathrm{2}}}{\sum A_i{\widehat{f_i}}^{\mathrm{2}}}}},\end{array}$$

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 finite-volume methods show second-order convergence of errors for the square cell mesh, with the finite-volume scheme showing significantly lower errors. With the PWL basis functions the variational method shows second-order convergence, while with Wachspress basis functions the variational method shows a varying order of convergence with near-second-order convergence at low resolution but with the order decreasing to first with higher resolutions. This suggests there is a small source of first-order error when using the Wachspress basis functions. For the hexagonal cell mesh, the finite-volume method shows only first-order convergence with significantly higher error than the variational method. This is because the line integral around vertices for the finite-volume 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 finite-volume integrals around squares or hexagons, opposite edges of the polygon lead to cancellation of first-order 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.

Figure 6 Scaling of L₂ error norm against grid resolution for the model calculation of the (a) u and (b) v components of stress divergence from an analytical strain rate tensor field for a planar test case with a regular mesh. Grid resolution is the mean distance between cell centers. Solid lines denote hexagonal cell meshes, whereas dashed lines signify square cell meshes. Ideal first- and second-order scaling gradients are also shown as dotted lines.

In summary, for regular planar meshes, for the gradient operators, the finite-volume and averaged PWL and averaged finite-volume 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 finite-volume scheme, while for square cell meshes the order of convergence is the same between the variational and finite-volume schemes, with the finite-volume 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{ϵ}}}_{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.

Figure 7 Orthographic projection of one hemisphere of the unit sphere test case described in Sect. 3.1.2. Panels (a) and (b) show the analytical forms of velocity used in the numerical schemes and used to calculate the desired analytical strain rates (shown in e–g) and the desired analytical stress divergence (shown in c and d). Panels (h), (i), (l), (m), (p), and (q) show the error in calculation of the ${\dot{\mathit{ϵ}}}_{\mathrm{11}}$ component of the strain rate for the (h, i) Wachspress variational scheme, the (l, m) PWL variational scheme, and the (p, q) finite-volume (FV) scheme. Panels (j), (k), (n), (o), (r), and (s) show the error in calculation of the u component of the stress divergence for the (j, k) Wachspress variational scheme, the (n, o) PWL variational scheme, and the (r, s) finite-volume (FV) scheme. Detail is shown around one of the pentagonal mesh cells for each scheme in (i), (k), (m), (o), (q), and (s). The error is the difference between the model and analytical values.

Errors generated for the strain rate component ${\dot{\mathit{ϵ}}}_{\mathrm{11}}$ for this test case are displayed in Fig. 7h, l, and p. These show that the error is significantly lower for the finite-volume 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 finite-volume scheme is enhanced around the 12 pentagonal cells found in the quasi-uniform 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 finite-volume scheme. The scaling of L₂ error norm with grid resolution is shown in Fig. 9a. L₂ 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 first-order convergence, and the finite-volume scheme shows second-order 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 finite-volume operator averaging show similar error scaling characteristics with close-to-second-order 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 first-order convergence for all strain rate operator methods. As for the hexagonal cell planar mesh, for both the L₂ and L_∞ norms, averaging the finite-volume scheme increases the error.

Figure 8 Frequency of grid cell error in the calculation of strain rate (a) and stress divergence (b) from an analytical velocity field for a spherical test case. The error is given as the absolute value of the difference between model and analytical fields. The results are plotted as histograms for the two variational schemes (using Wachspress and PWL basis functions and the alternative area formulation) and the finite-volume (FV) scheme.

Figure 9 Scaling of L₂ (a, b) and L_∞ (c, d) error norm against grid resolution for the model calculation of (a, c) the ${\dot{\mathit{ϵ}}}_{\mathrm{11}}$ component of the strain rate tensor from an analytical velocity field and (b, d) the u component of the stress divergence from an analytical strain rate tensor field on the unit sphere for a spherical test case. Grid resolution is the mean distance between cell centers. Only mesh cells where $\mathit{\theta }>\mathrm{20}{}^{\circ }$ are included. Ideal first- and second-order scaling gradients are also shown as dotted lines.

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 finite-volume scheme. Away from the pentagonal cells, the variational methods display smaller errors than the finite-volume scheme but show a significant enhancement of error at the pentagonal cells, something not found with the finite-volume 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₂ error norm against grid resolution is shown in Fig. 9b, where all methods display first-order 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 finite-volume method shows linear convergence for all grid resolutions.

Figure 10 Orthographic projection of one hemisphere of the stress divergence operator error for the unit sphere test case. Shown are results for the variational schemes using the alternate area denominator formulation. Panels (a) and (b) show results for the Wachspress variational scheme, while (c) and (d) show results for the PWL variational scheme. Detail is shown around one of the pentagonal mesh cells for each scheme in (b), (d).

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 MPAS-Seaice support quadrilateral grids, direct comparisons can be made between MPAS-Seaice and CICE. For idealized planar test cases it is possible to set up MPAS-Seaice 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 MPAS-Seaice 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

$$\begin{array}{}\text{(67)}& {\displaystyle }{u}_{\mathrm{a}}=\mathrm{5}-\mathrm{3}\sin {\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{L_x}}\sin {\displaystyle \frac{\mathit{\pi }y}{L_y}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\text{(68)}& {\displaystyle }{v}_{\mathrm{a}}=\mathrm{5}-\mathrm{3}\sin {\displaystyle \frac{\mathrm{2}\mathit{\pi }y}{L_y}}\sin {\displaystyle \frac{\mathit{\pi }x}{L_x}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}},\end{array}$$

while the ocean currents have the form

$$\begin{array}{}\text{(69)}& {\displaystyle }{u}_{\mathrm{o}}=\mathrm{0.1}{\displaystyle \frac{\mathrm{2}y-L_y}{L_y}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\text{(70)}& {\displaystyle }{v}_{\mathrm{o}}=-\mathrm{0.1}{\displaystyle \frac{\mathrm{2}x-L_x}{L_x}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}},\end{array}$$

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. Sea-ice 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 MPAS-Seaice. In this comparison MPAS-Seaice 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 sea-ice stresses are present (see Fig. 12e). The figure also shows that the three MPAS-Seaice 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 MPAS-Seaice using a regular hexagonal mesh. Here differences between the finite-volume 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 right-hand corner of the domain.

Figure 11 Forcing fields of (a) wind velocity and (b) ocean velocity used in a square domain test case.

Figure 12 (a) The u velocity component of CICE for the square domain test case using a quadrilateral mesh. (b, c, d) Difference ($\mathrm{\Delta }u/u_{max}$) between the u component of MPAS-Seaice and CICE in the square domain test case for the MPAS-Seaice Wachspress variational, PWL variational, and finite-volume (FV) schemes, respectively, and using a quadrilateral mesh. (e)–(h) As (a)–(d) but for the u component of the divergence of internal ice stress.

Figure 13 (a) The u velocity component of MPAS-Seaice with the Wachspress variational scheme for the square domain test case using hexagonal elements. (b, c) Difference ($\mathrm{\Delta }u/u_{max}$) between the u component of MPAS-Seaice using the Wachspress variational scheme and MPAS-Seaice using the PWL variational and finite-volume (FV) schemes, respectively, in the square domain test case using hexagonal cells. (d)–(f) As (a)–(c) but for the u component of the divergence of internal ice stress.

Figure 14 Principal stress components for the square domain test case and (dotted) the EVP yield curve. (a) CICE. (b–d) MPAS-Seaice on a quadrilateral mesh and using the Wachspress variational, PWL variational, and finite-volume (FV) schemes, respectively. (e–g) MPAS-Seaice on a hexagonal mesh and using the Wachspress variational, PWL variational, and finite-volume (FV) schemes, respectively.

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}}\cos \mathit{\theta },\mathrm{0})$, where $u_{\mathrm{0}}=\left(\mathrm{2}\mathit{\pi }R\right)/\left(\mathrm{12}\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 first-order upwind scheme to demonstrate the improvements in numerical diffusion gained by increasing the order of transport scheme. For both the IR scheme and the first-order 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 non-convergent, 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₂ 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 close-to-second-order accuracy (indicated by the dotted diagonal line) for the cosine bell and converges slightly below first-order 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 along-motion 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 upstream-based method.

Figure 15 Comparison of incremental remapping to a first-order upwind scheme for advection around the sphere. Panel (a) shows a cross section of a circular region of sea ice whose center lies on the Equator, with radius $R/\mathrm{3}$ (where R is the Earth's radius) and initial sea-ice concentration given by a cosine bell. Panel (b) shows a cross section of a slotted cylinder whose center lies on the Equator, with radius $R/\mathrm{2}$. The grid resolution is 120 km. The exact solution (which corresponds to the initial condition) is shown by a solid line, IR by long dashes, and upwind by short dashes.

Figure 16 Initial (a, d) and final (b, c, e, f) ice concentration contours at 120 km resolution for the cosine bell (a–c) and slotted cylinder (d–f) advection test cases. Results for the upwind advection scheme are shown in (b) and (e), while results from the incremental remapping scheme are shown in (c) and (f). Contours are at levels 0.05, 0.1, 0.3, 0.5, 0.7, 0.9, and 0.95 unless otherwise shown. The direction of transport is to the right.

Figure 17 Scaling of the L₂ error norm with grid resolution for the cosine bell (CB) and slotted cylinder (SC) advection tests shown in Fig. 15. IR scaling is shown by the solid lines, upwind by long-dashed lines. Theoretical linear and quadratic scaling are shown by short-dashed lines.

3.3 Column physics

To validate the column physics in MPAS-Seaice we make use of the fact that CICE and MPAS-Seaice share identical code. CICE and MPAS-Seaice were run with identical forcing and with dynamics disabled. Results from the two models were bit-for-bit identical, indicating a correct implementation of the column physics in MPAS-Seaice.

3.4 Global simulations

To validate the full MPAS-Seaice 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 MPAS-Seaice 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 6-hourly frequency are taken from the Coordinated Ocean-ice Reference Experiments (CORE) Corrected Inter-Annual 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, sea-surface 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 MPAS-Seaice 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 sea-ice extent between the MPAS-Seaice and CICE models and observational values for the years 1988 to 2007 inclusive. The observational sea-ice 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 sea-ice extent well represented in both hemispheres. The largest discrepancy occurs in the Southern Hemisphere, where austral summertime sea-ice extent is too low in both models.

Figure 18 Sea-ice extent (area with ice concentration greater than 15 %) by year for the Northern Hemisphere (a–c) and Southern Hemisphere (d–f), for MPAS-Seaice (a, d), CICE (b, e), and SSMI satellite observations (c, f). Northern Hemisphere total sea-ice volume by year for (g) MPAS-Seaice, (h) CICE, and (i) PIOMAS.

Figure 19 Spatial climatological maps for 1988 to 2007 of sea-ice concentration from (a–d) SSMI satellite observations processed with the NASATeam algorithm and (e–h) MPAS-Seaice. Differences between (i–l) MPAS-Seaice and (m–p) CICE ice concentration and SSMI observations. (q–t) Differences between MPAS-Seaice and CICE ice concentration. (a, e, i, m, q) Northern Hemisphere winter: January, February, and March. (b, f, j, n, r) Southern Hemisphere summer: December, January, and February. (c, g, k, o, s) Northern Hemisphere summer: July, August, and September. (d, h, l, p, t) Southern Hemisphere winter: June, July, and August.

Figure 19 shows a similar agreement, comparing sea-ice 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 sea-ice extent is well reproduced. In the Southern Hemisphere the sea-ice 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 MPAS-Seaice and CICE are explained by a number of differences in implementation between these models. Firstly, since MPAS-Seaice 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, MPAS-Seaice 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 sea-ice volume for the Northern Hemisphere is compared between the models and the Pan-Arctic 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 sea-ice thickness for the Northern Hemisphere in summer, autumn, and winter compared to observations of sea-ice 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 sea-ice 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 MPAS-Seaice and CICE have excess sea-ice thickness compared to ICESat observations in the Beaufort Sea and western Arctic basin and a deficit of sea-ice thickness in the Eurasian basin.

Figure 20 Spatial climatological maps of sea-ice thickness from ICESat (a–c) and MPAS-Seaice (d–f). Differences between (g–i) MPAS-Seaice and (j–l) CICE ice thickness and ICESat observations. (m–o) Differences between MPAS-Seaice and CICE ice thickness. (a, d, g, j, m) Northern Hemisphere summer. (b, e, h, k, n) Northern Hemisphere autumn. (c, f, i, l, o) Northern Hemisphere winter.

Finally, we perform simulations with MPAS-Seaice on a quasi-uniform 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 southern-hemispheric total sea-ice 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 quasi-uniform 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 basin-scale simulations. Compared to the Wachspress variational scheme, the PWL variational and finite-volume schemes have less effect on ice extent and a similar effect on ice volume as compared to the difference between the quasi-uniform SCVT and 1^∘ quadrilateral meshes. Differences are also strongly seasonal for these changes in operator, especially for total volume.

Figure 21 Difference (%) in hemispheric sea-ice extent (a, c) and volume (b, d) between global simulations using the Wachspress variational scheme with a quasi-uniform 30 km SCVT mesh and several other simulations. Results are shown for the Northern Hemisphere (a, b) and Southern Hemisphere (b, d).

4 Computational performance

Through the MPAS framework, MPAS-Seaice incorporates code parallelization through domain decomposition and message passing with the Message Passing Interface (MPI) library. To assess the computational performance of MPAS-Seaice when run on multiple processors, we perform a strong-scaling 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 time-varying 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 quasi-uniform SCVT mesh with cell-to-cell 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 MPAS-Seaice 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 mid-ocean. 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.

Figure 22 Strong-scaling performance characteristics of MPAS-Seaice for a global QU60km mesh (global) and a QU60km mesh with equatorial cells removed (polar). Simulated years per day of wall clock time (SYPD) is plotted against (a) the number of compute cores and (b) cells per core. Perfect strong scaling would be linear for panel (a) and inversely linear for panel (b). Three partitioning methods are employed as described in the main text: simple, region, and weight.

Load balancing is an issue with sea-ice 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 high-latitude cells will take longer to compute a time step than processors containing equatorial cells, which will have to wait periodically for high-latitude 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 high-latitude 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 high-latitude cells a larger weight than equatorial cells means fewer cells are included in high-latitude 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 sub-linear strong-scaling performance for MPAS-Seaice 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 near-linear 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. Sub-linear 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 Finite-Element/volumE Sea ice-Ocean Model (FESOM) (Koldunov et al., 2019).

Comparison of computational performance between MPAS-Seaice and CICE is non-trivial. While MPAS-Seaice 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 MPAS-Seaice 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 MPAS-Seaice and CICE we compare the SYPD achieved on Anvil for 10 d of the simulation described in Sect. 3.4 with both MPAS-Seaice and CICE using a fully global 1^∘ quadrilateral grid. Simulation throughput for MPAS-Seaice 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 MPAS-Seaice achieved approximately 70 % of the CICE throughput. The better computational performance of CICE is expected since the unstructured mesh in MPAS-Seaice necessitates less efficient memory access patterns. As a percentage of CICE model performance, the MPAS-Seaice velocity solver displayed the poorest throughput and the column physics the most competitive throughput.

Table 1 Model throughput (in SYPD) for various code sections for MPAS-Seaice and CICE. Results from 32 processors on Anvil with MPAS-Seaice using the “region” partitioning method and CICE using the “slenderX2” partitioning method.

5 Conclusions

We have described a new sea-ice model, MPAS-Seaice, 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 sea-ice 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 MPAS-Seaice 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 one-sided stencil of the variational strain rate operators results in poor error characteristics. The finite-volume 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 finite-volume stress divergence operator. This suggests that the variational scheme could be improved by modifying its strain rate operator to have a two-sided stencil. We investigated several averaging techniques to implement a two-sided stencil for the variational scheme, which resulted in improved error characteristics for the variational strain rate operator. For basin-scale sea-ice simulations, however, these alternate operators had only a small effect on simulation results.

MPAS-Seaice and CICE share the sophisticated suite of column physics and BGC originally developed in CICE, again allowing the rapid development of MPAS-Seaice. Global simulations with realistic forcing have validated MPAS-Seaice against similar simulations with CICE and against observations for sea-ice concentration, extent, and volume. MPAS-Seaice 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 sea-ice results from E3SM simulations. MPAS-Seaice shows power-law strong-scaling 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 MPAS-Seaice on variable-resolution meshes and examine more recent metrics for evaluating sea-ice dynamics, such as new statistical metrics of linear kinematic features (e.g., Hutter and Losch, 2020). Potential challenges with variable-resolution meshes include assessing the resolution invariance of sea-ice rheologies and developing resolution-invariant versions, efficiently using future heterogeneous computing architectures, as well as generating efficient domain partitions of highly variable-resolution meshes.