Implementation of a brittle sea ice rheology in an Eulerian, finite-difference, C-grid modeling framework: impact on the simulated deformation of sea ice in the Arctic

We have implemented the brittle Bingham–Maxwell sea ice rheology (BBM) into SI3, the sea ice component of NEMO. After discussing the numerical aspects and requirements that are specific to the implementation of a brittle rheology in the Eulerian, finite-difference, Arakawa C-grid framework, we detail the approach we have used. This approach relies on the introduction of an additional set of prognostic stress tensor components, sea ice damage, and sea ice velocity vector, following a grid point arrangement that expands the C-grid into the Arakawa E-grid. The newly implemented BBM rheology is first assessed by means of a set of idealized SI3 simulations at different spatial resolutions. Then, sea ice deformation rates obtained from simulations of the Arctic at a $\mathrm{1}/\mathrm{4}$° spatial resolution, performed with the coupled ocean–sea ice setup of NEMO, are assessed against satellite observations. For all these simulations, results obtained with the default current workhorse setup of SI3 are provided to serve as a reference. Our results show that using a brittle type of rheology, such as BBM, allows SI3 to simulate the highly localized deformation pattern of sea ice, as well as its scaling properties, from the scale of the model's computational grid up to the basin scale.

1 Introduction

Sea ice impacts the ocean and the atmosphere at both the local and the regional scale (Vihma, 2014; IPCC, 2022). In polar regions, the sea ice cover modulates the radiative and turbulent exchanges of heat, fresh water, gas, and momentum between the ocean and the atmosphere (e.g., Taylor et al., 2018, for a review). At the local scale, these fluxes are predominately controlled by the heterogeneity of the sea ice thickness, which itself is governed by the sea ice dynamics and the associated formation of leads and ridges. This stresses the relevance of accurately representing sea ice dynamics in simulations of the coupled multi-component Earth system, such as regional and global climate simulations, and even for short-term sea ice predictions.

The dynamical behavior of sea ice is controlled by processes that interact and evolve over a wide range of spatial and temporal scales. This multiscale nature of sea ice physics is fascinating and has triggered the curiosity of geophysicists since the early 1970s (Coon et al., 1974). More recently, scientific interest in sea ice dynamics has grown significantly due to the dramatic retreat and thinning of the Arctic sea ice cover. In addition, the abundance of new observations of sea ice kinematics, recorded by both in situ instruments (e.g., the buoy trajectories of the International Arctic Buoy Program, https://iabp.apl.uw.edu/data.html, last access: 30 July 2024) and satellites (e.g., the ice trajectories from the RADARSAT Geophysical Processor System; Kwok et al., 1998), has the potential to foster additional advancements in sea ice modeling.

The dynamics of sea ice is complex. For instance, Rampal et al. (2008) and Weiss et al. (2009) showed that the statistical properties of sea ice deformation are characterized by a coupled space–time multi-fractal scaling invariance, similar to what is observed for the deformation of the Earth's crust (Kagan and Jackson, 1991; Marsan and Weiss, 2010). The spatial and temporal scaling properties of sea ice deformation and their coupling provide evidence for the strong heterogeneity and intermittency that characterizes sea ice dynamics (Rampal et al., 2008).

Reproducing the discontinuous nature of sea ice – related to the presence of fractures and leads – in continuous sea ice models, as well as the complexity of the spatial patterns and temporal evolution of these features, poses a fundamental and major challenge (e.g., Bouchat et al., 2022; Hutter et al., 2022). As an effort to tackle this challenge, Dansereau et al. (2016) introduced the Maxwell elasto-brittle rheology (MEB). MEB was implemented into neXtSIM, a large-scale dynamical–thermodynamical Lagrangian finite-element sea ice model (Rampal et al., 2016), and used to evaluate the performance of this new rheology in a realistic pan-Arctic simulation. Wintertime sea ice deformations simulated by MEB in neXtSIM were first evaluated statistically against satellite observations, in terms of probability density functions (PDFs) and scaling invariance properties, in Rampal et al. (2016, 2019) and later in the two companion papers of Bouchat et al. (2022) and Hutter et al. (2022), showing satisfying results.

Recently, Ólason et al. (2022) introduced the brittle Bingham–Maxwell rheology (BBM) as an effort to address the incomplete treatment of the convergence of highly damaged sea ice in MEB. This deficiency of MEB results in the unrealistic representation of the ice thickness after a couple of years of model integration. Indeed, recent realistic BBM-driven multidecadal simulations performed with neXtSIM have been shown to reproduce (i) the scaling properties of sea ice deformation from the model grid cell up to the scale of the Arctic basin and (ii) the thickness pattern of the sea ice cover well when compared to observations (Ólason et al., 2022; Boutin et al., 2023).

Yet, performing coupled ocean–sea ice or Earth system CMIP-like simulations with neXtSIM in a numerically efficient manner remains challenging because the numerical coupling of neXtSIM to a third-party – generally Eulerian – GCM component requires the implementation of a relatively inefficient Lagrangian–Eulerian coupling strategy. Furthermore, the weak scalability capabilities of neXtSIM when run in parallel on more than a few tenths of processors and/or at spatial resolutions typically below 10 km have been shown to substantially hinder the scalability of coupled setups (Samaké et al., 2017). Thus, the implementation of BBM into an Eulerian CMIP-class sea ice model, such as SI3 of NEMO, has the potential to significantly benefit the sea ice, ocean, and climate modeling communities. First, it will facilitate the assessment of the sensitivity of the simulated sea ice dynamics to the type of rheology used in a modeling framework that these communities are familiar with. And second, the good scalability capabilities of NEMO (Tintó Prims et al., 2019) will allow performing realistic kilometer-scale simulations that use a brittle rheology.

As of today, a few efforts have been made to implement MEB in sea ice models comparable to SI3 in terms of discretization method and grid, such as the MIT general circulation model (Losch et al., 2010), or LIM, the former sea ice component of the NEMO modeling system (Rousset et al., 2015). And more recently, Plante et al. (2020) have successfully implemented MEB in the McGill sea ice model (Tremblay and Mysak, 1997; Lemieux et al., 2008, 2014). Overall, the work of these modeling groups has highlighted some challenging aspects that are specific to the implementation of a brittle rheology in a realistic Eulerian model that uses the finite-difference method on a staggered grid. As suggested by the work of Plante et al. (2020), when discretized on the Arakawa C-grid (Arakawa and Lamb, 1977), the same grid as used by SI3 (Vancoppenolle et al., 2023), brittle rheologies seem to be more prone to numerical instabilities than their viscous–plastic counterparts. In particular, they report that the stability of their MEB implementation is sensitive to the resort to spatial averaging, an interpolation technique that is traditionally used to relocate certain fields between the staggered points of the grid. Moreover, the need to advect the stress tensor, specific to brittle rheologies, poses another challenge when using the C-grid because it demands the advection of a scalar field, namely the shear element of the stress tensor, that is defined at the corner points of the grid cell.

In this paper, we propose a new discretization approach adapted to the numerical implementation of a brittle rheology in an Eulerian finite-difference, C-grid-based sea ice model. We describe how we have implemented this new approach into SI3 based on the expansion of the C-grid into an Arakawa E-grid.

As a first validation step of our BBM implementation, we discuss SI3 results obtained using the idealized test case setup of Mehlmann et al. (2021) at different horizontal resolutions. Then, as the second step, we compare the sea ice deformations obtained from realistic coupled ocean–sea ice simulations of the pan-Arctic against those constructed from satellite observations. To serve as a reference, the results of simulations run using the default workhorse setup of SI3 (based on the aEVP rheology of Kimmritz et al., 2016) are also included in both validation steps.

This paper is organized as follows. In Sect. 2, we summarize the equations of the sea ice dynamics model, discuss the aspects in which the numerical implementation of a brittle rheology may differ from that of a non-brittle viscous–plastic one, and detail the numerical aspects specific to our implementation of BBM into SI3. In Sect. 3, we describe the technical aspects of our SI3 simulations and discuss the results obtained with both the idealized and pan-Arctic configurations. In Sect. 4, we discuss some numerical aspects of our implementation and some limitations of the BBM rheology. Our conclusions are summarized in Sect. 5. Detailed nomenclature relating the acronyms and symbols used throughout the paper is outlined in Appendix A.

2 Model and implementation

2.1 Governing equations and constitutive law

The two-dimensional, vertically integrated, momentum equation for sea ice reads

$$\begin{array}{}\text{(1)}& m{\partial }_{\mathrm{t}}\mathit{u}=\mathbf{\nabla }\cdot \left(h\mathit{\sigma }\right)+A({\mathit{\tau }}_{\mathrm{a}}+\mathit{\tau })-mf\mathit{k}\times \mathit{u}-mg\mathbf{\nabla }H,\end{array}$$

where m is the mass of ice and snow per unit area, u is the ice velocity vector, h is the ice thickness, σ is the internal stress tensor, A is the sea ice fraction, τ_a is the wind stress vector, τ is the ocean current stress, f is the Coriolis frequency, k is the vertical unit vector, g is the acceleration of gravity, and H is the sea surface height. In the two-dimensional (plane stresses) case, the stress tensor is written as

$$\begin{array}{}\text{(2)}& \mathit{\sigma }=\left(\begin{array}{cc}{\mathit{\sigma }}_{\mathrm{11}}& {\mathit{\sigma }}_{\mathrm{12}}\\ {\mathit{\sigma }}_{\mathrm{12}}& {\mathit{\sigma }}_{\mathrm{22}}\end{array}\right).\end{array}$$

In general, a constitutive law relates σ to the strain rate tensor $\dot{\mathit{\epsilon }}$, defined as follows:

$$\begin{array}{}\text{(3)}& \dot{\mathit{\epsilon }}=\left(\begin{array}{cc}{\dot{\mathit{\epsilon }}}_{\mathrm{11}}& {\dot{\mathit{\epsilon }}}_{\mathrm{12}}\\ {\dot{\mathit{\epsilon }}}_{\mathrm{12}}& {\dot{\mathit{\epsilon }}}_{\mathrm{22}}\end{array}\right)\equiv \left(\begin{array}{cc}{\partial }_{x}u& \frac{\mathrm{1}}{\mathrm{2}}({\partial }_{y}u+{\partial }_{x}v)\\ \frac{\mathrm{1}}{\mathrm{2}}({\partial }_{y}u+{\partial }_{x}v)& {\partial }_{y}v\end{array}\right).\end{array}$$

As derived by Ólason et al. (2022) (their Eq. 20), the BBM constitutive model yields

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\partial }_{\mathrm{t}}\underset{\mathrm{¯}}{\mathit{\sigma }}=E\underset{\mathrm{¯}}{\mathit{K}}\cdot \underset{\mathrm{¯}}{\dot{\mathit{\epsilon }}}-\underset{\mathrm{¯}}{\mathit{\sigma }}{\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}\left(\mathrm{1}+\stackrel{\mathrm{̃}}{P}+{\displaystyle \frac{\mathit{\lambda }}{\mathrm{1}-d}}{\partial }_{\mathrm{t}}d\right)\\ & \text{ with }\underset{\mathrm{¯}}{\mathit{K}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}}\left(\begin{array}{ccc}\mathrm{1}& \mathit{\nu }& \mathrm{0}\\ \mathit{\nu }& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}-\mathit{\nu }\end{array}\right),\end{array}\end{array}$$

where E and λ are the elastic modulus and apparent viscous relaxation time of the ice, K is the elastic stiffness tensor, $\stackrel{\mathrm{̃}}{P}$ is a term introduced to prevent excessive ridging (see below), and d is the damage scalar: a variable that represents the density of fractures in the ice at the subgrid scale. The underbar notation indicates that the tensors are expressed in their Voigt form. E and λ are modulated by the sea ice concentration and damage as

$$\begin{array}{}\text{(5)}& {\displaystyle }E=& {\displaystyle }{E}_{\mathrm{0}}(\mathrm{1}-d)e^{-C(\mathrm{1}-A)},\text{(6)}& {\displaystyle }\mathit{\lambda }=& {\displaystyle }{\mathit{\lambda }}_{\mathrm{0}}\left[\right(\mathrm{1}-d)e^{-C(\mathrm{1}-A)}{]}^{\mathit{\alpha }-\mathrm{1}},\end{array}$$

where C is the compaction parameter constant and α is a constant exponent greater than 1. α fulfills the physical constraint that the relaxation time for the stress also decreases as damage increases and re-increases as the ice heals (i.e., damage decreases) because the material respectively loses and recovers the memory of reversible deformations (Dansereau et al., 2016).

The BBM constitutive model in Eq. (4) only differs from that of MEB through the inclusion of the term $\stackrel{\mathrm{̃}}{P}$: a threshold between reversible and permanent deformation regimes. As noted by Ólason et al. (2022), the inclusion of this term prevents the excessive convergence that is occurring in MEB simulations lasting longer than a season. For convergent stresses in the range $-P_{\text{max}}<{\mathit{\sigma }}_{\mathrm{I}}<\mathrm{0}$, the deformation is elastic; otherwise, it is visco-elastic. Ólason et al. (2022) interpret this threshold as the maximum pressure the ice can withstand before ridging. They consequently choose to let the ridging threshold, P_max, be proportional to the ice thickness to the power $\mathrm{3}/\mathrm{2}$ (Hopkins, 1998) and depend exponentially on the concentration (Hibler, 1979), i.e.,

$$\begin{array}{}\text{(7)}& \stackrel{\mathrm{̃}}{P}=\left\{\begin{array}{ll}\mathrm{0}& \text{if }{\mathit{\sigma }}_{\mathrm{I}}>\mathrm{0}\\ -\mathrm{1}& \text{if }-P_{\text{max}}<{\mathit{\sigma }}_{\mathrm{I}}<\mathrm{0}\\ \frac{P_{\text{max}}}{{\mathit{\sigma }}_{\mathrm{I}}}& \text{if }{\mathit{\sigma }}_{\mathrm{I}}<-P_{\text{max}}\end{array}\right.,\end{array}$$

where σ_I is the (isotropic) normal stress and P_max is the ridging threshold defined as

$$\begin{array}{}\text{(8)}& P_{\text{max}}=P_{\mathrm{0}}{\left({\displaystyle \frac{h}{h_{\mathrm{0}}}}\right)}^{\mathrm{3}/\mathrm{2}}{e}^{-C(\mathrm{1}-A)}.\end{array}$$

We follow Dansereau et al. (2016) and Ólason et al. (2022) in using a two-step approach to solve Eq. (4). As the first step, an initial estimate of $\underset{\mathrm{¯}}{\mathit{\sigma }}$, denoted as ${\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(\mathrm{i}\right)}$, is calculated assuming no change in damage:

$$\begin{array}{}\text{(9)}& {\partial }_{\mathrm{t}}{\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(\mathrm{i}\right)}=E\mathit{K}\cdot \underset{\mathrm{¯}}{\dot{\mathit{\epsilon }}}-\underset{\mathrm{¯}}{\mathit{\sigma }}{\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}(\mathrm{1}+\stackrel{\mathrm{̃}}{P}).\end{array}$$

Then, as the second step, the following test and adjustment are performed on the state of stress: if ${\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(\mathrm{i}\right)}$ is locally over-critical, i.e., located outside of the Mohr–Coulomb damage criterion (Fig. 1), an increment in ice damage, d_crit, is applied such that

$$\begin{array}{}\text{(10)}& \underset{\mathrm{¯}}{\mathit{\sigma }}=d_{\text{crit}}{\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(\mathrm{i}\right)},\end{array}$$

where ${\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(\mathrm{i}\right)}$ is the local value of the over-critical stress, and $\underset{\mathrm{¯}}{\mathit{\sigma }}$ is the corresponding post-failure (i.e., post-damage) stress. As discussed in Dansereau et al. (2016) and Plante and Tremblay (2021), d_crit is used to scale over-critical stresses back towards the Mohr–Coulomb damage criterion, assuming viscous relaxation to be negligible during the (comparatively very fast) damage process. The associated temporal evolution of the damage and adjustment of the stress state is given by

$$\begin{array}{}\text{(11)}& {\displaystyle }{\partial }_{\mathrm{t}}d& {\displaystyle }={\displaystyle \frac{\mathrm{1}-d_{\text{crit}}}{t_{\mathrm{d}}}}(\mathrm{1}-d),\text{(12)}& {\displaystyle }{\partial }_{\mathrm{t}}\underset{\mathrm{¯}}{\mathit{\sigma }}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}-d_{\text{crit}}}{t_{\mathrm{d}}}}{\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(\mathrm{i}\right)},\end{array}$$

where t_d is a characteristic timescale for damage propagation. In the BBM framework, d_crit is expressed as follows:

$$\begin{array}{}\text{(13)}& d_{\text{crit}}=\left\{\begin{array}{ll}\frac{c}{{\mathit{\sigma }}_{\text{II}}^{\left(\mathrm{i}\right)}+\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\left(\mathrm{i}\right)}}& \text{if }{\mathit{\sigma }}_{\mathrm{I}}^{\left(\mathrm{i}\right)}\ge -N\\ \frac{-N}{{\mathit{\sigma }}_{\mathrm{I}}^{\left(\mathrm{i}\right)}}& \text{otherwise},\end{array}\right.\end{array}$$

where c is the cohesion, and μ is the friction coefficient. The threshold N is used to prevent any numerical instability at very high normal stresses and is set large enough not to impact the solution noticeably. As suggested by Rampal et al. (2016), a slow restoring process is applied to the damage to account for the healing of ice under refreezing conditions. The rate of decrease of the damage associated with this refreezing is taken to be proportional to ΔT_h, the temperature difference between basal and surface ice:

$$\begin{array}{}\text{(14)}& {\partial }_{\mathrm{t}}d=-{\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{h}}}{k_{\text{th}}}},\end{array}$$

where k_th is the healing constant. This process can be decoupled from Eq. (11) due to the large separation of timescales between the healing and damaging processes.

Figure 1 Mohr–Coulomb yield envelope in the internal stress invariant coordinates (blue line). Illustration of how an over-critical stress state σ⁽ⁱ⁾ (initial estimate) is evolving (gray arrow) towards the corrected state σ when using the BBM rheology.

2.2 Numerical implementation: brittle versus viscous–plastic rheologies

To understand the extent to which the numerical implementation of a brittle rheology differs from that of a viscous–plastic one, let us first review the main differences between these rheologies and their respective classical numerical implementation.

First, the elasto-visco-brittle family of rheologies (MEB, BBM; Dansereau et al., 2016; Ólason et al., 2022) considers unfragmented sea ice to be an elastic and damageable solid. Fragmented sea ice is a visco-elastic material in which irreversible deformations dissipate the stresses. As opposed to viscous–plastic frameworks, elasticity is therefore a physical and non-negligible component of the model. It is modulated by the level of damage, d, which keeps the memory of the state of fragmentation of the sea ice cover. The combination of elasticity and damage, even if treated in an isotropic manner, naturally simulates a strong anisotropy and localization of the deformation, down to the nominal spatial and temporal scale (i.e., the grid resolution and time step of the model, respectively; Dansereau et al., 2016; Weiss and Dansereau, 2017; Rampal et al., 2019; Ólason et al., 2022). Therefore, all the mechanically related fields, such as damage, concentration, thickness, and velocity, tend to exhibit very sharp gradients. Second, in the BBM (as in the MEB) framework, a twofold approach is used to linearize the system of equations and solve the coupled constitutive and damage evolution equations: (i) an initial estimate, in which stress components are updated based on the constitutive law (Eq. 9), and (ii) a damage step in which the Mohr–Coulomb test is performed, resulting in a potential adjustment of local over-critical stresses and the associated increase in damage (Fig. 1, Eqs. 11 and 12). In viscous–plastic rheologies, which do not incorporate damage, no such twofold approach is necessary to solve the system of dynamical equations.

Figure 2 Point arrangement and staggering in a grid cell: (a) the C-grid as used in NEMO and (b) the E-grid. The letter d indicates the location of tracers, while the letters u and v indicate that of the i- and j-wise components of the velocity vector. Letters in brackets indicate the name of the grid points as referred to throughout the paper.

A third and major difference between the two types of model is that in brittle models, the stress tensor σ is a prognostic variable, while it is a diagnostic variable in viscous–plastic models. This implies that the implementation of a brittle rheology in an Eulerian framework, as opposed to that of a viscous–plastic rheology, should, in practice, consider the advection of σ, along with other – typically scalar – tracers (see Sect. 2.4). One could argue that, based on a scale analysis, the advection terms of the stress tensor components are somewhat negligible and that it is therefore acceptable to simply omit these terms (similarly to what is done for the ice velocity in the momentum equations). While these terms are indeed very small, we think that it is important to include them in the Eulerian implementation of a brittle model because of the strong interdependence that bounds the damage tracer and the stress tensor. This interdependence is the consequence of E and λ being a function of d (Eqs. 5 and 6) in the estimate of σ⁽ⁱ⁾ (Eq. 9). The damage, as a tracer that can live on for days, if not weeks, depending on the temperature conditions, has to be advected with the ice velocity, like any other tracer. As such, we think that the advection of the stress tensor is necessary to preserve the full spatial consistency between the damage tracer and the internal stress state, in particular in the case of simulations longer than a few days that involve significant sea ice displacements.

Finally, note that in their numerical implementation of BBM, Ólason et al. (2022) chose to solve the dynamics explicitly using a time step sufficiently small to account for the propagation of damage in the ice in a physically realistic manner. We follow the same approach in our implementation of BBM into SI3. Typically, this implies using a time step a few hundred times smaller (hereafter referred to as the dynamical time step, Δt) than that used for the thermodynamics and the advection (hereafter referred to as the advective time step, ΔT). This is implemented by means of a time-splitting approach. N_s, the number of (Δt-long) integrations to perform during one advective time step (ΔT), is imposed by ΔT and Δx, the horizontal resolution of the grid:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}}>\mathrm{2}{C}_{\mathrm{E}}\Rightarrow & N_{\mathrm{s}}>\mathrm{2}{C}_{\mathrm{E}}{\displaystyle \frac{\mathrm{\Delta }T}{\mathrm{\Delta }x}},\\ & \text{ with }{N}_{\mathrm{s}}\equiv {\displaystyle \frac{\mathrm{\Delta }T}{\mathrm{\Delta }t}}\text{ and }{C}_{\mathrm{E}}=\sqrt{{\displaystyle \frac{E}{\mathrm{2}(\mathrm{1}+\mathit{\nu }){\mathit{\rho }}_{\mathrm{i}}}}},\end{array}\end{array}$$

where C_E is the propagation speed of an elastic shear wave and ρ_i is the density of ice. Note that if ΔT is already constrained by Δx, as in NEMO, the choice of N_s becomes somewhat independent of the spatial resolution at which the model is run.

2.3 Numerical implementation of the BBM rheology

SI3 uses curvilinear coordinates on a fixed Eulerian mesh, and the spatial discretization is achieved by means of the finite-difference method (FD) on the Arakawa C-grid (Arakawa and Lamb, 1977). The use of the C-grid is justified based on numerical and practical grounds, as it ensures the exact collocation of ocean and sea ice horizontal velocity components, which simplifies the coupling with the ocean component of NEMO and prevents interpolation-related errors as well as extra computational load.

As shown in Fig. 2a, on the C-grid, tracers are defined at the cell centers, hereafter referred to as the T-point, while the x and y components of vectors are defined at the center of the right-hand and upper edges of each cell, respectively (hereafter U-point and V-point). The point located at the upper-right corner of each cell, known as the vorticity point, is referred to as the F-point. In the literature, this vorticity point is sometimes located at the bottom-left corner of the cell, with the U- and V-points possibly located at the left-hand and lower edges of the cell, and is sometimes referred to as the Z-point (Losch et al., 2010; Plante et al., 2020).

Regardless of the type of rheology considered, the main challenge posed by the C-grid is a consequence of the discretized FD expressions of the elements of the strain rate tensor $\dot{\mathit{\epsilon }}$ being staggered in space, with the trace elements ${\dot{\mathit{\epsilon }}}_{\mathrm{11}}$ and ${\dot{\mathit{\epsilon }}}_{\mathrm{22}}$ defined at the T-point and the shear rate ${\dot{\mathit{\epsilon }}}_{\mathrm{12}}$ defined at the F-point. Based on the constitutive law (Eq. 4), the same applies to the stress tensor σ. This staggering between the diagonal and off-diagonal elements of σ is appropriate when considering the discretization of the momentum equation (Eq. 1) because the discretized elements of the vector divergence of σ are then defined where they are needed: namely, at U- and V-points. However, this staggering becomes an issue whenever the parameterization of the constitutive law requires ${\dot{\mathit{\epsilon }}}_{\mathrm{12}}$ or σ₁₂ to be known at a T-point. This is the case, for instance, for the expression of the Δ parameter in elastic–viscous–plastic (EVP) models or that of the second stress invariant σ_II in MEB and BBM (i.e., Eq.13), as they require ${\dot{\mathit{\epsilon }}}_{\mathrm{12}}$ and σ₁₂, respectively, to be known at T-points. Moreover, in brittle rheologies, a value of d is required not only at the T-point, but also at the F-point to estimate E and λ (Eqs. 5 and 6) needed to update ${\mathit{\sigma }}_{\mathrm{12}}^{\left(\mathrm{i}\right)}$ (Eq. 9). On the C-grid, a common way to interpolate a scalar defined at F-points onto T-points is to simply use the average of this scalar on the four surrounding F- points and conversely to interpolate from T- to F-points. In the aEVP implementation of SI3 (Kimmritz et al., 2016), the problem posed by the staggering of tensor elements is overcome by using this averaging approach to interpolate the square of the shear rate ${\dot{\mathit{\epsilon }}}_{\mathrm{12}}$ from F- to T-points. Later on, the term $P/\mathrm{\Delta }$ is also interpolated from T- to F-points in order to estimate σ₁₂. In their implementation of MEB, Plante et al. (2020) also use this approach to interpolate the damage tracer at F-points. However, they report that using the same approach to estimate σ₁₂, and hence σ_II, at T-points when performing the Mohr–Coulomb test (Eq. 13) results in checkerboard instabilities. The solution they propose to prevent the occurrence of these instabilities is to introduce an additional σ₁₂ that is defined at T-points. This additional σ₁₂ is updated at each time step using – as an increment – the average of the four σ₁₂ increments computed at the surrounding F-points.

Note that based on the strong interdependence between the internal stress and the damage in brittle rheologies (Sect. 2.2), as well as the highly localized nature of the damage, we think that the use of the averaging technique to estimate d at the corner points of the C-grid cells should be avoided if possible. Indeed, by using such a technique, ${\mathit{\sigma }}_{\mathrm{12}}^{\left(\mathrm{i}\right)}$, as opposed to ${\mathit{\sigma }}_{\mathrm{11}}^{\left(\mathrm{i}\right)}$ and ${\mathit{\sigma }}_{\mathrm{22}}^{\left(\mathrm{i}\right)}$, is updated using values of E and λ that have been calculated using a poorer estimate of d at the F-point than what the rheology, together with Mohr–Coulomb test, would have produced at this location if explicitly solved on the F-point. This is because the four-point average of a variable such as the damage that is highly heterogeneous in space, even at the grid point scale, cannot provide a very accurate and reliable estimate of the local value.

Finally, with the C-grid, the implementation of the advection of σ₁₂ (F-point) in a way consistent (in terms of the advection scheme used) with that used for σ₁₁ and σ₂₂ (T-point) is somewhat challenging. That is because the advection of a scalar defined at the F-point, using the same scheme as that used for the advection of scalars at T-points, requires the existence of a u and a v at V- and U-points, respectively. These later considerations have prompted us to consider the use of a new spatial discretization approach for the implementation of BBM on the C-grid.

2.3.1 The E-grid approach

To avoid the interpolation of the damage and the stress components between the center and the corner points of the grid cell and allow the consistent advection of all the components of the stress tensor, an additional sea ice velocity vector, denoted as $(\widehat{u},\widehat{v})$, is introduced. As shown in Fig. 3b, the x component of this additional velocity, $\widehat{u}$, is defined at V-points, while its y component, $\widehat{v}$, is defined at U-points. Similarly, the damage tracer is also duplicated, with an additional occurrence at the upper-right corners of the grid cell, i.e., at F-points. This grid staggering arrangement corresponds to that of the Arakawa E-grid (Arakawa and Lamb, 1977; Janjić, 1984; Maier-Reimer et al., 1993).

Figure 3 Transition from (a) the conventional C-grid staggering as used in NEMO to (b) the E-grid staggering proposed in this study. T-centric (red) and F-centric (blue) cells. d is the damage tracer, u and v are the i- and j-wise components of the sea ice velocity vector, and σ_kl represents the components of the internal stress tensor. The $\widehat{x}$ notation indicates that variable x is specific to the F-centric grid. Note: the F-centric counterparts of u_i,j and v_i,j of the T-centric cell are ${\widehat{u}}_{i+\mathrm{1},j}$ and ${\widehat{v}}_{i,j+\mathrm{1}}$.

As suggested by Fig. 3b, the E-grid can be seen as a superposition of two C-grids, in which the cell center of the additional C-grid coincides with the upper-right corner of the original C-grid. For convenience, we will refer to these two grids as F-centric (additional) and T-centric (original), respectively.

In order to minimize the number of modifications and rewriting in the SI3 code, the idea was to restrict the use of this E-augmented C-grid to rheology and advection modules only. The rest of the code, which includes the thermodynamics, remains unmodified and relies entirely on the standard C-grid. As such, only rheology-specific tracers are defined in the E-grid fashion, i.e., at both T- and F-points. In our case, this applies only to the ice damage d and components of the internal stress tensor (even though components of a tensor cannot exactly be considered tracers when it comes to the advection; see Sect. 2.4). However, global tracers, such as ice concentration and thickness, which are updated within the thermodynamics module, remain defined at the T-point only. Consequently, these tracers are interpolated at the F-point within the rheology module whenever needed.

To summarize, in the proposed rheology-specific E-augmented C-grid approach, as shown in Fig. 3, the conventional C-grid model variables are augmented with (i) the u-velocity component at V-points and v-velocity component at U-points; (ii) the ice damage, σ₁₁, and σ₂₂ at F-points; and (iii) σ₁₂ at T-points. This approach implies that all the equations related to the dynamics, including constitutive and momentum equations, as well as the advection, have to be solved twice: on both the T- and F-centric grids. As detailed in Appendix B, the exact same discretization and numerical schemes can be used on both grids, with only the indices of the velocity components on the F-centric grid requiring particular attention: ${\widehat{u}}_{i+\mathrm{1},j}$ and ${\widehat{v}}_{i,j+\mathrm{1}}$ have to be used as the counterparts of u_i,j and v_i,j on the T-centric grid (Fig.3.b). This applies to the computation of the strain rate tensor (Sect. B2.1), the constitutive equation (Sect. B2.2), momentum equation (Sect. B3), the divergence of the stress tensor (Sect. B3.1), and the advection.

At this stage it is important to note that the doubling of the number of computational points implied by the transition to the E-grid in no way relates to an increase in the spatial resolution of the original C-grid. That is because the FD discretization of spatial derivatives on the E-grid (see Appendix B) still relies on the same local spatial increment, i.e., Δx, as that of the original C-grid, regardless of the subgrid considered (T- or F-centric).

2.3.2 The separation of solutions and how it is restrained

With the E-augmented C-grid approach, all rheology-specific prognostic variables are defined at the points where their value is required, and no interpolation is needed to solve the equations. It does, however, result in an apparent overdetermination, which allows the T- and F-centric solutions to evolve somewhat independently from one another. This separation of solutions rapidly degenerates into unrealistically noisy solutions as the spatial consistency of the fields between the two grids deteriorates.

This problem of grid separation has been known since the early adoption of the E-grid by the community (Arakawa, 1972; Mesinger, 1973; Janjić, 1974; Janjić and Mesinger, 1984; Mesinger and Popovic, 2010). It is mostly discussed in the context of the shallow-water equations and is often referred to as “(short) gravity wave decoupling” or “lattice separation”. Various treatments and methods have been proposed, from filtering approaches to more advanced ones such as the introduction of auxiliary velocity points midway between the neighboring tracer points (Mesinger, 1973; Janjić, 1974; Mesinger and Popovic, 2010). Interestingly, the E-grid was used in the Hamburg Large-Scale Geostrophic (LSG) model (Maier-Reimer et al., 1993) in order to achieve more accurate geostrophic balance, while retaining some advantages of the C-grid such as the straightforward discretization of the divergence. In their model, the problem of grid separation, already limited due to the use of a monthly time step, was overcome through adding horizontal viscosity and diffusion. Recently, Konor and Randall (2018) also mentioned the need to introduce a “horizontal mixing process” to avoid the “separation of solutions” when using the E-grid.

The cause of the separation of the two solutions resides in the weak coupling between the two grids, as they only exchange very little information. Specifically, in our case, the only exchange of information between the T- and F-centric grids occurs via the shear stress σ₁₂ and the ice velocity vector. The estimate of ${\mathit{\sigma }}_{\mathrm{12}}^{\left(\mathrm{i}\right)}$ of the T-centric grid (at F-point), based on Eq. (9), uses $\widehat{E}$, $\widehat{\mathit{\lambda }}$, and $\widehat{\stackrel{\mathrm{̃}}{P}}$ of the F-centric grid (at F-point) and conversely for ${\widehat{\mathit{\sigma }}}_{\mathrm{12}}^{\left(\mathrm{i}\right)}$. Similarly, the correction of ${\mathit{\sigma }}_{\mathrm{12}}^{\left(\mathrm{i}\right)}$ in Eq. (12), if occurring, uses ${\widehat{d}}_{\text{crit}}$ and ${\widehat{t}}_d$. For the velocity, the exchange of information occurs in the Coriolis term of Eq. (1) and through the advection of σ₁₂ via $\widehat{u}$ and $\widehat{v}$ (and that of ${\widehat{\mathit{\sigma }}}_{\mathrm{12}}$ via u and v). As suggested by results discussed later in this section, this exchange of information is not sufficient enough to prevent the decoupling of the solutions between the two grids. Hence, a numerical treatment is required to constrain the T- and F-centric solutions to remain spatially consistent with one another.

During the early phase of our development, we considered, implemented, and tested a variety of such treatments. As of now, only one has proven able to prevent the grid separation issue without leading to noisy and/or unrealistic solutions. This treatment, which operates sequentially on the T- and F-centric stress tensors at the dynamical time step level, is hereafter referred to as cross-nudging (CN). It consists of nudging each vertically integrated component of the T-centric stress tensor (σ) towards its F-centric counterpart ($\widehat{\mathit{\sigma }}$) interpolated at the relevant point under even time step integrations and conversely under odd time step integrations. This is written as follows:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}\left(\begin{array}{c}{\mathit{\sigma }}_{\mathrm{11}}\\ {\mathit{\sigma }}_{\mathrm{22}}\\ {\mathit{\sigma }}_{\mathrm{12}}\end{array}\right)=\left(\begin{array}{c}{\mathit{\sigma }}_{\mathrm{11}}\\ {\mathit{\sigma }}_{\mathrm{22}}\\ {\mathit{\sigma }}_{\mathrm{12}}\end{array}\right)& -{\displaystyle \frac{{\mathit{\gamma }}_{\text{cn}}}{N_{\mathrm{s}}}}\left(\begin{array}{l}{\mathit{\sigma }}_{\mathrm{11}}-{\stackrel{\mathrm{‾}}{\widehat{\mathit{\sigma }}}}_{\mathbf{11}}^T/h\\ {\mathit{\sigma }}_{\mathrm{22}}-{\stackrel{\mathrm{‾}}{\widehat{\mathit{\sigma }}}}_{\mathbf{22}}^T/h\\ {\mathit{\sigma }}_{\mathrm{12}}-{\stackrel{\mathrm{‾}}{\widehat{\mathit{\sigma }}}}_{\mathbf{12}}^F/{\stackrel{\mathrm{‾}}{h}}^F\end{array}\right)\\ & \text{ (even\hspace{0.17em}\hspace{0.17em}dynamical\hspace{0.17em}\hspace{0.17em}time step)}\\ \left(\begin{array}{c}{\widehat{\mathit{\sigma }}}_{\mathrm{11}}\\ {\widehat{\mathit{\sigma }}}_{\mathrm{22}}\\ {\widehat{\mathit{\sigma }}}_{\mathrm{12}}\end{array}\right)=\left(\begin{array}{c}{\widehat{\mathit{\sigma }}}_{\mathrm{11}}\\ {\widehat{\mathit{\sigma }}}_{\mathrm{22}}\\ {\widehat{\mathit{\sigma }}}_{\mathrm{12}}\end{array}\right)& -{\displaystyle \frac{{\mathit{\gamma }}_{\text{cn}}}{N_{\mathrm{s}}}}\left(\begin{array}{l}{\widehat{\mathit{\sigma }}}_{\mathrm{11}}-{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathbf{11}}^F/{\stackrel{\mathrm{‾}}{h}}^F\\ {\widehat{\mathit{\sigma }}}_{\mathrm{22}}-{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathbf{22}}^F/{\stackrel{\mathrm{‾}}{h}}^F\\ {\widehat{\mathit{\sigma }}}_{\mathrm{12}}-{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathbf{12}}^T/h\end{array}\right)\\ & \text{ (odd\hspace{0.17em}\hspace{0.17em}dynamical\hspace{0.17em}\hspace{0.17em}time step)},\end{array}\end{array}$$

where γ_cn is the CN coefficient, N_s is the time-splitting parameter (Eq. 15), the bar notation denotes the spatial interpolation from F- to T-points or T- to F-points (see Eq. A1 in Appendix A3), and stress components in bold are vertically integrated (Appendix A4). Each of the two tensors is “corrected” $N_{\mathrm{s}}/\mathrm{2}$ times during the course of one advective time step ΔT. The form of the term that modulates the nudging intensity, i.e., ${\mathit{\gamma }}_{\text{cn}}/N_{\mathrm{s}}$, ensures that the level of cross-nudging undergone by the two tensors under one ΔT is primarily controlled by γ_cn and remains somewhat independent of the choice of N_s.

Due to the strong damage–stress interdependence (Sect. 2.2), the CN is applied to σ⁽ⁱ⁾ rather than σ, i.e., after solving the constitutive equation (Eq. 9) but before computing d_crit and applying the stress correction (Eq. 12). Applying the CN after the stress correction stage (rather than before) may result in the use of a mix of (i) stress values that have been corrected through a local increase in damage and (ii) uncorrected stress values (with no increase in damage). This may lead to spatial inconsistencies between the post-CN stresses and the damage field.

Figure 4 Effect of using different values for the cross-nudging coefficient γ_cn on the simulated sea ice damage. Random snapshot of damage (at T-points, 13 January 1997) after 30 d of simulation using the specified value of γ_cn in a set of sensitivity experiments identical to SI3-BBM: (a) no cross-nudging, (b) γ_cn=0.1, (c) γ_cn=0.5, (d) γ_cn=1 as in SI3-BBM, (e) γ_cn=2, and (f) γ_cn=10.

The four-point spatial averaging (interpolation) used in the CN inevitably results in the introduction of a smoothing of the solution in space. As such, γ_cn is chosen to achieve the best compromise between the smoothing and the coupling of the T- and F-centric solutions. We have performed sensitivity tests with our pan-Arctic setup and we conclude, relying exclusively on the visual assessment of the simulated fields, that the right compromise is achieved when γ_cn is typically of the order of 1, with γ_cn=1 being the value used in our experiments. As illustrated in Fig. 4, with a value below 1, the solutions become increasingly noisy as γ_cn approaches zero. In particular, the damage field tends to exhibit strongly unrealistic straight-line features of high damage that are aligned along the x or y axis of the grid. Our results suggest that values of γ_cn typically above 2 lead to an excessive smoothing of the solutions (as shown, for example, for γ_cn=10 in Fig. 4f). The value of γ_cn appropriate for a given model setup is likely to be dependent on different factors that we have not identified yet. As such, we can only recommend that potential users of our implementation consider γ_cn a tuning parameter that should be adjusted for a given setup. However, simulations that we have performed at spatial resolutions of 1, 2, 4, and 10 km with the idealized test case discussed in Sect. 3.1 (not shown) suggest that γ_cn is only weakly influenced by the spatial resolution at which the model is run (values typically between 0.5 and 2 consistently yielding what we refer to as the best compromise).

2.4 Horizontal advection

In neXtSIM, the Lagrangian finite-element model used by Ólason et al. (2022), the advection occurs implicitly at each advective time step (also corresponding to the thermodynamics time step) through the ice-velocity-driven displacement of the mesh elements. As such, the rate of change of a prognostic scalar ϕ is $\dot{\mathit{\varphi }}\equiv {\partial }_{\mathrm{t}}\mathit{\varphi }$. In the present Eulerian context, however, the term relative to the horizontal advection has to be considered so that the rate of change of ϕ is now ${\partial }_{\mathrm{t}}\mathit{\varphi }+U{\partial }_x\mathit{\varphi }+V{\partial }_y\mathit{\varphi }$. In our implementation, as pointed out by Ólason et al. (2022), this advection term is computed and added to the trend of the prognostic scalar considered every advective time step. Thus, the sea ice velocity vectors U and V that we consider for the advection at the advective time step level represent the mean of the N_s successive velocity vectors (u and v) calculated under one time-splitting instance.

We use the second-order moment (SOM) advection scheme of Prather (1986) available in SI3 to advect the damage and the components of the stress tensors (considered to be scalars for now; see Sect. 2.4.1). Technically, the damage and stress tensor components defined at the T-point (d, σ₁₁, σ₂₂, and ${\widehat{\mathit{\sigma }}}_{\mathrm{12}}$) are advected using U and V defined at U- and V-points, respectively. Their F-point counterparts ($\widehat{d}$, ${\widehat{\mathit{\sigma }}}_{\mathrm{11}}$, ${\widehat{\mathit{\sigma }}}_{\mathrm{22}}$, and σ₁₂) are advected using $\widehat{U}$ and $\widehat{V}$ defined at V- and U-points, respectively. In practice, the exact same implementation of the advection scheme can be used to perform the advection at T- and F-points; the only difference is that for the advection of F-point scalars, the spatial indexing of the velocity components is staggered by one cell. Namely, ${\widehat{U}}_{i+\mathrm{1},j}$ and ${\widehat{V}}_{i,j+\mathrm{1}}$ have to be used in place of U_i,j and V_i,j (Fig. 3b).

As it is commonly done in sea ice models and justified by a scale analysis of the momentum equation, the term for the advection of momentum is neglected.

2.4.1 Advection of the internal stress tensor

In the Eulerian framework, the rate of change of a second-rank tensor must introduce additional terms to the material time derivative in order for the dynamics of the tensor to remain independent of the frame of reference (Oldroyd, 1950; Larson, 1988; Hinch and Harlen, 2021; Stone et al., 2023). These terms account for the effects of rotation and deformation of the medium on the evolution of the stress tensor and are gathered here in a symmetric tensor $\dot{\mathit{L}}$:

$$\begin{array}{}\text{(17)}& \dot{\mathit{\sigma }}\equiv {\displaystyle \frac{D\mathit{\sigma }}{Dt}}+\dot{\mathit{L}}\equiv {\partial }_{\mathrm{t}}\mathit{\sigma }+(\mathit{U}\cdot \mathbf{\nabla })\mathit{\sigma }+\dot{\mathit{L}}.\end{array}$$

As stressed by Snoeijer et al. (2020), one faces a “a somewhat unpleasant ambiguity” as two different formulations exist for $\dot{\mathit{L}}$. Both formulations are equally valid in terms of frame invariance, and so is any linear combination of the two. The first formulation is the upper-convected time derivative of σ, denoted as $\stackrel{▿}{\mathit{\sigma }}$,

$$\begin{array}{}\text{(18)}& \stackrel{▿}{\mathit{L}}=-(\mathbf{\nabla }\otimes \mathit{U}{)}^{⊺}\mathit{\sigma }-\mathit{\sigma }(\mathbf{\nabla }\otimes \mathit{U}),\end{array}$$

which, in component form, simplifies into

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\stackrel{▿}{L}}_{\mathrm{11}}& =-\mathrm{2}\left({\dot{\mathit{\epsilon }}}_{\mathrm{11}}{\mathit{\sigma }}_{\mathrm{11}}+{\partial }_{y}U{\mathit{\sigma }}_{\mathrm{12}}\right),\\ {\stackrel{▿}{L}}_{\mathrm{22}}& =-\mathrm{2}\left({\dot{\mathit{\epsilon }}}_{\mathrm{22}}{\mathit{\sigma }}_{\mathrm{22}}+{\partial }_{x}V{\mathit{\sigma }}_{\mathrm{12}}\right),\\ {\stackrel{▿}{L}}_{\mathrm{12}}& =-\left({\dot{\mathit{\epsilon }}}_{\mathrm{11}}+{\dot{\mathit{\epsilon }}}_{\mathrm{22}}\right){\mathit{\sigma }}_{\mathrm{12}}-{\partial }_{x}V{\mathit{\sigma }}_{\mathrm{11}}-{\partial }_{y}U{\mathit{\sigma }}_{\mathrm{22}}.\end{array}\end{array}$$

The second formulation is the lower-convected time derivative, denoted as $\stackrel{△}{\mathit{\sigma }}$,

$$\begin{array}{}\text{(20)}& \stackrel{△}{\mathit{L}}=\mathit{\sigma }(\mathbf{\nabla }\otimes \mathit{U}{)}^{⊺}+(\mathbf{\nabla }\otimes \mathit{U})\mathit{\sigma },\end{array}$$

with

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\stackrel{△}{L}}_{\mathrm{11}}& =\mathrm{2}\left({\dot{\mathit{\epsilon }}}_{\mathrm{11}}{\mathit{\sigma }}_{\mathrm{11}}+{\partial }_{x}V{\mathit{\sigma }}_{\mathrm{12}}\right),\\ {\stackrel{△}{L}}_{\mathrm{22}}& =\mathrm{2}\left({\dot{\mathit{\epsilon }}}_{\mathrm{22}}{\mathit{\sigma }}_{\mathrm{22}}+{\partial }_{y}U{\mathit{\sigma }}_{\mathrm{12}}\right),\\ {\stackrel{△}{L}}_{\mathrm{12}}& =\left({\dot{\mathit{\epsilon }}}_{\mathrm{11}}+{\dot{\mathit{\epsilon }}}_{\mathrm{22}}\right){\mathit{\sigma }}_{\mathrm{12}}+{\partial }_{y}U{\mathit{\sigma }}_{\mathrm{11}}+{\partial }_{x}V{\mathit{\sigma }}_{\mathrm{22}}.\end{array}\end{array}$$

These formulations of $\dot{L}$ are straightforward to implement in the model as they only involve multiplications between the components of tensors $\dot{\mathit{\epsilon }}$ and σ, which are all defined at both T- and F-points with the E-grid. Therefore, we have implemented both formulations in SI3. In our implementation, the standard advection trend for each tensor component, corresponding to the term (U⋅∇) σ in Eq. (17), is computed using the identical scheme to that used for regular scalar fields. The tensor-specific advection trend, $\dot{\mathit{L}}$, is computed according to Eq. (19) or (21). These two contributions are computed independently from one another using stress values that have not been updated yet by the advection process.

We have chosen to use the upper-convected formulation in both the idealized and pan-Arctic simulations presented in this paper. This choice is purely arbitrary and is not based on scientific considerations of any kind. It relies solely on the fact that the upper-convected formulation has been favored in the literature since Oldroyd, who introduced both formulations in his 1950 paper, found that his model would only realistically represent the flow around a rotating rod when using this formulation, as opposed to the lower-convected one (Hinch and Harlen, 2021). Nevertheless, two twin simulations of our reference BBM simulation (see Sect. 3.2) have been run, one using the traditional material derivative (i.e., $\dot{\mathit{L}}=\mathrm{0}$) and the second the lower-convected formulation. All the diagnostics and deformation statistics discussed later in this paper have been performed for these two additional simulations and no significant differences have been identified among the three options (PDFs of the total deformation for the reference and additional simulations are provided in Fig. C3 in Appendix C as an example).

Further work, involving, for example, the design of new idealized test cases, should be conducted to address this ambivalence and help identify which time-derivative formulation (or combination of them, such as the Gordon–Schowalter time-derivative discussed by Dansereau et al., 2016) is best adapted to sea ice rheology.

2.5 Construction of observed and simulated Lagrangian sea ice deformations

Our assessment of the NEMO pan-Arctic simulations relies on a multiscale statistical analysis of sea ice deformation rates constructed using observed and simulated Lagrangian sea ice trajectories during winter 1996–1997. Observed trajectories are taken from the RGPS (RADARSAT Geophysical Processor System Lagrangian trajectories) dataset of Kwok et al. (1998), while simulated trajectories are generated from the Eulerian sea ice velocities of SI3 by means of a sea ice particle tracker program.

The preprocessing and computing approach we use to construct sea ice deformations out of the raw RGPS Lagrangian trajectories is very similar to that used by Ólason et al. (2022), the main difference being that it relies on the tracking of quadrangles rather than triangles. To construct the SI3-derived synthetic version of these deformations, the tracking software seeds the identical points to those involved in the definition of the quadrangles selected for computing the RGPS deformation, respecting their initial position in space and time. These points are then tracked for about 3 d using the hourly averaged Eulerian sea ice velocities of SI3; the exact tracking duration used is that of the time interval between the two consecutive positions of the corresponding RGPS point.

Our period of interest, spanning 15 December 1996 to 20 April 1997, is divided into 3 d bins, which correspond to the nominal time resolution of the RGPS dataset.

2.5.1 Selection of RGPS trajectories

For each 3 d bin, an initial subset of the RGPS points is selected. Each point of this initial subset must satisfy the following requirements:

• The point has at least one position that occurs within the time interval of the bin; this position, or the earliest-occurring one if there is more than one occurrence, is selected and referred to as position no. 1.

• Position no. 1 is located at least 100 km away from the nearest coastline.

• The point has at least one upcoming position that occurs 3 d after position no. 1, with a tolerated deviation of ± 6 h, referred to as position no. 2 (in the event of more than one position satisfying this requirement, the position yielding the time interval the closest to 3 d is selected).

2.5.2 Quadrangulation of selected trajectories

A Delaunay triangulation is performed on this initial subset of points at position no. 1. Triangles whose areas differ by more than 25 % from the half of the nominal area of the quadrangles to be constructed (i.e., the square of the spatial scale under consideration) or with an angle below 5° or above 160° are excluded. Neighboring pairs of remaining triangles are then merged into quadrangles in order to transform the triangular mesh into a quadrangular mesh.

Aspiring quadrangles at position no. 2 are constructed by simply considering the exact same respective sets of four points as those defining quadrangles at position no. 1.

Finally, only points that define quadrangles that satisfy the following requirements at both position no. 1 and position no. 2 are retained:

• The square root of the area of the quadrangle falls within a ± 12.5 % range of agreement with the horizontal scale under consideration.

• The time position of each of the four points defining the vertices of the quadrangle should not differ from that of any of the other three points by more than 60 s.

• The thresholds for the minimum and maximum angles allowed are 40 and 140°, respectively.

2.5.3 Computation of deformation rates based on the quadrangles

For all quadrangles selected in a given 3 d bin, strain rates are computed based on their position no. 1 and no. 2 using the line-integral approximations (see, e.g., Lindsay and Stern, 2003, Eqs. 10–14).

The time location (date) assigned to a given deformation rate corresponds to the center of the time interval defined by position no. 1 and position no. 2 of each quadrangle. The spatial location of the deformation rates corresponds to the barycenter of the four vertices of the quadrangle considered at the center of this same time interval.

2.5.4 Construction of the simulated Lagrangian sea ice trajectories

To save computer resources, only the points from which valid RGPS deformation estimates were computed are retained. Each of these points is seeded using the same initialization date and location (bilinear interpolation) as its RGPS counterpart. It is then tracked during the same time interval of about 3 d (± 6 h) that separates the two consecutive records of the RGPS point considered. The tracking software uses a time step of 1 h and feeds on the hourly averaged simulated sea ice velocities of the SI3 experiments. Note that only the conventional C-grid velocities u and v of the T-centric cell are used to track the points ($\widehat{u}$ and $\widehat{v}$, available in the BBM-driven simulation, are not used).

3 Model evaluation

We use version 4.2.2 of the NEMO modeling system (Madec et al., 2022) as the basis for the development of the BBM rheology code extension and to perform both the idealized and coupled pan-Arctic simulations to be assessed. Since version 4, the default sea ice component of NEMO has been SI3 (Vancoppenolle et al., 2023).

3.1 Idealized simulations

To provide a first qualitative evaluation of the behavior of our BBM implementation, SI3 simulations were run on the idealized test case setup introduced by Mehlmann et al. (2021), including simulations using the default aEVP-driven SI3 setup for reference purposes. This test case, defined on a 512 km wide square domain, simulates a cyclone traveling in the northeastward direction over a thin layer of ice (h ≃ 0.3 m) that floats on an anticyclonically circulating ocean. This test case is well suited to illustrate the influence of the grid discretization on rheology-related processes such as the representation of linear kinematic features (LKFs) (Danilov et al., 2022, 2024), which makes it particularly relevant to our study. We use the identical setup and parameter values to those defined in Mehlmann et al. (2021) (see the “Code and data availability” section to access the SI3 namelists and forcing files). SI3 is run in stand-alone mode using SAS, the stand-alone surface module of NEMO. In SAS mode, SI3 uses a prescribed surface ocean state (current, height, temperature, and salinity) instead of being coupled to the ocean component of NEMO as in our pan-Arctic simulations (Sect. 3.2).

Figure 5 Sea ice total deformation (instantaneous) in the test case described by Mehlmann et al. (2021) after 48 h of simulation with SI3 using the default SI3 aEVP setup and the newly implemented BBM rheology (left-hand and middle column, respectively) and the corresponding PDFs (right-hand column) for simulations run at a spatial resolution of 2, 4, and 10 km (first, second, and third row, respectively).

The results of this test case, for both the aEVP and BBM rheology, are shown in Fig. 5. First we note that for the SI3 implementation of aEVP, the deformation fields obtained are in qualitative agreement with those of Mehlmann et al. (2021) (see, for instance, their Fig. 7). In the solutions obtained with BBM, we note the presence of a circular network of LKFs that contrast, by their arrangement, with the “spider-web-like” arrangement of the LKFs in the aEVP solution. In the BBM-driven simulation these LKFs are also simulated in the 10 km setup (Fig. 5h). The spatial pattern of the LKFs, particularly those accommodating the highest deformation, also look qualitatively different: apparently longer and with circular and concentric shapes with respect to the center of the forcing cyclone in the case of BBM, shorter and in radial alignment with respect to the forcing cyclone in the aEVP case.

We also find that the background deformation field is close to zero in the BBM solution, except along the LKFs, whereas the deformation looks more homogeneous in space in the aEVP solution. This can also be seen on the respective PDFs (Fig. 5c, f, and i) that exhibit different shapes and heavier tails in the BBM solution. We have verified that the aEVP solution is not too significantly impacted by the number of iterations used in the aEVP solver of SI3 by conducting the same aEVP experiments with a N_EVP = 1000 instead of N_EVP = 100 (Fig. C1).

Figure 6 Geographical extent, numerical grid, and actual local spatial resolution of the NANUK4 computational domain that is used in the experiments. For ease of visual representation of the grid cells, grid points have been subsampled by a factor of 4.

Finally, we note that the solutions do not contain any apparent numerical instabilities or noise for aEVP or BBM. The LKF-like features in the BBM solution at 10 km show a tendency to align horizontally, vertically, and diagonally with the grid. As of now, we are unable to provide an explanation of the mechanism responsible for these alignments, apart from their apparent connection with the use of a relatively coarse spatial resolution, as the solution obtained with the 4 km setup seems to be rid of them.

3.2 Coupled ocean–sea ice pan-Arctic simulations

The pan-Arctic simulations use SI3 coupled to the 3D ocean component of NEMO, named OCE. They are performed on the so-called NANUK4 regional configuration, which is an Arctic extraction of the standard global $\mathrm{1}/\mathrm{4}$° resolution NEMO gridded horizontal domain known as ORCA025 (Barnier et al., 2006). As such, and as shown in Fig. 6, the actual grid resolution of NANUK4 typically spans 10 up to 14 km in the central Arctic region. NANUK4 features two open lateral boundaries; the southernmost boundary is located at about 39° N in the Atlantic Ocean, and the second boundary is located south of the Bering Strait, at about 62° N in the Pacific Ocean. The vertical z-coordinate grid used for the ocean features 31 levels with a Δz of 10 m at the surface up to about 500 m at the deepest level, at a depth of 5250 m.

The hindcast nature of the simulations is achieved through the use of interannual atmospheric and oceanic forcings at the surface and at the two open boundaries, respectively. For the atmospheric forcing, both the ocean and the sea ice components receive, as surface boundary conditions, fluxes of momentum, heat, and fresh water at the air–sea and air–ice interface, respectively. These fluxes are computed every hour by means of bulk formulae using the hourly near-surface atmospheric state from the ERA5 reanalysis of the ECMWF (Hersbach et al., 2020) and the prognostic surface temperature of the relevant component (SST or ice surface temperature).

For the lateral boundary conditions of OCE, the 3D ocean is relaxed towards the monthly averaged 3D horizontal velocities, temperature, salinity, and SSH (2D) of the GLORYS2v4 ocean reanalysis (Ferry et al., 2012).

Table 1 Values of default aEVP-related parameters in SI3, as used in experiment SI3-default.

Table 2 Values of BBM-related parameters implemented in SI3, as used in experiment SI3-BBM.

Both OCE and SI3 use a time step of ΔT = 720 s, the advective time step. The coupling between these two components is also done at each advective time step. Our control simulation, named SI3-default, uses the default SI3 setup as provided in NEMO and thereby uses the aEVP rheology of Kimmritz et al. (2016). The second simulation, named SI3-BBM, only differs from SI3-default through the use of our implementation of the BBM rheology in place of aEVP and a higher value of the air–ice drag coefficient. Values of parameters relevant to both rheologies used in the two simulations are provided in Tables 1 and 2, respectively.

The two simulations are initialized on 1 December 1996 using the same restart data generated at the end of a 2-month spin-up performed with the SI3-default setup and run until 20 April 1997. This spin-up is initialized on 1 October 1996 by using the daily averaged ocean and sea ice data from the GLORYS2v4 reanalysis as an initial condition. More specifically, OCE is initialized at rest (no current) with the 3D temperature and salinity state of the reanalysis, while SI3 is initialized with the sea ice concentration and thickness. The 2-month spin-up we use is long enough to get the ocean velocities in the upper ocean into a good state with the given temperature and salinity fields. Our strategy implies that SI3-BBM is initialized on 1 December 1996, with a value of ice damage set to zero everywhere, which poses no issue as the time required to spin up the damage is very short (Bouillon and Rampal, 2015; Rampal et al., 2016), typically of the order of a few days. The analysis of the results is performed for the period 15 December 1996–20 April 1997, leaving the upper ocean and the sea ice cover in SI3-BBM 2 weeks to respond to the changed rheology, which should be ample.

For these simulations, the adjustable tuning parameters of SI3 are kept as close as possible to those of the reference configuration of NEMO (Tables 1 and 2). As such, the thermodynamic component uses five ice categories. Yet, the ice–atmosphere drag coefficient $C_{\mathrm{D}}^{\left(\mathrm{a}\right)}$ has been adjusted from 1.4 × 10⁻³ to 2 × 10⁻³ in SI3-BBM in order for the mean simulated deformation rate at the 10 km scale to be in agreement with that derived from the satellite observations against which we evaluate the model (see Sects. 2.5 and 3.3.3). In SI3-default, the default value of 1.4 × 10⁻³ satisfies this requirement and is left unchanged. For the time-splitting approach (Sect. 2.2), we use a dynamical time step of 7.2 s in SI3-BBM, which relates to a time-splitting parameter of N_s=100.

Figure 7 Maps of the sea ice total deformation rate at the 10 km spatial and 3 d temporal scale for the period centered about 24 December 1996, computed based on (a) RGPS Lagrangian data and (b, c) their synthetic counterparts constructed using the simulated sea ice velocities of SI3-BBM and SI3-default, respectively. Empty regions correspond to the absence of satellite data during the period concerned.

3.3 Comparison of simulated sea ice deformation statistics against satellite data

3.3.1 Probability density function of sea ice deformation rates

As illustrated by the maps of the 3 d total deformation rates (Fig. 7), RGPS clearly exhibits narrow and long features, commonly called linear kinematic features (LKFs; Kwok, 2001) along which the deformation is concentrated. Visually, LKFs simulated by SI3-BBM appear somewhat realistic in terms of both length and orientation, and the magnitude of the deformation rates along these LKFs is similar to that of RGPS. We note that SI3-default exhibits very smooth fields of deformation with almost no such localized features; this is consistent with the findings of recent studies that evaluate VP-driven sea ice simulations run at a grid resolution typically coarser than 5 km (e.g., Ólason et al., 2022; Bouchat et al., 2022).

Figure 8 PDFs of the (a) shear, (b) divergence, (c) convergence, and (d) total deformation rates at the 10 km spatial and 3 d temporal scale for RGPS data and their synthetic counterparts constructed using the simulated sea ice velocities of SI3-BBM and SI3-default. The light gray lines are for reference and correspond to a power law with an exponent of −3. Below each panel, the departure between the logarithm of the simulated and observed distributions is shown for each bin.

The PDFs of the total deformation rates (Fig. 8d) show that SI3-BBM exhibits a heavy tail similar to that of RGPS and that it can be approximated by a power law over the values corresponding to the last two percentiles of the RGPS distribution, although with slightly different exponents (−2.9 and −3, respectively). A look at the other invariants of the deformation (i.e., shear, divergence, and convergence rates in Fig. 8a–c) shows that SI3-BBM simulates large deformation events as seen in the observations, which suggests the ability of BBM to capture the heterogeneous character of sea ice deformation in this setup. In contrast, SI3-BBM is unable to reproduce the observed convergence over the full range of values present in the RGPS data (Fig. 8c). This limitation of the BBM rheology is further discussed in Sect. 4.2.

Our results suggest a propensity for SI3-default to underestimate the extreme values of deformation rates. This inadequacy of the model could very likely be mitigated by conducting a finer tuning of the parameters related to the viscous–plastic rheology, in particular through the better adjustment of the ratio between the ice compressive strength and the ice shear strength (Bouchat and Tremblay, 2017). Yet, conducting such a tuning is out of the scope of this paper. Interestingly, the underestimation of extreme deformation values set aside, SI3-default exhibits a power-law behavior similar to that of both observations and SI3-BBM, with similar exponents, except in convergence.

3.3.2 Time series of sea ice deformation rates

Following Ólason et al. (2022), we examine the 90th percentile of the total deformation (P90), chosen for its sensitivity to the high values that contribute to shaping the long tail of the PDFs of deformations. Technically, P90 is the value of deformation below which 90 % of deformation values in the frequency distribution fall. P90 is computed from each snapshot of deformation from mid-December 1996 to late April 1997 to evaluate the temporal evolution of the deformation. In addition to the 90th percentile, we also consider the 95th and 98th percentiles.

Figure 9 Time series of the (a) 90th, (b) 95th, and (c) 98th percentiles of the sea ice total deformation rate for winter 1996–1997 at the 10 km spatial and 3 d temporal scale for RGPS data and their synthetic counterparts constructed using the simulated sea ice velocities of SI3-BBM and SI3-default.

As illustrated in Fig. 9, SI3-BBM is in fairly good agreement with the observations, in particular for the P90 values (see Table 3). We note, however, that despite the ability of SI3-BBM to reproduce a variability similar to that observed, the higher the percentile value, the lower the agreement between the magnitudes. This suggests an inability of the BBM rheology to capture the most extreme deformation events.

Table 3 Bias, RMSE, and Pearson correlation of the deformation rate time series in Fig. 9 obtained between each simulation and RGPS.

We note that the biases and RMSEs are very similar between the two simulations. For P90 and P95, the values suggest fairly good agreement between the two simulations and the observations. The values for P98, however, highlight the inability of both models to reproduce extreme deformation events. This is in qualitative agreement with what Ólason et al. (2022) already reported. Yet, further investigation remains necessary to assess whether this is inherent to the BBM model or could be improved through the better adjustment of the rheological parameters.

3.3.3 Multi-fractal scaling analysis

The presence of heavy tails in the distributions shown in Fig 8 implies that one needs to consider higher moments than the mean to fully describe the statistics of the sea ice deformation (Sornette, 2006). Following Marsan et al. (2004), our multi-fractal scaling analysis should be restricted to the consideration of moments of order q > 0 because zero values exist in the deformation field. Moments of order q> 3 should also be disregarded because a transition occurs typically between q_c =  2.5 and q_c = 3 (Schertzer and Lovejoy, 1987). This is a consequence of the tails of the distributions for RGPS and SI3-BBM being well represented by a power law with an exponent of about −3, which would cause moments of order q>q_c to diverge (Savage, 1954).

We performed a multi-fractal spatial scaling analysis of the RGPS total deformation rates and their simulated counterparts, considering the moments q=1, 2, and 3 of the distributions. Both the observed and simulated statistics (mean, variance, and skewness) follow power laws (Fig. 10). In particular, the observed mean sea ice deformation rate $\langle \dot{\mathit{\epsilon }}\rangle$ is particularly well reproduced in SI3-BBM across the full range of spatial scales considered for this analysis and can be approximated by a power-law scaling, i.e., $\langle \dot{\mathit{\epsilon }}\rangle \sim L^{-\mathit{\beta }}$, where L is the spatial scale and β an exponent of about 0.15. We note that the atmospheric drag coefficient was used as the adjustment parameter in SI3 (Sect. 3.2), which led to the use of $C_{\mathrm{D}}^{\left(\mathrm{a}\right)}$ = 2 × 10⁻³ in SI3-BBM, while the $C_{\mathrm{D}}^{\left(\mathrm{a}\right)}$ = 1.4 × 10⁻³ used in SI3-default did not require any adjustment. Consistent with the results previously discussed, the higher moments, which characterize the largest and most extreme values of the distributions, remain underestimated in SI3-BBM compared to that derived from the observations. Indeed, the exponents of the power law that fits the SI3-BBM data (β = −0.6 and −1.34 for q = 2 and 3, respectively) are lower than those derived from RGPS data (β = −0.7 and −1.52). This indicates that SI3-BBM is not fully capturing the strength of the spatial scaling of sea ice deformation revealed by the observations or, in other words, that it fails to achieve the extremely high degree of spatial localization of the LKFs in the observations. Figure 10 suggests that the total deformation rates simulated by SI3-default cease to follow the expected power law for scales larger than 100 km. This is in line with published results (e.g., Hutter et al., 2018; Bouchat et al., 2022). Hutter et al. (2018) argue that the VP model needs approximately 10 grid cells to be able to resolve features, which suggests that the “effective resolution” of the model is 10 times coarser than that of the numerical grid on which it is run. This implies that one should instead consider fitting the deformation rates at a resolution 10 times coarser than that used by the model, i.e., 130 km in our case. This would yield power-law slopes that are in better agreement with those derived from observations. We argue that since sea ice deformation is a scale-invariant process at the geophysical scale, a sea ice model should be able to represent this scaling down to the model grid cell. Figure 10 suggests that our BBM implementation allows SI3 to achieve this despite the use of the Eulerian framework.

Figure 10 Spatial scaling analysis of the observed and simulated total deformation rate calculated over a 3 d timescale (all based on the motion of the same RGPS quadrangles) based on RGPS data and their synthetic counterparts constructed using the simulated sea ice velocities of SI3-BBM and SI3-default. Moments of order $q=\mathrm{1},\mathrm{2},\mathrm{3}$ of the distributions of the total deformation rate were calculated at scales spanning from 10 up to 640 km. The solid straight lines indicate the associated power-law scaling based on the least-square fit using values from 10 to 160 km. Values for 320 and 640 km are excluded due to excessive uncertainty resulting from the small sample size. Note: we used logarithmically spaced bins and applied an ordinary least-square method to the binned data in log–log space to get a reasonably accurate estimate of these power-law fits (Stern et al., 2018).

The simulated and observed structure functions (i.e., the dependence of the scaling exponents of the power law on the order of the moment) β(q) are shown in Fig. 11. The spatial scaling obtained from both the observations and SI3-BBM are multi-fractal because their structure functions is well approximated (in the sense of the least-square method) by a quadratic function of the type $\mathit{\beta }\left(q\right)=a{q}^{\mathrm{2}}+bq$. One should note that in the universal multi-fractal formalism, the structure functions are not required to be quadratic and can have a varying degree of nonlinearity (Lovejoy and Schertzer, 2007). A quadratic structure function, as obtained here, simply means that the process of sea ice deformation can be approximated by a lognormal multiplicative cascade model with a maximum degree of multi-fractality. The structure function of SI3-BBM shows a curvature a that has a magnitude comparable to that of RGPS, i.e., 0.15 versus 0.17. These values of curvature are in fair agreement with those obtained from Lagrangian simulations performed with neXtSIM and reported in previous studies: 0.14 in Rampal et al. (2016) and 0.11 in Rampal et al. (2019).

Figure 11 Structure functions β(q) for the RGPS data (black), SI3-BBM (blue), and SI3-default (red), where β indicates the exponent of the power-law fits indicated in Fig. 10 and q is the moment order.

4 Discussion

4.1 On the numerical implementation

The cross-nudging has a noteworthy analogy to the Asselin filter (Asselin, 1972) used when discretizing time derivatives of a prognostic variable by means of the leap-frog scheme (three time levels, centered, and second-order), in particular in the context of shallow-water equations. The goal of this Asselin filter is to subtly average the solutions of neighboring time levels to prevent the separation of trajectories between the even and odd time step levels (Marsaleix et al., 2012). As such, the cross-nudging can be seen as a sort of spatial and two-dimensional analog to the Asselin filter. Despite the crudeness of this approach, which tends to be problematic due to the unavoidable loss of conservation properties, the Asselin filter is still largely used in modern CMIP-class OGCMs like NEMO. Indeed, the ocean component of NEMO used in the simulations presented in this study still relies on it. As of now, our cross-nudging approach clearly lacks physical and numerical consistency, but it somehow allows demonstrating that the implementation of a brittle rheology, along with the advection of the internal stress tensor, is feasible onto an E-augmented C-grid, provided a method to prevent the separation of solutions is used. Nevertheless, we plan to further investigate the possibility of implementing approaches that are more physically and numerically consistent. For instance, an option is to apply the cross-nudging to the two invariants of the stress tensor (i.e., σ_I and σ_II) and the rate of internal work of the ice. This would introduce three equations for three invariant quantities, from which the three components of the stress tensor could be deduced afterward. Another option is to explore the possibility of deriving a numerical formulation inspired from that of Mesinger (1973) and Janjić (1974), in which auxiliary velocity (or stress) points are introduced midway between the neighboring tracer (or velocity) points.

Another critical requirement, this time stemming from the use of the Eulerian and finite-difference framework, has to do with the ability of the advection scheme to advect fields with as little numerical diffusion or dispersion as possible. This is particularly critical when using a brittle rheology like BBM, as most fields exhibit sharp gradients, often associated with linear kinematic features. We chose to use the scheme of Prather (1986), the dispersive scheme option of SI3, to favor the conservation of sharp gradients at the cost of potential noise and overshoots reminiscent of the Gibbs phenomenon. One could, however, consider the use of a different approach that would optimize the advection of sharp gradients, for instance a spatial discretization based on the discontinuous Galerkin method. This method has proven to be efficient and accurate in treating the advection of sea ice variables in the case of a brittle sea ice rheology such as MEB (Dansereau et al., 2017) but has not yet been tested in the context of large-scale, long-term sea ice simulations. This is the scope of our present work and future papers.

Table 4 Computational cost of 3 months (90 d) of pan-Arctic sea ice simulation at $\mathrm{1}/\mathrm{4}$° resolution with SI3 on the NANUK4 regional domain (31 vertical levels), with an advective time step of 720 s, run on 29 cores in parallel (with output data writing disabled to limit the influence of I/O). Default SI3 aEVP setup (N_EVP=100) versus newly implemented BBM rheology (N_s=100) for both stand-alone (SI3-SAS) and coupled (SI3-OCE) modes.

As discussed in Sect. 2.3.1, the use of the E-grid in the dynamics and advection modules of SI3 implies that equations specific to the momentum and the constitutive law are solved on both the T- and the F-centric grids. Moreover, with the need to advect the stress tensor and the damage tracer, specific to brittle rheologies, 2 × 4 additional scalar fields need to be advected. This inevitably leads to an increase in the computational cost of SI3. We have estimated this extra cost by comparing the wall-time length required to complete a 90 d simulation with each rheology using the same 29 cores in parallel on the same computer. Our results, summarized in Table 4, suggest that the increase in the computational cost associated with the use of BBM in place of aEVP is about 45 % when SI3 is used in a stand-alone mode (SAS). In standard coupled mode, with SI3 coupled to OCE, the BBM-related cost increase is about 20 %. This lower value is explained by the fact that by default, the coupling between OCE and SI3 is done sequentially. As such, the cost of SI3 simply adds up to that of OCE, and the cost of OCE is expected to be independent of the mode used (in our case: 113 and 114 CPU hours for SI3-default and SI3-BBM, respectively). Based on our results, the relative cost of SI3 in coupled mode is about 40 % when using the default aEVP setup and about 50 % with our BBM implementation.

4.2 On the simulated sea ice deformations

Based on comparisons against various types of observations, recent studies suggest that large-scale models using BBM can realistically simulate the dynamics and properties of sea ice (Ólason et al., 2022; Rheinlænder et al., 2022; Boutin et al., 2023; Regan et al., 2023). Yet, the deformation in convergence and the subgrid-scale processes related to sea ice ridging are not represented by BBM with the same degree of accuracy. The model overestimates the number of converging events with magnitudes of about 1 % to 5 % per day and underestimates the most extreme events (Fig. 8c and Ólason et al., 2022). As of now, parameter tuning, in particular that of the BBM-specific ridging threshold parameter P_max, has not helped to improve the agreement with the observed PDFs of convergence (not shown). Therefore, we conclude that some fundamental processes need to be reconsidered in BBM.

In Sect. 3.3.3, we find that the degree of multi-fractality of the deformation fields simulated by SI3-BBM is slightly lower than that obtained from the RGPS data. The fact that the deformation fields simulated by neXtSIM in Ólason et al. (2022) are in better agreement with RGPS in this regard suggests that this problem is linked to some numerical aspects of our BBM implementation rather than the BBM rheology itself. This is most likely the consequence of the introduction of some additional numerical dispersion and diffusion by the advection scheme and the cross-nudging treatment, respectively, as these two features are absent in neXtSIM. Moments of order 2 and 3 are expected to be more affected than the mean by an unwanted source of noise and diffusion, which might explain why SI3-BBM reproduces the mean across all scales remarkably well and why the power-law exponents for the variance and the skewness are underestimated. In this regard, the use of the finite-element method together with the discontinuous Galerkin method might prove to be a promising combination to simulate the multi-fractality of sea ice deformation even more accurately while remaining in the Eulerian and quadrilateral mesh framework.

5 Conclusions

The brittle Bingham–Maxwell rheology, known as BBM, has been successfully implemented into SI3, the CMIP-class, Eulerian finite-difference sea ice model of the NEMO modeling system. We have shown that our implementation, which features a prognostic ice damage tracer and a prognostic internal stress tensor, is able to realistically simulate sea ice deformation statistics on a pan-Arctic scale when compared to satellite observations.

Our implementation uses a new discretization approach that expands the C-grid of NEMO into the Arakawa E-grid in the parts of the code dedicated to sea ice dynamics. We have chosen to do so in order to (i) avoid resorting to the spatial averaging of prognostic fields, in particular the damage tracer, as an interpolation technique between the center and corner points of the grid cells, and (ii) allow the straightforward advection of the shear component of the stress tensor. However, by solving the dynamics on the E-grid, the issue of the grid separation is introduced. We have introduced a simple technique to prevent this grid separation in the form of the cross-nudging. This cross-nudging relies on the averaging of the components of the stress tensor and, as such, introduces a spatial smoothing of these components. Despite the fact that this aspect of our implementation is in contradiction to one of our initial motivations (i.e., avoid the use of spatial averaging), we think that our E-augmented C-grid approach is promising. This is because the damage tracer is never averaged, which we think is beneficial for the consistency of the brittle model, and the advection of the shear component of the stress tensor is straightforward and numerically consistent with that of the trace components.

For the advection of the stress tensor, we have chosen to use the upper-convected time derivative rather than its lower-convected counterpart, a combination of the two, or simply the standard material derivative. This choice, based on arbitrary considerations, has no significant impact on the deformation statistics presented in this paper. Both formulations are available in our implementation, which will allow SI3 users to further investigate on this matter, in particular by means of dedicated idealized test cases.

We carried out a statistical analysis of the sea ice deformation rates obtained from a set of realistic pan-Arctic coupled ocean–sea ice simulations for winter 1996–1997, performed with SI3 at a horizontal resolution of about 12 km. Based on a comparison with satellite observations, this analysis demonstrates that the use of the newly implemented BBM rheology results in simulated sea ice deformation statistics that are realistic. In particular, we show that the use of BBM allows simulating highly localized (nearly linear) kinematic features within the sea ice cover, along which the most substantial deformation rates are concentrated.

The observed non-Gaussian statistics of the sea ice deformation process are clearly present in the simulation that uses our newly implemented BBM rheology, except the most extreme values and more particularly those corresponding to the convergent mode of deformation. Since this drawback was already observed in the BBM-driven simulations of the Lagrangian sea ice model neXtSIM presented in Ólason et al. (2022), we think that it probably shows an intrinsic limitation of the current BBM rheological model, an issue that certainly merits investigating and fixing in the future. Finally, we show that the observed spatial scaling invariance property of sea ice deformation, and in particular its multi-fractal nature, is fairly well captured by the BBM-driven simulation but with a slightly lower degree of multi-fractality.

Appendix A: Nomenclature

A1 Table of symbols used in the text

Symbol	Definition	Units
m	mass of ice and snow per unitarea	kg m−2
ρi	density of sea ice	kg m−3
ρw	density of seawater	kg m−3
$\mathit{u}\equiv (u,v)$	sea ice velocity	m s−1
A	sea ice fraction	–
h	sea ice thickness	m
g	acceleration of gravity	m s−2
f	Coriolis frequency	s−1
k	vertical unit vector (z axis)	s−1
H	sea surface height	m
τ	ice–ocean stress	Pa
τa	wind (ice–atmosphere) stress	Pa
σ	internal stress tensor (2 × 2)	Pa
$\dot{\mathit{\epsilon }}$	strain rate tensor (2 × 2)	s−1
d	damage of sea ice	–
Δx	local resolution (size) of the gridmesh	m
C	compaction parameter	–
α	damage parameter (Dansereau, 2016)	–
E0,E	elastic modulus of undamaged and damaged ice	Pa
λ0,λ	apparent viscous relaxation time of undamagedand damaged ice	s
$\stackrel{\mathrm{̃}}{P}$	BBM-specific ridging term	–
Pmax	ridging threshold	Pa
P0	scaling parameter for Pmax	Pa
h0	reference ice thickness for Pmax	m
c	sea ice cohesion	Pa
ν	Poisson's ratio	–
σI	isotropic normal stress (firstinvariant of stress tensor)	Pa
σII	maximum shear stress (secondinvariant of stress tensor)	Pa
μ	internal friction coefficient	–
N	upper limit for compressivestress	Pa
CE	propagation speed of an elasticshear wave	m s−1
td	characteristic timescale for thepropagation of damage	s
dcrit	damage increment (Dansereau, 2016)	–
kth	healing constant for damage	K s
ΔTh	temperature difference betweenbottom and surface of ice	K

A2 Acronyms

NEMO	Nucleus for European Modeling of the Ocean
SI3	Sea-Ice modeling Integrated Initiative (sea icecomponent of NEMO)
OCE	3D ocean component of NEMO
SAS	Stand-alone surface module of NEMO (i.e.,SI3 stand-alone)
LIM	Louvain-La-Neuve sea-Ice Model
BBM	Brittle Bingham–Maxwell (rheology)
MEB	Maxwell elasto-brittle (rheology)
VP	Viscous–plastic (rheology)
FD	Finite difference (method)
CN	Cross-nudging (treatment)
PDF	Probability density function
LKFs	Linear kinematic features
GCM	General circulation model
OGCM	Ocean general circulation model
SST	Sea surface temperature
SSH	Sea surface height
RGPS	RADARSAT Geophysical Processor System(dataset)

A3 Notations related to the discretization on the E-grid

The bar + superscript notation refers to a spatial interpolation; $\stackrel{-X}{\mathit{\varphi }}$ is field ϕ interpolated onto X points.

Interpolation from T- to F-points, or conversely, is based on the average of the four nearest surrounding points (Fig. 3a).

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\stackrel{-\mathrm{T}}{\mathit{\varphi }}}_{i,j}=& \mathrm{1}/\mathrm{4}({\mathit{\varphi }}_{i,j}+{\mathit{\varphi }}_{i-\mathrm{1},j}+{\mathit{\varphi }}_{i-\mathrm{1},j-\mathrm{1}}+{\mathit{\varphi }}_{i,j-\mathrm{1}})\\ & \text{ (if }\mathit{\varphi }\text{ defined at F)}\\ {\stackrel{-\mathrm{F}}{\mathit{\varphi }}}_{i,j}=& \mathrm{1}/\mathrm{4}({\mathit{\varphi }}_{i+\mathrm{1},j+\mathrm{1}}+{\mathit{\varphi }}_{i,j+\mathrm{1}}+{\mathit{\varphi }}_{i,j}+{\mathit{\varphi }}_{i+\mathrm{1},j})\\ & \text{ (if }\mathit{\varphi }\text{ defined at T)}\end{array}\end{array}$$

Note: surrounding points located on land or open-boundary cells are excluded from the averaging.

For the interpolation from tracer (T or F) to velocity (U or V) points, or conversely, only the two nearest surrounding points are used.

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{\stackrel{-\mathrm{U}}{\mathit{\varphi }}}_{i,j}=& \mathrm{1}/\mathrm{2}({\mathit{\varphi }}_{i+\mathrm{1},j}+{\mathit{\varphi }}_{i,j})\text{ (if }\mathit{\varphi }\text{ defined at T)}\\ {\stackrel{-\mathrm{V}}{\mathit{\varphi }}}_{i,j}=& \mathrm{1}/\mathrm{2}({\mathit{\varphi }}_{i,j+\mathrm{1}}+{\mathit{\varphi }}_{i,j})\text{ (if }\mathit{\varphi }\text{ defined at T)}\\ {\stackrel{-\mathrm{U}}{\mathit{\varphi }}}_{i,j}=& \mathrm{1}/\mathrm{2}({\mathit{\varphi }}_{i,j}+{\mathit{\varphi }}_{i,j-\mathrm{1}})\text{ (if }\mathit{\varphi }\text{ defined at F)}\\ {\stackrel{-\mathrm{V}}{\mathit{\varphi }}}_{i,j}=& \mathrm{1}/\mathrm{2}({\mathit{\varphi }}_{i,j}+{\mathit{\varphi }}_{i-\mathrm{1},j})\text{ (if }\mathit{\varphi }\text{ defined at F)}\end{array}\end{array}$$

The “hat” notation $\widehat{x}$ refers to the F-centric counterpart of x, with x being a prognostic scalar or tensor (rank 1 or 2) defined in the T-centric grid (mind that if x is the element of a tensor, $\widehat{x}$ is not necessarily defined on F-points). Example: $\widehat{d}$ and ${\widehat{\mathit{\sigma }}}_{\mathrm{11}}$ are prognostic fields defined on F-points (the natural location for d and σ₁₁ on the C-grid is the T-point); similarly, ${\widehat{\mathit{\sigma }}}_{\mathrm{12}}$ is defined on T-points (the natural location for σ₁₂ on the C-grid is the F-point).

A4 Miscellaneous notations

$\underset{\mathrm{¯}}{\mathit{x}}$	symmetric tensor x expressed in its Voigt form
x(i)	initial estimate of variable x
at X	on the X points of the grid
σkl	vertically integrated components of tensor σ
	$\to {\mathit{\sigma }}_{\text{kl}}\equiv h{\mathit{\sigma }}_{\text{kl}}$ if σkl defined at T
	$\to {\mathit{\sigma }}_{\text{kl}}\equiv {\stackrel{\mathrm{‾}}{h}}^F{\mathit{\sigma }}_{\text{kl}}$ if σkl defined at F
$\stackrel{▿}{\mathit{x}}$	upper-convected time derivative of symmetric(rank 2) tensor x
$\stackrel{△}{\mathit{x}}$	lower-convected time derivative of symmetric(rank 2) tensor x

A5 Table of symbols related to the numerical implementation

Symbol	Definition	Units
ΔT	advective time step for the advection and the thermodynamics	s
Δt	dynamical time stepspecific to BBM (time splitting)	s
Ns	$\equiv \mathrm{\Delta }T/\mathrm{\Delta }t$, time-splittingparameter	–
k	time level index of timesplitting ($\mathrm{1}\le k\le N_{\mathrm{s}}$)	–
$A,h,d$	ice concentration, thickness, and damage of iceat T	–, m, –
$\stackrel{-\mathrm{F}}{A},\stackrel{-\mathrm{F}}{h}$	ice concentrationand thickness interpolated at F	–, m
$\widehat{d}$	damage of ice at F	–
$\dot{\mathit{\epsilon }}\equiv \left({\dot{\mathit{\epsilon }}}_{\mathrm{11}},{\dot{\mathit{\epsilon }}}_{\mathrm{22}},{\dot{\mathit{\epsilon }}}_{\mathrm{12}}\right)$	strain rate tensor (2 × 2)of the T-centric cell	s−1
$\widehat{\dot{\mathit{\epsilon }}}\equiv \left({\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{11}},{\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{22}},{\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{12}}\right)$	strain rate tensor (2 × 2)of the F-centric cell	s−1
$\mathit{\sigma }\equiv ({\mathit{\sigma }}_{\mathrm{11}},{\mathit{\sigma }}_{\mathrm{22}},{\mathit{\sigma }}_{\mathrm{12}})$	internal stress tensor(2×2) of the T-centric cell	Pa
$\widehat{\mathit{\sigma }}\equiv \left({\widehat{\mathit{\sigma }}}_{\mathrm{11}},{\widehat{\mathit{\sigma }}}_{\mathrm{22}},{\widehat{\mathit{\sigma }}}_{\mathrm{12}}\right)$	internal stress tensor(2×2) of the F-centric cell	Pa
$\stackrel{-\mathrm{U}}{A},\stackrel{-\mathrm{V}}{A}$	ice concentration interpolated at U and at V	m

$\stackrel{-\mathrm{U}}{h},\stackrel{-\mathrm{V}}{h}$	ice thickness interpolatedat U and at V	m
u,v	ice velocity at the Δt level(at U and at V)	m s−1
$\widehat{u},\widehat{v}$	ice velocity at the Δt level(at V and at U)	m s−1
U,V	ice velocity at the ΔTlevel (at U and at V)	m s−1
$\widehat{U},\widehat{V}$	ice velocity at the ΔTlevel (at V and at U)	m s−1
$C_{\mathrm{D}}^{\left(\mathrm{o}\right)}$	ice–ocean drag coefficient	–
τx,τy	ice–ocean stress at U andat V	Pa
${\widehat{\mathit{\tau }}}_x,{\widehat{\mathit{\tau }}}_y$	ice–ocean stress at V andat U	Pa
uo,vo	surface ocean current at Uand at V	m s−1
γcn	cross-nudging coefficient	–
$C_{\mathrm{D}}^{\left(\mathrm{a}\right)}$	ice–atmosphere drag coefficient	–
Δ𝚃x	T-centered Δx that connects two neighboringU-points	m
Δ𝚃y	T-centered Δy that connects two neighboringV-points	m
Δ𝙵x	F-centered Δx that connects two neighboringV-points	m
Δ𝙵y	F-centered Δy that connects two neighboringU-points	m
Δ𝚄x	U-centered Δx that connects two neighboringT-points	m
Δ𝚄y	U-centered Δy that connects two neighboringF-points	m
Δ𝚅x	V-centered Δx that connects two neighboringF-points	m
Δ𝚅y	V-centered Δy that connects two neighboringT-points	m

Appendix B: Algorithm and discretization

B1 Algorithm

Time-splitting loop (Δt)/for k=1 to N_s

• Compute elasticity $E,\widehat{E}$ and viscous relaxation time $\mathit{\lambda },\widehat{\mathit{\lambda }}$ as a function of damage $d^k,{\widehat{d}}^k$ and current sea ice concentration $A,\stackrel{-\mathrm{F}}{A}$ (Eqs. B5 and B6).

• Compute the normal stress invariant of σ^k and ${\widehat{\mathit{\sigma }}}^k$ → ${\mathit{\sigma }}_{\mathrm{I}}^k,{\widehat{\mathit{\sigma }}}_{\mathrm{I}}^k$ (Eq. B8).

• Compute $P_{\text{max}},{\widehat{P}}_{\text{max}}$ as a function of current sea ice thickness $h,\stackrel{-\mathrm{F}}{h}$ and concentration $A,\stackrel{-\mathrm{F}}{A}$ (Eq. B7).

• Compute $\stackrel{\mathrm{̃}}{P},\widehat{\stackrel{\mathrm{̃}}{P}}$ as a function of $P_{\text{max}},{\widehat{P}}_{\text{max}}$ and ${\mathit{\sigma }}_{\mathrm{I}},{\widehat{\mathit{\sigma }}}_{\mathrm{I}}$ (Eq. B9).

• Compute the three components of each strain rate tensor $\dot{\mathit{\epsilon }},\widehat{\dot{\mathit{\epsilon }}}$ based on sea ice velocities at time level k (Eqs. B1, B2, B3, and B4).

• Calculate an initial prognostic estimate of the stress tensors at time level $k+\mathrm{1}\to {\mathit{\sigma }}^{\left(\mathrm{i}\right)k+\mathrm{1}}$ and ${\widehat{\mathit{\sigma }}}^{\left(\mathrm{i}\right)k+\mathrm{1}}$ (Eq. B11).

• Apply cross-nudging between the vertically integrated σ^(i)k+1 and ${\widehat{\mathit{\sigma }}}^{\left(\mathrm{i}\right)k+\mathrm{1}}$ (Eq. 16).

• Apply a Mohr–Coulomb test on σ^(i)k+1 and ${\widehat{\mathit{\sigma }}}^{\left(\mathrm{i}\right)k+\mathrm{1}}$.

  • Compute the two invariants of σ^(i)k+1 and ${\widehat{\mathit{\sigma }}}^{\left(\mathrm{i}\right)k+\mathrm{1}}$ → ${\mathit{\sigma }}_{\mathrm{I}}^{\left(\mathrm{i}\right)k+\mathrm{1}},{\mathit{\sigma }}_{\text{II}}^{\left(\mathrm{i}\right)k+\mathrm{1}}$ and ${\widehat{\mathit{\sigma }}}_{\mathrm{I}}^{\left(\mathrm{i}\right)k+\mathrm{1}},{\widehat{\mathit{\sigma }}}_{\text{II}}^{\left(\mathrm{i}\right)k+\mathrm{1}}$ (Eq. B8).

  • Compute d_crit and ${\widehat{d}}_{\text{crit}}$ based on ${\mathit{\sigma }}_{\mathrm{I}}^{\left(\mathrm{i}\right)k+\mathrm{1}},{\mathit{\sigma }}_{\text{II}}^{\left(\mathrm{i}\right)k+\mathrm{1}}$ and ${\widehat{\mathit{\sigma }}}_{\mathrm{I}}^{\left(\mathrm{i}\right)k+\mathrm{1}},{\widehat{\mathit{\sigma }}}_{\text{II}}^{\left(\mathrm{i}\right)k+\mathrm{1}}$ (Eq. B12).

• Calculate a prognostic estimate of the stress tensors and damage at time level $k+\mathrm{1}\to {\mathit{\sigma }}^{k+\mathrm{1}}$, d^k+1 and ${\widehat{\mathit{\sigma }}}^{k+\mathrm{1}}$, ${\widehat{d}}^{k+\mathrm{1}}$.

  • Where $\mathrm{0}<d_{\text{crit}}<\mathrm{1}$ and/or $\mathrm{0}<{\widehat{d}}_{\text{crit}}<\mathrm{1}$ (over-critical stress state),
    → damage growth and stress adjustment (Eq. B13).

  • Elsewhere,
    → no damage growth and no stress adjustment (Eq. B14).

• Compute the divergence of the vertically integrated σ^k+1 and ${\widehat{\mathit{\sigma }}}^{k+\mathrm{1}}$ (Eqs. B16 and B17).

• Calculate a prognostic estimate of sea ice velocity at time level $k+\mathrm{1}\to u^{k+\mathrm{1}},v^{k+\mathrm{1}}$ and ${\widehat{u}}^{k+\mathrm{1}},{\widehat{v}}^{k+\mathrm{1}}$ (Eqs. B19 and B18).

SI3 (advective) time step (ΔT):

• BBM rheology (time-splitting loop above)

• Advection of generic SI3 prognostic tracers (A, h, etc.) at T-points using U and V

• Advection of d, σ₁₁, σ₂₂, and ${\widehat{\mathit{\sigma }}}_{\mathrm{12}}$ at T-points using U and V

• Advection of $\widehat{d}$, ${\widehat{\mathit{\sigma }}}_{\mathrm{11}}$, ${\widehat{\mathit{\sigma }}}_{\mathrm{22}}$, and σ₁₂ at F-points using $\widehat{U}$ and $\widehat{V}$

• Healing of damage (d and $\widehat{d}$) (Eq.14)

• Thermodynamics module of SI3 (update of A, h, etc.)

B2 Update of internal stress tensor in the T- and F-centric worlds

B2.1 Divergence, shear, and strain rate tensor of ice velocity

Following Hunke and Dukowicz (2002), here is how the components of the strain rate of the sea ice velocity vector are computed on the T- and F-centric grids based on the finite-difference method.

• Divergence rate $({\partial }_{x}u+{\partial }_{y}v)$:

  $$\begin{array}{}\text{(B1)}& \begin{array}{rl}{D}_{i,j}& ={\displaystyle \frac{\begin{array}{c}[{\mathrm{\Delta }}^{\mathtt{U}}yu{]}_{i,j}-[{\mathrm{\Delta }}^{\mathtt{U}}yu{]}_{i-\mathrm{1},j}\\ +[{\mathrm{\Delta }}^{\mathtt{V}}xv{]}_{i,j}-[{\mathrm{\Delta }}^{\mathtt{V}}xv{]}_{i,j-\mathrm{1}}\end{array}}{[{\mathrm{\Delta }}^{\mathtt{T}}x{\mathrm{\Delta }}^{\mathtt{T}}y{]}_{i,j}}},\\ {\widehat{D}}_{i,j}& ={\displaystyle \frac{\begin{array}{c}[{\mathrm{\Delta }}^{\mathtt{V}}y\widehat{u}{]}_{i+\mathrm{1},j}-[{\mathrm{\Delta }}^{\mathtt{V}}y\widehat{u}{]}_{i,j}\\ +[{\mathrm{\Delta }}^{\mathtt{U}}x\widehat{v}{]}_{i,j+\mathrm{1}}-[{\mathrm{\Delta }}^{\mathtt{U}}x\widehat{v}{]}_{i,j}\end{array}}{[{\mathrm{\Delta }}^{\mathtt{F}}x{\mathrm{\Delta }}^{\mathtt{F}}y{]}_{i,j}}}.\end{array}\end{array}$$

• Tension rate $({\partial }_{x}u-{\partial }_{y}v)$:

  $$\begin{array}{}\text{(B2)}& \begin{array}{rl}& T_{i,j}=\\ & {\displaystyle \frac{\begin{array}{c}\left(\right[u/{\mathrm{\Delta }}^{\mathtt{U}}y{]}_{i,j}-[u/{\mathrm{\Delta }}^{\mathtt{U}}y{]}_{i-\mathrm{1},j})[{\mathrm{\Delta }}^{\mathtt{T}}{y}^{\mathrm{2}}{]}_{i,j}\\ -\left(\right[v/{\mathrm{\Delta }}^{\mathtt{V}}x{]}_{i,j}-[v/{\mathrm{\Delta }}^{\mathtt{V}}x{]}_{i,j-\mathrm{1}})[{\mathrm{\Delta }}^{\mathtt{T}}{x}^{\mathrm{2}}{]}_{i,j}\end{array}}{[{\mathrm{\Delta }}^{\mathtt{T}}x{\mathrm{\Delta }}^{\mathtt{T}}y{]}_{i,j}}},\\ & {\widehat{T}}_{i,j}=\\ & {\displaystyle \frac{\begin{array}{c}\left(\right[\widehat{u}/{\mathrm{\Delta }}^{\mathtt{V}}y{]}_{i+\mathrm{1},j}-[\widehat{u}/{\mathrm{\Delta }}^{\mathtt{V}}y{]}_{i,j})[{\mathrm{\Delta }}^{\mathtt{F}}{y}^{\mathrm{2}}{]}_{i,j}\\ -\left(\right[\widehat{v}/{\mathrm{\Delta }}^{\mathtt{U}}x{]}_{i,j+\mathrm{1}}-[\widehat{v}/{\mathrm{\Delta }}^{\mathtt{U}}x{]}_{i,j})[{\mathrm{\Delta }}^{\mathtt{F}}{x}^{\mathrm{2}}{]}_{i,j}\end{array}}{[{\mathrm{\Delta }}^{\mathtt{F}}x{\mathrm{\Delta }}^{\mathtt{F}}y{]}_{i,j}}}.\end{array}\end{array}$$

• Shear rate $({\partial }_{y}u+{\partial }_{x}v)$:

  $$\begin{array}{}\text{(B3)}& \begin{array}{rl}& S_{i,j}=\\ & {\displaystyle \frac{\begin{array}{c}\left(\right[u/{\mathrm{\Delta }}^{\mathtt{U}}x{]}_{i,j+\mathrm{1}}-[u/{\mathrm{\Delta }}^{\mathtt{U}}x{]}_{i,j})[{\mathrm{\Delta }}^{\mathtt{F}}{x}^{\mathrm{2}}{]}_{i,j}\\ +\left(\right[v/{\mathrm{\Delta }}^{\mathtt{V}}y{]}_{i+\mathrm{1},j}-[v/{\mathrm{\Delta }}^{\mathtt{V}}y{]}_{i,j})[{\mathrm{\Delta }}^{\mathtt{F}}{y}^{\mathrm{2}}{]}_{i,j}\end{array}}{[{\mathrm{\Delta }}^{\mathtt{F}}x{\mathrm{\Delta }}^{\mathtt{F}}y{]}_{i,j}}},\\ & {\widehat{S}}_{i,j}=\\ & {\displaystyle \frac{\begin{array}{c}\left(\right[\widehat{u}/{\mathrm{\Delta }}^{\mathtt{V}}x{]}_{i,j}-[\widehat{u}/{\mathrm{\Delta }}^{\mathtt{V}}x{]}_{i,j-\mathrm{1}})[{\mathrm{\Delta }}^{\mathtt{T}}{x}^{\mathrm{2}}{]}_{i,j}\\ +\left(\right[\widehat{v}/{\mathrm{\Delta }}^{\mathtt{U}}y{]}_{i,j}-[\widehat{v}/{\mathrm{\Delta }}^{\mathtt{U}}y{]}_{i-\mathrm{1},j})[{\mathrm{\Delta }}^{\mathtt{T}}{y}^{\mathrm{2}}{]}_{i,j}\end{array}}{[{\mathrm{\Delta }}^{\mathtt{T}}x{\mathrm{\Delta }}^{\mathtt{T}}y{]}_{i,j}}}.\end{array}\end{array}$$

From the above, the three components of the 2D strain rate tensors are obtained.

$$\begin{array}{}\text{(B4)}& \begin{array}{r}{\left(\begin{array}{c}{\dot{\mathit{\epsilon }}}_{\mathrm{11}}\\ {\dot{\mathit{\epsilon }}}_{\mathrm{22}}\\ {\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{12}}\end{array}\right)}_{i,j}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{c}{D}_{i,j}+T_{i,j}\\ D_{i,j}-T_{i,j}\\ {\widehat{S}}_{i,j}\end{array}\right)\\ {\left(\begin{array}{c}{\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{11}}\\ {\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{22}}\\ {\dot{\mathit{\epsilon }}}_{\mathrm{12}}\end{array}\right)}_{i,j}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\begin{array}{c}{\widehat{D}}_{i,j}+{\widehat{T}}_{i,j}\\ {\widehat{D}}_{i,j}-{\widehat{T}}_{i,j}\\ S_{i,j}\end{array}\right)\end{array}\end{array}$$

B2.2 Update of the stress tensors

• Elasticity and viscous relaxation time of damaged ice

  $$\begin{array}{}\text{(B5)}& {\displaystyle }\begin{array}{rl}E& =E_{\mathrm{0}}(\mathrm{1}-d)e^{C(\mathrm{1}-A)}\\ \widehat{E}& =E_{\mathrm{0}}(\mathrm{1}-\widehat{d})e^{C(\mathrm{1}-\stackrel{-\mathrm{F}}{A})}\end{array}\text{(B6)}& {\displaystyle }\begin{array}{rl}\mathit{\lambda }& =\mathit{\lambda }\mathrm{0}\left[\right(\mathrm{1}-d)e^{C(\mathrm{1}-A)}{]}^{\mathit{\alpha }-\mathrm{1}}\\ \widehat{\mathit{\lambda }}& =\mathit{\lambda }\mathrm{0}\left[\right(\mathrm{1}-\widehat{d})e^{C(\mathrm{1}-\stackrel{-\mathrm{F}}{A})}{]}^{\mathit{\alpha }-\mathrm{1}}\end{array}\end{array}$$

Note that it is the averaged value of A at F-points, $\stackrel{-\mathrm{F}}{A}$, that is used in the equations for the F-centric grid.

• Ridging threshold

  $$\begin{array}{}\text{(B7)}& \begin{array}{rl}{P}_{\text{max}}& =P\mathrm{0}[h/h_{\mathrm{0}}{]}^{\mathrm{3}/\mathrm{2}}{e}^{C(\mathrm{1}-A)}\\ {\widehat{P}}_{\text{max}}& =P\mathrm{0}{\left[\stackrel{-\mathrm{F}}{h}/h_{\mathrm{0}}\right]}^{\mathrm{3}/\mathrm{2}}{e}^{C(\mathrm{1}-\stackrel{-\mathrm{F}}{A})}\end{array}\end{array}$$

Note that it is the averaged value of h at F-points, $\stackrel{-\mathrm{F}}{h}$, that is used in the second equation.

• Invariants of stress tensor

  $$\begin{array}{}\text{(B8)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{I}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}),{\mathit{\sigma }}_{\text{II}}=\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\widehat{\mathit{\sigma }}}_{\mathrm{12}}^{\mathrm{2}}}\\ & {\widehat{\mathit{\sigma }}}_{\mathrm{I}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\widehat{\mathit{\sigma }}}_{\mathrm{11}}+{\widehat{\mathit{\sigma }}}_{\mathrm{22}}),{\widehat{\mathit{\sigma }}}_{\text{II}}=\sqrt{{\left({\displaystyle \frac{{\widehat{\mathit{\sigma }}}_{\mathrm{11}}-{\widehat{\mathit{\sigma }}}_{\mathrm{22}}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}}\end{array}\end{array}$$

• $\stackrel{\mathrm{̃}}{P}$ term

  $$\begin{array}{}\text{(B9)}& \begin{array}{rl}& \stackrel{\mathrm{̃}}{P}=\left\{\begin{array}{ll}\frac{{\mathit{\sigma }}_{\mathrm{I}}}{-P_{\text{max}}}& \text{for }{\mathit{\sigma }}_{\mathrm{I}}<-P_{\text{max}}\\ -\mathrm{1}& \text{for }-P_{\text{max}}\le {\mathit{\sigma }}_{\mathrm{I}}<\mathrm{0}\\ \mathrm{0}& \text{for }{\mathit{\sigma }}_{\mathrm{I}}>\mathrm{0}\end{array}\right.\\ & \widehat{\stackrel{\mathrm{̃}}{P}}=\left\{\begin{array}{ll}\frac{{\widehat{\mathit{\sigma }}}_{\mathrm{I}}}{-{\widehat{P}}_{\text{max}}}& \text{for }{\widehat{\mathit{\sigma }}}_{\mathrm{I}}<-{\widehat{P}}_{\text{max}}\\ -\mathrm{1}& \text{for }-{\widehat{P}}_{\text{max}}\le {\widehat{\mathit{\sigma }}}_{\mathrm{I}}<\mathrm{0}\\ \mathrm{0}& \text{for }{\widehat{\mathit{\sigma }}}_{\mathrm{I}}>\mathrm{0}\end{array}\right.\end{array}\end{array}$$

• Multiplicator for stress update

  $$\begin{array}{}\text{(B10)}& \begin{array}{rl}\mathrm{\Omega }& ={\displaystyle \frac{\mathit{\lambda }}{\mathit{\lambda }+(\mathrm{1}+\stackrel{\mathrm{̃}}{P})\mathrm{\Delta }t}}\\ \widehat{\mathrm{\Omega }}& ={\displaystyle \frac{\widehat{\mathit{\lambda }}}{\widehat{\mathit{\lambda }}+(\mathrm{1}+\widehat{\stackrel{\mathrm{̃}}{P}})\mathrm{\Delta }t}}\end{array}\end{array}$$

• Initial update of stress tensor

  $$\begin{array}{}\text{(B11)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathrm{11}}^{\left(\mathrm{i}\right)k+\mathrm{1}}& =\mathrm{\Omega }\left[E\mathrm{\Delta }t{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}}\left({\dot{\mathit{\epsilon }}}_{\mathrm{11}}^k+\mathit{\nu }{\dot{\mathit{\epsilon }}}_{\mathrm{22}}^k\right)+{\mathit{\sigma }}_{\mathrm{11}}^k\right]\\ {\mathit{\sigma }}_{\mathrm{22}}^{\left(\mathrm{i}\right)k+\mathrm{1}}& =\mathrm{\Omega }\left[E\mathrm{\Delta }t{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}}\left(\mathit{\nu }{\dot{\mathit{\epsilon }}}_{\mathrm{11}}^k+{\dot{\mathit{\epsilon }}}_{\mathrm{22}}^k\right)+{\mathit{\sigma }}_{\mathrm{22}}^k\right]\\ {\widehat{\mathit{\sigma }}}_{\mathrm{12}}^{\left(\mathrm{i}\right)k+\mathrm{1}}& =\widehat{\mathrm{\Omega }}\left[\widehat{E}\mathrm{\Delta }t{\displaystyle \frac{\mathrm{1}-\mathit{\nu }}{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}}{\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{12}}^k+{\widehat{\mathit{\sigma }}}_{\mathrm{12}}^k\right]\\ {\widehat{\mathit{\sigma }}}_{\mathrm{11}}^{\left(\mathrm{i}\right)k+\mathrm{1}}& =\widehat{\mathrm{\Omega }}\left[\widehat{E}\mathrm{\Delta }t{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}}\left({\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{11}}^k+\mathit{\nu }{\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{22}}^k\right)+{\widehat{\mathit{\sigma }}}_{\mathrm{11}}^k\right]\\ {\widehat{\mathit{\sigma }}}_{\mathrm{22}}^{\left(\mathrm{i}\right)k+\mathrm{1}}& =\widehat{\mathrm{\Omega }}\left[\widehat{E}\mathrm{\Delta }t{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}}\left(\mathit{\nu }{\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{11}}^k+{\widehat{\dot{\mathit{\epsilon }}}}_{\mathrm{22}}^k\right)+{\widehat{\mathit{\sigma }}}_{\mathrm{22}}^k\right]\\ {\mathit{\sigma }}_{\mathrm{12}}^{\left(\mathrm{i}\right)k+\mathrm{1}}& =\mathrm{\Omega }\left[E\mathrm{\Delta }t{\displaystyle \frac{\mathrm{1}-\mathit{\nu }}{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}}{\dot{\mathit{\epsilon }}}_{\mathrm{12}}^k+{\mathit{\sigma }}_{\mathrm{12}}^k\right]\end{array}\end{array}$$

• Damage increment

  $$\begin{array}{}\text{(B12)}& \begin{array}{rl}{d}_{\text{crit}}& =\left\{\begin{array}{ll}\frac{c}{{\mathit{\sigma }}_{\text{II}}^{\left(\mathrm{i}\right)}+\mathit{\mu }{\mathit{\sigma }}_{\mathrm{I}}^{\left(\mathrm{i}\right)}}& \text{if }{\mathit{\sigma }}_{\mathrm{I}}^{\left(\mathrm{i}\right)}>-N\\ \frac{-N}{{\mathit{\sigma }}_{\mathrm{I}}^{\left(\mathrm{i}\right)}}& \text{otherwise}\end{array}\right.\\ {\widehat{d}}_{\text{crit}}& =\left\{\begin{array}{ll}\frac{c}{{\widehat{\mathit{\sigma }}}_{\text{II}}^{\left(\mathrm{i}\right)}+\mathit{\mu }{\widehat{\mathit{\sigma }}}_{\mathrm{I}}^{\left(\mathrm{i}\right)}}& \text{if }{\widehat{\mathit{\sigma }}}_{\mathrm{I}}^{\left(\mathrm{i}\right)}>-N\\ \frac{-N}{{\widehat{\mathit{\sigma }}}_{\mathrm{I}}^{\left(\mathrm{i}\right)}}& \text{otherwise}\end{array}\right.\end{array}\end{array}$$

• Update of damage and stress tensors

  • In regions where $\mathrm{0}<d_{\text{crit}}<\mathrm{1}$, the following applies.

    $$\begin{array}{}\text{(B13)}& \begin{array}{rl}& \begin{array}{c}{d}^{k+\mathrm{1}}=\begin{array}{c}{d}^k+(\mathrm{1}-d_{\text{crit}})\\ \times (\mathrm{1}-d^k)\mathrm{\Delta }t/t_{\mathrm{d}}\end{array}\\ {\underset{\mathrm{¯}}{\mathit{\sigma }}}^{k+\mathrm{1}}=\begin{array}{c}{\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(i\right)k+\mathrm{1}}-\left(\mathrm{1}-d_{\text{crit}}\right)\\ \times {\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(i\right)k+\mathrm{1}}\mathrm{\Delta }t/t_{\mathrm{d}}\end{array}\end{array}\text{with }{t}_{\mathrm{d}}={\mathrm{\Delta }}^{\mathtt{T}}x\sqrt{{\displaystyle \frac{\mathrm{2}(\mathrm{1}+\mathit{\nu }){\mathit{\rho }}_{\mathrm{i}}}{E}}}\\ & \begin{array}{c}{\widehat{d}}^{k+\mathrm{1}}=\begin{array}{c}{\widehat{d}}^k+(\mathrm{1}-{\widehat{d}}_{\text{crit}})\\ \times (\mathrm{1}-{\widehat{d}}^k)\mathrm{\Delta }t/{\widehat{t}}_d\end{array}\\ {\widehat{\underset{\mathrm{¯}}{\mathit{\sigma }}}}^{k+\mathrm{1}}=\begin{array}{c}{\widehat{\underset{\mathrm{¯}}{\mathit{\sigma }}}}^{\left(i\right)k+\mathrm{1}}-\left(\mathrm{1}-{\widehat{d}}_{\text{crit}}\right)\\ \times {\widehat{\underset{\mathrm{¯}}{\mathit{\sigma }}}}^{\left(i\right)k+\mathrm{1}}\mathrm{\Delta }t/{\widehat{t}}_d\end{array}\end{array}\text{with }{\widehat{t}}_d={\mathrm{\Delta }}^{\mathtt{F}}x\sqrt{{\displaystyle \frac{\mathrm{2}(\mathrm{1}+\mathit{\nu }){\mathit{\rho }}_{\mathrm{i}}}{\widehat{E}}}}\end{array}\end{array}$$

  • Elsewhere, the following applies.

    $$\begin{array}{}\text{(B14)}& \begin{array}{rl}& d^{k+\mathrm{1}}=d^k\\ & {\underset{\mathrm{¯}}{\mathit{\sigma }}}^{k+\mathrm{1}}={\underset{\mathrm{¯}}{\mathit{\sigma }}}^{\left(i\right)k+\mathrm{1}}\\ & {\widehat{d}}^{k+\mathrm{1}}={\widehat{d}}^k\\ & {\widehat{\underset{\mathrm{¯}}{\mathit{\sigma }}}}^{k+\mathrm{1}}={\widehat{\underset{\mathrm{¯}}{\mathit{\sigma }}}}^{\left(i\right)k+\mathrm{1}}\end{array}\end{array}$$

B3 Momentum equation

As opposed to aEVP for which SI3 uses the scheme of Kimmritz et al. (2016, 2017), here we chose to solve the equation for momentum (in both the T- and F-centric worlds) using the implicit scheme approach of Bouillon et al. (2009).

B3.1 Divergence of the vertically integrated stress tensor

• Definition

  $$\begin{array}{}\text{(B15)}& \left(\begin{array}{c}{\mathrm{\Delta }}_x\\ {\mathrm{\Delta }}_y\end{array}\right)\equiv \left(\begin{array}{c}\frac{\partial \left(h{\mathit{\sigma }}_{\mathrm{11}}\right)}{\partial x}+\frac{\partial \left(h{\mathit{\sigma }}_{\mathrm{12}}\right)}{\partial y}\\ \frac{\partial \left(h{\mathit{\sigma }}_{\mathrm{22}}\right)}{\partial y}+\frac{\partial \left(h{\mathit{\sigma }}_{\mathrm{12}}\right)}{\partial x}\end{array}\right)\end{array}$$

• Discretized in the T-centric cell

  $$\begin{array}{}\text{(B16)}& \begin{array}{rl}& {\mathrm{\Delta }}_x^{k+\mathrm{1}}{|}_{i,j}={\displaystyle \frac{{\left[{\mathit{\sigma }}_{\mathrm{11}}^{k+\mathrm{1}}h{\mathrm{\Delta }}^{\mathtt{T}}{y}^{\mathrm{2}}\right]}_{i+\mathrm{1},j}-{\left[{\mathit{\sigma }}_{\mathrm{11}}^{k+\mathrm{1}}h{\mathrm{\Delta }}^{\mathtt{T}}{y}^{\mathrm{2}}\right]}_{i,j}}{[{\mathrm{\Delta }}^{\mathtt{U}}x{\mathrm{\Delta }}^{\mathtt{U}}{y}^{\mathrm{2}}{]}_{i,j}}}\\ & +{\displaystyle \frac{{\left[{\mathit{\sigma }}_{\mathrm{12}}^{k+\mathrm{1}}\stackrel{-\mathrm{F}}{h}{\mathrm{\Delta }}^{\mathtt{F}}{x}^{\mathrm{2}}\right]}_{i,j}-{\left[{\mathit{\sigma }}_{\mathrm{12}}^{k+\mathrm{1}}\stackrel{-\mathrm{F}}{h}{\mathrm{\Delta }}^{\mathtt{F}}{x}^{\mathrm{2}}\right]}_{i,j-\mathrm{1}}}{[{\mathrm{\Delta }}^{\mathtt{U}}y{\mathrm{\Delta }}^{\mathtt{U}}{x}^{\mathrm{2}}{]}_{i,j}}}\text{(at U)}\\ & {\mathrm{\Delta }}_y^{k+\mathrm{1}}{|}_{i,j}={\displaystyle \frac{{\left[{\mathit{\sigma }}_{\mathrm{22}}^{k+\mathrm{1}}h{\mathrm{\Delta }}^{\mathtt{T}}{x}^{\mathrm{2}}\right]}_{i,j+\mathrm{1}}-{\left[{\mathit{\sigma }}_{\mathrm{22}}^{k+\mathrm{1}}h{\mathrm{\Delta }}^{\mathtt{T}}{x}^{\mathrm{2}}\right]}_{i,j}}{[{\mathrm{\Delta }}^{\mathtt{V}}y{\mathrm{\Delta }}^{\mathtt{V}}{x}^{\mathrm{2}}{]}_{i,j}}}\\ & +{\displaystyle \frac{{\left[{\mathit{\sigma }}_{\mathrm{12}}^{k+\mathrm{1}}\widehat{h}{\mathrm{\Delta }}^{\mathtt{F}}{y}^{\mathrm{2}}\right]}_{i,j}-{\left[{\mathit{\sigma }}_{\mathrm{12}}^{k+\mathrm{1}}\widehat{h}{\mathrm{\Delta }}^{\mathtt{F}}{y}^{\mathrm{2}}\right]}_{i-\mathrm{1},j}}{[{\mathrm{\Delta }}^{\mathtt{V}}x{\mathrm{\Delta }}^{\mathtt{V}}{y}^{\mathrm{2}}{]}_{i,j}}}\text{(at V)}\end{array}\end{array}$$

• Discretized in the F-centric cell

  $$\begin{array}{}\text{(B17)}& \begin{array}{rl}& {\widehat{\mathrm{\Delta }}}_x^{k+\mathrm{1}}{|}_{i,j}={\displaystyle \frac{{\left[{\widehat{\mathit{\sigma }}}_{\mathrm{11}}^{k+\mathrm{1}}\stackrel{-\mathrm{F}}{h}{\mathrm{\Delta }}^{\mathtt{F}}{y}^{\mathrm{2}}\right]}_{i,j}-{\left[{\widehat{\mathit{\sigma }}}_{\mathrm{11}}^{k+\mathrm{1}}\stackrel{-\mathrm{F}}{h}{\mathrm{\Delta }}^{\mathtt{F}}{y}^{\mathrm{2}}\right]}_{i-\mathrm{1},j}}{[{\mathrm{\Delta }}^{\mathtt{V}}x{\mathrm{\Delta }}^{\mathtt{V}}{y}^{\mathrm{2}}{]}_{i,j}}}\\ & +{\displaystyle \frac{{\left[{\widehat{\mathit{\sigma }}}_{\mathrm{12}}^{k+\mathrm{1}}h{\mathrm{\Delta }}^{\mathtt{T}}{x}^{\mathrm{2}}\right]}_{i,j+\mathrm{1}}-{\left[{\widehat{\mathit{\sigma }}}_{\mathrm{12}}^{k+\mathrm{1}}h{\mathrm{\Delta }}^{\mathtt{T}}{x}^{\mathrm{2}}\right]}_{i,j}}{[{\mathrm{\Delta }}^{\mathtt{V}}y{\mathrm{\Delta }}^{\mathtt{V}}{x}^{\mathrm{2}}{]}_{i,j}}}\text{(at V)}\\ & {\widehat{\mathrm{\Delta }}}_y^{k+\mathrm{1}}{|}_{i,j}={\displaystyle \frac{{\left[{\widehat{\mathit{\sigma }}}_{\mathrm{22}}^{k+\mathrm{1}}\stackrel{-\mathrm{F}}{h}{\mathrm{\Delta }}^{\mathtt{F}}{x}^{\mathrm{2}}\right]}_{i,j}-{\left[{\widehat{\mathit{\sigma }}}_{\mathrm{22}}^{k+\mathrm{1}}\stackrel{-\mathrm{F}}{h}{\mathrm{\Delta }}^{\mathtt{F}}{x}^{\mathrm{2}}\right]}_{i,j-\mathrm{1}}}{[{\mathrm{\Delta }}^{\mathtt{U}}y{\mathrm{\Delta }}^{\mathtt{U}}{x}^{\mathrm{2}}{]}_{i,j}}}\\ & +{\displaystyle \frac{{\left[{\widehat{\mathit{\sigma }}}_{\mathrm{12}}^{k+\mathrm{1}}h{\mathrm{\Delta }}^{\mathtt{T}}{y}^{\mathrm{2}}\right]}_{i+\mathrm{1},j}-{\left[{\widehat{\mathit{\sigma }}}_{\mathrm{12}}^{k+\mathrm{1}}h{\mathrm{\Delta }}^{\mathtt{T}}{y}^{\mathrm{2}}\right]}_{i,j}}{[{\mathrm{\Delta }}^{\mathtt{U}}x{\mathrm{\Delta }}^{\mathtt{U}}{y}^{\mathrm{2}}{]}_{i,j}}}\text{(at U)}\end{array}\end{array}$$

B3.2 Update of sea ice velocity

For clarity, we gather the contributions of the wind stress, the Coriolis term, and the SSH tilt vectors in a single vector term named (R_x,R_y). This is because these three terms do not present any particular challenge to express with respect to the existing implementation of aEVP. Note, however, that with the E-grid no spatial interpolation is needed to express the discretized Coriolis term.

Implicitness of the scheme is introduced through the use of sea ice velocity at level k+1 in the estimate of the basal ice–water stress vector (τ_x,τ_y).

$$\begin{array}{}\text{(B18)}& \begin{array}{rl}& {\mathit{\tau }}_x=Z_x^k(u_{\mathrm{o}}^k-u^{k+\mathrm{1}})\\ & \text{with }{Z}_x^k=\stackrel{-\mathrm{U}}{A}{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{D}}^{\left(\mathrm{o}\right)}\sqrt{(u_{\mathrm{o}}^k-u^k{)}^{\mathrm{2}}+(\stackrel{-\mathrm{U}k}{v_o}-{\widehat{v}}^k{)}^{\mathrm{2}}}\text{(at U)}\\ & {\mathit{\tau }}_y=Z_y^k(v_{\mathrm{o}}^k-v^{k+\mathrm{1}})\\ & \text{with }{Z}_y^k=\stackrel{-\mathrm{V}}{A}{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{D}}^{\left(\mathrm{o}\right)}\sqrt{(\stackrel{-\mathrm{V}k}{u_o}-{\widehat{u}}^k{)}^{\mathrm{2}}+(v_{\mathrm{o}}^k-v^k{)}^{\mathrm{2}}}\text{(at V)}\\ & {\widehat{\mathit{\tau }}}_x={\widehat{Z}}_x^k(\stackrel{-\mathrm{V}k}{u_o}-{\widehat{u}}^{k+\mathrm{1}})\\ & \text{with }{\widehat{Z}}_x^k=\stackrel{-\mathrm{V}}{A}{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{D}}^{\left(\mathrm{o}\right)}\sqrt{(\stackrel{-\mathrm{V}k}{u_o}-{\widehat{u}}^k{)}^{\mathrm{2}}+(v_{\mathrm{o}}^k-v^k{)}^{\mathrm{2}}}\text{(at V)}\\ & {\widehat{\mathit{\tau }}}_y={\widehat{Z}}_y^k(\stackrel{-\mathrm{U}k}{v_o}-{\widehat{v}}^{k+\mathrm{1}})\\ & \text{with }{\widehat{Z}}_y^k=\stackrel{-\mathrm{U}}{A}{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{D}}^{\left(\mathrm{o}\right)}\sqrt{(u_{\mathrm{o}}^k-u^k{)}^{\mathrm{2}}+(\stackrel{-\mathrm{U}k}{v_o}-{\widehat{v}}^k{)}^{\mathrm{2}}}\text{(at U)}\end{array}\end{array}$$

Then, the discretized equation of momentum yields the expression of the two velocity components at time level k+1.

$$\begin{array}{}\text{(B19)}& \begin{array}{rl}{u}^{k+\mathrm{1}}=& {\displaystyle \frac{\frac{{\mathit{\rho }}_{\mathrm{i}}\stackrel{-\mathrm{U}}{h}}{\mathrm{\Delta }t}{u}^k+Z_x{u}_{\mathrm{o}}^k+{\mathrm{\Delta }}_x^{k+\mathrm{1}}+R_x^k}{\frac{{\mathit{\rho }}_{\mathrm{i}}\stackrel{-\mathrm{U}}{h}}{\mathrm{\Delta }t}+Z_x}}\text{(at U)}\\ v^{k+\mathrm{1}}=& {\displaystyle \frac{\frac{{\mathit{\rho }}_{\mathrm{i}}\stackrel{-\mathrm{V}}{h}}{\mathrm{\Delta }t}{v}^k+Z_y{v}_{\mathrm{o}}^k+{\mathrm{\Delta }}_y^{k+\mathrm{1}}+R_y^k}{\frac{{\mathit{\rho }}_{\mathrm{i}}\stackrel{-\mathrm{V}}{h}}{\mathrm{\Delta }t}+Z_y}}\text{(at V)}\\ {\widehat{u}}^{k+\mathrm{1}}=& {\displaystyle \frac{\frac{{\mathit{\rho }}_{\mathrm{i}}\stackrel{-\mathrm{V}}{h}}{\mathrm{\Delta }t}{\widehat{u}}^k+{\widehat{Z}}_x\stackrel{-\mathrm{V}k}{u_o}+{\widehat{\mathrm{\Delta }}}_x^{k+\mathrm{1}}+{\widehat{R}}_x^k}{\frac{{\mathit{\rho }}_{\mathrm{i}}\stackrel{-\mathrm{V}}{h}}{\mathrm{\Delta }t}+{\widehat{Z}}_x}}\text{(at V)}\\ {\widehat{v}}^{k+\mathrm{1}}=& {\displaystyle \frac{\frac{{\mathit{\rho }}_{\mathrm{i}}\stackrel{-\mathrm{U}}{h}}{\mathrm{\Delta }t}{\widehat{v}}^k+{\widehat{Z}}_y\stackrel{-\mathrm{U}k}{v_o}+{\widehat{\mathrm{\Delta }}}_y^{k+\mathrm{1}}+{\widehat{R}}_y^k}{\frac{{\mathit{\rho }}_{\mathrm{i}}\stackrel{-\mathrm{U}}{h}}{\mathrm{\Delta }t}+{\widehat{Z}}_y}}\text{(at U)}\end{array}\end{array}$$

Appendix C: Additional figures

Figure C1 Sea ice total deformation (instantaneous) in the test case described by Mehlmann et al. (2021) after 48 h of simulation with SI3 using the default SI3 aEVP rheology with N_EVP=100 and N_EVP=1000 (left-hand and right-hand column, respectively) at 2, 4, and 10 km resolution (first, second, and third row, respectively).

Figure C2 Sea ice total deformation (instantaneous) in the test case described by Mehlmann et al. (2021) after 48 h of simulation with SI3 using the newly implemented BBM rheology with N_s=100 and N_s=200 (left-hand and right-hand column, respectively) at 2, 4, and 10 km resolution (first, second, and third row, respectively).

Figure C3 PDFs of the total deformation rates at the 10 km spatial and 3 d temporal scale for the synthetic counterparts of RGPS data constructed using the simulated sea ice velocities of three BBM-driven SI3 experiments that use different time-derivative formulations for the stress tensor: only material derivative (“Material”, black line with circle markers), upper-convected time derivative as in SI3-BBM (“Upper-C”, blue line with square markers), and lower-convected time derivative (“Lower-C”, red line with star markers).