A particle based simulation model for glacier dynamics

Introduction Conclusions References


Introduction
The formation and propagation of fractures underpins a wide range of important glaciological processes including crevasse formation, iceberg calving and rheological weakening of ice in shear margins and icefalls. Numerical simulations of glaciers, however, almost exclusively employ continuum methods, which treat ice as a continuous medium with uniform or smoothly varying properties. The difficulty of dealing with discontinuities in continuum models means that the effects of fracturing are routinely represented by simple parameters, such as depth of fracture penetration (Benn et al., 2007a, b;Nick et al., 2010), bulk "damage" Borstad et al., 2012), or ice softness (Vieli et al., 2006). While useful for many purposes, these approaches impose major limitations on the kinds of glacier behaviour that can be represented in prognostic models.
Iceberg calving and the fracture of ice remain intensively active topics of interest, which is a testament both to the difficulty of the work, and the long term monitoring required to quantify its statistical nature (Weertman, 1974(Weertman, , 1980Schulson, 2001). Calving constitutes up to 40-50 % of mass loss on marine terminating ice fronts (Burgess et al., 2005;Dowdeswell et al., 2008;Walter et al., 2010;Thomas et al., 2004) in the regions where it has been documented. Marine terminating glaciers and ice shelves account for almost all mass loss through calving in the case of Antarctica and about 50 % for Greenland (Rignot et al., 2011;Jacob et al., 2012). Calving glaciers are very variable, but two end J. A. Åström et al.: A particle based simulation model for glacier dynamics member types can be recognised: (i) glaciers with grounded termini, and (ii) floating ice shelves that are constrained only at their lateral margins. The two scenarios produce radically different types of calving: (i) small ice blocks that fall off the calving cliff in typically warm tidewater glacier settings, and (ii) large flat-topped bergs that can be tens of kilometres across from the colder ice shelves that fringe the polar ice sheets. At present ice sheet models do not typically incorporate calving as a function of atmospheric and oceanic forcing. Indeed, no formulation for calving has yet been agreed as suitable for models, though several have been proposed (Benn et al., 2007a;Nick et al., 2010;Bassis, 2011;Levermann et al., 2012), and indeed different ones may be suitable for different applications such as large scale models (e.g. Levermann et al., 2012) or basin-scale studies (Benn et al., 2007a) with floating ice tongues (e.g. Nick et al., 2010). Benn et al. (2007a) proposed a physically based model with the position of the calving front depending on the depth of penetration of surface crevasses, which in turn depends on the longitudinal strain rate. A modification suggested was to increase crevasse depth by the filling of crevasses by surface melt water, which is common occurrence in summer even on ice sheets, and certainly typical of many marine terminating smaller glaciers. Nick et al. (2010) introduced a further modification by including basal crevasses with a calving criterion when surface crevasses reach basal crevasses. Basal crevasses can penetrate much farther than surface air-filled crevasses, hence potentially triggering calving at a greater distance from the terminus. The existence of huge tabular icebergs, originating from floating ice shelves, provides ample motivation for incorporating this effect. An upper bound for the height of calving ice cliffs was suggested by Bassis and Walker (2012).
There is a long history of using particle models to simulate geophysical phenomena (Cundall and Strack, 1979;Jing, 1992;Gethin et al., 2001;Potyondy and Cundall, 2004), but usually the material behaviour studied with these models have been restricted to elastic and brittle properties. In this paper, we introduce a new, particle-based method for modelling ice flow, which allows elastic, viscous and brittle behaviour to be represented within a single framework. Although based on simple rules, a very wide range of glaciological phenomena emerge from the model, allowing detailed investigation of processes that are difficult or impossible to represent in continuum models. We first describe the model, then illustrate some of its potential applications, including calving events, the effects of variable basal friction, and threshold behaviour in sliding rates ("surging"). Further, we present an ice-deformation calculation comparing the results obtained with the particle model to the ones obtained with the FEM code Elmer/Ice (http://elmerice.elmerfem.org) for benchmarking. Fig. 1. The particles connected with a beam can (a) stretch when a force F is applied and (b) bend when a torque T is applied. Particles that overlap, i.e. come into contact (c) will experience repulsive forces. The amount of stretching and bending required for beam breaking is highly exaggerated as is the amount of particle overlap in the simulations.

Model
In our simulation model, a large ice-body is divided into discrete particles. A detailed theoretical description of the model is given in Appendix A. The typical diameter of the particles is in the present simulations of the order of 1 m. Initially, at the start of a simulation, particles are densely packed (close-packed) or deposited to form a large body, and the particles are assumed to be frozen together. The particles can either be arranged more or less randomly as in an amorphous solid, or as in a regular lattice. The frozen contacts between the particles are modelled by elastic massless beams which can break if stretched, sheared or bent beyond elastic threshold limits. In such a case the beam vanishes (see Fig. 1). Choosing a proper fracture criterion is far from trivial. A general elliptical criterion (Zhang and Eckert, 2005) that includes the 'classical' fracture criterion of Tresca, von Mises, Mohr-Coulomb and the maximum normal, i.e. hydrostatic pressure, stress criterion is rather useful. This criterion states that a material under tension fails at locations where in which σ I is the normal stress, or pressure determined by the first invariant of the stress tensor and σ I I is shear stress determined by the second invariant. σ 0 and τ 0 are material dependent constants. Instead of fracture stress thresholds, it is also possible to use thresholds for elastic strain. It is trivial to change between stress and strain via the relation σ = K , where σ is the stress tensor, K is stiffness tensor and is strain tensor. Also, energy based criteria can easily be formed. Then fracture takes place if the total elastic energy of a beam grows beyond a threshold. It is still an open question which criteria are best for glacier simulations. No mass is lost when beams break. If particles detach, they are able to flow past each other and thereby collide with other particles. The shape deformations of the particles are not calculated exactly. The contact forces in a collision are calculated as a function of overlapping of particles. The collisions are inelastic and kinetic energy is lost in every collision. This means that once all contacts are broken in an ice-block under local compression, it will display granular flow. In parts of an ice-block with extant connections the ice will continue to behave as an elastic solid. Under tension the ice is able to fracture via crack formation and propagation. The model should thus contain the necessary ingredients for simulating a visco-elastic material, like ice, that fractures, at least on a qualitative level. The equations of motion may vary slightly with the exact implementation for the interaction of the particles (cf. Appendix A), but can typically be written in the form: where M is the diagonal mass-matrix containing the masses (i.e. volume times density Vρ) and moments of inertia (ρ r 2 dV ) of the particles. Vectors r i andṙ i denote the position and velocity of particle i and r ij andṙ ij are the corresponding relative vectors for particles i and j . The diagonal damping matrix C contains damping coefficients for drag, i.e. drag force = (1/2)ρv 2 c D S, where v is velocity, S is cross-sectional area of the object to which the drag is applied, and c D is the Reynolds number dependent drag coefficient. The other damping matrix, C , contains coefficients for the inelastic collisions. The parameter γ ij is zero for particles not in contact and unity for particles in contact, and γ ij is unity for connected particles and zero otherwise. The stiffness matrix is denoted by K and F i is the sum of other forces acting on particle i. These forces may include gravitation (gρV ), buoyancy (gδρV ), where δρ is density difference between ice and fluid, atmospheric and hydrodynamic/static forces, etc. In order to simulate ice we use the Young's modulus E = 10 9 N m −2 , density ρ = 910 kg m −3 and fracture strain c = (1-5) ×10 −4 (Schulson, 1999). The damping coefficient for collisions is C = 10 5 N s m −1 . C represents the drag force on ice falling into water in a calving event. A typical value for this parameter is 10 3 kg s −1 . If the contacts between the particles are broken, the material consisting of only the particles behaves as a nearly incompressible fluid. If the contacts are not broken, the material consisting of particles and beams, deforms elastically under small deformations. In granular flow, the viscosity depends on various factors such as the packing density and the cohesion between particles in contact. These are two important parameters that affect the diffusion and the momentum transfer between colliding particles, which is the microscopic origin of viscosity in any material, including polycrystalline ice. One of the many contributions to viscosity of ice comes from grain-boundary sliding (Johari et al., 1995). In general, diffusion increases with temperature, which means more "liquid-like" for higher temperatures and more "solid-like" for lower temperatures. For ice, it therefore seems reasonable to model the viscous cohesion as a "melting-refreezing probability". This means that once the elastic beam that models the frozen contact between adjacent particles is stretched or bent, the probability for that beam to break becomes non-zero. In contrast, if the tensile strain of the beam reaches the fracture strain the contact always breaks. Also, when particles without a connecting beam are close to each other, a beam can be created with a small probability allowing the material to "refreeze". When combined, the two effects allow the material to undergo constant liquid-like deformation (or regelation in the case of ice bodies) as well as fracture. Notice that this method differs from fluid-like particle models such as smoothed particle hydrodynamics (Monghan, 1992(Monghan, , 2005. Furthermore, the melting-freezing probability can be adjusted to produce stress-dependent viscous flow obeying Glen's law for flow rate, D = A(T )σ n−1 e t D , where A(T ) is a temperature dependent Arrhenius factor, σ e is the second invariant of the deviatoric stress tensor, D is the strain-rate tensor, t D is the deviatoric stress tensor, and n ≈ 3. Correspondingly, dynamic viscosity is µ ≈ 1 Aσ 2 . Details of this derivation are given in Appendix A. The probabilities can be adjusted such that the desired viscosity can be acquired, that is the pre-factor, A, can be changed by adjusting the probabilities. Usefully this allows the temperature dependence of A in Glen's flow law to be incorporated in the model. Computational problems arise, however, from the fact that the time step length is limited by the rapid timescale of the brittle failure events to approximately 10 −4 s, while the relevant time scale for viscous flow of ice is much longer. To cover both relevant time scales in a single simulation is impractical. It is however possible to use lower viscosities and thereby higher strain rates and re-scale the simulation time to match ice behaviour as long as the viscous flow timescale remains slow compared with that for fracture events (Riikilä et al., 2013). This approach is somewhat similar to the plasticity model used by Timar et al. (2010). Another, simpler, approach to imitate viscous behaviour is to use a weak and short-range attraction force for particles that are close to being in contact, similar to cohesion models of wet granular materials (Huang et al., 2005). Both approaches seem to give fairly realistic results. The former approach is benchmarked against a continuum model below.
It is also possible to introduce friction between the particles. This would add another, tangential, interaction potential between grains in contact, although that would be in addition to existing interactions modelled which also include tangential forces. Moreover, as long as particles are connected to form larger blocks rather than just being individual particles, (which is almost always the case in the simulations), granular friction appears as a natural consequence of the surface roughness of the blocks. Because of these effects, and for www.the-cryosphere.net/7/1591/2013/ The Cryosphere, 7, 1591-1602, 2013 simplicity, we have not introduced an explicit friction force between the particles. A flow chart showing a representative algorithm of a simulation is presented in Fig. 2.

Figure 3 shows snapshots of a calving 30 m×30 m ice-block.
In this simulation the material is purely elastic, i.e. no viscous component is present. The Young's modulus and fracture strain are set rather low, 10 8 Pa and 10 −4 , respectively. This example is not intended to exactly mimic any particular real glacier. There is also a significant fraction of the beams missing to mimic damaged ice and there is an artificial crevasse at the top left of the ice block to initialise calving on the left side. The model block rests on a soft substrate that hinders the block from sliding. This "muddy sea floor" is modelled as a linear spring force prohibiting the block from sinking too deep and from moving sideways once it is stuck in the mud. The block is immersed in 20 m of water. The water is modelled here only as a buoyancy force. The simulation times of the snapshots are indicated. The time resolution in the simulation is 10 −4 s. In this case the initial configuration is unstable and as soon as the simulation starts at time t = 0, cracks appear and the ice-block calves. The duration of this single calving event is 10-20 s, which is realistic in comparison with similar events in nature. Figure 4 shows the fragment size distribution, n(s), from the simulations displayed in Fig. 3. n(s)ds is the number of fragments found in the size interval [s − ds/2, s + ds/2], and s is the number of particles in a fragment. This is, obviously, proportional to the volume of the fragment, i.e. roughly the same amount in m 3 . Note that the distribution in Fig. 4 is a relative probability distribution, which means that values of n(s) below unity can occur. Larger bin sizes have been used for larger fragments to avoid empty bins and the relative probabilities are calculated as fragments in one bin over the bin width. Results for three different simulations are shown and we distinguish between the size distributions early during the calving event and late during this process when the fragments have come to rest floating on the water. The fragment size distributions are compared with that of the universal crack-branching-merging model for fragmentation of brittle materials. This distribution is given by (Åström et al., 2004;Åström, 2006;Kekäläinen et al., 2007): where s is again the fragment size or volume, D is the dimension (D = 2 for the simulations here and D = 3 for real glaciers), s 0 is a parameter which depends on, e.g. the material and the fracture energy. The result shown in Fig. 3 is consistent with field data by Crocker (1993) and Savage et al. (2000) from Bonavista Bay on Newfoundland, and with the data of Dowdeswell and Forsberg (1992) from Kongsfjorden on Svalbard. The consistency shown means that the field data presented in these papers have, approximately, the same shape for the fragment size distribution as Eq. (3).
Next we turn to verifying the viscous behaviour without fracture. In order to verify the flow behaviour of the model a comparison with Elmer/Ice (Gagliardini and Zwinger, 2008;Zwinger and Moore, 2009) was made. In both simulations an ice block was placed on a flat surface with little/no friction and gravity as the only driving force. The result is displayed in Fig. 5. The deformation of the ice blocks are obviously quite similar. The main difference appears during the early times of the simulations. This is probably due to a partial jamming effect of the granular material. The particle model parameters are set such that the resulting viscosity is 10 5 times lower than in the Elmer model, leading to 10 5 times faster strain rate, i.e. strain rate is proportional to the inverse of viscosity. In a set of similar tests with varying shear rates, we also verified that the particle model can reproduce Glen's flow law (cf. Appendix A).
In order to further investigate the behaviour of our model, we simulated the dynamics of an ice-block on a slope. We chose a block of size 200 m×50 m on a 18 • slope. Again, the viscosity was roughly 10 5 lower than realistic values for ice, i.e. an Arrhenius factor ≈ 5×10 −19 s −1 Pa −3 . We would thus expect that the strain rates will be roughly 10 5 higher than realistic rates for ice deformation. It is thus, in a crude sense, possible to re-scale the simulation time, which is calculated in seconds, to approximately days (24 h = 0.864 × 10 5 s).
In order to mimic natural variation in bed "stickyness", we also introduced the possibility to locally switch between a high friction, i.e. a no-slip condition, and zero friction for the contact between the substrate and the ice-block. We anchored the base of the block, indicated by a red line in Figs. 6 and 7, by high friction against the substrate. We also included a pressure, corresponding to the over-burden pressure on the upper vertical edge of the ice-block (Figs. 6 and 7). This simulates roughly the pressure induced by a slab of ice with same thickness upstream. Finally, we investigate how the ice-block slides down-slope as a function of the fraction of the rest of the glacier being anchored to the substrate. Figure 6 displays the case when only the top part of the ice-block is anchored. The rest is free to slide down-slope without friction. It is evident from this figure that the anchored part is not enough to keep the ice-block in place. It breaks near the substrate and the entire block slides downslope. If the time is re-scaled as explained above, the velocity reaches roughly 5 m day −1 , which is comparable to observed rates (about 10-100 m d −1 , Cuffey and Paterson, 2010) for surging glaciers.
In the opposite limit, when there are additional frictional anchoring points on the lower part of the slope, the ice-block cannot move, but remains stuck. In this case only a smaller part of the block, near the lower edge, calves, fragments and flows a limited distance down the slope. As this layer of highly fragmented ice flows, it gets thinner and the force driving it down-hill decreases and it slows down. This is displayed in Fig. 7. www.the-cryosphere.net/7/1591/2013/ The Cryosphere, 7, 1591-1602, 2013 In order to quantify the above, Fig. 8 shows kinetic energy (E kin ) as a function of time for varying amounts of frictional contact on the lower slope as described in Figs. 6 and 7. The various lines in Fig. 8 correspond to two different phases, a surging phase, when the entire ice-block slides down the slope, and a quiescent phase, when only part of the front of the ice-block fractures and flows down-hill. For some of the simulated cases there appears to be a phase-transition during the simulation run time. In these cases the kinetic energy is initially in the quiescent phase and at some point in time the kinetic energy suddenly increases and rapidly approaches that of the surging phase. Sometimes surging does not appear for the entire block, and the kinetic energy only increases part-way towards the surging phase before stabilising.
Finally, to highlight the difference between the single calving event represented by Fig. 3, and the surging glaciers in Fig. 6, the fragment size distribution was calculated for the surging glacier simulation (Fig. 4). In this case the fragment size distribution was equivalent to that usually found for a grinding process (Åström, 2006).  5. Snapshots of a deforming ice block simulated with Elmer/Ice and our particle model. In the Elmer/Ice simulation (red markers) the snapshots are from time steps 0, 3, 5 and 7 yr and the particle model snapshots (red area) are from corresponding time steps. The size of the block is 30 m × 30 m and the time span of the Elmer/Ice simulation is roughly 10 8 s but only 10 3 s in the particle simulation.

Discussion
The new model we introduced in this paper is certainly not feasible to incorporate into ice sheet models given the extensive computing power required. It may however be used to investigate details of calving processes and relationships such as dependence of crevassing rate and fragment size on the water depth at the calving front, or the presence and influence of water in crevassses on fracture processes. The model also has considerable potential to test and improve parameterizations of fracturing and calving used in continuum models. The resolution of many models simply does not include small ice-cliff failure, and, since calving and fracturing are essentially discontinuous processes, introducing them into continuum models is problematic. Cuffey and Paterson (2010) summarised the situation as: most models either let ice shelves advance to the edges of the model grid, or assume that ice shelves terminate at a prescribed water depth (400 m typically). For marine-terminating glaciers that are not fully floating, most models either assume that calving rate increases with water depth, or constrain the ice front thickness instead of the calving rate. However, recent progress has been considerable in the field of parameterizing crevassing by weakening the ice rheology in a damage model.
Our discrete particle formulation may be seen as a complementary method to statistical continuum damage approaches that have been applied to ice shelves (e.g. Borstad et al., 2012) or to mountain glacier calving (e.g. Jouvet et al., 2011). This can be illustrated by, for example, Levermann et al. (2012) who formulated the vertically averaged ice strain rate tensor, which can be determined from the spatial derivative of the remotely sensed velocity pattern. His model can The Cryosphere, 7, 1591-1602, 2013 www.the-cryosphere.net/7/1591/2013/ then be tuned to observations of specific ice shelves, and no other observations are needed for the model to "carry itself forward" into the future. That is, it can "predict" calving without any other observation inputs (the necessary inputs would all come from the dynamic ice sheet model).
Over broad areas of an ice shelf, the viscosity is reduced by crevasses (e.g. along the flow units coming from different tributary ice streams and glaciers), or the ice may be strengthened by the presence of sub ice shelf freeze-on of ocean water or weakened by melting processes. The crevasses can be readily seen in imagery, and these images can be used to tune models for ice viscosity and fracture initiation stress in specific ice shelves or tributary ice streams to give similar patterns of both crevassing and velocities as observations (Albrecht and Levermann, 2012a, b). Above we demonstrated the importance of basal friction to the behaviour of the particle model. Fast outlet ice streams and surging glaciers are governed by the physics of basal sliding. In temperate glaciers (i.e. glaciers with temperatures at the pressure melting point) sliding behaviour is often tightly interlinked with basal hydrology. On the continental ice sheets, the fast flowing ice streams, and outlet glaciers owe their speed to basal sliding in addition to internal ice deformation. Schoof (2009) showed that a variety of friction laws converged on the Coulomb friction law in appropriate parametric limits which can usefully describe a plastic till rheology. The motion of ice streams appears to depend critically on the distribution and nature of regions of high drag ("sticky spots", Alley, 1993). It is not known what controls the present configuration of these features, though presumably they are related to the bed roughness and geometry either directly as a bedrock bump, or by routing water supply and till properties. Inverse methods can be used to determine the spatial variability of basal friction (Raymond and Gudmundsson, 2009;Morlighem et al., 2010;Petra et al., 2012;Arthern and Gudmundsson, 2010;Jay-Allemand et al., 2011;Schäfer et al., 2012), in an analogous way to the damage weakening of shear margins and crevasse damage on ice shelves. The cases when the ice-block is the entire time in either the surging or the quiescent phase the energy is indicated by discrete markers. For the blocks for which a phase transition appear, the kinetic energy is represented by continuous lines.
The discrete particle model we use here clearly suffers from lack of three dimensional geometry; hence it is presently limited to the testing or verification of the parameterization used in continuum models. It is plausible to incorporate basal friction laws that could mimic more accurately a plastic basal till with sticky spots. Our model prediction for the particle size distributions can be readily tested in observational data on fragmentation and calving.
Appendix A NUMFRAC -a particle based code for fracture mechanics in nonlinear materials A1 Elastic model NUMFRAC is a particle based code which has been designed to simulate fracture in nonlinear, e.g. visco-elastic or plastic, materials. The code has been used for modelling brittle fracture (Åström, 2006), cytoskeleton dynamics (Åström et al., 2010), and mechanics of fiber networks (Åström et al., 2012). A material is modelled as a set of particles that are connected by interaction potentials. The code allows for arbitrary arranged "packings" of particles of different size, kind and shape. For example, a rather loose random packing of spheres of different size and stiffness can be used as a model of an isotropic porous material with structure fluctuations. In contrast, a close packed hexagonal arrangement of identical spheres models a dense uniform anisotropic material. Geometry, interaction, and all other relevant parameters can be set separately for all particles and particle-particle interactions, Fig. 9. Model test result for Glen's flow law. Five constant shear stresses were applied to a test block and strain rate was measured and viscosity extracted. The melting-refreezing probabilities are set so that the pre-factor, A, in Glen's flow law is of the order 10 −10 , rather than 10 −24 , which would be the realistic range for ice: the simulated material in its shear dependent behaviour resembles Glen's flow law but has an overall viscosity that is about 14 orders of magnitude lower than that of ice. This corresponds to a viscosity roughly at the same order of magnitude as that of bitumen or asphalt. The simulation data is represented by red markers and fitting the function, x B /A gave A ≈ 1.5 × 10 −10 , B ≈ −2.0 as it should.
making the code an excellent tool for investigating mechanical behaviour of strongly disordered materials.
To initialise a simulation, the particles are arranged and packed to form the desired material. Thereafter the connections between particles are determined, e.g. by elastic beams, which then determine the potentials.
The equations of motion may vary from case to case. A simple example is given by: where M is a mass-matrix containing the masses and moments of inertia of the particles, r i ,ṙ i ,r i are the position, velocity and acceleration vectors of particle i, including rotations. r ij are the corresponding position vectors for all particles j that are connected to particle i. C is a damping matrix containing damping coefficient. K is the stiffness matrix and F i is the sum of other forces acting on particle i. These forces may include gravitation, buoyancy, atmospheric and hydrodynamic/static forces, etc. An example of an interaction potential between two particles is an Euler-Bernoulli (E-B) beam. The elastic energy of an E-B beam can be written as (1/2)kx 2 , where x is the displacement vector containing translational and rotational displacements of two connected mass points. If the components of the displacement vector, x 1 and x 4 are the displacements of the two connected mass points along the axis of the beam The Cryosphere, 7, 1591-1602, 2013 www.the-cryosphere.net/7/1591/2013/ that connects them, x 2 and x 5 the displacements in the perpendicular direction and x 3 and x 6 the rotations of the mass points, the stiffness matrix k for small deformations of a linear elastic E-B beam in two dimensions is: where I is the moment of inertia of the beam cross-section, L the length of the beam, E Young's modulus and A the crosssection area of the beam. In 3-D the matrix must be extended to a size of 12×12 entries. Euler-Bernoulli beams overestimate the bending/shear stiffness for short and fat beams. That is, beams with an aspect ratio smaller than or roughly equal to 20. For smaller aspect ratios, Timoshenko beams may be used.
Another efficient and easy potential which we have also used for the ice application below, is to define a tension stiffness which is a harmonic potential with respect to the distance between two particles, and another harmonic potential with respect to node rotations away from the axis that connects the mid-points of the two particles. These two stiffness moduli can be directly linked to the macroscopic Young's modulus and Poisson ratio of the material to be modelled. This potential is a type of shear-beam. Calving behaviour seems to be rather insensitive to which beam model is used. No significant differences in results could be detected. Timoshenko beams are, from a theoretical point of view, the most exact, but demand slightly more computations than Euler-Bernoulli beams or shear beams. For the results shown here we used the latter ones.
Once the interactions for the particle-particle connections are determined and the appropriate stiffness matrices are assembled, the global stiffness matrix, K, for the entire material body to be simulated is assembled by expanding, with zeros, the matrices for the individual particle-particle interaction matrices, and adding them together. Since the orientation is different for each interaction pair and the orientations change over time, the last term in Eq. (A1) on the left side, is typically nonlinear and can be obtained, for the present example, via: where dr ij is the time dependent incremental displacement vector for connected particles i and j , T is a time dependent rotation matrix for converting a global coordinate system to the orientation of the beam axis connecting i and j . In order to implement Eq. (A1) on a computer, the differential equation must be rewritten as a difference equation via: where r is now the global position vector containing all particle translations and rotations. t is time and t is the time incremental, i.e. time-step. This discrete form is easily implemented on a computer and given time and space-dependent boundary conditions, and r(t = 0), the time development r(t) can unambiguously be calculated.
For most fragmentation simulations there is further a need to determine a fracture criterion. If "softening" potentials are used, for example, the Lennard-Jones potential, there is no such need since the attractive force between mass points vanishes continuously. For harmonic potentials the fracture criterion must be defined explicitly.
Choosing a proper fracture criterion is far from trivial. A general elliptical criterion (Zhang and Eckert, 2005) that includes the "classical" fracture criterion of Tresca, von Mises, Mohr-Coulomb and the maximum normal, i.e. hydrostatic pressure, stress criterion is rather useful. This criterion states that a material under tension fails at locations where in which σ I is the normal stress, or pressure determined by the first invariant of the stress tensor and σ I I is shear stress determined by the second invariant. σ 0 and τ 0 are material dependent constants. Because both shear and normal stresses are easily defined for beams it is straightforward to implement this fracture criterion in the beam model. A fracture limit can either be defined for every beam separately or as a limit on the average stress over several beams. For a single beam, σ I can simply be set to the stress along the axis of the beam and shear σ I I the off-axis stress. For several beams (e.g. for all beams connected to a single mass point) the average stress tensor can be divided in a trace-less, isotropic and a diagonal part which then defines local shear and normal stress, respectively.

A2 Plastic and visco-elastic models
In order for the particle model to be able to model not only fracture of elastic material, but also plastic and visco-elastic materials more elaborate interaction potentials, as compared to purely elastic ones, must be used. The microscopic mechanism behind plasticity and visco-elasticity is that local interactions cannot just be broken but also reform in configurations different from the original one. This is the general principle behind irreversible material deformations like viscous www.the-cryosphere.net/7/1591/2013/ The Cryosphere, 7, 1591-1602, 2013 910 kg m −3 1000 kg m −3 0.1-1.0 GPa 10 −4 -10 −3 10 3 kg s −1 (10 −5 -10 −3 ) s flow and plastic deformation. In the particle model, plasticity can be introduced via an ordinary yield/fracture stress criterion combined with an additional criterion that allow new contacts between particles to be formed. These new contacts may be formed between any two particles that come close enough to each other as the material deforms. Close enough means here that the particles touch, or almost touch, each other. With this mechanism, the material stress will not only be dependent on the load as in an elastic material before any fracture has taken place, but also on the load history. A visco-elastic material differs from a plastic material in that the material stress depends not only on load and load history but also on time. This means that for example a constant, time independent, external displacement, stress will slowly relax to zero in a visco-elastic material.
It is rather straightforward to incorporate this behaviour in particle models. To demonstrate the procedure, Glen's flow law can be obtained as a result of the adopted shear-beams rheology via, so called, "melting-freezing probabilities". If these probabilities are set right, they can produce a stressdependent viscous flow obeying Glen's law for flow rate, D = A(T )σ n−1 e t D , where A(T ) is a temperature dependent Arrhenius factor, σ e is the second invariant of the deviatoric stress tensor, D is the strain-rate tensor, t D is the deviatoric stress tensor, and n ≈ 3. To derive Glen's law for our model in an approximate scalar form, we assume that melting events are random and uncorrelated which means that the probability for an event to appear with a time interval t obeys an exponential probability function P = 1−e −λ t . This function gives the probability of melting a beam during a timestep t, where λ is the rate of melting (melting events per second). Each time a beam is broken all elastic stress and strain, = σ/E, will be relaxed. Here , σ and E are strain, stress and Young's modulus of the beam. As long as the model material is under compression and the particles are close-packed, the model material will be practically incompressible. This implies that the deformation of the simulated material will be by shear flow only and of the form D ≈ t D E λ. By choosing λ = 2A U a 2 E 2 ≈ Aσ 2 E we obtain Glen's law D = A(T )σ 2 t D . Here U is the elastic energy of the beam and a is the diameter of the discrete blocks (and length of a beam). The elastic energy of a beam U is roughly equal to the elastic energy density σ 2 2E , multiplied by the effective volume of a beam a 2 . Correspondingly, dynamic viscosity is µ ≈ E λ ≈ 1 Aσ 2 . Figure 9 shows viscosity as a function of shear stress in a simulation series where a varying shear stress was applied to a test square.

A3 Application to glacier mechanics
The application to real world ice problems is still under development. Some general issues that are important input to the model are still poorly understood for glacier ice. The most important is probably the correlation and distribution of structural and mechanical disorder. Disorder, in this context, means the impurities, weaknesses, cracks that are abundant in natural ice. Such impurities have an affect on local stiffness moduli and fracture thresholds on the "particle scale", which is presently cubic-metre.
Since some of the model key parameters are only known with accuracies of an order of magnitude, there is no point to model other parameters to higher accuracy. Therefore we simply use, e.g. 910 3 kg m −3 for ice density, a 10 % density difference between ice and water for buoyancy. The Young's modulus of ice is a particularly difficult parameter. For pure crystalline ice the value of Young's modulus E is roughly 10 10 Pa (Schulson, 1999), while the effective stiffness S of polycrystalline glacier ice commonly drops to about 10 9 Pa (Vaughan, 1995). It is quite obvious that when the "quality" of ice weakens, stiffness, in particular when measured as tension, will decrease drastically. Eventually, as ice begins to melt or is highly damaged, it practically vanishes. The same kind of uncertainty is relevant for fracture strain of ice. For the present simulations we have used values [10 8 − 10 9 ] Pa for Young's modulus and [10 −4 − 10 −3 ] fracture strain. The damping coefficient for colliding ice-blocks is another parameter that is not known to us. It is, however, reasonable to expect that colliding square-metre ice blocks would not bounce back very efficiently, and therefore we use values similar to, or slightly less than that for critical damping of the harmonic potential. The water drag coefficient for a squaremetre object is roughly 10 3 kg s −1 . Typical simulation parameters are displayed in Table A1.
From a computational point of view the most problematic aspect of the simulations are the short time-spans that can be simulated. In practice, the only reasonable way to slightly reduce this problem is to use a viscosity that is lower than that of ice, but which still resembles Glen's flow law. The relaxation time near ice cliffs does not have a simple definition but is likely in the range of weeks or months. Since the calving time scale is of the order of a few seconds we can drastically reduce the viscous relaxation time without getting much interference between the two, i.e. the Deborah number defined as viscous relaxation time over typical calving time is of the order of 10 6 (Reiner, 1964). In other words, reducing the viscous timescale by 10 5 from days to seconds still leaves the viscous and the fracture timescales well separated.

A4 Computational implementation
The code has been constructed by the authors TR,JA,TT and does not use any commercial sub-routines or libraries. The code is written in Fortran and C++ and utilises MPI for parallel computing. The scalability to large glaciers, i.e. many particles, is close to perfect. The most severe limiting factor is the time-step, which is of the order 10 −4 − 10 −3 s. There is no useful way to formulate the time discretization in parallel and, at present, the code may, in the best case, simulate about one hour of glacier dynamics within 24 h of computing. The only reasonable way to speed up the computations is to calibrate the particle model such that the macroscopic viscosity of the model is considerably lower than that of ice as explained above.
Since computation of time cannot be made parallel, it is important, for computational efficiency to maximise the time-step without violating the stability of the computation. Since there are no explicit temperature fluctuations in the simulations and no external heat sources it is easy to check for stability by recording the total energy and check for energy conservation during simulations.
The particle size sets a very intuitive resolution limit for the simulation. At present the particles are of the order of one cubic metre, but the model is inherently scale invariant. As long as the equations of motion, Eq. (A1) remain unchanged, all parameters can be re-scaled without changing the simulation results. This means, the same simulation can be viewed as, e.g. a cubic-metre or cubic-millimetre simulation if the length unit in r, M, C, K and F are all changed accordingly. For example, if the length unit is decreased, M and F , which is dominated by gravity, are reduced in proportion to the volume of the particles, while K is reduced only in proportion to the surface area of the particles. This means that small glaciers are much more stable, under gravity, than large ones.