Calving glaciers and ice shelves

ABSTRACT Calving, or the release of icebergs from glaciers and floating ice shelves, is an important process transferring mass into the world’s oceans. Calving glaciers and ice sheets make a large contribution to sea-level rise, but large uncertainty remains about future ice sheet response to alternative carbon scenarios. In this review, we summarize recent progress in understanding calving processes and representing them in the models needed to predict future ice sheet evolution and sea-level rise. We focus on two main types of calving models: (1) discrete element models that represent ice as assemblages of particles linked by breakable bonds, which can explicitly simulate fracture and calving processes; and (2) continuum models, in which calving processes are parameterized using simple calving laws. With a series of examples using both synthetic and real-world ice geometries, we show how explicit models are yielding a detailed, process-based understanding of system physics that can be translated into predictive capability via improved calving laws. Graphical Abstract


Introduction
Understanding the processes responsible for the release of icebergs (calving) from glaciers and floating ice shelves has been a long-standing problem in glaciology. Despite early recognition of the potential instability of marineterminating glaciers [1,2] and some clear-sighted modelling studies [3,4], progress on the problem was sporadic until around 20 years ago. To a large degree, this reflects the complexity of calving processes. Calving can occur in lakes, fjords or other marine waters, as well as from terrestrial ice cliffs in polar and high-mountain regions. Calving glaciers may be fast or slow flowing, have grounded or floating termini, and occur in a wide range of temperate to polar environments. Additional obstacles to progress on the calving problem included the difficulty of making comprehensive observations in dangerous and remote locations, and the hitherto limited computer hardware and software available for detailed analysis. Recently, calving research has been spurred by increasing awareness of the sensitivity of marine-terminating glaciers to climate change, most dramatically illustrated by the disintegration of the northern Larsen ice shelves in 1995 and 2002, and flow acceleration and retreat of many calving glaciers in Greenland in the early twenty-first century [5,6]. The last decade has seen considerable progress in both observational and modelling capability, and a comprehensive understanding of calving has begun to emerge. Reliable predictive capability, if not quite within reach, is at least within sight.
In this paper, we review current understanding of ice fracture and calving processes from lake-and marine-terminating (or tidewater) glaciers, and ice shelves. We focus on recent advances in modelling calving, including explicit models of fracture processes and the simple calving laws required by large-scale ice sheet models. We conclude with a discussion of the future prospects for calving research.

Calving processes 2.1. Fracture
All calving is a consequence of the formation and growth of fractures. Mechanical deformation of ice occurs in contrasting regimes depending on the applied stress and the time scale of deformation. Over long time scales, polycrystalline glacier ice will flow as a non-linear fluid, through the aggregate effect of strains within and between grains [7,8]. Typical observed strain rates for viscous flow lie in the range 10 −11 to 10 −8 s −1 , but may be two orders of magnitude faster during transient events such as surges [9]. At higher strain rates (faster than 10 −3 s −1 ), viscosity can be largely ignored and ice behaves as an elastic-brittle material [10]. Between the viscous and elastic-brittle strain-rate regimes, ice deformation can occur via a complex mixture of creep and ductility combined with fracturing at many different length scales.
Brittle fracture can occur in three modes [11,12]: Mode 1: crack opening (tensile mode); Mode II: crack sliding (crack propagation parallel to the applied shear stress); and Mode III: crack tearing (crack propagation at right angles to the applied shear stress). In ice, tensile failure typically produces singular, sharp-sided cracks ( Figure 1). Once failure has occurred the ice is unable to support stress. In contrast, shear failure occurs where at least one of the principal stresses is compressive, and typically produces networks of microfractures that grow gradually in density and begin to interact. Finally, the microcracks coalesce and the intact ice between them collapses to form a compressive shear band in which further fragmentation occurs by grinding [7,13]. The shear band has finite strength, and can continue to support some level of applied stress. The stress or strain-rate thresholds for fracture initiation in ice are poorly defined, largely due to the highly variable properties of glacier ice, including temperature, presence of debris and other impurities, crystal anisotropy, pre-existing damage, loading history and other factors. Observations of the formation of tensile crevasses on glaciers indicate that the threshold for Mode I failure lies in the range 90-400 kPa [14]. For shear failure the threshold is somewhat higher. The Young's modulus of ice is roughly 10 GPa, and the failure strain in the brittle regime for laboratory samples is 10 −3 to 10 −4 [15], which gives a tensile fracture stress threshold of 1-10 MPa. The large discrepancy in failure stress for crevasse formation and laboratory-scale ice is partly a reflection of the size-and structure-dependent fracture strength, which is not uncommon for brittle and heterogeneous materials [16].
Fractures (crevasses) are very widespread on glaciers [11], but calving requires crack propagation sufficient to isolate ice blocks from the main body of the glacier. Crack propagation can be limited by a number of factors. For example, surface crevasses rarely penetrate deeper than a few tens of meters even under high-tensile stresses, because tensile deviatoric stress is increasingly offset by ice overburden pressure at depth. Pressure exerted by water within fractures provides an additional opening force, allowing surface crevasses to penetrate to greater depths [17]. Additionally, where ice is close to flotation, water-filled basal crevasses can penetrate almost through the full thickness of a glacier, potentially leading to calving [18,19]. High-tensile stresses at the ice surface are an insufficient condition for calving, and stress conditions throughout the ice mass must be conducive to continued crack growth. The process of crack growth is of course influenced by the presence of the crack itself; loss of strength in the cracked region means that any body forces acting on the glacier must be supported by the remaining intact ice, thus modifying the initial stress distribution. If the elastic strain energy released by fracturing exceeds the energy required to break more bonds at the crack tip, the fracture will propagate. If not, the crack will stabilize. Thus, the size, location and frequency of calving events reflect the sum of stresses arising from the interaction between the glacier and its environment, including such factors as buoyancy, undercutting by subaqueous melting, and boundary friction. The influence of these and other factors on calving styles is discussed in Section 2.3.

Damage
The development of micro-scale cracks and macro-scale crevasses alters the bulk properties of glacier ice and is collectively termed damage. The accumulation of damage is of great importance for many glaciological processes. For example, intense fracturing at the margins of ice streams and ice shelves reduces their ability to support stress, diminishing their ability to buttress inland ice and increasing rates of ice flow [20]. Damage (including the slow growth of sub-critical cracks) also increases the likelihood of bulk failure and calving [21][22][23]. An extreme end-member of damaged glacier ice is mélange, consisting of floating masses of calved icebergs that may be unbonded or cemented together by sea ice. Mélange can play an important role in suppressing calving by resisting the rotation and overturn of icebergs [24][25][26].
Effective treatment of damage requires solutions to two problems: (1) formulating an effective rheology for fractured ice and (2) determining functions to describe damage evolution. The rheology of ice damaged by fracture can be approximated as a viscoelastic fluid with a lower viscosity and Young's modulus than intact ice. This is typically done by invoking the so-called equivalence principle, that postulates that only the undamaged part of a damaged material can carry stress. Within such a meanfield approximation both viscosity η and Young's modulus E decrease linearly from their undamaged values η 0 and E 0 with increasing damage [22]. For completeness, an expression for the Poisson ratio ν d ð Þwould be required. As a first approximation ν ¼ ν 0 is independent of d. Bulk constitutive relations for damaged ice have been developed by Dansereau et al. [27], Borstad et al. [28] and Sun et al. [20], and a damage-mechanics-based calving model was introduced by Krug et al. [29]. Catastrophic failure (e.g. a calving event) occurs at some threshold value of damage, whereupon ice can no longer support the body forces acting on it.
Damage accumulates as a consequence of stress, allowing damage evolution functions to be defined [22,28]. A commonly used criterion is the Hayhurst criterion [22], which was designed to describe creep damage in ductile materials [30]. For uncorrelated diffuse damage a scalar representation of d may be enough, but when catastrophic failure occurs damage become highly correlated through crack propagation. This means the effective material will no longer be isotropic, and damage should be described by a higher-order tensor. Large deformations, often including significant rotations, will appear, which means that linear elasticity is no longer valid. Under such circumstances a mean-field representation of a material is usually no longer practical. Under some circumstances, particularly where ice is close to the pressure-melting point, healing can occur, reducing the level of damage. For glacier ice, few data are available to shed light on this important process.

Calving processes
Stresses sufficient to propagate fractures and detach ice-blocks can arise in several ways. Four basic calving scenarios were identified by Benn et al. [17]: (1) extension of the ice due to gradients in boundary friction; (2) tensile and shear stresses associated with force imbalance at ice fronts; (3) undercutting by subaqueous melting; and (4) torque arising from buoyant forces. Subsequent work has shown that this remains valid as a general framework, although our understanding of the processes associated with each scenario has increased substantially in the intervening decade. Here, we focus on recent research to illustrate the key processes and significance of these four fundamental calving scenarios.
The tractions at the basal and lateral boundaries of glaciers typically vary spatially, creating longitudinal and transverse stress gradients in the ice that can lead to fracture propagation (Figure 2(a)). In general, ice needs to be at or close to flotation to allow full-depth fracture penetration, because high water pressures in basal crevasses can offset the stabilizing effect of ice overburden pressure. Examples of calving due to longitudinal extension have been described by Walter et al. [31] at Columbia Glacier, Alaska. In 2008 and 2009, the retreating glacier terminus became buoyant, and the resulting loss of basal drag led to longitudinal stretching and the formation of full depth crevasses (rifts). These propagated over days to weeks, eventually releasing large bergs.
At the terminal cliffs of both glaciers and ice shelves, an imbalance always exists between outward-directed lithostatic pressure and backward-directed water pressure on the front. The resulting deviatoric stress increases with subaerial cliff height [32,33]. This stress is eccentrically distributed, being greatest at the base of the subaerial part of the terminal cliff. Tensile crevasses form at the glacier surface, whereas at depth ice can fail in shear if stresses are sufficiently high [34]. A stress threshold for shear failure of~1 MPa implies an upper bound on ice-cliff height of~110 m when both tensile and shear failure are considered, which closely matches observational data [35]. Buoyancy conditions imply that the total ice-front thickness (i.e. the subaerial plus submerged parts) cannot exceed~1000 m, and that catastrophic failure will occur in the absence of buttressing (e.g. support provided by mélange). This has major implications for the future stability of marine-based sectors of the Antarctic Ice Sheet, particularly Thwaites Glacier and adjacent catchments where the bed slopes inland to depths well in excess of 1 km [36]. If buttressing ice shelves are lost and ice fronts begin to retreat into deepening basins, Marine Ice Cliff Instability (MICI) could greatly accelerate rates of ice loss [37,38]. The factors controlling ice-cliff stability and consequent rates of retreat are not well understood, and the development of process-based models of MICI is an urgent priority.
Undercutting by subaqueous melting (melt-undercutting) encourages calving from tidewater and lake-calving glaciers by removing support from the subaerial part of the front (Figure 2(b)). This is now known to be a very important process wherever relatively warm water is brought into contact with glacier termini, including fjords in Alaska, Svalbard and Greenland [39][40][41][42][43]. Basal melting can also impact the stability of ice shelves, and is discussed separately below. Effective subaqueous melting requires tangential water motion across the ice front, which ensures efficient energy transfer. Wind-driven surface currents and wave action can encourage melt close the waterline in both marine and lacustrine environments [44][45][46]. In marine environments, rapid subaqueous melt is encouraged by the buoyant ascent of meltwater discharged from subglacial drainage systems, which can entrain warm ambient water and enhance heat flux to the ice [47,48]. Melt rates increase linearly with ambient water temperature and with meltwater discharge raised to a power between 0.3 and 0.75.
Melt-undercutting may be the dominant component of the total frontal ablation (frontal melt plus calving) on many tidewater glaciers [41,43]. The precise relationship between melt-undercutting and calving remains unclear, however. Depending on ice-front geometry and other factors, calved masses may be more or less extensive than the melted cavity. The former scenarioin which calving removes the entire overhanging part of the terminus plus additional icehas been dubbed the calving amplifier effect by O'Leary and Christoffersen [49], and analysis of stresses suggests that calving rates may be up to four times the subaqueous melt rate but decreasing with relative buoyancy [50]. However, the calving amplifier effect requires the presence of large melted cavities, and if frequent small calving events repeatedly remove parts of the overhang, the cavity may never become large enough for the amplifier to kick in. In this casewhich may be the norm for highly crevassed tidewater glacierslong-term calving rates may simply be directly paced by subaqueous melt rates [50,51].
Melt-undercutting also has non-local impacts on calving. For example, melt undercutting and enhanced calving in the vicinity of buoyant meltwater plumes can create embayments in the ice front, removing lateral support from the adjacent ice by breaking cross-glacier stress bridges [52]. This 'keystone effect' may thus allow local melting to affect calving losses across the whole terminus.
Under certain circumstances, the tongues of tidewater and lake-terminating glaciers can become super-buoyant. That is, the ice can be held below buoyant equilibrium and thus subjected to upward-directed buoyant forces. These forces can be relaxed by viscous flow and gradual uplift, or can lead to upward propagation of basal fractures and calving ( Figure 2(c,d)). Buoyancy driven calving was first described for lake-calving glaciers, where it can occur in response to ice-surface lowering through melt, or by lake-level rise [53,54]. More recently, it has been recognized that buoyancy driven calving is an important process on fast-flowing tidewater glaciers in Greenland [55][56][57][58]. As ice flows into deep water and approaches buoyancy, the eccentrically distributed stresses at the front cause forward bending. This initially prevents the ice from attaining buoyant equilibrium, and ice advance can increase the length of the super-buoyant tongue. This can become several hundreds of metres in length before buoyant forces are high enough to initiate calving, releasing very large bergs in dramatic events [55].
In nature, calving events can involve a mixture of the above scenarios, and different processes can occur in adjacent parts of the same glacier front. For example, Medrzycka et al. [59] found that on Rink Glacier, West Greenland, melt-undercutting is the main driver of calving near upwelling meltwater plumes at the glacier margins, whereas super-buoyancy dominates in the central part of the glacier.
Calving from ice shelves can occur by several mechanisms, ranging from the episodic release of huge tabular bergs to catastrophic disintegration events over periods of days to weeks [60]. For ice shelves in long-term steady state, infrequent tabular calving events can interrupt periods of ice front advance, forming multi-decadal cycles of ice-shelf growth and retreat [60,61]. Tabular bergs can be of immense size, up to~10 4 km 2 in area. Release of tabular bergs typically occurs following the slow propagation of full-depth fractures (rifts) across the shelf [62,63]. Rift propagation is largely driven by glaciological stresses, although rifts on thinner shelves (<~200 m) are also susceptible to oceanic swells [34]. A complicating factor in rift propagation is that shelf ice is typically inhomogeneous: ice of marine origin tends to be softer than meteoric ice, which inhibits propagation of rifts and can delay their progress across a shelf [64]. A recent well-studied example of rift-driven calving is the release of berg A68 from Larsen Ice Shelf C [65]. Following sporadic advance of a transverse rift since 2010, the 5,800 km 2 , 200 m thick berg was released between 10 and 12 July 2017. A similar large calving event preceded disintegration of Larsen B, and it is as yet unknown whether calving of A68 is part of a longterm cycle or a precursor of ice shelf break-up.
Disintegration events have affected numerous ice shelves in the northern Antarctic Peninsula in recent decades [66]. These appear to have been preconditioned by multiple factors, including thinning by surface and bottom melting, structural weakening, and loss of firn porosity by refreezing of surface meltwater [67][68][69][70]. Loss of firn porosity encourages accumulation of surface meltwater in ponds, which in turn can trigger rapid propagation of surface crevasses (hydrofracturing) [71]. Once initiated, fracturing of the shelf can spread rapidly in a chain reaction involving pond drainage, hydrostatic rebound, and ring-fracturing, culminating in the release and capsize of innumerable needle-shaped bergs [72,73]. The association of ice-shelf disintegration events in the Antarctic Peninsula with melt-pond formation emphasizes the need to understand ice surface hydrology, and the controls on storage and evacuation of surface melt, particularly with regard to the stability of the large ice shelves fringing the East and West Antarctic Ice Sheets [74,75].
Shelf disintegration, however, can also occur in colder environments with little or no evidence of surface melt [60]. These events involve intense crevassing and rifting along multiple lines of weakness, including basal crevasses originating near the grounding line, radial crevasses and fragmented shear margins [76] (Figure 3). Enhanced calving and fragmentation appears to be driven by complex sub-shelf processes, including pronounced and spatially varying melting in basal crevasses and subshelf channels, and interactions with damage evolution [77,78]. Prediction of the future evolution of Antarctic ice shelves requires that complex rifting and fragmentation processes can somehow be incorporated within prognostic models.

Calving models
Calving processes have been modelled using both analytical and numerical techniques. However, simulating calving events by setting and solving a set of equations is, in any practical sense, an impossible task. Thus, analytical approaches have chiefly been used to analyse stresses associated with specified glacier geometries, and to assess the implications for calving processes [32,33,49,58]. On the other hand, numerical techniques can be used to simulate, in detail, the entire calving process by representing ice as an assemblage of connected particles in discrete element models [79,80]. In models of ice sheet evolution, however, it is not possible to model every calving event in detail and calving losses must be represented using 'calving laws' based on empirical relationships or theoretical considerations. Here, we review explicit numerical models and the simple calving laws required to simulate ice sheet evolution.

Explicit models of fracture and calving
Ice fracturing processes can be simulated explicitly using discrete element models, which represent glaciers as assemblages of particles lined by breakable beams [79][80][81][82]. The most advanced model of this kind currently in use in glaciology is the Helsinki Discrete Element Model (HiDEM), which has been applied to a range of idealized and real-world problems [50,83,84]. The model code is freely available at: https://github.com/joea todd/HiDEM. Initial boundary conditions in HiDEM are set by the ice geometry and the depth of the adjacent water body. During a simulation, the dynamics of the ice is computed using a discrete version of Newton's equations of motion, using inelastic potentials to calculate the interactions of individual particles and beams. As the ice deforms elastically under its own weight, stresses on the beams increase. If stress reaches a failure threshold the beam breaks and particles (or aggregates) become disconnected but continue to interact as long as they are in contact. In a violent fragmentation event, a large number of interacting cracks and shock waves create a very complicated process that rapidly dissipates energy until the fragmenting body can find a new equilibrium and come to rest. The microscopic-scale interactions of beams and particles remain rather simple, but high complexity may arise from the large-scale dynamics of many coupled and nonlinearly interacting degrees of freedom.
In HiDEM, this approach is not limited to elastic-brittle materials, but can also encompass all kinds of creep, plasticity and viscosity [82]. All possible combinations can be modelled in terms of time, strain, strain-rate, or stress dependent interactions between particles. Small-scale deformations are handled by the interactions potentials, and large-scale (viscous) deformations by the destruction and reformation of potentials between particles that move past each other. However, for calving events, which typically occur on timescales of seconds to minutes, viscous processes can usually be neglected and it is sufficient to consider only brittle failure. In this case, breakage of beams is irreversible and consequent interactions are only solid body collisions.
The equations of motion can be written as: where m i , c i , d ij are diagonal mass, drag, and dissipation matrices, and k ij are stiffness matrices of connected elements. The first term on the left-hand side represents inertial forces, the second represents viscous drag, the third dissipative forces in collisions between elements, the fourth elastic forces between elements, and the fifth encompasses forces such as gravitation, interaction with the bed and buoyancy. As a 2D example we can write , the elements of d ij determine the energy dissipated in relative motion of elements in contact, and k ij is the stiffness matrix that when multiplied with the vector x ij ! ¼ x i ; y i ; θ i ; x j ; y j ; θ j À Á gives the elastic forces. The drag on a single isolated particle can be modelled fairly accurately as function of relative ice-water velocity. However, more accurate computations of the coupled interaction between hydrodynamics and agglomerates of particles would increase the computational demands tremendously, and are ignored in the current HiDEM model. Equation 1 can be written in a discrete form and time can be integrated forward once initial conditions (i.e. the initial vectors ðx i ; y i ; θ i Þ) are determined. It is something of an art form to design the exact setup of the discrete integration scheme to obtain a suitable compromise between numerical efficiency and accuracy. This has been investigated extensively for molecular dynamics, which has the same type of algorithm [85]. The final piece needed to set up a numerical model of fragmentation is to determine and implement a fracture threshold. This can be defined in terms of either a fracture stress, strain or elastic energy threshold. It is possible to implement fracture by reshaping the model to remove broken connections, but numerically more efficient is to make the elements of k ij dependent on stress/strain and elastic energy history. For example, fracture may be implemented by setting k ij ¼ 0 as soon as a fracture threshold is exceeded.

Fragment size distributions, FSDs
Fragmentation of elastic-brittle material is a process with very high entropymeaning that the number of possible ways a brittle material can break, each with more or less the same probability, is very high. This does not however mean that fragmentation is a process without repeatable characteristics, which can be used to compare the results of HiDEM simulations with observations. Not all calving events, of course, involve continuous and pervasive fragmentation. For example, tabular bergs released from floating ice tongues or shelves typically remain intact with minimal breakup following initial rift propagation. Most styles of calving do, however, involve complex fragmentation histories during and following berg release [83]. Failure of overhanging ice fronts following melt undercutting typically results in fragmentation of failed masses into many thousands of pieces [86], while capsize of deep, narrow bergs typically involves multiple interactions with adjacent ice masses, triggering further fragmentation [72]. Thus, the fragment size distribution (FSD) is often a very good measure for quantifying and categorizing fragmentation processes, yielding insights into calving mechanisms.
The FSD can be written as follows: n s ð Þ / f s ð Þds, where f(s)ds stands for density of fragments of size s in an interval ds. There are a few standard forms of the FSD that are relevant for calving. (1) Calving can be described as a Poisson process, which results in an exponential FSD, n s ð Þ / exp Às=s0 ð Þ, where s 0 is a crack density constant. (2) Calving debris and icebergs can be formed by branching and merging cracks, and n s ð Þ / s À 2DÀ1 ð Þ =D , where D is dimension. It is reasonable to set D = 2 for large ice sheets or shelves over length scales much larger than their thickness, and D = 3 for glaciers confined to valleys or fjords. 3) Where fracture is dominated by grinding, for example in mature shear zones, the FSD is typically n s ð Þ / s Àα . In 3D, typical grinding exponents are in the range / ,1:8 À 3.5 [80].
A simple 2D calving scenario is shown in Figure 4, representing a~100 m high, grounded terminus of a tidewater glacier resting on a uniform slope. The mechanical interaction of ice particles with water is represented by buoyancy and viscous drag only. The glacier terminus fractures as a result of its own weight, and initial crack propagation creates an FSD that is a result of a merging-branching process. As calving proceeds there is considerable grinding in shear bands, and consequently the FSD computed for late in the simulation is typical of grinding processes.

'Calving laws'
'Calving laws' are not physical laws based on rigorous and repeatable experiments, in the sense of Boyle's Law. Rather, they are simple empirical, theoretical, or hybrid functions intended to calculate bulk calving losses in numerical ice sheet models [87]. Calving laws can be divided into two broad categories: functions that predict calving rate and those that predict calving position. In principle, all calving laws could be recast in either form depending on timescale and other factors (J. Bassis, personal communication), although in practice the two approaches have been developed separately.
The earliest calving laws were rate laws, based on observed correlations between calving losses and environmental variables such as water depth at the terminus [88,89]. Correlations between calving rates and water depths are often remarkably strong within local populations of glaciers, but differ widely between regions, hinting at complex causal relationships [17]. Alley et al. [90] proposed an ice-shelf calving rate law based on correlation with longitudinal strain rates. This approach was extended to strain in two dimensions by Levermann et al. [91], who proposed the following relation for calving rate C: where e 1 and e 2 are orthogonal horizontal principal strain rates, and K 2 is an empirically determined proportionality factor. Morlighem et al. [92] took a similar approach to calculating calving from outlet glaciers. The physical reasons for the correlations underlying rate laws are often unclear. Van der Veen [93] argued that because (1) calving rate is defined as the difference between ice velocity and the rate of change in terminus position, and (2) ice velocity is typically much greater than rates of terminus position change, then correlations between e.g. water depth and calving rate may actually reflect causal relationships with ice velocity rather than calving. Similar arguments may apply to calving rate laws based on strain rates, rendering them of little use in non-steady state conditions.
Calving position laws predict the location of the calving front at any given time from the geometry and/or state of stress in the glacier. This approach was pioneered by van der Veen [93], who observed that the calving front of Columbia Glacier was typically located where ice thickness was 50 m greater than that required for flotation. Van der Veen proposed a height above buoyancy calving law, which was generalized by Vieli et al. [94] to define the position of the calving front as the point where ice thickness h is: where D W is water depth and q is a small fraction. A limitation of this calving law is that it cannot permit floating ice tongues or ice shelves, although these rarely form on most glaciers outside of Antarctica. A physical basis for the height above buoyancy calving law was proposed by Bassis and Walker [34], and is discussed below.
A calving position law was proposed by Benn et al. [17,87], based on the penetration of crevasses opened in response to tensile stresses. This crevasse depth calving law locates the glacier margin where surface crevasses penetrate to the waterline, or where surface and basal crevasses penetrate through the full thickness of the glacier [95]. The full-thickness criterion is probably more realistic for most calving scenarios. Early implementations of the calving law calculated crevasse depth from the Nye zero stress model, which assumes that a field of closely spaced crevasses will penetrate to a depth where the tensile stress in the along-flow direction is exactly balanced by the lithostatic pressure [96]. This has been generalized to include tensile stresses in all directions, such that crevasses are assumed to exist wherever the largest principal component of the Cauchy stress tensor, σ 1 > 0 [97,98]. The effect of water pressure P w in surface or basal crevasses can be included by setting σ 1 + P w > 0. Bassis and Walker [34] and Ma et al. [35] modified the crevasse depth calving law by allowing ice failure wherever shear stresses exceed a yield value, in addition to tensile crevasse penetration calculated from the Nye zero stress model. This predicts an upper bound on terminus ice thickness for given water depths, as a consequence of the large shear stresses associated with high, unsupported ice cliffs. Conversely, a lower bound on ice thickness is predicted from the penetration of surface and basal crevasses where ice approaches buoyancy ( Figure 5). Taken together, these bounds provide a physical justification for height-above-buoyancy calving laws, and an explanation why floating tongues cannot form on highly fractured outlet glaciers.
The presence of water in surface crevasses allows them to penetrate through greater thicknesses of ice, because water pressure counteracts the effect of ice overburden pressure on crevasse closure [17][18][19]. This makes crevasse-depth calving laws ideal for modelling ice-shelf disintegration through hydrofracturing, provided water depths can be calculated adequately [37]. However, the use of crevasse water depth as a tuning parameter in calving models has tended to over-emphasize glacier sensitivity to atmospheric forcing. This has been exacerbated by the fact that crevasse depth calving laws are insensitive to melt-undercutting, which leads to under-estimation of the impact of oceanographic forcing on calving losses [99,100].
An alternative approach to predicting the location of calving events is based on parameterization of damage. This has been applied to a wide range of problems, including the stability of hanging glaciers in high mountain terrain [22] and weakening of ice shelves [23]. Krug et al. [29] modelled tidewater glacier calving by combining a damage parameterization and a linear elastic fracture mechanics method of calculating crevasse depth. Most recently, Mercenier et al. [101] developed a calving rate law based on a stress threshold and a damage evolution function, tuned to fit data from Arctic glaciers.

Continuum models
Various calving laws have been implemented in continuum models to simulate the past and future behaviour of glaciers and ice sheets. An important issue is the representation of stress in these models. For computational reasons, most models of ice sheet evolution employ reduced stress formulations, such as the shallow ice approximation, that includes vertical shear but neglects longitudinal and transverse stresses [102], and the shallow shelf approximation, that includes longitudinal and transverse stresses but neglects vertical shear [103]. Models of the Antarctic Ice Sheet typically employ the former for ice-sheet interiors where basal drag accounts for most of the resisting forces, and the latter for ice stream and ice shelves where basal drag is small or zero [38,104]. Due to the intricate geometry of the Greenland Ice Sheet, calving model studies have tended to focus on individual outlet glaciers. Most studies have used reduced stress formulations implemented in 1D (flowline) [105,106] or 2D (map view) [92]. Recently, calving routines have been developed for the full Stokes continuum model Elmer/Ice in two and three dimensions [97,98], opening up new possibilities to explore the full range of calving processes in continuum models. However, it is currently impractical to solve the full Stokes equations for large (> 10 3 km) model domains over years of model time. Consequently, reduced stress models will continue to play an important role in the foreseeable future, particularly in view of the need for fast models for ensemble experiments [107]. Development of robust calving laws suitable for these platforms remains an important challenge.
For Antarctica, a strain rate-calving rate law implemented in the Potsdam Parallel Ice Sheet Model has been used to predict ice sheet response to greenhouse gas emissions [91,104,108]. These studies indicate that retreat of marine-based sectors of the ice sheet will likely contribute several metres to global sea level in the coming centuries to millennia. A different approach to the same problem was taken by Pollard et al. [37], who used a crevasse-depth calving law to simulate hydrofracturing of ice shelves, and an ice-cliff stability criterion together with a heuristic rate function to model ice-cliff instability. When tuned to fit model output to palaeo-sea-level targets, this model predicts very fast deglaciation of West Antarctica under high emission scenarios, with sea-level rise of~1 m by 2100 [38]. The calving laws used in these studies must be tuned to observations, introducing large uncertainties when extrapolated into the future. Clearly, narrowing of these uncertainties is a major priority.
For Greenland outlet glaciers, height-above-buoyancy and crevassedepth calving laws have been implemented in reduced stress ice-flowline models to explore a range of problems, including past, present and possible future ice-front fluctuations [94,105,106]. In addition, Enderlin et al. [109,110] used a similar model and synthetic geometries to explore the sensitivity of Greenland outlet glaciers to input uncertainties, including bed topography. Morlighem et al. [92] used a strain-rate-calving rate law in the Ice Sheet System Model (ISSM) to investigate the impact of submarine melting on Store Glacier. All of these models use vertically integrated stress fields, and assume buoyant equilibrium for floating parts of glacier termini. Consequently, they cannot explicitly represent melt-undercutting and buoyancy driven calving processes. These shortcomings are overcome in the full Stokes model Elmer/Ice, illustrated by studies of Store Glacier in 2D and 3D [97,98]. The results of the 3D model closely match observed spatial and seasonal patterns of calving, including detachment of large tabular bergs from buoyant parts of the terminus at the onset of the melt season ( Figure 6).

Discrete element models
Because they explicitly simulate fracture and calving processes, discrete element models can play an important role in efforts to understand calving processes and their relationship with environmental conditions. An early example is the use of a 2D discrete element model by Bassis and Jacobs [80] to investigate the influence of glacier geometry on calving style. More recently, discrete element models have been used to investigate the effects of mélange backstress on calving by Robel [111] and Burton et al. [26]. The latter study is particularly useful, as it combines model simulations with field data and laboratory experiments to quantify buttressing forces. The Helsinki Discrete Element Model (HiDEM) has been used to study a wide range of problems, using both synthetic and real glacier geometries. Some of these applications are discussed in the following sections. Figure 6. Calving on store glacier, Greenland, modelled in Elmer/Ice, including a large tabular berg on the near-side of the glacier front. This event immediately followed mélange break-up, and is similar in size, location and timing to observed events. Image by Joe Todd.

Synthetic geometries
Experiments with synthetic geometries show that many of the calving processes discussed in Section 2.2 emerge spontaneously from HiDEM under contrasting boundary conditions. In particular, variations in damage, basal drag, undercutting and buoyancy can produce a wide range of observed calving styles, including ice-cliff instability, buoyant uplift and overhang collapse following undercutting [50]. Figure 7 shows a matrix of model runs for varying undercut lengths U L , and water depth D W relative to that required for ice flotation D F . Each panel shows mean calving lengths modelled in HiDEM (numbers and vertical green lines) and the effective principal stress (EPS: σ 1 + P w ) calculated in Elmer/Ice. Note that the stress fields relate to intact ice at the start of the HiDEM simulations, and that stresses change significantly as cracks propagate. The results reveal the following patterns: (1) calving length increases with the height of the subaerial ice cliff; (2) calving length increases with undercut length; (3) the effect of undercutting decreases as relative buoyancy increases; and (4) super-buoyancy induces calving. Additional experiments (not shown) indicate that down-flow reduction in basal drag (e.g. associated with a waterpressure-dependent sliding law) can trigger calving by full-depth crevasse penetration.
These results indicate that calving processes form a continuous spectrum, depending on combinations of boundary conditions. This suite of behaviours arises from the iteration of a few simple rules governing the response of particles and inter-particle bonds to gravitational and buoyant forces. Thus, the complexity of the 'calving problem' arises from the wide range of possible glacier geometries and boundary conditions rather than the inherent complexity of the underlying processes.
The combination of Elmer/Ice and HiDEM allows calving processes to be analysed in considerable detail. For example, Benn et al. [50] investigated Greenland-style buoyant calving using Elmer/Ice to simulate the flow of ice into deep water, then importing glacier geometry for selected timesteps into HiDEM to simulate calving events. Calving of the super-buoyant tongue was associated with high-tensile stress at the ice base near the ungrounding point (Eϒϒ > 0.9 MPa) and compressive stress at the ice surface. This stress pattern is diagnostic of the torque created by large buoyant forces. The zero stress crevasse criterion, which permits crevassed ice where EPS > 0, predicts the penetration of basal crevasses only about half-way through the ice under these conditions, not enough to trigger calving according to the crevasse-depth calving law. This is because the Elmer/Ice stress field is a snapshot of conditions in intact material, and does not consider the effect that fracturing will have on the remaining ice. HiDEM, of course, simulates this for every time-step as fractures propagate. With this limitation in mind, the Elmer/Ice snapshots nevertheless provide valuable insights into the conditions required for calving, allowing identification of diagnostic stress patterns, and guiding the design of improved calving laws.

Real-world geometries
HiDEM can be used to study real-world fracturing and calving processes in considerable detail. Required input data include digital elevation models of the ice surface, the bed (i.e. rock or sediment) on which the ice rests, and (for floating ice shelves or tongues) the ice base. To improve simulation results it is also useful to use inversion techniques to estimate the effective friction between the moving ice and the bed, utilizing observed surface velocities [84]. The necessary data (ice velocities, surface and bed DEMs) are increasingly available for many regions. Here we discuss several examples that illustrate the capabilities of this approach.
The first example illustrates the use of HiDEM to simulate surface fracturing on the Basin 3 Glacier that drains the Austfonna ice cap in Svalbard. This glacier underwent a major surge following autumn 2012, and Gong et al. [112] used HiDEM in conjunction with Elmer/Ice and a simple hydrological model to test the hypothesis that surface water routed to the bed played a major role in surge activation and propagation. Modelled crevasse patterns compare well with observations (Figure 8), and when these are used to define water input points, the hydrological model predicts water flow towards regions of observed fast flow.
The second example focuses on Kronebreen, a fast-flowing (~3 m day −1 ) tidewater glacier in Svalbard, that has been the subject of detailed studies of ice dynamics, hydrology and calving [42,84,113,114]. At Kronebreen,   calving rates show a very strong correlation with fjord water temperature over seasonal timescales, and a correlation with ice velocities at times when velocity is high [42]. Vallot et al. [84] used HiDEM in conjunction with other models to investigate the influence of melt-undercutting and basal friction on calving. The results showed that undercutting was necessary to replicate observed calving rates during the summer melt season, and that low basal friction resulted in increased strain rates and calving near the end of the summer (Figure 9). These results explain the observed correlations, and show how HiDEM can elucidate the processes underlying observed glacier behaviour. The third example we present here focuses on fracturing of the Totten ice shelf in eastern Antarctica [115]. Figure 10 shows a comparison of observed fracture patterns with HiDEM results using the present-day shelf thickness and location of grounded zones as boundary conditions. The model captures well the tensile cracks that appear as a result of terminus spreading. In contrast, many of the transverse fractures (basal crevasses and rifts) are missing in the computed image. These fractures originate far upstream. A HiDEM computation of the grounding region ( Figure 10(d)) reveals that these fractures emerge as tensile fractures near the ice-base mainly over re-grounding islands. In nature, these basal fractures are transported down-stream by ice flow, and widen and grow by melt and stretching until they reach the terminus region where they precondition calving events (Figure 10(c)). These temporal effects are not currently included in the HiDEM simulation, which simulates elastic-brittle relaxation of the ice from an initial geometry. Realistic modelling of calving of the shelf, and investigation of future ice-shelf stability, must therefore take its long-term evolution into account.

Conclusions and future prospects
There have been major advances in both observational and modelling capability in recent years. Modelling tools such as HiDEM and Elmer/Ice in particular have allowed detailed analysis of calving events, and a comprehensive understanding of calving processes is beginning to emerge. Important problems, however, remain to be solved.
Many details of calving processes need to be worked out, requiring targeted observations and detailed modelling studies. Key problems include the following: (1) quantifying the links between plume dynamics, melt-undercutting and calving; (2) understanding divergent glacier behaviour under near-buoyant and super-buoyant conditions; (3) measuring and modelling the properties of damaged ice, including mélange; (4) quantifying the factors that precondition and trigger ice shelf collapse; and (5) understanding ice-cliff instability mechanisms, including the controls on rates of ice-cliff retreat. For these and other problems, the central challenge is to develop a detailed process-based understanding of system physics that can be translated into predictive capability.
Defining robust, flexible calving laws is a major priority. The strikingly consistent patterns observed in nature (e.g. relationships between ice-cliff heights and buoyancy conditions; Figure 4) suggest that simple calving laws can indeed reliably predict ice-front locations, despite the complexity of calving processes. But the optimum form of such calving laws still remains an open question.
It is clear that calving glaciers and ice sheets make a large contribution to sea-level rise, but much uncertainty remains about future ice sheet response to alternative carbon futures. The foremost source of uncertainty is the West Antarctic Ice Sheet, parts of which appear to have already embarked on a process of irreversible retreat [36]. Future rates of ice loss will depend critically on ice shelf response to atmospheric and oceanic forcing, and what happens once buttressing ice shelves are lost and ice-cliff instability kicks in [37,116]. Potentially, this will entail behaviour well beyond the observed range, involving complex interactions between fracturing, flow of damaged ice, calving, and marine processes. Meeting the challenge of predicting this behaviour will require an integrated approach, drawing upon all available and developing tools, from high-resolution explicit models to the simple calving laws required for long-term simulations of ice-sheet evolution.