Combining damage and fracture mechanics to model calving

Introduction Conclusions References


Introduction
The discharge of ice from the Greenland and Antarctic ice sheets increased sharply in recent decades (Shepherd et al., 2012), due to either intensified submarine melting, or to an increase in the rate of iceberg calving.Recent observations show that despite some regional differences, the ice loss is almost equally distributed between these two sink terms (Rignot et al., 2010;Depoorter et al., 2013).In 2013, ice loss in Antarctica due to iceberg calving has been estimated at 1321 ± 144 gigatonnes per year (Depoorter et al., 2013) and in Greenland between 2000 and 2005, at 357 gigatonnes per year (Rignot and Kanagaratnam, 2006).These figures could be higher, as destabilization of the front can exert strong positive feedback on glacier dynamics.Indeed, the abrupt collapse of the front can affect the equilibrium of the whole glacier, leading to both upstream thinning and acceleration due to loss of buttressing, which in turn further increases ice discharge (Gagliardini et al., 2010).The collapse of Larsen B ice shelf in 2002 (Scambos et al., 2004) or the disintegration of the floating tongue of Jakobshavn Isbrae, on the west coast of the Greenland ice sheet in the same year (Joughin et al., 2008a) are two examples on how perturbations near the ice front can affect the upstream and even the grounded part of adjacent glaciers.To accurately predict future changes in the ice sheet, improved understanding and representation of the processes that occur at the front are necessary, especially those concerning iceberg calving.
Most studies dealing with iceberg calving follow the approach proposed by Benn et al. (2007a, b).These authors introduced a criterion suggesting that calving occurs when a crevasse penetrates the glacier below the water level.The computation of the crevasse penetration depth is based on the work of Nye (1957), and depends on the equilibrium between longitudinal stretching (opening term) and cryostatic pressure (closing term).This so-called "crevasse depth" criterion has been applied to individual marine-terminated glaciers in Greenland and Antarctica, and enabled the successful reproduction of variations in the front (Nick et al., 2010;Otero et al., 2010;Nick et al., 2013;Cook et al., 2014).Although the model of Nick et al. (2010) accounts for basal crevasse propagation, which improves the ability of a model to reproduce observed behaviour, it is based on an instantaneous stress balance combined with an empirical criterion for calving.Consequently, it cannot account for the physical processes related to fracture propagation, such as the stress concentration at the tip of crevasses, and does not take into account the role of the stress history on the accumulation of Published by Copernicus Publications on behalf of the European Geosciences Union.
J. Krug et al.: Combining damage and fracture mechanics to model calving fractures, which may influence both the time of occurrence and the amplitude of calving.
Another approach to modelling calving uses discreteelement models (Bassis and Jacobs, 2013;Åström et al., 2013).This type of model has produced interesting results in terms of calving processes and iceberg size distribution.However, coupling such models with finite-difference or finite-element glacier or ice-sheet models is not straightforward, and their high computation cost limits their usefulness.
In the last few years, some authors have used continuum damage mechanics to represent both the development from micro-defects in the ice to the development of macro-scale crevasses, and their effects on the viscous behaviour of the ice while keeping a continuum approach (Duddu and Waisman, 2012;Borstad et al., 2012;Albrecht and Levermann, 2014).Initially developed for metal deformation (Kachanov, 1958), damage mechanics has recently been applied to ice dynamics to study the appearance of a single crevasse (Pralong et al., 2003;Pralong and Funk, 2005;Duddu and Waisman, 2013) or to average crevasse fields (Borstad et al., 2012).Linear elastic fracture mechanics (van der Veen, 1998a, b) has been used to describe the rapid propagation of surface and bottom crevasses through the ice.This approach has rarely been used in ice-sheet modelling.The reason is that a realistic representation of crevasses requires a high mesh refinement that is usually difficult to achieve when modelling large ice masses.
Here we consider an approach that combines damage mechanics and fracture mechanics.The proposed physically based calving model can cover both the accumulation of damage as the ice is transported through the glacier, and the critical fracture propagation that characterizes calving events.The slow development of damage represents evolution of viscous ice at long timescales, while the use of fracture mechanics makes it possible to consider calving events that occur at shorter timescales, during which the ice can be considered as a purely elastic medium.The description of the physics used is presented in Sect.2, and covers the beginning of damage and its development, fracture propagation and its arrest criterion.In Sect.3, the results of sensitivity tests on Helheim Glacier are discussed and suggestions are made for further research.

Ice flow and rheology
We consider an incompressible, isothermal and gravitydriven ice flow in which the ice exhibits a non-linear viscosity.The ice flow is ruled by Stokes equations (i.e.Navier-Stokes equations without an inertial term), for the momentum and the mass balance: where σ represents the Cauchy stress tensor, g the gravity force vector, ρ i the density of ice and u the velocity vector.The Cauchy stress tensor can be expressed as a function of the deviatoric stress tensor S and the isotropic pressure p, with σ = S − p I and p = −tr(σ )/3.Ice rheology is represented by a non-linear Norton-Hoff type flow law called Glen's flow law: This equation links the deviatoric stress S to the strain rate ε.
The effective viscosity η is written as where I 2 ε2 represents the square of the second invariant of the strain rate tensor, A is the fluidity parameter and E is an enhancement factor, usually varying between 0.58 and 5.6 in ice-flow models (Ma et al., 2010), depending on the fabric and on the stress state.

Boundary conditions
The upper surface is defined as a stress-free surface.In the coordinate system (x, y, z), it obeys the following equation: where z s refers to the elevation of the upper surface, and (u s , v s , w s ) are the surface velocities.The surface mass balance a s is prescribed as a vertical component only.As we disregard any effect of atmospheric pressure, normal and tangential stresses at the surface are zero: Subscripts n and t i respectively refer to normal (pointing outward) and tangential directions.Similar to the upper free surface, the bottom surface evolution is described by where (u b , v b , w b ) are the basal velocities, and a b represents the vertical component of the basal mass balance (melting or accretion).At the bed, the glacier can be either grounded or floating.The grounded part of the glacier undergoes shearing stress, which is represented by a non-linear Weertman-type friction law: where C and m = 1/3 are respectively the friction coefficient and exponent.u b is the norm of the sliding velocity with n b the normal outward-pointing unit vector to the bedrock.Where the glacier is floating, the free surface is forced by an external sea pressure normal to the surface: where ρ w is the water density, l w the sea level, and z b is the elevation of the bottom surface.The position between the grounded and floating part of the basal boundary, i.e. the grounding line, is part of the solution and is computed solving a contact problem following Durand et al. (2009) and Favier et al. (2012).The inverse method described in Jay-Allemand et al. ( 2011) is used to infer the basal friction coefficient C by reducing the mismatch between observed and modelled surface velocities.
The ice front is defined as a third free-surface, which can also undergo melting.In analogy to the other free surfaces this gives where (u f , v f , w f ) are the frontal velocities, and a f is the frontal mass balance.The ice front is exposed to the water pressure below sea level and is stress-free above sea level.These Neumann conditions read The physical and numerical parameters used in this paper are listed in Table A1.Some boundary conditions are specific to the 2-D flow line application applied in Sect.3, and are described in detail in Sect.3.2.

Continuum damage mechanics model
Continuum damage mechanics (CDM) was introduced by Kachanov (1958) to quantify the degradation of mechanical properties resulting from the nucleation of internal defects such as micro-cracks or voids.In this case, the internal defects must be small compared to the representative volume element over which damage is considered.In the case of ductile failure, when the propagation of a macro-crack is associated with the nucleation of defects (voids) ahead of the crack tip in a so-called fracture process zone (FPZ), damage mechanics has been used to describe the propagation of the macro-crack itself.So far, there is no evidence for such ductile failure in ice.Consequently, in the present work, we used damage mechanics to deal with the effect of a field of crevasses on ice flow and not to describe the propagation of an individual crevasse.Here, as stated by Lemaitre et al. (1988), we consider that CDM describes the evolution of phenomena in the medium from a virgin state to the triggering of macroscopic fractures.In this approach, the material is always considered as a continuous material, even when the level of damage is very high.Slow deformation is typically encountered in glaciology, when ice flows slowly under its own weight following the surface slope.CDM has been successfully used in ice-flow models to deal with some glaciological problems such as the flow acceleration of large damaged areas or the opening of crevasses in hanging glaciers (Xiao and Jordaan, 1996;Pralong et al., 2003;Pralong and Funk, 2005;Jouvet et al., 2011;Borstad et al., 2012).
The principle of CDM models is based on the use of a damage variable, usually denoted D, which represents the degradation of mechanical properties (stiffness, viscosity, etc) resulting from a population of defects whose effect is averaged at mesoscale (i.e. the grid size, a few metres).When considering an anisotropic approach, damage must be represented as a second-order tensor (Murakami and Ohno, 1981;Pralong and Funk, 2005).However, following Pralong and Funk (2005), here we consider isotropic damage as a first approximation.In this case, the state variable D is a scalar quantity that varies between 0 and 1.While the ice is considered undamaged for D = 0, full damage is attainable when D approaches 1.
To describe the effect of damage on the ice flow, an effective deviatoric part of the Cauchy stress tensor is introduced: This effective stress can be understood as the original force applied to an effective undamaged area only.Using the equivalence principle of Lemaitre et al. (1988), strain is affected by the damage only through an effective stress entering the constitutive equation (see Eq. 3) in place of S, and thus changes the expression of the Cauchy stress tensor σ in the Stokes equation (see Eq. 1).Thereby, there is no need to define an effective strain rate.
Our approach models the non-linear viscous flow of ice.The viscous behaviour of the ice is described using Glen's flow law, which links the deviatoric part of the Cauchy stress tensor to the strain rate tensor.Consequently, when accounting only for viscous deformation, damage to the ice only affects the deviatoric part of the Cauchy stress tensor, and not the cryostatic pressure.

Damage evolution
Damage is a property of the material at the mesoscale.It is therefore advected by the ice flow, and evolves over time depending on the stress field.To take this evolution into account, we recommend an advection equation: Compared with Pralong and Funk (2005)'s mathematical formalism, this equation does not involve the spin tensor, as the damage defined here is isotropic and results in a scalar variable.The right-hand side represents a damage source term f (χ ) that can be written as a function of a damage criterion χ and a numerical parameter B, henceforth referred to as damage enhancement factor: In Sect.3, we present some experiments to test the sensitivity of the CDM model to the damage enhancement factor.The choice of the damage criterion is crucial for the representation of damage increase, and its physical expression is a critical step in the formulation of a CDM model.Commonly used criteria are the Coulomb criterion (Vaughan, 1993), the von Mises criterion (Albrecht and Levermann, 2012), and the Hayhurst criterion (Pralong and Funk, 2005;Duddu andWaisman, 2013, 2012).However, these criteria are not necessarily applicable to damage to ice: the Coulomb criterion is used for a representation of frictional process under compressive loading (e.g.Weiss and Schulson, 2009).The von Misès criterion, expressed in terms of deviatoric stresses only, is independent of cryostatic "pressure" (either positive or negative) and symmetric in tension and compression.As the ice does not behave in the same way under tensile or compressive stresses (Schulson and Duval, 2009), we believe this criterion is not suited to describe ice fracturing processes.The Hayhurst criterion (Hayhurst, 1972;Gagliardini et al., 2013a), which involves the maximum principal stress, cryostatic pressure and the von Misès stress invariant, was designed to describe creep damage in ductile materials and allows for damage under uniaxial compression.However, as we want to describe crevasse opening under tension only, we believe that this criterion is not suitable for this work.
Instead, we use a pure-tensile criterion, described as a function of the maximum principal Cauchy stress σ I .This choice is consistent with the fact that we want to describe crevasse opening under pure tension, and is supported by Rist et al. (1999), who argued that crevasses tend to open normal to the direction of maximum tensile stress.This criterion is also able to represent the broad range of crevasses observed on glaciers, including splaying crevasses.The implementation of another criterion in the model would be straightforward, and would be an interesting parameter to investigate in future work.The criterion reads as follows: Here σ th represents a stress threshold for damage initiation.
The corresponding envelope of the damage criterion is represented in the space of Mohr circle in Fig. 1.The stress threshold for damage initiation σ th corresponds to the overload that must be applied in order to reach the ice strength and initiate degradation.To account for subgrid-scale heterogeneity, we introduce some noise on σ th : σ th = σ th ± δ σ th , where δ σ th σ th follows a standard normal distribution with a standard deviation of 0.05.σ th is the mean stress threshold, and usually reaches several tens of kilo-Pascals (Pralong and Funk, 2005).To avoid negative values for σ th arising from the sub-grid-scale heterogeneity, the lower bound σ th > 0 is applied.The sensitivity of the model to this parameter is discussed in Sect.3.This formulation of the damage criterion implies some limitations to the calving model.In particular, it does not account for shear and compressive failure mechanisms.However, it remains consistent with the approach of Benn et al. (2007b), who stated that it is the longitudinal stretching associated with longitudinal velocity gradients that exert firstorder control on the development of crevasses in glaciers.Moreover, it is consistent with the fracture mechanics approach explained in Sect.2.3, which considers crevasses opening in pure tension only.

Viscosity modification
As pointed out by Pralong et al. (2003) and Pralong and Funk (2005) in the case of Alpine hanging glaciers, ice flow is altered by the accumulation of micro-defects in the ice: damage softens the ice and accelerates the creep.This softening is taken into account by introducing effective stress in Glen's law.
Introducing the effective deviatoric stress tensor S and taking into account the equivalence principle, as described in Sect.2.2, Eq. ( 3) gives When the expression of the effective deviatoric stress tensor given in Eq. ( 12) is replaced by the deviatoric stress given in Eq. ( 8) it becomes By identification with Eq. ( 3), the enhancement factor can be linked E with the damage D, such that For undamaged ice (meaning D = 0), E = 1, the flow regime is unchanged.When the damage increases (D > 0), E > 1, the viscosity of the ice is reduced, and the velocity of the flow consequently increases.This formulation of the enhancement factor is consistent with the expected behaviour, and has already been used in previous studies, such as Borstad et al. (2012).The damage then evolves under the effect of the stress field, where the ice undergoes a critical tensile stress σ th in the direction of the principal stress, and exerts positive feedback on the velocity field.

Fracture mechanics
Continuum damage mechanics can be used to deal with the degradation of ice viscosity with increasing damage at long timescales.It can be understood as a way to simulate subcritical crevasse nucleation and propagation (Weiss, 2004) at a mesoscale, and their role in creep enhancement.However, calving events are triggered by rapid propagation of pre-existing fractures at very short timescales, and speeds can reach a significant fraction of the speed of sound.Thus, this process cannot be represented by viscous rheology (Weiss, 2004).Instead, at such short timescales, the medium should be considered as elastic.In these conditions, linear elastic fracture mechanics (LEFM) is a useful tool to account for these features and matches observed crevasse depths (Mottram and Benn, 2009).The application of LEFM to the penetration of surface crevasses, originally introduced by Smith (1976), has since been used by several authors (Rist et al., 1999;van der Veen, 1998a, b;Nath and Vaughan, 2003).
Here, a LEFM model is combined with the previously described damage model to achieve the formulation of a calving law taking into account fast and slow processes controlling glacier fracturing.In LEFM, three modes of crack propagation can be considered: Mode I opening, Mode II sliding and Mode III tearing.In the following, only the opening mode (Mode I) is considered.

LEFM theory
The key physical parameter of LEFM is the stress intensity factor K. van der Veen (1998a) proposed an expression for K I in an idealized case where the opening stress is constant in the vertical direction.For the opening mode I considered here, in the coordinate system (x, y, z), where σ xx is the horizontal component of the Cauchy stress tensor, d is the crevasse depth and β is a parameter that depends on the geometry of the problem.The crack is considered to propagate vertically.In the ideal case introduced by van der Veen (1998a), fracture propagation is a function of the difference between the opening deviatoric stress S xx , resulting from horizontal velocity gradients, and the cryostatic pressure (creep closure) σ p = ρ i g(H − z), corresponding to the weight of the ice, hence However, when considering real cases, the opening term S xx (and so σ xx ) is not constant over depth z or the lateral coordinate y.Consequently, the appropriate formula to calculate the stress intensity factor for an arbitrary stress profile σ xx (y, z) applied to the crack is given by the weight functions method (Labbens et al., 1974): This formula relies on the use of the superposition principle: in the case of linear elasticity, the value of the stress intensity factor at the tip of the crack can be seen as the sum of contributions of all individual point loads along the crack length.In our case, instead of considering the value of the along-flow component of the deviatoric stress tensor at the tip of the crack, we multiplied it by the weight function β(z, d, H ) at each vertical coordinate and integrated it over the initial crevasse depth (Labbens et al., 1974).In this way, the effect of an arbitrary stress profile on the stress intensity factor can be taken into account.
In LEFM theory, a fracture is able to propagate downward in the ice if the stress intensity factor is greater than the fracture toughness K Ic .Toughness is a property of the material and depends to a great extent on ice porosity.Several experiments have been carried out to link the value of K Ic to porosity (Fischer et al., 1995;Rist et al., 1996;Schulson and Duval, 2009).Inferred values ranged from 0.1 to 0.4 MPa m 1/2 .In our calving model, we use a constant value of 0.2 MPa m 1/2 .The sensitivity of the model to this value is discussed in Sect.3.3.4.
The weight function β(z, d, H ) depends on the geometry of the crevasse, and consequently on the problem considered.Among the weight functions for various crack and notched geometries that have been proposed, we chose the one corresponding to an edge crack in an infinitely wide plane plate   A1 for a list of the parameters.(Glinka, 1996).A complete description of the weight function and an illustration of the geometry is given in Fig. 2 and Appendix A.

Critical damage contour and fracture initiation
From Eq. ( 16), it is easy to understand that the triggering of crevasse propagation requires a combination of both sufficient tensile stress and sufficiently deep initial crevasse to exceed fracture toughness.In our model, the size of preexisting flaws is dictated by a contour of critical damage on the near surface of the glacier, where damage reaches a critical value D c .For application to the LEFM theory, we consider that the depth of this damaged layer corresponds to the initial crevasse depth d (see Fig. 3).One should bear in mind that the value of D c is also a threshold that needs to be set.The sensitivity of the model to this parameter is tested in Sect.3.
Compared to the work of van der Veen (1998a, b), we do not consider the presence of water-filled crevasses for the initiation of crack propagation, nor the formation of basal crevasses.These questions are briefly addressed in Sect.3.4.Without these features, however, our model is sufficient to provide a lower bound for crevasse propagation.

Fracture arrest
Once the conditions for fracture initiation are fulfilled, we consider that the crevasse propagates vertically.In van der Veen (1998b), crevasses propagate downward as long as the equation K I ≥ K Ic is satisfied, thus assuming that K I = K Ic represents both a crack propagation and a crack arrest criterion.Such an arrest criterion is probably misleading, as the stress intensity factor at arrest, though non-zero, is always lower than the stress intensity factor at propagation (Ravi-Chandar and Knauss, 1984), mostly because dynamic effects have to be taken into account for the arrest condition.Therefore, following Ravi-Chandar and Knauss (1984), we use a crevasse arrest criterion K I < K Ia , with K Ia = α K Ic and 0 < α < 1.The value of α for ice is unknown.In the following, we arbitrarily set α to 0.5.The sensitivity to this value is discussed in Sect.3.3.4.In this simplified LEFM framework, calving would theoretically occur only if K I remains larger than K Ia down to the bottom of the glacier.However, as a result of the cryostatic pressure and hydrostatic pressure, K I becomes negative before reaching d = H . Several authors have proposed alternative criteria to overcome this inconsistency.Benn et al. (2007a) proposed a first-order approach considering that calving of the aerial part of the glacier occurs when a surface crevasse reaches the sea level.This criterion is supported by two observations.First, Motyka (1997) showed that calving of the aerial part occurs when the crevasse reaches sea level, usually followed by the calving of the subaqueous part.Secondly, a surface crevasse reaching sea level may be filled with water, if a connection exists to the open sea (Benn et al., 2007b).In this case, the crevasse will propagate further downward.Indeed, the water adds a supplementary pressure ρ w g d w , where d w represents the height of the water in the crevasse, which is equal to the distance between sea level and the crack tip.Added to the tensile opening stress, this supplementary hydrostatic pressure counterbalances the cryostatic pressure and/or the ocean water back pressure.Consequently, the opening full stress S xx dominates the force balance and one would expect the crevasse to propagate downward.We calculated K I in the case where a crevasse fills with water on reaching sea level.The resulting stress intensity factor becomes positive (up to one kilometre upstream from the front) over the entire glacier depth, thus supporting Benn's criterion.This parameterization was successfully applied by Nick et al. (2010) on an idealized geometry, and we chose to use the same criterion.Thereby, the stress intensity factor is computed at a depth d f equal to sea level.
This framework has two consequences.First, the stress profile σ xx used to calculate K I | d f for the arrest criterion is estimated before the propagation of the crevasse.This propagation modifies σ xx , but this effect is not considered here.Second, if the condition nothing happens in the model whereas one would expect some brittle crevasse propagation down to d f .In other words, our model considers LEFM to describe calving but not to simulate crevasse propagation upstream of the calving front.
The calving model described in the previous sections is summarized in Fig. 4.
The CDM and LEFM models were implemented in the finite element ice sheet/ice flow model Elmer/Ice.Further information regarding Elmer/Ice can be found in Gagliardini et al. (2013b).One specific improvement of the model was incorporating an automatic re-meshing procedure when a calving event occurs.The mesh is rescaled, and the variables are interpolated onto the new geometry.The corresponding method is described in supplementary material in Todd and Christoffersen (2014).In addition, we implemented a horizontal interpolation scheme for the damage and the stress field, allowing calving to occur between nodes, thus reducing the horizontal mesh dependency.

A case study
The model was applied to Helheim Glacier, a fast-flowing well-monitored glacier in Southeast Greenland.The abundance of observations there allowed us to compare and constrain our model parameters against past glacier evolution.For our model development, we focus on a two-dimensional flowline problem.

Data sources
As stated by Andresen et al. (2011), Helheim Glacier's front position remained within a range of 8 km over the last 80 years.During this period, the glacier underwent several advances and retreats.In particular, Helheim underwent a major retreat between 2001 and 2005, before creating a floating tongue which advanced again between 2005 and 2006 (Howat et al., 2007).In the last decade, it has been intensively surveyed and studied (Luckman et al., 2006;Joughin et al., 2008b;Nick et al., 2009;Bevan et al., 2012;Cook, 2012;Bassis and Jacobs, 2013), and is therefore an interesting case study to calibrate a calving model.
As we focused on changes in the front and on representing calving, we only needed a bedrock topography covering the last 10 km up to the front to capture any basal features that could influence the behaviour of the front.For this reason, we chose to follow the work by Nick et al. (2009) and to use their data set in which the last 15 km of the glacier are well represented.In this data set, the initial front position corresponds to the May 2001 pre-collapse geometry.In addition, we chose to consider the glacier as isothermal near the ice front, where we set a constant temperature of −4.6 • C (following Nick et al., 2009).It should be noted here that a constant negative temperature may be inconsistent with the  calving criterion, as it will freeze any water in crevasses, and thus prevent the glacier from establishing a hydraulic connection with the proglacial-water body.However, this constant value was chosen in order to produce a velocity field consistent with observations, and, for sake of simplicity, we did not take the vertical variations in temperature into account.A constant surface mass balance a s was taken from Cook (2012) who fitted direct observations from stakes placed on the glacier between 2007 and 2008, which are assumed to represent the annual surface mass balance.

Flow line and numerical specifications
Following our notation system, the ice flows along the horizontal "x direction" and the "z direction" is the vertical axis.
The geometry covers the last 30 km of the glacier, with an initial thickness varying between 900 m at the inlet boundary, and 700 m at the front.Using the metric from Nick et al. (2009), the beginning of the mesh was located at kilometre 319, and the front at kilometre 347 (see Fig. 5).This geometry was discretized through a structured mesh of 10 500 quadrilateral elements, refined on the upper surface and at the front.The element size decreased from 150 to 33 m approaching the front in the horizontal direction, and from 3.3 to 68.0 m from top to bottom.Sensitivity tests were carried out to optimize the mesh size in the vertical direction, and glacier behaviour appeared to converge with increasing refinement, since only negligible differences appeared when refining more than 3. allowed for accurate fracture initiation and damage advection, as well as relatively rapid serial computation.We used an Arbitrary Lagrangian Eulerian (ALE) method to account for ice advection and mesh deformation, with a time step of 0.125 day.Sensitivity to the length of the time step was also undertaken.Below the chosen value, no more significant difference appeared in the behaviour of the glacier, or in the size and frequency of the calving event.
The specific boundary conditions used for the 2-D application are described below, using the boundary conditions presented in Sect.2.1.2.
The basal friction C was inferred from the 2001 surface velocity data presented in Howat et al. (2007).
For melting at the calving front, we chose an ablation function that increased linearly with depth, from zero at sea level to 1 m day −1 at depth.This constant melting parameterization was inspired by the work of Rignot et al. (2010) on four West Greenland glaciers.
The inflow boundary condition (x = 319 km) did not coincide with the ice divide.As we only considered the last kilometres before the front, we assumed that ice velocity is constant at the upstream boundary, and we imposed a vertically constant velocity profile of u x = 4000 m a −1 , in agreement with the observed surface velocity in Howat et al. (2007).However, this was not the case for the inlet flux, which can vary depending on changes in inlet boundary thickness.
When dealing with a 2-D flowline representation of the ice flow, we had to take some three-dimensional effects into account.In particular, lateral friction along the rocky-margins of the glacier can play a significant role by adding resistive stress to the flow.Here, we modified the gravitational force using a lateral friction coefficient k, as proposed by Gagliardini et al. (2010).
This coefficient depends on Glen's flow law parameters A and n, as well as on the channel width W , which was taken from Nick et al. (2009).The lateral friction coefficient was applied over the whole length of the glacier, and covered the entire lateral ice surface.It did not account for the tributary glacier, which merges with the principal stream near the centre of the flow segment.
Even if the velocity and the surface topography were known and correspond to the state of the glacier observed in May 2001, some relaxation was necessary to obtain a stable steady state (Gillet-Chaulet et al., 2012).Thus, for 8 years, we let the geometry adjust to the prescribed boundary conditions and to the previously inverted basal friction coefficient.

Model calibration: sampling strategy
The model was calibrated by varying the three parameters σ th , B, and D c in a given range.In order to distinguish the effect of the different damage parameters, we used a Latin hypercube sampling (LHS) method with a combination of 16 distinct parameters for (σ th , B), that we reproduced for three values of D c .
The stress threshold σ th was estimated to be in the range of [0.01, 0.20] MPa.The lower bound was near that given in Pralong and Funk (2005).If σ th > 0.20 MPa, the modelled stresses are not high enough to reach the damage threshold, and damage, when it occurs, is too weak.
The damage enhancement factor B was chosen within [0.5, 2.0] MPa −1 a −1 .The explanation for this wide range is the following: when the damage criterion χ is positive, the corresponding increase in damage releases the stress level in the ice.In reality, this process happens continuously, such that the stress level cannot exceed the edge of the envelope defined by the green line in Fig. 1.However, as our model deals with a finite time step, the stress level may jump in the "damaging" area (area shaded grey in Fig. 1).The role of "damaging" is then to "push" the stress back to the edge of the envelope: the rate of this stress displacement is set by the value of B. Ultimately, this parameter should be set such that the stress is displaced right to the edge of the damaging envelope.However, this value is difficult to evaluate, especially because stress relaxation occurs through non-linear viscous flow.
The critical damage value D c has already been inferred in several studies (Pralong and Funk, 2005;Borstad et al., 2012;Duddu and Waisman, 2013).Following these studies, we used three values of D C (0.4; 0.5; 0.6).

Model calibration: spin-up
Damage can be produced anywhere in the glacier.As we needed to obtain a steady state for the damage field, we had to let the damage created upstream be advected to the front.The model was therefore spun-up for 8 years: the front was kept fixed at its initial position, without submarine frontal melting, or calving.After 8 years, the front was released, and frontal melting and calving were activated.

Model response
Over the last century, Helheim Glacier has probably undergone several advance and retreat cycles (Andresen et al., 2011).Observations of sand deposits showed that variations in the front did not exceed 10 km.Knowledge of the potential trigger mechanisms for these cycles is still poor: according to Joughin et al. (2008b) and Andresen et al. (2011), the recent retreat may have been caused by higher summer air and ocean temperatures.Yet, the role of the calving activity under these changes remains unclear.For these reasons, we did not try to reproduce the exact chronology of the most recent retreat of Helheim Glacier.Instead, we studied the dynamic behaviour of the model with respect to the different sets of parameters, and tried to distinguish between unrealistic and realistic behaviours.To distinguish between different model behaviours, we ran the simulations for 10 years using the parameter sets.The model response was then split into the three classes illustrated in Fig. 6.The first class of model behaviour corresponds to when calving does not often happen.The glacier advances too far and builds up a floating tongue several kilometres in length (blue line, Fig. 6).The second class illustrates prolific calving, leading the front to retreat far upstream (yellow line).The last class of behaviour is consistent with observations, where the front is punctuated by regular or irregular calving events, forcing the glacier to remain within an acceptable range of values (red line).
This threefold classification of glacier behaviour was generalized to the 48 simulations.To eliminate aberrant behaviour, we ran a sanity check by considering plausible sets of parameters as those that lead to a simulated front position within the range [340, 350 km].The results are presented in Fig. 7, in the parameter space (B, σ th , D c ).The three classes can again be distinguished -blue triangles, yellow circles and red diamonds, representing respectively the case where the front exceeds 350 km, the case where the front retreats less than 340 km, and the case where the front remains between 340 and 350 km.
The steady advance of the front with no, or with only a few calving events (blue triangles) can be explained first by considering the parameters (B, σ th ).These underlying simulations are characterized by a low value of B and/or a high value of σ th .This means that either damage incrementation was too low, or the stress threshold was too high to trigger damage.In these cases, damage production may be not sufficient to reach the calving criterion D = D c .In addition, when σ th is too high, damage may only increase in the area where the traction is very high, i.e. on top of bumps, in the immediate vicinity of the surface.As a consequence, the damage never reaches a sufficient depth to trigger calving.
In contrast, the too rapid retreat of the front (yellow circles) can be explained as follows: when B is high and/or σ th is low, triggering damage is easy, and increments quickly, leading to a major damage at the glacier surface.As a consequence, the criterion D > D c can be more easily reached, leading to a too rapid sequence of calving events.
The chosen range for σ th produces acceptable results.When related to the tensile strength of snow and ice, it agrees with constraints from previous studies, especially those described in Borstad (2011), who inferred the tensile strength of snow and firn at the surface to be in the range of 10-50 kPa.Vaughan (1993) found values in the range of 100-400 kPa in the vicinity of crevasses, which is slightly higher than our range.However, deeper knowledge of the firn density at Helheim surface would probably help further constrain this parameter.
Finally, parameter D c is used to control whether the fracture propagates.If it is low, the conditions for crevasse propagation are easily reached (as soon as the criterion on LEFM is fulfilled).If D c is too high, damage may never reach sufficient depth to initiate fracture propagation (see Fig. 7).As a consequence, this value controls the location of the calving front.However, it is tightly linked with (B, σ th ) through the production of damage upstream and must be chosen in the same range as the level of damage at the front.
However, discretization into three classes of parameters should not mask the continuous behaviour of the glacier, depending on the value of (σ th , B, D c ), and the boundaries between classes are not abrupt.Thus, simulations may present a front position out of the range [340, 350 km] (and are considered to be unrealistic), but still have a general dynamics that is not far from realistic behaviours.

Realistic behaviour
The acceptable parameter combinations are located on the 3-D diagonal represented by the red diamonds in Fig. 7.During these simulations, the front remained between the two limits and behaved consistently with observed values.The simulation corresponding to the parameter set σ th = 0.11 MPa, B = 1.30MPa −1 , and D c = 0.50 (red star in Fig. 7), can be seen as an example of consistent behaviour.As represented in Fig. 8a and b, damage increases in the area where the traction is high enough to exceed σ th , mostly over bumps.This is consistent with observations of glaciers flowing over slope ruptures (Pralong and Funk, 2005) at a rate high enough to reach the critical damage value D c at the front and trigger calving (see Fig. 8c and d).Mottram and Benn (2009) investigated crevasse depth in the vicinity of the front of a freshwater calving glacier in Iceland, and showed that most crevasses were no deeper than 10 m.In our model experiment, damage reached depths of 5 to 15 m before calving occurred.As described previously, the value of d must be large enough to account for critical crevasse propagation, and be consistent with observations.Focusing on this parameter set, we highlighted the following main features of the glacier front dynamics.First, at short timescales, the calving activity can be described through a small-size/high-frequency event distribution.This activity is characterized by quasi-periodic retreats of the front of limited extent (between 50 and 150 m), associated with a period varying from about 3 to 15 days (and only 1 day for very rapid front retreat, see Sect.3.3.5).These features can be seen in the inset in Fig. 6).Additionally, some outliers (calving events larger than 500 m) interrupted the dynamics of the front from time to time.
The glacier front dynamics is also characterized by a larger-amplitude/lower-frequency oscillation (see Fig. 6), with a period of ∼ 500/600 days and a 500 m amplitude.This asymmetric cycle exhibits a slow advance but a faster retreat, mostly due to a rapid succession of calving events (but not related to the outliers).Despite the complex interactions and feedback between damage parameters and ice flow, we have a possible explanation for these features.Damage reduces ice viscosity, accelerates ice flow, and adjusts glacier geometry.This softening also reduces the stress level as well as the damage rate.In addition, the damage is advected until it reaches conditions in which crevasse propagation can trigger calving events (probably because of bending stresses at the front, or processes referred to as second-order processes in Benn et al., 2007b).This calving event results in an immediate increase in the tensile stress in the vicinity of the front (due to the unbalanced forces at the new ice cliff).The subsequent increase in damage where the geometry readjusts may trigger a "cascade" of calving events leading to a major (cumulative) retreat of the glacier front over a limited timescale (a few days/weeks), until the front reaches a position where the damage is too low to trigger calving.Repeating this process several times leads to the calving cycles highlighted in Fig. 6.To investigate the strength of this feature, we ran this simulation (red curve in Fig. 6) eight times, using only different local disorders.The Fourier spectrum of the position of the front (fast Fourier transform visible in Fig. 9) exhibits a lowfrequency peak for a period (500-600 days) approximately equal to the time that ice needs to be advected between the two main bumps visible in Fig. 8, which is likely related to the above-mentioned "cascade" mechanism.This feature was apparent with different parameters sets (B, σ th , D c -not shown).In Fig. 9, the smooth wide peak visible at around 3-15 days characterizes high-frequency calving events.This observation, combined with the mechanism proposed above, could explain the cycles in the position of the front, and also suggests that the surface geometry (driven by the bedrock topography) controls the calving dynamics by modifying the rate at which damage primarily develops in the ice.Thus, in the experiments presented in this paper, the slow cycles in Fig. 6 are not related to variability in external forcing (which remained constant), but are the result of internal glacier dynamics.This leads us to conclude that a variation of several kilometres in the position of the front is not necessarily the consequence of an external forcing, although it may well be triggered by it.

Sensitivity analysis
The analysis described in the previous section for one parameter set was extended to the 12 simulations that successfully passed the sanity check, and showed that the model's response remained qualitatively unchanged even when the damage parameters changed.Altogether, we www.the-cryosphere.net/8/2101/2014/The Cryosphere, 8, 2101-2117, 2014 identified 6534 distinct calving events.Characteristic calving event sizes emerged from the distribution, similar to those visible in the inset in Fig. 6, which are illustrated in Fig. 10.When Cook et al. (2014)'s model was applied on Helheim Glacier, the result was a mean calving event size of 450 m, with a frequency of around 14 days.On Store Glacier, with the same zero-forcing conditions, Todd and Christoffersen (2014) obtained a mean size of 80 m, and a frequency of around 8 days.Our 50-150 m range and our 3-15 days period is thus within the range of these two modelling studies.
However, it should be recalled that this plot cannot be interpreted as an iceberg's size distribution.Indeed, one has to distinguish between the front retreat and the size of the resulting iceberg(s) which may differ considerably, as the calved portion of ice may fragment into many icebergs and/or capsize.Although the distribution of the distance the front retreats could be an interesting parameter to calibrate the model, it would require continuous tracking of the front position of the actual glacier, as discrete determination of the position (Joughin et al., 2008b) may bias the estimation of the retreat of single events, particularly for the small population sizes.We are not aware of the existence of such a data set for Helheim Glacier.
In addition, we estimated the influence of the other parameters discussed above.The results showed that the model is only slightly sensitive to the K Ic parameter.As ice toughness increased, it became more and more difficult to trigger fracture propagation.However, changing the value of ice toughness did not alter the general behaviour of the system.Varying α did not have a significant impact.Indeed, the arrest criterion was computed at a depth equal to sea level.At this depth, the stress intensity factor is greater than the ice toughness (not shown), and thus, higher than K Ia , whatever the value of α.
Finally, model sensitivity to the initial heterogeneous distribution of micro-defects introduced in Sect.2.2.1 was tested.The standard deviation of the distribution of stress threshold δσ th σ th varied within the range [0.005; 0.2].The results were robust under these changes, although the spatial and temporal variability of the front position was modified (not shown).The higher the standard deviation, the greater the variation in the position of the front.

Outlook and further improvements
Several improvements could be made to the degree of physics included in the model.Their implementation should be straightforward, as soon as the knowledge of their underlying properties is sufficient to better constrain the resulting physical parameters.

Water-filled crevasses
Among the improvements that could be made to the model, the presence of surface melt water in crevasses is the most obvious.This has a major effect on their propagation, as it represents a supplementary hydrostatic force in the stress balance tending to favour crevasse opening.van der Veen (1998b) showed that a water-filled crevasse is able to propagate through the whole thickness of the glacier and trigger calving as soon as the lever of water reaches around 10-20 m below the surface.
The implementation of this feature in our modelling framework is straightforward.To do so, we just have to add a water pressure term in the expression of the Cauchy stress tensor, depending on the water level in the crevasse d w : Then, the stress intensity factor is computed following Eq.( 16), using σ xx (y, z) instead of σ xx (y, z).However, adding this feature requires knowledge of the water depth in the crevasse, or the rate of water input, which are both currently poorly constrained and difficult to measure.

Basal crevasses
At the current stage of development, the model does not include the representation of basal crevasses.These are cited as a possible explanation for the production of large tabular icebergs, as they require a full-thickness fracture.This basal propagation is only possible if the glacier is near or at flotation, such that the tensile stress is high enough to trigger a fracture (van der Veen, 1998a).Nick et al. (2010) included basal crevasses in a physically based calving law, by adding the pressure of water filling the crevasses in the estimation of the penetration length.Once the penetration of top and basal crevasses covers the full thickness of the glacier, calving occurs.In our experiments, we did not observe a sufficient tensile stress to generatedamage at the base of the glacier, and work is currently underway to study this process.Technically, capturing the fracture propagation at the base of the glacier relies on integrating the effect from of the water pressure at glacier base: where H p is the piezometric head, which represents "the height above the base to which water in a borehole to the bed will rise" (van der Veen, 1998a).But this implementation requires knowledge of the value of H p , which is difficult to estimate if the glacier is grounded.Considering H p = s l could be a reliable upper bound in the vicinity of the front.
In its present state, however, our framework still proposes a lower bound for calving event size and front retreat.

Crevasse shielding
Considering the case of highly crevassed glacier surfaces and closely spaced crevasses, our fracture propagation framework could not rely on a parameterization that only considers lone crevasse propagation.In this case, problems arise from the fact that the stress concentrations at crevasse tips are reduced by the presence of neighbouring crevasses (shielding effect).The consequence is that it would become harder for a crevasse to propagate downward under the same stress conditions.van der Veen (1998b) proposed a parameterization that could improve our model.In the case of a constant tensile stress profile S xx , the stress intensity factor is modified as a function of the distance between neighbouring crevasses l and the crevasse depth d: where S xx is the tensile stress, and D(L) is a weight function that depends on L. However, this kind of parameterization would require estimation of l, which cannot be obtained from our modelling framework at present.Once an estimation of l is obtained in another way, the introduction of such shielding effect in Eq. ( 16) will be straightforward.
The major effect of this development would be a supplementary delay in the time and the position of calving events, probably resulting in calving events occurring nearer to the calving front, where the tensile stress is higher, and as such resulting in smaller size distribution of icebergs.

Conclusions
In this work, we combined continuum damage mechanics and linear elastic fracture mechanics to propose a processbased calving model.This model is able to reproduce the slow development of small fractures leading to the appearance of macroscopic crevasse fields over long timescales, while considering ice as a viscous material.It also allows for the elastic behaviour of breaking ice, consistent with the critical crevasse propagation that triggers calving events, characterized by very short timescales.The model was applied to Helheim Glacier, which allowed us to constrain the parameter space for the most important free parameters.Within these constraints, the simulated ice front led to cyclic calving events, and the position of the glacier front remained in the range observed over the last century.
The three decisive parameters are the damage initiation threshold σ th , the effect of damage on viscosity, quantified by the enhancement factor B, and the critical damage variable D c .The first two parameters must be in a range that allows sufficient damage to occur upstream before reaching the front by advection.Then, the maximum value of damage should be close to D c up to a certain depth in order to trigger calving.
It should be borne in mind that this sensitivity test is based on the response of one specific glacier to a poorly known external forcing and with limited observations.Under these conditions, we show that some sets of parameters definitely produce reasonable behaviour, but these parameters should not be transferred to another configuration without being sure that the response of the model is realistic.Despite this limitation, this calving model based on realistic physical approaches gives reliable results and could be easily implemented in classical ice-flow models.
The calving process described in this paper is immediately driven by the variation in longitudinal stretching associated with horizontal velocity gradients, producing a first-order control of the calving rate, as stated by Benn et al. (2007b).Local aspects, involving undercutting or force imbalance at ice cliffs are described as second-order calving processes.Using this model, further work could be undertaken to increase our knowledge of these second-order phenomena.As stated by Glinka (1996), the weight function for the computation of the stress intensity factor depends on the specific geometry of the initial crack.For an edge crack in a finite width plate, the weight function is given by .
The weight function depends on three numerical parameters, polynomial functions of the ratio d : H , namely

Figure 1 .
Figure1.Damage envelope in the space of principal stresses.σ I and σ II represent the first and the second principal stress, respectively, while σ th is the stress threshold.The shaded area corresponds to stress conditions in the calving model that will trigger damage.

Figure 2 .
Figure 2. Crevasse shape.See TableA1for a list of the parameters.
Figure 3. (a) Shape of a grounded glacier and (b) close-up of one crevasse.The red curve illustrates the contour of critical damage D = D c , for which we computed the along-flow component of the Cauchy stress tensor σ xx multiplied by the weight function, and integrated over the crevasse depth d.

Figure 4 .
Figure 4. Algorithm of the calving model where t refers to the time step.The blue shape indicates the area of CDM application, where the ice is viscous, and the orange shape corresponds to the LEFM domain of application, where the ice is elastic, representing fracture propagation and calving events.

Figure 5 .
Figure 5. Glacier location and geometry.(a) Location on the Greenland ice sheet (red dot).(b) Close-up of the Helheim Glacier front and the flow line (red curve) studied here.(c) Mesh constructed for this segment of the flow line.The starting position corresponds to the front position at 347 km (May 2001).The blue line indicates sea level.(d) Close-up of the upper surface of the glacier at the calving front.

Figure 6 .Figure 7 .
Figure 6.Positions of the calving front over the 10-year period.Each colour corresponds to a set of parameters σ th , B, and D c .Blue represents an advance with almost no calving (σ th = 0.15 MPa, B = 0.84 MPa −1 , and D c = 0 .6); yellow represents a major retreat (σ th = 0.10 MPa, B = 1.75 MPa −1 , and D c = 0.4) ; red represents a behaviour that is consistent with observations (σ th = 0.11 MPa, B = 1.30MPa −1 , and D c = 0.50).The inset shows a close-up of the consistent case.

Figure 8 .
Figure 8. State of Helheim Glacier after 365 days of simulation for the set of parameter (σ th = 0.11 MPa, B = 1.30MPa −1 , and D c = 0.50) corresponding to the red star in Fig. 7. χ is the damage criterion, which is positive where damage will occur.(a) Damaging areas of Helheim Glacier and (b) close-up of the red rectangle; (c) damage field and close-up of the red rectangle (d).

Figure 9 .
Figure 9. Frequential response of variation in the calving front in the eight simulations using the set of parameters (σ th = 0.11 MPa, B = 1.30MPa −1 , D c = 0.50), computed over 10 years.Gray lines represent the frequency spectrum for each of the eight simulations, and the solid black line shows the corresponding mean.

Figure 10 .
Figure10.Distribution of the size of the calving events corresponding to the 12 realistic simulations.Different modes can be related to the different combinations of damage parameters.The axis representing the calving event size is truncated at 200 m.However, a few events, up to 2000 m were recorded.

Table A1 .
Physical and numerical parameters.Tunable parameters are in italic.