the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# CICE on a C-grid: new momentum, stress, and transport schemes for CICEv6.5

### Jean-François Lemieux

### William H. Lipscomb

### Anthony Craig

### David A. Bailey

### Elizabeth C. Hunke

### Philippe Blain

### Till A. S. Rasmussen

### Mats Bentsen

### Frédéric Dupont

### David Hebert

### Richard Allard

This article presents the C-grid implementation of the CICE sea ice model, including the C-grid discretization of the momentum equation, the boundary conditions (BCs), and the modifications to the code required to use the incremental remapping transport scheme. To validate the new C-grid implementation, many numerical experiments were conducted and compared to the B-grid solutions. In idealized experiments, the standard advection method (incremental remapping with C-grid velocities interpolated to the cell corners) leads to a checkerboard pattern. A modal analysis demonstrates that this computational noise originates from the spatial averaging of C-grid velocities at corners. The checkerboard pattern can be eliminated by adjusting the departure regions to match the divergence obtained from the solution of the momentum equation. We refer to this novel approach as the edge flux adjustment (EFA) method. The C-grid discretization with edge flux adjustment allows for transport in channels that are one grid cell wide – a capability that is not possible with the B-grid discretization nor with the C-grid and standard remapping advection. Simulation results match the predicted values of a novel analytical solution for one-grid-cell-wide channels.

- Article
(3082 KB) - Full-text XML
- BibTeX
- EndNote

CICE (Hunke et al., 2023) is a dynamic and thermodynamic sea ice model used for a variety of applications, such as climate modeling (e.g., DeRepentigny et al., 2020), sub-seasonal sea ice forecasting (e.g., Barton et al., 2021), and short-term sea ice forecasting (e.g., Smith et al., 2021). Since 2017, the model has been developed by the CICE Consortium, a group of institutions from the USA, Canada, Denmark, and Poland.

Earlier versions of CICE used the Arakawa B-grid (Arakawa and Lamb, 1977, i.e., the horizontal velocity components *u* and *v* are co-located at cell corners) for the spatial discretization. This co-location of *u* and *v* simplifies the treatment of the Coriolis term and of the off-diagonal part of the water stress term. Another interesting characteristic of the B-grid is the straightforward implementation of no-slip and no-outflow boundary conditions.

Recently, many CICE users have requested a C-grid capability, in which the *u* component is defined on east and west cell edges and the *v* component on north and south edges. The C-grid also has appealing characteristics. First, it allows for straightforward coupling with C-grid ocean models (e.g., NEMO and HYCOM; Madec and the NEMO System Team, 2008; Metzger et al., 2014) and atmospheric models (e.g., GEM; McTaggart-Cowan et al., 2019). Second, as opposed to a B-grid, it can represent transport along channels that are only one grid cell wide. Finally, the C-grid discretization better represents inertial–plastic compressive waves (Bouillon et al., 2009). For these reasons, members of the CICE Consortium decided to implement a finite-difference C-grid capability in CICE, which is presented here.

Other widely used continuum-based sea ice models with a C-grid discretization include the Sea Ice modelling Integrated Initiative (SI^{3}; Vancoppenolle et al., 2023), the sea ice component of the Massachusetts Institute of Technology general circulation model (MITgcm; Losch et al., 2010), and the Sea Ice Simulator version 2 (SIS2; Adcroft et al., 2019). In terms of dynamics, these sea ice models have many similarities (e.g., they include the elastic–viscous–plastic scheme) but differ in the way sea ice (and snow) transport is implemented. SI^{3} uses the Prather (1986) transport scheme, which is based on the conservation of second-order moments. The default approach for transport in the MITgcm is a second-order scheme with a superbee flux limiter. Other options for transport in the MITgcm include upwind schemes and a fifth-order scheme with weighted essentially non-oscillatory (WENO) limiters. The default transport method for sea ice and snow in SIS2 is a non-directionally split first-order upwind scheme, and SIS2 also includes options for directionally split piecewise parabolic, directionally split piecewise linear, and directionally split piecewise constant methods.

When we started working on the C-grid discretization, CICE already included a first-order upwind transport scheme. Upwind transport is naturally suited for a C-grid discretization since edge velocities can be directly used to calculate fluxes on the four sides of a quadrilateral grid cell. Although upwind transport has some desirable characteristics (e.g., it is conservative, stable, computationally efficient, and sign-preserving), it is very diffusive. Sharp features such as the ice edge are quickly smeared out when simulating transport with an upwind scheme.

The default advection algorithm for the B-grid in CICE is the remapping transport scheme (Dukowicz and Baumgardner, 2000; Lipscomb and Hunke, 2004). Although remapping is fundamentally a B-grid transport scheme, we decided to adapt it for the C-grid discretization. An important reason for doing this is to be able to reuse a significant part of the code and therefore speed up the implementation. The attractive properties of remapping also motivated our choice. Like upwind, incremental remapping is conservative, stable, and sign-preserving. Although the geometric calculations (computing departure regions for each cell edge) are relatively expensive, the method scales well when there are many tracers, as is the case in CICE. Also, with a linear reconstruction of scalar fields, remapping is much less diffusive than the upwind scheme and therefore better preserves sharp features (Lipscomb and Hunke, 2004).

In our first implementation of remapping for the C-grid, the C-grid velocities on cell edges were interpolated to the corners, and remapping was then applied the same way as for the B-grid. This approach, however, does not simulate transport in one-grid-cell-wide channels. Also, some idealized tests showed the presence of numerical noise (a checkerboard pattern) in fields such as sea ice concentration.

We therefore introduce a novel method for adapting the remapping transport scheme for the C-grid discretization. We refer to this approach as the edge flux adjustment (EFA) method. The EFA method is necessary to eliminate the checkerboard numerical noise. Another crucial advantage of the EFA method is that it allows for transport in one-grid-cell-wide channels – a capability that is not possible with the B-grid discretization. As such, our implementation of remapping with the EFA method should be seen as a hybrid B- and C-grid transport scheme.

The main contribution of this work is (1) the introduction of the novel EFA method for the remapping transport scheme. Other contributions include (2) a detailed description and validation of the CICE C-grid spatial discretization, including the formulation of boundary conditions, (3) the derivation of analytical solutions for one-grid-cell-wide channels, and (4) the interpretation of the checkerboard pattern associated with the interpolation of C-grid velocities when using the standard remapping scheme.

This article is structured as follows. Section 2 introduces the momentum and stress equations (Sect. 2.1) and the equations for transport (Sect. 2.2). Section 3 briefly introduces the C-grid spatial and temporal discretizations; more details can be found in Appendix A. Section 4 describes our initial implementation of incremental remapping for the C-grid along with its weaknesses. These weaknesses are corrected with the edge flux adjustment method, which is explained in Sect. 5. Section 6 gives an overview of the different tests used to validate the new C-grid discretization. Concluding remarks and future work are given in Sect. 7. Appendix B presents a modal analysis of the remapping checkerboard pattern. Appendix C describes some modifications to improve the robustness of the remapping method. Finally, Appendix D introduces a novel analytical solution for one-grid-cell-wide channels.

## 2.1 Momentum equation and rheology

The two-dimensional sea ice momentum equation is given by

where *m* is the combined mass per square meter of sea ice and snow; $\mathit{u}=u\mathit{i}+v\mathit{j}$ is the horizontal velocity vector with components *u* and *v*; ** i**,

**, and**

*j***are unit vectors, respectively, aligned with the**

*k**x*,

*y*, and

*z*axes of the coordinate system;

*f*is the Coriolis parameter;

*τ*_{a}is the air stress;

*τ*_{w}is the water (or ocean) stress;

*τ*_{b}is a seabed stress which represents the effect of grounded pressure ridges; ∇⋅

**is the rheology term with horizontal stress components**

*σ**σ*

_{11}=

*σ*

_{xx},

*σ*

_{22}=

*σ*

_{yy}, and

*σ*

_{12}=

*σ*

_{xy};

*g*

_{e}is the Earth's gravitational acceleration; and ∇

*H*

_{0}is the sea surface tilt.

When using the B-grid discretization in CICE, there are different approaches for representing sea ice rheology and for solving the momentum equation: the viscous–plastic (VP; Hibler, 1979) rheology, which involves an implicit solution; the elastic–viscous–plastic framework (EVP; Hunke, 2001), which is based on the VP rheology but relies on an explicit method; the revised EVP approach (rEVP; Lemieux et al., 2012; Bouillon et al., 2013; Kimmritz et al., 2015), and the elastic–anisotropic–plastic model (EAP; Tsamados et al., 2013). In the current C-grid implementation, only the EVP and rEVP approaches are available. The EVP implementation is presented below.

Before describing the EVP equations for the internal stresses, we list a few modifications that were made recently to improve the flexibility of the VP (B-grid only) and (r)EVP (B-grid and C-grid) approaches. First, following König Beatty and Holland (2010), the yield curve can include tensile strength (see Lemieux et al., 2016, for details about the implementation in CICE). Tensile strength improves the simulation of landfast ice in regions of deep water (Lemieux et al., 2016). Second, the stresses are formulated in terms of viscosities, as introduced by Hibler (1979). Although only the elliptical yield curve is currently available, the formulation with viscosities offers more flexibility for defining other yield curves (e.g., Zhang and Rothrock, 2005). Finally, the current implementation includes the plastic potential approach of Ringeisen et al. (2021). Due to the normal flow rule, the standard VP rheology tends to simulate fracture angles that are too wide compared to observations (Ringeisen et al., 2019; Hutter et al., 2022). This problem can be remedied with the use of the plastic potential, which defines post-fracture deformations (or flow rule) independently from the yield curve. Following Ringeisen et al. (2021), the plastic potential is also defined by an elliptical curve.

Given these latest innovations, the EVP equations for the internal stresses are given by

where ${\mathit{\sigma}}_{\mathrm{1}}={\mathit{\sigma}}_{\mathrm{11}}+{\mathit{\sigma}}_{\mathrm{22}}$; ${\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{11}}-{\mathit{\sigma}}_{\mathrm{22}}$; *p* is the replacement pressure (defined below); *ζ* and $\mathit{\eta}={e}_{G}^{-\mathrm{2}}\mathit{\zeta}$ are, respectively, the bulk and shear viscosities; *e*_{G} is the plastic potential's ellipse ratio of major to minor axes; and *T*_{d} is a damping timescale for elastic waves (Hunke, 2001). *T*_{d} is defined as *T*_{d}=*E*_{0}Δ*t*, where $\mathrm{0}<{E}_{\mathrm{0}}<\mathrm{1}$ is a parameter and Δ*t* is the advective time step. ${D}_{\mathrm{d}}={\dot{\mathit{\u03f5}}}_{\mathrm{11}}+{\dot{\mathit{\u03f5}}}_{\mathrm{22}}$, ${D}_{\mathrm{t}}={\dot{\mathit{\u03f5}}}_{\mathrm{11}}-{\dot{\mathit{\u03f5}}}_{\mathrm{22}}$, and ${D}_{\mathrm{s}}=\mathrm{2}{\dot{\mathit{\u03f5}}}_{\mathrm{12}}$ are the divergence, the horizontal tension, and the shearing strain rate defined from the components of a symmetric strain rate tensor. To improve the damping of elastic waves, Eqs. (3) and (4) follow the formulation proposed by Bouillon et al. (2013) (see also Koldunov et al., 2019, for details).

The bulk viscosity, *ζ*, is given by

where *P* is the ice strength, *k*_{t} is a parameter between 0 and 1 that determines tensile strength (König Beatty and Holland, 2010) and Δ is a deformation (i.e., strain rate) associated with the elliptical yield curve and expressed as $\mathrm{\Delta}={\left[{D}_{d}^{\mathrm{2}}+\frac{{e}_{F}^{\mathrm{2}}}{{e}_{G}^{\mathrm{4}}}({D}_{t}^{\mathrm{2}}+{D}_{s}^{\mathrm{2}})\right]}^{\mathrm{1}/\mathrm{2}}$, with *e*_{F} being the elliptical yield curve axis ratio. When Δ tends toward zero, *ζ* tends toward infinity. To prevent this singularity, the denominator Δ in Eq. (5) is replaced by Δ^{*}. There are two approaches in the code to define Δ^{*}. By default, the capping approach of Hibler (1979) is used. In this case, ${\mathrm{\Delta}}^{*}=max[\mathrm{\Delta},{\mathrm{\Delta}}_{\mathrm{min}}]$, where Δ_{min} is a small deformation. A second approach with ${\mathrm{\Delta}}^{*}=(\mathrm{\Delta}+{\mathrm{\Delta}}_{\mathrm{min}})$ allows for a smoother formulation (Kreyscher et al., 2000). Finally, the replacement pressure, *p*, ensures that stresses are zero in the absence of deformations:

The ice strength can be calculated with the approach of Hibler (1979) (referred to as H79) or with the Rothrock (1975) parameterization (referred to as R75). With H79, *P* is given by

where *P*^{*} and *C*^{*} are parameters, *a* is the sea ice concentration, and $\stackrel{\mathrm{\u203e}}{h}$ is the mean thickness. Details about the more complicated R75 approach can be found in Rothrock (1975) and Lipscomb et al. (2007).

## 2.2 Transport equation

Roach et al. (2018) recently introduced a joint probability distribution, *f*(*h*,*r*), of sea ice thickness, *h*, and floe size, *r*, in CICE. For simplicity, the transport equation is introduced here by only considering the sea ice thickness distribution, *g*(*h*). Further information about horizontal transport in CICE can be found in Lipscomb and Hunke (2004).

The evolution of *g*(*h*) is given by

where *f*_{t} is the rate of thermodynamic ice growth and *ψ* is a mechanical redistribution term. CICE solves this equation using an operator-splitting approach; the change in *g*(*h*) from one time level to the next is computed in three steps, where only one term on the right hand side is non-zero at each step. The last two terms in Eq. (8) are handled by the column physics model in CICE called Icepack. Here, we only consider the change in *g*(*h*) due to the horizontal transport term $-\mathrm{\nabla}\cdot \left(g\mathit{u}\right)$.

The spatial discretization of the new C-grid implementation is based on finite differences using curvilinear coordinates on a fixed Eulerian mesh. It differs from the B-grid implementation in the way the stresses, strain rates, and other terms in Eqs. (1)–(4) are discretized. Our implementation mostly follows the work of Bouillon et al. (2009), Bouillon et al. (2013), and Kimmritz et al. (2016). In this section, we give a brief overview of the C-grid spatial and temporal discretizations of the momentum equation. Appendix A describes in detail the C-grid spatial discretization of the air stress (Sect. A1), the seabed stress (Sect. A2), and the rheology term, as well as the time-stepping of the internal stresses (Sect. A3); Sect. A4 presents the time-stepping of the momentum equation. For both B- and C-grid implementations, the momentum equation is first advanced in time (with the ice thickness distribution held fixed), which is followed by the transport equation using the newly calculated sea ice velocity field.

Figure 1 shows where variables are defined for the C-grid discretization. Scalar variables such as ice thickness and ice concentration are defined at the tracer point, *T*. Unlike the B-grid, where both velocity components are co-located at the corners (the *U* points), the C-grid *u* component is at the midpoint of the east (*E*) edge, while the *v* component is at the midpoint of the north (*N*) edge. In the derivations below, a variable such as *u*_{E} refers to the *u* component of velocity evaluated on the east edge (Fig. 1), and the same idea applies for variables with subscript *N*, *U*, or *T*, which are, respectively, evaluated at the north edge, the northeast corner, or the tracer point. The land–ocean mask is originally defined at the *T* point and referred to as *M*_{T}, with *M*_{T}=1 for ocean cells and *M*_{T}=0 for land cells. Other masks (*M*_{U}, *M*_{E}, and *M*_{N}) are defined based on *M*_{T}. For example, *M*_{U}=1 only if the four surrounding cells are ocean cells in the *M*_{T} mask.

For the momentum equation, the most complex part of the C-grid implementation is the spatial discretization of the rheology term (Sect. A3). This involves the calculations of strain rates at the *U* (Sect. A3.1) and *T* (Sect. A3.2) points; of *ζ*, *η*, and the replacement pressure at the *T* points (Sect. A3.3); and of *η* at the *U* points (Sect. A3.5). No-slip and no-outflow boundary conditions are applied at the land–ocean boundaries using ghost velocities (see Sect. A3.1 for details). Following Bouillon et al. (2013), ${D}_{\mathrm{s}}^{\mathrm{2}}$ at a *T* point is obtained from a spatial average of ${D}_{s}^{\mathrm{2}}$ from the four neighboring *U* points. This is done to enhance numerical stability. Following Kimmritz et al. (2016), the code includes two methods for calculating *η* at the *U* points; the default method averages *η*_{T} from the neighboring ocean cells while the other approach is based on an averaging of the ice strength. Section A3.4 and A3.6 describe the time-stepping of the stresses at the *T* and *U* points. The calculation of the *x* and *y* components of ∇⋅** σ** (Sect. A3.7) at the

*E*points (

*F*

_{1E}) and the

*N*points (

*F*

_{2N}) requires

*σ*

_{1T},

*σ*

_{2T}, and

*σ*

_{12U}. Also,

*σ*

_{12T}is calculated in order to diagnose normalized internal stresses (Lemieux and Dupont, 2020) at the

*T*points.

Considering only the horizontal transport term, we discretize Eq. (8) in terms of partial ice concentration. The transport equation for thickness category *n* is given by

where *a*_{n} is the partial ice concentration for category *n*. Snow volume, ice volume, snow enthalpy, and ice enthalpy also need to be transported using equations of the form

where *h*_{n} is the snow or ice thickness for category *n* and *q*_{n} is the enthalpy.

The incremental remapping scheme solves Eqs. (9)–(11) in a unified way. Given the velocity field ** u**, departure regions are computed for each grid cell. Then the quantities

*a*

_{n},

*a*

_{n}

*h*

_{n}, and

*a*

_{n}

*h*

_{n}

*q*

_{n}are integrated over each departure region so that volume and internal energy are transferred conservatively between cells.

Our initial implementation of remapping for the C-grid discretization mostly followed that in Lipscomb and Hunke (2004). The differences are related to the interpolation of C-grid velocity components to the *U* points, as described in Sect. 4.1 below.

## 4.1 Interpolation of C-grid velocity components and computation of departure regions

Departure regions in the remapping transport scheme are defined by approximating backward trajectories using corner velocities (*U* points). In our first implementation of remapping for the C-grid, C-grid velocities were interpolated to the corners, and remapping was used in the same way as for the B-grid. The velocity components *u*_{E} and *v*_{N} are interpolated to the *U* points as

where *A*_{E} and *A*_{N} are grid cell areas evaluated at the *E* and *N* points, respectively, and the multiplication by *M*_{U}(*i*,*j*) ensures that the no-slip and no-outflow boundary conditions (BCs) are respected.

To improve the accuracy of the estimated departure regions, the midpoints of the backward trajectories are computed first. Then, velocity components are bilinearly interpolated to these midpoints. Finally, these interpolated velocities are used to calculate the departure points defining the departure regions (Lipscomb and Hunke, 2004).

Panel (a) in Fig. 2 shows an example of a departure region on the north edge of cell (*i*,*j*). The departure region is a quadrilateral defined by the left (cl) and right (cr) corner points and the left (dl) and right (dr) departure points.

## 4.2 Weaknesses of the initial approach for remapping using C-grid velocity components

We identified two notable weaknesses with this initial C-grid discretization and remapping implementation. First, a C-grid discretization offers the possibility of representing transport in one-grid-cell-wide channels, but our initial implementation did not do so. Because of the no-slip and no-outflow BCs, the *U* velocities are zero on both sides of such a channel, in which case the departure regions have zero area. Second, we found a numerical problem: in some idealized tests, we observed a checkerboard pattern in fields such as sea ice concentration. Panel (a) in Fig. 3 shows an example of this pattern, which indicates the presence of a spurious numerical mode (or modes). This numerical noise is not present when using the upwind scheme (not shown).

Appendix B presents a modal analysis of a simplified set of perturbed equations (momentum and transport). We show that a stationary wave explains the formation of the checkerboard pattern. This stationary wave is a consequence of the spatial averaging used to obtain *u*_{U} and *v*_{U} (Eqs. 12 and 13) for the remapping scheme.

## 5.1 Description of the method

To eliminate the checkerboard pattern, we introduce a novel approach, which we refer to as the edge flux adjustment (EFA) method. With this method, the remapping scheme still calculates departure regions using the *U*-point velocities. However, the C-grid velocity components are then used to adjust the edge fluxes (or departure regions.) Following the example in Fig. 2, the departure region is adjusted so that the total flux area *A*_{tot} is equal to *v*_{N}Δ*x*_{N}Δ*t*, where Δ*t* is the advective time step and *v*_{N}<0 in the example shown in Fig. 2. (Note that *A*_{tot} has the same sign as *v*_{n}). This is done by shifting the midpoint, dm, of the line segment connecting the departure points, dl and dr. Geometrically speaking, this creates an additional departure triangle (the orange triangle in Fig. 2b).

The code computes the coordinates of the shifted dm in a non-dimensional coordinate system. In this coordinate system, the points cl and cr have the coordinates $(-\mathrm{0.5},\mathrm{0})$ and (0.5,0). The non-dimensional flux area, denoted as ${\stackrel{\mathrm{\u0303}}{A}}_{\mathrm{tot}}$, is equal to ${A}_{\mathrm{tot}}/{A}_{N}$, where *A*_{N} is the grid cell area evaluated at the *N* point. Given the non-dimensional coordinates (*x*_{dl},*y*_{dl}) and (*x*_{dr},*y*_{dr}) for the departure points, the non-dimensional coordinates $({x}_{\mathrm{dm}}^{*},{y}_{\mathrm{dm}}^{*})$ for the shifted departure midpoint are obtained as follows:

where ${x}_{\mathrm{dm}}=({x}_{\mathrm{dl}}+{x}_{\mathrm{dr}})/\mathrm{2}$ and ${y}_{\mathrm{dm}}=({y}_{\mathrm{dl}}+{y}_{\mathrm{dr}})/\mathrm{2}$ are the initial departure midpoint coordinates and *α* is given by

The adjusted area flux can finally be computed using the points $\mathrm{cl}=({x}_{\mathrm{cl}},{y}_{\mathrm{cl}})$, $\mathrm{dl}=({x}_{\mathrm{dl}},{y}_{\mathrm{dl}})$, $\mathrm{cr}=({x}_{\mathrm{cr}},{y}_{\mathrm{cr}})$, $\mathrm{dr}=({x}_{\mathrm{dr}},{y}_{\mathrm{dr}})$, and ${\mathrm{dm}}^{*}=({x}_{\mathrm{dm}}^{*},{y}_{\mathrm{dm}}^{*})$.

The initial departure region shown in Fig. 2a is confined to the central region (i.e., cell (*i*,*j*) and/or cell $(i,j+\mathrm{1}))$. When part of the departure region (i.e., a triangle) is located, for example, in a corner cell (e.g., the northwestern cell $(i-\mathrm{1},j+\mathrm{1})$, not shown), the area of this corner triangle is subtracted from *A*_{tot} before adjusting the central portion of the departure region (the part lying in cell $(i,j+\mathrm{1})$). In this case, we identify the point (0,*y*_{i}) where the segment joining dl and dr intersects the left edge of cell $(i,j+\mathrm{1})$. The departure point dl is reset to (0,*y*_{i}) before finding the shifted midpoint, dm^{*}. Thus, the new triangle with vertices dl, dr, dm^{*} is always located within the central region.

The calculation of $({x}_{\mathrm{dm}}^{*},{y}_{\mathrm{dm}}^{*})$ as described above is the most common case, with both initial departure points on the same side of the edge (i.e., *y*_{dr}*y*_{dl}≥0). We give another less common example in which *y*_{dr} and *y*_{dl} have opposite signs. Then, there are two departure triangles: one on the left, with vertices cl, dl, ip, and one on the right, with vertices cr, dr, ip, where ip denotes the point (*x*_{i},0) at which the departure segment intersects the *x* axis.

In the case shown (Fig. 4), the non-dimensional coordinate *x*_{i} of the intersection point is greater than zero. There is a similar case not discussed here, with *x*_{i}<0.

The strategy is to fix the right-hand triangle (the “lone triangle”) while modifying the left-hand triangle. Given the required total area flux, *A*_{tot}=*v*_{N}Δ*x*_{N}Δ*t*, and the area *A*_{r△} of the right-hand triangle (in yellow), the EFA method first calculates the remaining area flux, ${A}_{l}={A}_{\mathrm{tot}}-{A}_{\mathrm{r}\mathrm{\u25b3}}$. The middle departure point is reset to $(\frac{{x}_{\mathrm{dl}}+{x}_{i}}{\mathrm{2}},\frac{{y}_{\mathrm{dl}}}{\mathrm{2}})$, and then dm is shifted so that the non-dimensional area of the green and orange triangles is equal to ${\stackrel{\mathrm{\u0303}}{A}}_{l}={A}_{l}/{A}_{N}$. The shifted middle departure point has non-dimensional coordinates

where in this case *α* is given by

When a part of the departure region is a corner triangle lying outside cells (*i*,*j*) and $(i,j+\mathrm{1})$, it is treated as described above for the most common case.

The EFA method ensures that the divergence (associated with edge fluxes) implied by the remapping is consistent with the divergence calculated by the dynamical solver (i.e., the EVP solver; see Eq. A24). Figure 3b shows that the EFA method prevents the formation of the checkerboard pattern. This pattern originates from an interaction between the solution of the momentum equation and the standard remapping scheme. To support this conclusion, we conducted the following experiment: velocity fields obtained with the C-grid discretization and the use of the EFA method in the remapping scheme were first stored. In a second simulation, both the dynamics (i.e., EVP) and EFA method are turned off and the stored velocity fields are used for transport in the remapping scheme. In this case, no checkerboard pattern develops, and the concentration field is very similar to the one shown in Fig. 3b.

The EFA code was added to the standard CICE remapping algorithm described in Lipscomb and Hunke (2004). Although the remapping algorithm for the B-grid is very robust, we found that some approximations in the computation of departure regions can be problematic with the new EFA method for the C-grid discretization. Long-term C-grid simulations with the EFA method exhibited rare failures on non-uniform grids. The remapping code in CICE includes many consistency checks to ensure that the solution is physically sound – for example, that transport does not lead to negative area or mass. The code failed on rare occasions with negative area and mass values close to land or the ice edge. These negative values were a result of approximations in the area of the departure region. Appendix C describes the changes in the remapping algorithm that were required to fix this robustness issue. Note that many of these modifications slightly alter the departure areas for both the B-grid and the C-grid (with the EFA method) discretizations. With these changes to the remapping, the algorithm is now robust for both B-grid and C-grid simulations on non-uniform grids.

## 5.2 Transport through one-grid-cell-wide channels

Interestingly, the EFA method also remedies the other weakness of our initial implementation: the absence of transport in one-grid-cell-wide channels. Although the initial departure regions have zero area (because the *U*-point velocities are zero along the channel edges), the departure regions are adjusted based on C-grid velocity components (e.g., the *v*_{N} component would be non-zero for a north–south channel). In this case, the departure region is simply defined by (*x*_{cl},*y*_{cl}), (*x*_{cr},*y*_{cr}), and $({x}_{\mathrm{dm}}^{*},{y}_{\mathrm{dm}}^{*})$, with ${x}_{\mathrm{dm}}^{*}=\mathrm{0}$.

A minor code modification was required nevertheless. Before this modification, the edge fluxes were calculated only when at least one of the two departure points was displaced from its corner. Given the non-displaced departure points for the north edge (i.e., ${u}_{U}(i-\mathrm{1},j)={v}_{U}(i-\mathrm{1},j)=\mathrm{0}$ and ${u}_{U}(i,j)={v}_{U}(i,j)=\mathrm{0}$), edge fluxes are now computed whenever $\left|{v}_{N}\mathrm{\Delta}{x}_{N}\mathrm{\Delta}t\right|>\mathrm{0}$. Similarly for the east edge (i.e., ${u}_{U}(i,j)={v}_{U}(i,j)=\mathrm{0}$ and ${u}_{U}(i,j-\mathrm{1})={v}_{U}(i,j-\mathrm{1})=\mathrm{0}$), edge fluxes are now computed whenever $\left|{u}_{E}\mathrm{\Delta}{y}_{E}\mathrm{\Delta}t\right|>\mathrm{0}$. Given that the remapping with the EFA method uses velocity components at the *U* point and at the *E* and *N* points and that it allows for transport in one-grid-cell-wide channels, it should be viewed as a hybrid B-and-C-grid approach.

To test the new capability for transport, we implemented an idealized configuration with a long one-grid-cell-wide channel. The initial ice conditions are *a*=0.5 and $\stackrel{\mathrm{\u203e}}{h}=\mathrm{1}$ m over a length of five grid cells (80 km), with *a*=0 and $\stackrel{\mathrm{\u203e}}{h}=\mathrm{0}$ m elsewhere. The wind blows from the west, and the fields are analyzed after 30 d of simulation. As expected, there is no transport with the B-grid discretization (not shown). Using the C-grid discretization, both the upwind method and the remapping with the EFA method lead to transport in the channel (Fig. 5). As shown by Lipscomb and Hunke (2004) for more complex configurations, remapping is much less diffusive than upwind.

We validated the new C-grid implementation with the EFA method using several approaches, including (1) symmetry tests, (2) thorough comparison of C-grid simulations with the default B-grid simulations, (3) comparison of C-grid simulations with analytical solutions, and (4) diagnostics of simulated stress states. This section gives an overview of the different tests. For all the simulations presented here, we used version 6.5.0 of CICE, which includes our C-grid modifications. Unless otherwise specified, we used default values for most physical and numerical parameters. In the experiments described below, we mostly modified the strength parameterization, *P*^{*}, Δ_{min}, the number of EVP subcycles (*n*_{evp}), and *E*_{0}, in order to broaden the variety of tests (e.g., ice strength parameterization) and in some cases to improve the numerical convergence of the EVP method (*n*_{evp} and *E*_{0}). See Table 1 for experiment details.

We conducted many idealized tests to verify the symmetry of simulated fields. Figure 6 shows an example. The thickness fields after 14 d are perfectly (bit-for-bit) symmetrical for the four oblique (i.e., northeast, southeast, southwest, and northwest) wind directions. Similar tests with the wind blowing from the west, east, north, and south also give symmetrical results, but with small differences (the maximum difference is 4 × 10^{−4} m). Changing the capping method to the smoother formulation (i.e., capping_method = 'sum') leads to bit-for-bit symmetrical results.

We also evaluated more realistic simulations with B-grid runs as references. We compared C-grid runs on 1 and 3° global grids to B-grid runs initialized and forced by the same fields. The Japanese 55-year Reanalysis (JRA-55; Kobayashi et al., 2015) is used for the atmospheric forcing fields, while the ocean forcing was derived from a Community Earth System Model (CESM) simulation (Kay et al., 2015). Panels (a) and (b), respectively, in Fig. 7 show the total simulated sea ice volume for the Northern and Southern hemispheres for a 1° B-grid simulation with remapping transport (reference, orange), a 1° B-grid simulation with upwind transport (blue), and a 1° C-grid simulation with remapping transport (dashed violet). Only the C-grid simulation uses the EFA method. Compared to the reference B-grid simulation, changing the grid discretization has a smaller impact on the total volume than changing the advection scheme. This is particularly evident in the Southern Hemisphere.

This conclusion is further supported by spatial maps of sea ice thickness. The monthly mean ice thickness in December 2009 for a 1° C-grid simulation with remapping is qualitatively correct, with the thickest ice found north of Greenland and in the Canadian Arctic Archipelago (Fig. 8). When compared to the reference simulation (B-grid with remapping), the ice is thinner in the regions of thick ice (Fig. 9a), although these differences are generally less pronounced than those for a B-grid with upwind compared to the reference (Fig. 9b). The same is true for the Southern Hemisphere (not shown).

The checkerboard pattern was not visible in our initial 1 and 3° C-grid simulations without EFA. Instead, the error manifested as much thicker ice in convergent regions. We speculate this was due to a feedback between excessive ridging associated with the checkerboard pattern in velocities (and hence divergence and convergence) and ice growth in open water formed through the ridging process.

The most complex part of the discretization of the momentum equation is the rheology term. A crucial test is to verify whether the simulated internal stresses are inside (viscous) or on (plastic) the elliptical yield curve. To improve the numerical convergence of the EVP solver, the number of subcycling iterations, *n*_{evp}, was increased from 240 (default) to 1200 and the elastic damping parameter, *E*_{0}, was reduced from 0.36 (default) to 0.09 (Exp5 in Table 1). Following the approach of Lemieux and Dupont (2020), we plotted the states of stress in stress-invariant coordinates of a snapshot after 5 d of a 1° C-grid simulation. Figure 10 confirms that the solution is viscous–plastic.

Analytical solutions are useful tools for verifying a numerical implementation. We derived novel analytical solutions for a one-grid-cell-wide channel with cyclic boundary conditions; see Appendix D for the derivation. As sea ice conditions are assumed constant in space and time, these solutions cannot be used to verify the simulated transport. Nevertheless, they provide steady-state analytical velocity values. Following the assumptions of this analytical solution with *a*=0.8 and $\stackrel{\mathrm{\u203e}}{h}=\mathrm{0.8}$ m for an east–west channel, the ice should be in the plastic (viscous) regime for wind speeds higher (lower) than ${u}_{\mathrm{a}*}=\mathrm{3.176}$ m s^{−1}. For *u*_{a}=4 m s^{−1} (i.e., ${u}_{\mathrm{a}}>{u}_{\mathrm{a}*}$), the steady-state analytical solution, *u*_{E}, calculated independently using a Python code is 0.04095 m s^{−1}. The steady-state solution, *u*_{E}, obtained with the new C-grid implementation matches the analytical solution up to the 12th digit. Similarly, for *u*_{a}=1.5 m s^{−1} (i.e., ${u}_{\mathrm{a}}<{u}_{\mathrm{a}*}$), the steady-state analytical solution *u*_{E} is much smaller and is equal to 7.135 × 10^{−6} m s^{−1}. The new C-grid steady-state solution, *u*_{E}, matches the analytical solution up to the ninth digit (a difference of less than 10^{−16} m s^{−1}). The same results are found when using a north–south channel. This experiment is referred to as Exp6 in Table 1.

We have designed and implemented a C-grid version of the CICE sea ice model. The C-grid spatial discretization is based on a finite-difference approach and follows the work of Bouillon et al. (2009, 2013) and Kimmritz et al. (2016). This article describes the finite-difference spatial discretization of the momentum equation, the implementation of no-slip and no-outflow boundary conditions, and the use of the remapping transport scheme with C-grid velocities.

The most notable contribution of this work is a novel method for the remapping referred to as the edge flux adjustment (EFA) method. Preliminary results from idealized experiments showed that the new C-grid discretization for the momentum equation and the use of the standard remapping transport scheme could produce checkerboard patterns in fields such as ice concentration. This numerical noise is not present when using the upwind transport scheme. A modal analysis of a simplified set of perturbed equations (i.e., momentum and transport with spatial discretization) shows that a stationary wave is responsible for the checkerboard pattern. This stationary wave results from the interpolation of C-grid velocity components to the *U* points for use with remapping, which is fundamentally a B-grid scheme (Dukowicz and Baumgardner, 2000; Lipscomb and Hunke, 2004). This interpolation can be viewed as a spatial averaging. Many authors (e.g., Batteen and Han, 1981) have demonstrated that spatial averaging can lead to checkerboard patterns when solving the shallow-water equations, which are similar to our simplified set of equations. The checkerboard pattern is eliminated by the EFA method.

The EFA method uses C-grid velocity components at their natural locations to modify the departure regions calculated by the remapping so that the implied divergence in the remapping is consistent with the divergence calculated by the dynamical solver. We also introduced some modifications to the calculation of the area of departure regions to increase the robustness of remapping with the EFA method on non-uniform grids.

A C-grid discretization offers the possibility of representing transport in one-grid-cell-wide channels. Because of the no-slip and no-outflow boundary conditions, the *U*-point velocities at the channel edges are zero, and there is therefore no transport when using the standard remapping. The EFA method, however, allows for transport in these channels by creating departure regions with non-zero area based on the C-grid velocities. As such, the remapping with the EFA method should be seen as a hybrid B-and-C-grid transport scheme. Moreover, remapping with the EFA method is much less diffusive than the upwind method for idealized channel tests.

For VP sea ice models, there are few existing analytical solutions due to the complexity of the rheology. We derived novel analytical solutions for one-grid-cell-wide channels and showed that with several simplifications (uniform ice conditions, constant wind, cyclic boundary conditions, and transport turned off), the sea ice velocity can be obtained analytically for both plastic and viscous regimes. The steady-state values simulated by CICE match the analytical ones.

We also conducted multiyear 1° global simulations to compare the C-grid solution (using the EFA method) to the reference B-grid solution (with standard remapping). We ran an additional B-grid simulation with upwind transport. Compared to the reference B-grid run, the C-grid discretization has a smaller impact on the total volume (and spatial differences) than changing the advection scheme from remapping to upwind, especially in the Southern Hemisphere.

Ongoing work within CICE Consortium modeling centers includes coupling this new CICE C-grid implementation to ocean models, such as the Modular Ocean Model (MOM6) and the Nucleus for European Modelling of the Ocean (NEMO), and to atmospheric models, such as the Global Environmental Multiscale Model (GEM).

The spatial discretization is presented below in the same order as in the code.

## A1 Air stress at the *E* and *N* points

For both B-grid and C-grid implementations, the air stress is calculated at the *T* point and then interpolated to the required locations. For the C-grid, *τ*_{ax} at the *E* point and *τ*_{ay} at the *N* point are weighted averages of the values at the neighboring *T* points and are given by

where *A*_{E}, *A*_{N}, and *A*_{T} are cell areas evaluated at the *E*, *N*, and *T* points.

## A2 Seabed stress at the *E* and *N* points

The seabed stress components are ${\mathit{\tau}}_{bxE}=-{C}_{bE}{u}_{E}$ and ${\mathit{\tau}}_{byN}=-{C}_{bN}{v}_{N}$, where the *C*_{b} coefficients are calculated as in Lemieux et al. (2016) or following the probabilistic approach of Dupont et al. (2022). For both approaches, *C*_{bE} and *C*_{bN} are written as

where *T*_{bE} and *T*_{bN} are factors that characterize the maximum possible seabed stress and *u*_{0} is a small velocity parameter that ensures a smooth transition between the static and dynamic regimes of the seabed stress. Velocity *v*_{E} is obtained by interpolating *v*_{N} to the *E* point, and *u*_{N} is found similarly by interpolating *u*_{E} to the *N* point. These interpolations are currently different than in the C-grid implementation of Lemieux et al. (2015) but will be made consistent in a subsequent version of CICE.

As opposed to the denominators in Eqs. (A3) and (A4), factors *T*_{bE} and *T*_{bN} do not vary during the EVP subcycling. They are therefore calculated before the subcycling. Following the approach of Lemieux et al. (2016), *T*_{bE} and *T*_{bN} are

where ${\stackrel{\mathrm{\u203e}}{h}}_{E}=max\left[{\stackrel{\mathrm{\u203e}}{h}}_{T}\right(i,j),{\stackrel{\mathrm{\u203e}}{h}}_{T}(i+\mathrm{1},j\left)\right]$; ${a}_{E}=max\left[{a}_{T}\right(i,j),{a}_{T}(i+\mathrm{1},j\left)\right]$; ${h}_{cE}={a}_{E}{h}_{\mathrm{w}E}/{k}_{\mathrm{1}}$, ${h}_{\mathrm{w}E}=min\left[{h}_{\mathrm{w}T}\right(i,j),{h}_{\mathrm{w}T}(i+\mathrm{1},j\left)\right]$; ${\stackrel{\mathrm{\u203e}}{h}}_{N}=max\left[{\stackrel{\mathrm{\u203e}}{h}}_{T}\right(i,j),{\stackrel{\mathrm{\u203e}}{h}}_{T}(i,j+\mathrm{1}\left)\right]$; ${a}_{N}=max\left[{a}_{T}\right(i,j),{a}_{T}(i,j+\mathrm{1}\left)\right]$; and ${h}_{cN}={a}_{N}{h}_{\mathrm{w}N}/{k}_{\mathrm{1}}$, ${h}_{\mathrm{w}N}=min\left[{h}_{\mathrm{w}T}\right(i,j),{h}_{\mathrm{w}T}(i,j+\mathrm{1}\left)\right]$. ${\stackrel{\mathrm{\u203e}}{h}}_{E}$ (${\stackrel{\mathrm{\u203e}}{h}}_{N}$), *h*_{cE} (*h*_{cN}), and *h*_{wE} (*h*_{wN}) are the mean ice thickness (or ice volume), the critical thickness, and the water depth at the *E* (*N*) point calculated using values at the *T* point. *k*_{1}, *k*_{2}, and *α*_{b} are three parameters of the seabed stress parameterization (Lemieux et al., 2016). The *T*_{b} factors are set to zero when the water depth is greater than 30 m.

When the seabed stress is computed based on the probabilistic approach, the calculation of the *T*_{b} factors is more complicated than with the method of Lemieux et al. (2016). Details can be found in Dupont et al. (2022). With the probabilistic approach, the *T*_{b} factor is first calculated at the *T* points, and *T*_{bE} and *T*_{bN} are then given by

## A3 Discretization of rheology

As opposed to the variational method used for the B-grid (Hunke and Dukowicz, 1997, 2002), our C-grid spatial discretization is based on finite differences. With this approach, the discretization of ∇⋅** σ** requires the calculation of

*σ*

_{1}and

*σ*

_{2}at the

*T*points and of

*σ*

_{12}at the

*U*points. The stresses are calculated in three steps: the computation of strain rates, the computation of viscosities and replacement pressure, and finally the time-stepping of the stresses from subcycle

*k*to

*k*+1. The subsections below explain how this is done for the

*T*and

*U*points. The sequence of computations follows that in the code.

The spatial discretization requires the computation of strain rates and components of the rheology term in curvilinear coordinates. The strain rates are given by

where *D*_{d} is the divergence, *D*_{t} is the tension, *D*_{s} is the shear strain rate, *ξ*_{1} and *ξ*_{2} are the non-dimensional coordinates, and *h*_{1} and *h*_{2} are scale factors, referred to as Δ*x* and Δ*y*.

In curvilinear coordinates, the *x* and *y* components of the divergence of the stress tensor are, respectively, as follows:

which are Eqs. (20) and (21) in Hunke and Dukowicz (2002).

### A3.1 Strain rates at the *U* points

To solve the momentum equation, shear stress, *σ*_{12}, is needed at the *U* points, including at land–ocean boundaries. This implies that strain rates and shear viscosities must be computed at these locations. In the code, strain rates at the *U* points are first calculated. The reason for doing this is to follow what is suggested in Bouillon et al. (2013); to enhance numerical stability, ${D}_{\mathrm{s}T}^{\mathrm{2}}$ in Δ_{T} is a weighted average of the ${D}_{\mathrm{s}U}^{\mathrm{2}}$ around it.

As described in Sect. A3.5, there are two methods for calculating *η* at the *U* points. Following Kimmritz et al. (2016), we refer to these methods as C1 and C2. The C1 method requires only *D*_{sU}, while the C2 method requires *D*_{dU}, *D*_{tU}, and *D*_{sU}. C1 is the default method, but for completeness, we explain here how *D*_{dU}, *D*_{tU}, and *D*_{sU} are computed.

To ease the implementation of the boundary conditions (BCs) in the code, strain rates at the *U* points are calculated differently than at the *T* points. To do so, the derivatives in Eqs. (A9), (A10), and (A11) are expanded. First, for the divergence, Eq. (A9), we expand the derivatives and write

The discretized form of divergence at the *U* point (*i*,*j*) is therefore

where ${u}_{N}^{*}(i+\mathrm{1},j)$, ${u}_{N}^{*}(i,j)$, ${v}_{E}^{*}(i,j+\mathrm{1})$, and ${v}_{E}^{*}(i,j)$ are modified versions of ${u}_{N}(i+\mathrm{1},j)$, *u*_{N}(*i*,*j*), ${v}_{E}(i,j+\mathrm{1})$, and *v*_{E}(*i*,*j*) to take the BCs into account. This is explained below. The velocity components interpolated to the *U* points follow Eqs. (12) and (13), while (the unmodified) *v*_{E} and *u*_{N} are calculated as

where ${A}_{N\mathrm{tot}}={\sum}_{k=\mathrm{0}}^{\mathrm{1}}{\sum}_{l=-\mathrm{1}}^{\mathrm{0}}{A}_{N}(i+k,j+l)$ and ${A}_{E\mathrm{tot}}={\sum}_{k=-\mathrm{1}}^{\mathrm{0}}{\sum}_{l=\mathrm{0}}^{\mathrm{1}}{A}_{E}(i+k,j+l)$. For example, to clarify the notation, ${A}_{N\mathrm{tot}}={\sum}_{k=\mathrm{0}}^{\mathrm{1}}{\sum}_{l=-\mathrm{1}}^{\mathrm{0}}{A}_{N}(i+k,j+l)={A}_{N}(i,j-\mathrm{1})+{A}_{N}(i,j)+{A}_{N}(i+\mathrm{1},j-\mathrm{1})+{A}_{N}(i+\mathrm{1},j)$.

Similarly, to calculate *D*_{t} at the *U* points, we expand the derivatives in Eq. (A10) and write the tension as

The discretized form of this equation is

Finally, for the shear strain rate, we write Eq. (A11) as

with the following discretized form:

Away from the land–ocean boundaries, ${u}_{E}^{*}(i,j)={u}_{E}(i,j)$, ${u}_{E}^{*}(i,j+\mathrm{1})={u}_{E}(i,j+\mathrm{1})$, ${v}_{N}^{*}(i,j)={v}_{N}(i,j)$, etc. However, at ocean–land boundaries, no-slip and no-outflow BCs are implemented by setting ${u}_{U}(i,j)={v}_{U}(i,j)=\mathrm{0}$ and using ghost values for the other terms. As an example, consider *v*_{N}(*i*,*j*) and ${v}_{N}(i+\mathrm{1},j)$ when the *T* cells at $(i+\mathrm{1},j)$ and $(i+\mathrm{1},j+\mathrm{1})$ are land cells. We need *v*_{N}(*i*,*j*) and ${v}_{N}(i+\mathrm{1},j)$ to calculate the $\partial v/\partial {\mathit{\xi}}_{\mathrm{1}}$ term in *D*_{sU}. As *v*_{N}(*i*,*j*) is in the ocean, ${v}_{N}^{*}(i,j)={v}_{N}(i,j)$. However, as ${v}_{N}(i+\mathrm{1},j)$ is on land, it is not defined and must be formulated using the BCs. We assume that *v*_{N} varies linearly at the ocean–land boundary. We therefore write ${v}_{N}=mx+b$, where *m* is the slope and *b* is the value of *v*_{N} at *x*=0, which is defined at the ocean–land boundary. Using the no-outflow condition implies that *b*=0. Given that ${h}_{\mathrm{1}N}(i,j)/\mathrm{2}$ (the distance between the ocean–land boundary and the *N*(*i*,*j*) point) and ${h}_{\mathrm{1}N}(i+\mathrm{1},j)/\mathrm{2}$ (the distance between the ocean–land boundary and the $N(i+\mathrm{1},j)$ point), it is easy to show that

where in the case of a uniform Cartesian grid, ${v}_{N}^{*}(i+\mathrm{1},j)$ is simply *v*_{N}(*i*,*j*) multiplied by −1. To take into account all the possible cases, the mask at the *N* point (*M*_{N}) is used for the final formulation of ${v}_{N}^{*}(i+\mathrm{1},j)$:

which reduces to ${v}_{N}^{*}(i+\mathrm{1},j)={v}_{N}(i+\mathrm{1},j)$ away from the ocean–land boundary (i.e., all four *T* cells are ocean cells).

### A3.2 Strain rates at the *T* points

Using Eq. (A9), a finite-difference approximation of the divergence at the *T* point is given by

Similarly, using Eq. (A10), the tension at the *T* point is given by

Following Bouillon et al. (2013), ${D}_{\mathrm{s}T}^{\mathrm{2}}$ is obtained as a weighted average of the neighboring ${D}_{\mathrm{s}U}^{\mathrm{2}}$:

where *A*_{U}(*i*,*j*) is the cell area evaluated at the *U* point and ${A}_{U\mathrm{tot}}={\sum}_{k=-\mathrm{1}}^{\mathrm{0}}{\sum}_{l=-\mathrm{1}}^{\mathrm{0}}{A}_{U}(i+k,j+l)$. At the *T* point, the strain rate, Δ_{T}, for the viscosities is then calculated as

### A3.3 Viscosities and replacement pressure at the *T* points

Using Δ_{T}(*i*,*j*) as calculated in Eq. (A27), *ζ*_{T}(*i*,*j*), with the capping approach of Hibler (1979), is obtained as follows:

where the ice strength, *P*_{T}, is also evaluated at the *T* point. Similarly, the replacement pressure at the *T* point is

If *ζ*_{T} and *p*_{T} are regularized with the smoother approach as in Kreyscher et al. (2000), the denominator $max\left[{\mathrm{\Delta}}_{T}\right(i,j),{\mathrm{\Delta}}_{\mathrm{min}}]$ in Eqs. (A28) and (A29) is replaced by $\left({\mathrm{\Delta}}_{T}\right(i,j)+{\mathrm{\Delta}}_{\mathrm{min}})$. The approach of Hibler (1979) can be used by setting capping_method = 'max' in the namelist, while the smoother formulation is used by setting capping_method = 'sum'. Finally, the shear viscosity at the *T* point is simply ${\mathit{\eta}}_{T}(i,j)={e}_{G}^{-\mathrm{2}}{\mathit{\zeta}}_{T}(i,j)$.

### A3.4 Time-stepping of the stresses at the *T* points

For our C-grid implementation, only *σ*_{1} and *σ*_{2} are required at the *T* points for time-stepping the velocity components using the momentum equation. Nevertheless, *σ*_{12} is also computed at the *T* points in order to calculate normalized stresses (Lemieux and Dupont, 2020) as diagnostics. Following Eqs. (2)–(4), the stresses at the *T* points are time-stepped from subcycle *k* to *k*+1 as follows:

where Δ*t*_{e} is the subcycling time step.

It is straightforward to solve the equations above for ${\mathit{\sigma}}_{\mathrm{1}T}^{k+\mathrm{1}}(i,j)$, ${\mathit{\sigma}}_{\mathrm{2}T}^{k+\mathrm{1}}(i,j)$, and ${\mathit{\sigma}}_{\mathrm{12}T}^{k+\mathrm{1}}(i,j)$. Note that *ζ*_{T}, *η*_{T}, *p*_{T}, and strain rates in the equations above are calculated with a velocity field at iteration *k*.

### A3.5 Viscosities at the *U* points

With our C-grid implementation, only the shear viscosity, *η*, is needed at the *U* points. Two methods in the code can be used to calculate *η*_{U}. The default method (visc_method = 'avg_zeta' in the namelist) is a weighted spatial average of the values at the *T* points. This is the C1 method of Kimmritz et al. (2016) and the same method used in Bouillon et al. (2013). With the C1 method, *η*_{U} is obtained from a weighted average of the *η*_{T} values in ocean cells around the *U* points. This can be concisely written as

where ${A}_{T\mathrm{tot}}={\sum}_{k=\mathrm{0}}^{\mathrm{1}}{\sum}_{l=\mathrm{0}}^{\mathrm{1}}{A}_{T}(i+k,j+l){M}_{T}(i+k,j+l)$.

The second method (visc_method = 'avg_strength' in the namelist) relies on a weighted spatial average of the ice strength values at the surrounding ocean *T* points. This is the C2 method of Kimmritz et al. (2016) and also the method used in Bouillon et al. (2009). The ice strength at the *U* point is given by

where *A*_{Ttot} is the same as for Eq. (A33) above.

Given that Δ_{U}(*i*,*j*) = ${\left[{D}_{\mathrm{d}U}^{\mathrm{2}}(i,j)+{e}_{F}^{\mathrm{2}}{e}_{G}^{-\mathrm{4}}\left({D}_{\mathrm{t}U}^{\mathrm{2}}\right(i,j)+{D}_{\mathrm{s}U}^{\mathrm{2}}(i,j\left)\right)\right]}^{\mathrm{1}/\mathrm{2}}$, the shear viscosity at the *U* point with capping_method = 'max' is given by

With capping_method = 'sum', it is given by

### A3.6 Time-stepping of the stresses at the *U* points

Using *η*_{U} and *D*_{sU}, the shear stress at the *U* point is advanced in time from subcycle *k* to subcycle *k*+1 according to the following equation:

which can easily be solved for ${\mathit{\sigma}}_{\mathrm{12}U}^{k+\mathrm{1}}(i,j)$. Note that *η*_{U} and *D*_{sU} in the equation above are calculated with a velocity field at iteration *k*.

### A3.7 Divergence of the stress tensor

Once the stresses at the *T* and *U* points have been advanced in time from *k* to *k*+1, the components of the rheology term can be calculated. Equations (A12) and (A13) introduced earlier can be rewritten as

Using finite differences, the discretized formulation of *F*_{1} at the *E* points is

while the discretized formulation of *F*_{2} at the *N* points is

## A4 Time-stepping of the momentum equation

When using the B-grid discretization, the sea ice momentum equation in CICE can be solved either explicitly with the EVP (or revised EVP or EAP) approach or implicitly with a Picard solver (similar to the one described in Lemieux et al., 2008). For now, only the EVP and revised EVP approaches are implemented for the C-grid discretization.

As this subsection describes the time-stepping, the grid indices (*i*,*j*) are omitted to simplify the description. Hence, *u*_{E}(*i*,*j*) and *v*_{N}(*i*,*j*) are here referred to as *u*_{E} and *v*_{N}. Neglecting the advection of momentum and introducing the EVP time-stepping, the momentum equations for the *u*_{E} and *v*_{N} components are

where the interpolated quantities, *v*_{E} and *u*_{N}, are calculated using Eqs. (A16) and (A17). The terms *m*_{E} and *m*_{N} are

All the terms in Eq. (A42) are evaluated at the *E* points, while all the terms in Eq. (A43) are evaluated at the *N* points. The seabed stress components are ${\mathit{\tau}}_{bxE}=-{C}_{bE}{u}_{E}^{k+\mathrm{1}}$ and ${\mathit{\tau}}_{byN}=-{C}_{bN}{v}_{N}^{k+\mathrm{1}}$, where the *C*_{b} coefficients are calculated as in Lemieux et al. (2016) or following the probabilistic approach of Dupont et al. (2022). Decomposing the water stress term, Eqs. (A42) and (A43) can be written as

with

where *ρ*_{w} is the water density, *C*_{dw} is the ocean drag coefficient, *θ*_{w} is the turning angle, sin *θ*_{w} terms have a hemisphere-dependent sign, *u*_{w} and *v*_{w} are the near-surface water velocity components (evaluated at either the *E* or *N* point), and *a*_{E} and *a*_{N} (the concentration at the *E* and *N* points) are given by

In a coupled framework, for example, *u*_{wE}, *v*_{wE}, *u*_{wN}, and *v*_{wN} could come from a C-grid ocean model.

As opposed to what is done for the B-grid, the Coriolis term and part of the water stress are explicit (i.e., at iteration *k*) because *u*_{E} and *v*_{N} are not co-located. After introducing ${l}_{E}=\frac{{m}_{E}}{\mathrm{\Delta}{t}_{e}}+{C}_{\mathrm{w}E}\mathrm{cos}{\mathit{\theta}}_{\mathrm{w}}+{C}_{bE}$, ${l}_{N}=\frac{{m}_{N}}{\mathrm{\Delta}{t}_{e}}+{C}_{\mathrm{w}N}\mathrm{cos}{\mathit{\theta}}_{\mathrm{w}}+{C}_{bN}$, ${r}_{E}={m}_{E}{f}_{E}\pm {C}_{\mathrm{w}E}\mathrm{sin}{\mathit{\theta}}_{\mathrm{w}}$, and ${r}_{N}={m}_{N}{f}_{N}\pm {C}_{\mathrm{w}N}\mathrm{sin}{\mathit{\theta}}_{\mathrm{w}}$, Eqs. (A46) and (A47) become

which can be easily solved for ${u}_{E}^{k+\mathrm{1}}$ and ${v}_{N}^{k+\mathrm{1}}$.

The explicit approach for the off-diagonal C-grid terms (as described above) is the same as that used by Kimmritz et al. (2016). Note that for the C-grid, the semi-implicit approach of Bouillon et al. (2009) could be used to solve for *u*^{k+1} and *v*^{k+1} (see their Eqs. 34 and 35).

We conducted many numerical experiments to understand and simplify the conditions that lead to the checkerboard pattern. The goal of this simplification is to allow us to perform a modal analysis and identify the cause of this spurious mode. From the original experiment, with results shown in Fig. 3a, we further simplify the problem by forcing the *v* component of velocity and the shear viscosity, *η*, to be zero. Having *η*=0 is equivalent to setting *e*_{G} to infinity. The ice strength is parameterized according to Hibler (1979). We neglect the snow. We also assume that the concentration is close to 1 and that the ice is in a single thickness category. In experiments for which the wind pushes the ice against a wall, the checkerboard starts to develop close to the wall. This was verified in a few idealized numerical experiments (not shown). We further noticed that when the checkerboard develops close to the wall, there is convergence and the ice is in the plastic regime. Considering the plastic regime and the absence of shear stress, the stress, *σ*=*σ*_{11}, is expressed as

which becomes $\mathit{\sigma}=-P$ because $\mathrm{\Delta}=\left|{D}_{\mathrm{d}}\right|$ and *D*_{d}<0 (convergence).

We write the momentum and transport equations as

where $P={P}^{*}\stackrel{\mathrm{\u203e}}{h}$ as we assume that the concentration is close to 1.

We linearize these equations around *h*_{0} and *u*_{0}. That is, $\stackrel{\mathrm{\u203e}}{h}={h}_{\mathrm{0}}+{h}^{\prime}$ and $u={u}_{\mathrm{0}}+{u}^{\prime}$, where *h*^{′} and *u*^{′} are small perturbations. As $\left|{u}_{\mathrm{0}}\right|\gg \left|{u}^{\prime}\right|$, we write the water stress as ${\mathit{\rho}}_{\mathrm{w}}{C}_{\mathrm{dw}}\left|{u}_{\mathrm{0}}\right|({u}_{\mathrm{0}}+{u}^{\prime})$. Because Eqs. (B2) and (B3) are valid for the base state, *h*_{0} and *u*_{0}, and if we neglect ${h}^{\prime}{u}^{\prime}$ terms, we have

Because the ice is compact and subject to a no-outflow boundary condition, it is reasonable to assume that the base state, *u*_{0}, is small close to the wall and that we can neglect products of *u*_{0} with *h*^{′} and with *u*^{′}. Hence, we neglect the terms $\mathit{\rho}{h}^{\prime}\partial {u}_{\mathrm{0}}/\partial t$, ${\mathit{\rho}}_{\mathrm{w}}{C}_{\mathrm{dw}}\left|{u}_{\mathrm{0}}\right|{u}^{\prime}$, and $\partial \left({h}^{\prime}{u}_{\mathrm{0}}\right)/\partial x$ and finally have

Equations (B6) and (B7) are similar to the one-dimensional shallow-water equations. Many authors have studied these equations and described the checkerboard patterns that depend on the spatial discretization (Schoenstadt, 1980; Batteen and Han, 1981; Walters and Carey, 1983; Le Roux et al., 2005).

We assume solutions of the form ${u}^{\prime}=\widehat{u}{e}^{-i\mathit{\omega}t}$ and ${h}^{\prime}=\widehat{h}{e}^{-i\mathit{\omega}t}$, where *i* is the unit imaginary number and *ω* is the frequency. Following Batteen and Han (1981), we adopt a semi-discrete approach; we only analyze the impact of the spatial discretization and do not consider the time discretization. We first obtain

We write $\widehat{u}=\stackrel{\mathrm{\u0303}}{u}{e}^{i(kx+ly)}$ and $\widehat{h}=\stackrel{\mathrm{\u0303}}{h}{e}^{i(kx+ly)}$, where $\stackrel{\mathrm{\u0303}}{u}$ and $\stackrel{\mathrm{\u0303}}{h}$ define the amplitudes and *k* and *l* the wavenumbers, and we conduct the analysis for a uniform Cartesian grid with grid cells of size Δ*x*×Δ*y*. The origin of our *x* and *y* coordinate system is at the *T* point of a grid cell – that is, the *T* point is at (0,0), while the *E* and *U* points are, respectively, at $(\frac{\mathrm{\Delta}x}{\mathrm{2}},\mathrm{0})$ and $(\frac{\mathrm{\Delta}x}{\mathrm{2}},\frac{\mathrm{\Delta}y}{\mathrm{2}})$. Evaluating Eq. (B8) at the *E* point, we obtain

which can be rearranged as

If the standard remapping (our initial implementation) is used, the departure regions are defined by trapezoids in our simple one-dimensional problem. The shape of these trapezoids depends on the *U*-point velocities, which are calculated from the average C-grid velocities as ${u}_{U}(i,j)=\frac{\left[{u}_{E}\right(i,j)+{u}_{E}(i,j+\mathrm{1}\left)\right]}{\mathrm{2}}$. The area of the trapezoid on the east edge is therefore proportional to $\frac{\left[{u}_{E}\right(i,j-\mathrm{1})+\mathrm{2}{u}_{E}(i,j)+{u}_{E}(i,j+\mathrm{1}\left)\right]}{\mathrm{4}}$. In this case, considering the (perturbed) fluxes for both edges in our simple one-dimensional problem, Eq. (B9) can be written as

which becomes

Using Eq. (B11) to replace ${h}_{\mathrm{0}}\stackrel{\mathrm{\u0303}}{u}$ in Eq. (B13), we obtain the dispersion relation:

Considering the smallest possible wavelength in the *y* direction (*λ*=2Δ*y*), wavenumber *l* is then $l=\mathit{\pi}/\mathrm{\Delta}y$. With that value of *l*, we have *ω*=0 in Eq. (B14), which means that this wave does not propagate: it is a stationary wave, which explains the presence of the checkerboard pattern. This [1+cos (*l*Δ*y*)] term characterizes the spurious divergence associated with the interpolation of velocities to the *U* points. Note that the smallest wavelength in the other direction, *λ*=2Δ*x*, is not a problem because ${\mathrm{sin}}^{\mathrm{2}}\left(\frac{k\mathrm{\Delta}x}{\mathrm{2}}\right)={\mathrm{sin}}^{\mathrm{2}}\left(\frac{\mathit{\pi}}{\mathrm{2}}\right)=\mathrm{1}$.

On the other hand, if the EFA method is used, the fluxes are based on rectangles defined by the *u*_{E} velocity components. Given the fluxes on the west and east edges, Eq. (B9) can be written as

which can be rearranged as follows:

Using Eq. (B11) to replace ${h}_{\mathrm{0}}\stackrel{\mathrm{\u0303}}{u}$ in Eq. (B16), we find the dispersion relation

Compared to Eq. (B14), the dispersion relation associated with the EFA method (Eq. B17) does not have the [1+cos (*l*Δ*y*)] term. As for the ${\mathrm{sin}}^{\mathrm{2}}\left(\frac{k\mathrm{\Delta}x}{\mathrm{2}}\right)$ term in Eq. (B14), the smallest wavelength, *λ*=2Δ*x*, does not create a stationary wave.

Long-term C-grid simulations showed that the novel EFA method with the original remapping algorithm (Lipscomb and Hunke, 2004) sometimes failed on non-uniform grids. These rare failures were due to negative area and mass values close to land or the ice edge.

These negative values were a result of approximations in the area of the departure regions. As explained in Sect. 5.1, the points that define the departure triangles and the shifted departure midpoints are calculated in non-dimensional coordinates. Once the triangles have been found, their areas are scaled to the true grid dimensions, with an area factor, *A*_{f}, assigned to each triangle. This factor is simply an approximation of the grid cell area at a certain location. Triangle areas *A*_{△} are calculated as follows:

where (*x*_{1},*y*_{1}), (*x*_{2},*y*_{2}), and (*x*_{3},*y*_{3}) are the non-dimensional coordinates of the three triangle vertices.

To enhance the robustness of the remapping, the new code modifies some of the area factors. We show two examples to summarize the problems and solutions. In the first example (Fig. C1), we assume that the ocean cell (*i*,*j*) has no ice in category *n* before the transport step. We examine the transport calculation for that category. We assume that cells $(i-\mathrm{1},j)$ and $(i,j+\mathrm{1})$ are land cells, while cells $(i+\mathrm{1},j)$ and $(i,j-\mathrm{1})$ are ocean cells. This means that cell (*i*,*j*) can have fluxes only across its east and south edges. We finally assume that cell $(i+\mathrm{1},j)$ has ice in category *n*. On the south edge (Fig. C1a), the shifted middle departure point (dm^{*}) is in the same cell (*i*,*j*) as the initial middle departure point. This reflects the fact that ${v}_{N}(i,j-\mathrm{1})$ has the same sign as ${v}_{U}(i,j-\mathrm{1})$ (i.e., ${v}_{N}(i,j-\mathrm{1})<\mathrm{0}$ and ${v}_{U}(i,j-\mathrm{1})<\mathrm{0}$). Less commonly, the departure region on the east edge is as shown in Fig. C1b, with ${u}_{E}(i,j)>\mathrm{0}$, ${u}_{U}(i,j-\mathrm{1})<\mathrm{0}$, and $\left|{u}_{E}\right(i,j\left)\right|>\left|{u}_{U}\right(i,j-\mathrm{1}\left)\right|$. In that case, the initial middle departure point, dm (not shown), in cell $(i+\mathrm{1},j)$ is shifted to dm^{*} in cell (*i*,*j*).

On the east edge of cell (*i*,*j*), the orange triangle represents an area flux out of the cell, while the yellow triangle is an incoming flux. On the south edge, all three triangles represent outgoing fluxes. As *a*_{n}=0 in cell (*i*,*j*), the fluxes associated with the orange and the dark blue triangle are zero. The only triangles that matter are the yellow triangles associated with the south and the east edges. The triangle associated with the south edge has vertices cr, dr, ip_{y}, while the one associated with the east edge has vertices cr, dr, ip_{x}. In non-dimensional coordinates, the incoming area flux across the east edge is larger than the outgoing flux across the south edge. This should not lead to a negative net flux for cell (*i*,*j*). However, if the area factor, *A*_{f}, is different for the two yellow triangles, the outgoing area can exceed the incoming area, leading to negative ice area in cell (*i*,*j*).

As described in Sect. 5.1, the EFA method uses the cell area evaluated at the midpoint of the edge to calculate the non-dimensional area of the departure region. Both yellow and orange triangles in Fig. C1b therefore have ${A}_{\mathrm{f}}={A}_{E}(i,j)$, where *A*_{E}(*i*,*j*) is the cell area evaluated at the *E* point. For the south edge, the dark blue and orange triangles have ${A}_{\mathrm{f}}={A}_{N}(i,j-\mathrm{1})$ (Fig. C1a). Because the yellow triangle is not in the central region of the south edge, it is not part of the adjustment process. In the previous version of the remapping, the area factor assigned to this triangle was ${A}_{\mathrm{f}}={A}_{U}(i,j-\mathrm{1})$. On highly deformed grids, ${A}_{U}(i,j-\mathrm{1})$ can be larger than *A*_{E}(*i*,*j*), resulting in negative fluxes.

To improve the robustness of the remapping for C-grid simulations, we modified the code by assigning the area factor *A*_{E}(*i*,*j*) to the yellow triangle associated with the south edge. Similar modifications are required for triangles referred to as TL (top-left), BL (bottom-left), TR (top-right), and BR (bottom-right) in the code. These modifications apply with or without the EFA method.

Differing area factors are less of an issue for the B-grid. Considering the same example, the departure region for the east edge would be defined by points cr, dr, and cl (Fig. C1b). The area flux into cell (*i*,*j*) would be much larger than the area flux out of the cell (yellow triangle in Fig. C1a), and there would be no negative fluxes.

We omit the details here, but a similar problem can arise with lone triangles (e.g., the yellow triangle in Fig. 4). This triangle is now assigned ${A}_{\mathrm{f}}={A}_{N}(i,j)$ instead of *A*_{U}(*i*,*j*) as was done before, in the original code. All the lone triangles now use *A*_{f} at the center of the edge that they border.

As a result of the code modifications described above, another change was required to prevent negative areas. For simplicity, we omit the EFA triangles in this explanation; see Fig. C2. Here, we assume that the ocean cells (*i*,*j*), $(i-\mathrm{1},j)$, and $(i+\mathrm{1},j)$ do not have ice in category *n*. Moreover, there are ocean cells to the north, with *a*_{n}>0, while the southern boundary is a coastline (i.e., land cells). We look at the fluxes for category *n*.

In rare situations, the segment joining dl and dr crosses two edges to form two corner triangles on the north edge of cell (*i*,*j*), as shown in Fig. C2a. Since *a*_{n}=0 in this cell, the departure region inside this cell associated with the north edge does not contribute to the total flux. This region, defined by the points cl, dl, ip_{y}, and cr, is shown in dark blue in Fig. C2a. Similarly, since *a*_{n}=0 in cell $(i+\mathrm{1},j)$, the green triangle for the north edge (Fig. C2a) and the green triangle for the east edge (Fig. C2b) do not contribute to the total flux. The two triangles that matter are the yellow ones. The one for the north edge, defined by cr, dr, ip_{x}, represents a flux out of cell (*i*,*j*), while the one for the east edge defined by cl, dl, ip_{y} corresponds to an incoming flux. In non-dimensional coordinates, the incoming area flux is greater than the outgoing flux. But if the two triangles have different area factors as described above, the net flux can be negative.

With the code changes described above, the yellow triangle for the east edge uses ${A}_{\mathrm{f}}={A}_{N}(i+\mathrm{1},j)$. To ensure robustness (i.e., positive areas) with these changes, the code now uses ${A}_{\mathrm{f}}={A}_{N}(i+\mathrm{1},j)$ for the yellow triangle associated with the north edge (Fig. C2a). This triangle is referred to as TR1 in the code. The green triangle, known as BR1, uses the area factor *A*_{E}(*i*,*j*). Similar modifications are required for the analogous TL1, BL2, and BR2 triangles.

We consider a uniform Cartesian grid with an east–west oriented one-grid-cell-wide channel by applying cyclic boundary conditions. The wind is constant and blows from the west. We further simplify the problem by assuming that the ocean is at rest and that the sea surface tilt term, the turning angle, and the Coriolis parameter are zero. The ice conditions are considered constant along the channel, with *a* being the concentration and $\stackrel{\mathrm{\u203e}}{h}$ the mean thickness. For this test, these fields do not change in time since transport, redistribution, and thermodynamics are turned off. Finally, we do not consider the plastic potential and simply set ${e}_{G}={e}_{F}=e$. With these simplifications, *v*=0 and the *u* momentum equation becomes

Because of the cyclic boundary conditions, $\frac{\partial {\mathit{\sigma}}_{\mathrm{11}}}{\partial x}=\mathrm{0}$. At steady state, the *u* momentum equation therefore becomes

Discretizing Eq. (D2) at the *E* point, we obtain

At steady state, the shear stresses are given by

where ${D}_{S}=\partial u/\partial y$ because *v* is zero. Given the ellipse parameter *e*, *η* is expressed as

where the capping formulation ${\mathrm{\u25b3}}^{*}=max[\mathrm{\u25b3},{\mathrm{\u25b3}}_{\mathrm{min}}]$ is used. Because ${\dot{\mathit{\u03f5}}}_{\mathrm{11}}$ and ${\dot{\mathit{\u03f5}}}_{\mathrm{22}}$ are zero, $\mathrm{\u25b3}={e}^{-\mathrm{1}}\left|{D}_{S}\right|$.

With ${\mathit{\tau}}_{\mathrm{a}E}(i,j)=a{\mathit{\rho}}_{\mathrm{a}}{C}_{\mathrm{da}}{u}_{\mathrm{a}}^{\mathrm{2}}$ and ${\mathit{\tau}}_{\mathrm{w}E}(i,j)=-a{\mathit{\rho}}_{\mathrm{w}}{C}_{\mathrm{dw}}{u}_{E}^{\mathrm{2}}(i,j)$, Eq. (D3) can be written as

where *ρ*_{a} is the air density, *C*_{da} the air drag coefficient, and *u*_{a} the surface wind velocity.

With the no-slip boundary condition, we can approximate the shear strain rate:

which means that ${D}_{\mathrm{s}U}(i,j)<\mathrm{0}$ and ${D}_{\mathrm{s}U}(i,j-\mathrm{1})=-{D}_{\mathrm{s}U}(i,j)$.

We want to solve Eq. (D6) for *u*_{E}(*i*,*j*). For simplicity, we drop (*i*,*j*) – that is, ${u}_{E}(i,j)={u}_{E}$. With strong winds, the ice is in the plastic regime – that is, ${\mathrm{\u25b3}}^{*}=\mathrm{\u25b3}={e}^{-\mathrm{1}}\left|{D}_{\mathrm{s}}\right|$. We can write Eq. (D6) as

The transition between the plastic and viscous regimes occurs for a wind velocity ${u}_{\mathrm{a}}={u}_{\mathrm{a}*}$. At this transition, △=△_{min}, which leads to a sea ice velocity of $e{\mathrm{\u25b3}}_{\mathrm{min}}\mathrm{\Delta}y/\mathrm{2}$. Replacing *u*_{E} in Eq. (D9) by that expression gives

Solving for *u*_{a*}, we get

If ${u}_{\mathrm{a}}>{u}_{\mathrm{a}*}$, the ice is in the plastic regime and *u*_{E} can be found by solving Eq. (D9):

where the first term is the free-drift velocity and the second term, which is there due to the rheology, slows down the ice. In the plastic regime, the shear stresses *σ*_{12U}(*i*,*j*) and ${\mathit{\sigma}}_{\mathrm{12}U}(i,j-\mathrm{1})$ are, respectively, $-{e}^{-\mathrm{1}}P/\mathrm{2}$ and ${e}^{-\mathrm{1}}P/\mathrm{2}$.

On the other hand, if the wind is weak (i.e, ${u}_{\mathrm{a}}<{u}_{\mathrm{a}*}$), the ice is in the viscous regime. In this case, ${\mathrm{\u25b3}}^{*}={\mathrm{\u25b3}}_{\mathrm{min}}$ and Eq. (D6) becomes

which can be rewritten as

The solution of Eq. (D14) is thus

The CICE code is available on GitHub at https://github.com/CICE-Consortium/CICE (last access: 4 September 2024). The simulations for this article were done using release 6.5.0 which can be obtained at https://github.com/CICE-Consortium/CICE/releases/tag/CICE6.5.0 (last access: 4 September 2024) and on Zenodo at https://doi.org/10.5281/zenodo.10056499 (Hunke et al., 2023). Release 6.5.0 includes Icepack 1.4.0. The atmospheric forcing fields (JRA55) and CESM oceanic forcing fields used for the global simulations can be found on Zenodo at https://doi.org/10.5281/zenodo.8118239 (CICE-Consortium, 2020) and https://doi.org/10.5281/zenodo.4660188 (CICE-Consortium, 2021).

JFL led the design and implementation of the C-grid discretization with contributions from AC, DAB, TASR, PB, and EH. AC implemented the test cases with contributions from DAB, EH, DH, JFL, PB, and RA. MB and WHL designed and implemented the first version of the edge flux adjustment method, and JFL and WHL improved the robustness of the method. JFL, FD, and PB investigated the remapping checkerboard pattern. JFL derived the analytical solution for the one-grid-cell-wide channel. JFL wrote the article with contributions from all the co-authors.

The contact author has declared that none of the authors has any competing interests.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

We would like to thank Carolin Mehlmann, Martin Losch, Thomas Richter, and Sergey Danilov for interesting discussions about sea ice model C-grid and CD-grid implementations and Mitchell Bushuk for his help with the SIS2 model. We are also grateful to two anonymous reviewers who provided very useful critiques and comments.

David A. Bailey and William H. Lipscomb were supported by the NSF National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under cooperative agreement no. 1852977. Till A. S. Rasmussen was funded by the Danish state through the National Centre for Climate Research (NCKF). Mats Bentsen was supported by the Research Council of Norway project INES (grant no. 270061). Elizabeth Hunke was supported by the Earth System Model Development program of the US Department of Energy Office of Biological and Environmental Research. Anthony Craig was funded through a National Oceanic and Atmospheric Administration contract in support of the CICE Consortium. David Hebert and Richard Allard were supported by the Office of Naval Research ESPC National Coordination and Infrastructure Support program.

This paper was edited by Christopher Horvat and reviewed by two anonymous referees.

Adcroft, A., Anderson, W., Balaji, V., Blanton, C., Bushuk, M., Dufour, C. O., Dunne, J. P., Griffies, S. M., Hallberg, R., Harrison, M. J., Held, I. M., Jansen, M. F., John, J. G., Krasting, J. P., Langenhorst, A. R., Legg, S., Liang, Z., McHugh, C., Radhakrishnan, A., Reichl, B. G., Rosati, T., Samuels, B. L., Shao, A., Stouffer, R., Winton, M., Wittenberg, A. T., Xiang, B., Zadeh, N., and Zhang, R.: The GFDL Global Ocean and Sea Ice Model OM4.0: Model Description and Simulation Features, J. Adv. Model. Earth Sy., 11, 3167–3211, https://doi.org/10.1029/2019MS001726, 2019. a

Arakawa, A. and Lamb, V. R.: Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, Methods in Computational Physics: Advances in Research and Applications, edited by: Chang, J., 17, 173–265, https://doi.org/10.1016/B978-0-12-460817-7.50009-4, 1977. a

Barton, N., Metzger, E. J., Reynolds, C. A., Ruston, B., Rowley, C., Smedstad, O. M., Ridout, J. A., Wallcraft, A., Frolov, S., Hogan, P., Janiga, M. A., Shriver, J. F., McLay, J., Thoppil, P., Huang, A., Crawford, W., Whitcomb, T., Bishop, C. H., Zamudio, L., and Phelps, M.: The Navy's Earth System Prediction Capability: A New Global Coupled Atmosphere-Ocean-Sea Ice Prediction System Designed for Daily to Subseasonal Forecasting, Earth and Space Science, 8, e2020EA001199, https://doi.org/10.1029/2020EA001199, e2020EA001199 2020EA001199, 2021. a

Batteen, M. L. and Han, Y.-J.: On the computational noise of finite-difference schemes used in ocean models, Tellus, 33, 387–396, https://doi.org/10.1111/j.2153-3490.1981.tb01761.x, 1981. a, b, c

Bouillon, S., Ángel Morales Maqueda, M., Legat, V., and Fichefet, T.: An elastic–viscous–plastic sea ice model formulated on Arakawa B and C grids, Ocean Model., 27, 174–184, https://doi.org/10.1016/j.ocemod.2009.01.004, 2009. a, b, c, d, e

Bouillon, S., Fichefet, T., Legat, V., and Madec, G.: The elastic-viscous-plastic method revisited, Ocean Model., 71, 2–12, https://doi.org/10.1016/j.ocemod.2013.05.013, 2013. a, b, c, d, e, f, g, h

CICE-Consortium: CICE JRA55 gx1 forcing by year. 2020.03.20, Zenodo [data set], https://doi.org/10.5281/zenodo.8118239, 2020. a

CICE-Consortium: CESM climatological 2D ocean forcing gx1 2021.04.02, Zenodo [data set], https://doi.org/10.5281/zenodo.4660188, 2021. a

DeRepentigny, P., Jahn, A., Holland, M. M., and Smith, A.: Arctic Sea Ice in Two Configurations of the CESM2 During the 20th and 21st Centuries, J. Geophys. Res.-Oceans, 125, e2020JC016133, https://doi.org/10.1029/2020JC016133, 2020. a

Dukowicz, J. K. and Baumgardner, J. R.: Incremental Remapping as a Transport/Advection Algorithm, J. Comput. Phys., 160, 318–335, https://doi.org/10.1006/jcph.2000.6465, 2000. a, b

Dupont, F., Dumont, D., Lemieux, J.-F., Dumas-Lefebvre, E., and Caya, A.: A probabilistic seabed–ice keel interaction model, The Cryosphere, 16, 1963–1977, https://doi.org/10.5194/tc-16-1963-2022, 2022. a, b, c

Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, 1979. a, b, c, d, e, f, g

Hunke, E., Allard, R., Bailey, D. A., Blain, P., Craig, A., Dupont, F., DuVivier, A., Grumbine, R., Hebert, D., Holland, M., Jeffery, N., Lemieux, J.-F., Osinski, R., Rasmussen, T., Ribergaard, M., Roach, L., Roberts, A., Turner, M., Winton, M., and Worthen, D.: CICE-Consortium/CICE: CICE Version 6.5.0, Zenodo [code], https://doi.org/10.5281/zenodo.10056499, 2023. a, b

Hunke, E. C.: Viscous-plastic sea ice dynamics with the EVP model: linearization issues, J. Comput. Phys., 170, 18–38, 2001. a, b

Hunke, E. C. and Dukowicz, J. K.: An elastic-viscous-plastic model for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867, 1997. a

Hunke, E. C. and Dukowicz, J. K.: The Elastic-Viscous-Plastic Sea Ice Dynamics Model in General Orthogonal Curvilinear Coordinates on a Sphere - Incorporation of Metric Terms, Mon. Weather Rev., 130, 1848–1865, 2002. a, b

Hutter, N., Bouchat, A., Dupont, F., Dukhovskoy, D., Koldunov, N., Lee, Y. J., Lemieux, J.-F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Ólason, E., Rampal, P., Rasmussen, T., Talandier, C., Tremblay, B., and Wang, Q.: Sea Ice Rheology Experiment (SIREx): 2. Evaluating Linear Kinematic Features in High-Resolution Sea Ice Simulations, J. Geophys. Res.-Oceans, 127, e2021JC017666, https://doi.org/10.1029/2021JC017666, 2022. a

Kay, J. E., Deser, C., Phillips, A., Mai, A., Hannay, C., Strand, G., Arblaster, J. M., Bates, S. C., Danabasoglu, G., Edwards, J., Holland, M., Kushner, P., Lamarque, J.-F., Lawrence, D., Lindsay, K., Middleton, A., Munoz, E., Neale, R., Oleson, K., Polvani, L., and Vertenstein, M.: The Community Earth System Model (CESM) Large Ensemble Project: A Community Resource for Studying Climate Change in the Presence of Internal Climate Variability, B. Am. Meteorol. Soc., 96, 1333–1349, https://doi.org/10.1175/BAMS-D-13-00255.1, 2015. a

Kimmritz, M., Danilov, S., and Losch, M.: On the convergence of the modified Elastic-Viscous-Plastic method for solving the sea ice momentum equation, J. Comput. Phys., 296, 90–100, https://doi.org/10.1016/j.jcp.2015.04.051, 2015. a

Kimmritz, M., Danilov, S., and Losch, M.: The adaptive EVP method for solving the sea ice momentum equation, Ocean Model., 101, 59–67, https://doi.org/10.1016/j.ocemod.2016.03.004, 2016. a, b, c, d, e, f, g

Kobayashi, S., Ota, Y., Harada, Y., Ebita, A., Moriya, M., Onoda, H., Onogi, K., Kamahori, H., Kobayashi, C., Endo, H., Miyaoka, K., and Takahashi, K.: The JRA-55 Reanalysis: General Specifications and Basic Characteristics, J. Meteorol. Soc. Jpn. Ser. II, 93, 5–48, https://doi.org/10.2151/jmsj.2015-001, 2015. a

Koldunov, N. V., Danilov, S., Sidorenko, D., Hutter, N., Losch, M., Goessling, H., Rakowsky, N., Scholz, P., Sein, D., Wang, Q., and Jung, T.: Fast EVP Solutions in a High-Resolution Sea Ice Model, J. Adv. Model. Earth Sy., 11, 1269–1284, https://doi.org/10.1029/2018MS001485, 2019. a

König Beatty, C. and Holland, D. M.: Modeling landfast sea ice by adding tensile strength, J. Phys. Oceanogr., 40, 185–198, https://doi.org/10.1175/2009JPO4105.1, 2010. a, b

Kreyscher, M., Harder, M., Lemke, P., and Flato, G. M.: Results of the Sea Ice Model Intercomparison Project: Evaluation of sea ice rheology schemes for use in climate simulations, J. Geophys. Res., 105, 11299–11320, 2000. a, b

Le Roux, D., Sène, A., Rostand, V., and Hanert, E.: On some spurious mode issues in shallow-water models using a linear algebra approach, Ocean Model., 10, 83–94, https://doi.org/10.1016/j.ocemod.2004.07.008, 2005. a

Lemieux, J.-F. and Dupont, F.: On the calculation of normalized viscous–plastic sea ice stresses, Geosci. Model Dev., 13, 1763–1769, https://doi.org/10.5194/gmd-13-1763-2020, 2020. a, b, c

Lemieux, J.-F., Tremblay, B., Thomas, S., Sedláček, J., and Mysak, L. A.: Using the preconditioned Generalized Minimum RESidual (GMRES) method to solve the sea-ice momentum equation, J. Geophys. Res., 113, C10004, https://doi.org/10.1029/2007JC004680, 2008. a

Lemieux, J.-F., Knoll, D. A., Tremblay, B., Holland, D. M., and Losch, M.: A comparison of the Jacobian-free Newton Krylov method and the EVP model for solving the sea ice momentum equation with a viscous-plastic formulation: a serial algorithm study, J. Comput. Phys., 231, 5926–5944, https://doi.org/10.1016/j.jcp.2012.05.024, 2012. a

, Lemieux, J.-F., Tremblay, L. B., Dupont, F., Plante, M., Smith, G. C., and Dumont, D.: A basal stress parameterization for modeling landfast ice, J. Geophys. Res.-Oceans, 120, 3157–3173, https://doi.org/10.1002/2014JC010678, 2015. a

Lemieux, J.-F., Dupont, F., Blain, P., Roy, F., Smith, G. C., and Flato, G. M.: Improving the simulation of landfast ice by combining tensile strength and a parameterization for grounded ridges, J. Geophys. Res.-Oceans, 121, 7354–7368, https://doi.org/10.1002/2016JC012006, 2016. a, b, c, d, e, f, g

Lipscomb, W. H. and Hunke, E. C.: Modeling sea ice transport using incremental remapping, Mon. Weather Rev., 132, 1341–1354, 2004. a, b, c, d, e, f, g, h, i

Lipscomb, W. H., Hunke, E. C., Maslowski, W., and Jakacki, J.: Ridging, strength, and stability in high-resolution sea ice models, J. Geophys. Res., 112, C03S91, https://doi.org/10.1029/2005JC003355, 2007. a

Losch, M., Menemenlis, D., Campin, J.-M., Heimbach, P., and Hill, C.: On the formulation of sea-ice models. Part 1: Effects of different solver implementations and parameterizations, Ocean Model., 33, 129–144, https://doi.org/10.1016/j.ocemod.2009.12.008, 2010. a

Madec, G. and the NEMO System Team: NEMO Ocean Engine Reference Manual, Zenodo, https://doi.org/10.5281/zenodo.8167700, 2023. a

McTaggart-Cowan, R., Vaillancourt, P. A., Zadra, A., Chamberland, S., Charron, M., Corvec, S., Milbrandt, J. A., Paquin-Ricard, D., Patoine, A., Roch, M., Separovic, L., and Yang, J.: Modernization of Atmospheric Physics Parameterization in Canadian NWP, J. Adv. Model. Earth Sy., 11, 3593–3635, https://doi.org/10.1029/2019MS001781, 2019. a

Metzger, E. J., Smedstad, O. M., Thoppil, P. G., Hurlburt, H. E., Cummings, J. A., Wallcraft, A. J., Zamudio, L., Franklin, D. S., Posey, P. G., Phelps, M. W., Hogan, P. J., Bub, F. L., and DeHaan, C. J.: US Navy Operational Global Ocean and Arctic Ice Prediction Systems, Oceanography, 27, 32–43, https://doi.org/10.5670/oceanog.2014.66, 2014. a

Prather, M. J.: Numerical advection by conservation of second-order moments, J. Geophys. Res.-Atmos., 91, 6671–6681, 1986. a

Ringeisen, D., Losch, M., Tremblay, L. B., and Hutter, N.: Simulating intersection angles between conjugate faults in sea ice with different viscous–plastic rheologies, The Cryosphere, 13, 1167–1186, https://doi.org/10.5194/tc-13-1167-2019, 2019. a

Ringeisen, D., Tremblay, L. B., and Losch, M.: Non-normal flow rules affect fracture angles in sea ice viscous–plastic rheologies, The Cryosphere, 15, 2873–2888, https://doi.org/10.5194/tc-15-2873-2021, 2021. a, b

Roach, L. A., Horvat, C., Dean, S. M., and Bitz, C. M.: An Emergent Sea Ice Floe Size Distribution in a Global Coupled Ocean-Sea Ice Model, J. Geophys. Res.-Oceans, 123, 4322–4337, https://doi.org/10.1029/2017JC013692, 2018. a

Rothrock, D. A.: The Energetics of the Plastic Deformation of Pack Ice by Ridging, J. Geophys. Res., 80, 4514–4519, 1975. a, b

Schoenstadt, A. L.: A transfer function analysis of numerical schemes Used to simulate geostrophic adjustment, Mon. Weather Rev., 108, 1248–1259, 1980. a

Smith, G. C., Liu, Y., Benkiran, M., Chikhar, K., Surcel Colan, D., Gauthier, A.-A., Testut, C.-E., Dupont, F., Lei, J., Roy, F., Lemieux, J.-F., and Davidson, F.: The Regional Ice Ocean Prediction System v2: a pan-Canadian ocean analysis system using an online tidal harmonic analysis, Geosci. Model Dev., 14, 1445–1467, https://doi.org/10.5194/gmd-14-1445-2021, 2021. a

Tsamados, M., Feltham, D., and Wilchinsky, A. V.: Impact of a new anisotropic rheology on simulations of Arctic sea ice, J. Geophys. Res., 118, 91–107, https://doi.org/10.1029/2012JC007990, 2013. a

Vancoppenolle, M., Rousset, C., Blockley, E., Aksenov, Y., Feltham, D., Fichefet, T., Garric, G., Guémas, V., Iovino, D., Keeley, S., Madec, G., Massonnet, F., Ridley, J., Schroeder, D., and Tietsche, S.: SI3, the NEMO Sea Ice Engine (4.2release_doc1.0), Zenodo, https://doi.org/10.5281/zenodo.7534900, 2023. a

Walters, R. A. and Carey, G. F.: Analysis of spurious oscillation modes for the shallow water and Navier-Stokes equations, Comput. Fluids, 11, 51–68, https://doi.org/10.1016/0045-7930(83)90013-0, 1983. a

Zhang, J. and Rothrock, D. A.: Effect of sea ice rheology in numerical investigations of climate, J. Geophys. Res., 110, C08014, https://doi.org/10.1029/2004JC002599, 2005. a

- Abstract
- Introduction
- Model equations
- Spatial and temporal discretizations of the momentum equation
- Initial approach for remapping using C-grid velocity components
- The edge flux adjustment method
- Validation of the C-grid implementation
- Conclusions
- Appendix A: Spatial and temporal discretization of the momentum equation
- Appendix B: Modal analysis of the remapping checkerboard pattern
- Appendix C: Improved robustness of remapping
- Appendix D: C-grid analytical solution for a one-grid-cell-wide channel
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Model equations
- Spatial and temporal discretizations of the momentum equation
- Initial approach for remapping using C-grid velocity components
- The edge flux adjustment method
- Validation of the C-grid implementation
- Conclusions
- Appendix A: Spatial and temporal discretization of the momentum equation
- Appendix B: Modal analysis of the remapping checkerboard pattern
- Appendix C: Improved robustness of remapping
- Appendix D: C-grid analytical solution for a one-grid-cell-wide channel
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References