the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Simulating intersection angles between conjugate faults in sea ice with different viscous–plastic rheologies

### Martin Losch

### L. Bruno Tremblay

### Nils Hutter

Recent high-resolution pan-Arctic sea ice simulations show
fracture patterns (linear kinematic features or LKFs) that are typical of
granular materials but with wider fracture angles than those observed in
high-resolution satellite images. Motivated by this, ice fracture is
investigated in a simple uni-axial loading test using two different
viscous–plastic (VP) rheologies: one with an elliptical yield curve and a
normal flow rule and one with a Coulombic yield curve and a normal flow rule
that applies only to the elliptical cap. With the standard VP rheology, it is
not possible to simulate fracture angles smaller than 30^{∘}. Further,
the standard VP model is not consistent with the behavior of granular
material such as sea ice because (1) the fracture angle increases with ice
shear strength; (2) the divergence along the fracture lines (or LKFs) is
uniquely defined by the shear strength of the material with divergence for
high shear strength and convergent with low shear strength; (3) the angle of
fracture depends on the confining pressure with more convergence as the
confining pressure increases. This behavior of the VP model is connected to
the convexity of the yield curve together with use of a normal flow rule. In
the Coulombic model, the angle of fracture is smaller ($\mathit{\theta}=\mathrm{23}{}^{\circ}$)
and grossly consistent with observations. The solution, however, is unstable
when the compressive stress is too large because of non-differentiable
corners between the straight limbs of the Coulombic yield curve and the
elliptical cap. The results suggest that, although at first sight the large-scale patterns of LKFs simulated with a VP sea ice model appear to be
realistic, the elliptical yield curve with a normal flow rule is not
consistent with the notion of sea ice as a pressure-sensitive and dilatant
granular material.

- Article
(6438 KB) - Full-text XML
- BibTeX
- EndNote

Sea ice is a granular material, that is, a material that is composed of ice floes of different sizes and shapes (Rothrock and Thorndike, 1984; Overland et al., 1998). In most large-scale models, sea ice is treated as a viscous–plastic continuum. It deforms plastically when the internal stress becomes critical in compression, shear, or tension; it deforms as a very viscous (creep) flow when the internal stress is relatively small (e.g., Hibler, 1979; Zhang and Hibler, 1997; Hunke and Dukowicz, 1997). The corresponding highly nonlinear sea ice momentum equations can be solved with modern numerical solvers to reproduce, in a qualitative way, observed linear patterns of sea ice deformation within reasonable computing time (Hutchings et al., 2004; Lemieux et al., 2010; Losch et al., 2010; Hutter et al., 2018). These linear kinematic features (LKFs) are places of large shear and divergence (Kwok, 2001). Leads that open along LKFs are responsible for an emergent anisotropy of such models, affecting the subsequent dynamics, mass balance, and the heat and matter exchanges between the ocean, ice, and atmosphere. It is therefore important to investigate whether sea ice fracture is represented accurately in continuum sea ice models.

The sea ice dynamics are complicated because of sharp spatial changes in material properties associated with discontinuities (e.g., along sea ice leads or ridges) and heterogeneity (spatially varying ice thickness and concentration). The sea ice momentum equations are difficult to solve numerically because of the nonlinear sea ice rheology. Since the first sea ice dynamics model, the elastic–plastic sea ice model based on data collected during the Arctic Ice Dynamics Joint Experiment (AIDJEX; Coon et al., 1974), several approaches to modeling sea ice have been developed. Sea ice has been modeled as an incompressible fluid (Rothrock, 1975), a viscous–plastic (VP) material (Hibler, 1979), an elastic–viscous–plastic (EVP) material (Hunke, 2001), a granular material (Tremblay and Mysak, 1997), an elastic anisotropic plastic (EAP) medium (Wilchinsky and Feltham, 2006), an elastic–decohesive medium (Schreyer et al., 2006), an elasto–brittle (EB) material (Rampal et al., 2016), and a Maxwell(viscous)–elastic–brittle (MEB) material (Dansereau et al., 2016). The actual number of approaches to sea ice modeling in the community, however, is much smaller. For example, 30 out of 33 global climate models in CMIP5 use some form of the standard VP rheology (Stroeve et al., 2014).

In spite of its success, the standard VP rheology is not undisputed. Coon et al. (2007) critically reviewed the assumptions behind current modeling practice since the original model of Coon et al. (1974), namely the zero tensile strength (ice is a highly fractured material) and isotropy assumptions of the sea ice cover and the rheological model. Originally, Coon et al. (1974) assumed sea ice to have cracks in all directions, justifying isotropic ice properties and isotropic rheologies. The use of continuum models such as the standard VP model for high-resolution simulations (grid spacings of 1–10 km) is also debated since the grid size approaches a typical floe size and clearly violates the continuum assumption. For instance, recent high-resolution simulations using the VP model used spatial resolution of approximately 500 m for a regional domain (Wang et al., 2006) and 1 km for a pan-Arctic domain (Hutter et al., 2018). It can be argued that if the mode of deformation of a single floe is similar to that of an aggregate of floes, a given rheology developed for a continuum can still be applicable at spatial resolutions of the order of the floe size (Overland et al., 1998; Feltham, 2008, Appendix C), but the validity of a given flow rule across scales is not clear. At any scale, the assumption of viscous creep for small deformations is not physical, and an elastic model would be appropriate for low stress states. The long viscous timescale, compared to the synoptic timescale of LKFs, of order 30 years (Hibler, 1979), however, allows viscous deformation to be viewed as a small numerical regularization with few implications for the dissipation of mechanical energy from the wind or ocean current (Bouchat and Tremblay, 2014), and the ice model can be considered an ideal plastic material. Tsamados et al. (2013) included anisotropy explicitly in a VP model and show that it improved the representation of ice thickness and ice drift compared to an EVP model. Alternative VP rheologies were never widely used in the community. These include a Coulombic yield curve with a normal flow rule (Hibler and Schulson, 2000), a parabolic lens and a tear drop (Pritchard, 1975), a diamond-shaped yield curve with normal flow rules (Zhang and Rothrock, 2005), a Mohr–Coulomb yield curve with a double-sliding deformation law (Tremblay and Mysak, 1997), or a curved diamond (Wang, 2007).

Previously, fracture lines (LKFs) in the pack ice were explained by brittle fracture (Marko and Thomson, 1977). Similar fracture patterns were also observed, from the centimeter scale in the lab to hundreds of kilometers in satellite observations (Schulson, 2004; Weiss et al., 2007). The scale invariance of the fracture processes at the floe scale has not been shown. This may come from a lack of observations at both high spatial and temporal resolution. Based on satellite observations (e.g., RADARSAT Geophysical Processor System, RPGS, or Advanced Very-High-Resolution Radiometer, AVHRR) and in situ internal ice stress measurements (e.g., from the Surface Heat Budget of the Arctic Ocean, SHEBA, experiment), Weiss et al. (2007) proposed to model winter sea ice as a material that undergoes brittle failure with subsequent inelastic deformation by sliding along LKFs. This idea was formalized with an additional parameterization to simulate damage associated with brittle fracture in an elasto–brittle (EB) and Maxwell–elasto–brittle (MEB) model (Girard et al., 2011; Rampal et al., 2016; Dansereau et al., 2016). We note that subsequent deformation in this model is considered to be elastic deformation (EB) or visco–elastic deformation (MEB) instead of plastic. That is, in the EB and MEB approaches, the material does not weaken when fracture occurs, but rather the Young's modulus is reduced, leading to larger elastic deformation for the same stress. From the scaling behavior of simulated sea ice deformation fields of EVP models (with 12 km grid spacing), it was found that the heterogeneity and the intermittency of deformation in the VP model are not consistent with Radarsat Geophysical Processor System (RGPS) data (Girard et al., 2009). In contrast, VP models were shown to be indeed capable of simulating the probability density functions (PDFs) of sea ice deformations and some of the scaling characteristics over the whole Arctic in agreement with the same observations, either with sufficient resolution (Spreen et al., 2017; Hutter et al., 2018) or with tuned shear and compressive strength parameters (Bouchat and Tremblay, 2017).

High-resolution sea ice models simulate LKF patterns in pack ice, where they
appear as lines of high deformation (Hutchings et al., 2005; Hutter et al., 2018). Previously fractured ice will be weaker and will affect
future sea ice deformation fields. The weakening associated with shear
deformation results from divergence and a reduction in ice concentration
along the LKFs. This mechanism introduces an anisotropy in high-resolution
simulations that is similar to observations with comparable spatial
resolution. Lead characteristics, including intersection angles between LKFs,
were studied a number of times (Lindsay and Rothrock, 1995; Hutchings et al., 2005; Wilchinsky et al., 2010; Bröhan and Kaleschke, 2014; Wang et al., 2016; Hutter et al., 2019). These studies show that VP models produce
LKFs with various confinements, scales, resolutions, and forcings. From
observations with different instruments (Landsat, Seasat/SAR, areal
photographs, AVHRR), typical fracture angles between intersecting LKFs of
$(\mathrm{15}\pm \mathrm{15}){}^{\circ}$ emerge at scales from 1 to 100 km
(Erlingsson, 1988;
Walter and Overland, 1993). Hutter et al. (2019) present an
LKF tracking algorithm and show that fracture angles (half of the
intersection angles) between LKFs in RGPS data follow a broad distribution
that peaks around 20^{∘}, in line with previous assessments
(e.g., Walter and Overland, 1993). Hutter et al. (2019) also show
that the distribution of fracture angles in a VP simulation with 2 km
grid spacing is biased, with a high modal value of 45^{∘} and with too
few small intersection angles between 15 and 25^{∘}. The observed
bias motivates the present investigation of the dependence of fracture angles
in different VP rheologies and model settings, that is, scale, resolution,
boundary conditions, model geometry, and variability in initial ice thickness
field.

The simulation of fractures in sea ice models has been studied in idealized model geometries before. Hibler and Schulson (2000) investigated the effect of embedded flaws – that favor certain angles of fractures – in idealized experiments using a Coulombic yield curve. Hutchings et al. (2005) showed that LKFs can be simulated with an isotropic VP model using an idealized model geometry. The shape of the elliptical yield curve (ratio of shear to compressive strength) in the standard VP model determines if ice arches can form in an idealized channel experiment (Hibler et al., 2006; Dumont et al., 2009). Pritchard (1988) investigated the yield curve's mathematical characteristics and derived angles between the principal stress directions and characteristics directions that depend on the tangent to the yield curve. These results show that stress states exist in plastic materials where no LKFs form and were later used to build a yield curve (Wang, 2007). To build an anisotropic rheology, Wilchinsky et al. (2010) used a discrete element model (DEM) in an idealized model domain and showed clear diamond-shaped fracture patterns. Idealized experiments are also used to investigate new rheologies, for example, the Maxwell–elastic–brittle (MEB) rheology (Dansereau et al., 2016) or the material-point method (MPM) (Sulsky et al., 2007), or to study the theoretical framework explaining the fracture angles (e.g., Dansereau et al., 2017, with a Mohr–Coulomb yield criterion in an MEB model). Recently, Heorton et al. (2018) compared simulated fractures by the EVP and EAP models using an idealized model geometry and wind forcing and showed that the anisotropic model creates sharper deformation features. To the best of our knowledge, the dependency of the fracture angles in sea ice on the shape of the yield curve using high-resolution models has not yet been investigated. This is another motivation of this study.

In this paper, we simulate the creation of a pair of conjugate faults in an ice floe with two different VP rheologies in an idealized experiment at a spatial resolution of 25 m. We explore the influence of various parameters of the rheologies and the model geometry (scale, resolution, confinement, boundary conditions, and heterogeneous initial conditions). The remainder of this paper is structured as follows. Section 2 presents the experimental setup: the VP framework (Sect. 2.1), the definition of the yield curve (Sect. 2.2), and the description of the idealized experiment (Sect. 2.3). Section 3 presents the results: first the control simulation is presented (Sect. 3.1), then we explore the sensitivity of the setup in Sect. 3.2 to scale, resolution and longer run time (Sect. 3.2.1), modified boundary conditions and lateral confinement (Sect. 3.2.2), and to heterogeneity in initial conditions (Sect. 3.2.3). Finally, we consider the effects of two different yield curves with different flow rules in Sect. 3.3: the elliptical (Sect. 3.3.1) and the Coulombic yield curve (Sect. 3.3.2). Discussion and conclusions follow in Sects. 4 and 5.

## 2.1 Viscous–plastic model

We use the Massachusetts Institute of Technology general circulation model
(MITgcm; Marshall et al., 1997) with its sea ice package that
allows for the use of different rheologies (Losch et al., 2010).
All thermodynamic processes have been turned off for our experiments. The
initial sea ice conditions, mean (grid cell averaged) thickness *h* and
fractional sea ice cover *A*, are advected by ice drift velocities with a
third-order flux limiter advection scheme (Hundsdorfer et al., 1995).
Ice drift is computed from the sea ice momentum equations

where *ρ* is the ice density, *h* is the grid cell averaged sea ice
thickness, ** u** is the velocity field,

*f*is the Coriolis parameter,

**is the vertical unit vector,**

*k*

*τ*_{air}is the surface air stress,

*τ*_{ocean}is the ocean drag, ∇

*ϕ*(0) is the gradient of sea surface height, and

**is the vertically integrated internal ice stress tensor. The form of**

*σ***defines the rheology. In the case of the standard VP model described in Hibler (1979), the components of**

*σ***are defined as**

*σ*
where *δ*_{ij} is the Kronecker delta, and summation over equal indices
is implied. *η* and *ζ* are the shear and bulk viscosities,
${\dot{\mathit{\u03f5}}}_{ij}$ is the strain rate tensor defined as

and *P* is the maximum compressive stress defined as a function of the ice
strength parameter *P*^{⋆}, mean sea ice thickness *h*, and the sea ice
concentration *A*:

where *C*^{⋆} is a free parameter.

The stress tensor ** σ** is often expressed in terms of
principal stresses

*σ*

_{1}and

*σ*

_{2}or stress invariants

*σ*

_{I}and

*σ*

_{II}. The principal stresses

*σ*

_{1}and

*σ*

_{2}are the principal components or eigenvalues of the stress tensor on a sea ice element. Eigenvalues always exist because the stress tensor is by definition symmetric. The principal stresses

*σ*

_{1}and

*σ*

_{2}can be expressed as a function of

*σ*

_{ij}as

This change of coordinates can then be represented as a rotation of the
coordinates by *ψ* (Fig. A1). This angle is
(Tremblay and Mysak, 1997)

Any linear combination of the principal stresses consists of stress invariants. One
common set of stress invariants is the mean normal stress (*σ*_{I})
and the maximal shear stress (*σ*_{II}). They can be written as

## 2.2 Yield curve

The VP rheology was originally developed to simulate ice motion on a basin scale (e.g., Arctic Ocean, Southern Ocean) (Hibler, 1979). In this model, stochastic elastic deformation is parameterized as highly viscous (creep) flow (Hibler, 1977). Ice is set in motion by surface air and basal ocean stresses moderated by internal ice stress. When the internal sea ice stress reaches a critical value in compression, tension, or shear, sea ice fails and relatively large plastic deformation takes place. Internal ice stress below these thresholds leads to highly viscous (creep) flow that parameterizes the bulk effect of many small reversible elastic deformation events. The timescale of viscous deformation is so high (≃30 years) that viscous deformation can be seen as regularization for better numerical convergence in the case of small deformation. Plastic deformations are relatively large and irreversible. Viscous deformations are negligibly small; in contrast to elastic deformation, they are also irreversible. The yield criterion is expressed as a 2-D envelope either in principal stress space or stress invariant space with a normal flow rule. The constitutive equations (Eq. 2) are derived assuming that the principal axes of stress coincide with the principal axes of strain. The stress state on the yield curve together with the normal flow rule therefore determines the relative importance of divergence (positive or negative) and shear strain rate at a point. The magnitude of the deformation is such that the stress state remains on the yield cure during plastic deformation.

In this study, we use two different yield curves: an elliptical yield curve
(Hibler, 1979) and a Coulombic yield curve
(Hibler and Schulson, 2000). The elliptical yield curve is used in
conjunction with a normal flow rule, while the Coulombic yield curve uses a
normal flow rule on the elliptical cap and a flow rule normal to the
truncated ellipse for the same first principal stress
(Hibler and Schulson, 2000, Appendix A). For the elliptical yield curve
(Fig. 1, black line), *η* and *ζ* are given by

with the abbreviation

In this abbreviation, the strain rate invariants are the divergence
${\dot{\mathit{\u03f5}}}_{\mathrm{I}}={\dot{\mathit{\u03f5}}}_{\mathrm{11}}+{\dot{\mathit{\u03f5}}}_{\mathrm{22}}$, and the maximum
shear strain rate ${\dot{\mathit{\u03f5}}}_{\mathrm{II}}=\sqrt{({\dot{\mathit{\u03f5}}}_{\mathrm{22}}-{\dot{\mathit{\u03f5}}}_{\mathrm{11}}{)}^{\mathrm{2}}+\mathrm{4}{\dot{\mathit{\u03f5}}}_{\mathrm{12}}^{\mathrm{2}}}$. $e=a/b$ is the
ellipse aspect ratio with the semi-major half-axes *a* and *b* (shown in blue
in Fig. 1). The ellipse aspect ratio *e* defines the shear
strength ${S}^{\star}={P}^{\star}/\mathrm{2}e$ of the material as a fraction of its
compressive strength (Bouchat and Tremblay, 2017). For the Coulombic yield
curve (Fig. 1, red curve), the shear viscosity *η* is capped on
the two straight limbs:

where *μ* is the slope of the Mohr–Coulomb limbs (Fig. 1), and *c*
is the cohesion value (the value of *σ*_{II} for *σ*_{I}=0) defined
relative to the tensile strength by $c=\mathit{\mu}\cdot {T}^{\star}$, where *T*^{⋆} is
defined as a fraction of *P*^{⋆}.

The theoretical angle of fracture *θ* can be calculated from the Mohr's
circle of stress and yield curve written in the local (reference) coordinate
system (Ip et al., 1991; Pritchard, 1988; Hibler and Schulson, 2000). Details are described in Appendix Appendix B. For a Mohr–Coulomb yield
criterion, *θ* follows immediately from the internal angle of friction
or the material shear strength. An instructive analogue is the slope of a
pile of sand on a table. Moist sand has a higher shear strength, and hence the
slope angle can be steeper (i.e., the angle *θ* is smaller).

## 2.3 Idealized experiment

An idealized compressive test is used to investigate the modes of sea ice fracture (Fig. 2). This experiment is standard in engineering (Schulson, 2004; Weiss et al., 2007). The numerical configuration is inspired by Herman (2016) and similar to the one shown in Dansereau et al. (2016). All experiments presented below use the same setup unless specified otherwise. The values of parameters and constants are presented in Table 1.

The model domain is a rectangle of size 10 km×25 km,
except for the experiments presented in Sect. 3.2.1
and 3.2.2. An ice floe of
size 8 km×25 km, surrounded by 1 km of open
water on the eastern and western sides, is compressed with a linearly (in
time) increasing strain rate from the north against a solid southern
boundary. The eastern and western strips of open water avoid interesting
dynamics being confounded by the choice of lateral boundary conditions along
the open boundaries. We use a no-slip condition for the southern boundary,
constraining lateral ice motion. Note that the results presented below are
not sensitive to the choice of boundary condition on the eastern and western
boundaries. Because the simulation time and the ice velocities are small, the
Coriolis force in the momentum equations are neglected. Ocean and sea ice are
initially at rest. The only term left in the momentum equation
(Eq. 1) that is relevant for our experiment is the stress
divergence term, ∇⋅*σ*. The ice floe has a uniform
concentration of 100 % and a thickness of 1 m. The spatial resolution
of the model is 25 m. The angle of fracture is measured with the
angle measuring tool of the GNU Image Manipulation Program (GIMP,
https://www.gimp.org/, Version 2.8.22, last access: 4 April 2019). All angles measured in this study have an error
of approximately 1^{∘}. The finite size of the grid spacing
widens the deformation line, and the fracture spreads over several pixels
because of the obliquity of the fracture. Automatic algorithms for measuring
LKF intersection angles are described in Linow and Dierking (2017) and
Hutter et al. (2019).

We solve the nonlinear sea ice momentum equations with a Picard or fixed point iteration with 1500 nonlinear or outer-loop iterations. Within each nonlinear iteration, the nonlinear coefficients (drag coefficients and viscosities) are updated and a linearized system of equations is solved with a line successive (over-)relaxation (LSR) (Zhang and Hibler, 1997). The linear iteration is stopped when the maximum norm of the updates is less than ${\mathit{\u03f5}}_{\mathrm{LSR}}={\mathrm{10}}^{-\mathrm{11}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, but we also limit the number iterations to 1500. Typically, 1500 nonlinear iterations are required to reach a state close enough to the converged solution. Note that this criterion is much stricter than that proposed by Lemieux and Tremblay (2009) – this is so because of slow convergence due to the highly nonlinear rheology term and the high spatial resolution.

On the open eastern and western boundaries, we use von Neumann boundary conditions for velocity, thickness and concentration, and ice can escape the domain without any restrictions:

where E and W denote the eastern and western boundaries, respectively. Strain is applied to the ice at the northern boundary by prescribing a velocity that increases linearly with time:

where *a*_{v} is the prescribed acceleration, and N denotes the northern boundary.

We use simple uni-axial loading experiments to investigate the creation of pair of conjugate faults and their intersection angle. After presenting the results of simulations with the default parameters (Sect. 3.1), we explore the effects of experimental choices: confining pressure, choice of boundary conditions (i.e., von Neumann versus Dirichlet), domain size, and spatial resolution and inhomogeneities (i.e., localized weakness) in the initial thickness and concentration field (Sect. 3.2). Finally, we study the behavior of two viscous–plastic rheologies with different yield curves and compare these dependencies to what we can infer from smaller- and larger-scale measurements from laboratory experiment and RGPS observations (Sect. 3.3).

## 3.1 Uni-axial compressive test – default parameters

With default parameters (Table 1), a diamond-shaped fracture appears in the shear strain rate and divergence fields after a few seconds of integration (Fig. 3). After one time step (or 0.1 s), the stress states already lie on the yield curve, and the fracture is readily seen in the deformation fields (divergence and shear). We iterate for a total of 20 s in order for the signal to be apparent in the thickness and concentration fields. We do this to more clearly show the link between position of the stress states on the yield curve and the resulting deformation defined by the normal flow rule in the standard VP rheology of Hibler (1979). The shear deformation (${\dot{\mathit{\u03f5}}}_{\mathrm{II}}$) shows where the ice slides in friction and deforms plastically. From Fig. 3, the simulated intersection angle is $\mathit{\theta}=(\mathrm{34}\pm \mathrm{1}){}^{\circ}$.

After a few time steps, the ice thickness decreases, particularly along the
LKFs (Fig. 3c) where divergence is maximal. Note that the loading
axis in our simple 1-D experiment is also the second principal axis, and
consequently the stress states are migrating along the *σ*_{2} axis as the
strain rate at the northern boundary increases. Fracture occurs after plastic
failure when the stress state reaches the yield curve and the ice starts to
move in divergence. This occurs in the half of the ellipse closer to the
origin (for *e*>1) where the normal to the flow rule points in the
direction of positive divergence (or first strain rate invariant) (see
Fig. 4). This explains the simulated divergent flow field
and lower ice thickness particularly along LKFs.

## 3.2 Sensitivity experiments

In this section, we test the sensitivity of the standard VP model simulation (Sect. 3.1) to the choice of resolution, scale, and run time (Sect. 3.2.1), boundary conditions and confinement pressure (Sect. 3.2.2), and heterogeneity in the initial sea ice mass field (Sect. 3.2.3).

### 3.2.1 Domain size, spatial resolution, and length of integration

The angle of intersection between a pair of conjugate faults does not change with domain size and spatial resolution (Fig. 5). This is expected because non-dimensionalizing the divergence of the internal ice stress term (the only term that remains in this simple uni-axial test experiment) by setting ${u}^{\prime}=u/U$, ${x}^{\prime}=x/L$ gives the same equations in non-dimensional form, irrespective of the initial ice thickness or spatial resolution. Consequently, the control and sensitivity experiments are scale-independent, and the behavior of the standard VP model can be readily compared with results from RGPS, AVHRR, or laboratory experiments.

Continuing the integration to 2700 s (45 min), compared to 20 s in the reference simulation, leads to the creation of smaller diamond-shaped ice floes due to secondary and tertiary fracture lines (Fig. 6). The openings are visible in the thickness and concentration fields with thinner, less concentrated ice in the lead. In this longer experiment, the sea ice also ridges, for instance, at the center of the domain, where the apex of the diamonds fails in compression. There is also some thicker ice at the northern boundary induced by the specified strain rate at the northern boundary. The fracture pattern and presence of secondary and tertiary fracture lines are in line with results from laboratory experiments (Schulson, 2004) and with AVHRR and RGPS observations.

In the following, we always show results after 5 s of integration because our main focus is on the initial fracture of the ice, that is, the instant when the ice breaks for the first time under compression.

### 3.2.2 Boundary conditions and geometry

The dynamics responsible for the ice fracture and location of the fracture (presented above) take place far away from the eastern and western boundaries and therefore do not depend on the choice of the corresponding boundary conditions. We now investigate the sensitivity of the results to the choice of boundary condition at the southern boundary. To this end, we force the fracture line to intersect the southern boundary by reducing the domain size to 10 km×10 km with an ice floe of 8 km×10 km in the interior. In this case, the fracture develops from corner to corner, and the angle is solely determined by the geometry of the ice floe, that is, $\mathit{\theta}=\mathrm{arctan}({l}_{x}/{l}_{y})$ (Fig. 7b). With a free-slip boundary condition at the southern boundary, the fracture angle is similar to the one from the control simulation (Fig. 7a). That is, the no-slip condition concentrates the stress to the corner of the ice floe touching the boundary and predetermines the fracture location. A free-slip boundary condition is therefore considered more physical in such idealized experiments where fractures lines can extend from one boundary to another. This result can have implications for simulation of LKFs in the Arctic that would extend from one boundary to another, for instance in the Beaufort Sea.

No-slip or free-slip boundary conditions have little impact on the fracture angle in the larger domain used in the control run simulation because the LKFs always only touch one boundary and end in open water (results not shown). With the free-slip boundary conditions, the stresses and strains are only different south of the diamond fracture pattern because ice can move along the southern boundary, and the second fracture cannot form.

We now explore the effect of confining pressure on the eastern and western
boundaries on the angle of fracture when using a (convex) elliptical yield
curve with a normal flow rule. To do so, we replace the open boundaries to
the east and the west with solid walls and the open water gaps with ice of
thicknesses *h*_{c}. Note that the ice strength is linearly related to the ice
thickness (Eq. 4). Therefore the normal stress at the edge of
the floe is completely defined by the thickness of the surrounding ice.

With an increasing lateral confinement pressure (i.e., an increasing ice
thickness *h*_{c} next to the main floe), all stress states are moved to higher
compressive stresses (blue curve in Fig. 4), and the
fracture angle increases (Fig. 8). In this case, the stress
states are again migrating in a direction parallel to the *σ*_{2} axis but
with a non-zero *σ*_{1} value. The stress states of the ice along the
fracture are therefore located in a region of higher compressive stresses on
the yield curve where the divergence is reduced or even changes sign. With
increasing confinement, the stress states of the ice floe move to more
negative values of *σ*_{1} along a line of constant *σ*_{2} (blue
line in Fig. 4) with deformation moving towards more
convergent states. Between *h*_{c}=0.2 and *h*_{c}=0.3, the regime changes
from lead opening to ridging, as the fracture angle increases to values above
45^{∘}. This is inconsistent with the behavior of a granular
material where the angle of fracture is independent of confining pressure in
uni-axial loading laboratory experiment.

### 3.2.3 Effects of the heterogeneity

So far, all initial conditions have been homogeneous in thickness and concentration within the ice floe. In practice, sea ice (in a numerical model but also in reality) is not homogeneous. A local weakness in the initial ice field is likely the starting point of a crack within the ice field (e.g., Herman, 2016, her Fig. 5c). Local failures raise the stress level in adjacent grid cells, and a crack can propagate. Note that the crack propagation in an “ideal” plastic model such as the VP model is instantaneous, and this propagation is not seen between time steps. As a consequence, lines of failure will likely develop between local weaknesses. The location of weaknesses in the ice field together with the ice rheology (yield curve and flow rule) both determine the fracture angles (Hibler and Schulson, 2000; Aksenov and Hibler, 2001). The influence of previous leads on subsequent lead creation have been studied with a discrete element model (Wilchinsky et al., 2011) and has been used to constrain new anisotropic rheologies that include the effects of embedded anisotropic leads (Wilchinsky and Feltham, 2011, 2012).

To illustrate this behavior, we start new simulations from an initial ice
field with two areas of zero ice thickness and zero ice concentration, hence
weaker ice (Fig. 9a). After 5 s these simulations
yield fracture patterns that are dramatically different from those of the
control run simulation (Sect. 3.1): the
fracture lines now start and terminate at the locations of the weak ice
areas. Still, changing the shear strength of the ice (by changing *e*)
changes the fracture pattern (Fig. 9b and c). With *e*=1, the
angles are much wider than with *e*=2, which is consistent with the general
dependence of fracture angles on *e* (see Sect. 3.3.1). Our simulations cannot lead to conclusive
statements about the relative importance of heterogeneity of initial
conditions and yield curve parameters for the fracture pattern, but we can
state that both affect the simulations in a way that requires treating them
separately to avoid confounding effects. Details are deferred to a dedicated
study.

## 3.3 Effects of the yield curve on the fracture angle

### 3.3.1 Elliptical yield curve

Keeping ${P}^{\star}=\mathrm{27.5}$ kN m^{−1} at its default value, the maximal
shear strength ${S}^{\star}={P}^{\star}/\mathrm{2}e$ is varied by changing the ellipse
ratio *e*. Scaling the absolute values of *P*^{⋆} and *S*^{⋆} while
keeping *e* constant does not change the fracturing pattern as the tangent to
the ellipse stays the same (not shown). Changing the ellipse aspect ratio *e*
has a large effect on the fracture angle. The fracture angle decreases
monotonically as the shear strength of the material (or *e*) decreases, from
61^{∘} for *e*=0.7 to 32^{∘} for *e*=2.6. This is clearly
inconsistent with the behavior of a granular material; in the sand castle
analogue this would correspond to a dry sand castle with steeper walls than a
moist sand castle. From the simple schematic of Fig. 4, it
becomes clear that with increasing *e*, the intersection of the *σ*_{2}
axis with the yield curve gradually migrates from the left side of the
ellipse to the right, where the normal to the yield curve points increasingly
towards convergent motion. We present a theoretical explanation for the
sensitivity of the fracture angle to the shear strength of the material (*e*,
for the ellipse) in Appendix Appendix B by
rewriting the elliptical yield curve in local coordinates in the fracture
plane (*σ*,*τ*) instead of principal or stress invariant coordinates.
The fracture angle is then determined from the slope of the tangent to the
yield curve in local coordinates, and this angle follows from the Mohr's
circle (see, for instance, Popov, 1976).

Bouchat and Tremblay (2017) suggest a smaller ellipse aspect ratio (e.g.,
*e*=0.7) to obtain a closer match with RADARSAT-derived distribution of
deformation rates in pan-Arctic simulations at 10 km resolution. From
Figs. 10 and 11, the corresponding fracture angle
is $\mathit{\theta}=(\mathrm{61}\pm \mathrm{1}){}^{\circ}$, that is, much larger than that is derived
from RADARSAT images. *e* also changes the distribution of the stress states
on the yield curve. As the stress state migrates along the principal stress
*σ*_{2} until it reaches the yield curve in our uni-axial compressive
test, the stress states are in the second half of the ellipse for *e*<1 and
the resulting deformation is in convergence (or ridging). The ice thickness
increases due to ridging along the shear lines (Fig. 11). In a
longer simulation with *e*=0.7 (not shown) the ice does not open but only
ridges, with thicker ice building up within the ice floe. This is in strong
contrast to the results with *e*=2.0 presented in
Sect. 3.2.1, where the initial floe breaks up and
separate floes form.

### 3.3.2 Coulombic yield curve

In this section, we replace the elliptical yield curve with a Coulombic yield
curve (Hibler and Schulson, 2000). This yield curve consists of a
Mohr–Coulomb failure envelope – two straight limbs in principal or stress
invariant space with a slope *μ* – capped by an elliptical yield curve
for high compressive stresses. Note that the flow rule applies only to the
elliptical cap in this yield curve. For the two straight limbs, the yield
curve is normal to the truncated ellipse with the first stress invariant
*σ*_{I}. For a Mohr–Coulomb yield curve, the fracture angle depends
directly on the slope of the Mohr–Coulomb limb of the yield curve.
Appendix Appendix A provides a theoretical
explanation of how the angle of fracture depends on the internal angle of
friction.

The slope of the Mohr–Coulomb limbs of the Coulombic yield curve *μ* is
varied between 0.3 and 1.0 (corresponding to an internal angle of friction
*ϕ*=arcsin (*μ*) of 17.5 to 90^{∘}) to study how the fracture
angle depends on the shear strength of the material. In all experiments with
the Coulombic yield curve, we use a tensile strength of 5 *%* of
*P*^{⋆} and an ellipse ratio *e*=1.4, following
Hibler and Schulson (2000). The tensile strength is introduced mainly for
numerical reasons. With zero tensile strength, the state of stress in a
simple uni-axial compressive test with no confinement pressure is tangential
to the yield curve at the origin (failure in tension) and on the two straight
limbs (failure in shear) simultaneously, resulting in a numerical
instability. With tensile stress (or confinement pressure) included, the
state of stress reaches the yield only on the two limbs of the yield curve
(see Fig. 12a).

For the Coulombic yield curve, there are two distinct regimes of failure.
When the *σ*_{2} axis intersects the yield curve on the two straight
limbs, which happens for our configuration for angles of friction $\mathit{\varphi}<\mathrm{45}{}^{\circ}$ (Fig. 12a, left hand side for *μ*=0.7 or
$\mathit{\varphi}=\mathrm{44}{}^{\circ}$), the angle of fracture $\mathit{\theta}=\mathit{\pi}/\mathrm{4}-\mathit{\varphi}/\mathrm{2}$ as per
standard theory (Appendix Appendix A). When the
*σ*_{2} axis intersects the yield curve on the elliptical cap, which
happens for $\mathit{\varphi}>\mathrm{45}{}^{\circ}$ (Fig. 12c, for *μ*=0.95 or
$\mathit{\varphi}=\mathrm{72}{}^{\circ}$), we observed a discontinuity in the fracture angle
associated with the non-differentiable corner in the yield curve. Note that
this corner cannot be removed (by changing the *P*^{⋆} and *e* of the
elliptical cap) as the two straight Mohr–Coulomb limbs are defined as a
truncation of the ellipse. For $\mathit{\varphi}\approx \mathrm{45}{}^{\circ}$ in our
configuration, the numerical solver has difficulties reaching convergence
because of the non-differentiable corner in the yield curve between the
elliptical cap and the two straight limbs (Fig. 12b
for *μ*=0.8 or $\mathit{\varphi}=\mathrm{53}{}^{\circ}$). Finally for very small angles *ϕ*, a
large number of fractures, as opposed to single well defined fracture lines,
appear because of the weakness of the material in shear. This behavior is not
something that is typically observed in a uni-axial compressive test of a
granular material which generally have higher shear resistance. Note that the
value of *ϕ* that is characteristic of the individual regimes depends on
the amount of tensile strength.

Our idealized experiments using the VP rheologies resolve fracture lines as described by Hutchings et al. (2005) and akin to observations (Kwok, 2001). The fracturing of the ice floe creates smaller floes in a manner that appears realistic, for example, compared to Landsat-7 images (Schulson, 2004, Fig. 2). At the high resolution of 25 m the original interpretation of the continuity assumption, namely that each grid cell should represent a distribution of floes (Coon et al., 1974), is no longer valid, but we show that the fracture angle is independent of resolution and scale as expected. Instead, the emerging discontinuities and the polygonal diamond shape of the fracture lines that appear as floes spanning many grid cells are a consequence of the mathematical characteristics of the VP model (Pritchard, 1988). Diamond-shaped floes are observed in the Arctic Ocean (Erlingsson, 1988; Walter and Overland, 1993) and also modeled using a discrete element model (DEM) in an idealized experiment (Wilchinsky et al., 2010). The elastic anisotropic plastic (EAP) rheology assumes predominately diamond-shaped floes in sea ice (Wilchinsky and Feltham, 2006). A sea ice model with EAP creates sharper fractures than a model with the elastic–viscous–plastic (EVP; Hunke and Dukowicz, 1997) rheology (Heorton et al., 2018). The authors concluded that the anisotropic model may improve the fracturing process for sea ice, especially by creating areas of oriented weaknesses, and particularly at coarse resolution where the fracture is not resolved by the grid spacing. In the experiments presented here, the VP rheologies lead to sharp and anisotropic fracture lines without any additional assumptions.

We explored some experimental choices to separate their effects from those of
the rheology parameters. The fracture angles do not depend on the spatial
resolution and domain size as expected in our idealized numerical experiment
setup (Sect. 3.2.1, Fig. 5). The
maximum viscosities in the VP model are very high, and consequently, the VP
model can be considered an ideal plastic material (i.e., a model with an
elastic component that has an infinite elastic wave speed). For this reason,
fracture in a VP model occurs almost instantaneously. Observed timescales of
fracture are on the order of 10 s for 60 m floe diameters
(Dempsey et al., 2012, Fig. 6b), and from typical
elastic wave speeds of 200–2000 m s^{−1}, large cracks of order
1000 km can form in minutes to hours (Marsan et al., 2012).

In our setup, the no-slip boundary condition has little effect on the fracture pattern, but our results suggest that in basin-wide simulations the choice of boundary conditions affects the fracture depending on the geometry and stress direction. The no-slip condition appears to be unphysical. It acts to concentrate the stress on the corners of the floe and forces the fracture to occur at this location. This should motivate a more thorough investigation of the boundary conditions for LKFs that form between one shoreline and another. Similar results were obtained from analytical solutions in idealized geometry for the Mohr–Coulomb yield curve with a double sliding deformation law (Sirven and Tremblay, 2014).

The confining pressure (i.e., thin ice imposed on the side of the domain) changes the distribution of stress within the domain. This results in different deformation patterns (shear and divergence) and different fracture angles because the yield curve is convex and uses a normal flow rule. From this we can conclude that by surrounding our floe with open water, we get the most acute angles from the rheology in this uni-axial compression setup. This is not consistent with the behavior of typical granular material for which an angle of fracture is independent of the confining pressure (Hutter and Rajagopal, 1994). Details of a heterogeneous ice cover also affect the fracture pattern. LKFs link the weaknesses in the ice cover, but the pattern still depends on the preferred fracture angles implied by the model rheology. In summary, we are confident that our choice of parameters allows us to isolate the effects of the rheology and the yield curve on the fracturing process.

In granular material, large shear resistance is linked to contact normals
between floes that oppose the shear motion and lead to dilatation
(Balendran and Nemat-Nasser, 1993). In our experiments, increasing shear strength
in the standard VP model (reducing the ellipse aspect ratio *e*) does not
decrease but increases the fracture angle. This is in contrast to the
behavior of granular material where larger shear strength leads to lower
fracture angles – think of a moist sand castle versus a dry sand castle. In
addition, high shear strength in the VP model with the elliptical yield curve
leads to convergence along the fracture plane, whereas observations (e.g.,
RADARSAT-derived deformation fields) show a range of positive and negative
divergence along LKFs – in accordance with laboratory tests of granular
material that show a variable internal angle of friction that depends on the
distribution of the contact normals between individual floes
(Hutter and Rajagopal, 1994). Inspection of the stress states in the 2-D stress
plane suggests that the intersection of the yield curve with the *σ*_{2}
axis has an important role in the fracture process. This intersection point
appears to determine the fracture angle. In fact, the angle is determined
from the intersection of the Mohr's circle of stress with the yield curve to
give a theoretical relationship between the fracture angle and the ellipse
ratio *e*. With our experiments, we were able to confirm this relationship
empirically.

Arctic-wide simulations improve metrics of sea ice concentration, thickness,
and velocity by decreasing the value of *e* of the standard elliptical yield
curve, that is, by adding shear and bi-axial tensile and compressive strength
(Miller et al., 2005; Ungermann et al., 2017). The representation of
sea ice arches improves with smaller *e* (Dumont et al., 2009), as do
LKF statistics (Bouchat and Tremblay, 2017). Our results, however, show that
this makes the fracture angles larger, which is in stark contrast to what we
expect to be necessary to improve the creation of LKFs in sea ice models.

The fracture angle and the sea ice opening and ridging depending on the
deformation states are consistent with the theory of the yield curve analysis
developed in Pritchard (1988) and the Mohr's circle
framework that we present in
Appendix Appendix B. Interestingly, a change
of ice maximum compressive strength *P*^{⋆} with a constant *e* has no
influence on the LKF creation, although *P*^{⋆} is usually thought of as
the principal parameter of sea ice models in climate simulations
(e.g., Schmidt et al., 2014). The effects of bi-axial tensile
strength *T*^{⋆} on fracture processes require further investigation,
especially given the fact that the assumption of zero tensile strength is
being challenged (Coon et al., 2007). The ice strength parameter
*C*^{⋆} (the parameter governing the change of ice strength depending on
ice concentration; Eq. 4) was not studied here, although it
appears to be an important tuning parameter, and it also helps to improve
basin-wide simulations (Ungermann et al., 2017). The simulations
presented in this study are not realistic and cannot be compared directly to
observations of ice floe fracture. For instance, our idealized ice floe is
homogeneous, while sea ice is known to feature some weaknesses like thermal
cracks or melt ponds.

With the Coulombic yield curve, the simulated fracture angle can be smaller
than for the elliptical yield curve. For *μ*=0.7 ($\mathit{\varphi}=\mathrm{44}{}^{\circ}$)
theory predicts ${\mathit{\theta}}_{\mathrm{MC}}=\mathrm{22.8}{}^{\circ}$
(Appendix Appendix B). The simulated fracture
angle with *μ*=0.7 of $\mathit{\theta}=\mathrm{23.5}{}^{\circ}$ is close to the $\simeq \mathrm{20}{}^{\circ}$ described in Hibler and Schulson (2000).
Erlingsson (1988) developed a different Mohr–Coulomb
theory linking the internal angle of friction and the fracture angle. This complex
theory takes into account the fractal (or self-similar) nature of sea ice. It
gives different results but is inadequate for a single ice floe simulated as
presented here. Based on the results of Pritchard (1988),
Wang (2007) used observed fracture patterns to design a
curved diamond yield curve. But this yield curve also contains a
non-differentiable point, which will be problematic for numerical reasons.
The Coulombic yield curve used here uses a normal flow, and consequently
divergence will always be present along shear lines. In situ measurements,
however, show that the deformations follow a non-normal flow rule
(Weiss et al., 2007), and large-scale observations show both divergence and
convergence (ridges) along LKFs (Stern et al., 1995). There are
alternative flow rules still to be explored, for example, a double-sliding law
with (Ip et al., 1991) or without dilatation included
(Balendran and Nemat-Nasser, 1993; Tremblay and Mysak, 1997).

Motivated by the observation that the intersection angles in a 2 km Arctic-wide simulation of sea ice are generally larger than in the RGPS dataset (Hutter et al., 2019), the fracturing of ice under compression was studied with two VP rheologies in a highly idealized geometry and with very small grid spacing of 25 m. The main conclusions are given in the following.

In our experimental configuration with uni-axial compression, fracture angles
below 30^{∘} are not possible in a VP model with an elliptical
yield curve. Observations suggest much lower values. We find an empirical
relationship between the fracture angle and the ellipse ratio *e* of the
elliptical yield curve that can be fully explained by the convexity of the
yield curve (Appendix Appendix B). In contrast
to expectations, increasing the maximum shear strength in the sea ice model
increases the fracture angle. Along a fracture line, there can be both
divergence and convergence depending on the shear strength of the ice, linked
to the flow rule. The simulated ice opens and creates leads with an ellipse
ratio *e*>1 (shear strength is smaller than compressive strength) and ridges
for *e*<1 (shear strength is larger than compressive strength).

With a modified Coulombic yield curve, the fracture angle can be decreased to
values expected from observations, but the non-differentiable corner points
of this yield curve lead to numerical (convergence) issues and, for some
values of the coefficient of internal friction *μ*, to fracture patterns
that are difficult to interpret. At these corner points, two different slopes
meet and give two non-unique solutions for fracture angles and deformation
directions. We recommend avoiding non-differentiable yield curves (with a
normal flow rule) in viscous–plastic sea ice models.

More generally, the model produces diamond-shaped fracture patterns. Later the ice floe disintegrates and several smaller floes develop. The fracturing process in the ice floe in our configuration is independent of the experiment resolution and scale but sensitive to boundary conditions (no-slip or free-slip). The fracture angle in the VP model is also sensitive to the confining pressure. This is not consistent with the notion of sea ice as a granular material. Unsurprisingly, the yield curve plays an important role in fracturing sea ice in a numerical model as it governs the deformation of the ice as a function of the applied stress.

The idealized experiment of a uni-dimensional compression is useful to explore the effects of the yield curve because all other parameters are controlled. Historically, the discrimination between the different yield curves was not possible because of the scarcity of sea ice drift data. Model comparisons to recent sea ice deformation datasets, such as from RADARSAT, imply that we would need to increase the shear strength with the ellipse in the standard VP rheology to match observations (Bouchat and Tremblay, 2017). We find that this increases the fracture angles, creating a dilemma. Therefore, the high-resolution idealized experiment presented in this work provides a framework to investigate and discriminate different rheologies – a yield curve and a flow rule.

If Arctic-wide sea ice simulations with a resolution of 25 m are not feasible today because of computational cost, we can still imagine small experiments being useful for process modeling on small scales when local and high-resolution observations (e.g., wind, ice velocities) are available. For example, such process modeling studies could be used to constrain the rheology with data from the upcoming MOSAiC campaign (Dethloff et al., 2016) that will provide a full year of sea ice observations in pack ice. Such simulations would also need to take into account the effects of heterogeneous ice cover and wind patterns, with potentially convergent and divergent wind forcing. Most climate models use the standard VP rheology (Stroeve et al., 2014) or one of its variants (e.g., EVP). Results presented here, however, imply that a more physical yield curve with a (possibly non-associative) flow rule is required. Such a yield curve would have to be continuous in all representations, differentiable without corners, have some cohesion, and be consistent with available observations of fracture angles in convergent and divergent flow.

No datasets were used in this article. All simulation data have been obtained with the MITgcm (http://mitgcm.org, last access: 4 April 2019). Model configuration and code modifications are described in detail in the paper. Additionally they are available on GitHub (https://github.com/dringeis/MITgcm/tree/obcs_seaice_cont+mc, last access: 4 April 2019).

Below, we derive a relationship between the fracture angle and the internal
angle of friction for a Mohr–Coulomb yield criterion for completeness. We
consider an arbitrary piece of a 2-D medium (Fig. A1a) that is
subject to stresses in physical stress space *σ*_{ij} ($i=\mathrm{1},\mathrm{2}$).
Computing the change of coordinates as described in Eq. (5), we
can consider the principal stresses (*σ*_{1},*σ*_{2}) applied on the
medium (Fig. A1b). From the force balance, the normal stress
*σ* and the shear stress *τ* on a plane at an angle *θ* from the
principal stress axis can be written as (see Fig. A1b and
Popov, 1976)

where d*A* is the area of the friction plane on which the stresses are
applied (in 2-D it is just a line). The second trigonometric term comes from
the fact that this surface is tilted compared to the direction of stresses
*σ*_{1} and *σ*_{2}. Using THE angle sum and difference identities of
trigonometry, we can write the stresses *σ* and *τ* in terms of the
principal stresses *σ*_{1} and *σ*_{2} as

In terms of the stress invariants *σ*_{I} and *σ*_{II}, this gives

The Mohr–Coulomb failure criterion can be written in the fracture plane stress space (see Fig. A2) as

where *ϕ* is the internal angle of friction, and *c* the cohesion when no
stresses are applied (Verruijt, 2018). Substituting
Eqs. (A1) and (A1) in Eq. (A1), we get

and after multiplying both sides by cos (*ϕ*),

By geometrical construction (see Fig. A2), the MC criterion is satisfied when (see also Verruijt, 2018, Sect. 20.4)

so that Eq. (A3) becomes

from which we get the classical result of material deformation physics:

A yield curve can be defined in the local stress (*σ*_{ij}), principal
stress (*σ*_{1,2}), or stress invariant (*σ*_{I,II}) spaces. The
latter gives the center and radius of the Mohr's circle of stress defining
all equivalent stress states (*σ*,*τ*) for all angles with respect to
a reference coordinate system. This allows the translation of the elliptical
yield curve from the standard principal or stress invariant space to a local
stress coordinate system (*σ*_{ij}). In this sense, we can plot the
yield curve in (*σ*,*τ*) space as the envelope of all Mohr's circles
for each point on the elliptical yield curve defined in stress invariant
coordinates (see Fig. B1 for an illustration with the
elliptical yield curve). In the following, we refer to this envelope of all
Mohr's circles as the reconstructed yield curve. The tangent to this curve
can be expressed as (Fig. B2)

We can then express the fracture angle for stress states on the yield curve envelope by placing Eq. (B1) in Eq. (A5):

This is the same relation presented (Pritchard, 1988) and
used previously (Wang et al., 2006) but is obtained within the
(*σ*,*τ*) stress space.

## B1 Elliptical yield curve

From the previous equations, some implications about the elliptical yield
curve immediately follow. As shown in Fig. 4, in a
uni-directional compressive setup the slope of a tangent to the yield curve
changes with the ellipse ratio. The convexity of the ellipse implies that the
ratio $\frac{\mathit{\tau}}{\mathit{\sigma}}=\mathrm{tan}\left(\mathit{\varphi}\right)$ of shear strength *τ* to
compressive strength *σ* becomes smaller with smaller *e*. If we compute
the slope of the tangent to the elliptical yield curve at the intersection
point between the yield curve and the *σ*_{2} axis, we get

Inserting this relationship into Eq. (B2) gives the angle of
fracture for the uni-axial compressive experiment with an ellipse ratio *e*:

Note that a yield curve in (*σ*_{I,II}) space with a tangent slope above
unity does not have a Mohr's circle that can be tangent to the yield curve in
(*σ*,*τ*) space (orange circle on in Fig. B1). This
implies that no angle of fracture can be derived for these stress states.
This is the case for the elliptical yield curve for low and high compressive
stresses. It is still unclear what happens in the VP model for stress states
on the yield curve that have a tangent with a slope higher than unity
(see also Pritchard, 1988). Note also that for some
(*σ*_{I},*σ*_{II}) states, the ice will actually fail in tension, as
the reconstructed yield curve with a few points in the first and fourth
quadrant.

The shear and bulk viscosities are symmetrical about the center of the ellipse. This implies that they are equal for divergence and convergence. Clearly this is not physical since, for shear deformations where ice floes continue to interact with one another (termed the quasi-static flow regime (Babić et al., 1990), divergent flow counterintuitively should have more ice–ice interactions and higher viscosities than convergent flow – because divergent flow is the result of a higher number of contact normals opposing the shear. When the divergence is large and floes no longer interact, the shear and bulk viscosities are still symmetrical about the center of the ellipse. While this is nonphysical, it does lead to more numerical stability because the extra viscosity or dissipation of energy regularizes the problem. We also note that a yield curve with a tangent that has a slope smaller than 1 (in absolute value) in the first and fourth quadrant (positive first principal stress) is unphysical because it would lead to a diamond-shaped pair of ice fracture, even in a uni-axial tensile test, which is inconsistent with laboratory experiments (Cox and Richter-Menge, 1985; Menge and Jones, 1993). We conclude that adding tensile strength to the elliptical yield curve may not be physical. The behavior of the elliptical yield curve in uni-axial tensile tests will be explored elsewhere.

## B2 Coulombic yield curve

Applying Mohr's circle to the Coulombic yield curve explains why the
non-differentiable corners in the yield curve lead to numerical problems
(Fig. B3). The tangent does not vary smoothly, and the
reconstructed yield curve in the failure plane (*σ*,*τ*) becomes
discontinuous (Fig. B3, red line). As shown in
Sect. 3.3.2, when the stress states fall on only one of
the two parts (ellipse or limb) the conjugate faults form as expected. Using
Eq. (B3), with *μ* as the slope of the Mohr–Coulomb limbs of
the Coulombic yield curve, the fracture angle is given by

which is identical to Eq. (A5).

The Coulombic yield curve with an internal angle of friction of 1 (*μ*=1)
and no cohesion (*c*=0) (also called the truncated ellipse method,
TEM; Hibler and Schulson, 1997, Appendix) only has one possible solution with an
angle of fracture equal to 0^{∘} (i.e., conjugate pairs of fracture are not
possible). Zero cohesion implies that the ice will deform, even for nearly no
stress. This yield curve also appears unphysical to us.

DR designed the experiments, ran the simulations, and interpreted the results with the help of ML and LBT. NH contributed to the discussion on LKFs in simulations and observations. DR prepared the paper with contributions from all co-authors.

The authors declare that they have no conflict of interest.

We would like to thank Jennifer Hutchings and Harry Heorton for their helpful
reviews on this paper, as well as Amélie Bouchat and Mathieu Plante for
useful discussion during this work. This project was 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” (IRTG 1904 ArcTrain). The authors
would like to thank the Isaac Newton Institute for Mathematical Sciences for
support and hospitality during the program “Mathematics of sea ice
phenomena” when work on this paper was undertaken. This work was supported by EPSRC
grant numbers EP/K032208/1 and EP/R014604/1. This work is a contribution to
the Canadian Sea Ice and Snow Evolution (CanSISE) network funded by the
Natural Sciences and Engineering Research Council (NSERC) of Canada, the
Marine Environmental Observation Prediction and Response (MEOPAR) Network, and
the NSERC Discovery Grants Program awarded to L. Bruno Tremblay.

The article processing charges for this open-access

publication were covered by a Research

Centre of the Helmholtz Association.

This paper was edited by Daniel Feltham and reviewed by Harry Heorton and Jennifer Hutchings.

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, UTAM Symposium, Fairbanks, Alaska, USA, 13–16 June 2000, Kluwer Academic Publishers, 2001. a

Babić, M., Shen, H. H., and Shen, H. T.: The stress tensor in granular shear flows of uniform, deformable disks at high solids concentrations, J. Fluid Mech., 219, 81–118, https://doi.org/10.1017/S0022112090002877, 1990. 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, https://doi.org/10.1016/0022-5096(93)90049-L, 1993. a, b

Bouchat, A. and Tremblay, B.: Energy dissipation in viscous-plastic sea-ice models, J. Geophys. Res.-Oceans, 119, 976–994, https://doi.org/10.1002/2013JC009436, 2014. 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, https://doi.org/10.1002/2017JC013020, 2017. a, b, c, d, e

Bröhan, D. and Kaleschke, L.: A Nine-Year Climatology of Arctic Sea Ice Lead Orientation and Frequency from AMSR-E, Remote Sensing, 6, 1451–1475, https://doi.org/10.3390/rs6021451, 2014. 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, b

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, b, c, d

Cox, G. F. N. and Richter-Menge, J. A.: Tensile Strength of Multi-Year Pressure Ridge Sea Ice Samples, J. Energ. Resour.-ASME, 107, 375–380, https://doi.org/10.1115/1.3231204, 1985. a

Dansereau, V., Weiss, J., Saramito, P., and Lattes, P.: A Maxwell elasto-brittle rheology for sea ice modelling, The Cryosphere, 10, 1339–1359, https://doi.org/10.5194/tc-10-1339-2016, 2016. a, b, c, d

Dansereau, V., Démery, V., Berthier, E., Weiss, J., and Ponson, L.: Fault orientation in damage failure under compression, arXiv:1712.08530 [cond-mat, physics:physics], available at: http://arxiv.org/abs/1712.08530 (last access: 4 April 2019), 2017. a

Dempsey, J. P., Xie, Y., Adamson, R. M., and Farmer, D. M.: Fracture of a ridged multi-year Arctic sea ice floe, Cold Reg. Sci. Technol., 76–77, 63–68, https://doi.org/10.1016/j.coldregions.2011.09.012, 2012. a

Dethloff, K., Rex, M., and Shupe, M.: Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC), EGU General Assembly Conference Abstracts, 18, available at: https://ui.adsabs.harvard.edu/#abs/2016EGUGA..18.3064D/abstract (last access: 4 April 2019), 2016. 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, b

Erlingsson, B.: Two-dimensional deformation patterns in sea ice, J. Glaciol., 34, 301–308, 1988. a, b, c

Feltham, D. L.: Sea Ice Rheology, Annu. Rev. Fluid Mech., 40, 91–112, https://doi.org/10.1146/annurev.fluid.40.111406.102151, 2008. a

Girard, L., Weiss, J., Molines, J. M., Barnier, B., and Bouillon, S.: Evaluation of high-resolution sea ice models on the basis of statistical and scaling properties of Arctic sea ice drift and deformation, J. Geophys. Res.-Oceans, 114, C08015, https://doi.org/10.1029/2008JC005182, 2009. a

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, https://doi.org/10.3189/172756411795931499, 2011. 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, 20170349, https://doi.org/10.1098/rsta.2017.0349, 2018. a, b

Herman, A.: Discrete-Element bonded-particle Sea Ice model DESIgn, version 1.3a – model description and implementation, Geosci. Model Dev., 9, 1219–1241, https://doi.org/10.5194/gmd-9-1219-2016, 2016. 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

Hibler, W. D.: A Dynamic Thermodynamic Sea Ice Model, J. Phys. Oceanogr., 9, 815–846, https://doi.org/10.1175/1520-0485(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 sea-ice fracture and flow in numerical investigations of climate, Ann. Glaciol., 25, 26–32, https://doi.org/10.3189/S0260305500190019, 1997. a

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, h, i

Hibler, W. D., Hutchings, J. K., and Ip, C. F.: sea-ice arching and multiple flow States of Arctic pack ice, Ann. Glaciol., 44, 339–344, https://doi.org/10.3189/172756406781811448, 2006. a

Hundsdorfer, W., Koren, B., vanLoon, M., and Verwer, J. G.: A Positive Finite-Difference Advection Scheme, J. Comput. Phys., 117, 35–46, https://doi.org/10.1006/jcph.1995.1042, 1995. a

Hunke, E. C.: Viscous–Plastic Sea Ice Dynamics with the EVP Model: Linearization Issues, J. Comput. Phys., 170, 18–38, https://doi.org/10.1006/jcph.2001.6710, 2001. a

Hunke, E. C. and Dukowicz, J. K.: An Elastic–Viscous–Plastic Model for Sea Ice Dynamics, J. Phys. Oceanogr., 27, 1849–1867, https://doi.org/10.1175/1520-0485(1997)027<1849:AEVPMF>2.0.CO;2, 1997. a, b

Hutchings, J. K., Jasak, H., and Laxon, S. W.: A strength implicit correction scheme for the viscous-plastic sea ice model, Ocean Model., 7, 111–133, https://doi.org/10.1016/S1463-5003(03)00040-4, 2004. a

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, K. and Rajagopal, K. R.: On flows of granular materials, Continuum Mech. Therm., 6, 81–139, https://doi.org/10.1007/BF01140894, 1994. 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, 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/tc-13-627-2019, 2019. a, b, c, d, e

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

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., no. 94 in Solid Mechanics and Its Applications, Springer Netherlands, 315–322, https://doi.org/10.1007/978-94-015-9735-7_26, 2001. a, b

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

Lemieux, J.-F., Tremblay, B., Sedláček, J., Tupper, P., Thomas, S., Huard, D., and Auclair, J.-P.: Improving the numerical convergence of viscous-plastic sea ice models with the Jacobian-free Newton–Krylov method, J. Comput. Phys., 229, 2840–2852, https://doi.org/10.1016/j.jcp.2009.12.011, 2010. a

Lindsay, R. W. and Rothrock, D. A.: Arctic sea ice leads from advanced very high resolution radiometer images, J. Geophys. Res.-Oceans, 100, 4533–4544, https://doi.org/10.1029/94JC02393, 1995. a

Linow, S. and Dierking, W.: Object-Based Detection of Linear Kinematic Features in Sea Ice, Remote Sensing, 9, 493, https://doi.org/10.3390/rs9050493, 2017. 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, https://doi.org/10.1016/j.ocemod.2009.12.008, 2010. a, b

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

Marsan, D., Weiss, J., Larose, E., and Métaxian, J.-P.: Sea-ice thickness measurement based on the dispersion of ice swell, J. Acoust. Soc. Am., 131, 80–91, https://doi.org/10.1121/1.3662051, 2012. a

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, https://doi.org/10.1029/96JC02775, 1997. a

Menge, J. A. R. and Jones, K. F.: The tensile strength of first-year sea ice, J. Glaciol., 39, 609–618, https://doi.org/10.3189/S0022143000016506, 1993. a

Miller, P. A., Laxon, S. W., and Feltham, D. L.: Improving the spatial distribution of modeled Arctic sea ice thickness, Geophys. Res. Lett., 32, L18503, https://doi.org/10.1029/2005GL023622, 2005. 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

Popov, E. P.: Mechanics of Materials, 2nd edn., Prentice Hall, Englewood Cliffs, N.J., USA, 1976. a, b

Pritchard, R. S.: An Elastic-Plastic Constitutive Law for Sea Ice, J. Appl. Mech., 42, 379–384, https://doi.org/10.1115/1.3423585, 1975. a

Pritchard, R. S.: Mathematical characteristics of sea ice dynamics models, J. Geophys. Res.-Oceans, 93, 15609–15618, https://doi.org/10.1029/JC093iC12p15609, 1988. a, b, c, d, e, f, g

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/tc-10-1055-2016, 2016. a, b

Rothrock, D. A.: The steady drift of an incompressible Arctic ice cover, J. Geophys. Res., 80, 387–397, https://doi.org/10.1029/JC080i003p00387, 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

Schmidt, G. A., Kelley, M., Nazarenko, L., Ruedy, R., Russell, G. L., Aleinov, I., Bauer, M., Bauer, S. E., Bhat, M. K., Bleck, R., Canuto, V., Chen, Y.-H., Cheng, Y., Clune, T. L., Genio, A. D., Fainchtein, R. d., Faluvegi, G., Hansen, J. E., Healy, R. J., Kiang, N. Y., Koch, D., Lacis, A. A., LeGrande, A. N., Lerner, J., Lo, K. K., Matthews, E. E., Menon, S., Miller, R. L., Oinas, V., Oloso, A. O., Perlwitz, J. P., Puma, M. J., Putman, W. M., Rind, D., Romanou, A., Sato, M., Shindell, D. T., Sun, S., Syed, R. A., Tausnev, N., Tsigaridis, K., Unger, N., Voulgarakis, A., Yao, M.-S., and Zhang, J.: Configuration and assessment of the GISS ModelE2 contributions to the CMIP5 archive, J. Adv. Model. Earth Syst., 6, 141–184, https://doi.org/10.1002/2013MS000265, 2014. a

Schreyer, H. L., Sulsky, D. L., Munday, L. B., Coon, M. D., and Kwok, R.: Elastic-decohesive constitutive model for sea ice, J. Geophys. Res.-Oceans, 111, C11S26, https://doi.org/10.1029/2005JC003334, 2006. a

Schulson, E. M.: Compressive shear faults within arctic sea ice: Fracture on scales large and small, J. Geophys. Res.-Oceans, 109, C07016, https://doi.org/10.1029/2003JC002108, 2004. a, b, c, d

Sirven, J. and Tremblay, B.: Analytical Study of an Isotropic Viscoplastic Sea Ice Model in Idealized Configurations, J. Phys. Oceanogr., 45, 331–354, https://doi.org/10.1175/JPO-D-13-0109.1, 2014. 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, https://doi.org/10.5194/tc-11-1553-2017, 2017. a

Stern, H. L., Rothrock, D. A., and Kwok, R.: Open water production in Arctic sea ice: Satellite measurements and model parameterizations, J. Geophys. Res.-Oceans, 100, 20601–20612, https://doi.org/10.1029/95JC02306, 1995. a

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/tc-8-1839-2014, 2014. a, b

Sulsky, D., Schreyer, H., Peterson, K., Kwok, R., and Coon, M.: Using the material-point method to model sea ice dynamics, J. Geophys. Res.-Oceans, 112, C02S90, https://doi.org/10.1029/2005JC003329, 2007. 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/1520-0485(1997)027<2342:MSIAAG>2.0.CO;2, 1997. a, b, c, d

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

Ungermann, M., Tremblay, L. B., Martin, T., and Losch, M.: Impact of the Ice Strength Formulation on the Performance of a Sea Ice Thickness Distribution Model in the Arctic, J. Geophys. Res.-Oceans, 122, 2090–2107, https://doi.org/10.1002/2016JC012128, 2017. a, b

Verruijt, A.: An Introduction to Soil Mechanics, Theory and Applications of Transport in Porous Media, Springer International Publishing, available at: http://www.springer.com/gp/book/9783319611846 (last access: 4 April 2019), 2018. a, b

Walter, B. A. and Overland, J. E.: The response of lead patterns in the Beaufort Sea to storm-scale wind forcing, Ann. Glaciol., 17, 219–226, https://doi.org/10.3189/S0260305500012878, 1993. a, b, c

Wang, K.: Observing the yield curve of compacted pack ice, J. Geophys. Res.-Oceans, 112, C05015, https://doi.org/10.1029/2006JC003610, 2007. a, b, c

Wang, K., Leppäranta, M., and Kõuts, T.: A study of sea ice dynamic events in a small bay, Cold Reg. Sci. Technol., 45, 83–94, https://doi.org/10.1016/j.coldregions.2006.02.002, 2006. a, b

Wang, Q., Danilov, S., Jung, T., Kaleschke, L., and Wernecke, A.: Sea ice leads in the Arctic Ocean: Model assessment, interannual variability and trends, Geophys. Res. Lett., 43, 7019–7027, https://doi.org/10.1002/2016GL068696, 2016. a

Weiss, J., Schulson, E. M., and Stern, H. L.: Sea ice rheology from in-situ, satellite and laboratory observations: Fracture and friction, Earth Planet. Sc. Lett., 255, 1–8, https://doi.org/10.1016/j.epsl.2006.11.033, 2007. a, b, c, d

Wilchinsky, A. V. and Feltham, D. L.: Anisotropic model for granulated sea ice dynamics, J. Mech. Phys. Solids, 54, 1147–1185, https://doi.org/10.1016/j.jmps.2005.12.006, 2006. a, b

Wilchinsky, A. V. and Feltham, D. L.: Modeling Coulombic failure of sea ice with leads, J. Geophys. Res.-Oceans, 116, C08040, https://doi.org/10.1029/2011JC007071, 2011. a

Wilchinsky, A. V. and Feltham, D. L.: Rheology of Discrete Failure Regimes of Anisotropic Sea Ice, J. Phys. Oceanogr., 42, 1065–1082, https://doi.org/10.1175/JPO-D-11-0178.1, 2012. a

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, C10002, https://doi.org/10.1029/2009JC006043, 2010. a, b, c

Wilchinsky, A. V., Feltham, D. L., and Hopkins, M. A.: Modelling the reorientation of sea-ice faults as the wind changes direction, Ann. Glaciol., 52, 83–90, https://doi.org/10.3189/172756411795931831, 2011. 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, b

Zhang, J. and Rothrock, D. A.: Effect of sea ice rheology in numerical investigations of climate, J. Geophys. Res.-Oceans, 110, C08014, https://doi.org/10.1029/2004JC002599, 2005. a