the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Implementation of a brittle sea ice rheology in an Eulerian, finitedifference, Cgrid modeling framework: impact on the simulated deformation of sea ice in the Arctic
Laurent Brodeau
Pierre Rampal
Einar Ólason
Véronique Dansereau
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, finitedifference, Arakawa Cgrid 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 Cgrid into the Arakawa Egrid. 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.
 Article
(9491 KB)  Fulltext XML
 BibTeX
 EndNote
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 multicomponent Earth system, such as regional and global climate simulations, and even for shortterm 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 multifractal 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 elastobrittle rheology (MEB). MEB was implemented into neXtSIM, a largescale dynamical–thermodynamical Lagrangian finiteelement sea ice model (Rampal et al., 2016), and used to evaluate the performance of this new rheology in a realistic panArctic 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 BBMdriven 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 CMIPlike simulations with neXtSIM in a numerically efficient manner remains challenging because the numerical coupling of neXtSIM to a thirdparty – 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 CMIPclass 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 kilometerscale 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 finitedifference method on a staggered grid. As suggested by the work of Plante et al. (2020), when discretized on the Arakawa Cgrid (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 Cgrid 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 finitedifference, Cgridbased sea ice model. We describe how we have implemented this new approach into SI3 based on the expansion of the Cgrid into an Arakawa Egrid.
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 panArctic 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 nonbrittle 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 panArctic 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.1 Governing equations and constitutive law
The twodimensional, vertically integrated, momentum equation for sea ice reads
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 twodimensional (plane stresses) case, the stress tensor is written as
In general, a constitutive law relates σ to the strain rate tensor $\dot{\mathit{\epsilon}}$, defined as follows:
As derived by Ólason et al. (2022) (their Eq. 20), the BBM constitutive model yields
where E and λ are the elastic modulus and apparent viscous relaxation time of the ice, K is the elastic stiffness tensor, $\stackrel{\mathrm{\u0303}}{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
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 reincreases 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{\u0303}}{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 viscoelastic. Ó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.,
where σ_{I} is the (isotropic) normal stress and P_{max} is the ridging threshold defined as
We follow Dansereau et al. (2016) and Ólason et al. (2022) in using a twostep approach to solve Eq. (4). As the first step, an initial estimate of $\underset{\mathrm{\xaf}}{\mathit{\sigma}}$, denoted as ${\underset{\mathrm{\xaf}}{\mathit{\sigma}}}^{\left(\mathrm{i}\right)}$, is calculated assuming no change in damage:
Then, as the second step, the following test and adjustment are performed on the state of stress: if ${\underset{\mathrm{\xaf}}{\mathit{\sigma}}}^{\left(\mathrm{i}\right)}$ is locally overcritical, i.e., located outside of the Mohr–Coulomb damage criterion (Fig. 1), an increment in ice damage, d_{crit}, is applied such that
where ${\underset{\mathrm{\xaf}}{\mathit{\sigma}}}^{\left(\mathrm{i}\right)}$ is the local value of the overcritical stress, and $\underset{\mathrm{\xaf}}{\mathit{\sigma}}$ is the corresponding postfailure (i.e., postdamage) stress. As discussed in Dansereau et al. (2016) and Plante and Tremblay (2021), d_{crit} is used to scale overcritical 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
where t_{d} is a characteristic timescale for damage propagation. In the BBM framework, d_{crit} is expressed as follows:
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:
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.
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 elastoviscobrittle 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 viscoelastic material in which irreversible deformations dissipate the stresses. As opposed to viscous–plastic frameworks, elasticity is therefore a physical and nonnegligible 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 overcritical 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.
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 σ^{(i)} (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 timesplitting approach. N_{s}, the number of (Δtlong) integrations to perform during one advective time step (ΔT), is imposed by ΔT and Δx, the horizontal resolution of the grid:
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 finitedifference method (FD) on the Arakawa Cgrid (Arakawa and Lamb, 1977). The use of the Cgrid 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 interpolationrelated errors as well as extra computational load.
As shown in Fig. 2a, on the Cgrid, tracers are defined at the cell centers, hereafter referred to as the Tpoint, while the x and y components of vectors are defined at the center of the righthand and upper edges of each cell, respectively (hereafter Upoint and Vpoint). The point located at the upperright corner of each cell, known as the vorticity point, is referred to as the Fpoint. In the literature, this vorticity point is sometimes located at the bottomleft corner of the cell, with the U and Vpoints possibly located at the lefthand and lower edges of the cell, and is sometimes referred to as the Zpoint (Losch et al., 2010; Plante et al., 2020).
Regardless of the type of rheology considered, the main challenge posed by the Cgrid 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 Tpoint and the shear rate ${\dot{\mathit{\epsilon}}}_{\mathrm{12}}$ defined at the Fpoint. Based on the constitutive law (Eq. 4), the same applies to the stress tensor σ. This staggering between the diagonal and offdiagonal 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 Vpoints. However, this staggering becomes an issue whenever the parameterization of the constitutive law requires ${\dot{\mathit{\epsilon}}}_{\mathrm{12}}$ or σ_{12} to be known at a Tpoint. 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 σ_{12}, respectively, to be known at Tpoints. Moreover, in brittle rheologies, a value of d is required not only at the Tpoint, but also at the Fpoint to estimate E and λ (Eqs. 5 and 6) needed to update ${\mathit{\sigma}}_{\mathrm{12}}^{\left(\mathrm{i}\right)}$ (Eq. 9). On the Cgrid, a common way to interpolate a scalar defined at Fpoints onto Tpoints is to simply use the average of this scalar on the four surrounding F points and conversely to interpolate from T to Fpoints. 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 Tpoints. Later on, the term $P/\mathrm{\Delta}$ is also interpolated from T to Fpoints in order to estimate σ_{12}. In their implementation of MEB, Plante et al. (2020) also use this approach to interpolate the damage tracer at Fpoints. However, they report that using the same approach to estimate σ_{12}, and hence σ_{II}, at Tpoints 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 σ_{12} that is defined at Tpoints. This additional σ_{12} is updated at each time step using – as an increment – the average of the four σ_{12} increments computed at the surrounding Fpoints.
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 Cgrid 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 Fpoint than what the rheology, together with Mohr–Coulomb test, would have produced at this location if explicitly solved on the Fpoint. This is because the fourpoint 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 Cgrid, the implementation of the advection of σ_{12} (Fpoint) in a way consistent (in terms of the advection scheme used) with that used for σ_{11} and σ_{22} (Tpoint) is somewhat challenging. That is because the advection of a scalar defined at the Fpoint, using the same scheme as that used for the advection of scalars at Tpoints, requires the existence of a u and a v at V and Upoints, respectively. These later considerations have prompted us to consider the use of a new spatial discretization approach for the implementation of BBM on the Cgrid.
2.3.1 The Egrid 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 Vpoints, while its y component, $\widehat{v}$, is defined at Upoints. Similarly, the damage tracer is also duplicated, with an additional occurrence at the upperright corners of the grid cell, i.e., at Fpoints. This grid staggering arrangement corresponds to that of the Arakawa Egrid (Arakawa and Lamb, 1977; Janjić, 1984; MaierReimer et al., 1993).
As suggested by Fig. 3b, the Egrid can be seen as a superposition of two Cgrids, in which the cell center of the additional Cgrid coincides with the upperright corner of the original Cgrid. For convenience, we will refer to these two grids as Fcentric (additional) and Tcentric (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 Eaugmented Cgrid to rheology and advection modules only. The rest of the code, which includes the thermodynamics, remains unmodified and relies entirely on the standard Cgrid. As such, only rheologyspecific tracers are defined in the Egrid fashion, i.e., at both T and Fpoints. 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 Tpoint only. Consequently, these tracers are interpolated at the Fpoint within the rheology module whenever needed.
To summarize, in the proposed rheologyspecific Eaugmented Cgrid approach, as shown in Fig. 3, the conventional Cgrid model variables are augmented with (i) the uvelocity component at Vpoints and vvelocity component at Upoints; (ii) the ice damage, σ_{11}, and σ_{22} at Fpoints; and (iii) σ_{12} at Tpoints. 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 Fcentric 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 Fcentric 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 Tcentric 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 Egrid in no way relates to an increase in the spatial resolution of the original Cgrid. That is because the FD discretization of spatial derivatives on the Egrid (see Appendix B) still relies on the same local spatial increment, i.e., Δx, as that of the original Cgrid, regardless of the subgrid considered (T or Fcentric).
2.3.2 The separation of solutions and how it is restrained
With the Eaugmented Cgrid approach, all rheologyspecific 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 Fcentric 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 Egrid 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 shallowwater 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 Egrid was used in the Hamburg LargeScale Geostrophic (LSG) model (MaierReimer et al., 1993) in order to achieve more accurate geostrophic balance, while retaining some advantages of the Cgrid 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 Egrid.
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 Fcentric grids occurs via the shear stress σ_{12} and the ice velocity vector. The estimate of ${\mathit{\sigma}}_{\mathrm{12}}^{\left(\mathrm{i}\right)}$ of the Tcentric grid (at Fpoint), based on Eq. (9), uses $\widehat{E}$, $\widehat{\mathit{\lambda}}$, and $\widehat{\stackrel{\mathrm{\u0303}}{P}}$ of the Fcentric grid (at Fpoint) 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 σ_{12} 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 Fcentric 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 Fcentric stress tensors at the dynamical time step level, is hereafter referred to as crossnudging (CN). It consists of nudging each vertically integrated component of the Tcentric stress tensor (σ) towards its Fcentric 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:
where γ_{cn} is the CN coefficient, N_{s} is the timesplitting parameter (Eq. 15), the bar notation denotes the spatial interpolation from F to Tpoints or T to Fpoints (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 crossnudging 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 σ^{(i)} 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 postCN stresses and the damage field.
The fourpoint 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 Fcentric solutions. We have performed sensitivity tests with our panArctic 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 straightline 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 finiteelement 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 icevelocitydriven 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 timesplitting instance.
We use the secondorder 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 Tpoint (d, σ_{11}, σ_{22}, and ${\widehat{\mathit{\sigma}}}_{\mathrm{12}}$) are advected using U and V defined at U and Vpoints, respectively. Their Fpoint counterparts ($\widehat{d}$, ${\widehat{\mathit{\sigma}}}_{\mathrm{11}}$, ${\widehat{\mathit{\sigma}}}_{\mathrm{22}}$, and σ_{12}) are advected using $\widehat{U}$ and $\widehat{V}$ defined at V and Upoints, respectively. In practice, the exact same implementation of the advection scheme can be used to perform the advection at T and Fpoints; the only difference is that for the advection of Fpoint 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 secondrank 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}}$:
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 upperconvected time derivative of σ, denoted as $\stackrel{\u25bf}{\mathit{\sigma}}$,
which, in component form, simplifies into
The second formulation is the lowerconvected time derivative, denoted as $\stackrel{\u25b3}{\mathit{\sigma}}$,
with
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 Fpoints with the Egrid. 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 tensorspecific 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 upperconvected formulation in both the idealized and panArctic 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 upperconvected 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 lowerconvected 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 lowerconvected 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 timederivative formulation (or combination of them, such as the Gordon–Schowalter timederivative 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 panArctic 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 SI3derived 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 earliestoccurring 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 lineintegral 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 Cgrid velocities u and v of the Tcentric cell are used to track the points ($\widehat{u}$ and $\widehat{v}$, available in the BBMdriven simulation, are not used).
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 panArctic 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 aEVPdriven 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 rheologyrelated 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 standalone mode using SAS, the standalone 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 panArctic simulations (Sect. 3.2).
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 “spiderweblike” arrangement of the LKFs in the aEVP solution. In the BBMdriven 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).
Finally, we note that the solutions do not contain any apparent numerical instabilities or noise for aEVP or BBM. The LKFlike 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 panArctic simulations
The panArctic simulations use SI3 coupled to the 3D ocean component of NEMO, named OCE. They are performed on the socalled 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 zcoordinate 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 nearsurface 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).
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 SI3default, uses the default SI3 setup as provided in NEMO and thereby uses the aEVP rheology of Kimmritz et al. (2016). The second simulation, named SI3BBM, only differs from SI3default 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 2month spinup performed with the SI3default setup and run until 20 April 1997. This spinup 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 2month spinup 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 SI3BBM 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 SI3BBM 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^{−3} to 2 × 10^{−3} in SI3BBM 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 SI3default, the default value of 1.4 × 10^{−3} satisfies this requirement and is left unchanged. For the timesplitting approach (Sect. 2.2), we use a dynamical time step of 7.2 s in SI3BBM, which relates to a timesplitting parameter of N_{s}=100.
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 SI3BBM 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 SI3default exhibits very smooth fields of deformation with almost no such localized features; this is consistent with the findings of recent studies that evaluate VPdriven sea ice simulations run at a grid resolution typically coarser than 5 km (e.g., Ólason et al., 2022; Bouchat et al., 2022).
The PDFs of the total deformation rates (Fig. 8d) show that SI3BBM 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 SI3BBM 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, SI3BBM 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 SI3default 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, SI3default exhibits a powerlaw behavior similar to that of both observations and SI3BBM, 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 midDecember 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.
As illustrated in Fig. 9, SI3BBM 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 SI3BBM 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.
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 Multifractal 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 multifractal 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 SI3BBM 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 multifractal 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 SI3BBM across the full range of spatial scales considered for this analysis and can be approximated by a powerlaw 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^{−3} in SI3BBM, while the ${C}_{\mathrm{D}}^{\left(\mathrm{a}\right)}$ = 1.4 × 10^{−3} used in SI3default 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 SI3BBM compared to that derived from the observations. Indeed, the exponents of the power law that fits the SI3BBM 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 SI3BBM 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 SI3default 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 powerlaw slopes that are in better agreement with those derived from observations. We argue that since sea ice deformation is a scaleinvariant 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.
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 SI3BBM are multifractal because their structure functions is well approximated (in the sense of the leastsquare 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 multifractal 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 multifractality. The structure function of SI3BBM 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).
4.1 On the numerical implementation
The crossnudging has a noteworthy analogy to the Asselin filter (Asselin, 1972) used when discretizing time derivatives of a prognostic variable by means of the leapfrog scheme (three time levels, centered, and secondorder), in particular in the context of shallowwater 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 crossnudging can be seen as a sort of spatial and twodimensional 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 CMIPclass 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 crossnudging 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 Eaugmented Cgrid, 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 crossnudging 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 finitedifference 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 largescale, longterm sea ice simulations. This is the scope of our present work and future papers.
As discussed in Sect. 2.3.1, the use of the Egrid 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 Fcentric 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 walltime 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 standalone mode (SAS). In standard coupled mode, with SI3 coupled to OCE, the BBMrelated 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 SI3default and SI3BBM, 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 largescale 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 subgridscale 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 BBMspecific 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 multifractality of the deformation fields simulated by SI3BBM 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 crossnudging 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 SI3BBM reproduces the mean across all scales remarkably well and why the powerlaw exponents for the variance and the skewness are underestimated. In this regard, the use of the finiteelement method together with the discontinuous Galerkin method might prove to be a promising combination to simulate the multifractality of sea ice deformation even more accurately while remaining in the Eulerian and quadrilateral mesh framework.
The brittle Bingham–Maxwell rheology, known as BBM, has been successfully implemented into SI3, the CMIPclass, Eulerian finitedifference 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 panArctic scale when compared to satellite observations.
Our implementation uses a new discretization approach that expands the Cgrid of NEMO into the Arakawa Egrid 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 Egrid, the issue of the grid separation is introduced. We have introduced a simple technique to prevent this grid separation in the form of the crossnudging. This crossnudging 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 Eaugmented Cgrid 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 upperconvected time derivative rather than its lowerconvected 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 panArctic 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 nonGaussian 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 BBMdriven 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 multifractal nature, is fairly well captured by the BBMdriven simulation but with a slightly lower degree of multifractality.
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)  – 
E_{0},E  elastic modulus of undamaged and damaged ice  Pa 
λ_{0},λ  apparent viscous relaxation time of undamagedand damaged ice  s 
$\stackrel{\mathrm{\u0303}}{P}$  BBMspecific ridging term  – 
P_{max}  ridging threshold  Pa 
P_{0}  scaling parameter for P_{max}  Pa 
h_{0}  reference ice thickness for P_{max}  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 
C_{E}  propagation speed of an elasticshear wave  m s^{−1} 
t_{d}  characteristic timescale for thepropagation of damage  s 
d_{crit}  damage increment (Dansereau, 2016)  – 
k_{th}  healing constant for damage  K s 
ΔT_{h}  temperature difference betweenbottom and surface of ice  K 
A2 Acronyms
NEMO  Nucleus for European Modeling of the Ocean 
SI3  SeaIce modeling Integrated Initiative (sea icecomponent of NEMO) 
OCE  3D ocean component of NEMO 
SAS  Standalone surface module of NEMO (i.e.,SI3 standalone) 
LIM  LouvainLaNeuve seaIce Model 
BBM  Brittle Bingham–Maxwell (rheology) 
MEB  Maxwell elastobrittle (rheology) 
VP  Viscous–plastic (rheology) 
FD  Finite difference (method) 
CN  Crossnudging (treatment) 
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 Egrid
The bar + superscript notation refers to a spatial interpolation; $\stackrel{X}{\mathit{\varphi}}$ is field ϕ interpolated onto X points.
Interpolation from T to Fpoints, or conversely, is based on the average of the four nearest surrounding points (Fig. 3a).
Note: surrounding points located on land or openboundary 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.
The “hat” notation $\widehat{x}$ refers to the Fcentric counterpart of x, with x being a prognostic scalar or tensor (rank 1 or 2) defined in the Tcentric grid (mind that if x is the element of a tensor, $\widehat{x}$ is not necessarily defined on Fpoints). Example: $\widehat{d}$ and ${\widehat{\mathit{\sigma}}}_{\mathrm{11}}$ are prognostic fields defined on Fpoints (the natural location for d and σ_{11} on the Cgrid is the Tpoint); similarly, ${\widehat{\mathit{\sigma}}}_{\mathrm{12}}$ is defined on Tpoints (the natural location for σ_{12} on the Cgrid is the Fpoint).
A4 Miscellaneous notations
$\underset{\mathrm{\xaf}}{\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{\u203e}}{h}}^{F}{\mathit{\sigma}}_{\text{kl}}$ if σ_{kl} defined at F  
$\stackrel{\u25bf}{\mathit{x}}$  upperconvected time derivative of symmetric(rank 2) tensor x 
$\stackrel{\u25b3}{\mathit{x}}$  lowerconvected 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 
N_{s}  $\equiv \mathrm{\Delta}T/\mathrm{\Delta}t$, timesplittingparameter  – 
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 Tcentric 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 Fcentric 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 Tcentric 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 Fcentric 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 
u_{o},v_{o}  surface ocean current at Uand at V  m s^{−1} 
γ_{cn}  crossnudging coefficient  – 
${C}_{\mathrm{D}}^{\left(\mathrm{a}\right)}$  ice–atmosphere drag coefficient  – 
Δ^{𝚃}x  Tcentered Δx that connects two neighboringUpoints  m 
Δ^{𝚃}y  Tcentered Δy that connects two neighboringVpoints  m 
Δ^{𝙵}x  Fcentered Δx that connects two neighboringVpoints  m 
Δ^{𝙵}y  Fcentered Δy that connects two neighboringUpoints  m 
Δ^{𝚄}x  Ucentered Δx that connects two neighboringTpoints  m 
Δ^{𝚄}y  Ucentered Δy that connects two neighboringFpoints  m 
Δ^{𝚅}x  Vcentered Δx that connects two neighboringFpoints  m 
Δ^{𝚅}y  Vcentered Δy that connects two neighboringTpoints  m 
B1 Algorithm
Timesplitting 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{\u0303}}{P},\widehat{\stackrel{\mathrm{\u0303}}{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 crossnudging 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}}$.

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 (timesplitting loop above)

Advection of generic SI3 prognostic tracers (A, h, etc.) at Tpoints using U and V

Advection of d, σ_{11}, σ_{22}, and ${\widehat{\mathit{\sigma}}}_{\mathrm{12}}$ at Tpoints using U and V

Advection of $\widehat{d}$, ${\widehat{\mathit{\sigma}}}_{\mathrm{11}}$, ${\widehat{\mathit{\sigma}}}_{\mathrm{22}}$, and σ_{12} at Fpoints 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 Fcentric 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 Fcentric grids based on the finitedifference 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.
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 Fpoints, $\stackrel{\mathrm{F}}{A}$, that is used in the equations for the Fcentric 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 Fpoints, $\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{\u0303}}{P}$ term
$$\begin{array}{}\text{(B9)}& \begin{array}{rl}& \stackrel{\mathrm{\u0303}}{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{\u0303}}{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{\u0303}}{P})\mathrm{\Delta}t}}\\ \widehat{\mathrm{\Omega}}& ={\displaystyle \frac{\widehat{\mathit{\lambda}}}{\widehat{\mathit{\lambda}}+(\mathrm{1}+\widehat{\stackrel{\mathrm{\u0303}}{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})\phantom{\rule{0.25em}{0ex}}\mathrm{\Delta}t/{t}_{\mathrm{d}}\end{array}\\ {\underset{\mathrm{\xaf}}{\mathit{\sigma}}}^{k+\mathrm{1}}=\begin{array}{c}{\underset{\mathrm{\xaf}}{\mathit{\sigma}}}^{\left(i\right)k+\mathrm{1}}\left(\mathrm{1}{d}_{\text{crit}}\right)\\ \times {\underset{\mathrm{\xaf}}{\mathit{\sigma}}}^{\left(i\right)k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}\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})\phantom{\rule{0.25em}{0ex}}\mathrm{\Delta}t/{\widehat{t}}_{d}\end{array}\\ {\widehat{\underset{\mathrm{\xaf}}{\mathit{\sigma}}}}^{k+\mathrm{1}}=\begin{array}{c}{\widehat{\underset{\mathrm{\xaf}}{\mathit{\sigma}}}}^{\left(i\right)k+\mathrm{1}}\left(\mathrm{1}{\widehat{d}}_{\text{crit}}\right)\\ \times {\widehat{\underset{\mathrm{\xaf}}{\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{\xaf}}{\mathit{\sigma}}}^{k+\mathrm{1}}={\underset{\mathrm{\xaf}}{\mathit{\sigma}}}^{\left(i\right)k+\mathrm{1}}\\ & {\widehat{d}}^{k+\mathrm{1}}={\widehat{d}}^{k}\\ & {\widehat{\underset{\mathrm{\xaf}}{\mathit{\sigma}}}}^{k+\mathrm{1}}={\widehat{\underset{\mathrm{\xaf}}{\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 Fcentric 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\phantom{\rule{0.25em}{0ex}}{\mathit{\sigma}}_{\mathrm{11}}\right)}{\partial x}+\frac{\partial \left(h\phantom{\rule{0.25em}{0ex}}{\mathit{\sigma}}_{\mathrm{12}}\right)}{\partial y}\\ \frac{\partial \left(h\phantom{\rule{0.25em}{0ex}}{\mathit{\sigma}}_{\mathrm{22}}\right)}{\partial y}+\frac{\partial \left(h\phantom{\rule{0.25em}{0ex}}{\mathit{\sigma}}_{\mathrm{12}}\right)}{\partial x}\end{array}\right)\end{array}$$ 
Discretized in the Tcentric 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\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{T}}{y}^{\mathrm{2}}\right]}_{i+\mathrm{1},j}{\left[{\mathit{\sigma}}_{\mathrm{11}}^{k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}h\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{T}}{y}^{\mathrm{2}}\right]}_{i,j}}{[{\mathrm{\Delta}}^{\mathtt{U}}x\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{U}}{y}^{\mathrm{2}}{]}_{i,j}}}\\ & +{\displaystyle \frac{{\left[{\mathit{\sigma}}_{\mathrm{12}}^{k+\mathrm{1}}\stackrel{\mathrm{F}}{h}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{F}}{x}^{\mathrm{2}}\right]}_{i,j}{\left[{\mathit{\sigma}}_{\mathrm{12}}^{k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}\stackrel{\mathrm{F}}{h}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{F}}{x}^{\mathrm{2}}\right]}_{i,j\mathrm{1}}}{[{\mathrm{\Delta}}^{\mathtt{U}}y\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{U}}{x}^{\mathrm{2}}{]}_{i,j}}}\text{(atU)}\\ & {\mathrm{\Delta}}_{y}^{k+\mathrm{1}}{}_{i,j}={\displaystyle \frac{{\left[{\mathit{\sigma}}_{\mathrm{22}}^{k+\mathrm{1}}h\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{T}}{x}^{\mathrm{2}}\right]}_{i,j+\mathrm{1}}{\left[{\mathit{\sigma}}_{\mathrm{22}}^{k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}h\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{T}}{x}^{\mathrm{2}}\right]}_{i,j}}{[{\mathrm{\Delta}}^{\mathtt{V}}y\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{V}}{x}^{\mathrm{2}}{]}_{i,j}}}\\ & +{\displaystyle \frac{{\left[{\mathit{\sigma}}_{\mathrm{12}}^{k+\mathrm{1}}\widehat{h}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{F}}{y}^{\mathrm{2}}\right]}_{i,j}{\left[{\mathit{\sigma}}_{\mathrm{12}}^{k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}\widehat{h}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{F}}{y}^{\mathrm{2}}\right]}_{i\mathrm{1},j}}{[{\mathrm{\Delta}}^{\mathtt{V}}x\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{V}}{y}^{\mathrm{2}}{]}_{i,j}}}\text{(atV)}\end{array}\end{array}$$ 
Discretized in the Fcentric 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}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{F}}{y}^{\mathrm{2}}\right]}_{i,j}{\left[{\widehat{\mathit{\sigma}}}_{\mathrm{11}}^{k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}\stackrel{\mathrm{F}}{h}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{F}}{y}^{\mathrm{2}}\right]}_{i\mathrm{1},j}}{[{\mathrm{\Delta}}^{\mathtt{V}}x\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{V}}{y}^{\mathrm{2}}{]}_{i,j}}}\\ & +{\displaystyle \frac{{\left[{\widehat{\mathit{\sigma}}}_{\mathrm{12}}^{k+\mathrm{1}}h\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{T}}{x}^{\mathrm{2}}\right]}_{i,j+\mathrm{1}}{\left[{\widehat{\mathit{\sigma}}}_{\mathrm{12}}^{k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}h\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{T}}{x}^{\mathrm{2}}\right]}_{i,j}}{[{\mathrm{\Delta}}^{\mathtt{V}}y\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{V}}{x}^{\mathrm{2}}{]}_{i,j}}}\text{(atV)}\\ & {\widehat{\mathrm{\Delta}}}_{y}^{k+\mathrm{1}}{}_{i,j}={\displaystyle \frac{{\left[{\widehat{\mathit{\sigma}}}_{\mathrm{22}}^{k+\mathrm{1}}\stackrel{\mathrm{F}}{h}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{F}}{x}^{\mathrm{2}}\right]}_{i,j}{\left[{\widehat{\mathit{\sigma}}}_{\mathrm{22}}^{k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}\stackrel{\mathrm{F}}{h}\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{F}}{x}^{\mathrm{2}}\right]}_{i,j\mathrm{1}}}{[{\mathrm{\Delta}}^{\mathtt{U}}y\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{U}}{x}^{\mathrm{2}}{]}_{i,j}}}\\ & +{\displaystyle \frac{{\left[{\widehat{\mathit{\sigma}}}_{\mathrm{12}}^{k+\mathrm{1}}h\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{T}}{y}^{\mathrm{2}}\right]}_{i+\mathrm{1},j}{\left[{\widehat{\mathit{\sigma}}}_{\mathrm{12}}^{k+\mathrm{1}}\phantom{\rule{0.25em}{0ex}}h\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{T}}{y}^{\mathrm{2}}\right]}_{i,j}}{[{\mathrm{\Delta}}^{\mathtt{U}}x\phantom{\rule{0.25em}{0ex}}{\mathrm{\Delta}}^{\mathtt{U}}{y}^{\mathrm{2}}{]}_{i,j}}}\text{(atU)}\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 Egrid 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}).
Then, the discretized equation of momentum yields the expression of the two velocity components at time level k+1.

The NEMO source code used to perform the experiments is based on the official release 4.2.2 of NEMO; it is available on Zenodo with the DOI https://doi.org/10.5281/zenodo.11551599 (The NEMO team, 2024).

New and modified Fortran90 source files relative to our implementation of the BBM rheology in version 4.2.2 of NEMO/SI3 are available on Zenodo with the DOI https://doi.org/10.5281/zenodo.11581840 (Brodeau, 2024a).

The Python software used to seed and build Lagrangian trajectories out of the SI3 hourly sea ice velocities is named
sitrack
; the version used to perform the present study is available on Zenodo with the DOI https://doi.org/10.5281/zenodo.10457918 (Brodeau, 2024b). 
The Python software used to compute the RGPS and modelbased sea ice deformation rates based on quadrangles and perform the scaling analysis is named
mojito
; the version used to perform the present study is available on Zenodo with the DOI https://doi.org/10.5281/zenodo.10457924 (Brodeau, 2924c). 
Model data produced and analyzed in this study, namely SI3 hourly output files for panArctic simulations SI3BBM and SI3default, are available on Zenodo with the DOI https://doi.org/10.5281/zenodo.11582103 (Brodeau, 2024d).

Model data as well as setup, forcing, and configuration files for the “cyclone” test case experiments (Mehlmann et al., 2021) discussed and analyzed in this study are available on Zenodo with the DOI https://doi.org/10.5281/zenodo.11615982 (Brodeau, 2024e).
This study is based on the original idea of PR and EÓ, who suggested using the Egrid. The implementation of the BBM rheology into SI3 has been carried out by LB, with the help of PR. The numerical simulations and statistical analyses have been performed by LB. LB also prepared all the figures. Interpretations of the results were provided by LB, PR, and EÓ. LB and PR led the writing of the manuscript, incorporating input from EÓ. Suggestions and improvements on the text have also been provided by VD.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We would like to express our gratitude to one anonymous reviewer and Mathieu Plante for conducting a thorough and constructive review. Their comments and suggestions have greatly contributed to improving and clarifying this manuscript.
This research received support through Schmidt Sciences, LLC, via the SASIP project (grant G2466154).
This paper was edited by Riccardo Farneti and reviewed by Mathieu Plante and two anonymous referees.
Arakawa, A.: Design of the UCLA general circulation model, Tech. Rept. No. 7, 116 pp., University of California at Los Angeles, CA, 90024, 1972. a
Arakawa, A. and Lamb, V. R.: Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, in: Methods in Computational Physics: Advances in Research and Applications, edited by: Chang, J., Elsevier, https://doi.org/10.1016/b9780124608177.500094, ISBN 9780124608177; ISSN: 00766860, pp. 173–265, 1977. a, b, c
Asselin, R.: Frequency Filter for Time Integrations, Mon. Weather Rev., 100, 487–490, https://doi.org/10.1175/15200493(1972)100<0487:fffti>2.3.co;2, 1972. a
Barnier, B., Madec, G., Penduff, T., Molines, J.M., Treguier, A.M., Sommer, J. L., Beckmann, A., Biastoch, A., Böning, C., Dengg, J., Derval, C., Durand, E., Gulev, S., Remy, E., Talandier, C., Theetten, S., Maltrud, M., McClean, J., and Cuevas, B. D.: Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddypermitting resolution, Ocean Dynam., 56, 543–567, 2006. a
Bouchat, A. and Tremblay, B.: Using seaice deformation fields to constrain the mechanical strength parameters of geophysical sea ice, J. Geophys. Res.Oceans, 122, 5802–5825, https://doi.org/10.1002/2017jc013020, 2017. a
Bouchat, A., Hutter, N., Chanut, J., Dupont, F., Dukhovskoy, D., Garric, G., Lee, Y., Lemieux, J.F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Ólason, E., Rampal, P., Rasmussen, T., Talandier, C., Tremblay, B., and Wang, Q.: Sea Ice Rheology Experiment (SIREx), part I: Scaling and statistical properties of seaice deformation fields, J. Geophys. Res.Oceans, 127, e2021JC017667, https://doi.org/10.1029/2021jc017667, 2022. a, b, c, d
Bouillon, S. and Rampal, P.: Presentation of the dynamical core of neXtSIM, a new sea ice model, Ocean Model., 91, 23–37, https://doi.org/10.1016/j.ocemod.2015.04.005, 2015. a
Bouillon, S., Maqueda, M. Á. M., Legat, V., and Fichefet, T.: An elastic–viscous–plastic sea ice model formulated on Arakawa B and C grids, Ocean Model., 27, 174–184, https://doi.org/10.1016/j.ocemod.2009.01.004, 2009. a
Boutin, G., Ólason, E., Rampal, P., Regan, H., Lique, C., Talandier, C., Brodeau, L., and Ricker, R.: Arctic sea ice mass balance in a new coupled ice–ocean model using a brittle rheology framework, The Cryosphere, 17, 617–638, https://doi.org/10.5194/tc176172023, 2023. a, b
Brodeau, L.: brodeau/bbm4si3: 1.1.0 (1.1.0), Zenodo [code], https://doi.org/10.5281/zenodo.11581840, 2024a. a
Brodeau, L.: brodeau/sitrack: 1.0 (1.0), Zenodo [code], https://doi.org/10.5281/zenodo.10457918, 2024b. a
Brodeau, L.: brodeau/mojito: 1.0 (1.0), Zenodo [code], https://doi.org/10.5281/zenodo.10457924, 2024. a
Brodeau, L.: SI3 seaice hourly output for SI3BBM and SI3default PanArctic simulations (Brodeau et al., 2024, final version), Zenodo [data set], https://doi.org/10.5281/zenodo.11582103, 2024d. a
Brodeau, L.: NEMO/SASSI3 setup, forcing and output for SI3BBM and SI3default idealized “cyclone” simulations (Brodeau et al., 2024, final version), Zenodo [data set], https://doi.org/10.5281/zenodo.11615982, 2024e. a
Coon, M. D., Maykut, S. A., Pritchard, R. S., Rothrock, D. A., and Thorndike, A. S.: Modeling the pack ice as an elasticplastic material, AIDJEX Bull., 24, 1–105, 1974. a
Danilov, S., Mehlmann, C., and Fofonova, V.: On discretizing seaice dynamics on triangular meshes using vertex, cell or edge velocities, Ocean Model., 170, 101937, https://doi.org/10.1016/j.ocemod.2021.101937, 2022. a
Danilov, S., Mehlmann, C., Sidorenko, D., and Wang, Q.: CDtype discretization for sea ice dynamics in FESOM version 2, Geosci. Model Dev., 17, 2287–2297, https://doi.org/10.5194/gmd1722872024, 2024. a
Dansereau, V.: A MaxwellElastoBrittle model for the drift and deformation of sea ice, PhD thesis, Université Grenoble Alpes, Grenoble, France, 2016. a, b
Dansereau, V., Weiss, J., Saramito, P., and Lattes, P.: A Maxwell elastobrittle rheology for sea ice modelling, The Cryosphere, 10, 1339–1359, https://doi.org/10.5194/tc1013392016, 2016. a, b, c, d, e, f, g
Dansereau, V., Weiss, J., Saramito, P., Lattes, P., and Coche, E.: Ice bridges and ridges in the MaxwellEB sea ice rheology, The Cryosphere, 11, 2033–2058, https://doi.org/10.5194/tc1120332017, 2017. a
Ferry, N., Barnier, B., Garric, G., Haines, K., Masina, S., Parent, L., Storto, A., Valdivieso, M., Guinehut, S., and Mulet, S.: NEMO: the modelling engine of global ocean reanalyses, Mercator Ocean Quarterly Newsletter, 46, 46–59, 2012. a
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., MuñozSabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049, https://doi.org/10.1002/qj.3803, 2020. a
Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, https://doi.org/10.1175/15200485(1979)009<0815:adtsim>2.0.co;2, 1979. a
Hinch, J. and Harlen, O.: Oldroyd B, and not A?, J. NonNewton. Fluid, 298, 104668, https://doi.org/10.1016/j.jnnfm.2021.104668, 2021. a, b
Hopkins, M. A.: Four stages of pressure ridging, J. Geophys. Res.Oceans, 103, 21883–21891, https://doi.org/10.1029/98jc01257, 1998. a
Hunke, E. C. and Dukowicz, J. K.: The Elastic–Viscous–Plastic Sea Ice Dynamics Model in General Orthogonal Curvilinear Coordinates on a Sphere—Incorporation of Metric Terms, Mon. Weather Rev., 130, 1848–1865, https://doi.org/10.1175/15200493(2002)130<1848:tevpsi>2.0.co;2, 2002. a
Hutter, N., Losch, M., and Menemenlis, D.: Scaling Properties of Arctic Sea Ice Deformation in a HighResolution ViscousPlastic Sea Ice Model and in Satellite Observations, J. Geophys. Res.Oceans, 123, 672–687, https://doi.org/10.1002/2017jc013119, 2018. a, b
Hutter, N., Bouchat, A., Dupont, F., Dukhovskoy, D., Koldunov, N., Lee, Y. J., Lemieux, J.F., Lique, C., Losch, M., Maslowski, W., Myers, P. G., Ólason, E., Rampal, P., Rasmussen, T., Talandier, C., Tremblay, B., and Wang, Q.: Sea Ice Rheology Experiment (SIREx): 2. Evaluating Linear Kinematic Features in HighResolution Sea Ice Simulations, J. Geophys. Res.Oceans, 127, e2021JC017666, https://doi.org/10.1029/2021jc017666, 2022. a, b
IPCC: The Ocean and Cryosphere in a Changing Climate: Special Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK, https://doi.org/10.1017/9781009157964.005, pp. 203–320, 2022. a
Janjić, Z. I.: A Stable Centered Difference Scheme Free of TwoGridInterval Noise, Mon. Weather Rev., 102, 319–323, https://doi.org/10.1175/15200493(1974)102<0319:ASCDSF>2.0.CO;2, 1974. a, b, c
Janjić, Z. I.: Nonlinear Advection Schemes and Energy Cascade on SemiStaggered Grids, Mon. Weather Rev., 112, 1234–1245, https://doi.org/10.1175/15200493(1984)112<1234:nasaec>2.0.co;2, 1984. a
Janjić, Z. I. and Mesinger, F.: Finite difference methods for the shallow water equations on various horizontal grids, in: Numerical Methods for Weather Prediction, vol. 1, Seminar 1983, ECMWF, Shinfield Park, Reading, p. 29–101, 1984. a
Kagan, Y. Y. and Jackson, D. D.: LongTerm Earthquake Clustering, Geophys. J. Int., 104, 117–134, https://doi.org/10.1111/j.1365246x.1991.tb02498.x, 1991. a
Kimmritz, M., Danilov, S., and Losch, M.: The adaptive EVP method for solving the sea ice momentum equation, Ocean Model., 101, 59–67, https://doi.org/10.1016/j.ocemod.2016.03.004, 2016. a, b, c, d
Kimmritz, M., Losch, M., and Danilov, S.: A comparison of viscousplastic sea ice solvers with and without replacement pressure, Ocean Model., 115, 59–69, https://doi.org/10.1016/j.ocemod.2017.05.006, 2017. a
Konor, C. S. and Randall, D. A.: Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations – Part 1: Nonhydrostatic inertia–gravity modes, Geosci. Model Dev., 11, 1753–1784, https://doi.org/10.5194/gmd1117532018, 2018. a
Kwok, R.: Deformation of the Arctic Ocean Sea Ice Cover between November 1996 and April 1997: A Qualitative Survey, in: IUTAM Symposium on Scaling Laws in Ice Mechanics and Ice Dynamics, edited by: Dempsey, J. P. and Shen, H. H., Springer Netherlands, Dordrecht, 315–322, ISBN 9789401597357, https://doi.org/10.1007/9789401597357_26, p. 315–322, 2001. a
Kwok, R., Schweiger, A., Rothrock, D. A., Pang, S., and Kottmeier, C.: Sea ice motion from satellite passive microwave imagery assessed with ERS SAR and buoy motions, J. Geophys. Res.Oceans, 103, 8191–8214, https://doi.org/10.1029/97jc03334, 1998. a, b
Larson, R. G.: Constitutive Equations for Polymer Melts and Solutions, Elsevier, https://doi.org/10.1016/c20130042843, ISBN 9780409901191, 1988. a
Lemieux, J., Tremblay, B., Thomas, S., Sedláček, J., and Mysak, L. A.: Using the preconditioned Generalized Minimum RESidual (GMRES) method to solve the seaice momentum equation, J. Geophys. Res.Oceans (1978–2012), 113, C10004, https://doi.org/10.1029/2007jc004680, 2008. a
Lemieux, J.F., Knoll, D. A., Losch, M., and Girard, C.: A secondorder accurate in time IMplicit–EXplicit (IMEX) integration scheme for sea ice dynamics, J. Comput. Phys., 263, 375–392, https://doi.org/10.1016/j.jcp.2014.01.010, 2014. a
Lindsay, R. W. and Stern, H. L.: The RADARSAT Geophysical Processor System: Quality of Sea Ice Trajectory and Deformation Estimates, J. Atmos. Ocean. Tech., 20, 1333–1347, https://doi.org/10.1175/15200426(2003)020<1333:trgpsq>2.0.co;2, 2003. a
Losch, M., Menemenlis, D., Campin, J.M., Heimbach, P., and Hill, C.: On the formulation of seaice models. Part 1: Effects of different solver implementations and parameterizations, Ocean Model., 33, 129–144, https://doi.org/10.1016/j.ocemod.2009.12.008, 2010. a, b
Lovejoy, S. and Schertzer, D.: Nonlinear Dynamics in Geosciences, vol. 59 of Nonlinear Dynamics in Geosciences, Springer New York, New York, NY, https://doi.org/10.1007/9780387349183_18, pp. 311–337,2007. a
Madec, G., BourdalléBadie, R., Chanut, J., Clementi, E., Coward, A., Ethé, C., Iovino, D., Lea, D., Lévy, C., Lovato, T., Martin, N., Masson, S., Mocavero, S., Rousset, C., Storkey, D., Müeller, S., Nurser, G., Bell, M., Samson, G., Mathiot, P., Mele, F., and Moulin, A.: NEMO ocean engine, Zenodo, https://doi.org/10.5281/zenodo.1464816, 2022. a
MaierReimer, E., Mikolajewicz, U., and Hasselmann, K.: Mean Circulation of the Hamburg LSG OGCM and Its Sensitivity to the Thermohaline Surface Forcing, J. Phys. Oceanogr., 23, 731–757, https://doi.org/10.1175/15200485(1993)023<0731:mcothl>2.0.co;2, 1993. a, b
Marsaleix, P., Auclair, F., Duhaut, T., Estournel, C., Nguyen, C., and Ulses, C.: Alternatives to the Robert–Asselin filter, Ocean Model., 41, 53–66, https://doi.org/10.1016/j.ocemod.2011.11.002, 2012. a
Marsan, D. and Weiss, J.: Space/time coupling in brittle deformation at geophysical scales, Earth Planet. Sc. Lett., 296, 353–359, https://doi.org/10.1016/j.epsl.2010.05.019, 2010. a
Marsan, D., Stern, H., Lindsay, R., and Weiss, J.: Scale dependence and localization of the deformation of arctic sea ice, Phys. Rev. Lett., 93, 178501, https://doi.org/10.1103/physrevlett.93.178501, 2004. a
Mehlmann, C., Danilov, S., Losch, M., Lemieux, J. F., Hutter, N., Richter, T., Blain, P., Hunke, E. C., and Korn, P.: Simulating Linear Kinematic Features in ViscousPlastic Sea Ice Models on Quadrilateral and Triangular Grids With Different Variable Staggering, J. Adv. Model. Earth Sy., 13, e2021MS002523, https://doi.org/10.1029/2021ms002523, 2021. a, b, c, d, e, f, g, h
Mesinger, F.: A method for construction of secondorder accuracy difference schemes permitting no false twogridinterval wave in the height field, Tellus, 25, 444–458, https://doi.org/10.1111/j.21533490.1973.tb00629.x, 1973. a, b, c
Mesinger, F. and Popovic, J.: Forward–backward scheme on the B/E grid modified to suppress lattice separation: the two versions, and any impact of the choice made?, Meteorol. Atmos. Phys., 108, 1–8, https://doi.org/10.1007/s0070301000801, 2010. a, b
Ólason, E., Boutin, G., Korosov, A., Rampal, P., Williams, T., Kimmritz, M., Dansereau, V., and Samaké, A.: A new brittle rheology and numerical framework for largescale seaice models, J. Adv. Model. Earth Sy., 14, e2021MS002685, https://doi.org/10.1029/2021MS002685, 2022. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s
Oldroyd, J. G.: On the formulation of rheological equations of state, P. Roy. Soc. Lond. A Mat., 200, 523–541, https://doi.org/10.1098/rspa.1950.0035, 1950. a
Plante, M. and Tremblay, L. B.: A generalized stress correction scheme for the Maxwell elastobrittle rheology: impact on the fracture angles and deformations, The Cryosphere, 15, 5623–5638, https://doi.org/10.5194/tc1556232021, 2021. a
Plante, M., Tremblay, B., Losch, M., and Lemieux, J.F.: Landfast sea ice material properties derived from ice bridge simulations using the Maxwell elastobrittle rheology, The Cryosphere, 14, 2137–2157, https://doi.org/10.5194/tc1421372020, 2020. a, b, c, d
Prather, M. J.: Numerical advection by conservation of secondorder moments, J. Geophys. Res., 91, 6671, https://doi.org/10.1029/jd091id06p06671, 1986. a, b
Rampal, P., Weiss, J., Marsan, D., Lindsay, R., and Stern, H.: Scaling properties of sea ice deformation from buoy dispersion analysis, J. Geophys. Res., 113, C3, https://doi.org/10.1029/2007jc004143, 2008. a, b
Rampal, P., Bouillon, S., Ólason, E., and Morlighem, M.: neXtSIM: a new Lagrangian sea ice model, The Cryosphere, 10, 1055–1073, https://doi.org/10.5194/tc1010552016, 2016. a, b, c, d, e
Rampal, P., Dansereau, V., Ólason, E., Bouillon, S., Williams, T., Korosov, A., and Samaké, A.: On the multifractal scaling properties of sea ice deformation, The Cryosphere, 13, 2457–2474, https://doi.org/10.5194/tc1324572019, 2019. a, b, c
Regan, H., Rampal, P., Ólason, E., Boutin, G., and Korosov, A.: Modelling the evolution of Arctic multiyear sea ice over 20002018, The Cryosphere, https://doi.org/10.5194/egusphereegu2315844, 2023. a
Rheinlænder, J. W., Davy, R., Ólason, E., Rampal, P., Spensberger, C., Williams, T. D., Korosov, A., and Spengler, T.: Driving Mechanisms of an Extreme Winter Sea Ice Breakup Event in the Beaufort Sea, Geophys. Res. Lett., 49, e2022GL099024, https://doi.org/10.1029/2022gl099024, 2022. a
Rousset, C., Vancoppenolle, M., Madec, G., Fichefet, T., Flavoni, S., Barthélemy, A., Benshila, R., Chanut, J., Levy, C., Masson, S., and Vivier, F.: The LouvainLaNeuve sea ice model LIM3.6: global and regional capabilities, Geosci. Model Dev., 8, 2991–3005, https://doi.org/10.5194/gmd829912015, 2015. a
Samaké, A., Rampal, P., Bouillon, S., and Ólason, E.: Parallel implementation of a Lagrangianbased model on an adaptive mesh in C++: Application to seaice, J. Comput. Phys., 350, 84–96, https://doi.org/10.1016/j.jcp.2017.08.055, 2017. a
Savage, L. J.: The Foundations of Statistics, Wiley publications in statistics, Wiley, 294 pp., https://books.google.fr/books?id=wqzV4GoYMgEC (last access: 5 August 2024), 1954. a
Schertzer, D. and Lovejoy, S.: Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes, J. Geophys. Res.Atmos., 92, 9693–9714, https://doi.org/10.1029/jd092id08p09693, 1987. a
Snoeijer, J. H., Pandey, A., Herrada, M. A., and Eggers, J.: The relationship between viscoelasticity and elasticity, P. Roy. Soc. AMath Phy., 476, 2243, https://doi.org/10.1098/rspa.2020.0419, 2020. a
Sornette, D.: Power Law Distributions, Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/3540331824_4, pp. 93–121, 2006. a
Stern, H. L., Schweiger, A. J., Stark, M., Zhang, J., Steele, M., and Hwang, B.: Seasonal evolution of the seaice floe size distribution in the Beaufort and Chukchi seas, Elementa: Science of the Anthropocene, 6, 48, https://doi.org/10.1525/elementa.305, 2018. a
Stone, H. A., Shelley, M. J., and Boyko, E.: A note about convected time derivatives for flows of complex fluids, Soft Matter, 19, 5353–5359, https://doi.org/10.1039/d3sm00497j, 2023. a
Taylor, P., Hegyi, B., Boeke, R., and Boisvert, L.: On the Increasing Importance of AirSea Exchanges in a Thawing Arctic: A Review, Atmosphere, 9, 41, https://doi.org/10.3390/atmos9020041, 2018. a
The NEMO team: NEMO release version 4.2.2, Zenodo [code], https://doi.org/10.5281/zenodo.11551599, 2024. a
Tintó Prims, O., Castrillo, M., Acosta, M. C., MulaValls, O., Sanchez Lorente, A., Serradell, K., Cortés, A., and DoblasReyes, F. J.: Finding, analysing and solving MPI communication bottlenecks in Earth System models, J. Comput. Sci.Neth., 36, 100864, https://doi.org/10.1016/j.jocs.2018.04.015, 2019. a
Tremblay, L.B. and Mysak, L. A.: Modeling sea ice as a granular material, including the dilatancy effect, J. Phys. Oceanogr., 27, 2342–2360, https://doi.org/10.1175/15200485(1997)027<2342:msiaag>2.0.co;2, 1997. a
Vancoppenolle, M., Rousset, C., Blockley, E., Aksenov, Y., Feltham, D., Fichefet, T., Garric, G., Guémas, V., Iovino, D., Keeley, S., Madec, G., Massonnet, F., Ridley, J., Schroeder, D., and Tietsche, S.: SI3, the NEMO Sea Ice Engine, Zenodo [code], https://doi.org/10.5281/ZENODO.7534900, 2023. a, b
Vihma, T.: Effects of Arctic Sea Ice Decline on Weather and Climate: A Review, Surv. Geophys., 35, 1175–1214, https://doi.org/10.1007/s1071201492840, 2014. a
Weiss, J. and Dansereau, V.: Linking scales in sea ice mechanics, Philos. T. Roy. Soc. A, 375, 20150352, https://doi.org/10.1098/rsta.2015.0352, 2017. a
Weiss, J., Marsan, D., and Rampal, P.: Space and time scaling laws induced by the multiscale fracturing of the Arctic sea ice cover, in: IUTAM Symposium on scaling in solid mechanics, 25–29 June 2007, Cardiff, UK, 10, https://doi.org/10.1007/9781402090332, 2009. a
 Abstract
 Introduction
 Model and implementation
 Model evaluation
 Discussion
 Conclusions
 Appendix A: Nomenclature
 Appendix B: Algorithm and discretization
 Appendix C: Additional figures
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Model and implementation
 Model evaluation
 Discussion
 Conclusions
 Appendix A: Nomenclature
 Appendix B: Algorithm and discretization
 Appendix C: Additional figures
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References