## Abstract

The drift and deformation of sea ice floating on the polar oceans is caused by the applied wind and ocean currents. Over ocean basin length scales the internal stresses and boundary conditions of the sea ice pack result in observable deformation patterns. Cracks and leads can be observed in satellite images and within the velocity fields generated from floe tracking. In a climate sea ice model the deformation of sea ice over ocean basin length scales is modelled using a rheology that represents the relationship between stresses and deformation within the sea ice cover. Here we investigate the link between emergent deformation characteristics and the underlying internal sea ice stresses using the Los Alamos numerical sea ice climate model. We have developed an idealized square domain, focusing on the role of sea ice rheologies in producing deformation at spatial resolutions of up to 500 m. We use the elastic anisotropic plastic (EAP) and elastic viscous plastic (EVP) rheologies, comparing their stability, with the EAP rheology producing sharper deformation features than EVP at all space and time resolutions. Sea ice within the domain is forced by idealized winds, allowing for the emergence of five distinct deformation types. Two for a low confinement ratio: convergent and expansive stresses. Two about a critical confinement ratio: isotropic and anisotropic conditions. One for a high confinement ratio and isotropic sea ice. Using the EAP rheology and through the modification of initial conditions and forcing, we show the emergence of the power law of strain rate, in accordance with observations.

This article is part of the theme issue ‘Modelling of sea-ice phenomena’.

## 1. Introduction

Sea ice floating on the polar oceans is composed of many individual floes and floe aggregates. These floes are continually breaking apart, sliding against one another and thermodynamically healing. These interactions result in discontinuous drift and deformation fields. Sea ice is continually under the influence of external wind and ocean forcing that, along with the internal ice stresses resulting from floe interaction and deformation and Coriolis acceleration, make up the sea ice force balance [1,2].

Accurately representing the sea ice force balance is essential for modelling the drift and thickness distribution of sea ice over seasonal and climate time scales. The resulting sea ice drift is important in the momentum transfer from the atmosphere, through sea ice, to the ocean below, driving basin-scale currents in the Arctic Ocean.

The relation between stress and deformation of sea ice is known as sea ice rheology. When modelling the drift and deformation of sea ice within a numerical climate model, a numerical implementation of the rheology of sea ice over ocean dynamical length scales is needed. In recent years, new anisotropic rheologies have been developed that take inspiration from observable oriented features in satellite observations, for example, the elastic anisotropic plastic (EAP) rheology used in this paper [3] and the Maxwell elasto-brittle rheology [4]. The EAP rheology has been shown to produce significantly different sea ice drift and thickness distribution patterns in the Arctic Ocean compared to the isotropic elastic viscous plastic (EVP) rheology [2].

The deformation and force balance of 10 cm samples of sea ice have been investigated in laboratory settings [5]. Typically, these samples are compressed in two dimensions until they fail. The ratio of the magnitude of stress in two dimensions, the stress confinement ratio *R*_{int} (as described in §2a(i)), is found to control the mechanical failure of the ice. For *R*_{int} below a critical confinement ratio *R*_{crit}, ice fails with Coulombic lines of shear deformation, lines that are inclined to the principal components of the stress tensor. For *R*_{int} > *R*_{crit} there is a different failure mechanism, with either a single unoriented feature or the disintegration of the sample [6]. These deformation characteristics are argued to be present in sea ice over all length scales [7,8], which is the underlying physical principle behind the EAP rheology presented in §2a(i).

The drift and deformation of sea ice have been observed both from floating arrays of remote buoys [9], and from comparing images from orbiting satellites [10]. The distribution of deformation from these sources is proposed to follow a power law distribution [9,11]. This distribution describes the probability density of the magnitude of strain rate over the entire sea ice pack, with likelihood of a strain rate proportional to that strain rate to a power *r*_{c} < 0. This distribution has been used to investigate how well climate sea ice models represent the deformation and thus the rheology of sea ice over climate length and time scales [12], with higher-resolution models typically better able to reproduce the power law distribution.

In this paper, we document our study into the link between observable deformation features and the underlying internal stress characteristics. We use a state of the art sea ice climate model, the Los Alamos sea ice model (CICE), with the EAP and EVP rheologies on an idealized domain, documented in §2a,b. We use idealized wind forcing in order to impose stress states within the sea ice that produce particular deformation patterns, §2c. We explore the model responses at high resolution, §2d. We explore the link between stress confinement in the wind and internal ice stresses and the resulting deformation patterns, §3b. In the case of EAP, we look at the effect of the anisotropic alignment of the sea ice structure upon deformation characteristics, particularly the deformation characteristics of reorienting winds and sea ice structural alignment, §3c.

## 2. Model configuration

### (a) Model equations

The numerical sea ice model CICE calculates the drift and deformation of sea ice using the vertically integrated, horizontal momentum balance [13]:
2.1The left-hand side of the equation represents the rate of change of momentum, with **u** the sea ice drift velocity and *m* the mass of sea ice per unit area, balanced by, in order: the Coriolis acceleration −*mf*_{c}**k** × **u**, with *f*_{c} the Coriolis parameter and **k** the unit vector normal to the sea surface; the ice concentration *C* weighted applied stress from the atmosphere *τ*^{a} and ocean *τ*^{o}; gravitational acceleration from the ocean tilt ** S** set to zero for this study; and the divergence of the internal ice stress tensor ∇ ·

**.**

*σ*The atmospheric drag is calculated from the 10 m wind speed *U*_{a} with *τ*^{a} = *ρ*_{a}*C*_{a}|*U*^{a}|*U*^{a}, where *ρ*_{a} is the air density and *C*_{a} is a drag parameter. Ocean stresses are typically calculated from the difference between ice and ocean velocities along with a turning angle. In this study, we have a still ocean and no turning angle with *τ*^{o} = − *ρ*_{o}*C*_{o}|** u**|

**|, where**

*u**ρ*

_{o}is the ocean density and

*C*

_{o}an ocean drag coefficient.

#### (i) Elastic anisotropic plastic rheology

To represent floe-scale interactions within a continuum model, Wilchinsky & Feltham [3] developed the EAP rheology that sums together the forces arising between many diamond-shaped floes within an arbitrary area of sea ice cover. This rheology was implemented into the Los Alamos sea ice model CICE by Tsamados *et al.* [14], with further investigations into its role in the sea ice force balance presented by Heorton *et al.* [2]. The orientation of ice floes within the area is recorded using a model variable that changes to represent the breaking and healing of ice floes. When considering the interaction between individual floes within an arbitrary area, one is able to derive the plastic sliding and ridging stresses that are combined to give the full local ice stress:
2.2The local sliding, *σ*^{f}_{s}, and ridging, *σ*^{f}_{r}, stresses are obtained for the floe orientations (given by unit vector ** r**) and their relative motion (determined from strain rate ) as described by Wilchinsky & Feltham [3].

*P*

_{(r,s)}(

*h*) are the ridging and sliding ice strengths for ice of thickness

*h*. The individual floe stresses are expanded over an arbitrary area containing many floe alignments represented by

*ψ*(

**) =**

*r**ψ*( −

**) (with ) and total internal stress given by where the sliding ice strength**

*r**P*

_{s}=

*kP*

_{r}with

*k*a constant.

The structure tensor, used to represent floe orientation, is similarly defined as
2.3which is a tensor of unit trace, i.e. *A*_{11} + *A*_{22} = 1 with eigenvalues *A*_{1}, *A*_{2} = 1 − *A*_{1} associated with eigenvector *A*_{1}, the principal component of ** A** giving the direction of local anisotropic alignment

**. The largest eigenvalue 0.5 <**

*r**A*

_{1}< 1 indicates the level of local anisotropy, with

*A*

_{1}= 1 being fully anisotropic and

*A*

_{1}= 0.5 fully isotropic. The structure tensor changes in time due to local forcing with 2.4where D

^{c}/D

^{c}

*t*is the co-rotational time derivative accounting for advection and rigid body rotation and the

*F*_{( )}terms represent various processes that realign ice floes.

*F*_{frac} determines the ice floe reorientation due to fracture, and explicitly depends upon sea ice stress direction but not its magnitude. Following Wilchinsky & Feltham [3], we use four failure modes defined by the internal stress confinement ratio *R*_{int} = *σ*_{1}/*σ*_{2}, where *σ*_{1,2} are the principal components of the stress tensor: (i) biaxial tension causing longitudinal splitting; (ii) uniaxial compression/tension causing axial splitting; (iii) biaxial compression with a low confinement ratio causing in-plane shear rupture; and (iv) biaxial compression with a large confinement ratio causing out-of-plane shear rupture. Modes (i), (iv) cause no realignment of floes and modes (ii), (iii) align floes with ** r** parallel to the direction of greatest compressive stress (

*σ*

_{2}). The formulation for

*F*_{frac}is 2.5where

**is a tensor that results in the major principal axis of**

*S***aligning with the principal direction of**

*A***associated with**

*σ**σ*

_{1};

*k*

_{f}reflects the rate of fracture formation in the sea ice cover and, as with [14], we choose the reference value of

*k*

_{f}= 10

^{−3}s

^{−1}which gives 90% anisotropic alignment within 6 h of case (ii) or (iii) occurring. The value of

*R*

_{c}= 0.3 used is in accordance with the laboratory-scale observations of Golding

*et al.*[6] and Iliescu & Schulson [5]. The thermodynamic healing of the sea ice structure tensor is turned off (

*F*_{iso}=

**0**) to allow us to focus on floe reorientation due to fracture

*F*_{frac}.

#### (ii) Elastic viscous plastic rheology

The EVP rheology [15] is a numerical implementation that elastically approximates the viscous plastic (VP) rheology of Hibler [16]. In this rheology, the deformation of the sea ice is modelled as plastic for high stress states and as a highly viscous fluid for low stress states to ease the numerical complexity of distinguishing between plastic and non-plastic deformation (see [13] for further description). This is represented through the stress tensor with
2.6aand
2.6bwhere *η* and *ζ* are the shear and bulk viscosities, *e* defines the elliptical plastic yield curve aspect ratio, is a scaling factor representing the magnitude of local strain and *p* is the ice strength as discussed below.

### (b) Model domain, boundary and initial conditions

We use a square grid with constant Coriolis acceleration of *f*_{c} = 1.46 × 10^{−4} s^{−1} to make the model applicable to the polar regions. The size of the domain is *d* = 2000 km square. To simplify the model dynamics, we use a single thickness category and the ice strength parametrization of Hibler [16] with where *p** = 2700 N m^{−2} and *c* = 20 are constants. Tests showed that using the alternative strength parametrization of Rothrock [17] and five thickness categories made no qualitative difference to our simulation results. The ice strength parametrization used is the same as that used by Hutchings *et al.* [18], though we have differences in the way the strength and thus internal stresses are initialized, as discussed below.

To simulate the stress characteristics within the continuous sea ice pack special boundary conditions are needed (figure 1). Boundary conditions for open or closed boundaries cause either too little or too much stress and thus no deformation features, due to dissipation of stress or the ‘locking up’ of the sea ice pack, respectively. For this reason, we use the boundary conditions used in Hutchings *et al.* [18] to produce realistic deformation characteristics on a square grid. For a buffer region of *d*_{bc} ≈ 50 km from the ice edge (tuned for each model resolution to give realistic deformation characteristics, figure 1), we force the sea ice to be in a quasi-free drift state with spatially uniform stress and strain. After stress equations are solved within the elastic sub-cycle (described fully within Hunke *et al.* [19]), giving realistic stress conditions across the whole domain, we stop the internal ice stresses from dissipating out of the domain by imposing the sea ice drift speed within the buffer region (*d*_{bc}) by balancing the wind forcing and a linear ocean drag with
2.7where *U*^{o}_{g} = 0 represents the still ocean. As experienced by Hutchings *et al.* [18] this boundary condition gives uniform stress and strain rate within the buffer region as seen in the later figure in §3a. The model simulations presented in this paper are numerically stable throughout. Over a transition region of size *d*_{trans} = 300 km from the domain edge, the uniform strain rate changes into observable features. The inner region beyond *d*_{trans} is the area analysed within the results sections (figure 1).

The model is initialized using data from an Arctic-wide CICE run using the restart configuration. There are many variables required to restart a model run. Most we take from a single point in the Arctic-wide run with 2% noise, with idealized: sea ice velocity near zero, sea ice concentration of 0.999, a constant thin snow cover of 0.1 m and constant idealized sea ice structure tensor and thickness, which we discuss in §3. The rest of the CICE model is unmodified with the thermodynamic and thickness distribution equations solved for. Our model set-up is in contrast to Hutchings *et al.* [18] who solve the VP rheology over a uniform ice cover. So as to allow discrete kinematic features to form, they also introduce noise to the otherwise continuous ice field, but into the initialization of ice strength only.

### (c) Idealized forcing

In order to gain insight into the role of sea ice rheology in producing sea ice deformation features, we perform simulations with idealized atmospheric and oceanic forcing. To induce near constant internal stress conditions over the idealized domain, we impose a wind forcing similar to that used both by Wilchinsky *et al.* [20,21] for discrete element simulations of sea ice and by Hutchings *et al.* [18]. In all model runs, the ocean velocity *U*^{o}_{g} is set to zero. We construct wind forcing by imposing the gradient of the wind stress, with *τ*^{a}_{yy}/*τ*^{a}_{xx} = *R*_{wind} a constant and the maximum wind speed **U**^{a}_{max}. The wind velocity field can be generated from and the boundary conditions of (*u*^{a}, *v*^{a}) = (**U**^{a}_{max}, 0)|^{y=d/2}_{x=0} and (*u*^{a}, *v*^{a}) = ( − **U**^{a}_{max}, 0)|^{y=d/2}_{x=d}, measuring (*x*, *y*) from the southwest corner of the model domain with
2.8This wind pattern has winds heading eastwards from the western edge of the domain and westwards from the eastern edge of the domain. The winds diverge out of the northern and southern edges for *R*_{wind} < 0 (see the left-hand side of figure 1), and head southwards from the northern edge of the domain and north from the southern edge of the domain for *R*_{wind} > 0 (right-hand side of figure 1). The velocity field is symmetric about the lines *x* = *d*/2 and *y* = *d*/2. The model is initialized with a zero wind field (*u*^{a}, *v*^{a}) = **0** that increases linearly to the idealized wind forcing set over 6 h of model time.

### (d) Model stability

#### (i) Resolution

The model resolution has been tested at 10 km, typical of modern high-resolution global climate models, 2 km, typical of very high-resolution regional models, and 500 m. The time resolution is decreased for finer spatial resolution to allow the equations of motion to be solved, with a time step of 600 s at 10 km, 30 s at 2 km and 5 s at 500 m. For all cases, we use 200 elastic sub-cycles within the CICE numerical solver. These time steps give model convergence at each resolution with the model results not changing by more than 1% for shorter time steps or increased sub-cycles. The model responses have been widely tested to confirm a correct response. For example, for the 2 km EVP runs presented in figure 2, the same result is achieved for a model time step of 5 s and 1000 sub-cycles.

The deformation patterns for the EAP rheology in figure 2 show good correlation for increasing resolution. The linear deformation features have a similar shape for the three resolutions. At increasing resolution, finer details can be observed: finer parallel slip lines for *R*_{wind} = 0.4 and the appearance of comb cracks at 2 km and 500 m resolution for *R*_{wind} = − 0.4. When using the EVP rheology (figure 2) only the 500 m resolution under winds with *R*_{wind} = − 0.4 gives linear deformation features after only 6 h of model time.

We have repeated many of the runs presented in this paper with the wind field rotated by 45° about the centre of the domain to test the numerical sensitivity of our chosen square grid. The characteristic lines of deformation and the principal component of the structure tensor rotate with the forcing field as seen in the rotating wind experiments in the penultimate figure in this paper.

#### (ii) Model run length

The sea ice model has been run for 10 days of model time at 2 km resolution under winds with *R*_{wind} = − 0.8 (shown in figure 3), chosen to produce high rates of deformation and minimal mechanical thickening at the centre of the domain. We test whether the emergent deformation characteristics are stable or change with time.

For the EAP rheology oriented lines of deformation are observable from 3 h of model time, and the sea ice is strongly aligned (*A*_{1} > 0.95) from after 10 h of model time. From 12 h to 10 days, the features have little variation. All differences between these two time points are due to thickening sea ice, primarily in areas of ridging.

With the EVP rheology, the model takes longer to stabilize. There are no linear deformation features at 6 h and few at 12 h. Identifiable evenly spaced lines of deformation are apparent from 3 days of model run and remain until the end of the run. Using a higher resolution may give more identifiable features (see figure 2), but is computationally beyond the scope of this study.

## 3. Results

### (a) Reference runs with the EAP rheology, initially isotropic

Here we present the stress–strain rate relationship for two reference runs with an initial isotropic sea ice cover (*A*_{1} = 0.5): the case of *R*_{wind} = − 0.8 and divergent shear deformation (figure 4) and of *R*_{wind} = 0.8, both shown in figure 5. The wind speed increases linearly over the first 6 h of model run, with the wind stress and total strain rate being proportional to the square of the wind speed; see dark red and dark blue lines in figure 5. The internal stress magnitude increases to a steady value by 3 h of model run. The internal stress confinement *R*_{int} remains near constant for *R*_{wind} = 0.8 from 3 h onwards. For *R*_{wind} = − 0.8, *R*_{int} is increasing from 3 h to the end of the run as the sea ice becomes increasingly anisotropically aligned ().

In figure 4 we show the shear and divergent strain rate and internal stress for the run with *R*_{wind} = − 0.8. The internal stresses are near spatially uniform with lines of reduction in stress corresponding to lines of high deformation, and reduced or increased stress near the domain boundaries. The lines of deformation have high shear strain rate and either divergent or compressive failure. For all the runs with higher *R*_{wind}, divergent deformation is rare and the lines of shear only correspond to compressive failure. In the centre of the domain, far enough from the imposed boundary conditions, there is near uniform stress and the lines of deformation are evenly spaced. The confinement of the internal stress is less than the critical confinement ratio *R*_{int} < *R*_{crit} for the whole run, causing the sea ice structure to become anisotropically aligned as indicated by the overlaying quadrilaterals in figure 4 being diamond shaped. The principal component of the sea ice structure tensor (*A*_{1}, as indicated by the major axis of the diamonds) is aligned to the principal component of the wind stress gradient, not the wind direction. This is in agreement with Heorton *et al.* [2] who found highly variable angles between the direction of wind stress and sea ice structure within the central Arctic.

For the convergent wind field with *R*_{wind} = 0.8 shown in figure 6 the compressive stress and strain rate dominates. The centre of the domain is dominated by a region of near constant compressive stress and deformation with little shear stress. Radiating away from the central region there are slip lines of shear and compressive deformation that correspond to reduced compressive and increased shear stresses. The majority of the sea ice remains isotropic, with the sea ice within slip lines becoming anisotropically aligned (not shown).

### (b) Varying the internal stress confinement *R*_{int}, initially isotropic

Here we investigate the relationship between the internal stress confinement ratio *R*_{int} and the deformation characteristics of sea ice displayed in figures 6 and 5. Since *R*_{int} is an emergent property we cannot directly impose it and so we perform a set of simulations with different values of imposed *R*_{wind}, which directly determines *R*_{int}. We calculate *R*_{int} = *σ*_{1}/*σ*_{2} from the principal components or eigenvectors of the internal ice stress tensor *σ* for each grid cell point, with the values plotted in figure 5 taken from within the inner region *d*_{trans} from the domain edge.

#### (i) Stress and strain characteristics

As with the runs presented in §3a, the wind speed increases to **U**^{a}_{max} = 15 m s^{−1} over the first 6 h of the model run. The magnitude of internal stress reaches its maximum from 3 h of model run whereas the strain rate increases with the wind speed to its maximum at 6 h of model run (see figure 5). The runs with higher *R*_{wind} have greater magnitude of internal stress due to greater confinement. The magnitude of strain rate appears to be equal for the wind confinements *R*_{wind} of equal magnitude (+0.8, − 0.8, for example), due to the symmetry of equations (2.8), resulting in the total applied wind stress for these runs being equal.

In the EAP rheology, the runs that have *R*_{int} < *R*_{crit} become anisotropically aligned (see the *A*_{1} and *R*_{int} panels of figure 5). For increasing alignment within the sea ice cover *R*_{int} coverges to *R*_{crit} = 0.3.

The experiments with *R*_{wind} < 0.1 have curved lines for the probability density function (pdf) of shear strain rate (figure 5), indicating the deformation taking place over many low strain rate events. For *R*_{int} > *R*_{crit} there is a complex curve showing the deformation being collected into fewer, high strain rate features. The combined pdf for all runs (black dashed line) approaches a linear relationship and thus the power law of sea ice deformation suggested by Girard *et al.* [12] and Rampal *et al.* [9].

#### (ii) Emergent failure modes

For an initial isotropic cover of sea ice, varying *R*_{wind} reveals two failure modes of sea ice. The runs with a large positive *R*_{wind} in the range 0.1 < *R*_{wind} < 0.8 produce shear strain rate concentrated into near parallel lines that diverge towards the edge of the domain; see figure 6. The high shear deformation coincides with convergent deformation. The lines are perpendicular to the principal component of the wind stress gradient. This mode we attribute to the out-of-plane shear failure of sea ice, ridging for example, that inspired failure mode (iv) in equation (2.5).

The second failure mode happens for *R*_{wind} < 0.1 and features diagonally intersecting lines of high shear strain rate. The number of lines and incident angle between them both appear to increase for decreasing *R*_{wind}; for example the run with *R*_{wind} = − 0.8 has the greatest number of intersecting lines of shear that intersect at an angle of ≈ 90°. Hutchings *et al.* [18], who performed a fuller analysis of linear feature intersection, also found the same relationship of increasing feature intersection angle for decreasing *R*_{wind}. We attribute this failure mode to in-plane shear rupture with the sea ice pack, the inspiration for mode (iii) in equation (2.5). Although mode (ii) can also cause the anisotropic alignment observable in this case, uniaxial compression represented by a negative confinement ratio is little represented within the histogram of confinement in figure 5.

The wind stress confinement *R*_{wind} is not the same as the sea ice internal stress confinement below it. Figure 5 shows an emergent bimodal distribution of *R*_{int}. For *R*_{wind} ≤ 0.0, *R*_{int} peaks near, but strictly less than, *R*_{crit}. For *R*_{wind} > 0.1 and increasing, *R*_{int} has a peak that approaches 1.0. The peak at *R*_{crit} corresponds to observable Coulombic lines of shear and anisotropy, and the peak at *R*_{int} = 1.0 corresponds to parallel lines of shear and isotropy. The experiment with *R*_{wind} = 0.1 (yellow line) displays both peaks, as the domain is split into isotropic regions with parallel deformation features and anisotropic regions with Coulombic deformation features, as seen in figure 6. The bimodal deformation modes of sea ice have been observed in laboratory experiments. Golding *et al.* [6], for example, found a clear bimodal distribution of slip line intersection angle about *R*_{crit} when performing laboratory experiments with imposed stress confinement. A laboratory experiment where the internal stresses could be observed in comparison to imposed external stresses, if possible, would produce a distribution of stress confinement ratio that can be used to contrast with those we show for the EAP and EVP rheologies in this paper.

### (c) The role of alignment and realignment

When the sea ice internal stress confinement is less than the critical confinement ratio *R*_{int} < *R*_{crit}, the sea ice structure tensor becomes aligned (). The internal stress characteristics and link between the confinement of wind stress *R*_{wind} and internal stress *R*_{int} are a function of the directional alignment of the structure tensor (the direction of *A*_{1}). To investigate these relationships we perform the runs described above in §3b but with a pre-aligned sea ice structure tensor and contrast with the EVP rheology. To link the different sea ice structures together, we investigate longer runs with changing wind patterns.

#### (i) Pre-aligned initial conditions

For anisotropically aligned sea ice, pre-aligned to the expected alignment for the given wind field (principal component *A*_{1} aligned with the *x* axis), the parallel slip line failure mode for *R*_{int} > *R*_{crit} shown in figure 6 does not occur. The intersecting Coulombic deformation patterns as shown in figure 6 do occur for *R*_{int} < *R*_{crit} but the sea ice internal stress confinement ratio *R*_{int} is much more closely concentrated around the critical confinement ratio *R*_{crit} for all the runs to produce figure 7. This causes the Coulombic slip lines to only occur for *R*_{wind} ≤ 0.2 for this pre-alignment. For *R*_{int} > *R*_{crit} (for winds with *R*_{wind}≥0.2 in this case) the increased shear stress of aligned sea ice and the converging sea ice combine, resulting in little shear deformation and no identifiable features. For the case of *R*_{int} ≈ *R*_{crit} a new deformation pattern is observed as shown in the last figure in this paper. This failure mode has bands of shear failure that are similarly oriented to the Coulombic slip lines but are much wider and have a blurry appearance. In some cases the bands may be formed of closely packed comb cracks perpendicular to the overall band of shear, though more investigation is needed to confirm this. This failure mode is also observed for the *R*_{wind} = 0.0 run in figure 6 if continued for 2 days, at which point *A*_{1} ≈ 1 and conditions are very similar to the pre-aligned case. Apart from the visual features this aligned failure mode has similar characteristics to the Coulombic failure mode: high shear stress, *R*_{int} ≈ *R*_{crit} and a similar pdf of strain rate magnitude.

When the model is initiated with sea ice aligned perpendicular to the expected alignment for the applied wind field (principal component *A*_{1} aligned with the *y* axis), then there are two failure modes encountered. For *R*_{wind} > 0.0 the sea ice remains anisotropic as in the initial conditions and the increased shear stresses compared to the isotropic case result in no obvious deformation characteristics. For these cases (‘warm’ coloured lines in figure 7*c*), *R*_{int} remains closer to *R*_{crit} than in the isotropic case, with increasing *R*_{wind} resulting in decreased *R*_{int}. For *R*_{wind} ≤ 0.0, the sea ice is able to realign with *R*_{int} going lower than *R*_{crit} within the run time. These runs all go through an ≈ 6 h period of realignment where the structure tensor ** A** and internal stress confinement change rapidly. During this period and afterwards, Coulombic slip lines occur though they move orientation and location much more than in the runs illustrated in figure 6. As the internal stress state changes, divergent weakening occurs, giving the high shear strain rates in the second peak on the pdf in figure 7. For the two failure modes in this experiment there is a bimodal distribution in the histogram for

*R*

_{int}. Those that realign appear to have strictly

*R*

_{int}<

*R*

_{crit}and likewise those that do not realign have

*R*

_{int}>

*R*

_{crit}for all cases.

In comparison to these runs, when using EVP rheology there is no emergent relationship between *R*_{int} and any critical confinement ratio. There is a different peak in the histogram of *R*_{int} in figure 7 for each experiment with varying *R*_{wind}. The pdf of strain rate follows the same curved profile for each EVP experiment.

#### (ii) Changing wind fields: alternating the wind stress confinement *R*_{wind}

The top three panels edged in pink and the line plot within figure 8 illustrate a simulation where the wind forcing alternates. The wind speed increases linearly from still (*u*^{a}, *v*^{a} = 0) conditions at *t* = 0 to a wind field with **U**^{a}_{max} = 15 m s^{−1} and *R*_{wind} = 0.4 at *t* = 6 h with the wind field then held constant until *t* = 12 h. From *t* = 12 to 18 h, the wind pattern is smoothly varied to *R*_{wind} = − 0.4, with the wind speed held constant. From *t* = 24 to 30 h, the wind field is interpolated back to have *R*_{wind} = 0.4 and then held constant until *t* = 36 h, as shown in the line plot in figure 8. The first wind condition causes the isotropic parallel slip lines as found in the *R*_{wind} = 0.4 run in figure 6 with the majority of the ice remaining isotropic and *R*_{int} > *R*_{crit} (blue and red lines in figure 8). As the wind field changes to have *R*_{wind} < 0, the internal stress changes to have *R*_{int} < *R*_{crit} and the sea ice begins to anisotropically align with intersecting slip lines forming as seen in the *R*_{wind} = − 0.4 run in figure 6. As the sea ice becomes anisotropically aligned, *R*_{int} approaches *R*_{crit}, where it remains even after the wind field returns to have *R*_{wind} = 0.4, as the the stress state within the sea ice cannot return the sea ice to isotropy. The deformation field no longer has parallel slip lines and bears a resemblance to the *R*_{int} ≈ *R*_{crit} anisotropic mode in figure 9. The *R*_{wind} = − 0.4 wind field causes heterogeneous irreversible change to the sea ice alignment.

#### (iii) Changing wind fields: rotating the wind direction

The final five panels edged in green in figure 8 show snapshots of a rotating run where the wind field with **U**^{a}_{max} = 15 m s^{−1} and *R*_{wind} = − 0.4 rotates by 45° anticlockwise about the centre of the grid during *t* = 12–18, 24–30, 36–42 and 48–54 h to return to the wind conditions at 6 h, as this wind field has horizontal and vertical symmetry (see green arrows on figure 8). For each forcing arrangement, intersecting Coulombic shear slip lines occur diagonal to, and the sea ice structure becomes aligned parallel to, the principal component of the wind stress gradient at that point in time (see deformation patterns and overlaying quadrilaterals in figure 8). During the transition between the forcing arrangements, large slip lines occur (see animation in the electronic supplementary material) that are more widely spaced and have higher local strain rates. These higher strain rates occur due to the reorientation of the sea ice structure tensor and result in the linear profile in pdf of strain rate magnitude (see green line in the plot for pre-aligned vertical structure in figure 7). In comparison the pdf for the alternating run also shows some linearity, and the rotating run with the EVP rheology has increased high strain rates compared with the runs with constant wind forcing (green line in the bottom row of figure 7).

## 4. Concluding remarks

We have performed idealized numerical experiments using a sea ice model to investigate the link between applied wind and internal sea ice stress conditions and observable deformation characteristics. We have used a sea ice model commonly used within global climate models (the Los Alamos sea ice model, CICE) with minimal adaptation so it runs with idealized initial conditions on a square grid. We run the model at 10 km, 2 km and 500 m spatial resolutions using both the anisotropic EAP and isotropic EVP rheologies. This set-up enables us to illustrate the sea ice dynamical phenomena that can be expected to occur within state of the art high-resolution sea ice climate models and is suitable for testing and comparing all sea ice model rheologies.

We have successfully imposed internal sea ice stress states using winds fields with a constant confinement of the wind stress gradient (*R*_{wind}, see §2c). We have discovered an emergent bimodal relationship between *R*_{wind} and the internal stress confinement *R*_{int}. This result links the EAP rheology and laboratory experiments such as Golding *et al.* [6] where a bimodal relationship between imposed stress confinement and fracture alignment is observed, thus successfully continuing the hierarchy of ice failure over all length scales [7]. With the anisotropic EAP rheology, our numerical experiments with varying *R*_{wind} and initial alignment have revealed five characteristic failure modes as illustrated in figure 9. For *R*_{int} > *R*_{crit}, there are two failure modes. When the sea ice is isotropic, parallel shear slip lines can form. However, when the sea ice is anisotropic, its increased shear strength resists the formation of shear slip lines and the sea ice deforms compressively over large length scales. For *R*_{int} < *R*_{crit}, the sea ice becomes anisotropically aligned due to the mechanics of the rheology in equation (2.5) and there is only one distinctive failure mode, intersecting shear slip lines. A further failure mode occurs for diverging sea ice, where *R*_{wind} = 0.4 with the wind on a bearing from the centre to the edge of the domain. In this case, the sea ice fails in divergence, *R*_{int} < 0 and the sea ice becomes anisotropically aligned.

Previous attempts at characterizing the rheology of sea ice over basin length scales have focused on the distribution of strain rate [11], observable from satellite observation of sea ice drift and an emergent property of sea ice models [12]. Observations show that the pdf of strain magnitude follows a power law when collecting the data from over basin and seasonal length and time scales [9]. In this paper, we investigate the distribution of strain rate for constant and changing stress conditions. We discover that both the EAP and EVP have fewer high strain rate events compared to the expected power law when considering an individual model run with constant stress conditions. However, when using the EAP rheology and either forcing the sea ice structure to rapidly realign or using changing wind conditions, a power law emerges (see figure 7, discussed in §3c). When considering previous studies into the time scaling of sea ice deformation [22–24], this behaviour is expected. This result suggests that for uniform wind conditions, calm periods between Arctic winter storms for example, the observed pdf of strain rate magnitude may not follow a power law. However, due to the scale of the Arctic and that the deformation of sea ice results from internal stresses, including at its boundaries, the likelihood of low stress conditions throughout the Arctic may be low.

In this paper, we have shown clear contrasting features of the EAP and EVP rheologies when using them as part of the currently available CICE model set-up. The EAP rheology consistently gives observable deformation features for all the model resolutions studied, but only for specific situations when using the EVP rheology (§2d). The EAP is also able to produce a power law for the probability of strain rate magnitude and has an emergent critical confinement ratio for the internal sea ice stresses (figure 8).

## Data accessibility

This article has no additional data.

## Authors' contributions

H.D.B.S.H. adapted the CICE model, designed and carried out the experiments, performed the data analysis and drafted the manuscript. D.L.F. provided expert knowledge and experimental direction. M.T. performed preliminary experiments and model development. All authors read and approved the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by EPSRC grant number EP/K032208/1. This paper was produced as part of the NERC project number NE/K011510/1.

## Acknowledgments

The authors would like to thank Jennifer Hutchings for her advice on running CICE in a square domain. We would also like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the 'Mathematics of Sea Ice Phenomena' programme when work on this paper was undertaken.

## Footnotes

One contribution of 15 to a theme issue ‘Modelling of sea-ice phenomena’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4161083.

- Accepted June 27, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.