A particle based simulation model for glacier dynamics

This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. A particle based simulation model for glacier dynamics Åström, Jan; Riikilä, Timo; Tallinen, Tuomas; Zwinger, Thomas; Benn, D.; Moore, J.C.; Timonen, Jussi Åström, J., Riikilä, T., Tallinen, T., Zwinger, T., Benn, D., Moore, J.C., & Timonen, J. (2013). A particle based simulation model for glacier dynamics. Cryosphere, 7(5), 1591-1602. https://doi.org/10.5194/tc-7-1591-2013 2013 The Cryosphere, 7, 1591–1602, 2013 www.the-cryosphere.net/7/1591/2013/ doi:10.5194/tc-7-1591-2013 © Author(s) 2013. CC Attribution 3.0 License. The Cryosphere Open Access A particle based simulation model for glacier dynamics J. A. Åström1, T. I. Riikilä2, T. Tallinen2, T. Zwinger1, D. Benn3,4, J. C. Moore5,6,7, and J. Timonen2 1CSC – IT Center for Science, P.O. Box 405, 02101, Esbo, Finland 2Department of Physics, University of Jyväskylä, P.O. Box 35 (YFL), 40014, Jyväskylä, Finland 3Department of Geology, University Centre in Svalbard, 9171 Longyearbyen, Norway 4School of Geography and Geosciences, University of St Andrews, Fife, KY16 8ST, UK 5State Key Laboratory of Earth Surface Processes and Resource Ecology, College of Global Change and Earth System Science, Beijing Normal University, Beijing, China 6Arctic Centre, University of Lapland, PL122, 96100 Rovaniemi, Finland 7Department of Earth Sciences, Uppsala University, Villavägen 16, Uppsala, 75236, Sweden Correspondence to: J. A. Åström ([email protected]) Received: 29 January 2013 – Published in The Cryosphere Discuss.: 6 March 2013 Revised: 15 August 2013 – Accepted: 19 August 2013 – Published: 8 October 2013 Abstract. A particle-based computer simulation model was developed for investigating the dynamics of glaciers. In the model, large ice bodies are made of discrete elastic particles which are bound together by massless elastic beams. These beams can break, which induces brittle behaviour. At loads below fracture, beams may also break and reform with small probabilities to incorporate slowly deforming viscous behaviour in the model. This model has the advantage that it can simulate important physical processes such as ice calving and fracturing in a more realistic way than traditional continuum models. For benchmarking purposes the deformation of an ice block on a slip-free surface was compared to that of a similar block simulated with a Finite Element full- Stokes continuum model. Two simulations were performed: (1) calving of an ice block partially supported in water, similar to a grounded marine glacier terminus, and (2) fracturing of an ice block on an inclined plane of varying basal friction, which could represent transition to fast flow or surging. Despite several approximations, including restriction to twodimensions and simplified water-ice interaction, the model was able to reproduce the size distributions of the debris observed in calving, which may be approximated by universal scaling laws. On a moderate slope, a large ice block was stable and quiescent as long as there was enough of friction against the substrate. For a critical length of frictional contact, global sliding began, and the model block disintegrated in a manner suggestive of a surging glacier. In this case the fragment size distribution produced was typical of a grinding process. 1 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” (Jouvet et al., 2011; 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, 1980; Schulson, 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 Published by Copernicus Publications on behalf of the European Geosciences Union. 1592 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 Fis applied and (b) bend when a torque Tis 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. 2 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 1m. 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 σ2 I σ2 0 +σ2 II τ2 0 ≥1,(1) in which σIis the normal stress, or pressure determined by the first invariant of the stress tensor and σII is shear stress determined by the second invariant. σ0and τ0are 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, Kis 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 The Cryosphere, 7, 1591–1602, 2013 www.the-cryosphere.net/7/1591/2013/ J. A. Åström et al.: A particle based simulation model for glacier dynamics 1593 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: M¨ ri+C˙ ri+X j γij C0˙ rij +X j γ0 ij Krij =Fi,(2) where Mis the diagonal mass-matrix containing the masses (i.e. volume times density Vρ) and moments of inertia (ρRr2dV) of the particles. Vectors riand ˙ ridenote the position and velocity of particle iand rij and ˙ rij are the corresponding relative vectors for particles iand j. The diagonal damping matrix Ccontains damping coefficients for drag, i.e. drag force = (1/2)ρv2cDS, where vis velocity, S is cross-sectional area of the object to which the drag is applied, and cDis the Reynolds number dependent drag coefficient. The other damping matrix, C0, contains coefficients for the inelastic collisions. The parameter γij is zero for particles not in contact and unity for particles in contact, and γ0 ij is unity for connected particles and zero otherwise. The stiffness matrix is denoted by Kand Fiis 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=109Nm−2, density ρ=910kgm−3and fracture strain c=(1–5)×10−4(Schulson, 1999). The damping coefficient for collisions is C0=105Nsm−1.Crepresents the drag force on ice falling into water in a calving event. A typical value for this parameter is 103kgs−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, 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 etD, where A(T ) is a temperature dependent Arrhenius factor, σeis the second invariant of the deviatoric stress tensor, Dis the strain-rate tensor, tDis 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−4s, 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 1594 J. A. Åström et al.: A particle based simulation model for glacier dynamics Fig. 2. A flow chart showing the algorithm in a schematic fashion. Details may vary. 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. 3 Results Figure 3 shows snapshots of a calving 30m×30m 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, 108Pa 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 20m 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−4s. 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–20s, 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)dsis the number of fragments found in the size interval [s−ds/2,s +ds/2], and sis the number of particles in a fragment. This is, obviously, proportional to the volume of the fragment, i.e. roughly the same amount in m3. 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): n(s) ∝s−(2D−1)/D exp(−s/s0), (3) where sis again the fragment size or volume, Dis the dimension (D=2 for the simulations here and D=3 for real The Cryosphere, 7, 1591–1602, 2013 www.the-cryosphere.net/7/1591/2013/ J. A. Åström et al.: A particle based simulation model for glacier dynamics 1595 Fig. 3. Snapshots of a calving ice-block. The size of the block is 30m×30m. The block rests on a soft substrate that efficiently hinders the block from sliding, thus resembling, e.g. a muddy sea floor. The block is immersed in 20m of water. The simulation times of the snapshots are indicated. glaciers), s0is 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 105 times lower than in the Elmer model, leading to 105times 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 200m×50m on a 18◦slope. Again, the viscosity was roughly 105lower than realistic values for ice, i.e. an Arrhenius factor ≈5×10−19 s−1Pa−3. We would thus expect that the strain rates will be roughly 105higher 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 (24h=0.864×105s). 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 5mday−1, which is comparable to observed rates (about 10–100md−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 1596 J. A. Åström et al.: A particle based simulation model for glacier dynamics Fig. 4. (A) The fragment size distribution for simulations like the one in Fig. 3. The figure displays the results for three different simulations and both the distribution early during the calving event and later when the fragments have come to rest floating on the water. The line is the distribution predicted from Eq. (3). (B) The fragment size distribution n(s) for a surging glacier. In this case n(s) ∝s−2.3. In order to quantify the above, Fig. 8 shows kinetic energy (Ekin) 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). Fig. 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 7yr and the particle model snapshots (red area) are from corresponding time steps. The size of the block is 30m×30m and the time span of the Elmer/Ice simulation is roughly 108s but only 103s in the particle simulation. 4 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 (400m typically). For marine-terminating glaciers that are notfullyfloating, most modelseitherassume thatcalvingrate 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/ J. A. Åström et al.: A particle based simulation model for glacier dynamics 1597 Fig. 6. Snapshots of a 200m×50m ice-block on a 18◦slope. The red line marks a high friction contact. The colours of the ice-block represent elastic areas where the elastic beams between particles have not been broken (green) and areas where beams are no longer intact (gray). The pressure of an upstream ice slab on the slope above the ice-block is not shown in the figure. 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 Fig. 7. Similar as in Fig. 6, but with several “frictional anchors”, indicated by red markers, also on the lower part of the slope. 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. www.the-cryosphere.net/7/1591/2013/ The Cryosphere, 7, 1591–1602, 2013 1598 J. A. Åström et al.: A particle based simulation model for glacier dynamics Fig. 8. The kinetic energy (Ekin) as function of time for “surging” simulations exemplified by Figs. 6 and 7 for 7 different simulations. 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, xB/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: M¨ ri+C˙ ri+X j Krij =Fi,(A1) where Mis a mass-matrix containing the masses and moments of inertia of the particles, ri,˙ ri,¨ riare the position, velocity and acceleration vectors of particle i, including rotations. rij are the corresponding position vectors for all particles jthat are connected to particle i.Cis a damping matrix containing damping coefficient. Kis the stiffness matrix and Fiis 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)kx2, where xis the displacement vector containing translational and rotational displacements of two connected mass points. If the components of the displacement vector, x1and x4are 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/