Discrete-Element bonded-particle Sea Ice model DESIgn , version 1 . 3 . Model description and implementation

Abstract. This paper presents theoretical foundations, numerical implementation and examples of application of the two-dimensional Discrete-Element bonded-particle Sea Ice model – DESIgn. In the model, sea ice is represented as an assemblage of objects of two types: disk-shaped "grains" and semi-elastic bonds connecting them. Grains move on the sea surface under the influence of forces from the atmosphere and the ocean, as well as interactions with surrounding grains through direct contact (Hertzian contact mechanics) and/or through bonds. The model has an experimental option of taking into account quasi-three-dimensional effects related to the space- and time-varying curvature of the sea surface, thus enabling simulation of ice breaking due to stresses resulting from bending moments associated with surface waves. Examples of the model's application to simple sea ice deformation and breaking problems are presented, with an analysis of the influence of the basic model parameters ("microscopic" properties of grains and bonds) on the large-scale response of the modeled material. The model is written as a toolbox suitable for usage with the open-source numerical library LIGGGHTS. The code, together with full technical documentation and example input files, is freely available with this paper and on the Internet.


Introduction
Sea ice cover in polar and subpolar seas is a complex assemblage of ice blocks of various sizes, thicknesses, ages, structures and properties resulting from their genesis, typically consisting of multiple cycles of partial melting, (re)freezing and mechanical deformation resulting from the action of external agents (wind, waves, solar radiation, etc.) and from interactions with surrounding ice.In favorable conditions, the ice blocks may join (freeze) to form larger blocks (ice floes), behaving like semi-rigid bodies, so that the deformation of ice is localized, limited to narrow shear and compression zones.This type of ice cover is characteristic of the compact, central Arctic ice pack.Close to the ice edge, extensive breaking primarily caused by ocean waves produces ice behaving as a polydisperse granular material composed of individual, relatively small floes of various diameters.In all cases, many important aspects of sea ice dynamics are directly related to its discrete, discontinuous nature.Consequently, although some of the large-scale effects of those processes can be parametrized in continuum sea ice models, mechanisms underlying them can be investigated and understood only by means of models that properly take into account the fundamental physics.
Probably the most appreciable way in which the granular nature of sea ice influences other processes, including its own behavior, is through the floe size (mean and the floe-size distribution, FSD).Examples include mechanical weakening of ice after storms resulting from its fragmentation (Holt and Martin, 2001;Asplin et al., 2012Asplin et al., , 2014), manifesting itself for example in the intensity of inertial ice motion (Gimbert et al., 2012;Haller et al., 2014) and influencing large-scale changes in sea ice extent (Kohout et al., 2014); strong contribution of the form drag due to floe edges to the total surface drag coefficient at medium ice concentrations (Steele et al., 1989;Lu et al., 2011;Lüpkes et al., 2012); lateral melting rates dependent on the total floe perimeter within a given area (Steele, 1992); spatial light distribution under ice, with effects analogous to those observed under melt ponds (Frey et al., 2011); and wave propagation and attenuation in the marginal ice zone (MIZ; e.g., Williams et al., 2013a, b).The range of these mechanisms is much wider than this list suggests (see further Sect.2.2) and is only beginning to be appreciated as new observational and computational techniques provide new insights into sea ice physics and dynamics.
The Discrete-Element, bonded-particle Sea Ice model -DESIgn -presented here has been developed as a tool for studying the processes mentioned above at the floe level, with hope that it will help to deepen our understanding of ice dynamics at different scales, and possibly to develop parametrizations of relevant processes for continuum models.
This paper presents the theoretical background, underlying assumptions and equations of the model, its numerical implementation and examples of applications to sea ice problems.The model is an extension of the earlier versions described in Herman (2013a, b).The paper is accompanied by the model code and a full technical documentation, so that it can be freely used and modified by anyone, as described in the last section.A full, very detailed description of all equations used in the model is provided, even if some of them are standard in discrete-element models (DEMs) -for the sake of completeness and in order to make it easier for users inexperienced in DEMs to configure and run their own simulations.The main purpose of the modeling results presented here is the verification of the model's consistency and its new features rather than validation against observational data, which will be presented in further works.
The model is two-dimensional (2-D), but it enables one to take into account some wave-related effects, i.e., stresses resulting from flexural moments acting on sea ice when surface waves are present.It can be applied to a wide range of sea ice types, although it is worth stressing here that the word "granular" in the present context describes macroscopic, largescale ice properties, i.e., the fact that it is composed of individual ice floes.It does not refer in any way to smallerscale material structures at the level of (groups of) ice crystals.Also, in view of specific assumptions underlying the model (e.g., the above-mentioned two-dimensionality), it is not suitable for early stages of sea ice formation, like frazil and grease ice; pancake ice can be regarded as a rough lower limit of the model validity in terms of floe size.On the other hand, although the model was designed for sea ice, it can be applied to other 2-D materials composed of disk-shaped grains.
The paper is structured as follows.The next section contains a short review of previous attempts to account for the granular nature of sea ice in numerical models of its dynamics.Section 3 begins with the general concept and underlying assumptions of the model proposed here, followed by the presentation of the equations.The mechanics of grains and bonds is discussed in Sects.4 and 5, followed in Sect.6 by the presentation of the types of internal forcing implemented in the present model version.The numerical implementation of the model is described in Sect.7. Modeling results illustrating the most important aspects of the bond-related model behavior are presented in Sect.8, followed by a discussion of possible directions of the further model development in Sect.9.
2 Modeling granular effects in sea ice dynamicsa brief review

Sea ice floe-size distribution
From the point of view of processes analyzed here, one of the most important properties of sea ice is the FSD.Since the seminal paper of Rothrock and Thorndike (1984), observations from various parts of the world showed that typical FSDs are very wide and can be well approximated by a power law (e.g., Lytle et al., 1997;Holt and Martin, 2001;Paget et al., 2001;Toyota and Enomoto, 2002;Inoue et al., 2004;Toyota et al., 2006;Lu et al., 2008;Steer et al., 2008;Toyota et al., 2011).However, only a few attempts have been made to create models that would explain the observed variation of the FSD shapes, including the values of the exponent of the distribution and deviations from the power law occurring in certain situations.The validity and range of applicability of the proposed statistical models (Herman, 2010;Toyota et al., 2011;Perovich and Jones, 2014) remain to be assessed.Power laws are produced by a very wide range of models, including very simple models of breaking, making the selection of a proper model and the identification of mechanisms that are important in any particular real-world situation a challenge.

Parametrization of granular effects in continuum models
A number of parametrizations have been developed to improve the performance of continuum sea ice models in situations when granular effects have a significant influence on the large-scale sea ice behavior.One group of parametrizations is related to the so-called collisional rheology, describing stress in fragmented sea ice due to inelastic collisions between floes, relevant especially in the MIZ.In the existing collisional rheology models, stress is calculated from the average collision frequency and momentum transferred during collisions, assuming uniform distribution of the floes on the sea surface and either constant (Shen et al., 1984(Shen et al., , 1986;;Leppäranta et al., 1989) or variable floe sizes (Lu et al., 1989).
Although the models proved useful in reproducing some aspects of sea ice flow in the MIZ (Feltham, 2005), collisional rheology is very rarely taken into account in ice modeling, presumably partly because the rather unrealistic assumptions on which it is based tend to underestimate the contribution of collisional effects to the total stress.Tremblay and Mysak (1997) developed a parametrization of dilatation effects, i.e., flow-induced expansion of granular materials accompanying their slow shear deformation and related to the directional distribution of contacts between neighboring grains.Their model relates the effective angle of friction of sea ice -its macroscopic property -to the microscopic angle of friction characterizing the material.Implemented in a continuum sea ice model of the Arctic, the parametrization permits ice concentrations lower than 1 in the central ice pack, mimicking the lead-opening processes.Steele (1992) proposed a parametrization of the lateralmelting rate in sea ice composed of separate floes of a given mean diameter.A number of studies concentrate on parametrization of floe-related effects on the effective skin and drag coefficients over fragmented sea ice (Steele et al., 1989;Lu et al., 2011;Lüpkes et al., 2012, and references therein), as well as on various aspects of sea ice-wave interactions in the MIZ.For example, Dumont et al. (2011) formulated the first combined strain and stress failure criteria for ice breaking in a flexural mode due to the action of waves.Their results were used by Williams et al. (2013a, b) to formulate a wave-ice interaction model for the MIZ, in which waves break the ice, thus determining the maximum floe size and the FSD, and in turn the FSD is used to estimate the wave attenuation term in the wave energy balance equation.
Following the well-established theory describing the evolution of the ice-thickness distribution in continuum sea ice models, analogous equations for the floe-size distribution (Zhang et al., 2015) and for the joint floe-size and thickness distribution (Horvat and Tziperman, 2015) have been proposed recently, thus providing a general framework within which advection, lateral freezing and melting, as well as ridging and fragmentation processes can be parametrized, making continuum models more suitable for the MIZ.

Discrete-element methods in sea ice modeling
Although continuum models remain a standard tool for simulating sea ice dynamics, especially at large scales, a number of discrete-element models have been developed in recent decades in which sea ice is represented as an assemblage of interacting objects ("particles"); although the models share the same underlying idea, they differ in terms of the shape and properties of their building blocks, details of the contact mechanics formulations, parametrization of physical processes not explicitly accounted for in the model (e.g., ridging), as well as numerical algorithms used to solve the model equations.Some models combine an Eulerian, grid-based approach typical for continuum models with a Lagrangian, particle-based approach used in DEMs -examples include particle-in-cell (PIC) models (Flato, 1993;Huang and Savage, 1998), distributed-mass/discrete-floe (DMDF) models (Rheem et al., 1997;Fujisaki et al., 2007) or smoothedparticle hydrodynamics (SPH) models (Gutfraind andSavage, 1997a, b, 1998;Li et al., 2014).In PIC and DMDF, ice floes are represented by individual particles, advected in a Lagrangian manner based on velocities obtained from momentum equations solved on a fixed grid.By contrast, in SPH models Lagrangian particles represent sets of discrete ice floes.The SPH models with Mohr-Coulomb rheology, implemented within the viscous-plastic approach of Hibler, proved particularly useful for sea ice problems with strong deformation zones and/or complicated geometry (Gutfraind andSavage, 1997a, b, 1998).It is worth stressing that the authors use a DEM model, conceptually very similar to the one proposed here (Savage, 1995;Sayed et al., 1995), to verify their SPH model -they treat DEM as "a very useful tool for the determination of the most appropriate rheology to describe ice as a continuum".Their DEM model is 2-D, based on disk-shaped particles, and takes into account contact forces as well as air and water drag as external forcing to calculate the motion of each particle in the system.
In a series of papers by Hopkins and coworkers (Hopkins and Hibler III, 1991;Hopkins, 1992Hopkins, , 1996;;Hopkins et al., 2004;Hopkins and Thorndike, 2006), the ice pack of the DEM model consists of blocks of thick, multi-year ice and thin ice filling spaces between them, with both types of ice approximated as convex polygons.The model contains a ridging parametrization and, in the newer version that includes elastic joints that may fail under tension or compression, enables one to simulate crack generation and propagation in compact sea ice and produces realistic floe-area distributions in computations of sustained deformation in the Arctic basin.Hopkins (1994) and Hopkins et al. (2004) developed a joint-particle DEM model of ridging due to flexural failure of a floating ice sheet.The same model with an incorporated shear-rupture mechanism was used by Wilchinsky et al. (2010) to study sea ice fragmentation under convergent wind stresses, and by Wilchinsky et al. (2011) in a study on the influence of changing wind direction on the evolution of the crack patterns in sea ice.Hopkins and Tuhkuri (1999) used a 3-D DEM model based on disk-shaped particles with a circular edge to study floe rafting and underturning during compression; the same model has been applied to simulations of pancake ice on waves (Hopkins and Shen, 2001;Sun and Shen, 2012).A DEM model of sea ice has also been used by Xu et al. (2012) to study sea ice breaking and floe formation in the MIZ due to the action of waves.Although their model is based on the same joint-particle principles as the one proposed in this paper, it is used in 2D-V (2-D vertical); i.e., the initial ice floe representing MIZ, composed of a 1 m thick layer of small, bonded particles, is subject to flexural strain induced by the oscillatory motion of prescribed amplitude and frequency.In a very recent work, Rabatel et al. (2015) proposed a DEM sea ice model in which floes are defined as rigid, polygonal objects subject to oceanic and atmospheric drag, and collisions are treated as linear complementarity problems.The model reproduces certain aspects of "granular" sea ice dynamics, e.g., formation of floe clusters due to dissipation of kinetic energy during collisions.Yulmetov et al. (2016) developed a similar DEM based on polygonal elements and applied it to simulate forces acting on an iceberg towed through a broken sea ice field.
Finally, the works by Herman (2011Herman ( , 2012Herman ( , 2013a, b, c) , b, c) present results obtained with DEMs of increasing complexity that have led to the development of the model presented here: from a simple event-driven molecular-dynamics type of model suitable for low-concentration sea ice (Herman, 2011(Herman, , 2012)), to a DEM including a contact mechanics model suitable for studying deformation of a compact ice cover (Herman, 2013a, b, c).The common focus of these works has been on strong polydispersity (i.e., heterogeneity of sizes) of sea ice floes and the role that it plays in sea ice dynamics, including cluster formation of ice floes and stress transmission in ice under shear deformation.
3 Bonded-particle discrete-element sea ice model

General model concept
The sea ice model proposed in this work describes the motion and interactions of two types, or classes, of objects.Following the nomenclature used in similar models of rocks and soils, we will refer to those two types of objects as "grains" or "particles", and "bonds".(In the documentation of the numerical libraries on which the numerical model is based, the terms "atoms" and "bonds" are used independently of the type of those objects; see Sect.7.) Together they form a heterogeneous granular material with large-scale macroscopic behavior depending on the mechanical properties of its individual parts interacting at a "microscopic" scale.In rock mechanics, where bonds represent the cement filling spaces between grains of sedimentary or crystalline rocks, such models have been shown to reproduce a wide range of observed rock behavior, including crack formation and propagation, damage accumulation, and dilational effects (Potyondy and Cundall, 2004;Cho et al., 2007;Asadi et al., 2012;Bahaaddini et al., 2013).Similarly, bonded-particle models are successfully used to study the behavior of cohesive soils, as well as calving and fracturing of glaciers (Åström et al., 2013).
In the context of sea ice, we define the particles as diskshaped sea ice blocks moving within a 2-D space representing the sea surface.The bonds form when the neighboring disks freeze together.In other words, bonds represent new, usually thinner ice filling cracks, leads and other open spaces between thicker ice blocks.Throughout this paper, sets of bonded particles will be called floes -similarly as in Hopkins and Thorndike (2006), who treated ice floes as groups of frozen discrete elements of their model.Depending on the context and the purpose of a particular model configuration, the ice disks may be interpreted as elementary "building blocks" of sea ice floes, or as floes themselves (Fig. 1).In the first case, useful for example in studies of floe formation and evolution of the floe-size distribution (FSD), it is reasonable to base the simulation on a very large number of disks with sizes much smaller than the expected floe size.In the second case, the FSD may be regarded as given, and its influence on the sea ice behavior analyzed in a way similar to that in Herman (2011, 2012, 2013a, b), whose model is a special case of the one presented here.
There are two essentially independent mechanisms of interactions between neighboring particles.The first, based on repulsive and frictional forces between particles, requires that they are in direct contact with each other.The second requires that the particles are connected with an elastic bond.Crucially, whereas forces are transmitted in both cases, bonds are also able to transmit momentum.Another substantial difference between the two interaction types results from the fact that bonds have a certain tensile strength.Consequently, if they are present in a material, it attains a tensile strength at a macroscopic level as well.
From the point of view of sea ice modeling, the bondedparticle approach has a number of important advantages and provides an opportunity to overcome several drawbacks of the existing models.In particular, the model is suitable for arbitrary ice concentration A, from a compact ice cover at A 1 to very loose ice composed of freely drifting floes at A 1. The same equations describe short pairwise collisions and semi-permanent contact between floes.Also, the model is suitable for both a rapid flow regime with relatively large ice velocities and for quasi-static deformation in a compact ice pack.Furthermore, non-cylindrical floes can be modeled as bonded assemblages of disk-shaped elementary particles, so that the model equations are still solved for simple, centrally symmetric objects, without the need to calculate and store their orientation -which is not the case in models based, e.g., on polygonal elements, which require complicated algorithms for calculating grain overlap, force momenta and orientation (see, e.g., Hopkins, 1992Hopkins, , 2004)).
In typical soil and rock applications of bonded-particle models, the simulations are initiated with fully bonded material, and the bonds are allowed to break but not to recover during a simulation (see, e.g., Potyondy and Cundall, 2004;Bahaaddini et al., 2013).In the present version of the sea ice model, new bonds may be created at selected time instances between pairs of grains located closer to each other than a specified distance (see the technical documentation).

Basic definitions and assumptions
Let us consider an ensemble of N disk-shaped ice grains, each with a constant mass density ρ, and with variable thickness h i and radii r i , for i = 1, . .., N .To make a clear distinction between the horizontal and vertical dimensions, we adopt the following notation in a Cartesian coordinate system: ] the translational velocity of its mass center, and ] its angular velocity.As already mentioned, the model is two-dimensional, with some effects related to surface waves included.Precisely, this means that z i ≡ 0 and u z,i ≡ 0, and the distance between disks i and j is calculated as δ ij = x i −x j −(r i +r j ), i.e., within the x 1 x 2 plane.In the absence of waves, the vertical axes of the disks are perpendicular to the horizontal sea surface.If waves are present, z i and u z,i are still equal to zero; the only change is that ω i = 0 and the deviation of the disks' axes from the vertical axis are described by a tilt [θ i , 0] = [θ 1,i , θ 2,i , 0].For details of the wave-related effects, see further Sect.6.3.
The horizontal space S available to the particles may be unlimited or bounded by rigid walls representing shorelines, concrete structures, etc.; in a general case, they may change their position over time.The motion of the disks within S is influenced by (i) the body and surface forces from the ocean and the atmosphere and (ii) the particle-particle and particlewall interaction forces.For the sake of simplicity and conciseness, the particle-wall forces are not included in the equations formulated below (in the model implementation, they are treated in exactly the same way as their particleparticle counterparts, as the contact model used is equally valid for both kinds of contact).
Let C i (t) denote the set of all disks in contact with disk i at a certain time instance t.For each j ∈ C i , the respective pairwise contact force will be denoted with F c,ij .It is convenient to express this force as a sum of two components, normal and tangential to the plane of contact: F c,ij,n and F c,ij,t .The normal component contributes to the translational motion of disks i and j , and the tangential component to their rotation.The respective torque M c,ij equals where r ij is a vector of length r i pointing from the center of disk i to the contact point with disk j .Analogously, let B i (t) denote the set of all disks that are bonded to disk i at time t, and F b,ij and M b,ij the bondtransmitted force and torque, respectively, calculated for all j ∈ B i .Again, F b,ij is a sum of two components, F b,ij,n and F b,ij,t , but the total torque M b,ij not only results from F b,ij,t , but also includes bending and twisting components, as detailed in Sect.5.2.
Finally, the net external force acting on disk i, F e,i , can be obtained from the density of surface and body forces, f s,i and f b,i (measured in N m −2 and N m −3 , respectively): where S i and V i denote the total surface area and volume of the disk.Analogously, the net moment of those forces is where r denotes the distance from the disk center.

Momentum equations for grains
Let us define a unit vector pointing upward, n = [0, 0, 1], and a 3 × 3 projection matrix H from the 3-D space to the x 1 x 2 plane: H 11 = H 22 = 1 and all remaining elements are zero.By definition, the horizontal position x i and tilt θ i are related to u i and ω i , respectively, through The general form of the equation describing the translational motion of the ith grain is where m i = πρh i r 2 i is the mass of the grain and t denotes time.Analogously, the angular-momentum equations can be written in the form of Euler's rotation equations for a rigid body: In the grain-fixed frame of reference, in which these equations are formulated, the inertia tensor is diagonal and its principal values are In accordance with the assumptions formulated in the previous section, the torques due to contact forces have only vertical components and hence do not appear in Eq. ( 7).
Each of the normal and tangential forces F c,ij,n and F c,ij,t is a sum of two terms.The two normal forces are the repulsive force and the viscoelastic damping force; the tangential force is a sum of the shear force and the damping force.The damping components depend on the relative velocity between the interacting particles.At each time step, the amplitude of the tangential force is limited through the Coulomb friction law: where µ is the static yield criterion.
The details of the formulation of the contact forces depend on the selected contact model and on the geometry of the interacting particles, as described in Sect.S1 of the Supplement, which also contains the derivation of the full set of equations used in the present sea ice model.Both alternative formulations available are based on the Hertzian contact mechanics.Additional details of the Hertzian model in a general context can be found, e.g., in Brilliantov et al. (1996), Zhang and Makse (2005), Schwager (2007) and Zhou (2011).In sea ice studies, (elements of) this model were used, e.g., by Hopkins and Thorndike ( 2006), Fortt and Schulson (2011) and Xu et al. (2012).
Notably, in the quasi-3-D version of the model, the contact forces are calculated in the same way as in the 2-D version, i.e., without taking into account the tilt between grains -which amounts to the assumption that the orientations of the axes of symmetry of neighboring grains do not deviate significantly from each other.

Bond properties
In the following, the bonds are identified by the pair of indices of grains that they connect.Each bond, cuboid in shape, is characterized by the following set of properties: thickness h ij ; length b ij ; width 2R ij (measured in the direction perpendicular to the line connecting the centers of grains i and j ; see Fig. 2); Young's modulus E b ; ratio of the normal to shear stiffness λ ns ; tensile strength σ t,max ; compressive strength σ c,max ; shear strength τ max .In the present model formulation, where λ R ∈ (0, 1], λ b ∈ (0, 1] -similarly to E b , σ t,max , σ c,max , τ max -are global parameters, common for all bonds of a given type (the model enables one to define a number of bond types with different properties).When λ b = 1, elastic deformation is calculated as if it were distributed across grains (as, e.g., in the sea ice model of Hopkins et al., 2004); when λ b → 0, it is limited to narrow zones at the grains' boundaries.Note that the distance δ between the grains' edges (Fig. 2) is not included in the calculation of b ij .The normal and shear stiffness, k n,ij and k t,ij , depends on E b and the bond's length: The moments of inertia of the bond are

Bond mechanics
The forces and torques acting on the grains connected with a bond result from the (finite) relative displacement and rotation of those grains; they can be decomposed into axial, tangential, bending and twisting components (Obermayr et al., 2013).Similarly to the history effects in the contact model, the force increment during a time period t can be calculated based on differences between the grains' linear and angular velocities, u ij = u j − u i and ω ij = ω j − ω i .
The components of the force due to the relative displacement are calculated from a linear elastic material law, in which the force is proportional to the displacement, given by t u ij : where S ij = 2R ij h ij is the cross-sectional area of the bond, γ d is a damping coefficient (preventing spurious oscillations of the forces), and u ij,n and u ij,t denote components of u ij normal and tangential to the plane of contact (or, equivalently, parallel and perpendicular to the bond axis).In the 2-D model, both forces act in the horizontal plane and thus F b,ij,t contributes to the grains' rotation around the z axis.
Analogously, twisting and bending moments due to changes in the relative orientation of grains i and j are given by where ω ij,n and ω ij,t denote the normal and tangential components of ω ij , respectively.In total, M b,ij,t in Eqs. ( 7) and ( 8) is given by In the 2-D version of the model, the twisting moment M b,ij,tw ≡ 0 and only the z component of the bending moment M b,ij,bn can be different from zero.Note that in Eqs. ( 14)-( 17), the forces and moments have to be rotated after each time step due to changes in the orientation of the bond.

Breaking criteria
According to the classical beam theory, the shear stress acting on the bond can be calculated as The normal stress reaches its maximum value at the bond peripheries.It has two components, one from the bending moment resulting from the relative rotation of the grains (letters "C" and "T" in Fig. 3), and one from the normal force.
The bending moment produces tension and compression on the opposite sides of the bond, which may be enhanced or reduced by the normal force depending on its sign.Thus, the maximum tensile and compressive normal stress can be written as In the present version, the bond breaks if at least one of the stress components (Eqs.19-21) acting on that bond exceeds the bond strength, i.e., if

External forcing
In terms of the formulation of forces acting on the grains, the model is very flexible and enables one to specify any combination of forces that may be space-and time-varying and depend on the properties of the individual grains (e.g., their mass or size).To make the configuration of the model more (a) convenient, formulae describing the forces most relevant to the motion of sea ice on the sea surface have been implemented in the code and the corresponding forces can be activated easily by means of simple commands described in the User's Guide.These forces include the Coriolis force and the skin and form drag due to the wind and surface current.

The Coriolis force
The Coriolis force acting on the ith grain, F C,i , is given as where the Coriolis parameter f = 2 Z sin φ, Z denotes the angular velocity of the Earth, and the latitude φ can be constant or spatially variable.The net torque due to the Coriolis force: M C,i ≡ 0.

Wind and surface currents
In most real-world situations, the dominating surface forces acting on sea ice floes are the atmospheric and oceanic skin drag, τ ha,i and τ hw,i , and body drag, τ va,i and τ vw,i : where ρ a and ρ w denote the air and water density, the grain.The skin drag acts on the upper and lower surfaces of the grains, respectively, with the surface area of both equal to π r 2 i .The atmospheric and oceanic body drags, τ va,i and τ vw,i , act on the grain's edges above and below the water line, respectively.Here we assume that the vertical area exposed to τ va,i equals π r i h f,i , and the area exposed to τ vw,i equals π r i (h i − h f,i ), where h f,i = h i (ρ w − ρ)/ρ w is the grain's freeboard.In the case of deformed ice, additional drag acting on the slopes of ridges and keels could be taken into account by modifying these expressions.In the present form, the net forces from the atmosphere and the ocean integrated over the respective surface areas become The net torque of F a,i equals zero (a direct consequence of the assumption that the air-ice stress does not depend on the grain's motion).Because of the dependence of F w,i on the disk's velocity relative to the water, the z component of the torque associated with this force, M z,w,i , is different from zero and has a damping effect on the disk's rotation: 6.3 Surface waves

General idea
In the MIZ, as well as in regions with low ice concentrations, sea ice is affected by surface waves (wind waves and, most importantly, swell).Flexural stresses related to the curvature of the sea surface are one of -or presumably thedominant factors leading to floe breakup (e.g., Dumont et al., 2011;Williams et al., 2013a).However, these stresses, being related to forces and torques acting out of the horizontal x 1 x 2 plane, cannot be taken into account in a 2-D model.A question emerges whether -and how -some of the waveinduced effects can be included in the model without introducing full three-dimensionality, i.e., without having to abandon the obvious advantage of calculating the floe-floe distances and solving the equations within a 2-D plane.
The wave-related effects available in the present version of the model are under development and should be treated as a starting point for more advanced models.At present, two wave-related processes have been implemented: forces due to the oscillating surface current, and a net moment of buoyancy forces due to the time-varying sea surface slope and curvature.These two mechanisms tend to be relevant in different conditions.The alternating convergence and divergence associated with oscillatory motion of the sea surface, and the resulting tensile and compressive stress, influence pancake ice formation and dynamics (Shen et al., 2004).It is generally significant for relatively small floes, whereas larger ones are affected by the curvature of the sea surface (Dumont et al., 2011).
Obviously, there are a number of other wave-related effects that are very important but not included in the present model.Most crucially, the wave properties have to be prescribed -they may be spatially and temporarily variable, but are unaffected by the ice, which means that, for example, wave scattering and reflection at the floes' edges cannot be taken into account.Similarly, although the tilt of the grains is calculated, this is not the case for their vertical movement.Whereas some of these additional aspects of wave-ice interactions can be relatively easily implemented in the type of model described here, others, like for example rafting and ridging during compression, would require either a fully 3-D computation (like in Hopkins and Tuhkuri, 1999) or some kind of parametrization -see the last section for more discussion.

Oscillating surface current
Let us suppose that the sea surface elevation ξ = ξ(x, t) is given as a superposition of N wv propagating harmonic deepwater waves: where a n , ω n and φ n denote wave amplitude, frequency and phase, respectively, of the nth component and k n = [k 1,n , k 2,n , 0] is its wavenumber vector.The instantaneous horizontal velocity u wv,n of water particles at the sea surface (z = 0) associated with that component equals After summing the contribution from all elementary waves and averaging over the surface of the ith grain, we obtain cos ϕ n dS and finally The oscillatory current is effective only for small grains and long waves (k n r i → 0) when sin(k n r i )/(k n r i ) → 1.At the opposite extreme, when the grain diameter is large in comparison to the wavelength, the influence of the oscillations cancels out.Equation ( 33) has a computationally convenient form, as it is a product of a time-dependent term evaluated at the position of the grain's center and terms that are functions only of the grain's size.
If the oscillating current is to be taken into account in the model, u wv,i is added to the formula for the current-induced force (Eq.29), u w = ūw +u wv,i , where ūw denotes the slowly varying component of the total current.

Wave-induced horizontal torque
As mentioned above, the presence of waves and the associated space-and time-varying slope of the sea surface induces torque and rotation around the horizontal axes of the ice grains.In the following, an assumption is made that the x components of torque are produced by the unbalanced buoyancy forces acting on a disk if its upper surface is not parallel to the local sea surface, as shown in Fig. 4. It is also assumed for simplicity that exactly half of the disk experiences an excess of buoyancy, the other half an excess of gravity (see also Dumont et al., 2011).As defined in Sect.3.2, θ i = [θ 1,i , θ 2,i , 0] denotes the tilt of the disk; i.e., θ 1,i and θ 2,i are angles between the z axis and the projection of the symmetry axis of the disk onto the x 1 z and x 2 z planes, respectively.Furthermore, let θ s,i = [θ s,1,i , θ s,2,i , 0] denote the mean sea surface slope "under" the disk: tan The unbalanced part of the (vertical) buoyancy force acting on an elementary volume V (r) at the horizontal distance r = r[cos θ, sin θ, 0] from the grain center (dashed area in Fig. 4) is F wv,z,i (r) = ρg h i (r)rdθ dr, where h i (r) = β i • r, β i = tan(θ i − θ s,i ), and ρ equals ρ or ρ w for the emerged or submerged part of the grain.Thus, the total torque, integrated over the disk's volume, is It should be added as part of M e,i to the right-hand side of Eq. ( 7); as the vertical component of M wv,i equals zero, it does not contribute to Eq. ( 8).
For unbonded grains, their angular momentum resulting from Eq. ( 34), and thus variation of their tilt, depends both on their size and the wave characteristics: they decrease with increasing grain radius and wavelength.

The numerical model
The numerical model is based on two libraries designed for effective simulation of large systems of objects interacting through a variety of short-or long-range forces: The code of the DESIgn model has the form of a toolbox that -thanks to the modular, easily extendable nature of the above libraries -can be incorporated into the standard LIGGGHTS program in a straightforward manner, as described in the attached documentation.Importantly, many changes to the model configuration, including specification of additional forcing types, can be made via configuration files, without modification and recompilation of the code.All information on the availability of the code can be found in Sect."Code availability" at the end of this paper.Details regarding the numerical aspects of the model can be found in the User's Guide (available in the "doc" folder in the attached material) and in the documentation of LAMMPS/LIGGGHTS.DESIgn uses the standard methods of the solution of the governing equations implemented in LIGGGHTS.Therefore, only the most important facts regarding the numerical aspects of the model are given here.
The momentum equations are integrated in time using the energy-conserving velocity Verlet solver suitable for finitesize grains and taking into account not only the position and velocity of the center of mass of the particles, but also their angular velocity.Due to numerical stability and so-called energy drift issues, the computations require very small time steps.For sea ice simulations, in which grains typically have diameters of 10 0 -10 2 m, time steps of 10 −4 -10 −3 s are necessary, making the model very expensive computationally.Therefore, realistic model configurations require parallel computations.

Modeling results
All simulations described in Herman (2011Herman ( , 2012Herman ( , 2013a, b, c), b, c), obtained with older, LAMMPS-based versions of the model, can be reproduced with its present version, with proper model settings.Therefore, the results presented here concentrate on the new model features related to the functioning of bonds and to the influence of waves on sea ice.They have been obtained with very simple configurations and can be treated as test cases that verify the model behavior in clearly defined, easily interpretable situations.Both groups of calculations presented below (Sects.8.1 and 8.2) were performed for the so-called reference model setting, as well as for a number of settings differing from the reference ones in terms of one selected parameter, so that the influence of that parameter on the model behavior could be analyzed.The model parameters used in the two reference simulations are listed in Tables 1-3.Simulations similar to those presented below are useful for calibration purposes, as they help to find relationships between the "microscopic" properties of grains and bonds and the desired "macroscopic" properties of sea ice such as its compressive or tensile strength, which in real sea ice varies strongly with changing temperature, ice porosity, microstructure, etc.For a summary of material properties of sea ice, see, e.g., Schulson (1999) and Petrovic (2003).

Sea ice breakup under plane stress
In the first set of simulations, a rectangular sample of compact sea ice (densely packed grains with uniform size distribution, fully connected with their neighbors) is subject to a prescribed uniaxial tensile, uniaxial compressive, or shear strain.The strain rate is obtained by setting to zero the velocity of the grains located at the lower boundary of the domain, and moving the grains located at the upper boundary with a specified velocity until terminal failure (see Herman, 2013c, for a similar model configuration without bonds).In all cases the strain rate increased linearly in time from ε 0 to the maximum value ε e .In the great majority of cases, macroscopic failure of the sample occurred before the end of the simulation.Figure 5 shows examples of damage patterns resulting from compressive, tensile and shear deformation.The evolution of the global maximum normal stress and the shear stress under compressive and shear strain is shown in Figs. 6 and 7, respectively.Analogous plots for tensile strain simulations are shown in Fig. S3 in the Supplement.In all cases, the initial increase in strain results in a fast, approximately linear increase in stress.The rate of the stress accumulation in the material depends on its properties.In particular, it increases with increasing mean bond thickness h m and bond Young's modulus E b (panels a-d in Figs. 6 and 7 and Fig. S3).By contrast, the width of the bond-thickness distribution δ h -that is, the spatial inhomogeneity of the bond thickness -hardly influences the slope of the stress curves.However, it understandably does influence the final damage pattern: for a given value of h m , higher δ h implies a larger number of thin, weak bonds distributed throughout the material, and thus more potential spots where breaking can be initialized.Consequently, a more complex damage pattern develops, with a larger number of fracture zones with more ragged surfaces.This can be seen even in the simplest configurations, like those in Fig. 5a, b: under tension, a uniform material tends to break along an approximately straight line, whereas large δ h in situation in panel b resulted in two competing fractures propagating in opposite directions.Easier initiation of breaking is also responsible for slightly lower macroscopic strength of samples with larger δ h (panels a, b in Figs. 6 and 7 and Fig. S3).
The role of λ ns , defined in Eq. ( 13), is more complex, as its influence on the macroscopic material behavior depends on the type of deformation.(It is worth stressing that λ ns determines the shear stiffness k t ; it has no influence on the normal stiffness k n .)Under uniaxial compression, higher k t stabilizes the bonds and contributes to the overall strength of the material.Under shear strain, varying λ ns manifests itself in very different fracture patterns.Small λ ns , i.e., large k t , produces fracture zones aligned approximately with the direction of motion of the upper boundary of the sample, as in Fig. 5f.The dominating bond breaking mechanism in this case is related to the shear-stress criterion (Eq.19) in zones of high velocity gradient.By contrast, large λ ns and small k t are conducive to breaking related to extensive normal stresses acting on bonds.As a result, tensile fractures develop, penetrating deep into the sample, as for example in Fig. 5e.In all simulations initialized with ε 0 = 0, the initial gradual buildup of stress was related to only isolated breaking of single bonds.This phase was followed by a rapid, avalanchelike increase of bond breaking rates, manifesting itself in ter-minal failure of the material and the associated sudden drop in internal stress as the fraction of broken bonds increased (Figs. 8,9).Under compression, damage has the character of a single, clearly defined event; under shear, damage tends to occur through a large number of events of comparable magnitude.However, the course of stress buildup and the associated deformation of the sample is obviously strongly related to the history of strain, which can be illustrated by varying ε 0 .In the case of compressive strain, putting ε 0 close to ε e amounts to suddenly hitting the modeled sample, which results in very different failure patterns than those described above, with wide damage zones -see Fig. 5d and, for more extreme examples with rapidly increasing strain rates, Fig. S4.
A very important aspect of the model, mentioned in the introduction, is that both types of interactions -those due to bonds and those due to a direct contact between grainsact mutually and contribute to the overall properties of the modeled sea ice.In the simulations discussed above, the majority of grains were in direct contact with their neighbors.Even though, macroscopically, the stress due to pairwise in- teractions tends to be more than an order of magnitude lower than stress transmitted through bonds (e.g., it never exceeded 2-3 % in the reference model run), it in many ways influences the sea ice behavior.In particular, pairwise interactions play a crucial role in the development of the fracture zones, as, obviously, after bond breaking they are the only type of grain-grain interactions present.Thus, friction between grains influences sliding along deformation zones, as well as their width.grain-grain forces are strong in regions subject to compression and along fractures (bottom right in Fig. 11b), whereas bond-transmitted forces are strongest in regions under tension.

Wave-induced sea ice fragmentation
In the second group of simulations, compact, undamaged sea ice is subject to flexural stresses generated by regular deep-water surface waves with a prescribed, spatially uniform amplitude a, period T and length L. All computations are performed until no bond breaking occurs.For the purpose of the basic verification of the model concept described in Sect.6.3, simple one-dimensional (1-D) simulations are performed first (Sect.8.2.1), in which a set of N linearly arranged identical grains bonded to their neighbors are subject to a prescribed propagating wave.This 1-D setting is very close to a 2-D one in which a unidirectional wave propagates through a regular matrix of identical grains, each bonded to its four neighbors.Because in this case the stresses acting on bonds oriented parallel to the wave crests, related to the Poisson effect, are close to zero, the model produces long parallel stripes with a breaking pattern exactly the same as that obtained with a 1-D model version.Therefore, in Sect.8.2.2, devoted to 2-D simulations, only results obtained with an irregular arrangement of grains are analyzed, with a focus on the influence of this factor on the obtained floe sizes and shapes.

One-dimensional simulations
Let us consider a set of N identical grains with radius r and thickness h, arranged regularly (with spacing 2r) along a line parallel to the wave propagation direction.In this case, in the absence of other forcing, the problem reduces to 2N + N b equations, where N b denotes the number of bonds: These equations follow directly from Eqs. ( 5), (7), and ( 17) with an assumption of the damping coefficient γ d = 1.Some indices have been omitted for the sake of simplicity, and it is assumed that I x 1 ,i = I g for all grains, and I x 1 ,i,i+1 = I b and k n,i,i+1 = k n for all bonds.If the initial floe is freely floating on the sea surface, the outermost grains (i = 1 and i = N) have only one neighbor, N b = N − 1 and M 0,1 = M N,N +1 = 0.If one end of a floe is "frozen", e.g., connected to the coast or the landfast ice that remains motionless, we set N b = N and assume ω N+1 = θ N +1 = 0 and M 0,1 = 0.The wave forcing is calculated directly from Eq. ( 34), and the breaking criteria  reduce to a single one: where again h i,i+1 = h b for all bonds.
In the reference simulation summarized in Table 2, the initial floe has length equal to n 0 L, where n 0 = 100 and L is the wavelength, and with the selected grain radius of 3 m, there are roughly 26 grains per wavelength.Initially, no breaking is allowed (hence the infinite bond strength in Table 2) and n 0 is varied to analyze the response of floes of various sizes to the wave forcing.Figures 12 and 13 show the range of variability of the grains' tilt (separated into the rigid and flexural components) and the stress acting on bonds for three selected cases: small (n 0 = 1 and n 0 = 5) freely floating floes (Fig. 12a-d), and large (n 0 = 100) floes both freely floating (Fig. 12e, f)  and connected to landfast ice (Fig. 13a, b).In all cases, three floe lengths, (n 0 − 1 2 )L, n 0 L and (n 0 + 1 2 )L, are considered, because the motion of floes with lengths being an exact multiple of the wavelength tends to be different from the motion of similar floes that do not fulfill this criterion.Whereas in general, as expected, the amplitude of the rigid floe motion increases with decreasing floe size, and for small floes (n 0 < 10) exceeds the amplitude of the flexural motion by at least an order of magnitude (compare continuous and dashed lines in Fig. 12a, c), the rigid motion remains negligible for floes with n 0 ∈ N, and their flexural response is dominated by bending in the middle of the floe.Notably, the smallest floe considered, that of size 0.5L (green lines in Fig. 12a), follows the slope of the sea surface, which in the case con-sidered reaches ±1.15 • , almost exactly, even with a slight inertia-related "overshoot".An analysis of periodograms of the grains' tilts in the cases presented in Figs. 12 and 13 shows that for grains within large floes the power spectra have one dominating frequency -that of the external forcing -whereas additional peaks are present in spectra representing grains within small floes, corresponding to the normal modes of the floes.
Understandably, if only one end of the floe is allowed to move freely, the above-mentioned influence of the exact ratio of the floe size to the wavelength is much less pronounced (Fig. 13a).
For large floes, the amplitude of their flexural motion is largest in the vicinity of their free ends: in the cases tested,  (e, f).In all cases, both ends of the floe could move freely (see text).The distance along the floes is measured relative to the floes' centers.In (a, c, e), continuous lines show the range of variability (i.e., the minimum and maximum values recorded during a simulation) of the flexural components of the floe's motion; analogously, the dashed lines show the range of variability of the rigid motion.In (b, d, f), the continuous and dashed lines show the maximum and the standard deviation of the stress, respectively.Note different ranges of the vertical axes in the panels.Note that in (a), the dashed and continuous red lines cannot be distinguished due to very small values.
it has a roughly exponential profile within the distance of 10-20 wavelengths from the free floe's ends (Figs.12e and  13a).This is extremely important from the point of view of stress acting on the ice (Figs.12f and 13b) and, conse-quently, on the pattern of breaking and the final distribution of floe sizes.Figure 14 shows an example of the breaking history of a large floe like that shown in red in Fig. 13, obtained with bond strength σ t,max = 0.33 × 10 5 Pa, so that  the wave-induced stresses exceed this value close to the floe edge, but not in its inner parts.In this case, breaking starts at the floe boundary, gradually progresses deeper into the ice, and continues until no floes larger than 0.5L remain (inset in Fig. 14a).The pattern of dots in Fig. 14 shows that breaking events are clustered and tend to occur in series spanning a few wavelengths and a small fraction of the wave period.Indeed, almost 40 % of breaking events occur within the distance of one wavelength from the previous one, and the histogram of distances between subsequent breaking events has a clear maximum at ∼ 0.25L (not shown).The final pattern of floe sizes is not perfectly regular (Fig. 14b), but there is one clearly dominating floe size seen in the histogram in Fig. 14a.To the right of this maximum, the shape of the floesize distribution is approximately exponential.This, together with a clearly defined upper limit on floe sizes, suggests that this mechanism alone cannot produce heavy-tailed distributions often observed in the MIZ.

1236
A. Herman: Discrete-element bonded-particle sea ice model Two aspects of these results are worth stressing.First, this breaking pattern, progressing from the ice edge towards inner regions, has been obtained with a spatially constant wave amplitude.In a more realistic configuration, with wave amplitude decreasing with the distance from the ice edge due to attenuation, this effect would be even stronger.Second, with a constant wave amplitude, progressive breaking is possible only if the wave-induced stress acting on the ice far from its edge remains smaller than its strength.In real-world situations, such "tuning" of the wave steepness to the ice strength is presumably rare in stationary settings, but can be expected to occur frequently during periods of wave amplitude increasing in time, e.g., when the onshore wind strengthens or when swell from distant locations arrives at the ice edge with gradually increasing amplitude.Thus, breaking similar to that shown in Fig. 14 may be more frequent than the oversimplified model setting used here suggests.On the other hand, obviously, the great majority of combinations of the model parameters produce either very strong, almost instantaneous fragmentation of the whole initial floe (with the majority of bonds destroyed and a clearly defined maximum floe size) or almost no fragmentation (with only a very small fraction of bonds broken).

Two-dimensional simulations
As already mentioned in the introduction to this section, a 2-D model initialized with a regular, square matrix of identical grains produces floes in the form of long stripes parallel to the wave crests.However, in order for the model to be applicable to more general conditions, e.g., with waves coming from different directions, it is desirable that the model produces realistic results (in terms of both floe sizes and shapes) when it is initialized with randomly distributed grains of different sizes, and thus contains bonds with a range of spatial orientations not aligned with the wave direction.
As can be expected, most aspects of the behavior of the 2-D model are fully analogous to those of the 1-D model.In particular, the dependence of the amplitude of the flexural and rigid motions of the floes on their size, described in the previous section, is very similar in 1-D and 2-D.As previously, the model has been run for a reference run (Table 3) and for a set of configurations with one selected parameter varied.Obviously -as in the 1-D case -most combinations of the coefficients produce either almost no breaking or very intense breaking (like that shown in Fig. 15b).However, within a narrow range of the model parameters between these two "regimes", breaking results in floes with rather unrealistic shapes, with "branches" of connected grains stretching in different directions (some of these features can be seen in Fig. 15a) or even floes with "holes" in the middle, filled with loose grains or smaller floes.This tendency for "branching" floes to survive breaking is related to another rather peculiar aspect of these configurations: a wide range of sizes of the floes.For example, as shown in Fig. 16, for wave pe- It was not possible to obtain a wide distribution of floe sizes in the 1-D model, which supports the notion that this feature is directly linked to irregular floe shapes in the 2-D model.Also, it must be remembered that these features tend to disappear if some irregularity is introduced to the model (e.g., if the forcing is specified as a superposition of more than one elementary wave with different directions and phases), so that it may seem irrelevant in more realistic model settings than those analyzed here, but nevertheless it signals an undesired feature of the model.The tendency to produce  non-convex floes decreases with increasing shear stiffness and decreasing shear strength of bonds (relative to their normal stiffness and compressive strength, respectively), but it remains to be investigated whether adjusting these parameters accordingly does not produce any other negative effects.In any case, the general conclusion is that the orientations of bonds in the model do determine the resulting geometry of the floes.In particular, although the floes tend to be elongated in the direction of the wave crests (Fig. 15), it is not possible to obtain long stripes even with perfectly unidirectional wave forcing.

Discussion and further perspectives
The modeling results presented in this paper have been limited to very simple configurations, with the goal of testing the basic model features and analyzing the influence of the model parameters on its behavior.Further work is necessary to verify the model -its underlying assumptions, numerical algorithms, etc. -in a wider range of configurations.Importantly, the modular structures of LIGGGHTS in general and of the sea ice toolbox in particular make it relatively easy to modify and/or replace some parts should better or alternative solutions become available.In the same way, the model may be extended with new features, including time variability of the properties of grains and bonds (e.g., changes in their thickness due to thermodynamic processes), other contact models or bond breaking criteria, or a model of wavesea ice interactions in which the wave characteristics are not prescribed but modified based on the properties of the ice cover.As noted earlier in Sect.6.3.1, the wave-induced flexural breaking of ice is just one of a number of breaking mechanisms that may be of different importance in different conditions.In particular, it is desirable to include at least some of those mechanisms in DESIgn, e. account buoyancy forces acting on the ice on the vertically moving sea surface.This is planned for the next version of the model that will be available soon.It is also worth mentioning that the model can be applied to simulate the motion of icebergs (solitary or embedded in sea ice) and any other objects floating on the sea surface, provided suitable forcing is formulated.On the other hand, the basic assumption regarding the two-dimensional nature of contact forces between grains makes it difficult to implement in the model processes that are inherently three-dimensional, like for example rafting and ridging.Some parametrizations of these processes will be necessary to make the model applicable to situations when ridging cannot be disregarded.
The fact that DESIgn is a toolbox of LIGGGHTS offers a possibility to relatively easily combine this model with a wide range of functionalities offered by LIGGGHTS.In particular, the mesh-free SPH method available in the opensource LIGGGHTS version could be used to simulate sea ice in a manner similar to Gutfraind andSavage (1997a, b, 1998).Also, coupling the DEM sea ice model with a continuum model of the ocean and/or atmosphere is possible thanks to the CFD-DEM engine (Goniva et al., 2012), where CFD stands for computational fluid dynamics (see also Zhao and Shan, 2013).The tool provides solvers for both unresolved and resolved coupling.In the first method, the grains are taken into account in the Navier-Stokes equations by means of a volume fraction they occupy in each computational cell.The second method, suitable for relatively large grains, resolves the fluid motion around each grain.The CFD-DEM could be used, e.g., to solve the equations of the ocean mixed layer under fragmented ice, or to simulate the wave-ice interactions as a fully coupled problem.The above-mentioned resolved coupling would be particularly useful for that purpose, allowing for simulation of 3-D water motion around floating ice floes.Among more technical issues, meshes with complex geometry can be used in simulations with coastlines and other boundaries.
Among the challenges related to the usage of DESIgn (and other, similar models) in realistic sea ice problems, validation with observational data undoubtedly belongs to the most urgent ones.Although substantial progress has been made in recent years in terms of availability of high-resolution remote-sensing data, the temporal and spatial resolution of those data is still too low to capture some rather subtle effects described in Herman (2011Herman ( , 2012Herman ( , 2013a, b, c) , b, c) and in this paper.Moreover, they only provide information on the position of objects within images, from which average velocities between subsequent snapshots can be calculated with the help of motion-tracking algorithms.Validation at shorter timescales and in terms of other variables requires carefully planned and usually expensive in situ or laboratory measurements, e.g., with accelerometers and stress sensors placed on (a number of) ice floes within a given area, combined with meteorological and oceanographic observations.Because of very high costs of such undertakings, numerical models may be particularly useful in planning of experiments, suggesting processes worth investigating and methods suitable for their analysis.

Figure 1 .
Figure 1.Two approaches to the usage of grains and bonds in DE-SIgn: without bonds, when ice floes are identical to grains and the FSD is prescribed (a); with bonds, when ice floes are assemblages of bonded grains and the FSD is obtained as a solution (b).

Figure 2 .
Figure2.Geometry of two grains, i and j , connected with a semielastic bond.

Figure 3 .
Figure 3. Forces and torques acting on an elastic bond connecting grains i and j : top view (a) and side view (b).Letters C and T denote points where maximum compressive and tensile stress, respectively, occurs due to the bending moment caused by ω ij,z in (a) and ω ij,t in (b).Gray curved arrow in (b) denotes the twisting moment associated with ω ij,n .
u a and u w are the wind and current velocities, and C ha , C hw , C va , and C vw -drag coefficients at the respective surfaces of www.geosci-model-dev.net/9

Figure 5 .
Figure 5. Example damage patterns obtained in simulations of an initially compact sample under uniaxial tensile (a, b), uniaxial compressive (c, d), and shear (e, f) strain.Thick gray lines show the bonds between grains.Model parameters that differed from the reference run are given with each panel.See Figs.S3-S5 for more images.

Figure 6 .
Figure 6.Amplitude of the maximum normal(a, c, e) and shear (b, d, f) stress due to bonded interactions in under uniaxial compressive strain, with variable model parameters.

Figure 7 .
Figure 7.As in Fig. 6 but in simulations under shear strain.

Figure 8 .
Figure 8. Temporal evolution of the fraction of broken bonds and the rate of bond breaking (a); relationships between the fraction of broken bonds and the global normal and shear stress in simulations under uniaxial compression (reference model settings).

Figure 9 .
Figure 9.As in Fig. 8 but in simulation under shear strain.

1233Figure 10 .
Figure 10.Instantaneous normal stress (color scale; in Pa) acting on individual grains in simulation under uniaxial compressive strain shortly after terminal failure of the material (t = 600 s; see Figs. 6 and 8).Panel (b) shows a fragment within the rectangle marked in (a), with bonds between grains illustrated with thick blue lines.Reference model settings.

Figure 11 .
Figure 11.Instantaneous shear stress (color scale; in Pa) acting on individual grains due to bond interactions (a) and pairwise interactions (b) in simulation under shear strain at t = 2000 s (see Figs. 7 and 9).Reference model settings.Note different color scales in (a) and (b).

Figure 12 .
Figure 12. Results of 1-D simulations of sea ice response to waves (without breaking): tilt of the grains (a, c, e) and stress acting on bonds (b, d, f), for floe lengths equal to (n 0 − 1 2 )L (green), n 0 L (red) and (n 0 + 1 2 )L (blue); n 0 = 1 in (a, b), n 0 = 5 in (c, d), and n 0 = 100 in(e, f).In all cases, both ends of the floe could move freely (see text).The distance along the floes is measured relative to the floes' centers.In (a, c, e), continuous lines show the range of variability (i.e., the minimum and maximum values recorded during a simulation) of the flexural components of the floe's motion; analogously, the dashed lines show the range of variability of the rigid motion.In (b, d, f), the continuous and dashed lines show the maximum and the standard deviation of the stress, respectively.Note different ranges of the vertical axes in the panels.Note that in (a), the dashed and continuous red lines cannot be distinguished due to very small values.

Figure 13 .
Figure13.As in Fig.12e, f, i.e., for n 0 = 100, but for a floe with the right end "frozen" (see text).The distance along the floes is measured relative to the floes' right edges.

Figure 14 .
Figure 14.Wave-induced breaking of an ice floe with initial size 100L, "frozen" at its right end (x = 0).In the main panel in (a), individual breaking events are shown with dots in function of their position and time.The inset in (a) shows the histogram of the floe sizes (l/L) after breaking.In (b), a fragment of the initial floe, 20L long, is shown, with breaking positions marked with vertical lines.

1237Figure 15 .
Figure 15.Fragments of the model domain at the end of the reference simulation, with T = 25 s (a; wavelength L = 976 m), and a simulation with T = 23 s (b; wavelength L = 826 m), showing the pattern of floes.The wave propagation direction is from top to bottom.Note that the whole area of both panels is covered with ice; i.e., the large white spaces are large floes, not open water.

Figure 16 .
Figure 16.Rank-order statistics of floe sizes (a: number of grains in a floe; b: floe surface area) obtained in simulations with different wave periods T .The dashed line in (b) marks the area of the largest individual grain in the ensemble.Results of the reference run are shown with crosses.

Table 1 .
Physical and numerical model parameters used in the reference simulations in Sect.8.1.

Table 2 .
Physical and numerical model parameters used in the reference simulations in Sect.8.2.1.

Table 3 .
Physical and numerical model parameters used in the reference simulations in Sect.8.2.2.