Articles | Volume 15, issue 6
The Cryosphere, 15, 2873–2888, 2021
The Cryosphere, 15, 2873–2888, 2021

Research article 24 Jun 2021

Research article | 24 Jun 2021

Non-normal flow rules affect fracture angles in sea ice viscous–plastic rheologies

Non-normal flow rules affect fracture angles in sea ice viscous–plastic rheologies
Damien Ringeisen1,2, L. Bruno Tremblay3, and Martin Losch1 Damien Ringeisen et al.
  • 1Alfred-Wegener-Institut, Helmholtz-Zentrum für Polar- und Meeresforschung, Bremerhaven, Germany
  • 2MARUM – Center for Marine Environmental Sciences, Leobener Str. 8, 28359, Bremen, Germany
  • 3Department of Atmospheric and Oceanic Sciences, McGill University, Montréal, Canada

Correspondence: Damien Ringeisen (


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 uni-axial 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, post-fracture 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 uni-axial loading tests. To this end, an elliptical plastic potential – which defines the post-fracture deformations, or flow rule – is introduced independently of the elliptical yield curve. As a consequence, the post-fracture 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 root-mean-square error below 1.3. The new rheology allows for angles of fracture smaller than 30 in uni-axial 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 uni-axial 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 post-fracture deformation from the shape of the yield curve. The orientation of the post-fracture 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 non-elliptical plastic potential would be required to change the orientation of deformation and to match deformation statistics derived from satellite measurements.

1 Introduction

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 (Badgley1961), 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 Hibler2001). 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; Erlingsson1988; Anderson1942; Schall and van Hecke2010). 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) (Handin1969). The dominant deformation process along LKFs is shear. Sometimes this shear is associated with non-zero 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; Hibler1977), elastic–plastic (EP; Coon et al.1974), elastic–anisotropic–plastic (EAP; Tsamados et al.2013), or Maxwell elasto-brittle (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 post-fracture deformation, called the flow rule. The flow rule is normal to the plastic potential (Drucker and Prager1952). 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 (Hibler1979; Zhang and Rothrock2005) and some with a non-normal flow rule (Ip et al.1991; Tremblay and Mysak1997; Hibler and Schulson2000; Wang2007). 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 non-normal flow rules are often used in other scientific fields to model granular geo-materials (Vermeer and De Borst1984; Mánica et al.2018).

LKFs have been studied in satellite observations (Stern et al.1995; Kwok2001; Schulson and Hibler2004; 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 un-deformed 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 Tremblay2017; Hutter et al.2018), albeit at higher resolution than Maxwell elasto-brittle 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 (Hibler1979) 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 well-studied 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 ϕ (Coulomb1776; Mohr1900):

(1) θ C = π 4 - ϕ 2 .

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 Mysak1997, 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 sub-grid 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

(2) θ R = π 4 - δ 2 .

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; Vardoulakis1980)

(3) θ A = π 4 - 1 4 ( ϕ + δ ) ,

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: θCθθR, with δ<ϕ in sands. If ϕ=δ, then θR=θC=θ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 Nemat-Nasser1993; Vardoulakis and Graf1985; Bolton1986). 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 (Vermeer1990). 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 Schulson2000; Hutchings et al.2005; Wang2006). 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 uni-axial 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 Thomson1977; 15±1.5 in Erlingsson1988; 17 to 18 in Cunningham et al.1994) and a clear peak in the distribution of angles between 20 and 25 (Hutter and Losch2020). In addition, uni-axial 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 Nemat-Nasser1993; Alshibli and Sture2000; Han and Drescher1993; Desrues and Hammad1989). Note that some of these last experiments are tri-axial tests, and that bi-axial tests of 2D granular sea ice might yield different results as sea ice can “escape” in the vertical direction. Bi-axial 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 raft-like 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 (Erlingsson1988; Marko and Thomson1977; 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 Losch2020). 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 Tremblay2017). 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 non-normal flow rule on fracture angles. We use the non-normal flow rule as a means of separating the state of stress (at failure) and the post-fracture deformation. To this end, we study the non-normal 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 non-normal flow rules. Following Ringeisen et al. (2019), we test the new model in simple uni-axial 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 Sea ice model and rheology

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:

(4) ρ h u t = - ρ h f k × u + τ a + τ o - ρ h ϕ s + σ ,

where ρ is the ice density, h is the grid-cell-averaged 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 ϵ˙ 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 ϵ˙ij=(ui)/(xj). The constitutive equations then have the following form:

(5) σ = f ( ϵ ˙ , χ ) .

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 (Hibler1979):

(6) P = P h e - C * ( 1 - A ) ,

where C is a free parameter (typically C=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 χ:

(7) F ( σ i j , χ ) = 0 .

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.

Figure 1Schematic yield curve F (blue) and plastic potential G (red) for a normal (a) and non-normal (b) flow rule. The flow rule (orange) for a given stress on the yield curve is normal to the plastic potential (red) for the same σI. Note that the stress and strain invariant axes are assumed to coincide.


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 (ϵ˙I, ϵ˙II), as shown in orange in Fig. 1, indicates if the grid cell deforms in convergence (ϵ˙I<0) or divergence (ϵ˙I>0) and shear (ϵ˙II). Just as with the yield curve, the plastic potential can be written as

(8) G ( σ i j , χ ) = 0 .

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

(9) G σ i j ( σ i j , χ ) = λ ϵ ˙ i j ,

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 (ϵ˙11, ϵ˙22, ϵ˙12, ϵ˙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 ϵ˙ij=1/2((ui)/(xj)+(uj)/(xi)). 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 non-linear 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

(10) ζ max = 1 2 Δ min P = 2.5 × 10 8 s P

and ηmax=(ζmax)/eG2 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. Δmin=2×10-9s-1 (Hibler1979, 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 non-normal flow rule

We now build a rheology with an elliptical yield curve and a non-normal 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=ktP (König Beatty and Holland2010), the yield curve's ellipse ratio eF, and the plastic potential's ellipse ratio eG – 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

(12) σ II , X = 1 e X P 2 k t - σ I 2 - σ I P ( 1 - k t ) .

Following Hibler (1977, 1979), we derive the constitutive equations σij:

(13) σ i j = 2 η ϵ ˙ i j + ζ - η ϵ ˙ k k δ i j - P ( 1 - k t ) 2 δ i j ,

where the shear and bulk viscosities η and ζ are defined by

(14) ζ = P ( 1 + k t ) 2 Δ and η = ζ e G 2 = P ( 1 + k t ) 2 e G 2 Δ


(15) Δ = ( ϵ ˙ 11 - ϵ ˙ 22 ) 2 + e F 2 e G 4 ( ϵ ˙ 11 - ϵ ˙ 22 ) 2 + 4 ϵ ˙ 12 2 .

Figure 2 shows an example of a yield curve and plastic potential, with the resulting flow rule. For eG>eF, the absolute value of the divergence is smaller and the shear strain rate is larger compared to a normal flow rule (eG=eF) and vice versa for eG<eF.

Figure 2Elliptical yield curve with a non-normal flow rule, a yield curve ellipse aspect ratio eF=2 (blue), and a plastic potential ellipse aspect ratio eG=4 (red). The gray and orange arrows show the normal and non-normal flow rules, respectively.


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 uni-axial compressive test (Ringeisen et al.2019) to materials with a non-associated 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 post-failure 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 3Link between fracture angle and yield curve: (a) arbitrary yield curve F (blue) and plastic potential G (red) in stress invariant space. The plastic potential and yield curve intersect at a stress state p for illustration purposes only. The red arrow is perpendicular to G but non-normal to the yield curve F. The tangent to the plastic potential G at point p has a slope μG=tan (γG) and intersects the σI axis at point q (thin red line). For reference, the normal and tangent to the yield curve F are shown as a thin blue arrow and line. Dashed gray lines show the principal stress axes. (b) Mohr's circle for the fracture state p in (a) (for normal in blue and for non-normal flow rule in red) in the fracture plane of reference (σ, τ) of center σIp and radius σIIp. The thin red line is the tangent to Mohr's circle passing through the point q on the σ axis. By this geometrical construction, sin(δ)=tan(γG)=μG (only valid for |μG|1). δ is called the dilatancy angle. Again for comparison with (a), the blue lines and dots depicts the case of the normal flow rule (i.e., G=F). When considering the plastic potential, the angle of fracture is written from the dilatancy angle δ as θ=π/4-δ/2. By comparison again, δ=ϕ in the case of a normal flow rule and the Roscoe theory (θ=π/4-δ/2) reduces to the Mohr–Coulomb theory (θ=π/4-ϕ/2).


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 uni-axial compression, we (A) determine the stress state on the yield curve for uni-axial compression as a function of the yield curve ellipse ratio eF and (B) compute the slope of the plastic potential at that stress state as a function of the plastic potential ellipse ratio eG. Finally, we combine (2) and (B) to compute the theoretical prediction for the fracture angle as a function of ellipse ratios eG and eF. Figure 3 shows the geometrical construction that links the angle of dilatancy δ to the slope of the plastic potential tan (γG):

(16) sin ( δ ) = tan ( γ G ) = - σ II , G σ I .

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 sin(δ)=ϵ˙I/ϵ˙II, because the ratio of ϵ˙I to ϵ˙II is equal to the slope of the plastic potential -σII,G/σ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

(17) θ R = π 4 - δ 2 .

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 non-normal flow rule in uni-axial compression along the y axis. In this case, σ11=σ12=0; σ22<0; and the principal stresses and stress invariants can be written as


From Eq. (23), the maximum shear stress σII,Fp in the fracture plane in uni-axial compression can be expressed as

(24) σ II , F p ( σ I p ) = - σ I p ,

where p indicates the stress state at the fracture. Figure 4 shows the stress trajectory in principal stress space for uni-axial compression. It also shows how the flow rule changes for the same stress state when using two different elliptical plastic potentials.

Figure 4Trajectory of maximum normal stress (red arrow) in a uni-axial loading test experiment in a material with two different elliptical yield curves (blue) and plastic potentials (dashed orange and dash-dotted teal). The orange and teal arrows show the flow rule normal to the plastic potential of the same color for the same stress state. For eG<eF, the ratio of divergence to shear increases. The opposite is true for eG>eF. A similar figure in principal stress space is presented in Ringeisen et al. (2019).


In the following, we use the normalized stress invariants σI=σI/P and σII=σII/P to simplify the notation. The slope of the yield curve or the plastic potential depends only on eF and eG and not on P. Substituting Eq. (12), σI, and σII in Eq. (24), we obtain

(25) σ II p = - σ I p = 1 e F k t - σ I p ( σ I p + 1 - k t ) ,

and solve the first stress invariant σIp on the fracture plane in uni-axial compression:

(26) σ I p = ( k t - 1 ) - ( 1 - k t ) 2 + 4 k t ( 1 + e F 2 ) 2 ( 1 + e F 2 ) .

The slope of the tangent at σIp to the plastic potential is given by the derivative of Eq. (12):

(27) σ II , G σ I ( σ I p ) = 1 2 e G - 2 σ I p - 1 + k t k t - σ I p ( σ I p + 1 - k t ) .

Substituting Eq. (26) into Eq. (27) yields


or, for zero tensile strength (kt=0),

(29) σ II , G σ I σ I p , k t = 0 = 1 2 e F e G ( 1 - e F 2 ) .

The fracture angle can finally be written as a function of eG and eF from Eq. (18):

(30) θ e , n n ( e F , e G ) = 1 2 arccos 1 2 e F e G ( e F 2 - 1 ) .

As expected, for eF=eG=e, we recover the fracture angle derived in Ringeisen et al. (2019):

(31) θ e , n ( e ) = 1 2 arccos 1 2 1 - 1 e 2 .
3 Experimental setup and numerical scheme

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 Lx=10 km and Ly=25 km. At y=0, we use a closed, solid boundary with a no-slip condition (i.e., u=v=0). At x=0 and Lx, we use Neumann boundary conditions:

(32) .

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=Ly, 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 av is -5×10-4ms-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.

Figure 5Model domain with a solid wall at y=0 (red), Dirichlet boundary conditions where u is 0 at y=0, and prescribed velocities at y=Ly. Open boundaries at x=0,Lx (green) with Neumann boundary conditions. For the conservation of mass, ice thickness, and concentration equations (h, A), Neumann boundary conditions are used on all boundaries. θ is the measured fracture angle with respect to the vertical; the blue line represents an LKF.


The non-linear momentum equation, Eq. (4), is integrated using a Picard solver with 15 000 non-linear (or outer-loop) iterations (Losch et al.2010). For the linearized problem within each non-linear iteration, we use a line successive (over-)relaxation (LSR) method (Zhang and Hibler1997), with a tolerance criterion of |uk-uk-1|max<10-11ms-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 non-linear iterations, but it does so as the non-linear system approaches a converged solution. We chose a very small tolerance and residual norm for the solution of the linear and non-linear problem in order to simulate a clean fracture with a well-defined 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 Tremblay2009). We expect that numerical sea ice models are computationally more challenging with a non-normal flow rule than with a normal flow rule. The non-normality of the flow rule relative to the yield curve introduces more complexity because Drucker's stability postulate is not satisfied (Vermeer and De Borst1984; Balendran and Nemat-Nasser1993). This particular uni-axial 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;, last access: 3 June 2021). We estimated the accuracy as ±1 (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; Ip1993), because it has no influence on the angle of fracture (not shown).

4 Results

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 (kt=0).

Figure 6Diamond-shaped fracture pattern in the shear deformation field ϵ˙II for eF=2.0 and three different values of eG after 5 s of simulation. For the non-normal flow rule (a, c), there are primary and secondary fracture lines, in contrast to the normal flow rule (b) where a single pair of fracture lines are simulated. The fracture angles are 29.92±1.28 for eG=1.4, 34.3±0.25 for eG=2.0, and 40.7±0.94 for eG=4.0. The error corresponds to 2 standard deviations (2σ) of the measured fracture angles.


Figure 6 shows the fracture pattern for the standard yield curve ellipse ratio eF=2.0 and three values of the plastic potential ellipse ratio: eG=1.4, 2.0, and 4.0. The fractures form a diamond shape, which is similar to the shapes observed at large scales (Erlingsson1988) and in laboratory experiments (Schulson2001) 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 (eG=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 diamond-shape floes are simulated. With the non-normal flow rule (eG=1.4 and eG=4) we make three observations:

  1. 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 non-normality of the flow rule (the ratio of divergence to shear strain rate differs from that of the shear to normal stress). For instance, the L2 norm of the residual (R) in the non-linear equations decreases by 4 orders of magnitude for the normal flow rule compared with 2 orders of magnitude for the non-normal flow rule for the same number of non-linear iterations (15 000), specifically to R=8×10-4 for eG=eF=2 and to R=6×10-2 for eG=1.4 and eF=2. Note that a Jacobian-free 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).

  2. The width and activity of the LKFs is also affected by the flow rule. With eG=1.4, the lines are thinner, the shear along the LKFs is smaller, and there is little shear between the fracture lines. With eG=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 eG=1.4, the deformation at the fracture is mainly in divergence, while for eG=4.0, the deformation is mainly in shear and there is more stress transmitted to the ice in between the fracture lines.

  3. The fracture angle changes as the plastic potential changes. The angles are wider with eG=4 than eG=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 eF=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 eG. The choices of yield curve ellipse ratios eF 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 eG changes for each of the four values of eF. There is a clear dependence of the fracture angles on the relative eccentricity of the plastic potential and yield curve. For eG>eF, the shear strain rate increases along the LKFs (see Fig. 6c), and the fracture angles tend toward 45 as eG increases, in agreement with the theory (Eq. 18). For eG<eF, the flow rule implies more divergence (for eF>1 or convergence for eF<1) and less shear along the LKFs (see Fig. 6a), and the fracture angles move away from 45 as eG decreases. More generally, for eF<1, the fracture angle increases with increasing convergence along the LKFs as eG decreases. For eF>1, the fracture angle decreases with increasing divergence as eG decreases. For eF=1 (a circular yield curve), the fracture angles are independent of eG 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 (eG).

The colored, dashed lines in Fig. 7a show the fracture angles θe,nn(eF,eG) predicted by Eq. (30). The coefficient of determination r2 and the root-mean-square error (RMSE) between the simulated angle of fracture and theoretical predictions are 0.97 and 0.37 for eF=0.7, 0.95 and 1.22 for eF=2.0, and 0.97 and 0.47 for eF=4.0. The RMSE is 0.37 for eF=1.0, with r2 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 (eG=eF; Eq. 31).

For completeness, Fig. 7b also show the theoretical predictions for a constant plastic potential ellipse ratio eG for varying yield curve ellipse ratios eF. The fracture angles become smaller as eF increases. Yield curve ellipse ratios smaller than eF=1 do not create fracture angles below 45.

Figure 7(a) Fracture angles as a function of the plastic potential ellipse ratio eG for different yield curve ellipse ratios (eF=0.7, 1.0, 2.0, and 4.0). The markers with ranges are the mean and 2 standard deviations of the fracture angles. The dashed lines show the Roscoe angle (Eq. 30). The arrows mark the Coulomb angles as a function of eF, which are constant with respect to eG. Colors indicate the value of eF for lines and markers. The r2 values between theory and modeled angles for eF=0.7, 2.0, and 4.0 are 0.97, 0.95, and 0.97. (b) Roscoe fracture angle computed from Eq. (30) as a function of eF with a constant eG, for illustration. As eF changes, both the stress state and the flow rule change (see Fig. 4), resulting in more complex behavior. The dotted black line for the normal flow rule (eF=eG) is drawn for reference.


5 Discussion

The idealized experiments using the elliptical yield curve with a non-normal 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 eF<1 does not allow fracture angles smaller than 45 in uni-axial compression. To reduce the fracture angles with yield curve ellipse ratios of eF>1, plastic potential ellipse ratios eG smaller than the yield curve ellipse ratio are required; that is, eG<eF. 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 post-fracture 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 (Roscoe1970).

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:

  1. 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 (Coulomb1776). This behavior is linked to the specific shape of the elliptical yield curve with a maximum shear stress at P/2 and an ascending and a descending part. In principle, we can decrease the fracture angle with increasing shear strength (eF decreasing) by decreasing eG but only if eF>1, but then the flow rule is far from normal, making the numerical convergence difficult.

  2. 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 Drescher1993; Desrues and Hammad1989; Alshibli and Sture2000). 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/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; Schulson2002). 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 (Schulson2002).

  3. 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 eG 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 eF. 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.

  4. The fracture angles in the standard VP models are larger than observed. Using a non-normal flow rule allows us to change the fracture angle in uni-axial 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 non-associated flow rules. Our goal is to provide a reference for the future development of viscous–plastic rheologies with non-normal 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.

Figure 8Alternative yield curves and flow rules: the Mohr–Coulomb yield curve with shear (non-normal) flow rule (blue; Ip et al.1991), the modified Coulombic yield curve with normal (elliptic part) and non-normal (linear part) flow rule (orange; Hibler and Schulson2000), and the teardrop yield curve with a normal flow rule (red; Zhang and Rothrock2005). The elliptical yield curve with eF=2.0 is shown for reference (thin black line). P is the compressive ice strength, and T is the tensile ice strength.


Non-normal flow rules can be combined with the Mohr–Coulomb family of yield curves. For a Mohr–Coulomb yield curve with a double-sliding 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 eGeF; hence δ≃0 and θ=45 (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 Mysak1997) 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; Kaus2010; Mancktelow2006). 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 Schulson2000). 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 non-differentiable 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 non-normal 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 near-normal, 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 non-normal 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 (Rothrock1975; Zhang and Rothrock2005) 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 non-normal 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 normal-flow-rule 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 time-varying 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 time-step, 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 fast-changing 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 Tremblay2017; 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 Schulson2009); the distribution of fracture angles should agree with observations (Marko and Thomson1977; Erlingsson1988; 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 time-varying distribution of the contact normals (Balendran and Nemat-Nasser1993), floe size distributions (Horvat and Tziperman2017; 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 high-spatial-resolution 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 post-fracture deformation in a 2D granular material such as sea ice. This interplay is likely different than for the well-studied case of a solid, homogeneous, 3D block of ice (e.g., Schulson2002). Sea ice floating on the ocean surface can escape vertically when it forms ridges under confined compression (Hopkins1994). 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 laboratory-scale experimental results are available (Weiss et al.2007; Schulson et al.2006b; Weiss and Schulson2009). The sea ice flow size distribution varies in time (summer/winter) and space (marginal ice zone/central Arctic) (Rothrock and Thorndike1984). These variations change the mechanical properties (e.g., distribution of contact normals) and thermodynamic properties (e.g., lateral melt) of sea ice (Horvat and Tziperman2017). Designing more appropriate rheologies for improved high-resolution climate models and more accurate sea ice prediction systems requires consolidated observations of these still-unclear physical processes.

6 Conclusions

The flow rule, which dictates the post-fracture 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 non-normal 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 uni-axial 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 (Roscoe1970; Tremblay and Mysak1997). This new rheology partially solves issues raised in an earlier study (Ringeisen et al.2019). The use of a plastic potential or non-normal 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 non-normal 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 high-resolution 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 high-resolution 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.

Code availability

The sea ice rheology used in this paper is implemented in the sea ice package of the MITgcm, Mid 2021 release (, Campin et al.2021; Marshall et al.1997; Losch et al.2014).

Data availability

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.

Author contributions

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.

Competing interests

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.

Financial support

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 open-access publication were covered by the University of Bremen.

Review statement

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,, 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., Al-Ani, Q. a. J. L., and Assadi, A.: Plastic deformation and failure in granular media, Géotechnique, 27, 53–74,, 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: (last access: 3 June 2021), 1961. a

Balendran, B. and Nemat-Nasser, S.: Double sliding model for cyclic deformation of granular materials, including dilatancy effects, J. Mech. Phys. Solids, 41, 573–612,, 1993. a, b, c, d

Bolton, M. D.: The strength and dilatancy of sands, ICE Publishing, Géotechnique, 36, 65–78,, 1986. a

Bouchat, A. and Tremblay, B.: Using sea-ice deformation fields to constrain the mechanical strength parameters of geophysical sea ice, J. Geophys. Res.-Oceans, 122, 5802–5825,, 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,, 2006. a

Campin, J.-M., Heimbach, P., Losch, M., Forget, G., Adcroft, A., Dussin, R., et al.: MITgcm/MITgcm: checkpoint67z (Version checkpoint67z), Zenodo,, 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,, 2007. a

Coon, M. D., Maykut, A., G., Pritchard, R. S., Rothrock, D. A., and Thorndike, A. S.: Modeling The Pack Ice as an Elastic-Plastic 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,, 1994. a, b, c

Dansereau, V., Weiss, J., Saramito, P., and Lattes, P.: A Maxwell elasto-brittle rheology for sea ice modelling, The Cryosphere, 10, 1339–1359,, 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ńsk-Sobieszewo, 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,, 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,, 2009. a

Erlingsson, B.: Two-dimensional 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 sea-ice mechanics based on elasto-brittle rheology, Ann. Glaciol., 52, 123–132,, 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,, 2010. a

Han, C. and Drescher, A.: Shear Bands in Biaxial Tests on Dry Coarse Sand, Soil and Foundations, 33, 118–132,, 1993. a, b

Handin, J.: On the Coulomb–Mohr failure criterion, Journal of Geophysical Research (1896–1977), 74, 5343–5348,, 1969. a

Heorton, H. D. B. S., Feltham, D. L., and Tsamados, M.: Stress and deformation characteristics of sea ice in a high-resolution, anisotropic sea ice model, Philos. T. R. Soc. A, 376, 20170 349,, 2018. a, b

Hibler, W. D.: A viscous sea ice law as a stochastic average of plasticity, J. Geophys. Res., 82, 3932–3938,, 1977. a, b, c, d

Hibler, W. D.: A Dynamic Thermodynamic Sea Ice Model, J. Phys. Oceanogr., 9, 815–846,<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,, 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,, 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,, 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,, 2005. a, b, c, d

Hutter, N. and Losch, M.: Feature-based comparison of sea ice deformation in lead-permitting sea ice simulations, The Cryosphere, 14, 93–113,, 2020. a, b

Hutter, N., Martin, L., and Dimitris, M.: Scaling Properties of Arctic Sea Ice Deformation in a High-Resolution Viscous-Plastic Sea Ice Model and in Satellite Observations, J. Geophys. Res.-Oceans, 123, 672–687,, 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,, 2019. a, b, c

Ip, C. F.: Numerical investigation of different rheologies on sea-ice 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 sea-ice 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,, 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,, 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 High-Resolution Sea Ice Model, J. Adv. Model. Earth Sy., 11, 1269–1284,, 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,, 2001. a

König Beatty, C. and Holland, D. M.: Modeling Landfast Sea Ice by Adding Tensile Strength, J. Phys. Oceanogr., 40, 185–198,, 2010. a

Lemieux, J.-F. and Tremblay, B.: Numerical convergence of viscous–plastic sea ice models, J. Geophys. Res.-Oceans, 114, C05009,, 2009. a

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

Losch, M., Fuchs, A., Lemieux, J.-F., and Vanselow, A.: A parallel Jacobian-free Newton–Krylov solver for a coupled sea ice-ocean model, J. Comput. Phys., 257, 901–911,, 2014. a, b

Mancktelow, N. S.: How ductile are ductile shear zones?, GeoScienceWorld, Geology, 34, 345–348,, 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,, 1977. a, b, c

Marshall, J., Adcroft, A., Hill, C., Perelman, L., and Heisey, C.: A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers, J. Geophys. Res.-Oceans, 102, 5753–5766,, 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,, 2018. a

Nguyen, A. T., Menemenlis, D., and Kwok, R.: Arctic ice-ocean simulation with optimized model parameters: Approach and assessment, J. Geophys. Res.-Oceans, 116, C04 025,, 2011. a

Nguyen, A. T., Kwok, R., and Menemenlis, D.: Source and Pathway of the Western Arctic Upper Halocline in a Data-Constrained Coupled Ocean and Sea Ice Model, American Meteorological Society, J. Phys. Oceanogr., 42, 802–823,, 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,, 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 elasto-brittle rheology, The Cryosphere, 14, 2137–2157,, 2020. a

Rampal, P., Bouillon, S., Ólason, E., and Morlighem, M.: neXtSIM: a new Lagrangian sea ice model, The Cryosphere, 10, 1055–1073,, 2016. a

Rampal, P., Dansereau, V., Olason, E., Bouillon, S., Williams, T., Korosov, A., and Samaké, A.: On the multi-fractal scaling properties of sea ice deformation, The Cryosphere, 13, 2457–2474,, 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,, 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 Ocean-Sea Ice Model, J. Geophys. Res.-Oceans, 123, 4322–4337,, 2018. a

Roscoe, K. H.: The Influence of Strains in Soil Mechanics, Géotechnique, 20, 129–170,, 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,, 1975. a

Rothrock, D. A. and Thorndike, A. S.: Measuring the sea ice floe size distribution, J. Geophys. Res.-Oceans, 89, 6477–6486,, 1984. a

Schall, P. and van Hecke, M.: Shear Bands in Matter with Granularity, Annu. Rev. Fluid Mech., 42, 67–88,, 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,, 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,, 2004. a

Schulson, E. M., Fortt, A. L., Iliescu, D., and Renshaw, C. E.: Failure envelope of first-year Arctic sea ice: The role of friction in compressive fracture, John Wiley & Sons, Ltd., J. Geophys. Res.-Oceans, 111, C11S25,, 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. post-terminal failure, Acta Mater., 54, 3923–3932,, 2006b. a

Spreen, G., Kwok, R., Menemenlis, D., and Nguyen, A. T.: Sea-ice deformation in a coupled ocean–sea-ice model and in satellite remote sensing data, The Cryosphere, 11, 1553–1573,, 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,, 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,, 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,<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,, 2013. a, b

Vardoulakis, I.: Shear band inclination and shear modulus of sand in biaxial tests, Int. J. Numer. Anal. Met., 4, 103–119,, 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,, 1985. a

Vermeer, P. A.: The orientation of shear bands in biaxial tests, Géotechnique, 40, 223–236,, 1990. a, b

Vermeer, P. A. and De Borst, R.: Non-associated plasticity for soils, concrete and rock, Delft University of Technology, Heron, 29, 1984.  a, b

Wang, K.: Pack ice as a two-dimensional granular plastic: a new constitutive law, Ann. Glaciol., 44, 317–320,, 2006. a

Wang, K.: Observing the yield curve of compacted pack ice, J. Geophys. Res.-Oceans, 112, C05015,, 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,, 2009. a, b

Weiss, J., Schulson, E. M., and Stern, H. L.: Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction, Earth Planet. Sci. Lett., 255, 1–8,, 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,, 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 sea-ice models, J. Comput. Phys., 340, 519–533,, 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,, 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,, 2005. a, b, c

Short summary
Deformations in the Arctic sea ice cover take the shape of narrow lines. High-resolution sea ice models recreate these deformation lines. Recent studies have shown that the most widely used sea ice model creates fracture lines with intersection angles larger than those observed and cannot create smaller angles. In our work, we change the way sea ice deforms post-fracture. This change allows us to understand the link between the sea ice model and intersection angles and create more acute angles.