Momentum equations for ice floes

The general form of the linear momentum equation of the ith floe is (e.g. Reference LeppärantaLepparanta, 2005):

where is the mass of the floe, u _i is the velocity of its mass center, t is time, f is the Coriolis parameter, k is a unit vector pointing vertically upward and Si and Vi are the surface and volume of the floe. Analogously, the angular momentum equation can be written as

where w_i denotes the angular velocity of floe i, r the distance from the floe’s center and r _ij a vector pointing from its center to the contact point with floe j.

External forcing acting on the floes

In this work it is assumed that there are no external body forces and that four terms contribute to the total surface forcing, :

The atmospheric and oceanic skin drag, τ _ha,i and τh_w,i, act on the upper and lower surface of the floe, respectively (with surface area of both equal to π_i ² ). The atmospheric body drag, τ_{va, i,} acts on the floe’s edges above the waterline and, in the case of deformed ice, on the slopes of ridges. Here we assume that the vertical area exposed to τ_vai equals πr_ihf _i , where is the floe’s freeboard and ρ_w is water density. Analogously, the oceanic body drag, τνw, i, acts on the floe’s edges below the waterline and on the slopes of keels, and the area exposed to it equals π _ri ( h_i - hfi).

The terms on the right-hand side of Eqn (3) are calculated as

where ρ_a denotes the air density, u _a and u _w are the wind and current velocities, respectively, and C _ha, C _hw, C _va and C _vw are drag coefficients at the horizontal and vertical surfaces of the floe.

Thus, the final form of Eqns (1) and (2) is (see also Reference LeppärantaLepparanta, 2005):

As analyzed by Reference HermanHerman (2011, Reference Herman2012) (for a simplified case without the Coriolis force and without rotational motion of the floes), the response of the floes to the forcing is size-dependent, with larger floes generally having larger equilibrium velocity directed at wider angles to the wind, and reacting more slowly to changes in forcing. In situations with wide, heavy-tailed floe-size distributions, this fact is responsible for a number of important phenomena, including formation of force networks at concentrations below the jamming transition, discussed in this work.

Floe/floe contact mechanics

Each of the normal and tangential forces, , and in Eqns (6) and (7) is a sum of two terms:

The two normal forces are the Hertzian contact force and the damping force . Analogously, the components of the tangential force are the shear force and the damping force . According to Hertzian mechanics, the forces between the interacting objects depend on the distance, 8, between their edges: the forces equal zero for δ ≤ 0 and increase proportionally to δ ≤ 9 for 8 < 0. Details of the Hertzian model in a general context are given by, for example, Reference Brilliantov, Spahn, Hertzsch and PoschelBrilliantov and others (1996), Reference Zhang and MakseZhang and Makse (2005), Reference SchwagerSchwager (2007) and Reference ZhouZhou (2011). In sea-ice studies, (elements of) this model were used, for example, by Reference Hopkins and ThorndikeHopkins and Thorndike (2006) and Reference Fortt and SchulsonFortt and Schulson (2011). From the point of view of this study, the Hertzian contact model has two advantages. Firstly, it is suitable for polydisperse systems and, secondly, the coefficients in the formulations of the terms of Eqn (8) are functions of well-defined material properties of the particles (in this case, sea ice). The material properties of sea ice used in the simulations were obtained by Reference SchulsonSchulson (1999).

The numerical model

Replacement of the hard-disk approach used by Reference HermanHerman (2011, Reference Herman2012) by the more advanced collision model described at the start of this section has some very important (predominantly positive) consequences for the numerical implementation of the model. Whereas hard-disk algorithms progress from one collision to the next, which makes them extremely computationally expensive in situations with high collision rates, the present approach allows the use of a constant time-step, Δt, limited by the characteristic collision duration in the analyzed system. Thus, the model run times are, firstly, predictable and, secondly, to a large degree independent of the details of a particular model configuration.

The numerical model has been built based on the LAMMPS library (Reference PlimptonPlimpton, 1995; http://lammps.sandia. gov/), one of the classical molecular dynamics programs, designed for large systems of particles interacting through a variety of short- or long-range forces. The code is highly parallel, based on neighbor lists and effective spatial decomposition techniques. For the purpose of this work, LAMMPS has been extended to disk-shaped particles in 2-D domains. For all N particles, the Newton equations of motion, Eqns (6) and (7), are solved by means of the velocity-Verlet integrator.

A list of the model parameters used in the simulations is given in Table 1.

Jamming phase transition in pure convergence

Before we proceed to more complex situations, it is useful to analyze the time evolution of a system of floes subjected to pure convergence, without wind. As we will see below, this scenario allows us to define two clearly separated regimes, or states, of sea ice.

Figure 1a shows the evolution of the internal pressure, p, with increasing ice concentration, A. Initially, floe/floe collisions are sporadic, and the corresponding levels of internal stress are close to zero. Increasing A leads to decreasing floe/floe distances and higher collision rates. The relationship between p and A is approximately exponential, until a certain limiting value of A, denoted AJ, is reached, marking the onset of a rapid, hyperexponential increase of pressure. As in other granular systems (Reference Pica Ciamarra, Nicodemi and ConiglioPica Ciamarra and others, 2010), the exact numerical value of AJ depends on the geometrical and material properties of the particles (in this case, exponent a of the FSD, which determines the packing efficiency of floes on the sea surface, and the elastic and damping coefficients that influence the repulsive force between interacting floes), as well as on the path along which the system is driven towards the jammed state. However, in the analyzed cases the value of A _J was found to be relatively stable, between 0.918 and 0.921. It is worth noting that, due to the high packing efficiency of floes with wide size distributions, AJ is higher than typically obtained in monodisperse granular materials.

Fig. 1. Jamming transition in sea ice under pure convergence for dA/dt = 0.04 h–1: (a) evolution of pressure, p(A); (b) rank-order distribution of floe/floe distances, 6, below and above AJ, in the three situations marked with dots in (a).

Qualitative differences between the states below and above AJ manifest themselves in marked differences in properties of the networks of contact forces between the floes. Before the onset of jamming, there is still enough free space that the repulsive forces generated during collisions make it impossible for permanent floe/floe contact to develop. Snapshots of the unjammed system show forces only between pairs of floes undergoing a collision at that particular time instant (Fig. 2a). Accordingly, floe/floe distances tend to be small, but larger than zero (Fig. 1b), the mean contact number η _c ≈ 1 (binary collisions) and the non-rattler fraction f _nr ≈ 0 (Fig. 3). When the system approaches AJ, formation of ‘clusters’ of floes with larger numbers of neighbors (Fig. 2b) leads to increasing p, η _c and f_nr. At AJ, the force network percolates the whole system (Fig. 2c), the distances between a large number of neighboring floes drop to zero (Fig. 1b) and r] _c exceeds 3, which is the minimum value necessary to stabilize the positions of the floes (the number of neighbors of individual floes is roughly proportional to their perimeters). Obviously, under pure convergence η _a and τ remain close to zero (not shown).

Fig. 2. Simulated distribution of instantaneous forces between ice floes under pure convergence, at A equal to (a) 0.908, (b) 0.914 and (c) 0.920, i.e. in situations corresponding to those marked with dots in Figure 1. For each floe i, a line is drawn from its center to the center of floe j=Ci(t). For clarity the plots show only fragments of the modeled system.

Fig. 3. Evolution of the mean contact number, η _c, and the fraction of non-rattler floes, f_nr, with ice concentration, A, during the pure-convergence simulation.

The jamming transition thus separates two very different regimes. Below AJ, in the unjammed state, floe/floe collisions are the dominant source of the internal stress. Moreover, because the floes have some space to move independently of their neighbors, the size-dependent response of individual floes to external forcing plays an important role in the macroscopic behavior of sea ice (Reference HermanHerman, 2011, Reference Herman2012). Above AJ, owing to the existence of a system-size force network, the ice possesses a certain stiffness, or robustness, and can resist certain levels of stress without undergoing irreversible deformation. (Of course, the results become unrealistic at some A 3> AJ, when the load exceeds the ice strength and the floes break, raft on top of each other, etc.)

It is worth noting that even in the jammed state, the fraction of force-bearing floes, f_nr, is still much lower than 1 (Fig. 3). The percolated force network has an openwork structure, with patches of unjammed (usually relatively small) floes surviving among force-bearing floes up to A far above AJ (Fig. 2c). The consequences of this force-network property on point measurements of stress in the ice are briefly discussed in the final section.

Force networks under uniform wind

In order to illustrate the role of the size-dependent response of the ice floes to external forcing, let us analyze a simple case of ice under a uniform wind with ua=10ms1 and a constant ice concentration, A (varied between model runs). The simulations presented below were performed within square domains with periodic boundaries. Tests with other types of boundary conditions and in domains of different sizes showed that the results were not sensitive to the formulation of the boundary conditions.

Snapshots of the system in a quasi-equilibrium state for two selected values of A (Fig. 4) show that wind-forced sea ice has features resembling those present in the transitional phase preceding the jamming transition in the case analyzed above (Fig. 2b). Although the force networks do not percolate the entire model domain, ‘locally jammed’ subregions develop with above-average density of the contact and force networks. The similarity manifests itself in the values of pressure (cf. Figs 5a and 1a) and r] _c (Figs 5b and 3). However, in many ways the wind-forced scenarios are different from the pure-convergence case. Importantly, with decreasing A and increasing role of the size-dependent response of individual floes, the force and contact networks become more localized and more asymmetric. For A<∼0.84, characteristic bands of densely packed floes develop, with higher-than-average p and rj _c, separated by regions of locally lower ice concentration and sporadic, binary collisions (Fig. 4a). Contrary to the convergence case, the system has significant values of the shear stress, r, and anisotropy, ?] _a (Fig. 5), which increases with decreasing A (from almost zero at A = AJ). Thus, wind-forced sea ice possesses some features resembling the shear-jammed states in 2-D granular materials, described by Reference Zhang, Majmudar, Tordesillas and BehringerZhang and others (2010) and Reference Bi, Zhang, Chakraborty and BehringerBi and others (2011).

Fig. 4. Simulated distribution of instantaneous forces between ice floes subject to a constant wind with u_ax = 10ms^_1 and u_ay = 0, at A equal to (a) 0.82 and (b) 0.90. For clarity, the plots show only fragments of the modeled system.

Fig. 5. Statistical measures of the force and contact networks in simulations with constant wind: (a) pressure, p, and shear stress, r; (b) mean contact number, rj _c; and (c) anisotropy measure, η_a.

Force networks in a marginal ice zone

Finally, a set of simulations representing a simplified MIZ was performed, assuming that the floes along the lower boundary of the model domain, representing the fast-ice edge, are stationary (‘frozen’). Periodic boundary conditions were applied to the left and right boundaries. As previously, the wind velocity ua=10ms1 in all simulations; the wind direction relative to the ice edge was kept constant in each individual model run, but varied between simulations. All simulations were performed until a quasi-stationary state was reached.

The modeled across-ice-edge (y-direction) profiles of along-ice-edge ice velocity, u_x(y), are similar to those produced by other models (Reference FelthamFeltham, 2005). At first, ice velocity changes smoothly across the width of the MIZ, but this stage is short-lived. The final velocity profile is characterized by a narrow region of strong velocity gradients, separating stationary ice adjacent to the model boundary from the outer zone, moving with an approximately constant speed, dependent on the alongshore wind velocity (Fig. 6). It should be noted that the motionless inner zone is not limited to the floes arbitrarily declared as ‘frozen’, but extends into the MIZ and comprises a strip of jammed floes.

Fig. 6. Along-ice-edge ice velocity of individual floes, u, _x(y), simulated with three different wind directions (in all cases, ua=10ms and u_a,y denotes the cross-shore component of u _a).

Probably one of the most striking features of the model behavior is the erratic character of the internal stress. Although many general, global characteristics of the modeled MIZ remain stable in time, the area-averaged stress undergoes irregular changes, with uneven cycles of stress build-up and sudden drops (Fig. 7), corresponding to localized jamming events. Unloading of jams can also be seen in Figure 6: during local rearrangements of floes some of them (usually smaller ones) acquire speeds larger than the equilibrium speed corresponding to their position within the MIZ. Importantly, as shown in Figure 7b, although the changes of p and τ are correlated (with a correlation coefficient of 0.68), their instantaneous ratio varies considerably from close to 1 up to ∼12.

Fig. 7. (a) Fragments of simulated time series of pressure, p, and shear stress, r, and (b) their ratio, p/r, in the MIZ under stationary wind forcing (ua=10ms1 and uay=5ms)

Interestingly, the pressure levels, although more than two orders of magnitude higher than in the simple wind-forced cases analyzed above, still remain within the range corresponding to the transitional regime just below AJ (Figs 7 and 1a). Also, independently of the wind direction, nc=1.6 statistical properties of the contact networks remain relatively stable in all analyzed cases. Importantly, they are not sensitive to details of the model configuration. In particular, changes in the material properties of the ice (provided that they remain within reasonable limits) influence the model behavior in the initial phase of the simulations, but their significance is very limited after a dynamic equilibrium is reached under a specified forcing.

Finally, it is worth noticing that the MIZ simulations described above were initialized with floes having random positions and sizes drawn from a single, prescribed FSD (Table 1), whereas in real-world MIZs the average floe size and the shape parameters of the FSD depend on the distance from the ice edge (e.g. Reference Lu, Li, Zhang and DongLu and others, 2008; Reference Toyota, Haas and TamuraToyota and others, 2011). In the idealized cases analyzed here, if the simulations are carried out for a sufficiently long time, the largest/smallest floes tend to migrate towards the inner/outer sections of the MIZ, suggesting the existence of an inherent sorting mechanism transferring largest floes to regions of lower average flow velocities (although I do not suggest that a similar mechanism is responsible for the observed changes of the FSD properties across MIZs).