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

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.

1 Introduction

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 Experimental setup

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

$$\begin{array}{}\text{(1)}& \mathit{\rho }h{\displaystyle \frac{\partial \mathit{u}}{\partial t}}=-\mathit{\rho }hf\mathit{k}\times \mathit{u}+{\mathit{\tau }}_{\mathrm{air}}+{\mathit{\tau }}_{\mathrm{ocean}}-\mathit{\rho }h\mathrm{\nabla }\mathit{\varphi }\left(\mathrm{0}\right)+\mathrm{\nabla }\cdot \mathit{\sigma },\end{array}$$

where ρ is the ice density, h is the grid cell averaged sea ice thickness, u is the velocity field, f is the Coriolis parameter, k is the vertical unit vector, τ_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

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{ij}=\mathrm{2}{\mathit{\eta }}_{ij}{\dot{\mathit{\epsilon }}}_{ij}+\left(\mathit{\zeta }-\mathit{\eta }\right){\dot{\mathit{\epsilon }}}_{kk}{\mathit{\delta }}_{ij}-{\displaystyle \frac{P}{\mathrm{2}}}{\mathit{\delta }}_{ij},\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\dot{\mathit{ϵ}}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right),\end{array}$$

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:

$$\begin{array}{}\text{(4)}& P=P^{\star }h{e}^{-C^{\star }(\mathrm{1}-A)},\end{array}$$

where C^⋆ is a free parameter.

The stress tensor σ is often expressed in terms of principal stresses σ₁ and σ₂ or stress invariants σ_I and σ_II. The principal stresses σ₁ and σ₂ 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 σ₁ and σ₂ can be expressed as a function of σ_ij as

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}+\sqrt{({\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}{)}^{\mathrm{2}}+\mathrm{4}{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}}\right),\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}-\sqrt{({\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}{)}^{\mathrm{2}}+\mathrm{4}{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}}\right).\end{array}$$

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

$$\begin{array}{}\text{(7)}& \tan \left(\mathrm{2}\mathit{\psi }\right)={\displaystyle \frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{12}}}{{\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}}}.\end{array}$$

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

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{I}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{\mathrm{11}}+{\mathit{\sigma }}_{\mathrm{22}}),\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{II}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{({\mathit{\sigma }}_{\mathrm{11}}-{\mathit{\sigma }}_{\mathrm{22}}{)}^{\mathrm{2}}+\mathrm{4}{\mathit{\sigma }}_{\mathrm{12}}^{\mathrm{2}}}.\end{array}$$

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

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }={\displaystyle \frac{P}{\mathrm{2}\mathrm{\Delta }}},\text{(11)}& {\displaystyle }& {\displaystyle }\mathit{\eta }={\displaystyle \frac{\mathit{\zeta }}{e^{\mathrm{2}}}},\end{array}$$

with the abbreviation

$$\begin{array}{}\text{(12)}& {\displaystyle }\mathrm{\Delta }=\sqrt{{\dot{\mathit{ϵ}}}_{\mathrm{I}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{e^{\mathrm{2}}}}{\dot{\mathit{ϵ}}}_{\mathrm{II}}^{\mathrm{2}}}.\end{array}$$

In this abbreviation, the strain rate invariants are the divergence ${\dot{\mathit{ϵ}}}_{\mathrm{I}}={\dot{\mathit{ϵ}}}_{\mathrm{11}}+{\dot{\mathit{ϵ}}}_{\mathrm{22}}$, and the maximum shear strain rate ${\dot{\mathit{ϵ}}}_{\mathrm{II}}=\sqrt{({\dot{\mathit{ϵ}}}_{\mathrm{22}}-{\dot{\mathit{ϵ}}}_{\mathrm{11}}{)}^{\mathrm{2}}+\mathrm{4}{\dot{\mathit{ϵ}}}_{\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:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{MC}}=min\left\{\mathit{\eta },{\displaystyle \frac{\mathrm{1}}{{\dot{\mathit{ϵ}}}_{\mathrm{II}}}}\left[\mathit{\mu }\left({\displaystyle \frac{P}{\mathrm{2}}}-\mathit{\zeta }\cdot {\dot{\mathit{ϵ}}}_{kk}\right)-c\right]\right\},\end{array}$$

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^⋆.

Figure 1 Elliptical yield curve (black) with ellipse aspect ratio $e=a/b=\mathrm{2}$. Coulombic yield curve (red) and elliptical capping with internal angle of friction (μ). Both e and μ are measures of the shear strength of the material. The normal flow rule applies only to the elliptical part of the yield curves. For the two straight limbs of the Coulombic yield curve, the flow is normal to the truncated ellipse (dashed-dotted line) with the same first stress invariant. Note that the axes σ₁, σ₂ and σ_I, σ_II do not have the same scale.

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.

Figure 2 Model domain with a solid wall on the southern (red) boundary (Dirichlet boundary conditions with u=0) and prescribed southward velocities on the northern orange boundary (u=0, $v=a_v\cdot t+v_{\mathrm{i}}$; Eq. 9) and open boundaries to the east and the west (green) with von Neumann boundary conditions. θ is the measured fracture angle with the blue line representing an LKF.

Table 1 Model parameters of the reference simulation.

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{ϵ}}_{\mathrm{LSR}}={\mathrm{10}}^{-\mathrm{11}}\mathrm{m}{\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:

$$\begin{array}{}\text{(14)}& {\left.{\displaystyle \frac{\partial u}{\partial x}}\right|}_{\mathrm{E},\mathrm{W}}={\left.{\displaystyle \frac{\partial v}{\partial x}}\right|}_{\mathrm{E},\mathrm{W}}={\left.{\displaystyle \frac{\partial A}{\partial x}}\right|}_{\mathrm{E},\mathrm{W}}={\left.{\displaystyle \frac{\partial h}{\partial x}}\right|}_{\mathrm{E},\mathrm{W}}=\mathrm{0},\end{array}$$

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:

$$\begin{array}{}\text{(15)}& {\left.v\right|}_{\mathrm{N}}\left(t\right)=a_v\cdot t+v_{\mathrm{i}};{\left.u\right|}_{\mathrm{N}}=\mathrm{0};{\left.{\displaystyle \frac{\partial A}{\partial y}}\right|}_{\mathrm{N}}={\left.{\displaystyle \frac{\partial h}{\partial y}}\right|}_{\mathrm{N}}=\mathrm{0},\end{array}$$

where a_v is the prescribed acceleration, and N denotes the northern boundary.

3 Results

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{ϵ}}}_{\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 }$.

Figure 3 (a) First and (b) second strain invariants, (c) ice thickness anomaly ($\mathrm{\Delta }h=h-\mathrm{1}$), and (d) stress states in normalized stress invariant space along with the elliptical yield curve after 5 s of integration. The first and second strain invariants represent the divergence and maximum shear strain rate, respectively. The modeled angle of fracture 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 σ₂ 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.

Figure 4 Schematic of stress states and failure in principal stress space. Black arrows show how stresses move from zero at the beginning of loading towards the yield curve until failure. Red points show the stress states at failure – the intersection point between the second principal axis 2 (in red) and the elliptical yield curve – for different ellipse ratios e=2, 1, 0.7. The red arrows show the direction of deformation with a normal flow rule. The blue points and arrows show the case when the ice floe is confined and the loading will lead to extra stress in the direction of σ₁.

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).

Figure 5 Maximum shear strain rate (second strain invariant) after 10 s of integration for the default domain size and Δx=100 m (a) and 500 m (b) and for the default Δx and a doubled domain size of 20 km×50 km (c). Note that for the case of the double domain (c), the southward velocity at the northern boundary was also doubled to keep the deformation rate constant and that this simulation is limited to 2 s for numerical efficiency.

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.

Figure 6 Sea ice thickness (a), concentration (b), maximum shear strain rate (c), and divergence (d) after 45 min of integration (2700 s) in a uni-axial loading test. To make these longer simulations possible, both nonlinear and linear iterations are limited to 150 per time step. Results show the development of secondary fracture lines in all fields after the first fracture line has formed.

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.

Figure 7 Maximum shear strain rate after 5 s of integration in a reduced size domain (8 km×10 km) with free-slip (a) and no-slip (b) boundary conditions. Note that the no-slip boundary condition forces the fracture to occur at the corner of the domain, leading to a larger angle of $\mathit{\theta }=\mathrm{39}{}^{\circ }$ versus $\mathrm{34}\pm \mathrm{1}{}^{\circ }$ in the control experiment. This suggests that the choice of boundary conditions in current sea ice models needs to be revisited.

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 }=\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 σ₂ axis but with a non-zero σ₁ 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 σ₁ along a line of constant σ₂ (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.

Figure 8 Maximum shear strain rates (left) and stress state in stress invariant space (right) after 5 s of integration for different confinement pressure: h_c=0.05 m (a) and h_c=0.3 m (b). Note how stress states with divergent strain rates (a) migrate left towards convergent strain rates (b).

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).

Figure 9 Sea ice thickness with two ice-free areas (a) and maximum shear strain rates for two different ellipse aspect ratios (b, c) after 5 s of integration. The position of the ice weaknesses determines the location and angle of the fracture lines, and also the rheology parameter e has an entirely different effect. The main fractures lines are at angles of 25 and 34^∘ for e=2.0 and 57.6^∘ for e=1.0.

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⁻¹ 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 σ₂ 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 σ₂ 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.

Figure 10 Fracture angles as a function of ellipse aspect ratio e with constant P^⋆ (red, bottom scale; Sect. 3.3.1). The theoretical relationship ${\mathit{\theta }}_{\mathrm{th},\mathrm{ell}}=\frac{\mathrm{1}}{\mathrm{2}}\arccos \left[\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}-\frac{\mathrm{1}}{e^{\mathrm{2}}}\right)\right]$ (dashed black curve; Eq. B3 in the Appendix) fits the modeled angles almost perfectly with R²=0.9995 and $\sqrt{\mathrm{VAR}}=\mathrm{0.089}$. The simulated fracture angles for the Coulombic yield curve as a function of the slope of the Mohr–Coulomb limbs (blue, top scale; Sect. 3.3.2) fit the theoretical relationship ${\mathit{\theta }}_{\mathrm{th},c}=\frac{\mathrm{1}}{\mathrm{2}}\arccos \left(\mathit{\mu }\right)$ only for μ≤0.7 (black line; Eq. B4 in the Appendix). The errors bars mean that there was more than one unique fracture line: for a small μ, the ice breaks easily along the lateral edges of the floe. For μ>0.7 ($\mathit{\varphi }=\mathrm{44}{}^{\circ }$), the ambiguity appears because the stress states are both on the linear limbs and on the elliptical cap. For μ≥0.9 (blue line), the fracture angle is the same as for the ellipse for e=1.4.

Figure 11 Maximum shear strain (a), ice thickness anomaly (b), divergence (c), and stress state in stress invariant space (d) after 5 s of integration for a smaller ellipse aspect ration (e=0.7 compared to e=2 in the reference run in Sect. 3.1). Compared to the control run on Fig. 3, the angle of fracture is larger ($\mathit{\theta }=(\mathrm{61}\pm \mathrm{1}){}^{\circ }$), the stress states are in the second half of the ellipse (with strain rates pointing into the convergent direction), and there is convergence along the fracture lines (b) in agreement with the schematic in Fig. 4.

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.

Figure 12 Maximum shear strain (top) and stress state in stress invariant space (bottom) for different internal angles of friction. (a) μ=0.7 or $\mathit{\varphi }=\mathrm{44}{}^{\circ }$, (b) μ=0.85 or $\mathit{\varphi }=\mathrm{58}{}^{\circ }$ and (c) μ=0.95 or $\mathit{\varphi }=\mathrm{72}{}^{\circ }$ after 5 s of integration. The angles of fracture are θ=23, (28±2) and 41^∘. Figure 10 illustrates how θ depends on μ for a Coulombic yield curve.

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 σ₂ 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 σ₂ 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.

4 Discussion

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⁻¹, 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 σ₂ 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).

5 Conclusions

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.

Appendix A: Fracture angle

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 (σ₁,σ₂) 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)

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }\mathit{\sigma }\mathrm{d}A={\mathit{\sigma }}_{\mathrm{2}}\sin \left(\mathit{\theta }\right)\sin \left(\mathit{\theta }\right)\mathrm{d}A+{\mathit{\sigma }}_{\mathrm{1}}\cos \left(\mathit{\theta }\right)\cos \left(\mathit{\theta }\right)\mathrm{d}A,\text{(A2)}& {\displaystyle }& {\displaystyle }\mathit{\tau }\mathrm{d}A=-{\mathit{\sigma }}_{\mathrm{2}}\cos \left(\mathit{\theta }\right)\sin \left(\mathit{\theta }\right)\mathrm{d}A+{\mathit{\sigma }}_{\mathrm{1}}\cos \left(\mathit{\theta }\right)\sin \left(\mathit{\theta }\right)\mathrm{d}A,\end{array}$$

where dA 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 σ₁ and σ₂. Using THE angle sum and difference identities of trigonometry, we can write the stresses σ and τ in terms of the principal stresses σ₁ and σ₂ as

$$\begin{array}{}\text{(A3)}& {\displaystyle }& {\displaystyle }\mathit{\sigma }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}})+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}})\cos \left(\mathrm{2}\mathit{\theta }\right),\text{(A4)}& {\displaystyle }& {\displaystyle }\mathit{\tau }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}})\sin \left(\mathrm{2}\mathit{\theta }\right).\end{array}$$

In terms of the stress invariants σ_I and σ_II, this gives

$$\begin{array}{}\text{(A5)}& {\displaystyle }& {\displaystyle }\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{I}}+{\mathit{\sigma }}_{\mathrm{II}}\cos \left(\mathrm{2}\mathit{\theta }\right),\text{(A6)}& {\displaystyle }& {\displaystyle }\mathit{\tau }={\mathit{\sigma }}_{\mathrm{II}}\sin \left(\mathrm{2}\mathit{\theta }\right).\end{array}$$

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

$$\begin{array}{}\text{(A7)}& \mathit{\tau }=-\tan \left(\mathit{\varphi }\right)\mathit{\sigma }+c,\end{array}$$

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

$$\begin{array}{}\text{(A8)}& {\mathit{\sigma }}_{\mathrm{II}}\sin \left(\mathrm{2}\mathit{\theta }\right)=-\tan \left(\mathit{\varphi }\right){\mathit{\sigma }}_{\mathrm{I}}-\tan \left(\mathit{\varphi }\right){\mathit{\sigma }}_{\mathrm{II}}\cos \left(\mathrm{2}\mathit{\theta }\right)+c,\end{array}$$

and after multiplying both sides by cos (ϕ),

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{II}}\left[\sin \left(\mathrm{2}\mathit{\theta }\right)\cos \left(\mathit{\varphi }\right)+\cos \left(\mathrm{2}\mathit{\theta }\right)\sin \left(\mathit{\varphi }\right)\right]=\\ \text{(A9)}& {\displaystyle }& {\displaystyle }-{\mathit{\sigma }}_{\mathrm{I}}\sin \left(\mathit{\varphi }\right)+c\cos \left(\mathit{\varphi }\right).\end{array}$$

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

$$\begin{array}{}\text{(A10)}& {\mathit{\sigma }}_{\mathrm{II}}=-{\mathit{\sigma }}_{\mathrm{I}}\sin \left(\mathit{\varphi }\right)+c\cos \left(\mathit{\varphi }\right),\end{array}$$

so that Eq. (A3) becomes

$$\begin{array}{}\text{(A11)}& \sin \left(\mathrm{2}\mathit{\theta }\right)\cos \left(\mathit{\varphi }\right)+\cos \left(\mathrm{2}\mathit{\theta }\right)\sin \left(\mathit{\varphi }\right)=\sin (\mathrm{2}\mathit{\theta }+\mathit{\varphi })=\mathrm{1},\end{array}$$

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

$$\begin{array}{}\text{(A12)}& \mathrm{2}\mathit{\theta }+\mathit{\varphi }={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\Rightarrow \mathit{\theta }={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}-{\displaystyle \frac{\mathit{\varphi }}{\mathrm{2}}}.\end{array}$$

Figure A1 Stress state in physical stress space (a) and in an arbitrary coordinate system oriented at an angle θ with respect to the principal stress axes (b). The principal stresses are the eigenvalues of the stress tensor in an arbitrary coordinate system, and the angle ψ is derived from the rotation matrix composed of the two eigenvectors. Note that in the study above there is no shear stress (σ12=0, so principal axes and physical axes are aligned (ψ=0).

Figure A2 Mohr's circle of stress (black) with the Mohr–Coulomb yield criterion (red) of the angle of internal friction ϕ (red) and cohesion c in (σ,τ) space. From Eq. (A5), the deformation is created with an angle θ that can be represented in Mohr's circle (blue).

Appendix B: Fracture angle and yield curve

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)

$$\begin{array}{}\text{(B1)}& \sin \left(\mathit{\varphi }\right)=\tan \left(\mathit{\gamma }\right)=\mathit{\mu }={\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{II}}}{\partial {\mathit{\sigma }}_{\mathrm{I}}}}.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }\mathit{\theta }\left({\mathit{\sigma }}_{\mathrm{I}}\right)& {\displaystyle }={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\arcsin \left({\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{II}}}{\partial {\mathit{\sigma }}_{\mathrm{I}}}}\left({\mathit{\sigma }}_{\mathrm{I}}\right)\right)\\ \text{(B2)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\arccos \left({\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{II}}}{\partial {\mathit{\sigma }}_{\mathrm{I}}}}\left({\mathit{\sigma }}_{\mathrm{I}}\right)\right).\end{array}$$

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 }}=\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 σ₂ axis, we get

$$\begin{array}{}\text{(B3)}& {\left.{\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{II}}}{\partial {\mathit{\sigma }}_{\mathrm{I}}}}\right|}_{{\mathit{\sigma }}_{\mathrm{1}}=\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{e^{\mathrm{2}}}}\right).\end{array}$$

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

$$\begin{array}{}\text{(B4)}& {\mathit{\theta }}_{\mathrm{th},\mathrm{ell}}\left(e\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\arccos \left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{e^{\mathrm{2}}}}\right)\right].\end{array}$$

Figure B1 Illustration of the Mohr's circle applied to the elliptical yield curve (black ellipse) in σ,τ space, some examples of Mohr's circles (blue), and the reconstructed yield curve (red) in the fracture plane space. The orange Mohr's circle illustrates the case in which no fracture lines exists, for $\left|\mathit{\mu }\right|>\mathrm{1}$.

Figure B2 Mohr's circle of stress with an arbitrary yield curve (black line) in the fracture plane reference. tan (γ)=μ is the tangent to the yield curve, and ϕ is the internal angle of friction as described in Appendix Appendix A. We note that $\sin \left(\mathit{\varphi }\right)=\tan \left(\mathit{\gamma }\right)=\mathit{\mu }$ (for $\left|\mathit{\mu }\right|\le \mathrm{1}$). For a slightly different Mohr's circle (grey), the blue and red tangents meet in the same point on the σ axis.

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

$$\begin{array}{}\text{(B5)}& {\mathit{\theta }}_{\mathrm{th},c}\left(\mathit{\mu }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\arccos \left(\mathit{\mu }\right),\end{array}$$

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.

Figure B3 Mohr's circle applied to the Coulombic yield curve (in black) in σ,τ space, the Mohr's circle for the cusps between the elliptical cap and the Mohr–Coulomb linear limbs (blue circle), and the yield curve in (σ,τ) space (red). We can see the effect of combining two regimes; for the same Mohr's circle, two different angles coexist (red circles) and are apart from each other.