Research article 24 Jun 2021
Research article  24 Jun 2021
Nonnormal flow rules affect fracture angles in sea ice viscous–plastic rheologies
 ^{1}AlfredWegenerInstitut, HelmholtzZentrum für Polar und Meeresforschung, Bremerhaven, Germany
 ^{2}MARUM – Center for Marine Environmental Sciences, Leobener Str. 8, 28359, Bremen, Germany
 ^{3}Department of Atmospheric and Oceanic Sciences, McGill University, Montréal, Canada
 ^{1}AlfredWegenerInstitut, HelmholtzZentrum für Polar und Meeresforschung, Bremerhaven, Germany
 ^{2}MARUM – Center for Marine Environmental Sciences, Leobener Str. 8, 28359, Bremen, Germany
 ^{3}Department of Atmospheric and Oceanic Sciences, McGill University, Montréal, Canada
Correspondence: Damien Ringeisen (damien.ringeisen@awi.de)
Hide author detailsCorrespondence: Damien Ringeisen (damien.ringeisen@awi.de)
The standard viscous–plastic (VP) sea ice model with an elliptical yield curve and a normal flow rule has at least two issues. First, it does not simulate fracture angles below 30^{∘} in uniaxial compression, in contrast with observations of linear kinematic features (LKFs) in the Arctic Ocean. Second, there is a tight, but unphysical, coupling between the fracture angle, postfracture deformation, and the shape of the yield curve. This tight coupling was identified as the reason for the overestimation of fracture angles. In this paper, these issues are addressed by removing the normality constraint on the flow rule in the standard VP model. The new rheology is tested in numerical uniaxial loading tests. To this end, an elliptical plastic potential – which defines the postfracture deformations, or flow rule – is introduced independently of the elliptical yield curve. As a consequence, the postfracture deformation is decoupled from the mechanical strength properties of the ice. We adapt Roscoe's angle theory, which is based on observations of granular materials, to the context of sea ice modeling. In this framework, the fracture angles depend on both yield curve and plastic potential parameters. This new formulation predicts accurately the results of the numerical experiments with a rootmeansquare error below 1.3^{∘}. The new rheology allows for angles of fracture smaller than 30^{∘} in uniaxial compression. For instance, a plastic potential with an ellipse aspect ratio smaller than 2 (i.e., the default value in the standard viscous–plastic model) can lead to fracture angles as low as 22^{∘}. Implementing an elliptical plastic potential in the standard VP sea ice model requires only small modifications to the standard VP rheology. The momentum equations with the modified rheology, however, are more difficult to solve numerically. The independent plastic potential solves the two issues with VP rheology addressed in this paper: in uniaxial loading experiments, it allows for smaller fracture angles, which fall within the range of satellite observations, and it decouples the angle of fracture and the postfracture deformation from the shape of the yield curve. The orientation of the postfracture deformation along the fracture lines (convergence and divergence), however, is still controlled by the shape of the plastic potential and the location of the stress state on the yield curve. A nonelliptical plastic potential would be required to change the orientation of deformation and to match deformation statistics derived from satellite measurements.
Sea ice plays a significant role in the energy budget of the climate system and therefore has a strong influence on future climate projections. Sea ice dynamics are located primarily along narrow lines of deformation, called linear kinematic features (LKFs), where floes slide along and grind against each other. LKFs can form in divergence, creating stretches of open water or leads, or in convergence, creating piles of ice or ridges. LKFs in the Arctic sea ice cover influence the Earth system in many ways: heat and moisture exchange take place primarily over open water (Badgley, 1961), and salt rejection during ice formation in leads creates dense water and influences the thermohaline circulation (Nguyen et al., 2011, 2012; Itkin et al., 2015). Locally, the ice strength depends on the sea ice state (e.g., thickness, concentration, and damage), which in turn is affected by sea ice fracture with thermodynamic growth in opening leads and with local dynamical growth during ridge formation. One observable and quantifiable feature of LKFs in Arctic sea ice is the intersection angles between individual LKFs. The LKFs have an influence on the local ice strength, emergent anisotropy and future deformation in the pack ice, and therefore sea ice mass balance (Aksenov and Hibler, 2001). Reproducing the LKFs patterns, density, and orientation is important for accurate sea ice and climate projections at high resolution.
LKFs are ubiquitous features of granular media, and sea ice is often described as such a granular material (Overland et al., 1998; Erlingsson, 1988; Anderson, 1942; Schall and van Hecke, 2010). Similarly to the crumbling of rocks, sea ice also exhibits brittle fracture, as floes break into smaller pieces. Brittle behavior adds a level of complexity because it implies that models must represent both the dynamics of intact ice (brittle – fracture or elastic regime) and the dynamics of a fractured system (granular – friction or plastic regime) (Handin, 1969). The dominant deformation process along LKFs is shear. Sometimes this shear is associated with nonzero divergence, and this divergence along shear bands is referred to as dilatancy (Stern et al., 1995). Granular matter theory can explain the dilatancy along LKFs. In this work, we consider sea ice a granular material and focus on the dynamics of the fractured system. We use the term fracture as the failure of a compact assemblage of floes and define the fracture angle as half of the angle between intersecting LKFs.
Different rheological models assume different material behavior before and after fracture. Common sea ice rheological models are, for example, viscous–plastic (VP; Hibler, 1977), elastic–plastic (EP; Coon et al., 1974), elastic–anisotropic–plastic (EAP; Tsamados et al., 2013), or Maxwell elastobrittle (MEB; Dansereau et al., 2016). In these different rheological models, various stress–strain(rate) relationships, or constitutive equations, can be defined. In the following, we refer to models with different constitutive equations as different rheologies. We focus on the VP rheological model. A specific VP rheology is defined by a yield curve and plastic potential. The yield curve defines the stress criteria for the transition from small viscous deformations (creep) to the large plastic deformations (friction). The plastic potential determines the ensuing postfracture deformation, called the flow rule. The flow rule is normal to the plastic potential (Drucker and Prager, 1952). The plastic potential can be independent of, or equal to, the yield curve. In the latter case, the flow rule is also normal to the yield curve and is called a normal flow rule or associated flow rule. Several yield curves have been used in sea ice VP models, some with a normal flow rule (Hibler, 1979; Zhang and Rothrock, 2005) and some with a nonnormal flow rule (Ip et al., 1991; Tremblay and Mysak, 1997; Hibler and Schulson, 2000; Wang, 2007). We reiterate that two plastic models with the same yield curve but with different flow rules are referred to as two different rheologies, as they behave differently during the creation of LKFs.
The viscous–plastic rheology is an appropriate continuum rheology for modeling sea ice as a granular material because it includes (1) a yield condition for plastic deformation and (2) a flow rule that allows the representation of the divergent and convergent motion along shear lines, that is, the dilatancy observed in granular media. Continuum plastic flow models with normal or nonnormal flow rules are often used in other scientific fields to model granular geomaterials (Vermeer and De Borst, 1984; Mánica et al., 2018).
LKFs have been studied in satellite observations (Stern et al., 1995; Kwok, 2001; Schulson and Hibler, 2004; Weiss et al., 2007) and numerical models (Spreen et al., 2017; Hutter et al., 2018). In VP sea ice models, LKFs are represented as narrow zones of plastic deformation in a background field of nearly undeformed ice (viscous creep) (Hutchings et al., 2005). This behavior has been argued to be the reason for low temporal intermittency and spatial localization in VP models, leading to a spatial and temporal scaling of LKFs that is different from that of observations (Rampal et al., 2016). LKFs emerge clearly in plastic flow models at high resolution (Hutchings et al., 2005; Hutter et al., 2018; Koldunov et al., 2019). VP models reproduce observed intermittency and spatial localization even without brittle fracture dynamics (Bouchat and Tremblay, 2017; Hutter et al., 2018), albeit at higher resolution than Maxwell elastobrittle models (e.g., Rampal et al., 2019).
New models have been designed to represent sea ice fracture, for example, brittle models with a damage parameter that keeps the memory of previous fracture (Dansereau et al., 2016; Girard et al., 2011) or anisotropic viscous–plastic rheologies models (Tsamados et al., 2013; Heorton et al., 2018). Still, as of today, the viscous–plastic rheology with an elliptical yield curve and normal flow rule (Hibler, 1979) is the de facto standard rheology in global climate models. For example, of the 33 global climate models of the Climate Model Intercomparison Project 5 (CMIP5), 30 use the VP rheology with an elliptical yield curve and normal flow rule (Stroeve et al., 2014). Below, we refer to this rheology as the standard VP rheology.
The orientation of LKFs is a wellstudied subject in the field of engineering and granular materials (LKFs are called shear bands in this field). Two classical solutions coexist and set two limit angles for the orientation of fractures: the Coulomb angle (static behavior) and the Roscoe angle (dynamic behavior). The Coulomb angle of fracture θ_{C} between the fracture line and the first principal stress is determined by the Mohr–Coulomb criterion. It is a function only of the internal angle of friction ϕ (Coulomb, 1776; Mohr, 1900):
Roscoe (1970) challenged this concept by considering the case of dilatant material and found from experiments with sand that the dilatancy angle δ is the main parameter determining the orientation of shear bands (see Fig. 6 in Tremblay and Mysak, 1997, for a definition of the dilantancy angle δ in the context of sea ice modeling.). Dilatancy refers to divergence along shear bands or LKFs. This divergence is a function of the distribution of contact points between individual floes at the subgrid scales. A positive angle of dilatancy is associated with contact points that (on average) oppose the macroscopic shear motion and create divergence along the shear band, while negative dilatancy is associated with a closing of the shear line (ridging in the case of sea ice). The Roscoe angle of fracture is defined as
A general theory derived from experiments with sand that takes into account both the angle of friction and the angle of dilation combines the Coulomb and Roscoe angles as (Arthur et al., 1977; Vardoulakis, 1980)
where θ_{A} is called the Arthur angle. Tremblay and Mysak (1997) used this general theory to design their sea ice rheology. Vermeer (1990) proposed a theoretical framework based on the grain size and showed that the angle of fracture in most experiments falls between the two extremes: ${\mathit{\theta}}_{\mathrm{C}}\le \mathit{\theta}\le {\mathit{\theta}}_{\mathrm{R}}$, with δ<ϕ in sands. If ϕ=δ, then ${\mathit{\theta}}_{\mathrm{R}}={\mathit{\theta}}_{\mathrm{C}}={\mathit{\theta}}_{\mathrm{A}}$, and the flow rule is normal to the yield curve. In other words, for ϕ=δ, the principal axes of stress and the principal axes of strain are coaxial. This condition, however, is not generally satisfied for granular materials: experiments with sand have shown differences between ϕ and δ of the order of 30^{∘} (Balendran and NematNasser, 1993; Vardoulakis and Graf, 1985; Bolton, 1986). Note that both mechanisms, friction and dilatancy, are not radically different: a larger dilatancy angle implies a larger grain size and more contact normals opposing the flow, hence more friction (Vermeer, 1990). The concept of the internal angle of friction can be used to link the orientation of LKFs to plastic rheologies with a normal flow rule (Ringeisen et al., 2019, Appendix B). So far, only the yield curve has been considered when investigating the orientation of LKFs in the viscous–plastic model (as in, e.g., Hibler and Schulson, 2000; Hutchings et al., 2005; Wang, 2006). Therefore, it is unknown which of the three angles (Coulomb, Roscoe, Arthur) provides the most accurate prediction for this case.
The fracture angles with the standard VP rheology cannot be smaller than 30^{∘} in uniaxial compression, even by changing the ellipse aspect ratio e (Ringeisen et al., 2019). In contrast, observations show fracture angles generally below 30^{∘} (e.g., 14^{∘} in Marko and Thomson, 1977; $\mathrm{15}\pm \mathrm{1.5}{}^{\circ}$ in Erlingsson, 1988; 17 to 18^{∘} in Cunningham et al., 1994) and a clear peak in the distribution of angles between 20 and 25^{∘} (Hutter and Losch, 2020). In addition, uniaxial loading compression experiments with lateral confinement (achieved via the addition of thinner ice surrounding the ice slab) presented in Ringeisen et al. (2019) showed that (1) the angle of fracture is a function of the slope of the yield curve in stress invariant space, (2) the ellipse aspect ratio determines the divergence along the LKFs, and (3) the fracture angle is a function of the confining pressure. These three properties of the standard VP rheology do not agree with the theory and observations of granular media behavior, namely that shear band orientations and divergent or convergent motion at the slip lines are a function mainly of the shear strength of the material and orientation of the contact normals (or dilatancy angle) and that the confining pressure has only a limited effect (Balendran and NematNasser, 1993; Alshibli and Sture, 2000; Han and Drescher, 1993; Desrues and Hammad, 1989). Note that some of these last experiments are triaxial tests, and that biaxial tests of 2D granular sea ice might yield different results as sea ice can “escape” in the vertical direction. Biaxial tests on sea ice samples show that small confinements lead to coulombic shear fault fractures with a similar internal friction coefficient and similar fracture angle. However, larger confinements lead to a spalling raftlike behavior with a broader range of fracture angles (Schulson et al., 2006a). The fracture angles are similar in different regions of the Arctic with different background stress conditions (Erlingsson, 1988; Marko and Thomson, 1977; Cunningham et al., 1994). This observation supports the hypothesis that the angle of fracture in sea ice is independent of the confining pressure. Finally, the distribution of intersection angles simulated with the standard VP rheology at high resolution has a peak of the distribution at larger angles than the RADARSAT Geophysical Processor System (RGPS) dataset (Hutter et al., 2019; Hutter and Losch, 2020). The unphysical behavior of the standard VP rheology (points 1, 2, 3 presented above and the distribution of angles) is connected to the shape of the yield curve in conjunction with a normal flow rule.
The flow rule has the advantage that it can be observed with remote sensing methods, in contrast to observing stress which requires in situ measurements. The ratio of shear to divergence along the shear bands or LKFs allows the inference of the dilatancy angle of granular material. Observations of sea ice drift in the Arctic show that most of the deformation takes place in shear with some divergence (Stern et al., 1995). The distribution of the ratio of divergence and convergence can be reproduced by modifying the ellipse aspect ration e of the standard VP rheology (Bouchat and Tremblay, 2017). Separating the link between the fracture angle and the flow rule from the yield curve is necessary to design VP rheologies that are consistent with observed sea ice deformation.
In this paper, we investigate the effects of a nonnormal flow rule on fracture angles. We use the nonnormal flow rule as a means of separating the state of stress (at failure) and the postfracture deformation. To this end, we study the nonnormal flow rule in the context of the standard VP rheological model using a similar shape for the plastic potential (i.e., an ellipse) because (1) the ellipse is widely used in the community and (2) its behavior is well documented (compared to other models), providing a solid basis for comparison. For these two reasons, we use the elliptical yield curve despite the fact that it is not the most appropriate yield curve to model sea ice as a granular material. This paper provides a new generalized theoretical framework for any viscous–plastic material with normal or nonnormal flow rules. Following Ringeisen et al. (2019), we test the new model in simple uniaxial loading experiments where the relationship between fracture angle and flow rule can be easily identified.
The paper is structured as follows. Section 2 describes the model (Sect. 2.1), the new rheology (Sect. 2.2), and a general theory linking the fracture angles and a general flow rule (Sect. 2.3). Sects. 3 and 4 describe the idealized experimental setup and the results. Section 5 discusses these results and their implications for current and future rheologies. Conclusions follow in Sect. 6.
2.1 Building the sea ice VP constitutive equations
We consider sea ice a 2D viscous–plastic material. The ice velocities are calculated from the sea ice momentum equations:
where ρ is the ice density, h is the gridcellaveraged sea ice thickness, u is the ice drift velocity field, f is the Coriolis parameter, k is the vertical unit vector, τ_{a} is the surface air stress, τ_{o} is the ocean drag, ∇ϕ_{s} is acceleration from the gradient of sea surface height, and σ is the vertically integrated internal ice stress tensor defined by the sea ice VP constitutive equations. The constitutive equations define the vertically integrated stress tensor σ as a function of the strain rate tensor $\dot{\mathit{\u03f5}}$ and the state variables χ (e.g., ice thickness, ice strength, ice concentration). The components of the strain rate tensor are computed from the velocities as ${\dot{\mathit{\u03f5}}}_{ij}=(\partial {u}_{i})/(\partial {x}_{j})$. The constitutive equations then have the following form:
In the sea ice VP model the stresses are independent of the strain rates for large deformation events (the plastic states with stresses on the yield curve) and they depend on the strain rates for small deformations (the viscous states with stresses inside the yield curve). It is this set of equations that defines the rheology of sea ice and determines the fracture pattern and the opening or closing along the fractures.
One of the state variables in the model is the maximum compressive strength P. This variable represents the maximum compressive stress that sea ice can bear in uniform compression before ridging. We use the following simple standard relationship (Hibler, 1979):
where C^{⋆} is a free parameter (typically ${C}^{\star}=\mathrm{20}$), h is the mean ice thickness, A is the fractional sea ice area cover in a grid cell, and P^{⋆} is the ice strength of 1 m of ice at 100 % concentration (A=1).
The yield curve represents the stress states for which sea ice deforms plastically while enclosing the stress states for which sea ice slowly deforms viscously. We express the yield curve as a function of the stresses σ_{ij} and the state variables χ:
The yield curve can be represented in principal stress (σ_{1} and σ_{2}) or stress invariant space (σ_{I} and σ_{II}). Figure 1 shows an arbitrary yield curve in stress invariant space. Although Eq. (7) determines if the deformation is plastic or viscous, it does not determine how the ice will deform after fracture. In order to obtain a closed system of equations, we define a plastic potential that defines the flow rule.
The plastic potential determines the direction of deformation for stress states on the yield curve. The flow rule represents the direction of deformation in the grid cell. The orientation of the flow rule in the coordinate system (${\dot{\mathit{\u03f5}}}_{\mathrm{I}}$, ${\dot{\mathit{\u03f5}}}_{\text{II}}$), as shown in orange in Fig. 1, indicates if the grid cell deforms in convergence (${\dot{\mathit{\u03f5}}}_{\mathrm{I}}<\mathrm{0}$) or divergence (${\dot{\mathit{\u03f5}}}_{\mathrm{I}}>\mathrm{0}$) and shear (${\dot{\mathit{\u03f5}}}_{\text{II}}$). Just as with the yield curve, the plastic potential can be written as
The direction of the deformation, called the flow rule, is perpendicular to the plastic potential. This is shown in red in Fig. 1b and mathematically expressed by
where λ>0 is the unknown flow rate. The flow rule is applied for stress states on the yield curve at the same compressive stresses (orange arrows in Fig. 1b). If the plastic potential and the yield curve are the same (G=F), the flow rule is called an associated or normal flow rule, as the flow rule is also perpendicular to the yield curve (see Fig. 1a).
Using Eqs. (7) and (9), we can write a system of five equations (four from Eq. 9 and one from Eq. 7) for five unknowns (σ_{11}, σ_{22}, σ_{12}, σ_{12}, λ). Solving this system of equations allows us to write the constitutive equations for the sea ice model as function of the components of the strain rate (${\dot{\mathit{\u03f5}}}_{\mathrm{11}}$, ${\dot{\mathit{\u03f5}}}_{\mathrm{22}}$, ${\dot{\mathit{\u03f5}}}_{\mathrm{12}}$, ${\dot{\mathit{\u03f5}}}_{\mathrm{21}}$) and the state variables χ.
After deriving these constitutive equations, we assume that the stress and strain rate tensors are symmetric; that is, σ_{12}=σ_{21} and ${\dot{\mathit{\u03f5}}}_{ij}=\mathrm{1}/\mathrm{2}\left(\right(\partial {u}_{i})/(\partial {x}_{j})+(\partial {u}_{j})/(\partial {x}_{i}\left)\right)$. The symmetry follows from ignoring the rotation in an isotropic medium. Note that we first need to solve this system of equations without using the symmetry condition; the symmetry condition is only invoked at the end. Applying the symmetry before solving the system of equation would change the nature of the initial tensor, and the resulting constitutive equations would be different.
An ideal plastic model, with the stresses independent of the strain rates, has a singularity because the nonlinear viscosities tend to infinity as the strain rates tend to zero. Hibler (1977) solved this issue with a regularization that limits the value of the bulk and shear viscosities ζ and η to a maximum value. When the viscosities are capped to their maximum values, the stresses are linearly related to the strain rates and the material behaves as a viscous material. VP sea ice models typically cap the viscosity at
and ${\mathit{\eta}}_{\text{max}}=\left({\mathit{\zeta}}_{\text{max}}\right)/{e}_{G}^{\mathrm{2}}$ to regularize the momentum equations. When this regularization is in effect, ζ and η are independent of the deformation field (Δ) and the stress divergence reduces to harmonic viscosity with constant coefficients. ${\mathrm{\Delta}}_{\text{min}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{9}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$ (Hibler, 1979, 1977) translates to a deformation timescale of almost 16 years. Therefore, viscous deformations are slow and negligible with respect to the plastic deformations that operate on a synoptic timescale, and sea ice VP rheologies can be considered ideal plastic. The viscous behavior can be seen as a consequence of regularizing the viscosities rather than of an implementation of a type of physical behavior.
2.2 Elliptical yield curve with nonnormal flow rule
We now build a rheology with an elliptical yield curve and a nonnormal flow rule; that is, we use a plastic potential G that is different from the yield curve F. By doing this, we will change the orientation of the flow rule, without changing the yield stress state (see Figs. 2 and 4 for some examples). We use a different, but still elliptical, plastic potential for simplicity; this choice requires only minor modifications to a typical VP sea ice model. We define the yield condition F and the plastic potential G as a function of the state variables χ – the ice compression strength P, the ice tensile strength T=k_{t}P (König Beatty and Holland, 2010), the yield curve's ellipse ratio e_{F}, and the plastic potential's ellipse ratio e_{G} – by
for $X=F,G$ for the yield curve or the plastic potential. Using Eq. (11), we write σ_{II} as a function of σ_{I} as
Following Hibler (1977, 1979), we derive the constitutive equations σ_{ij}:
where the shear and bulk viscosities η and ζ are defined by
with
Figure 2 shows an example of a yield curve and plastic potential, with the resulting flow rule. For e_{G}>e_{F}, the absolute value of the divergence is smaller and the shear strain rate is larger compared to a normal flow rule (e_{G}=e_{F}) and vice versa for e_{G}<e_{F}.
2.3 Linking fracture and flow rule
In this section, we generalize the theory linking the rheological model and the fracture angles in a simple uniaxial compressive test (Ringeisen et al., 2019) to materials with a nonassociated flow rule. To this end, we follow the theory of Roscoe (1970), where the angle of fracture depends uniquely on the angle of dilatancy of a granular material. Based on laboratory experiments, Roscoe (1970) states that the velocity characteristics (the postfailure deformation) seem to be a better predictor than the stress characteristics (the stress at failure) for the orientation of shear bands in granular materials. The Roscoe angles can then be compared to the Coulomb angles, as defined in Ringeisen et al. (2019), and to the results from the idealized experiments in Sect. 4.
Figure 3 illustrates the case of an arbitrary yield curve with an arbitrary plastic potential. To adapt the Roscoe angles to sea ice modeling, we proceed as follows: (1) the stress state on the yield curve (point p in Fig. 3a) defines the position and size of Mohr's circle at the fracture (blue circle in Fig. 3b) and (2) the slope of the plastic potential determines the point on the Mohr's circle where deformation takes place; that is, the slope directly predicts the fracture angle θ as a function of the dilatancy angle δ (as per Roscoe theory, Fig. 3b). For the special case of uniaxial compression, we (A) determine the stress state on the yield curve for uniaxial compression as a function of the yield curve ellipse ratio e_{F} and (B) compute the slope of the plastic potential at that stress state as a function of the plastic potential ellipse ratio e_{G}. Finally, we combine (2) and (B) to compute the theoretical prediction for the fracture angle as a function of ellipse ratios e_{G} and e_{F}. Figure 3 shows the geometrical construction that links the angle of dilatancy δ to the slope of the plastic potential tan (γ_{G}):
Note that the minus sign above was included in the derivative of the yield curve function in Eqs. (B1) and (B2) of Ringeisen et al. (2019). This equation agrees with the definition of Roscoe (1970) that $\mathrm{sin}\left(\mathit{\delta}\right)={\dot{\mathit{\u03f5}}}_{\mathrm{I}}/{\dot{\mathit{\u03f5}}}_{\text{II}}$, because the ratio of ${\dot{\mathit{\u03f5}}}_{\mathrm{I}}$ to ${\dot{\mathit{\u03f5}}}_{\text{II}}$ is equal to the slope of the plastic potential $\partial {\mathit{\sigma}}_{\text{II},G}/\partial {\mathit{\sigma}}_{\mathrm{I}}$, as the flow rule is perpendicular to the plastic potential. Figure 3 also shows the normal flow rule, which, in agreement with the coulombic theory, would lead to different fracture angles (light blue lines). From Fig. 3, the fracture angle can be written as
Substituting Eq. (16) in the equation above, the relationship between the fracture angle and the plastic potential becomes
We calculate the fracture angles for the elliptical yield curve with a nonnormal flow rule in uniaxial compression along the y axis. In this case, ${\mathit{\sigma}}_{\mathrm{11}}={\mathit{\sigma}}_{\mathrm{12}}=\mathrm{0}$; σ_{22}<0; and the principal stresses and stress invariants can be written as
From Eq. (23), the maximum shear stress ${\mathit{\sigma}}_{\text{II},F}^{p}$ in the fracture plane in uniaxial compression can be expressed as
where p indicates the stress state at the fracture. Figure 4 shows the stress trajectory in principal stress space for uniaxial compression. It also shows how the flow rule changes for the same stress state when using two different elliptical plastic potentials.
In the following, we use the normalized stress invariants ${\mathit{\sigma}}_{\mathrm{I}}^{\prime}={\mathit{\sigma}}_{\mathrm{I}}/P$ and ${\mathit{\sigma}}_{\text{II}}^{\prime}={\mathit{\sigma}}_{\text{II}}/P$ to simplify the notation. The slope of the yield curve or the plastic potential depends only on e_{F} and e_{G} and not on P. Substituting Eq. (12), ${\mathit{\sigma}}_{\mathrm{I}}^{\prime}$, and ${\mathit{\sigma}}_{\text{II}}^{\prime}$ in Eq. (24), we obtain
and solve the first stress invariant ${\mathit{\sigma}}_{\mathrm{I}}^{p}$ on the fracture plane in uniaxial compression:
The slope of the tangent at ${\mathit{\sigma}}_{\mathrm{I}}^{p}$ to the plastic potential is given by the derivative of Eq. (12):
Substituting Eq. (26) into Eq. (27) yields
or, for zero tensile strength (k_{t}=0),
The fracture angle can finally be written as a function of e_{G} and e_{F} from Eq. (18):
As expected, for ${e}_{F}={e}_{G}=e$, we recover the fracture angle derived in Ringeisen et al. (2019):
Following Ringeisen et al. (2019), we load a rectangular ice floe of 8 km by 25 km with a uniform thickness of h=1 m and a uniform sea ice concentration of A=100 % (see Fig. 5). The numerical domain has the dimensions L_{x}=10 km and L_{y}=25 km. At y=0, we use a closed, solid boundary with a noslip condition (i.e., $u=v=\mathrm{0}$). At x=0 and L_{x}, we use Neumann boundary conditions:
On the left and right sides of the domain (x<1 km and x>9 km), we have open water between the ice floe and the boundary to ensure that the boundaries have no effect on the simulation. At y=L_{y}, we use a Dirichlet boundary condition for ice velocity (v the velocity in the y direction increasing linearly in time simulating an axial loading test) and a Neumann boundary condition for ice thickness and concentration:
where a_{v} is $\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{2}}$. The grid spacing of the domain is 25 m, and the time step is 0.1 s. For simplicity, the Coriolis parameter f is 0.
The nonlinear momentum equation, Eq. (4), is integrated using a Picard solver with 15 000 nonlinear (or outerloop) iterations (Losch et al., 2010). For the linearized problem within each nonlinear iteration, we use a line successive (over)relaxation (LSR) method (Zhang and Hibler, 1997), with a tolerance criterion of ${u}_{k}{u}_{k\mathrm{1}}{}_{max}<{\mathrm{10}}^{\mathrm{11}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, where k is the linear iteration index. We use an inexact approach with only a maximum of 200 linear iterations for the linearized equations; the linearized system does not reach the tolerance criterion for the first nonlinear iterations, but it does so as the nonlinear system approaches a converged solution. We chose a very small tolerance and residual norm for the solution of the linear and nonlinear problem in order to simulate a clean fracture with a welldefined fracture angle – for comparison with theory and observations. These criteria are much stricter than common recommendations for Arctic sea ice simulations (e.g., Lemieux and Tremblay, 2009). We expect that numerical sea ice models are computationally more challenging with a nonnormal flow rule than with a normal flow rule. The nonnormality of the flow rule relative to the yield curve introduces more complexity because Drucker's stability postulate is not satisfied (Vermeer and De Borst, 1984; Balendran and NematNasser, 1993). This particular uniaxial loading experiment is also complex to solve numerically because the forcing is localized on the boundary, in contrast to real geophysical system integrations where wind and ocean currents are acting over the entire surface of the ice.
The intersection angles between the LKFs are measured manually with the Measure Tool from the GNU Image Manipulation Program (GIMP; version 2.8.16; https://www.gimp.org/, last access: 3 June 2021). We estimated the accuracy as $\pm \mathrm{1}{}^{\circ}$ (Ringeisen et al., 2019). The first 5 s of simulations is used to define the sea ice fracture and calculate the fracture angle. Although the forced deformation is very slow, the stresses already reach the yield curve in the first time step (0.1 s). The fracture is created immediately, but because of the large viscosity of the viscous states with a deformation timescale of approximately 35 years, the fracture progression is not visible immediately (Ringeisen et al., 2019). Therefore we show the deformation after 5 s. During this 5 s there is no fundamental change other than the initial deformation becoming clearer. The angle of each fracture line is measured and used to compute the average fracture angle and the standard deviation. Note that the fracture angles do not depend on resolution, scale, geometry, or boundary conditions (see Ringeisen et al., 2019, their Sect. 3.2). We do not use a replacement pressure scheme (Ip et al., 1991; Ip, 1993), because it has no influence on the angle of fracture (not shown).
We study the evolution of the fracture angle θ when the plastic potential changes while the yield curve stays the same (see Fig. 4 for details). In this manner, the ice breaks for the exact same stress state but with a different flow rule. For simplicity, we test here the elliptical yield curve without tensile strength (k_{t}=0).
Figure 6 shows the fracture pattern for the standard yield curve ellipse ratio e_{F}=2.0 and three values of the plastic potential ellipse ratio: e_{G}=1.4, 2.0, and 4.0. The fractures form a diamond shape, which is similar to the shapes observed at large scales (Erlingsson, 1988) and in laboratory experiments (Schulson, 2001) and modeled with DEMs (Wilchinsky et al., 2010) or other continuous sea ice models (Ringeisen et al., 2019; Heorton et al., 2018). With a normal flow rule (e_{G}=2.0), single pairs of fracture lines with one unique fracture angle, large deformation along the LKFs, and smaller deformations (by several orders of magnitude) within diamondshape floes are simulated. With the nonnormal flow rule (e_{G}=1.4 and e_{G}=4) we make three observations:

Asymmetric secondary fracture lines appear, in contrast to the normal flow rule simulation. We attribute the asymmetry and presence of secondary fractures to the lack of full numerical convergence associated with the violation of Drucker's principle, or the nonnormality of the flow rule (the ratio of divergence to shear strain rate differs from that of the shear to normal stress). For instance, the L_{2} norm of the residual (R) in the nonlinear equations decreases by 4 orders of magnitude for the normal flow rule compared with 2 orders of magnitude for the nonnormal flow rule for the same number of nonlinear iterations (15 000), specifically to $R=\mathrm{8}\times {\mathrm{10}}^{\mathrm{4}}$ for ${e}_{G}={e}_{F}=\mathrm{2}$ and to $R=\mathrm{6}\times {\mathrm{10}}^{\mathrm{2}}$ for e_{G}=1.4 and e_{F}=2. Note that a Jacobianfree Newton–Krylov (JFNK) solver with a quadratic local numerical convergence does not perform better because the global convergence is poor with a combination of localized forcing and a high grid resolution (Losch et al., 2014; Williams et al., 2017).

The width and activity of the LKFs is also affected by the flow rule. With e_{G}=1.4, the lines are thinner, the shear along the LKFs is smaller, and there is little shear between the fracture lines. With e_{G}=4.0, the fracture lines are broader, the shear strain rate along the LKFs is higher, and there is more shear between the fracture lines. With e_{G}=1.4, the deformation at the fracture is mainly in divergence, while for e_{G}=4.0, the deformation is mainly in shear and there is more stress transmitted to the ice in between the fracture lines.

The fracture angle changes as the plastic potential changes. The angles are wider with e_{G}=4 than e_{G}=1.4. The effect of flow rule orientation on the fracture angles is discussed below.
We now present results from four sets of simulations with fixed yield curve ellipse ratios at e_{F}=0.7, 1.0, 2.0, and 4.0. For each of these, we test the sensitivity of the results to changes in the plastic potential ellipse ratio e_{G}. The choices of yield curve ellipse ratios e_{F} are the standard value of Hibler (1979), values suggested by Bouchat and Tremblay (2017) and Dumont et al. (2009), and an extreme value resulting in a very small shear strength and smaller fracture angles.
Figure 7a shows how the fracture angles evolve as the plastic potential ellipse ratio e_{G} changes for each of the four values of e_{F}. There is a clear dependence of the fracture angles on the relative eccentricity of the plastic potential and yield curve. For e_{G}>e_{F}, the shear strain rate increases along the LKFs (see Fig. 6c), and the fracture angles tend toward 45^{∘} as e_{G} increases, in agreement with the theory (Eq. 18). For e_{G}<e_{F}, the flow rule implies more divergence (for e_{F}>1 or convergence for e_{F}<1) and less shear along the LKFs (see Fig. 6a), and the fracture angles move away from 45^{∘} as e_{G} decreases. More generally, for e_{F}<1, the fracture angle increases with increasing convergence along the LKFs as e_{G} decreases. For e_{F}>1, the fracture angle decreases with increasing divergence as e_{G} decreases. For e_{F}=1 (a circular yield curve), the fracture angles are independent of e_{G} because the fracture takes place at the peak of the yield curve and the flow rule is not affected by changes in the plastic potential ellipse ratio (e_{G}).
The colored, dashed lines in Fig. 7a show the fracture angles ${\mathit{\theta}}_{e,nn}({e}_{F},{e}_{G})$ predicted by Eq. (30). The coefficient of determination r^{2} and the rootmeansquare error (RMSE) between the simulated angle of fracture and theoretical predictions are 0.97 and 0.37^{∘} for e_{F}=0.7, 0.95 and 1.22^{∘} for e_{F}=2.0, and 0.97 and 0.47^{∘} for e_{F}=4.0. The RMSE is 0.37^{∘} for e_{F}=1.0, with r^{2} being inapplicable. That is, the theory predicts the fracture angles accurately. This result shows that the flow rule plays a major role in the simulated fracture angle for a given rheology. The dashed black line show the change in the fracture angle with a normal flow rule (e_{G}=e_{F}; Eq. 31).
For completeness, Fig. 7b also show the theoretical predictions for a constant plastic potential ellipse ratio e_{G} for varying yield curve ellipse ratios e_{F}. The fracture angles become smaller as e_{F} increases. Yield curve ellipse ratios smaller than e_{F}=1 do not create fracture angles below 45^{∘}.
The idealized experiments using the elliptical yield curve with a nonnormal flow rule confirm that the type of deformation and the fracture angle are intimately linked with the shape of the plastic potential. We observe that, irrespective of the plastic potential elliptical aspect ratio, a yield curve ellipse ratio of e_{F}<1 does not allow fracture angles smaller than 45^{∘} in uniaxial compression. To reduce the fracture angles with yield curve ellipse ratios of e_{F}>1, plastic potential ellipse ratios e_{G} smaller than the yield curve ellipse ratio are required; that is, e_{G}<e_{F}. The idealized experiments show that with a plastic potential in a viscous–plastic model we can separate the yield criterion from the resulting deformation (flow rule). This allows decoupling the mechanical strength properties of the material (ice) from its postfracture behavior. The results illustrate clearly how the yield curve defines the stress for which the ice will deform, that is, the transition between viscous and plastic deformation, and how the relative shape of the plastic potential with respect to the yield curve defines both the fracture angle and the type of deformation (convergence or shear) along the fracture line. The resulting fracture angles are in excellent agreement with the Roscoe angle predictions (Roscoe, 1970).
Understanding the link between the rheology and fracture angle is necessary for choosing or designing a rheology that is capable of reproducing the observed intersection angles between pairs of LKFs and consequently the emerging anisotropy. An independent plastic potential may resolve several inconsistencies of the standard elliptical yield curve with a normal flow rule (discussed in Ringeisen et al., 2019), namely the following:

In the standard VP model with an elliptical yield curve and normal flow rule, adding shear strength increases the fracture angle, in contradiction to granular matter theory (Coulomb, 1776). This behavior is linked to the specific shape of the elliptical yield curve with a maximum shear stress at $P/\mathrm{2}$ and an ascending and a descending part. In principle, we can decrease the fracture angle with increasing shear strength (e_{F} decreasing) by decreasing e_{G} but only if e_{F}>1, but then the flow rule is far from normal, making the numerical convergence difficult.

Because of the elliptical shape of the yield curve, the angle of fracture in the standard VP model changes with confining pressure (Ringeisen et al., 2019, Sect. 3.2.2, Fig. 8) unlike in laboratory experiments with granular materials (e.g., sand) where the fracture angle is only weakly sensitive to the confining pressure (Han and Drescher, 1993; Desrues and Hammad, 1989; Alshibli and Sture, 2000). This behavior cannot be eliminated with an elliptical plastic potential, as the normal stress along the LKFs increases with confining pressure and the flow rule changes from divergence to convergence as one passes the maximum shear stress at $P/\mathrm{2}$ (Ringeisen et al., 2019). A different plastic potential function would change this behavior. However, this would make the model implementation and numerical convergence even more difficult. We note that a 3D granular material like sand cannot release stress by ridging as sea ice does. A 2D material, such as sea ice, can ridge and “escape to the third dimension” after fracture. Therefore, we expect a change in the fracture angles at large confinement. Laboratory experiments show this behavior and yield stresses in sea ice change above a critical confinement ratio (Golding et al., 2010; Schulson, 2002). It is still not clear whether these results can be extrapolated to the modeling of sea ice as a 2D medium at the geophysical scale, although several common features can be found (Schulson, 2002).

In the standard VP model with a normal flow rule, the divergence and convergence are set by the ellipse ratio of the yield curve and thus by the relative amounts of compressive and shear stress. The plastic potential ellipse ratio e_{G} changes the flow rule but does not change the sign of the divergence along the LKFs, which is solely determined by the yield curve ellipse ratio e_{F}. With the elliptical plastic potential, convergent motion remains convergent and only the ratio of shear to convergence changes. To change this behavior, a different shape of plastic potential is required, for example a teardrop plastic potential.

The fracture angles in the standard VP models are larger than observed. Using a nonnormal flow rule allows us to change the fracture angle in uniaxial compression to values below 30^{∘}. This is not possible with a normal flow rule (Ringeisen et al., 2019).
We discuss the elliptical yield curve here because it is the most commonly used one and its behavior is better documented than any other model in use in the community. This provides a known reference for studying the use of nonassociated flow rules. Our goal is to provide a reference for the future development of viscous–plastic rheologies with nonnormal flow rules rather than suggest a new VP rheology. Alternatives to the elliptic yield curve have been used before, for instance, the Mohr–Coulomb, the Coulombic yield curve, or the teardrop yield curves (Fig. 8). The concept of a plastic potential in conjunction with these yield curves may also prove useful in solving the issues described above. A detailed analysis of the simulations using the family of Mohr–Coulomb and teardrop yield curves is beyond the scope of this work and will be presented in a subsequent study. Below, we use the experience from our simulations to infer how alternative yield curves may address deficiencies in the standard VP rheology.
Nonnormal flow rules can be combined with the Mohr–Coulomb family of yield curves. For a Mohr–Coulomb yield curve with a doublesliding law (i.e., pure shear deformation; Ip et al., 1991), the Roscoe theory predicts a fracture angle of approximately 45^{∘} that is independent of the slope of the yield curve. This behavior can be mimicked using an elliptical yield curve and plastic potential by setting e_{G}≫e_{F}; hence δ≃0 and $\mathit{\theta}=\mathrm{45}{}^{\circ}$ (Fig. 7a). This contradicts the Coulomb theory, which predicts an angle of fracture that depends exclusively on the internal angle of friction (Eq. 1). Combining an angle of dilatancy with a Mohr–Coulomb yield curve (Tremblay and Mysak, 1997) would allow an angle of fracture depending on δ that is different with shear and divergence (δ>0) or with convergence (δ<0) along the LKFs. Such a fracture angle and divergence would be independent of the shear strength and the confining pressure in agreement with Roscoe's angle of fracture so that such a rheology could potentially solve all four issues in Sect. 5. It is also important to note that the Mohr–Coulomb yield curves do not satisfy the convexity requirements of Drucker's stability postulate. Mohr–Coulomb yield curves in plastic Earth mantle models lead to a variety of fracture angles corresponding to the Coulomb angle, the Roscoe angle, and the intermediate Arthur angles (Buiter et al., 2006; Kaus, 2010; Mancktelow, 2006). However, such geological models are designed for an incompressible medium, and making inferences for the compressible formulation of sea ice models is difficult.
The Coulombic yield curve uses the two straight limbs from the Mohr–Coulomb yield curve and an elliptical cap of the standard VP rheology for large compressive stresses (Hibler and Schulson, 2000). In this rheology, the flow rule over the two straight limbs is defined by the elliptical yield curve; that is, the ellipse serves as a plastic potential for the Mohr–Coulomb yield curve. The Coulombic yield curve leads to unrealistic and asymmetrical fracture lines (i) when the stress states fall onto the nondifferentiable intersection between the straight limbs and the elliptical cap (Ringeisen et al., 2019) and (ii) when the stress states fall onto the two straight limbs with the nonnormal flow rule. Note that straight and symmetric fracture lines in this rheology are only possible when all the stress states are on the Mohr–Coulomb limbs and the flow rule at the fracture line is nearnormal, that is, at the location where the normal to the elliptic plastic potential is nearly perpendicular to the limbs of the Mohr–Coulomb yield curve (Ringeisen et al., 2019). Hibler and Schulson (2000) have already inferred that the flow rule may have an effect on the angle of fracture, but the authors limited their case to the framework of flawed ice and did not consider Roscoe's theory of dilatancy. The rheology of Hibler and Schulson (2000) was tested in an idealized experiment that was more complex than ours (Hutchings et al., 2005), but the effect of using a nonnormal flow rule was not explored. The complexity of their setup may explain the observed difference between simulated and predicted angles. Note that the rheology in Hibler and Schulson (2000) was built by changing the shape of the yield curve a posteriori, while the rheology presented here solves the constitutive equations rigorously.
The teardrop yield curve with a normal flow rule (Rothrock, 1975; Zhang and Rothrock, 2005) is divergent for a wide range of normal stresses and for all practical purposes consists of a continuously differentiable version of the Coulombic yield curve. This asymmetry between divergent and convergent deformation for different normal stresses decreases the effect of confinement on the fracture angle – issue 2 in Sect. 5 – and reduces the fracture angle for any confinement pressure – issue 4. This yield curve does not address issue 1, because adding shear strength in the teardrop yield curve also increases the fracture angle.
As the main disadvantage of a nonnormal flow rule, we found that it leads to slower convergence of the numerical solver. Solving the momentum equation accurately requires more solver iterations, and failure to converge is more frequent than for standard normalflowrule rheologies. In our simulations, this numerical issue manifests itself by the presence of multiple and asymmetrical fracture lines despite the fact that our experiment geometry and forcing are exactly symmetrical. This asymmetry is not expected and is not found with normal flow rules. The fracture lines with a normal flow rule are symmetrical and come in pairs (Ringeisen et al., 2019). In practice, the numerical convergence issue will go unnoticed in simulations using realistic geometries and timevarying wind forcing. In these simulations, while the number of iterations typically used (O(10)) is much smaller than that required for full convergence, at each timestep, a new iteration typically uses the solution of the previous time step as the initial estimate. With this, together with slowly varying forcing in space and time, the number of solver iterations per forcing cycle is large, in contrast to the fastchanging forcing in this study (every time step). Whether this behavior (asymmetry and multiple fracture lines) will also be present in realistic simulation using spatially and temporally varying wind forcing remains to be tested.
The following criteria should be considered when building a new rheology: the spatial and temporal scaling of sea ice deformation should agree with observations (Bouchat and Tremblay, 2017; Hutter et al., 2018); the flow rule should reproduce the correct divergence along LKFs (Stern et al., 1995); the yield curve should include some tensile strength (Coon et al., 2007) and be Coulombic in nature in agreement with observed internal stress invariants from ice stress buoys (Weiss and Schulson, 2009); the distribution of fracture angles should agree with observations (Marko and Thomson, 1977; Erlingsson, 1988; Cunningham et al., 1994; Hutter et al., 2019); and the sea ice mechanical strength properties (i.e., yield curve) and deformation (i.e., flow rule for VP rheologies) should vary in time and space depending on additional variables or parameterizations, for example, the timevarying distribution of the contact normals (Balendran and NematNasser, 1993), floe size distributions (Horvat and Tziperman, 2017; Roach et al., 2018), or a damage parameter (Dansereau et al., 2016; Plante et al., 2020), as per observations and laboratory or numerical experiments (Overland et al., 1998; Hutter et al., 2019).
Although highspatialresolution observations from satellites are available from optical instruments (e.g., from the Landsat or Sentinel programs), higher temporal resolutions of sea ice deformation and flow size distributions are still unavailable. The new Sentinel constellation and in situ observations from the field program the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) may bridge this gap. There is also a knowledge gap in the interplay between yield stresses and the postfracture deformation in a 2D granular material such as sea ice. This interplay is likely different than for the wellstudied case of a solid, homogeneous, 3D block of ice (e.g., Schulson, 2002). Sea ice floating on the ocean surface can escape vertically when it forms ridges under confined compression (Hopkins, 1994). This behavior differs from laboratory tests with 3D granular material like sand that use axial symmetry. Generally, information about sea ice resistance in different configurations (e.g., confinement) and the resulting fracture angles and deformation (ridging or opening) is also still missing, although some laboratoryscale experimental results are available (Weiss et al., 2007; Schulson et al., 2006b; Weiss and Schulson, 2009). The sea ice flow size distribution varies in time (summer/winter) and space (marginal ice zone/central Arctic) (Rothrock and Thorndike, 1984). These variations change the mechanical properties (e.g., distribution of contact normals) and thermodynamic properties (e.g., lateral melt) of sea ice (Horvat and Tziperman, 2017). Designing more appropriate rheologies for improved highresolution climate models and more accurate sea ice prediction systems requires consolidated observations of these stillunclear physical processes.
The flow rule, which dictates the postfracture deformation, has a fundamental effect on the orientation of fractures lines in a viscous–plastic (VP) sea ice model. To test this, we added an elliptical plastic potential (allowing for a nonnormal flow rule) to the standard VP rheology with an elliptical yield curve, therefore modifying the flow rule without changing the yielding stress state. We tested this new rheology with numerical experiments in uniaxial compression using the standard VP model of Hibler (1979). The modeled fracture angles are in agreement with the Roscoe angle, a theory based on experiments with granular materials that includes an angle of dilatancy (Roscoe, 1970; Tremblay and Mysak, 1997). This new rheology partially solves issues raised in an earlier study (Ringeisen et al., 2019). The use of a plastic potential or nonnormal flow rule allows for the simulation of smaller fracture angles between pairs of linear kinematic features, in agreement with satellite observations. Because of the elliptical yield curve, the fracture angles still depend on the confinement pressure (Ringeisen et al., 2019), and the elliptical plastic potential only modifies the ratio of divergence relative to shear but not the direction of deformation at the fracture lines (convergence or divergence). The momentum equations for a rheological model with a nonnormal flow rule are more difficult to solve numerically and produce multiple lines of fractures that are asymmetrical (despite the geometrical symmetry of the problem), in contrast with a model with a normal flow rule. It is necessary to understand the effect of the flow rule on the fracture angle to design VP rheologies for highresolution sea ice modeling that reproduce both fracture angles and deformation along the fracture lines and the behavior of sea ice as a granular material.
Designing a rheology for highresolution simulations requires information about sea ice fracture angles and sea ice strength in a wide range of stress conditions (i.e., compression with or without confinement, pure shear, tension), yet this is unavailable at high temporal and spatial resolution. The observations of MOSAiC (Dethloff et al., 2016) in 2019–2020 may provide valuable data from continuous ice radar imaging, stress sensors, and arrays of drift buoys that will greatly help improve sea ice model dynamics.
The sea ice rheology used in this paper is implemented in the sea ice package of the MITgcm, Mid 2021 release (https://doi.org/10.5281/zenodo.4968496, Campin et al., 2021; Marshall et al., 1997; Losch et al., 2014).
No datasets were used in this article. All simulation data have been obtained with the model cited in the “Code availability” section below. The idealized model configuration is described in detail in Sect. 3.
DR designed the rheology and implemented the code changes with ML. DR ran the experiments. DR and BT designed the theory linking sea ice rheologies and granular matter theory. DR wrote the manuscript with contributions from BT and ML.
The authors declare that they have no conflict of interest.
The authors would like to thank Véronique Dansereau and Harry Heorton for their review and numerous comments that helped improve this paper. The authors also thank Jennifer Hutchings and Yevgeny Aksenov for their comments and involvement as editors. The authors are thankful to Mathieu Plante for his comments on this paper, as well as to Stephanie Deboeuf and Guillaume Ovarlez for discussions on granular materials rheologies.
This project has been supported by the Deutsche Forschungsgemeinschaft (DFG) through the International Research Training Group “Processes and impacts of climate change in the North Atlantic Ocean and the Canadian Arctic” (grant no. IRTG 1904 ArcTrain). This work is a contribution to the Natural Sciences and Engineering Research Council Discovery Grant awarded to L. Bruno Tremblay.
The article processing charges for this openaccess publication were covered by the University of Bremen.
This paper was edited by Yevgeny Aksenov and Jennifer Hutchings and reviewed by Véronique Dansereau and Harry Heorton.
Aksenov, Y. and Hibler, W. D.: Failure Propagation Effects in an Anisotropic Sea Ice Dynamics Model, in: IUTAM Symposium on Scaling Laws in Ice Mechanics and Ice Dynamics, edited by: Dempsey, J. P. and Shen, H. H., Solid Mechanics and Its Applications, 363–372, Springer, the Netherlands, 2001. a
Alshibli, K. A. and Sture, S.: Shear Band Formation in Plane Strain Experiments of Sand, J. Geotech. Geoenviron., 126, 495–503, https://doi.org/10.1061/(ASCE)10900241(2000)126:6(495), 2000. a, b
Anderson, E. M.: The dynamics of faulting and dyke formation with applications to Britain, Oliver and Boyd, 1942. a
Arthur, J. R. F., Dunstan, T., AlAni, Q. a. J. L., and Assadi, A.: Plastic deformation and failure in granular media, Géotechnique, 27, 53–74, https://doi.org/10.1680/geot.1977.27.1.53, 1977. a
Badgley, F. I.: Heat balance at the surface of the Arctic Ocean, in: Proceedings of the 29th Annual Western Snow Conference, Western Snow Conference, Spokane, Washington, available at: https://westernsnowconference.org/node/1205 (last access: 3 June 2021), 1961. a
Balendran, B. and NematNasser, S.: Double sliding model for cyclic deformation of granular materials, including dilatancy effects, J. Mech. Phys. Solids, 41, 573–612, https://doi.org/10.1016/00225096(93)90049L, 1993. a, b, c, d
Bolton, M. D.: The strength and dilatancy of sands, ICE Publishing, Géotechnique, 36, 65–78, https://doi.org/10.1680/geot.1986.36.1.65, 1986. 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, b, c, d
Buiter, S. J. H., Babeyko, A. Y., Ellis, S., Gerya, T. V., Kaus, B. J. P., Kellner, A., Schreurs, G., and Yamada, Y.: The numerical sandbox: comparison of model results for a shortening and an extension experiment, Analogue and Numerical Sandbox Models, Geol. Soc. Sp., 253, 29–64, https://doi.org/10.1144/GSL.SP.2006.253.01.02, 2006. a
Campin, J.M., Heimbach, P., Losch, M., Forget, G., Adcroft, A., Dussin, R., et al.: MITgcm/MITgcm: checkpoint67z (Version checkpoint67z), Zenodo, https://doi.org/10.5281/zenodo.4968496, 2021. a
Coon, M., Kwok, R., Levy, G., Pruis, M., Schreyer, H., and Sulsky, D.: Arctic Ice Dynamics Joint Experiment (AIDJEX) assumptions revisited and found inadequate, J. Geophys. Res.Oceans, 112, C11S90, https://doi.org/10.1029/2005JC003393, 2007. a
Coon, M. D., Maykut, A., G., Pritchard, R. S., Rothrock, D. A., and Thorndike, A. S.: Modeling The Pack Ice as an ElasticPlastic Material, AIDJEX Bulletin, 24, 1–106, 1974. a
Coulomb, C. A.: Sur une application des règles de maximis et minimis à quelques problèmes de statique, relatifs à l'architecture, Acad. Sci. Paris Mem. Math. Phys., 7, 343–382, 1776. a, b
Cunningham, G., Kwok, R., and Banfield, J.: Ice lead orientation characteristics in the winter Beaufort Sea, in: Proceedings of IGARSS '94 – 1994 IEEE International Geoscience and Remote Sensing Symposium, 3, 1747–1749, https://doi.org/10.1109/IGARSS.1994.399553, 1994. a, b, c
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
Desrues, J. and Hammad, W.: Shear banding dependency on mean stress level in sand, in: Proc. of the Int. Workshop on Numerical Methods for Localization and Bifurcation of Granular Bodies, GdańskSobieszewo, Poland, 25–30 September, 57–67, 1989. a, b
Dethloff, K., Rex, M., and Shupe, M.: Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC), EGU General Assembly Conference Abstracts, 18, https://ui.adsabs.harvard.edu/#abs/2016EGUGA..18.3064D/abstract, 2016. a
Drucker, D. C. and Prager, W.: Soil mechanics and plastic analysis or limit design, Q. Appl. Math., 10, 157–165, 1952. a
Dumont, D., Gratton, Y., and Arbetter, T. E.: Modeling the Dynamics of the North Water Polynya Ice Bridge, J. Phys. Oceanogr., 39, 1448–1461, https://doi.org/10.1175/2008JPO3965.1, 2009. a
Erlingsson, B.: Twodimensional deformation patterns in sea ice, J. Glaciol., 34, 301–308, 1988. a, b, c, d, e
Girard, L., Bouillon, S., Weiss, J., Amitrano, D., Fichefet, T., and Legat, V.: A new modeling framework for seaice mechanics based on elastobrittle rheology, Ann. Glaciol., 52, 123–132, https://doi.org/10.3189/172756411795931499, 2011. a
Golding, N., Schulson, E. M., and Renshaw, C. E.: Shear faulting and localized heating in ice: The influence of confinement, Acta Mater., 58, 5043–5056, https://doi.org/10.1016/j.actamat.2010.05.040, 2010. a
Han, C. and Drescher, A.: Shear Bands in Biaxial Tests on Dry Coarse Sand, Soil and Foundations, 33, 118–132, https://doi.org/10.3208/sandf1972.33.118, 1993. a, b
Handin, J.: On the Coulomb–Mohr failure criterion, Journal of Geophysical Research (1896–1977), 74, 5343–5348, https://doi.org/10.1029/JB074i022p05343, 1969. a
Heorton, H. D. B. S., Feltham, D. L., and Tsamados, M.: Stress and deformation characteristics of sea ice in a highresolution, anisotropic sea ice model, Philos. T. R. Soc. A, 376, 20170 349, https://doi.org/10.1098/rsta.2017.0349, 2018. a, b
Hibler, W. D.: A viscous sea ice law as a stochastic average of plasticity, J. Geophys. Res., 82, 3932–3938, https://doi.org/10.1029/JC082i027p03932, 1977. a, b, c, d
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, b, c, d, e, f, g
Hibler, W. D. and Schulson, E. M.: On modeling the anisotropic failure and flow of flawed sea ice, J. Geophys. Res.Oceans, 105, 17105–17120, https://doi.org/10.1029/2000JC900045, 2000. a, b, c, d, e, f, g
Hopkins, M. A.: On the ridging of intact lead ice, J. Geophys. Res.Oceans, 99, 16351–16360, https://doi.org/10.1029/94JC00996, 1994. a
Horvat, C. and Tziperman, E.: The evolution of scaling laws in the sea ice floe size distribution, J. Geophys. Res.Oceans, 122, 7630–7650, https://doi.org/10.1002/2016JC012573, 2017. a, b
Hutchings, J. K., Heil, P., and Hibler, W. D.: Modeling Linear Kinematic Features in Sea Ice, Mon. Weather Rev., 133, 3481–3497, https://doi.org/10.1175/MWR3045.1, 2005. a, b, c, d
Hutter, N. and Losch, M.: Featurebased comparison of sea ice deformation in leadpermitting sea ice simulations, The Cryosphere, 14, 93–113, https://doi.org/10.5194/tc14932020, 2020. a, b
Hutter, N., Martin, L., and Dimitris, M.: 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, c, d
Hutter, N., Zampieri, L., and Losch, M.: Leads and ridges in Arctic sea ice from RGPS data and a new tracking algorithm, The Cryosphere, 13, 627–645, https://doi.org/10.5194/tc136272019, 2019. a, b, c
Ip, C. F.: Numerical investigation of different rheologies on seaice dynamics, PhD thesis, Dartmouth College, New Hampshire, United States, 1993. a
Ip, C. F., Hibler, W. D., and Flato, G. M.: On the effect of rheology on seasonal seaice simulations, Ann. Glaciol., 15, 17–25, 1991. a, b, c, d
Itkin, P., Losch, M., and Gerdes, R.: Landfast ice affects the stability of the Arctic halocline: Evidence from a numerical model, J. Geophys. Res.Oceans, 120, 2622–2635, https://doi.org/10.1002/2014JC010353, 2015. a
Kaus, B. J. P.: Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation, Tectonophysics, 484, 36–47, https://doi.org/10.1016/j.tecto.2009.08.042, 2010. a
Koldunov, N. V., Danilov, S., Sidorenko, D., Hutter, N., Losch, M., Goessling, H., Rakowsky, N., Scholz, P., Sein, D., Wang, Q., and Jung, T.: Fast EVP Solutions in a HighResolution Sea Ice Model, J. Adv. Model. Earth Sy., 11, 1269–1284, https://doi.org/10.1029/2018MS001485, 2019. 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., Solid Mechanics and Its Applications, 315–322, Springer, the Netherlands, Dordrecht, https://doi.org/10.1007/9789401597357, 2001. a
König Beatty, C. and Holland, D. M.: Modeling Landfast Sea Ice by Adding Tensile Strength, J. Phys. Oceanogr., 40, 185–198, https://doi.org/10.1175/2009JPO4105.1, 2010. a
Lemieux, J.F. and Tremblay, B.: Numerical convergence of viscous–plastic sea ice models, J. Geophys. Res.Oceans, 114, C05009, https://doi.org/10.1029/2008JC005017, 2009. 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
Losch, M., Fuchs, A., Lemieux, J.F., and Vanselow, A.: A parallel Jacobianfree Newton–Krylov solver for a coupled sea iceocean model, J. Comput. Phys., 257, 901–911, https://doi.org/10.1016/j.jcp.2013.09.026, 2014. a, b
Mancktelow, N. S.: How ductile are ductile shear zones?, GeoScienceWorld, Geology, 34, 345–348, https://doi.org/10.1130/G22260.1, 2006. a
Marko, J. R. and Thomson, R. E.: Rectilinear leads and internal motions in the ice pack of the western Arctic Ocean, J. Geophys. Res., 82, 979–987, https://doi.org/10.1029/JC082i006p00979, 1977. a, b, c
Marshall, J., Adcroft, A., Hill, C., Perelman, L., and Heisey, C.: A finitevolume, incompressible Navier Stokes model for studies of the ocean on parallel computers, J. Geophys. Res.Oceans, 102, 5753–5766, https://doi.org/10.1029/96JC02775, 1997. a
Mohr, O.: Welche Umstände bedingen die Elastizitätsgrenze und den Bruch eines Materials, Zeitschrift des Vereins Deutscher Ingenieure, 46, 1572–1577, 1900. a
Mánica, M. A., Gens, A., Vaunat, J., and Ruiz, D. F.: Nonlocal plasticity modelling of strain localisation in stiff clays, Comput. Geotech., 103, 138–150, https://doi.org/10.1016/j.compgeo.2018.07.008, 2018. a
Nguyen, A. T., Menemenlis, D., and Kwok, R.: Arctic iceocean simulation with optimized model parameters: Approach and assessment, J. Geophys. Res.Oceans, 116, C04 025, https://doi.org/10.1029/2010JC006573, 2011. a
Nguyen, A. T., Kwok, R., and Menemenlis, D.: Source and Pathway of the Western Arctic Upper Halocline in a DataConstrained Coupled Ocean and Sea Ice Model, American Meteorological Society, J. Phys. Oceanogr., 42, 802–823, https://doi.org/10.1175/JPOD11040.1, 2012. a
Overland, J. E., McNutt, S. L., Salo, S., Groves, J., and Li, S.: Arctic sea ice as a granular plastic, J. Geophys. Res., 103, 21845–21868, https://doi.org/10.1029/98JC01263, 1998. a, b
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
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
Rampal, P., Dansereau, V., Olason, 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
Ringeisen, D., Losch, M., Tremblay, L. B., and Hutter, N.: Simulating intersection angles between conjugate faults in sea ice with different viscous–plastic rheologies, The Cryosphere, 13, 1167–1186, https://doi.org/10.5194/tc1311672019, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w
Roach, L. A., Horvat, C., Dean, S. M., and Bitz, C. M.: An Emergent Sea Ice Floe Size Distribution in a Global Coupled OceanSea Ice Model, J. Geophys. Res.Oceans, 123, 4322–4337, https://doi.org/10.1029/2017JC013692, 2018. a
Roscoe, K. H.: The Influence of Strains in Soil Mechanics, Géotechnique, 20, 129–170, https://doi.org/10.1680/geot.1970.20.2.129, 1970. a, b, c, d, e, f
Rothrock, D. A.: The energetics of the plastic deformation of pack ice by ridging, J. Geophys. Res., 80, 4514–4519, https://doi.org/10.1029/JC080i033p04514, 1975. a
Rothrock, D. A. and Thorndike, A. S.: Measuring the sea ice floe size distribution, J. Geophys. Res.Oceans, 89, 6477–6486, https://doi.org/10.1029/JC089iC04p06477, 1984. a
Schall, P. and van Hecke, M.: Shear Bands in Matter with Granularity, Annu. Rev. Fluid Mech., 42, 67–88, https://doi.org/10.1146/annurevfluid121108145544, 2010. a
Schulson, E. M.: Fracture of Ice on Scales Large and Small, in: IUTAM Symposium on Scaling Laws in Ice Mechanics and Ice Dynamics, edited by: Dempsey, J. P. and Shen, H. H., Solid Mechanics and Its Applications, 161–170, Springer, the Netherlands, 2001. a
Schulson, E. M.: Brittle Failure of Ice, GeoScienceWorld, Rev. Mineral. Geochem., 51, 201–252, https://doi.org/10.2138/gsrmg.51.1.201, 2002. a, b, c
Schulson, E. M. and Hibler, W. D.: Fracture of the winter sea ice cover on the Arctic ocean, C. R. Phys., 5, 753–767, https://doi.org/10.1016/j.crhy.2004.06.001, 2004. a
Schulson, E. M., Fortt, A. L., Iliescu, D., and Renshaw, C. E.: Failure envelope of firstyear Arctic sea ice: The role of friction in compressive fracture, John Wiley & Sons, Ltd., J. Geophys. Res.Oceans, 111, C11S25, https://doi.org/10.1029/2005JC003235, 2006a. a
Schulson, E. M., Fortt, A. L., Iliescu, D., and Renshaw, C. E.: On the role of frictional sliding in the compressive fracture of ice and granite: Terminal vs. postterminal failure, Acta Mater., 54, 3923–3932, https://doi.org/10.1016/j.actamat.2006.04.024, 2006b. a
Spreen, G., Kwok, R., Menemenlis, D., and Nguyen, A. T.: Seaice deformation in a coupled ocean–seaice model and in satellite remote sensing data, The Cryosphere, 11, 1553–1573, https://doi.org/10.5194/tc1115532017, 2017. a
Stern, H. L., Rothrock, D. A., and Kwok, R.: Open water production in Arctic sea ice: Satellite measurements and model parameterizations, John Wiley & Sons, Ltd., J. Geophys. Res.Oceans, 100, 20601–20612, https://doi.org/10.1029/95JC02306, 1995. a, b, c, d
Stroeve, J., Barrett, A., Serreze, M., and Schweiger, A.: Using records from submarine, aircraft and satellites to evaluate climate model simulations of Arctic sea ice thickness, The Cryosphere, 8, 1839–1854, https://doi.org/10.5194/tc818392014, 2014. 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, b, c, d, e
Tsamados, M., Feltham, D. L., and Wilchinsky, A. V.: Impact of a new anisotropic rheology on simulations of Arctic sea ice, J. Geophys. Res.Oceans, 118, 91–107, https://doi.org/10.1029/2012JC007990, 2013. a, b
Vardoulakis, I.: Shear band inclination and shear modulus of sand in biaxial tests, Int. J. Numer. Anal. Met., 4, 103–119, https://doi.org/10.1002/nag.1610040202, 1980. a
Vardoulakis, I. and Graf, B.: Calibration of constitutive models for granular materials using data from biaxial experiments, ICE Publishing, Géotechnique, 35, 299–317, https://doi.org/10.1680/geot.1985.35.3.299, 1985. a
Vermeer, P. A.: The orientation of shear bands in biaxial tests, Géotechnique, 40, 223–236, https://doi.org/10.1680/geot.1990.40.2.223, 1990. a, b
Vermeer, P. A. and De Borst, R.: Nonassociated plasticity for soils, concrete and rock, Delft University of Technology, Heron, 29, 1984. a, b
Wang, K.: Pack ice as a twodimensional granular plastic: a new constitutive law, Ann. Glaciol., 44, 317–320, https://doi.org/10.3189/172756406781811358, 2006. a
Wang, K.: Observing the yield curve of compacted pack ice, J. Geophys. Res.Oceans, 112, C05015, https://doi.org/10.1029/2006JC003610, 2007. a
Weiss, J. and Schulson, E. M.: Coulombic faulting from the grain scale to the geophysical scale: lessons from ice, J. Phys. D Appl. Phys., 42, 214 017, https://doi.org/10.1088/00223727/42/21/214017, 2009. a, b
Weiss, J., Schulson, E. M., and Stern, H. L.: Sea ice rheology from insitu, satellite and laboratory observations: Fracture and friction, Earth Planet. Sci. Lett., 255, 1–8, https://doi.org/10.1016/j.epsl.2006.11.033, 2007. a, b
Wilchinsky, A. V., Feltham, D. L., and Hopkins, M. A.: Effect of shear rupture on aggregate scale formation in sea ice, J. Geophys. Res.Oceans, 115, C10 002, https://doi.org/10.1029/2009JC006043, 2010. a
Williams, J., Tremblay, L. B., and Lemieux, J.F.: The effects of plastic waves on the numerical convergence of the viscous–plastic and elastic–viscous–plastic seaice models, J. Comput. Phys., 340, 519–533, https://doi.org/10.1016/j.jcp.2017.03.048, 2017. a
Zhang, J. and Hibler, W. D.: On an efficient numerical method for modeling sea ice dynamics, J. Geophys. Res.Oceans, 102, 8691–8702, https://doi.org/10.1029/96JC03744, 1997. a
Zhang, J. and Rothrock, D. A.: Effect of sea ice rheology in numerical investigations of climate, J. Geophys. Res.Oceans, 110, C08 014, https://doi.org/10.1029/2004JC002599, 2005. a, b, c