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

The standard viscous-plastic (VP) sea ice model with an elliptical yield curve and normal flow rule does not simulate fracture angles below 30◦ in uni-axial compression, in stark contrast with observations of Linear Kinematic Features (LKFs) in the Arctic Ocean. In this paper, we remove the normality constraint in the standard VP model and study its impact on the fracture angle in a simple uni-axial compressive loading test. To this end, we introduce a plastic potential independent of the yield curve that defines the post-fracture deformations or flow rule. The numerical experiments show that the fracture 5 angle strongly depends on the flow rule details. For instance, a plastic potential with an ellipse aspect ratio smaller than that of the standard ellipse gives fracture angles that are as low as 22◦. A newly adapted theory – based on one developed from observations of granular material – predicts numerical simulations of the fracture angles for plastic materials with a normal or non-normal flow rule with a root-mean-square error below 1.3◦. Implementing an elliptical plastic potential in the standard VP sea ice model requires only minor modifications. The modified rheology, however, takes longer to solve numerically for a 10 fixed level of numerical convergence. In conclusion, the use of a plastic potential addresses several issues with the standard VP rheology: the fracture angle can be reduced to values within the range of satellite observations and it can be decoupled from the exact shape of the yield curve. Furthermore, a different plastic potential function will be required to change the post-fracture deformation along the fracture lines (convergence or divergence) and to make the fracture angle independent on the confining pressure (as in observations). 15


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. Linear Kinematic Features (LKFs), narrow lines of deformation observed in the Arctic sea ice cover, emerge in high-resolution simulations (Kwok, 2001;Hutchings et al., 2005). LKFs in the Arctic sea ice cover influence the Earth 20 system in many ways: heat and matter exchange take place primarily over open water (Badgley, 1965). Salt rejection during ice formation in leads creates dense water and influences the thermohaline circulation (Nguyen et al., 2011(Nguyen et al., , 2012Itkin et al., Mysak, 1997, for a definition of the dilantancy angle δ in the context of sea ice modeling.). The Roscoe angle of fracture is A general theory derived from experiments with sand that takes into account both the angle of friction and the angle of dilation combines the Coulomb and Roscoe angles as (Arthur et al., 1977;Vardoulakis, 1980): (3) 60 Tremblay and Mysak (1997) used this general theory to design their sea ice rheology. Vermeer (1990) proposed a theoretical framework based the grain size and showed that the angle of fracture in most experiments lie between the two extremes: θ C ≤ θ ≤ θ R , with δ < φ in granular materials like sand. If φ = δ then θ R = θ C = θ A , and the flow rule is normal to the yield curve. In other words, the principal axes of stress and the principal axes of strain are coaxial. This condition, however, is not generally satisfied for granular materials (Balendran and Nemat-Nasser, 1993). Experiments with sand have shown differences 65 between φ and δ of the order of 30 • (Vardoulakis and Graf, 1985;Bolton, 1986). Note that both mechanisms, friction and dilatancy, are not radically different: a larger dilatancy angle implies a larger grain size, more contact normals, hence more friction. Ringeisen et al. (2019) used the theory of the internal angle of friction with a normal flow rule in their appendix B to link fracture angles to rheology.
The fracture angles with the standard sea ice rheology cannot be smaller than 30 • in uni-axial compression, even by changing 70 the ellipse aspect ratio e (Ringeisen et al., 2019). This minimum angle is in conflict with observations that report fracture angles (half of the intersection angles) generally below 30 • . Observations report fracture angles of 14 • (Marko and Thomson, 1977), 15 ± 1.5 • (Erlingsson, 1988), 17 • to 18 • (Cunningham et al., 1994), 24 • (µ = 0.9, Weiss and Schulson, 2009), and 20 • to 25 • (Hutter and Losch, 2020). Further, uni-axial compression experiments (Ringeisen et al., 2019) showed that: (1) the angle of fracture is a function of the gradient of shear strength to compressive strength set by the ellipse aspect ratio, (2) the ellipse 75 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 comply with the theory and observations of granular media behaviour, namely that shear band orientations and divergent/convergent motion at the slip lines are a function solely of the shear strength of the material and orientation of the contact normals (or dilatation angle). This unphysical behaviour of the standard VP rheology is connected to the shape of the yield curve in conjunction with a normal flow rule. 80 The flow rule has the advantage that it can be observed with remote sensing methods, contrary to stress which need in-situ measurements. The ratio of shear and divergence along the LKFs allows to infer the dilatancy angle. Observations show that most of the deformation takes place in shear, with 98% of deformations showing more shear than divergence or convergence (Stern et al., 1995). The ellipse ratio of the standard model can be modified to fit this distribution (Bouchat and Tremblay, 2017). Also, laboratory experiments with first-year ice showed that flow rules are non-normal during brittle failures (Weiss 85 et al., 2007). Separating the link between the fracture angle and the flow rule from the yield curve is necessary to design rheologies that are consistent with observed sea ice deformation.
This paper focuses on VP rheologies. Different models represent sea ice dynamics with different material properties, for example, Viscous-Plastic (VP, Hibler, 1977), Elastic-Plastic (EP, Coon et al., 1974), or Maxwell-Elasto-Brittle (MEB, Dansereau potential (flow rule) must be defined. The yield curve defines the stress criteria for the transition from small viscous deformations to large plastic deformations. The plastic potential determines the ensuing post-fracture deformation, called the flow rule.
The flow rule is normal to the plastic potential (Drucker and Prager, 1952). The plastic potential can be set independently, or be equal to the yield curve. In the latter case, the flow rule is also perpendicular to the yield curve and is called a normal-flow rule or associated flow rule. Several yield curves have been used in sea ice models, some with a normal flow rule (Hibler, 1979;95 Zhang and Rothrock, 2005) and some with a non-normal flow rule (Ip et al., 1991;Tremblay and Mysak, 1997;Hibler and Schulson, 2000;Wang, 2007). It is important to note that two models with the same material properties sharing the same yield curve but with different flow rules are two different rheologies.
In this paper, we investigate the effects of a non-normal flow rule on fracture angles and present a generalized framework to use with viscous-plastic material with any flow rules (normal or non-normal). To this end, we introduce a plastic potential 100 independent from yield curve. The new model is tested 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 (2.1), the new rheology (2.2), and a general theory linking the fracture angles and a general flow rule (2.3). The sections 3 and 4 describe the idealized experimental setup and the results.
Section 5 discusses these results and their implication on current and future rheologies. Conclusions follow in section 6. We consider sea ice as a 2D viscous-plastic material. The ice velocities are calculated from the sea ice momentum equations: 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, 110 k is the vertical unit vector, τ a is the surface air stress, τ o is the ocean drag, ∇φ(0) 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 link the stress tensor σ to the deformation tensor˙ and the state variables χ (ice thickness, ice strength, ice concentration, etc.). The components of the strain rate tensor are computed from the velocities as˙ ij = ∂ui ∂xj . The constitutive equations then have the form: In an ideal plastic model, the stresses are independent of the strain rates; in the 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. 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.
One of the state variables in the model is the maximal compressive strength P . This variable represents the maximal compressive stress that sea ice can bear in uniform compression before ridging. We use the simple standard relationship (Hibler, 1979): 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 125 cell, and P is the ice strength of 1 m ice at 100% concentration (A = 1). Some other state variables are a function of P ; for instance, the tensile strength T is usually defined as T = k t · P , where the tensile factor k t > 0 (König Beatty and Holland, 2010). Others are not, such as the ellipse aspect ratio (Hibler, 1979) or the internal angle of friction (Ip et al., 1991). The equation for the yield curve in a VP model is written in terms of the state variables.
For two-dimensional sea ice, stress is a rank two tensor; thus, it has four components. The yield curve represents the stress 130 states that are deforming plastically while enclosing the states of stress for slow viscous deformations. We express the yield curve as a function of the stresses σ ij and the state variables χ: The yield curve can be represented in principal stress (σ 1 and σ 2 ) or stress invariants space (σ I and σ II ). Figure 1 shows an arbitrary yield curve in stress invariants space. Although equation (7) determines if the deformation is plastic or viscous, it 135 does not determine how the ice will deform after fracture. In order to obtain a closed system of equations, we define a plastic potential that defines the flow rule.
The plastic potential determines the direction of deformation for stress states on the yield curve. Just as the yield curve, the plastic potential can be written as: The direction of the deformation, called the flow rule, is perpendicular to the plastic potential. This is shown in red on Fig. 1b and mathematically expressed by where λ > 0 is the unknown flow rate. The flow rule is applied for stress states on the yield curve at the same compressive stresses (orange arrows in Fig. 1b). If the plastic potential and the yield curve are the same (G = F ), the flow rule is called an 145 associative or normal flow rule, as the flow rule is also perpendicular to the yield curve (see Fig. 1a).
Using Eq. (7) and Eq. (9), we can write a system of 5 equations (four from Eq. 9 and one from Eq. 7) for 5 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 = 150 σ 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 changes 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 155 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.

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 160 different from the yield curve F . We use a different, but still elliptical plastic potential for simplicity: this choice only requires only minor modifications to a typical VP sea ice model. We define the yield condition F and the plastic potential G as a function of the state variables χ: the ice compression strength P , the ice tensile strength T = k t P (König Beatty and Holland, 2010), the yield curve's ellipse ratio e F , and the plastic potential's ellipse ratio e G by 165 for X = F, G for the yield curve or the plastic potential. Using Eq. (10), we write σ II as a function of σ I as: Figure 2. Elliptical 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.
Following Hibler (1977Hibler ( , 1979, we derive the constitutive equations σ ij : where the shear and bulk viscosities η and ζ are defined by: Figure 2 shows an example of yield curve and plastic potential, with the resulting flow rule. For e G > e F , the absolute value of the divergence is smaller and the shear strain rate is larger compared to a normal flow rule (e G = e F ) and vice versa for 175 e G < e F .

Linking fracture and flow rule
In this section, we generalize the theory linking the rheological model and the fracture angles in 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.
180 Figure 3 illustrates the case of an arbitrary yield curve with an arbitrary plastic potential. The figure shows the geometrical construction that links the angle of dilatancy δ to the slope of the plastic potential tan(γ G ): Figure 3. Link between fracture angle and yield curve: a) Arbitrary yield curve F (blue) and plastic potential G (red) in stress invariants space. The plastic potential and yield curve intersect at a stress state p for illustration purposes only. The orange 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. Mohr's circle of stress for this stress state (blue) has a radius of σ p II an centers on σ p I on the σI axis. Gray dashed lines show the principal stress axes. b) Mohr's circle for the fracture state p in a) in the fracture plane of reference (σ, τ ) of center σ p I and radius σ p II . The thin red line is the tangent to the Mohr's circle that passes through the point q on the σ axis. By this geometrical construction, sin(δ) = tan(γG) = µG (for |µG| ≤ 1).
δ is called the dilatancy angle. Again for comparison, the transparent blue lines show the corresponding construction for a normal flow rule from panel a).
Note that the minus sign above was included in the derivative of the yield curve in Eq. (B1) and (B2) of Ringeisen et al. (2019).
This equation agrees with the definition of Roscoe (1970) sin(δ) =˙ İ II , because the ratio of˙ I to˙ II is equal to the slope of 185 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: Substituting Eq. (15) 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 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. 21, the maximum shear stress σ p II,F in the fracture plane in uni-axial compression can be expressed as

200
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.
In the following, we use the normalized stress invariants σ I = σI P and σ II = σII P to simplify the notation. The slope of yield curve or the plastic potential depends only on e F and e G , but not on P . Substituting Eq. (11), σ I , and σ II in Eq. (22), we obtain, and solve the first stress invariant σ p I on the fracture plane in uni-axial compression .
The slope of the tangent at σ p I to the plastic potential is given by the derivative of Eq. (11): Substituting Eq. (24) into Eq. (25), yields or for zero tensile strength (k t = 0), The fracture angle can finally be written as a function of e G and e F from Eq. (17): θ e,nn (e F , e G ) = 1 2 arccos 1 2e F e G (e 2 F − 1) .
As expected, for e F = e G = e we recover the fracture angle derived in Ringeisen et al. (2019): 3 Experimental setup and numerical scheme On the left and right sides of the domain (x < 1 km and x > 9 km), we have open water between the ice floe and the boundary to ensure that the boundaries have no effect on the simulation. At (y = L y ), we use a Dirichlet boundary condition for ice velocity (v the velocity in y-direction increasing linearly in time simulating an axial loading test) and a Neumann boundary condition for ice thickness and concentration : v(t)| y=Ly = a v · t , u(t)| y=Ly = 0 ; ∂A ∂y y=Ly = ∂h ∂y y=Ly = 0 The non-linear momentum equations are integrated using a Picard solver with 15 000 non-linear (or outer-loop) iterations . For the linearized problem within each iteration, we use a line successive (over-)relaxation (LSR) method (Zhang and Hibler, 1997), with a tolerance criterion of |u k −u k−1 | max < 10 −11 m s −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 235 not reach the tolerance criterion for the first non-linear iterations, but 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 Tremblay, 2009). We expect modeling sea ice with a non-normal flow rule to be more challenging than with a normal flow rule. The non-coaxiality of the 240 deviatoric stress and strain rate introduces more complexity because Drucker's postulate for stability is not respected (Vermeer and De Borst, 1984;Balendran and Nemat-Nasser, 1993). 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 with the Measure Tool from the GNU Image Manipulation Program 245 (GIMP, version 2.8.16, gimp.org). The first 5 seconds of simulations are used to define the sea ice fracture and calculate the fracture angle. The angle of each fracture lines is measured and used to compute the average fracture angle and the standard deviation σ = √ VAR. Note that the fracture angles do not depend on resolution, scale, geometry, or boundary conditions (see Ringeisen et al., 2019, their Sec. 3.2 ). We do not use a replacement pressure scheme (Ip et al., 1991;Ip, 1993), because it has no influence with 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 (k t = 0).  3. The fracture angle changes as the plastic potential changes. The angles are wider with e G = 4 than e G = 1.4. The effect of flow rule orientation on the fracture angles is discussed below.

275
We now present results from four sets of simulations with fixed yield curve ellipse ratios at e F = 0.7, 1.0, 2.0, 4.0. For each of these, we test the sensitivity of the results to changes in the plastic potential ellipse ratio e G . The choice of yield curve ellipse ratios e F are: the standard value of (Hibler, 1979), values suggested by Bouchat and Tremblay (2017) and Dumont et al. (2009), and an extreme value resulting in a very small shear strength and smaller fracture angles. shear along the LKFs, and the fracture angles move away from 45 • as e G decreases. More generally, for e F < 1, the fracture angle increases with increasing convergence along the LKFs as e G decreases. For e F > 1, the fracture angle decreases with 285 increasing divergence as e G becomes smaller. For e F = 1 (a circular yield curve), the fracture angles are independent of e G because the fracture takes place at the peak of the yield curve and the flow rule is not affected by changes of the plastic potential ellipse ratio (e G ). Fig. 7 show the fracture angles θ e,nn (e F , e G ) predicted by Eq. (28). The coefficient of determination R 2 and the root-mean-square error (RMSE) between the simulated angle of fracture and theoretical predictions For completeness, Figure 7b also show the theoretical predictions for a constant plastic potential ellipse ratio e G as the yield 295 curve ellipse ratio e F changes. The modeled angles for e G = 4.0 are shown as an example. The fracture angles become smaller as the e F increase. Yield curve ellipse ratio smaller than e F = 1 do not create fracture angles below 45 • .

The coloured dashed lines in
The idealized experiments using the elliptical yield curve with a non-normal flow rule show that the type of deformation and the fracture angle are intimately linked with the shape of plastic potential. We observe that a yield curve ellipse ratio e F < 1 does 300 not allow for fracture angles smaller than 45 • in uni-axial compression when using an elliptical plastic potential, irrespective of the plastic potential elliptical aspect ratio. To reduce the fracture angles with yield curve ellipse ratios e F > 1, one needs to use plastic potential ellipse ratios e G smaller than the yield curve ellipse ratio, that is, e G < e F . The numerical experiments show that the use of a plastic potential in a viscous-plastic model allows separating the yield criterion from the resulting deformation (flow rule). This allows to decouple the mechanical strength properties of the material (ice) from its post-fracture behaviour.

305
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 type of deformation (convergence or shear) along the fracture line and the fracture angle. The resulting fracture angles are in excellent agreement with the Roscoe angle predictions (Roscoe, 1970).
Understanding the link between rheology and fracture angle is necessary for choosing or designing a rheology that is capable 310 of reproducing the observed angle of intersection between LKFs and consequently the emerging anisotropy. In principle, it may be possible to solve several inconsistencies of the standard elliptical yield curve with a normal flow rule Ringeisen et al.
(discussed in 2019) by using a plastic potential that is independent of the yield curve, namely: 1. In the standard VP model with elliptical yield curve and normal flow rule, adding shear strength increases the fracture angle, contrary to granular matter theory (Coulomb, 1773). This behaviour is linked to the shape of the elliptical yield 315 curve. In principle, we can decrease the fracture angle with increasing shear strength (e F decreasing) by decreasing e G , but only if e F > 1. When doing this, the flow rule becomes very non-normal, making the numerical convergence difficult.
2. The angle of fracture in the standard VP model changes with confining pressure unlike laboratory experiments with granular materials (e.g. sand) where the fracture angle is relatively insensitive to the confining pressure (Alshibli and 320 Sture, 2000). This behaviour cannot be eliminated with an elliptical plastic potential, as the normal stresses along the LKFs increases with confining pressure and the flow rule changes from divergence to convergence. A different plastic potential function would change this behaviour. However, this would make the model implementation and numerical convergence more difficult. However, 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 3rd dimension" after fracture. 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 e G changes the flow rule but does not change the sign of the divergence along the LKFs which is solely determined by the yield curve ellipse ratio e F . With the elliptical plastic potential, convergent motion remains convergent and only the ratio of shear to convergence changes. To change this behaviour, 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 • ; something that is not possible with a normal flow rule (Ringeisen et al., 2019).
Other yield curves have been used in previous studies; for instance, the Mohr-Coulomb, the Coulombic yield curve, or the 335 teardrop yield curves. The use of a plastic potential in conjunction with these yield curves may also prove useful in solving these issues. 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 apply knowledge gained from the simulations presented above, and discuss how alternative yield curves may address deficiencies in the standard VP rheology.
It is possible to include different (non-normal) flow rules with the Mohr-Coulomb family of yield curve. The  yield curve with a pure shear flow rule (Ip et al., 1991) would create fracture angle approximately equal to 45 • , independently of the slope of the yield curve. It corresponds to the case e G e F , hence δ 0 and θ = 45 • , as shown in Fig. 7). This contradicts the Coulomb angle theory, which predicts an angle of fracture that depends solely on the internal angle of friction (Eq. 1). Including an angle of dilatancy with a Mohr-Coulomb yield curve (Tremblay and Mysak, 1997) would allow for different angles of fracture with shear and divergence or convergence along the LKFs depending on δ. Such fracture angle and 345 divergence would be independent of the shear strength and the confining pressure in agreement with Roscoe's angle of fracture.
Such a rheology could potentially solve the all four issues. It is also important to note that the Mohr-Coulomb yield curves do not satisfy the convexity requirements of Drucker's postulate of stability. Mohr-Coulomb yield curves in plastic earth mantle models showed different fracture angles corresponding to the Coulomb angle, Roscoe angle, and the intermediate Arthur angles (Buiter et al., 2006;Kaus, 2010;Mancktelow, 2006). However, such geological models are usually incompressible, and 350 making inferences for the compressible formulation of sea ice models is difficult. The investigation of the fracture angles with Mohr-Coulomb yield curves is left for future work.
The Coulombic yield curve uses the two straight limbs from the Mohr-Coulomb yield curve and an elliptical cap of the standard VP rheology for large compressive stresses (Hibler and Schulson, 2000). In this rheology, the flow rule over the two straight limbs is defined by the elliptical yield curve; that is, the ellipse serves as a plastic potential for the  yield curve. The Coulombic yield curve leads to unrealistic and asymmetrical fracture lines i) when the stress states lie at the non-differentiable intersection between the straight limbs and the elliptical cap (Ringeisen et al., 2019), and ii) when the stress states lie on 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 located 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 ellipse plastic potential is nearly perpendicular to the limbs 360 of the Mohr-Coulomb yield curve (Ringeisen et al., 2019).
The teardrop yield curve with a normal flow rule (Rothrock, 1975;Zhang and Rothrock, 2005) is divergent for a wide range of normal stresses and for all practical purposes consists of a continuously differentiable version of the Coulombic yield curve. This asymmetry between divergent and convergent deformation for different normal stresses decreases the effect of confinement on the fracture angle -issue 2 -and reduces the fracture angle for any confinement pressure -issue 4. This yield 365 curve does not address issue 1, because adding shear strength in the tear-drop yield curve also increases the fracture angle.
The main disadvantage of a non-normal flow rule is the slower numerical convergence. 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 experiments are entirely symmetrical. The fracture lines with a normal flow rule are symmetrical and come in 370 pairs (Ringeisen et al., 2019). The poorer numerical convergence in practice will go unnoticed in high-resolution simulations using realistic geometries, since the number of iterations typically used (O(10)) is much smaller than that required for full convergence. On the upside, at each time-step, a new iteration typically use the solution from the previous timestep as initial conditions. 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. Whether this behaviour (asymmetry and multiple fracture lines) will 375 also be present in realistic simulation using spatially and temporally varying wind forcing remains to be tested.
The following criteria should be considered when building a new rheology. The spatial and temporal scaling of sea-ice deformation should agree with observations (Bouchat and Tremblay, 2017;Hutter et al., 2018); the flow rule should reproduce the correct divergence along LKFs (Stern et al., 1995); the yield curve includes some tensile strength (Coon et al., 2007) and reproduces observed distributions of internal stress when ice deforms (Weiss and Schulson, 2009); the distribution of 380 fracture angles should agree with observations (Marko and Thomson, 1977;Erlingsson, 1988;Cunningham et al., 1994;Hutter et al., 2019); the sea ice mechanical strength properties (yield curve) and deformation (flow rule) should vary in time and space depending on the time-varying distribution of the contact normals and floe size distributions, as per observations and laboratory or numerical experiments (Overland et al., 1998;Hutter et al., 2019;Horvat and Tziperman, 2017;Roach et al., 2018;Balendran and Nemat-Nasser, 1993;Dansereau et al., 2016).

385
Although high spatial resolution observations from satellite are available from optical instruments (e.g., from the Landsat or Sentinel programs), higher temporal resolution of sea ice deformation and flow size distributions are still lacking. So is the combined knowledge of the failure stresses and their associated deformation of sea ice as a 2D granular material. Sea ice floating on the ocean surface can "escape" vertically when ridging under confined compression (Hopkins, 1994). This behaviour differs from laboratory test that use axial symmetry and general knowledge about 3D granular material like sand. Generally, 390 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 laboratory scale experimental results are available Weiss et al. (2007); Schulson et al. (2006); Weiss and Schulson (2009). The sea ice flow size distribution varies in time (summer/winter) and space (marginal ice zone/central Arctic) (Rothrock and Thorndike, 1984). These variations change the mechanical properties (e.g., distribution of contact normals) and thermodynamic properties (e.g., lateral melt) of sea ice (Horvat and Tziperman, 2017).

395
Consolidated observations of these two physical processes are needed to design new rheologies for high-resolution climate modelling.
The flow rule, which dictates the post-fracture deformation, has a fundamental effect on the orientation of fractures lines in a viscous-plastic sea ice model. To test this, we added an elliptical plastic potential (allowing for a non-normal flow rule) to the 400 standard viscous plastic rheology with an elliptical yield curve. We tested this new rheology with numerical experiments in uni-axial compression using the standard viscous-plastic model of Hibler (1979). The modeled fracture angles are in agreement with the Roscoe angle, a theory based on experiments with granular materials that includes an angle of dilatancy (Roscoe, 1970;Tremblay and Mysak, 1997). This new rheology partially solves issues raised in an earlier study (Ringeisen et al., 2019). The use of a plastic potential or non-normal flow rule allows for the simulations of smaller fracture angles between pairs of Linear 405 Kinematic Features, in agreement with satellite observations. The fracture angles, however, still depend on the confinement pressure, and the elliptical plastic potential does not modify the direction of deformation at the fracture lines (convergence or divergence), only the ratio of divergence relative to shear. The momentum equations for a rheological model with a nonnormal flow rule are more difficult to solve numerically, and produce multiple lines of fractures that are asymmetrical (despite symmetry of the problem), in contrast with a model with a normal flow rule. Understanding the effect of the flow rule on the 410 fracture angle is necessary to design VP rheologies for high-resolution sea-ice modeling that both reproduce fracture angles and deformation along the fracture lines, and the behaviour of sea ice as a granular material.
Designing a rheology for high-resolution simulation requires information on 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 unavailable at high temporal and spatial resolution. The observations of the Multidisciplinary drifting Observatory for the Study of Arctic Climate 415 (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 modified version of MITgcm used in this study is available at https://github.com/dringeis/MITgcm/tree/ell_nnfr 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 of BT and ML.