Introduction
Sea ice is frozen ocean water that covers large areas of the polar oceans, especially in winter, when it can reach 15 × 106 km2 in the Arctic and 20 × 106 km2 in the Antarctic. As sea ice reflects up to 80% of incident solar radiation, insulates the ocean from the atmosphere and contributes to the oceanic energy and salt balance, it plays a significant role in determining polar and global climate. Sea ice is therefore an indispensable part of Global Climate Models (GCMs). Although GCMs use continuum isotropic sea-ice models, there is strong observational evidence of discontinuous anisotropic sea-ice deformation, concentrated not only in leads and ridges at the floe scale, but also in bands extending from tens to thousands of kilometres (Kwok, Reference Kwok2001). Such bands are usually referred to as slip lines (Reference ErlingssonErlingsson, 1991; Reference Overland, McNutt, Salo, Groves and LiOverland and others, 1998) or linear kinematic features (Reference KwokKwok, 2001). the slip lines frequently intersect at acute angles apparently independent of the spatial scale (Reference Walter, Overland and TuretWalter and others, 1995;Reference WeissWeiss, 2001) and delineate regions of approximately diamond-shaped floe aggregates.
The floe aggregates, which we also refer to as blocks, consist of floes frozen together. Reference Hopkins, Frankenstein and ThorndikeHopkins and others (2004) used a discrete-element model to study sea-ice fracture due to different imposed patterns of wind stress, allowing only compressive and tensile failure of inter-floe joints, and found that the formed aggregates had a roughly rectangular shape. This is in contrast to the observed diamond-shaped floe aggregates. It has been argued (e.g. Reference SchulsonSchulson, 2001; Weiss and Schulson, 2009), based upon a strong similarity between sea-ice leads and wing cracks observed in the laboratory, that slip lines form through a mechanism of shear rupture, i.e. failure/crack formation once a sufficient shear stress is reached. Reference Wilchinsky, Feltham and HopkinsWilchinsky and others (2010) incorporated a Mohr–Coulomb shear failure mechanism into the discrete-element model and found that such diamond-shaped aggregates can indeed be formed under uniaxial compression if the tensile strength is an order of magnitude lower than the compressive strength. This result provided an insight into the physics of formation of anisotropic sea-ice structures under a constant wind stress.
In more realistic situations the wind direction is expected to change and the formed cracks will freeze over. Reference SchulsonSchulson (2004) interpreted Landsat-7 images for the Beaufort Sea near the Canadian Arctic Archipelago and identified two separate sliding events along faults oriented in different directions that have presumably been formed by different, consecutive low-confinement compressive wind stresses. In this case the faults formed during the first event seem to have been inactive during the second event and therefore did not affect fault formation and sea-ice deformation during the latter event. Interpreting sea-ice images can provide only limited information on the deformation process, and modelling sea-ice fault reorientation under changing wind directions can provide a more comprehensive picture of the process. Moreover, understanding how changes in wind direction could affect the ice-block shape and orientation is important not only in terms of discrete models, but also for anisotropic continuum models of sea ice that incorporate the diamond-shaped structure of ice blocks (Wilchinsky and Feltham, 2006), as they require an evolution equation for the blocks’ orientation. Here we use the discrete-element model of sea ice developed by Reference HopkinsHopkins (1996) and Reference Hopkins, Frankenstein and ThorndikeHopkins and others (2004) with the incorporated shear rupture mechanism (Wilchinsky and others, 2010) to gain an insight into how changing the wind direction by 90˚ affects the shape and orientation of the sea-ice aggregates formed for different degrees of existing crack refreezing.
In the next section we introduce salient aspects of the discrete-element model incorporating the shear rupture mechanism, and in the following sections we present the simulations of new crack formation for different degrees of refreezing of cracks before changing the wind-stress direction, concentrating in particular upon the range of crack angles and failure modes.
Model Configuration
Sea-ice fracture is studied here using a two-dimensional discrete-element model (Reference HopkinsHopkins, 1996; Reference Hopkins, Frankenstein and ThorndikeHopkins and others, 2004) that has been extended by Reference Wilchinsky, Feltham and HopkinsWilchinsky and others (2010) to account for sea-ice shear failure, in addition to compressive and tensile failure used in the original version. Detailed descriptions of the model can be found in the above papers. the model parameters we use here are given in Table 1.
We consider an ice pack 400 km in size that consists of polygonal floes typically 4 km wide, whose shape is produced by a Voronoi tessellation of the domain in Cartesian coordinates x and y describing zonal and meridional directions, respectively. Initially the floes are connected by frozen joints. the mean thickness of the floes is 3 m and the new ice thickness determining the compressive strength of the joints between the floes is 0.25 m. the floes are considered to be rigid and any deformation at the joints is determined by mutual motion of the adjacent edges of neighbouring floes. the joints can sustain compressive,
Table 1. Model parameters
tensile and shear forces. If they break up, the edges can slide with the ice kinetic sliding friction coefficient taken as 0.32, as was shown by rafting experiments (Reference Hopkins and TuhkuriHopkins and Tuhkuri, 1999), which also falls into the observed range determined by edge-to-edge sliding of several metre-long floes when no stick and slip behaviour occurs (Reference Lishman, Sammonds, Feltham and WilchinskyLishman and others, 2009). Reference Hopkins and ThorndikeHopkins and Thorndike (2006) showed little sensitivity of the model to the kinetic sliding coefficient values. In the model the mutual floe displacement, δ, is split into components δ n and δ s, normal and tangential to the joint, respectively (Fig. 1). the displacement is distributed over the whole floe average length, L, so the depth-integrated normal and shear stresses at a joint are given by
Fig. 1. Adjacent floe interaction.
where E is the compressive elastic modulus, G = E/[2(1 + ν)] is the shear modulus where ν is the Poisson ratio, and h is the mean ice thickness. This also determines k n /k s = 2(1 + ν). In the model, the elastic modulus is 1 GPa and the Poisson ratio is 0.3, taken around the values given by Reference Evans and UntersteinerEvans and Untersteiner (1971). the elastic stress is combined with a viscous stress as a parallelly connected spring and dashpot, included to damp elastic waves and model inelasticity of the ice (Hopkins and others, 2004). the model allows for compressive, tensile and shear failure of the joints with the yield curve presented in Figure 2. the compressive failure stress (strength), σ c, is based on unconfined buckling (Reference Kovacs and SodhiKovacs and Sodhi, 1980) integrated over ice thickness,
Fig. 2. Failure criterion for the joints.
for h min = 0:25 m. the tensile strength of the joints, σ t, is a natural consequence of the assumed cohesion of the joints and linearity of the assumed Coulomb failure curve. It describes collapsing of tensile failure mode I range onto only one point of zero shear stress. Reference Weiss, Schulson and SternWeiss and others (2007) considered stresses recorded during the Surface Heat Budget of the Arctic Ocean (SHEBA) expedition and found only a relatively narrow range of shear stresses at mode I tensile failure at 50 kPa, determining σ t = 12:5 kPa m for h min = 0:25m here. Experiments have shown that the value for sea-ice tensile strength is highly sensitive to temperature and porosity (Reference Richter-Menge and JonesRichter-Menge and Jones, 1993), as well as the specimen width associated with the crack (Reference Dempsey, Adamson and MulmuleDempsey and others, 1999), which are uncertain in our model. In this modelling study we choose the tensile failure stress (strength), σ t, to be ten times smaller than the compressive strength, σ c. While the tensile and compressive strengths are not generally tied to each other, such a ratio between them is crucial for determining formation of diamond-shaped sea-ice aggregates under uniaxial compression (Reference Wilchinsky, Feltham and HopkinsWilchinsky and others, 2010), which corresponds to observations. Moreover this ratio is of the same order as those observed during the SHEBA experiment (Reference Weiss, Schulson and SternWeiss and others, 2007).
The floes in our model are connected by joints modelling the presence of thin ice whose fracture failure mechanism must be prescribed. Generally, laboratory experiments on fracture failure (Reference Schulson and NickolayevSchulson and Nickolayev, 1995; Reference SchulsonSchulson, 2001; Schulson and others, 2006) provide nearly linear yield surfaces of ice in terms of the standard first and second stress invariants (pressure and maximum shear stress). However, our model does not determine depth-integrated stresses acting on the line parallel to the joint. We therefore assume a Coulomb failure law for the joints that describes a similar linear relationship, but in terms of shear and normal tractions at the joint rather than the stress invariants. In this case, shear rupture occurs when the tangential (shear) component, σ s, of the traction at a joint reaches the Coulomb limit,
where σ n is the normal component of the traction at the joint, and μ is the shear rupture coefficient taken to be 0.6, which is within the observed failure envelope slope range. Note that since this failure mechanism describes fractural failure of sea ice and not friction between ice surfaces, μ may differ from the ice friction coefficient (Schulson and others, 2006).
The shear failure criterion (Equation (4)) can be rewritten as
or in terms of displacements as
where δ
c is the compressive failure displacement. As the pure tensile failure is included in Equation (4) for σ
s, for a given shear traction an inter-floe crack forms in a region of a joint where the normal traction reaches either σ
c or 
with only the latter being dependent on the shear traction. These strengths can be reached when the mutual floe displacement increases sufficiently, leading to a proportional increase of the elastic stress at joints (Fig. 1). When a crack extends beyond 95% of the joint length (edge length), the joint fails. After this the floe edges can separate, forming a lead or overlap making a ridge with the vertically integrated ridging stress determined by Reference HopkinsHopkins (1994) as
where l is the length of ice pushed into the ridge.
Sea-ice deformation in the model occurs under the action of a gradient of wind stress, , with stress varying linearly in the x (eastward) or y (northward) directions,
where ρ a is the air density, U wind is the eastward (northward) speed of wind blowing at the western (southern) end of the model pack, while an equal but opposite wind blows at the eastern (northern) end, L d is the domain size, and the tensor e determines the wind-stress direction, with only one nonzero component: either exx = -1 for the initial stage of the experiment where the wind blows zonally along x, or eyy = −1 after it rotates by 90˚ into a meridional direction. Here we use a constant wind speed, U wind = 4ms-1, which corresponds to the wind-stress gradient of 1.15×10-7 Pam-1, as it produces floe aggregates larger than the floe scale. the wind drag applied to a floe is
where A is the floe area and x is its position. the effect of subtracting 
from the position vector is to make the wind drag symmetrical about the centre of the domain. the action of the wind drag force gives rise to a momentum contribution aligned with the force. If the ice pack is homogeneous, the linear variation of the wind force leads to the same mutual displacement of all the floes and to a homogeneous stress distribution. Boundary effects are minimized by calculating the average stress in the central region and applying it to the boundaries every 10 s. All boundary blocks belonging to the same side of the domain then move simultaneously under the imposed stresses. the time-step was 
ensuring that the energy balance between the wind drag work, fracture energy, inelastic dissipation (including ridge building), frictional dissipation and water drag work has an error of <1%.
Simulations
The model starts from an initial configuration of floes at rest. In the first stage of the experiment the wind stress is applied zonally along the x-axis for 3000 s, which is sufficient to break the frozen ice pack into floe aggregates and attain a steady configuration (Fig. 3 for χ = 0). the ice blocks are roughly diamond-shaped, aligned with the compression direction and permeated with internal cracks. Our aim is to study how wind rotation affects the crack pattern. We restrict ourselves to considering a right-angle wind rotation, due to the shape effect of our square domain. If after the first stage the wind is immediately rotated by 90˚, then no further damage occurs. However, we expect the cracks to refreeze at least partially. As discussed by Reference Wilchinsky, Feltham and HopkinsWilchinsky and others (2010), the cracks which are internal to the blocks are likely to be inactive, as the main deformation would occur at the cracks delineating the blocks. Similarly, here we also assume that since they are not active the internal cracks will refreeze completely. the remaining cracks, delineating the block boundaries, are considered to refreeze partially, by a fraction, , of the crack length itself. Generally, since cracks form when the displacement (Equation (6)) has been attained, crack refreezing in the model requires reduction of floe displacement at each cracked joint, which is impossible to achieve without affecting the overlap regions between all adjacent floe edges, which would arbitrarily change the stress field in the ice pack. Therefore, we model refreezing by introducing the cracks formed in the first stage into the corresponding joints of the initial, undeformed ice-pack configuration and reducing the crack length by a fraction, . This configuration is then used as an initial state for the second-stage run where the same wind-stress gradient is applied meridionally (along the y-axis) for 4000 s. Since, as was mentioned earlier, applying the new, rotated wind stress onto the ice aggregates without refreezing does not have any significant effect, we refer to the configuration which is the result of the first-stage run, and which is used as the initial configuration for the second stage as given by χ = 0. Figure 3 compares the ice-cover secondary break-up after refreezing and wind rotation for various degrees of refreezing; it can also be seen that when the block internal cracks are not completely refrozen, the sea-ice cover sustains a high degree of damage. During sea-ice deformation, boundary blocks at the same side remain aligned as the boundary conditions make them move simultaneously. Inspecting the ice-cover failure after the first stage with wind along x, one can see that the internal cracks are preferentially aligned in the zonal direction and their intersection is therefore reduced. Imposing a wind stress rotated at a right angle produces internal cracks in a perpendicular direction, significantly increasing the probability of internal crack intersection. As a result, the ice cover breaks down into floe-scale blocks if the internal cracks are not removed after stage one.
Fig. 3. Ice cover for different refreezing factors, (including = 0), used as an initial configuration for the second-stage runs after partially refreezing the cracks delineating the blocks and completely refreezing the internal cracks. Different colours distinguish between different blocks. the ice cover which resulted from a χ = 0:05 run when no internal cracks were completely refrozen is shown with a high degree of fragmentation. Joints shown in black delineate blocks used in our statistical analysis. Cracks filtered out in our statistical analyses are shown in white. the cracks are filtered out to exclude damage zones by retaining only those cracks separating different ice blocks of which at least one is more than ten times the average floe area. We also remove cracks around the rectangular boundary floe joints. Cracks that surround blocks completely contained within another larger block are also filtered out.
Generally as the refreezing factor, increases, the blocks gradually reorient themselves from the zonal to the meridional direction. For medium values of the block shape is less distinct; however, a clear trend in the crack angle distribution change can be seen from the angular crack distribution in Figure 4. Later we denote angles relative to the compressive stress direction through , which are taken clockwise for consistency with Reference Ashby and HallamAshby and Hallam (1986), as in Figure 5, while their corresponding values in the original coordinate system through are taken anticlockwise in a standard way. In both cases, we consider only first quadrants, θ = π/2 ψ. − ψ At the limiting values of being 0 and 1, the cracks in our simulations (Fig. 4) are distributed bimodally around the critical flaw angle, ψc, determined by wing-crack formation around the corresponding most compressive stress directions (Reference Jaeger and CookJaeger and Cook, 1979; Reference Ashby and HallamAshby and Hallam, 1986; Reference SchulsonSchulson, 2004),
Fig. 4. The histograms of the normalized crack length distribution against the angle around the compression direction, x, for different refreezing factors, , under uniaxial compression. the bin width is 3˚, and the results are averaged over ten different initial Voronoi tessellations. the cracks are filtered as described in Figure 3. the dotted lines show the critical flaw angle, 
(Equation (10)), for compression along x, the dashed line shows θc for compression along y, and the dot-dashed line describes the angle, θr, at which a partially refrozen crack has the same stress intensity as the critical flaw when the wind stress is along y. PDF: probability density function.
Fig. 5. Wing-crack formation under uniaxial compression. the critical flaw is aligned at the angle determined by tan(2 ψc) = 1/μ (Reference Jaeger and CookJaeger and Cook, 1979; Reference Ashby and HallamAshby and Hallam, 1986; Reference SchulsonSchulson, 2004), where for ∅ = 0:392π the maximum stress intensity on an adjacent wing crack is attained.
The maxima, however, are not very prominent, presumably due to the limited number of discrete angles determined by the floe edges, along which cracks can only form. As the refreezing factor increases, the maxima gradually migrate between these two limiting values.
Discussion
Crack angles
In order to gain some insight into the nature of the crack reorientation mechanism, we first note that the critical flaw angle, ψc, is determined by maximizing the stress intensity at the adjacent wing cracks (Reference Jaeger and CookJaeger and Cook, 1979; Reference Ashby and HallamAshby and Hallam, 1986). In the wing-crack theory the critical flaw is subject to frictional sliding (Fig. 5), where the frictional stress is linearly proportional to the normal stress. In our sea-ice model, the situation is physically different: the joint is strained elastically until it reaches the failure criterion (Equation (4)). This initial elastic deformation is, however, analogous to the subfailure regime of frictional sliding, and the stress during failure of the joint is then similar to the stress state during frictional sliding. While our failure criterion (Equation (4)) also includes a tensile strength, σ t, it is constant and so does not affect the maximization procedure for the wing-crack stress intensity, and therefore does not affect the critical flaw angle, ψc. Because of this analogy, here, for simplicity, we consider the wing-crack theory in order to determine how the critical flaw angle would depend on the refreezing factor, χ.
Similarly to Reference Ashby and HallamAshby and Hallam (1986), we consider a two-dimensional situation (Fig. 5), but with only one principal stress, σ < 0. the shear stress along a flaw inclined at an angle, ψ, to the stress direction is
while the normal stress across it is
where we retain the local coordinate notation (x and y) for consistency with Reference Ashby and HallamAshby and Hallam (1986). the shear failure stress given in Equation (4) opposes the sliding motion at the crack, so the effective shear stress there becomes
The mode I stress intensity in a small adjacent wing crack maximized with regard to ∅ is (Reference Ashby and HallamAshby and Hallam, 1986)
where a is the critical crack length. Maximization of the mode I stress intensity with regard to gives the critical flaw angle (Equation (10)), and the maximum stress intensity at this flaw reaches K
I( ψc). In our model, the inter-floe joints transfer the stress through the refrozen parts of the cracks as well as through the broken parts where the ridging force (Equation (7)) acts normal to the cracks. the ridging stress at the beginning of the ridging process (l = 0) is 6.53 kPa, which is <4% of the compressive strength (Equation (3)) and therefore will be ignored here. the refrozen part of the crack constitutes a fraction, χ, of its length, while the remaining part remains stress-free as the ridging stress has been ignored. This leads to a 1/χ increase in the stress acting in the refrozen part of the crack relative to its unbroken state, which leads to the corresponding increase in the stress intensity. Generally, for a single, disconnected flaw, the effect of such a stress increase on the stress intensity would be absorbed by the surrounding unbroken joints. In particular, for an open flaw with zero shear traction, the stress intensity given by Equation (14) remains finite as → 0. However, since in our case the partly refrozen flaws are all connected as they delineate separate ice blocks formed during the initial sea-ice failure, the flaws subject to opening would also be only partly refrozen and therefore have the same stress concentration. Thus we can assume that the stress is increased uniformly, and the stress intensity can be obtained by Equation (13) with σ/χ substituting for σ. As the stress intensity at partially refrozen cracks increases, one would expect rupture to occur primarily in them at the critical angle, ψc, relative to the new wind direction, which we would use for now as our local coordinate system. However, the refrozen cracks are concentrated around the critical flaw angles determined by the initial wind direction, which would now be 
so there are not enough refrozen cracks available to fail at ψc. At the same time, when other angles are considered, the stress intensity reduces. In this case, two effects compete: the stress intensity is reduced when the crack angle deviates from the critical flaw angle, but it is increased as the refreezing factor, χ, decreases. If we find an angle, ψr, which determines the stress intensity at refrozen cracks equal to the stress intensity of new cracks generated at critical flaw angles,
then the refrozen cracks would preferentially not fail at angles > ψr, as the stress intensity there would be lower than the stress intensity at the new critical flaws directed at angle ψc. the tensile strength, σ t, can be cancelled from Equation (15) which then becomes an equation for ψr,
The solution can be found as
which is shown in Figure 6. the angle distribution maxima in the original coordinate system (θ= π/2 = ψ ) associated with the initial wind direction roughly fall between the new critical flaw angle, θc, shown by the dashed line, and the refrozen crack angle with the same stress intensity, θr, shown by the dash-dotted line in Figure 4. As the refreezing factor, χ, increases, ψr → ψc as the stress acting across the refrozen part of the crack decreases until it reaches the magnitude of stress at new cracks (Fig. 6). As ψr approaches ψc, the crack length distribution peaks become situated closer to ψc than ψr. For small the limiting angle for refrozen crack failure, ψr, relative to the most compressive stress direction is close to the initial critical flaw angle, 
(Fig. 6). the stress intensity at the refrozen cracks increases as becomes smaller relative to 
and since ψr and are close, there are enough refrozen cracks available for failure at <
r and the peak is situated roughly halfway between ψr and ψc. As increases, ψr gradually moves away from 
so there are fewer refrozen cracks available for failure, which makes failure of refrozen cracks less probable, and the crack length distribution peaks move away from ψr towards ψc.
Fig. 6. Solution for Equation (17) relative to the highest compression direction.
Failure modes
Given a shear stress at a joint, the failure displacements, 
or δ
c given by Equation (6), must be reached for the joint to fail. We normalize the tensile failure displacement to 
where δ
r can be interpreted as a failure regime parameter: for a tensile strength ten times smaller than the compressive strength it varies from -1 (pure compressive failure) through 0 (pure shear rupture) to 0.1 (pure tensile failure). Figure 7 shows histograms of the normalized crack
Fig. 7. The normalized histogram of the crack length distribution against the normalized tensile failure displacement, δr, for varying refreezing factors, χ, under uniaxial compression. the bin size is 0.1. the dashed line shows the δr values of Equation (18), 
determining the minimum failure energy. the presence of cracks due to purely compressive failure gives rise to a singular large value of the distribution function at the first histogram bin (δr = -1), whose value is shown as f (−1). the cracks are filtered as described in Figure 3.
length distribution during joint failure against δ
r, from which it can be seen that the most probable shear rupture regime is shifting from a compressive mode to a tensile one as the refreezing factor, χ, increases. the position of the peak of the tensile mode (shown by the vertical dashed line) was found by Reference Wilchinsky, Feltham and HopkinsWilchinsky and others (2010) to be determined by minimizing the joint failure energy per unit crack length, 
over all possible failure scenarios, which determines the corresponding failure displacements
This failure regime involves a tensile shear rupture of the joints through opening of the wing cracks. the dominance of shear failure under compression at small refreezing factors, χ, is explained by the fact that failure mainly occurs at refrozen cracks which are aligned at a less acute angle relative to the compressive stress direction. Figure 8 shows the crack angle distribution for cracks failing under compression or tension separately. For both high and low refreezing factors, the cracks that failed under compression are distributed between the critical flaw angle, θc, and the refrozen flaw-limiting angle, θr, given by Equations (10) and (17), which also corresponds to the minima of the crack length distribution under tension. the main difference is that in the χ = 0:05 case the fraction of cracks that failed under tension is ~16%, whereas it is ~41% for χ = 1.
Fig. 8. The histograms of normalized crack length distribution against the angle around direction x for cracks formed under compression and tension separately. the fractions of the cracks failed under tension are 15.75% for χ = 0:05 and 40.51% for χ = 1.
Conclusions
A discrete-element model was used to study how a change in wind direction can alter the orientation of sea-ice floe aggregates mainly formed through shear failure. the initial sea-ice aggregates were formed by imposing a 1.15×10- 7 Pam-1 wind stress gradient for 3000 s in the zonal direction. the formed cracks were then partially refrozen by reducing their length by a particular fraction, the refreezing factor, and the accumulated deformation was disregarded for simplicity. the partially healed sea-ice pack then underwent the same wind-stress gradient for 4000 s, but rotated by 90˚. In order to prevent high damage of the sea-ice pack due to intersection of internal cracks generated in the first and second stages, they were deemed to have refrozen completely. Our calculations show that, as the degree of crack refreezing increases, the crack angle gradually rotates from the critical flaw angle during wing-crack formation relative to the initial wind direction to the critical flaw angle determined by the new wind direction. the position of the maxima in the crack length angular distribution is bounded by the wing-crack theory critical flaw angle, ψc, and an angle, ψr, where the partially refrozen cracks have the same stress intensity. New cracks are less likely to form at angles much lower than ψc, while the partially refrozen cracks are less likely to fail at angles much higher than ψr relative to the new compression direction. Since the number of suitably aligned refrozen cracks decreases as ψr rotates away from the most likely refrozen crack angle as the refreezing factor increases, the secondary cracks become more aligned with ψc rather than ψr. For low refreezing factors, the prevailing mode of shear failure involves compression between the floe joints as the refrozen cracks which are the easiest to fail are at less acute angles relative to the compression direction, and for a 5% refrozen fraction only 16% of cracks fail under tension. For higher refreezing factors, the failure occurs at the wing-crack theory critical flaw angle (30˚) which is more acute relative to the compression direction, so the most likely shear failure occurs under tension at the wing cracks, and for completely refrozen cracks, or failure from the initial, unbroken configuration, 41% of all cracks fail under tension.
