A New Brittle Rheology and Numerical Framework for Large‐Scale Sea‐Ice Models

Abstract

We present a new brittle rheology and an accompanying numerical framework for large‐scale sea‐ice modeling. This rheology is based on a Bingham‐Maxwell constitutive model and the Maxwell‐Elasto‐Brittle (MEB) rheology, the latter of which has previously been used to model sea ice. The key strength of the MEB rheology is its ability to represent the scaling properties of simulated sea‐ice deformation in space and time. The new rheology we propose here, which we refer to as the brittle Bingham‐Maxwell rheology (BBM), represents a further evolution of the MEB rheology. It is developed to address two main shortcomings of the MEB rheology and numerical implementation we were unable to address previously: excessive thickening of the ice in model runs longer than about one winter and a relatively high computational cost. In the BBM rheology and numerical framework these shortcomings are addressed by demanding that the ice deforms under convergence in a purely elastic manner when internal stresses lie below a given compressive threshold. Numerical performance is improved by introducing an explicit scheme to solve the ice momentum equation. In this paper, we introduce the new rheology and numerical framework. Using an implementation of BBM in version two of the neXtSIM sea‐ice model (neXtSIMv2), we show that it gives reasonable long term evolution of the Arctic sea‐ice cover and very good deformation fields and statistics compared to satellite observations.

Introduction

The drift and deformation of sea ice is a key aspect of the over-all state of the ice cover. Large-scale drift redistributes ice, affecting where it forms, melts, and is collected, while small scale deformation opens up leads and builds ridges, which influence virtually all interactions between the atmosphere, ocean, and ice in ice-covered areas. The pan-Arctic drift and thickness distribution are relatively well observed (e.g., Colony & Thorndike, 1984; Kwok et al., 2013; Kwok & Rothrock, 2009; Ricker et al., 2017; Rothrock et al., 2008), while lead and ridge formation can be both directly observed at high resolution and linked to the Linear Kinematic Features (LKFs) observed from satellite (Kwok et al., 1998).

The drift and deformation of ice in a sea-ice model is determined by the solution of the momentum equation. This equation has several terms, with one of the most important ones being the internal stress term (e.g., Steele et al., 1997). The relationship between the internal stress and resulting deformation is referred to as a rheology and virtually all continuum, geophysical-scale sea-ice models used currently employ the viscous-plastic rheology (VP; Hibler, 1979) or the elastic-viscous-plastic rheology (EVP; Hunke & Dukowicz, 1997), which only addresses numerical issues with the VP. The VP rheology treats the ice as a continuum and assumes it deforms in a viscous manner with a high viscosity until the internal stress reaches a plastic threshold, determined by a yield curve which usually has an elliptic shape. Several important improvements have been made to the original VP rheology (such as Hunke & Dukowicz, 1997; Lemieux et al., 2010; Bouillon et al., 2013; Kimmritz et al., 2016), but the physical principles remain the same.

The VP rheology has enjoyed tremendous success and is used for time scales from days to centuries and spatial scales from tens of kilometres to basin scales. It is, however, not without its faults, both when it comes to the underlying assumptions (see in particular Coon et al., 2007) and the results produced by models that use it. There is generally a very large spread in key prognostic variables such as thickness, concentration, and drift in model inter-comparison studies—well beyond observed variability (Chevallier et al., 2016; Tandon et al., 2018). The sharp gradients in velocities, which are known as LKFs and are related to ridge and lead formation, are also poorly reproduced in any VP-based model running at a coarser resolution than about 2 km, a resolution that is an order of magnitude higher than the observational data (Hutter & Losch, 2020; Spreen et al., 2017). While it is not clear whether these shortcomings are due to the VP physics, numerics, or other factors (e.g., Bouchat et al., 2022; Hutter et al., 2022), modifying the model physics is a plausible avenue of investigation. Several authors have, therefore, suggested alternate approaches to the VP rheology, such as Tremblay and Mysak (1997); Wilchinsky and Feltham (2004); Schreyer et al. (2006); Girard et al. (2011); Dansereau et al. (2016).

The rheology presented here is the latest realisation of a branch of rheologies that traces its origin back to investigations of satellite observations obtained with the Radarsat Geophysical Processing System (RGPS, Kwok et al., 1998) and buoys trajectories from the International Arctic Buoy Program (IABP). Both data sets have proven to be a particularly rich source of information on sea-ice dynamics. For the sake of the current discussion, the most important result of the investigations of the RGPS data set is the discovery of the existence of a spatial scale invariance in the way sea ice deforms and of its associated fractal properties (e.g., Hutchings et al., 2011; Marsan et al., 2004; Oikkonen et al., 2017; Rampal et al., 2008; Weiss & Marsan, 2004). These observations indicate a possible way forward for the development of sea-ice rheological models: to be consistent with the observations the models must represent the propagation of fracturing and the associated spatial and temporal correlations in the sea-ice deformation field, and they must include a sub-grid-scale parameterization of the fracturing.

Sea-ice models using the VP rheology have been shown to capture the grid-scale propagation of fracturing for scales that are about an order of magnitude lager than the model resolution (Bouchat et al., 2022; Girard et al., 2011; Hutter & Losch, 2020; Spreen et al., 2017). This is witnessed by the fact that the models exhibit spatial scaling at these larger scales, albeit sometimes with the wrong power law exponent. The fact that they don't exhibit scaling at, or near the model resolution strongly indicates that they lack a good sub-grid-scale parameterization of fracturing.

It is important to consider the sub-grid-scale behavior because the triggering of fracture formation will always occur at scales much smaller than the model scale (possibly as small as the meter scale). This unresolved process must, therefore, be properly parameterized in order for the model to be physically consistent at the grid scale and, as much as possible, not resolution dependant. Given the observed scale invariance of sea-ice deformation and related quantities (e.g., Marsan et al., 2004; Rampal et al., 2009, 2008; Ólason et al., 2021) we can also assume that correctly capturing the small scale behavior will affect what happens at a larger scale.

Following these ideas and the work of Marsan et al. (2004), Weiss and Marsan (2004), Schulson (2004), Schulson and Hibler (2004), and Weiss et al. (2007), Girard et al. (2011) suggested the elasto-brittle (EB) rheology. This is a damage propagation model where the fracture density or damage at the sub-grid scale is parameterized using a single scalar variable which value is altered whenever the local stress exceeds the Mohr-Coulomb failure criterion. Girard et al. (2011) showed that the EB model could reproduce not only the observed spatial scaling, but also the localisation and other qualitative properties of the deformation field. Following this, Dansereau et al. (2016) then proposed a further development of the EB rheology in the form of the Maxwell-elasto-brittle (MEB) rheology. The MEB is a viscous-elastic rheology which allows the model to simulate also the large—and permanent—deformations occurring once the ice is fractured and fragmented. In parallel, Bouillon and Rampal (2015), Rampal et al. (2016), and Rampal et al. (2019) implemented and used the EB and MEB rheologies in the neXtSIM large-scale sea-ice model to evaluate these rheologies against observations over spatial and temporal scales spanning several orders of magnitudes.

Despite the very encouraging results of Dansereau et al. (2016), Dansereau et al. (2017), Rampal et al. (2019), and Ólason et al. (2021), the MEB rheology as proposed by Dansereau et al. (2016) and implemented in Rampal et al. (2019), leads to excessive convergence of highly damaged ice, causing it to pile up and become unrealistically thick, a problem not experienced by models using the VP rheology. Furthermore, in order to achieve acceptable numerical performance for longer simulations, Rampal et al. (2019) used a much longer time step than Dansereau et al. (2016) and did not use a fixed-point iteration scheme like Dansereau et al. (2016). This causes the model not to converge to the correct solution, impacts the damage propagation, and ultimately leads to a substantial dependence of model behavior on the time step. In this paper we present a new physical and numerical framework designed to address those issues, while retaining the main characteristics and results already obtained using MEB.

In the following we will first present the revised rheology and proposed numerical framework, discussing both the use of the Bingham-Maxwell constitutive model in a damage-propagation framework and the use of an explicit solver to improve the code's efficiency. We then evaluate this rheology and framework as implemented in the neXtSIM sea-ice model. We consider this a sufficiently substantial improvement of the model for it to now warrant the neXtSIMv2 moniker, which we will use hereafter to refer to neXtSIM with the Bingham-Maxwell rheology (BBM) rheology. In Section 3 we first evaluate model results against the RGPS observations, demonstrating the model's abilities in reproducing certain observed large-scale properties of sea-ice deformation. Thereafter, in Section 4, we demonstrate that this new framework gives very reasonable results in terms of large-scale drift and thickness distribution in a decade-long simulation of the Arctic ice cover. In Section 5 we then discuss the model's sensitivity to key parameters.

Model Description Motivation

Before describing in detail the modeling framework we discuss the rationale behind the change suggested to the MEB rheology and the new numerical implementation. These are the addition of a threshold for permanent deformation in compression and the use of an explicit solver, respectively.

Our motivation behind amending the MEB rheology is that neither the EB nor the MEB rheologies provide sufficient resistance to ice compression. This is because once damaged, the ice compresses readily allowing prevailing winds and currents to pile up very thick ice without any substantial resistance. For simulations lasting more than about a year this results in the formation of unrealistic, thick ice patches (thicker than 5 m, see Figure 1) of which the number and thickness increase over time. Our approach in addressing this is to replace the Maxwell constitutive model used in MEB with a Bingham-Maxwell constitutive model (e.g., Bingham, 1922; Cheddadi et al., 2008; Irgens, 2008; Saramito, 2021). Using this constitutive model in the context of sea ice was originally suggested by Dansereau (2016), although they suggested a different stress criterion. Schematically speaking, the Bingham-Maxwell constitutive model consists of a dashpot and a friction element in parallel, connected to a spring in series (Figure 2), with the friction element being the key distinguishing feature between MEB and BBM. The dashpot and spring still follow the same visco-elastic rheology coupled to a progressive damage mechanism as in Dansereau et al. (2016), while the condition we use for the friction element is that for −P_max < σ_N < 0 we have elastic behavior without permanent deformations, while otherwise we have both elastic and stress-dissipative behavior. Here σ_N is the mean normal stress in the ice and P_max is a compressive strength threshold. This setup is chosen to simulate ridging in high compression and a resistance to ridging when the compressive stress is below a threshold. Different formulations of the threshold are possible (including the one suggested by Dansereau, 2016, to represent friction between ice floes), but the one above is designed to treat compression and give the best results in both preventing excessive convergence and producing reasonable deformation results as discussed in the following sections.

Figure 1. Snapshot of simulated sea ice thickness distribution on 1 January 1999, after 4 years of simulation using the MEB rheology in neXtSIM.

Figure 2. Panel (a) A schematic of the Bingham-Maxwell constitutive model showing a dashpot and a friction element connected in parallel, with both connected to a spring in series. Panel (b) The yield criterion in the stress invariant plane {σN, τ}, as well as the elastic limit Pmax, and the ridging (I), elastic (II), and diverging (III) regimes.

The justification for using an explicit solver lies in the necessity to capture the propagation of damage while optimizing simulation times. Dansereau et al. (2016) introduced the concept of a characteristic time scale for damage evolution, t_d, as the time of propagation of (shear) elastic waves and used a semi-implicit scheme with a fixed-point (Picard) iteration scheme with a time step Δt ≥ t_d. Such a scheme is computationally demanding and Rampal et al. (2019) eventually used a semi-implicit solver, without a fixed-point iteration scheme, and Δt ≫ t_d, to reduce computational cost. As a result, their model results are dependent on the time-step length and the solution is most likely not fully converged. In opting for an explicit solver with a time-splitting scheme we update only rapidly changing variables (velocity, stress, and damage) at a short time step, while doing advection and thermodynamics at a longer time step. This is based on the fact that fracture formation happens at a speed similar to that of sound in the ice and is thus much faster than the sea ice drift speed. The use of an explicit solver is also inspired by the work of Hunke and Dukowicz (1997), who showed that in the case of the EVP model one can use a time step for the explicit solver determined by the elastic time scale and not the much shorter viscous time scale. This result also holds here (see Appendix A).

Using an explicit solver requires Δt < t_d to explicitly resolve the damage propagation. This time-step requirement is, however, not particularly imposing, as t_d ∝ Δx (see Appendix A) and there is considerable experience within the sea-ice modeling community in solving the sea-ice momentum equation explicitly in a computationally efficient manner. This was in fact the main goal of Hunke and Dukowicz (1997) in choosing an explicit solver for the EVP rheology. Moreover, typical values of t_d are similar to, or even larger, than values typically used for the elastic time scale of the EVP rheology. It is, therefore, reasonable to assume that the same sub-time stepping approach can be used here as in the EVP rheology. It is important to note that elasticity in the EVP rheology is not intended to be physical, but is introduced for numerical expediency and elastic waves in EVP should, therefore, be damped (e.g., Bouillon et al., 2013). Elasticity in BBM is, however, physical so there is no need to damp any resulting elastic waves.

The Brittle Bingham-Maxwell Constitutive Model

The EB and MEB rheologies are centered around the idea of damaging and damage propagation, and the BBM also relies on this concept, using the same damaging mechanism as MEB. The key difference between these rheologies lies in the constitutive model, with the EB using a damaging spring, MEB using a damaging Maxwell model, and the BBM being a damaging Bingham-Maxwell model. The Maxwell model consists of a dashpot and a spring in parallel, while the Bingham-Maxwell model consists of a dashpot and a friction element in parallel, connected in series with a spring (Figure 2). The inclusion of a friction element is thus the key difference between MEB and BBM. Here we will derive the constitutive model resulting from the use of a Bingham-Maxwell constitutive model with damage, link this to the damage mechanism, and then present the appropriate temporal discretization of the system.

Constitutive Model

The constitutive model used here is the Bingham-Maxwell model together with a dependence of elasticity and viscosity on damage. The Bingham-Maxwell model is a set up of a dashpot and friction element in parallel, connected in series with a spring (Figure 2). The condition we use for the friction element is defined in terms of the normal stress [Image Omitted. See PDF]as we aim to prevent excessive thickening. In divergent conditions (σ_N > 0), the stress in the friction element is 0 and only the dashpot is active. In this case the total stress is the same as the elastic stress and the viscous stress (σ = σ_E = σ_v) and the total displacement is the sum of the elastic and viscous displacements [Image Omitted. See PDF]In the range −P_max < σ_N < 0, the friction element is able to prevent any permanent deformation (ɛ_v = 0 and ɛ = ɛ_E) and we have a pure elastic behavior, with [Image Omitted. See PDF]For σ_N < −P_max, the friction element is no longer able to prevent permanent convergent deformation. We note that P_max is the key quantity introduced in the BBM rheology, compared to the MEB.

In a one-dimensional Bingham-Maxwell constitutive model (as in Figure 2, panel b) the friction element stress is constant (at P_max) and the viscous stress is related to the total stress by [Image Omitted. See PDF]which may be rewritten as [Image Omitted. See PDF]In the two dimensional case we use the normal stress σ_N as threshold to get [Image Omitted. See PDF]This ensures that the simulated ice retains some resistance to compression, even in a highly damaged state. Recalling Figure 2, we generalize the relationship between σ and σ_v as [Image Omitted. See PDF] [Image Omitted. See PDF]The threshold P_max thus separates the elastic and visco-elastic, or reversible and permanent deformation phases of the Bingham-Maxwell constitutive model. We assume that there is a relationship between the threshold P_max and ice thickness, which is related to the process of ridging, and so we have used the form [Image Omitted. See PDF]where h₀ = 1 m is a constant reference thickness, P a constant to parameterize P_max, following the results of Hopkins (1998), and C is a constant similar to the compaction parameter introduced by Hibler (1979). Different formulations for P_max may be considered, as briefly discussed in Section 5.

Brittle behavior is ensured by using a slightly modified version of the damaging mechanism of Dansereau et al. (2016). We write the elasticity E and viscosity η as a function of damage d and ice concentration A as [Image Omitted. See PDF] [Image Omitted. See PDF]where E₀ and η₀ are the undamaged elasticity and viscosity, and α > 1 is a constant. Undamaged ice has d = 0, while highly damaged ice has d → 1 and d = 1 is never reached. We use a different dependence of η on A compared to Dansereau et al. (2016), using e^{−Cα(1−A)}, instead of e^−C(1−A). This gives more realistic behavior at low and medium ice concentration, as discussed further in Section 5.

Following Dansereau et al. (2016), we can now apply the elastic stiffness tensor K to the time derivative of Equation 2 and multiply with the elasticity to get [Image Omitted. See PDF]We assume plane stress conditions, so the stiffness tensor operation $$\mathbf{K}:\dot{\varepsilon }$$ is [Image Omitted. See PDF]where ν is Poisson's ratio. As the elastic stress is, by definition of Equation 3 [Image Omitted. See PDF]its time derivative is [Image Omitted. See PDF]Calculating $$\dot{E}$$ from Equation 9 we get [Image Omitted. See PDF]noting that changes in concentration, A, are much slower and can be ignored (see Appendix B for details).

The viscous stress then relates to the viscous displacement as [Image Omitted. See PDF]and to the total stress by [Image Omitted. See PDF]The elastic stress is related to the total stress as [Image Omitted. See PDF]since the stress in each serially connected element must be equal to the total stress. By using Equations 7, 15–18 we can now rewrite Equation 11 as [Image Omitted. See PDF]or [Image Omitted. See PDF]where λ = η/E = λ₀(1 − d)^α−1 is the viscous relaxation time, with λ₀ the undamaged viscous relaxation time.

For the time rate of change of damage, $$\dot{d}$$ we have $$\dot{d} > 0$$ only when damaging occurs, otherwise $$\dot{d}=0$$. We will, therefore, link the $$-\sigma \dot{d}/(1-d)$$ term of Equation 20 to the damaging process below, noting that this term of the equation is zero when the stress is inside the failure envelope. Note also, that for $$\dot{d}=0$$ and $$\tilde{P}=0$$ the MEB constitutive law is recovered (Equation 4 of Dansereau et al., 2016).

Damaging and Healing

Damaging occurs in the BBM rheology whenever the simulated stress in a grid cell or element is outside the failure envelope, or yield curve. The failure envelope of the BBM rheology is the Mohr-Coulomb criterion: [Image Omitted. See PDF]where τ and σ_N are the stress invariants (shear and mean normal stress, respectively), μ is the internal friction coefficient and c is the cohesion (see Figure 2). Following Bouillon and Rampal (2015), we let the cohesion scale with model resolution, as [Image Omitted. See PDF]where c is the model cohesion, Δx is the distance between model node points, and c_ref is the cohesion at the reference length scale, l_ref. We here use the lab scale, l_ref = 10 cm as the reference length scale, where we know the cohesion to be of the order of 1 MPa (e.g., Schulson et al., 2006). In addition to the Mohr-Coulomb criterion we cap the yield curve at high compressive normal stress, as discussed below.

Stress states outside the failure envelope are not physical and so whenever the modeled stress states fall outside the criterion, the damage, d, is modified in order to bring the stresses back inside the yield curve. We note that Equation 20 can be written as [Image Omitted. See PDF]with the last term being [Image Omitted. See PDF]We now consider the case of damaging changing the stress from a stress state outside the yield curve, σ′, to a stress state on the failure envelope, σ, over a time t_d. We then have [Image Omitted. See PDF]and [Image Omitted. See PDF]Assuming that the damaging process is uniform over time t_d, we can write this as [Image Omitted. See PDF]Combining Equations 24 and 27 we get [Image Omitted. See PDF]We can assume that for stresses inside the yield curve d_crit = 0 at all times. Following Dansereau et al. (2016), we set the characteristic time scale of the propagation of damage to be [Image Omitted. See PDF]with c_E being the propagation speed of an elastic shear wave, ν being Poisson's ratio, ρ the ice density, and E the elasticity as before. We note that Equation 27 gives an equation for the change in stress due to damaging which allows us to simplify the time stepping described below.

The variable d_crit is the distance between the simulated stress and the yield curve. Here we use the formulation of Plante et al. (2020), but limiting on the compressive stress following (Bouillon & Rampal, 2015). This upper limit is there to counteract instabilities that set in at very high σ_N (as pointed out by Plante et al., 2020). This limit is a numerical construct and is chosen high enough so that it does not influence the results. We scale the limit using the same scaling relationship as for the cohesion, as the onset of instability at high compression is related to the value of cohesion. The resulting equation for the limit is [Image Omitted. See PDF]where N_ref is the limit at the reference length scale l_ref. The resulting equation for d_crit then reads [Image Omitted. See PDF]Using this formulation, stress states outside the yield curve have 0 < d_crit < 1. Since the elasticity and viscosity of the rheology depends on the damage, the damaging process as described above ensures that the stresses are always relaxed to within the yield curve over the time scale t_d.

Damaged ice must heal with time and this is done via a simple restoring term as originally introduced by Bouillon and Rampal (2015) and used in Rampal et al. (2016). [Image Omitted. See PDF]Here t_h is the characteristic time scale of healing, which we chose to depend on the temperature difference between the base of the ice and of the snow-ice interface, that is, t_h = k_th/ΔT, where k_th is a constant and ΔT is the temperature difference. This formulation is somewhat ad hoc, but it prevents melting ice from healing which improves thickness and concentration distribution in summer and has very limited effect in winter. The time scale of healing is much larger than that of damaging (t_h ≫ t_d), and so Equations 28 and 32 can be solved separately.

Temporal Discretization

The temporal discretization of Equation 20, using an implicit scheme, is relatively straightforward and very similar to that of Dansereau et al. (2016). Assuming no damage occurs, $$\dot{d}=0$$ and we write $$\dot{\sigma }$$ in terms of the previous time step and the current estimate, σⁿ and σ′ respectively, giving [Image Omitted. See PDF]where all variables are from the previous time step (n), and Δt is the time-step length. Rearranging gives [Image Omitted. See PDF]

If the stress σ′ is inside the failure envelope we set σⁿ⁺¹ = σ′ and continue. If the stress is outside the envelope, however, damaging occurs. In this case, damage is updated using the damage evolution in Equation 28, which should be discretized explicitly as [Image Omitted. See PDF]This can be rearranged as [Image Omitted. See PDF]The super-critical stress resulting from Equation 34 is then relaxed assuming a discretization of Equation 27 of the form [Image Omitted. See PDF]which can be rewritten as [Image Omitted. See PDF]This relaxation may also be replaced by a recalculation of the stress using the full Equation 20 and dⁿ⁺¹, but using Equation 38 is substantially more cost-efficient and the results are virtually identical.

An Explicit Solver for the Momentum Equation

The use of an explicit solver for the sea-ice momentum equation is well known within the sea-ice modeling community, and the current implementation follows very closely that of Hunke and Dukowicz (1997) and Danilov et al. (2015). There have been various improvements made to the EVP rheology of Hunke and Dukowicz (1997) in the last years (Bouillon et al., 2013; Kimmritz et al., 2016; Lemieux et al., 2012), attempting to find the best way of using a sub-time stepping scheme to converge the EVP solution to the implicit VP solution. In our case it is more appropriate to think not of sub-time stepping, but rather time-splitting, where the dynamic processes changing velocity and internal stress are resolved at a much shorter time step than advection and thermodynamic processes. Such time-splitting is well known in ocean models (e.g., Hallberg, 1997; Killworth et al., 1991) and the original EVP of Hunke and Dukowicz (1997) can also be considered as a time-splitting approach. We base our solver very closely on that of Hunke and Dukowicz (1997), it being a good fit for our purpose, and a widely adopted and well-understood method.

The momentum equation of sea ice is (e.g., Bouillon & Rampal, 2015; Connolley et al., 2004; Danilov et al., 2015) [Image Omitted. See PDF]where m = Aρh is the ice mass per unit area, $$\vec{u}$$ is the ice velocity, σ is the internal stress tensor, h is the ice slab thickness (not volume per unit area), ρ the ice density, $${\vec{\tau }}_{a}$$ and $${\vec{\tau }}_{w}$$ are the atmosphere and ocean stress terms, respectively, $${\vec{\tau }}_{b}=-{C}_{b}\vec{u}$$ is the basal stress term introduced in Lemieux et al. (2015), $$mf\vec{k}\times \vec{u}$$ is the Coriolis term, with vertical unit vector $$\vec{k}$$, and $$mg\vec{\nabla }\eta $$ is the ocean-tilt term. We write explicitly the integrated internal stress as σh following Sulsky et al. (2007) and Bouillon and Rampal (2015). On staggered grids, the tracers m, A, and h are generally collocated and not collocated with the velocities, so we prefer to divide Equation 39 with A to get [Image Omitted. See PDF]ignoring a factor of 1/A in front of the internal and basal stress terms. Those terms disappear very quickly with decreasing concentration, so the error incurred is very small (of the order of 0.1%). The correct dependence of these terms on A is also very uncertain.

The atmosphere and ocean stress terms are assumed to be quadratic, having the forms [Image Omitted. See PDF]and [Image Omitted. See PDF]respectively, where ρ_a and ρ_w are the atmosphere and ocean densities, C_a and C_w atmosphere and ocean drag coefficients, θ_a and θ_w the atmosphere and ocean turning angles, and $${\vec{u}}_{w}$$ is the ocean velocity.

The momentum equation, together with the drag terms, is then discretized in time (using Cartesian coordinates with i, j = 1, 2 implying x and y direction) as (Hunke & Dukowicz, 1997) [Image Omitted. See PDF]where ɛ_ijk is here the Levi-Civita symbol and $${c}^{\prime }={\rho }_{w}{C}_{w}\vert {\vec{u}}_{w}-{\vec{u}}^{k}\vert $$. This then gives a set of equations that can be solved for the velocity components to give [Image Omitted. See PDF] [Image Omitted. See PDF]with [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]Given a form for σ^k+1 and a spatial discretization, Equations 44 and 45 are easily solved to give the velocity components at each grid point or mesh node.

In this approach σ^k+1 depends on σ^k and $${\left({u}_{1},{u}_{2}\right)}^{k}$$, through $${\dot{\varepsilon }}^{k}$$ in Equation 34. A more correct temporal discretization of Equation 20 would use $${\dot{\varepsilon }}^{k+1}$$, but this is only available when solving the momentum equation implicitly. Using $${\dot{\varepsilon }}^{k}$$ and not $${\dot{\varepsilon }}^{k+1}$$ will result in an error in the estimate of σ^k+1, which in turn may lead to excessive damaging as well. We have not investigated the extent to which this affects the results, but a way to do so is to iterate over the equations for σ^k+1 and (44) and (45) until the velocity used to calculate σ^k+1 have converged to $${\left({u}_{1},{u}_{2}\right)}^{k+1}$$.

The spatial discretization of Equations 44 and 45 using finite differences was discussed by Hunke and Dukowicz (1997) for an Arakawa B-grid and by Bouillon et al. (2009) for both the Arakawa B and C-grids. As we have chosen to implement the new rheology in the finite element model neXtSIMv2, we have followed Danilov et al. (2015) for a discretization using the finite elements method, but there are no apparent impediments for a finite difference implementation of the new rheology.

In their implementation of the Finite Element sea-ice model, FESIM (version 2), Danilov et al. (2015) use nodal quadratures in all terms that do not involve spatial derivatives, in order to save computational time by not computing (unnecessary) mass matrices. The authors derive a weak formulation of the momentum Equation 40 by multiplying it with test functions, integrating over the domain, and integrating the internal stress term by parts to get a weak formulation. As the resulting lumped mass matrix $$\left({M}_{lm}^{L}\right)$$ is diagonal, its diagonal entries are simply one third of the sums of areas of triangles containing the vertex considered, A_c/3. Equations 44 and 45 can then be used directly, but with [Image Omitted. See PDF] [Image Omitted. See PDF]and [Image Omitted. See PDF] [Image Omitted. See PDF]where ∑_c(l) denotes the sum over all the elements adjacent to node l and ∑_m(c) denotes the sum over all the nodes belonging to element c. Note that in neXtSIMv2 the momentum equation is solved on the polar-stereographic plane and we do not include the metric factors as present in Danilov et al. (2015).

Implementation

The implementation of BBM into neXtSIMv2 that is used hereafter uses a time-splitting method wherein all equations except those related to the velocity, stress, and damage updates are solved using a long, main time step, Δt_m. This includes damage healing, according to Equation 32, thermodynamics, and advection. The velocity, stress, and damage fields (except for healing) are then updated at a higher frequency. The higher frequency time stepping simply consists of splitting the long time step into N_sub short dynamical time steps, Δt. The long time step is limited by the stability of the advection scheme, while the dynamical time step is limited by the stability of the BBM rheology. In neXtSIMv2, a single dynamical time step consists of the following: (Algorithm 1)

1 AlgorithmA Single Dynamical Time Step in the Implementation of BBM Into neXtSIMv2

1. Calculate σ_N and P_max according to Equations 1 and 8, respectively 2. Calculate σ′ according to Equation 34 3. Calculate d_crit according to Equation 31 4. if d_crit < 1 then 5.    Update damage according to Equation 36 6.    Update σⁿ⁺¹ according to Equation 38 7. else 8.    Set σⁿ⁺¹ = σ′ 9. end if 10. Calculate u₁ and u₂ using Equations 44 and 45 11. Update u₁ and u₂ on ghost-nodes of the parallelization sub-domains

In addition to the dynamical solver described here, thermodynamic growth is calculated using the approach of Winton (2000) and advection is done using the Lagrangian scheme of Samaké et al. (2017). We also use the two-class approach of Hibler (1979) for calculating ice growth in open water.

Evaluation of Simulated Deformation

Here we present a simplified evaluation of the simulated deformation. This evaluation was performed by qualitative visual analysis of deformation maps (see Figures 3 and 4), probability density functions, quantitative metrics including bias and root mean square error of deformation time series, and spatial scaling analysis. The goal of applying these metrics on the two model runs is to illustrate the sensitivity of the metrics to obviously different spatial patterns of deformation, rather than a comprehensive evaluation of the different rheologies.

Figure 3. Maps of sea ice divergence (day−1) for 2 February 2007 as observed by radarsat geophysical processing system and simulated by neXtSIMv2 using the brittle Bingham-Maxwell and mEVP rheologies.

Figure 4. Maps of sea ice shear (day−1) for 2 February 2007 as observed by radarsat geophysical processing system and simulated by neXtSIMv2 using the brittle Bingham-Maxwell and mEVP rheologies.

As explained in subsections below the metrics were computed for sea ice deformation from three sources of ice drift:

• SAR-based observations of ice drift from the RADARSAT Geophysical Processor System (RGPS, Kwok et al., 1998);

• neXtSIMv2 with the new BBM rheology (BBM);

• neXtSIMv2 with the mEVP rheology (Bouillon et al., 2009);

The main goal here is to compare BBM against observations. We include the mEVP simulations as a reference for the commonly used (E)VP models and we choose not to compare to results obtained with MEB, since we have already established that it is not suitable for longer simulations.

The model setup is similar to that in Rampal et al. (2019), except that here we use the BBM where they used MEB. In the two runs (BBM, mEVP) neXtSIMv2 is initialized on 15 November 2006 and runs until 30 April 2007. Atmospheric forcing is derived from the ERA5 reanalysis (Hersbach et al., 2020) and oceanic forcing from the TOPAZ4 reanalysis (Sakov et al., 2012). Initial sea ice thickness and concentration are set from a combination of data from OSISAF (Tonboe et al., 2016), TOPAZ4, and ICESAT (available at: https://icdc.cen.uni-hamburg.de/seaicethickness-satobs-arc.html, last access: August 2020), as described in Rampal et al. (2019). Initial sea ice damage is set to zero. In all three runs the explicit solver is used and the time step and spatial resolution are the same. The difference is in the rheological part of the model: BBM uses equations from Section 2.2 as they are, in mEVP we follow the implementation of Danilov et al. (2015) with minor changes discussed in Appendix C. We use model time steps of Δt_m = 900 s and Δt = 7.5 s, which is equivalent to 120 sub-iterations, for both BBM and mEVP. For the mEVP we set the α_mEVP and β_mEVP parameters to 500 following Koldunov et al. (2019). We also tested running the mEVP with 500 and 1,000 sub-iterations, but the differences in results are minor (see Appendix D). Table 1 lists the main model parameters and the values used here.

Table 1 Key Model Parameters and the Values Used in the Experiments Presented Here

Parameter	Symbol	Value
Ice–atmosphere drag coefficient	Ca	2.0 × 10−3
Ice–ocean drag coefficient	Cw	5.5 × 10−3
Undamaged elasticity	E0	5.96 × 108 Pa
Undamaged viscous relaxation time	λ0	1 × 107 s
Damage parameter	α	5
Scaling parameter for the riding threshold	P	10 kPa
Cohesion at the reference scale	cref	2 MPa
Poisson ratio	ν	1/3
Ice density	ρ	917 kg/m3
Maximum compressive stress at the reference scale	Nref	10 GPa
Temperature dependent healing time scale	kth	15 days/20 K
Main model time step	Δtm	900 s
Dynamical time step	Δt	7.5 s
Mean resolution	Δx	10 km
mEVP convergence parameters	αmEVP, βmEVP	500
mEVP ellipse aspect ratio	e	2
mEVP ice strength	P*	27.5 kN/m2
mEVP ice tensile strength	T*	0 kN/m2

Details on Computation of Sea Ice Deformation Rates

Sea-ice drift is computed from the RGPS data the same way as in Stern and Lindsay (2009), with “snapshots” of the sea-ice drift created from the Lagrangian displacement data. For a given target time the snapshot contains all observations of drift that start before this time, end after it and are separated by 3 days. Sea-ice drift from the model is computed similar to Rampal et al. (2019), with drifters in the model seeded at the location of the RGPS snapshot points, and these drifters then advected together with the model elements for the same duration as in the RGPS snapshot. Unlike in Rampal et al. (2019), the simulated trajectories are re-initialized every 3 days to exactly match the RGPS snapshots. The sea ice deformation components divergence (ɛ_div) and shear (ɛ_shear) formulation are exactly the same as in Rampal et al. (2019): [Image Omitted. See PDF] [Image Omitted. See PDF]where u_x, u_y, v_x and v_y are components of the ice drift velocity gradient.

Maps of divergence and shear rate computed from an example snapshot of RGPS-data based sea-ice drift for 2 February 2007 are compared against modeled results in Figures 3 and 4. Similar to maps in Rampal et al. (2019) and Marsan et al. (2004) the RGPS maps clearly show presence of narrow and long fractures in sea ice in the central Arctic, while the deformation field closer to the coast is more homogeneous. Visually the BBM maps appear quite realistic—length, width and orientation of fractures, as well as magnitude of deformation rates is similar to the RGPS observations. The mEVP maps, on the other hand, show very smooth fields of deformation with few obvious ice cracks.

Sea Ice Deformation Probability Distribution

Probability density functions (PDFs) were computed from all snapshots of sea ice deformation components for RGPS, BBM and mEVP and plotted in Figure 5. Comparison of PDFs shows that for both divergence and shear BBM fits very well with observations, yet slightly underestimating the highest shear values. High values of convergence (above 0.1 day⁻¹) (defined as negative values of divergence with opposite sign) are underestimated. mEVP, on the other hand overestimates very small deformations and significantly underestimates the main portion of the spectrum.

Figure 5. Probability density functions of three sea ice deformation components computed from all snapshots in 2007. Colors denote radarsat geophysical processing system observations (blue) and nextSIMv2 runs: brittle Bingham-Maxwell (orange) and mEVP (green) rheologies.

Sea-Ice Deformation Time Series

We have seen that both the spatial field and the PDFs are characterized by a small number of high deformation values. This is exemplified by the LKFs (Figures 3 and 4) and the long tail of the PDFs (Figure 5). To better analyze this, a metric sensitive to these high values should be used. The 90th percentile (denoted as P90) was selected as such a metric. P90 is the value of deformation below which 90% of deformation values in the frequency distribution fall. For evaluation of the temporal evolution of the deformation, P90 was computed from each snapshot of deformation in 2007. Values of P90 from RGPS and neXtSIMv2 were plotted and inter-compared using bias (b) and root mean square error (RMSE, e): [Image Omitted. See PDF] [Image Omitted. See PDF]where ϵ_N and ϵ_R are ice shear P90 values from neXtSIMv2 and RGPS and ⟨⟩ denotes averaging. The P90 time series (see Figure 6) show that while neither rheology can capture the highest peaks in deformation rates, the BBM results are clearly closer to RGPS, with a lower bias (b_BBM = 0.014, b_mEVP = 0.028) and RMSE (e_BBM = 0.012, e_mEVP = 0.016).

Figure 6. Time series sea ice shear P90 for 2007 as observed by radarsat geophysical processing system (blue) and simulated by neXtSIMv2 using the brittle Bingham-Maxwell (orange) and mEVP (green) rheologies.

It is noteworthy that the BBM rheology is able to instantaneously react to stronger forcing with rapidly increased deformation, and the timing of these periods of high deformation matches well with peaks in the observations. However, in the mEVP rheology deformation is lower, increases slower, and lags behind the observed rates. We expect both the P90 time series and the tail of the PDF presented in the following sub-section to be influenced by how well the atmospheric model represents extreme storms. This aspect is not investigated here.

Spatial Scaling Analysis

The spatial scaling analysis of the RGPS, BBM, and mEVP deformation distributions was performed similar to (Marsan et al., 2004). To form a distribution of the total deformation rate (ɛ_tot) at the nominal spatial scale of 10 km the triangular elements from RGPS and corresponding elements from BBM or mEVP runes were selected with the area between 40 and 60 km² (corresponding to initial RGPS triangles with sides 10 × 10 × 14 km). The shear and divergence components were computed on these triangles as described above and total deformation was computed as their geometric mean. On larger spatial scales (namely at 20, 40, 80, 160, 320, 640 and 1,000 km) the following procedure was used: the Arctic ocean was split by a grid with size equal to the analyzed spatial scale; area-weighted average of velocity gradients (u_x, u_y, v_x, v_y) from elements falling in each grid cell was computed; shear, divergence and total deformation rates were computed from the averaged velocity gradients. This procedure was repeated for 3-day fields of deformations acquired between 10 December 2006 and 10 May 2007.

The moments of distributions at each spatial scale were computed as $$\langle {\varepsilon }_{tot}^{q}\rangle $$ with order q = 1, 2 and 3. A power-law scaling function $$\langle {\varepsilon }_{tot}^{q}\rangle ={L}^{-\beta (q)}$$ was fitted for each moment using the least squares method. Moments, power-law functions and structure functions β(q) are plotted on Figure 7, where β indicates the exponent of the power-law fits and q is the moment order. The filled area indicate standard deviation from averaging moments through December 2006 - May 20.07.

Figure 7. Spatial scaling analysis of total deformation fields derived from the radarsat geophysical processing system (blue) and neXtSIMv2 runs using the brittle Bingham-Maxwell (orange) and mEVP (green) rheologies. (a) Moments of the distributions of the total deformation rate ɛtot calculated at a temporal scale of 3 days and space scales varying from 10 to 1,000 km. (b) Structure functions, where β indicates the exponent of the power-law fits and q is the moment order.

Evaluation of Simulated Thickness

One of the main motivation of the development of the BBM rheology was to be able to run long-term simulation without encountering the problem of excessive thickening that occurs with the MEB rheology as implemented by Rampal et al. (2019). In this section, we evaluate sea ice thickness in long-term simulations to ensure that BBM leads to reasonable values of the sea ice thickness, just like models using viscous-plastic based rheologies do (e.g., Zampieri et al., 2021, using mEVP).

Model Setup

We use a neXtSIMv2 setup very similar as the one used in Section 3, but with different initialization and simulation length. The model domain has been extended to encompass a larger part of the Eastern Greenland coast as well as the Barents and Kara seas (see Figure 8). Two simulations are run, one with the BBM rheology and one with the mEVP rheology. In the following, we refer to these two simulations as BBM and mEVP, respectively. The sea-ice rheology is the only difference between these two simulations. They are initialized on 1 January 1995 with ice conditions provided by PIOMAS (Schweiger et al., 2011) and are run over 20 years. Atmospheric forcings are provided by the hourly data set from the ERA5 reanalysis (Hersbach et al., 2020).

Figure 8. (a) Evolution of the 7 day running mean sea ice thickness over the domain for the mEVP and brittle Bingham-Maxwell rheology (BBM) simulations. Available data from the CS2-SMOS v2.2 product are also shown for comparison with their associated uncertainty in the shaded area. The corresponding spatial distribution for all the period covered by the CS2-SMOS v2.2 product between 2010 and 2016 is also presented for the mEVP (b) and BBM (c) simulations, as well as for the CS2-SMOS v2.2 product (d). The black solid line in (b, c, d) represents the 1.5 m sea ice thickness contour in each data set and the dashed contour line represents the borders of the model domain.

We also run 4 additional experiments using the BBM rheology to investigate the impact of the parameters P and the exponent of the thickness dependency of P_max in Equation 8. These experiments are initialized from the reference BBM simulation on 1 January 2000 and run for 5 years. The first two of them are similar to the BBM reference simulation with the exception of the value of P, set to 6 and 14 kPa. The third and fourth experiment use an exponent for the dependency of P_max on h equal to 1 and 2 respectively, instead of 3/2 in the reference simulation. The values of P in these two simulation have been adjusted to obtain the same value of P_max for h = 2 m.

Sea Ice Thickness Evaluation

For our evaluation, we compare the sea-ice thickness from the BBM and mEVP simulations to version 2.2 of the merged CS2-SMOS estimated sea thickness product (Ricker et al., 2017) (available at ftp://ftp.awi.de/sea_ice/product/cryosat2_smos/v202/nh/, last access March 2021). This product provides a 7 day averaged estimate of the pan-Arctic sea-ice thickness distribution. It is available daily during the freezing season, from mid-October to early April, starting from November 2010.

The evolution of the domain-averaged sea-ice thickness over the whole run for the two simulations is presented in Figure 8a. We used a 7 day running mean to be consistent with the CS2-SMOS estimated thickness when it is available. Here we can see that there is no spurious thickening of the sea ice in the BBM simulation, hence confirming it can be used for more than year-long simulations. The two simulations furthermore show very similar trend and inter-annual variability. The only difference is that ice is generally thicker in the BBM simulation, resulting in a positive offset of its associated curve compared to the mEVP one. The comparison with CS2-SMOS estimated thickness after 15 years of simulations show a reasonable agreement for the BBM simulation, despite a small negative bias. This negative bias is slightly larger for the mEVP simulation but can be reduced for either of these two simulations with an appropriate tuning of thermodynamical parameters.

We also check the sea ice thickness spatial distribution (Figures 8b–8d) for the overlapping period covered by the CS2-SMOS product and our simulations. In general, both simulations show distribution patterns similar to the observations, even though they underestimate the ice thickness. The extent of thick ice (represented by the 1.5 m contour in Figures 8b–8d) in the BBM simulation is however larger than in the mEVP simulation, showing a better agreement with the thick ice distribution in the CS2SMOS data set. This underestimation is particularly visible in places where ice is thicker than 2 m in the CS2-SMOS product. The underestimation of the sea ice thickness for thick ice and the overestimation of sea ice thickness for thin ice are a known problem of sea ice models (Schweiger et al., 2011). Note however that the BBM simulation seems to better reproduce the decreasing gradient of ice thickness from the northern coast of Greenland toward the North Pole than the mEVP one, in which thick ice is only found in a narrow band along the Greenland coast.

Our results show that the BBM rheology yields a reasonable sea-ice thickness magnitude and distribution when compared to observations in a way that is very similar to the results obtained with mEVP. Further studies should focus on the sea ice mass balance of a model using the BBM rheology to better understand how sea ice dynamics interact with thermodynamics.

Discussion

Given the role of spatial scaling analysis in the development of the EB and MEB models we have done a spatial scaling analysis of the BBM results as well. This shows that BBM closely follows the RGPS observations, both in terms of scaling and structure function. For P = 0 kPa we recover the MEB equations, as stated previously, and using this to run MEB within the new numerical framework shows only minor differences between the two in terms of scaling (not shown). This is consistent with previously published MEB results (e.g., Figure 3 in Rampal et al., 2019). The mEVP significantly underestimates all three moments indicating that the density distribution of deformations remain almost normal up to very small spatial scales, even if the model is run on a Lagrangian mesh. We note also that mEVP scaling results diverge significantly from the fit at the smallest scales. These results are consistent with the scaling analysis of approximately 10 km resolution (Eulerian) models performed by Bouchat et al. (2022). This shows that the source of the heterogeneity we see in the BBM runs is the model physics and not the Lagrangian advection scheme—although the advection scheme may help preserving this heterogeneity once formed.

The BBM adds to the MEB by introducing a new parameterization, which is that of the maximum pressure, P_max (see Equation 8). Here P_max is a threshold between the regimes of reversible and permanent deformations, which we interpret as the maximum pressure the ice can withstand before ridging. In Equation 8 we have chosen to use P ∝ h^3/2, leaving the constant of proportionality, P as a tunable parameter and the main new parameter of the rheology. The model results are reasonably sensitive to the value of this parameter. This is true for both the deformation patterns and the large-scale thickness distribution, both of which show a qualitatively continuous and monotonous response to changes in P for P > 0 kPa.

We explored manually the parameter space for P, and Figure 9 shows maps of shear rate for a given day and a range of values for P ∈ [0, 18] kPa, demonstrating the effect of P on the deformation patterns. Using P = 0 kPa we see that using BBM gives a qualitative improvement of the deformation patterns, compared to MEB. For P > 0 kPa there are also clear variations in the quality of the deformation patterns depending on P. For 0 < P ≲ 6 kPa the features are not as straight as expected, while for P ≳ 14 kPa they start to become too localized and intense with not enough deformation occurring between them. Modifying the cohesion (c_ref) also affects the deformation patterns; using a small value giving a large number of small, less intense features, while larger values give a smaller number of large, more intense features (not shown). A reasonable range for c_ref appears to be within 1 and 3 MPa. These comparisons are at the moment very qualitative, but we find that using the current tools we have at our disposal (such as scaling analysis and LKF detection) give either inconclusive results or require further development to be used to tune this new rheology against observed deformation.

Figure 9. Maps of sea ice shear for 2 February 2007 as simulated by neXtSIMv2 with the brittle Bingham-Maxwell rheology and P = 0, 2, 6, 10, 14, 18 kPa, in panels a, b, c, d, e, and f, respectively.

Using different values of P also affects the large-scale thickness distribution in the Arctic. Figure 10 shows how using P = 6 kPa and P = 14 kPa modifies the long term averaged thickness field, compared to P = 10 kPa. In it, we see a clear thickening by about 20 cm and thinning by about 10 cm for P = 6 kPa and P = 14 kPa, respectively. This is to be expected, as a lower P value allows the ice to ridge more readily and so the observed difference in thickness is due to an increase or decrease in ridging. We also don't expect the response to be symmetric around an optimal P value because P_max ∝ h^3/2 and not P_max ∝ h.

Figure 10. (a) January to March sea ice thickness climatology from 2000 to 2004 for the reference brittle Bingham-Maxwell rheology (BBM) simulation (P = 10 kPa and Pmax ∝ h3/2). Panels (b) and (c) show the difference for this same quantity between simulations using with P = 6 kPa (b) and P = 14 kPa (c) and the reference BBM simulation.

In addition to the sensitivity to the value of P we note that the formulation of P_max is not immediately obvious. Here we have chosen to relate the maximum pressure to ice thickness following Hopkins (1998). Other possible choices we explored were to use a constant, to use P_max ∝ h (similar to Hibler, 1979) and P_max ∝ h² (as per Rothrock, 1975). A dependence on the ice thickness is likely to be more complicated in reality, and other ice state parameters may have to be taken into account. Different formulations, such as relating P_max to the level of damage, are also possible, but were not explored here.

Using the different formulations of P_max listed above does not have a notable effect on the deformation patterns, but it does affect the large-scale thickness distribution. Figure 11 shows how using P_max ∝ h and P_max ∝ h² compares to the reference implementation with P_max ∝ h^3/2. In these experiments we chose the constant of proportionality such that P_max is the same in all three cases for 2 m thick ice. The figure shows a clear pattern of pivoting in the thickness anomalies between the different cases. For P_max ∝ h the ice that is thicker than 2 m in the reference experiment becomes even thicker, while for P_max ∝ h² it is thinner. The change in thickness is of the order of 20 cm. This behavior is expected, based on the model response to simply changing P in the reference implementation. Even though the difference between the different formulations is clear we still cannot conclusively determine which one gives the best results because uncertainties in observed ice thickness and unrelated model parameters are most likely larger than the signal we see here.

Figure 11. (a) Similar to Figure 10a. Panels (b) and (c) are also similar to Figures 10b and 10c but this time for two simulations with different dependencies of Pmax on (h) (b) Pmax ∝ h and (c) Pmax ∝ h2. Values of P in each simulation have been adjusted to obtain the same value of Pmax for h = 2 m as in the reference brittle Bingham-Maxwell rheology (BBM) simulation. The solid black line in each panel delimits the 2 m sea ice thickness contour in the BBM reference simulation.

Using the chosen set of parameters for the BBM, we see only minor differences between the thickness distribution and evolution of BBM and mEVP (Figure 8). This indicates a very strong influence of the atmospheric and oceanic forcing on the ice state—as is to be expected. We note, however, that the mean ice thickness using the BBM is slightly higher, and that this behavior can be reproduced with the mEVP by increasing the h₀ parameter of the Hibler (1979) two-category ice formation scheme. This shows that more ice is produced in leads using the BBM—which is also to be expected as that model clearly produces more openings (Figure 3). A plausible mechanism for this is that more ice is produced in a lead that opens, refreezes, and then closes mechanically, than would have been produced under level ice. A lead can only open if ice is either being ridged or exported down-stream, so this will also act to increase the mean ice thickness, except in the vicinity of export gates, such as the Fram Strait.

The difference between BBM and mEVP is much greater if we use the ice thickness scheme of Rampal et al. (2019), who added a dynamically inert thin, or young ice class (not shown). The role of ice formation in leads is, therefore, most likely underestimated using only the two categories of Hibler (1979) in this context, but further investigation of this is outside the scope of this paper.

In addition to proposing a new constitutive model, we here also propose a new relationship between the viscosity and sea-ice concentration in Equation 10. We introduced this change because with the original formulation of Dansereau et al. (2016) low-concentration ice behaved in a more rigid-like manner than what is readily observed. This was particularly evident in the Fram Strait and along the East Greenland coast where we saw arching during summer in the Fram Strait and the ice in the East Greenland Current was too loose and did not flow as close to the coast as can be seen in observations.

The original viscosity formulation of Dansereau et al. (2016) (who use e^−C(1−A), instead of e^{−Cα(1−A)}) is only an educated first guess when it comes to the relationship between viscosity and concentration (as they themselves point out). Our reformulation is motivated by the fact that the original formulation gives too viscous ice at low concentration, as well as the idea that there should be a relationship between damage and concentration, as for instance waves are more likely to break the ice into small floes where ice concentration is low (Boutin et al., 2021; Williams et al., 2017). Our equation for η can be re-written as $$\eta ={\eta }_{0}{\left[(1-d){e}^{-C(1-A)}\right]}^{\alpha }$$ to underline this connection.

Although our formulation gives reasonably good results, the connection between damage, floe-size distribution, and concentration should be investigated in substantially more detail still. One reason for further investigation is that the theoretical basis for the current formulation is probably weak and an in-depth study of the transition between the collisional and continuum regimes should yield a much better justified formulation. Another reason is that we have seen that the formulation of the relationship between viscosity and concentration affects the PDF of convergence (Figure 5), and the convergence PDF is still not as well reproduced by our model as the shear and divergence PDFs. There is, therefore, clearly room for improvement here, from both a theoretical and practical point of view. A possible way forwards here is to build on the work of Hibler (1977); Shen et al. (1986); Feltham (2005) who derive equations for the flow of ice in the marginal ice zone that resemble those of a viscous fluid. This could lead to a more realistic formulation of Equation 10 for the limits d → 0 and A → 0.

A final point to make is that of the numerical performance of the proposed system. In practical terms then the neXtSIMv2 implementation of mEVP and BBM differs only in the calculation of σ. The BBM routine to calculate σ is longer and more complex than the mEVP routine (about 65 lines vs. about 45 lines, with more loops) and takes about 4 times the time to execute. In the neXtSIMv2 implementation this means that solving the momentum equation using BBM takes about 25% longer than it takes using mEVP, when both use 120 sub-cycling steps in our 10 km resolution setup with a model time step of 900 s.

One way to speed up the BBM execution is to reduce the undamaged elasticity, E₀, which allows for a longer time step, or fewer sub-cycling steps (as per Equation A8). Reducing E₀ to quarter of the value used so far allows us to double the dynamical time step, or halve the number of sub-cycling steps. This makes the BBM 20% faster than mEVP. Reducing E₀ even further reduces the stability of the system, but we did not attempt to pinpoint the numerically optimum value for E₀ further. Reducing E₀ this way does not reduce the quality of the results presented in here, but we have yet to fully explore the effect of reducing E₀.

Summary and Conclusions

In this paper we present a new rheology and an accompanying numerical framework for large-scale sea-ice modeling. We refer to this rheology and framework as the brittle BBM. The BBM is a further development of the elasto-brittle (EB) and Maxwell-elasto-brittle (MEB) rheologies that have been used to simulate sea ice previously in large-scale models. The main motivation behind this new development is twofold: First, to address the missing physics in the MEB rheology related to the convergence mode of deformation, and that was responsible for allowing both unrealistic local (ridges) and basin-scale thickening of the sea ice cover over time. Second, to reduce the high numerical cost associated with the semi-implicit solver used for MEB in the neXtSIM model so far.

Following the work presented in this paper we can conclude the following:

• The BBM rheology provides a good distribution of deformation magnitude and temporal variability of the highest deformation rates. The maps of deformation rates are very realistic with sharp, well localized (down to the model grid scale) features.

• Using the BBM rheology we can simulate a realistic spatial ice thickness distribution and temporal evolution.

• Using an explicit solver to solve the underlying equations delivers numerical performance similar to that of the (m)EVP rheology.

A AppendixStability Analysis

We perform a von-Neumann stability analysis for the 1D case. We presume the motion and spatial variation only to happen in the x-direction, the coefficients to be constants and all forcing to be represented by τ. In 1D, the contribution of the elastic-stiffness tensor reduces to $$\mathbf{K}:{\dot{\varepsilon }}^{n}={\partial }_{x}{u}^{n-1}$$. Abbreviating σ = σ₁₁ and $${D}^{-1}=\dot{d}/(1-d)$$, and assuming h to be constant, the discretized equations (Equation 33 including the damage term as in 20, and the sea-ice momentum Equations 44 and 45) in 1D read [Image Omitted. See PDF] [Image Omitted. See PDF]with $$\chi {:=}{\left(1+\frac{{\Delta}t}{\lambda }\left(1+\tilde{P}\right)+\frac{{\Delta}t}{D}\right)}^{-1}$$. Given that $$-1\le \tilde{P}\le 0$$ (see Equation 7b), we always have χ ∈ (0, 1].

Assuming χ to be constant in x-direction, we eliminate σ from Equation A1–A2. Therefore, we first take the spatial derivative of Equation A2 to get an explicit representation of ∂σⁿ⁺¹/∂x: [Image Omitted. See PDF]replace this expression in Equation A1 and use Equation A1 at the previous time step to derive at [Image Omitted. See PDF]with $$\psi {:=}k{\Delta}t\sqrt{E/\rho }\in (0,\pi ]$$ and −k² being the eigenvalue of $${\partial }_{xx}^{2}$$ with k² ≤ π²/Δx². With the elastic wave speed $${c}_{E}{:=}\sqrt{E/\rho }$$ and the elastic timescale, which is equal to the damage propagation time t_d ≔ Δx/c_E, we have ψ = (Δxk)Δt/t_d.

To derive a formal stability condition, we study the amplification factor ξ = uⁿ⁺¹/uⁿ. The homogeneous Equation A4, where the forcing $$\frac{\tau {\Delta}t}{h\rho }(1-\chi )$$ is ignored, can be reformulated as: [Image Omitted. See PDF]which has the solutions [Image Omitted. See PDF]The formal stability constraint reads |ξ| ≤ 1, but bearing in mind that the underlying set of equations is highly nonlinear and in order to have a stable algorithm, the stronger constraint |ξ| < 1 should hold. The angle, ω, of ξ = |ξ| exp(iω) should also be sufficiently small to resolve oscillations that may occur during the time-stepping process (see also Kimmritz et al., 2015). For instance, ω = π/2 would provoke a change in sign in every second time step. Thus ω should ideally satisfy ω ≪ π/2. Figure A1 shows both, the maximum magnitude, max |ξ_1,2|, and the maximum angle, max(ω_1,2), in the χ, ψ space for the limits k = Δx⁻¹.

Figure A1. Stability regions of the simplified 1D case in the {χ, ψ}-plane. Contour lines show the maximum angle ω of ξ1,2 between 0 and π/2 and for π. The coloring depicts max |ξ1,2|, with max |ξ1,2| > 1 shaded gray. The dotted cyan lines are the functions ψ=χ−1−1 $\psi =\sqrt{{\chi }^{-1}}-1$ (where max(ω1,2) = 0) and ψ=χ−1+1 $\psi =\sqrt{{\chi }^{-1}+1}$ (where max(ω1,2) = π/2).

The values for max |ξ_1,2| and max(ω_1,2) fall into three main regions (see Figure A1):

The first region (gray area) collects unstable solutions where max |ξ_1,2| > 1. Solutions in this area occur, when a too large time step Δt fails to properly resolve the stress redistribution of undamaged or slightly damaged ice, or ice in or very near the elastic regime $$\left(\tilde{P}\approx -1\right)$$.

The second region (yellow lower left area) contains stable solutions with |ξ_1,2| close to 1 and no phase ω_1,2 = 0. It is characterized by $$\psi < \sqrt{{\chi }^{-1}}-1$$ (lower dotted cyan curve in Figure A1). In this case, the time step is small enough to resolve the stress redistribution without any phase changes in ξ, but error damping remains very small.

Solutions in the third region, lying between these two other regions in the {χ, ψ} plane, are stable and show faster damping of the error compared to solutions located in the lower left corner. They are, however, oscillatory as ω_1,2 > 0. Here the angles ω_1,2 are arranged in conjugate pairs (As in the EVP case, see Kimmritz et al., 2015), and so solutions in this third region have the real component $$\mathit{Re}\left({\xi }_{1,2}\right)=\frac{1}{2}\left(1+\chi -\chi {\psi }^{2}\right)$$ and the imaginary components $$\mathit{Im}\left({\xi }_{1,2}\right)=\pm \sqrt{\chi -{\left(1+\chi -\chi {\psi }^{2}\right)}^{2}/4}$$, resulting in max |ξ_1,2| being of the order of $$\sqrt{1/2}\left(1+\chi -\chi {\psi }^{2}\right)$$ as a conservative estimate. To ensure a stable solution we need ω < π/2, which means that ψ should be smaller than $$\sqrt{{\chi }^{-1}+1}$$ (upper dotted cyan curve in Figure A1). This condition is the most constraining when χ = 1, resulting in [Image Omitted. See PDF]This gives a global constraint on the time step Δt [Image Omitted. See PDF]

From Equation A8 we can immediately see that the stability of the BBM framework is determined by the horizontal resolution of the model and the propagation speed of damage. For practical purposes it is important to note that the time step scales with the horizontal resolution, that is, Δt ∝ Δx, and not the resolution squared, as one would expect from a purely viscous fluid. Second, the time step scales with the propagation time of damage, which in turn scales with the undamaged elasticity as $${t}_{d}\propto 1/\sqrt{E}$$. This means that one can increase the time step of the model if the elasticity is reduced, as noted in the discussion (Section 5).

B AppendixRelevance of Changes in Concentration to the Constitutive Equation

In Section 2.2.1 we derive the constitutive equations for the BBM rheology assuming that changes in concentration, A, are slow and can be ignored. This assumption can be justified by considering the full temporal derivative of E, derived from Equation 9: [Image Omitted. See PDF]to derive the time derivative of σ_E as [Image Omitted. See PDF]Now using Equation B2, together with 7 and 16, 17–18, we can derive the analogue of Equation 20 as [Image Omitted. See PDF]

If we assume the ice is not damaging, that is, $$\dot{d}=0$$, we see that for $$\dot{A}$$ to be negligible we must have [Image Omitted. See PDF]The largest values for divergence observed in the Arctic at 10 km resolution are about 10%/day, so for the inequality to hold for highly deforming ice (and with C = 20) we have [Image Omitted. See PDF]With λ = η/E and following Equations 9 and 10, the condition above holds for d ≳ 0.7 when A = 1 and A ≲ 0.7 when d = 0.

Comparing model fields of λ and divergence shows that the condition above also holds in general, in particular because damage must become quite high (≳ 0.7) before any deformation will occur. We have also implemented Equation B2, using $$\dot{A}=-\nabla \cdot \left(\vec{v}A\right)$$ in neXtSIMv2 and this gives results that are not significantly different from the ones we present in the paper's main text.

C AppendixThe mEVP Implementation

We choose to re-arrange slightly the mEVP equations in the neXtSIMv2 implementation, in order to have a more general code which requires only small changes to switch between mEVP, EVP, and MEB. In mEVP the momentum equation is generally written as (e.g., Danilov et al., 2015) [Image Omitted. See PDF]or [Image Omitted. See PDF]Here β is the mEVP damping parameter, n denotes the sub-time step number, u⁰ is the velocity before entering the sub-cycling, F_j = ∂σ_ij/∂x_i is the internal stress terms, and other terms are as before.

The right hand side of Equation C2 can be written as [Image Omitted. See PDF]With b ≔ β + 1, we now have [Image Omitted. See PDF]This is equivalent to using a modified time step [Image Omitted. See PDF]and an extra term in the equation of [Image Omitted. See PDF]With this, Equations 44 and 45 become (now using β from Hunke & Dukowicz, 1997) [Image Omitted. See PDF] [Image Omitted. See PDF]with α, β, τ_x, τ_y, and c′ as before. In the code it is trivial to switch between the normal and modified time steps and to include or not the additional term to efficiently switch between the mEVP and EVP time stepping.

D AppendixThe Effect of Using a Large Number of Sub-Iterations With mEVP

In addition to using 120 sub-iterations we also tested running the mEVP with 500 and 1,000 sub-iterations. The main impact is that with higher number of sub-iterations the deformation field becomes more localized (Figure D1), but since the number of features is very small, then the P90 value is lowered (Figure D2 and Section 3.3) and the magnitude of the three moments of the spatial scaling analysis is reduced (Figure D3 and Section 3.4). The effect of using a large number of sub-iterations on the PDFs is barely noticeable (not shown).

Figure D1. Maps of sea ice shear (day−1) for 2 February 2007 as simulated by neXtSIMv2 with the mEVP rheology and 120, 500 and 1,000 sub-iteration steps (panels a, b, and c, respectively).

Figure D2. Time series sea ice shear P90 for 2007 as simulated by neXtSIMv2 with the mEVP rheology and 120, 500 and 1,000 sub-iteration steps (N).

Figure D3. Spatial scaling analysis of total deformation fields as simulated by neXtSIMv2 with the mEVP rheology and 120, 500 and 1,000 sub-iteration steps (N). (a) Moments of the distributions of the total deformation rate ɛtot calculated at a temporal scale of 3 days and space scales varying from 10 to 1,000 km. (b) Structure functions, where β indicates the exponent of the power-law fits and q is the moment order.

Acknowledgments

This work was supported by Norges Forskningsråd (Grant Nos. 263044, 270061, 276730, and 325292), the Bjerknes Centre for Climate Research (through the AOI project), and the SASIP project which is supported by Schmidt Futures–a philanthropic initiative that seeks to improve societal outcomes through the development of emerging science and technologies. The authors would like to thank Camille Lique and the Pôle de Calcul et de Données Marines (PCDM) for providing DATARMOR computational resources, as well as acknowledging computational resources provided by UNINETT Sigma2 - the National Infrastructure for High Performance Computing and Data Storage in Norway. The authors would like to thank Mathieu Plante and an anonymous reviewer for their very helpful review of the paper. We would also like to thank Piotr Minakowski and Laurent Brodeau for their comments on the manuscript.

Data Availability Statement

The main model results used in this paper are available for download on Zenodo with a https://doi.org/10.5281/zenodo.6816983.