Introduction
The largest source of uncertainty in projections of future sea-level rise is the response of the Antarctic ice sheet to a warming climate (e.g. Pachauri and others, Reference Pachauri2014; Edwards and others, Reference Edwards2021). Unlike in Greenland, the majority of Antarctica's ocean-terminating glaciers extend into the ocean to form floating ice shelves (Cuffey and Paterson, Reference Cuffey and Paterson2010). Horizontal stress gradients in these fringing shelves, associated with pinning points and embayment walls, reduce the ice flux across the grounding line by buttressing the upstream glaciers (e.g. Dupont and Alley, Reference Dupont and Alley2005; Goldberg and others, Reference Goldberg, Holland and Schoof2009; Gudmundsson, Reference Gudmundsson2013; Pegler, Reference Pegler2018).
Antarctic ice shelves lose mass to the ocean primarily through ocean-driven submarine melting and iceberg calving (e.g. Liu and others, Reference Liu2015). Although complex, the thermodynamics associated with large-scale ice-shelf melt are relatively well understood (Jenkins and Holland, Reference Jenkins and Holland2007). However, both the small- and large-scale physical rifting and fracturing processes controlling iceberg calving have proven more difficult to understand and model (Benn and others, Reference Benn, Warren and Mottram2007b). Calving changes the geometry of ice shelves through the production of icebergs, which can separate the ice from pinning points or embayment walls. These detachments reduce the buttressing of ice shelves, increasing the ice flux across their grounding lines. Because much of the West Antarctic ice sheet is grounded below sea level and ice flux across the grounding line depends strongly on the ice thickness there, the loss of buttressing affects the stability of the grounding line itself (Schoof, Reference Schoof2007; Joughin and Alley, Reference Joughin and Alley2011; Gudmundsson, Reference Gudmundsson2013; Pegler, Reference Pegler2018; Martin and others, Reference Martin, Cornford and Payne2019).
Early attempts to model iceberg calving focused on calving from grounded margins and suggested a proportionality between water depth or ice thickness and calving rate (e.g. Brown and others, Reference Brown, Meier and Post1982; Warren, Reference Warren1992). Benn and others (Reference Benn, Hulton and Mottram2007a) and Amundson and Truffer (Reference Amundson and Truffer2010) later extended the ice-thickness-criterion calving model to support floating margins, and others have suggested heuristic relationships based on ice temperature (Reeh, Reference Reeh1968) or strain rate (Alley and others, Reference Alley2008; Levermann and others, Reference Levermann2012), but these are unlikely to be broadly applicable, both geographically and under a changing climate. Nye (Reference Nye1957) developed a process-based ‘zero-stress model’ that calculates the depth to which surface crevasses penetrate, and Weertman (Reference Weertman1973) derived an equivalent formulation for ice-shelf basal crevasses. However, the latter predicts that basal crevasses will rarely penetrate the entire ice thickness of freely floating ice shelves unless surface depressions become filled with water (van der Veen, Reference van der Veen1998a, Reference van der Veen1998b).
In areas with sufficiently large surface melt rates, networks of supraglacial lakes can form on an ice shelf. When these lakes fill crevasses, they drive downward penetration by applying an additional tensile stress at the crevasse tip (Weertman, Reference Weertman1973; van der Veen, Reference van der Veen1998b), a process known as hydrofracture. Combined with the elastic response of the ice shelf to the meltwater load, these lake networks can connect systems of crevasses throughout the ice. In this case, drainage of one supraglacial lake may cause a chain reaction among adjacent lakes that destabilizes the entire ice shelf, causing it to collapse (Banwell and others, Reference Banwell, MacAyeal and Sergienko2013). The Larsen B Ice Shelf underwent such a meltwater-induced disintegration in 2002, and Borstad and others (Reference Borstad2012) used the theory of continuum damage mechanics in a large-scale ice-dynamics model as a proxy for the ‘health’ of the ice shelf prior to the collapse (Borstad and others, Reference Borstad2016). Others have used damage mechanics (e.g. Pralong and Funk, Reference Pralong and Funk2005; Duddu and Waisman, Reference Duddu and Waisman2012; Sun and others, Reference Sun, Cornford, Moore, Gladstone and Zhao2017) and linear elastic fracture mechanics (e.g. Kenneally and Hughes, Reference Kenneally and Hughes2002) to model iceberg calving, but the computational costs are challenging when incorporating these complex, process-scale models within large-scale ice-sheet models in order to project ice-sheet evolution over centennial and millennial timescales.
Recently, Bassis and Ma (Reference Bassis and Ma2015) developed a process-based, crevasse-depth model that provides a link between the propagation of fractures through an ice shelf and basal melt, the existence of which has been suggested by previous studies (e.g. Liu and others, Reference Liu2015). Using a linear stability analysis, they deduced an evolution equation that represents crevasse depth as a continuum variable and governs the deepening and widening of crevasses in time. Here, we implement this crevasse evolution model as a passive tracer in the BISICLES ice flow model (Cornford and others, Reference Cornford2013) and apply it to a number of idealized ice-shelf geometries to compare it to analytic test cases and to demonstrate its capability in modeling ice-shelf extent in different environments with widely varying characteristics.
Model description
Damage mechanics
Our model represents ice-shelf crevasse penetration depths by the effective ‘damage’, a continuum variable we define as
where δh is the total crevasse depth and h is the local ice thickness, as shown schematically in Fig. 1. This damage variable is defined over the range [0, 1], where r = 0 in the absence of crevasses and r = 1 when crevasses extend through the full thickness of the ice shelf. This is similar to the model proposed by Sun and others (Reference Sun, Cornford, Moore, Gladstone and Zhao2017), which incorporates crevasse depth directly into the ice rheology. However, unlike other reported models, which use ad hoc definitions of damage laws, the crevasse depths in our model evolve based on a pseudo-plastic necking instability (Bassis and Ma, Reference Bassis and Ma2015). In contrast, Borstad and others (Reference Borstad2016) define damage through a viscosity weakening factor, which they find by way of inversion with assumed ice temperatures and recompute instantaneously based on computed strain rates, whereas Albrecht and Levermann (Reference Albrecht and Levermann2012) define damage as a fracture density field with a specified fracture growth rate that depends only on the strain and an empirical parameter to be fitted.
Fig. 1. Schematic diagram of the model system showing a longitudinal cross section of an ice shelf. The ice flux at the grounding line is held constant in time and directed along the x-direction. Damage r is defined as the local ratio between the crevasse depth δh and total ice thickness h.
We describe the necking instability using the results of a linear stability analysis, in which the strain-rate-weakening nature of ice magnifies large-scale ice-shelf thinning and local ductile deformation. In the long-wavelength limit appropriate for depth-integrated, two-dimensional (2-D), dynamic approximations to full-Stokes ice-sheet models, where the widths of basal crevasses are large compared to the ice thickness, damage evolves according to (Bassis and Ma, Reference Bassis and Ma2015)
where $\dot {\varepsilon }_{1}$
is the depth-averaged, largest principal strain rate and $\dot {m}$
is the total melt rate (which may be negative for accumulation/accretion). Although $\dot {m}$
may include contributions from both surface and basal melting, here we focus exclusively on the basal contribution, leaving exploration of surface crevasses and the process of hydrofracture for future work. We have also assumed that the melt rate within a crevasse is the same as the ambient melt rate of the basal surface of the ice, though observations have shown everything from accumulation of marine ice within crevasses (e.g. Fricker and others, Reference Fricker, Popov, Allison and Young2001) to enhanced melt (e.g. Dutrieux and others, Reference Dutrieux2014).
The dimensionless number S 0 is a ratio between hydrostatic pressure and the largest principal deviatoric stress within the ice shelf:
where τ 1 is the largest principal deviatoric stress. We use ρ i = 910 kg m−3 as the meteoric ice density, ρ w = 1028 kg m−3 as the density of sea water and g = 9.81 m s−2 as the acceleration due to gravity. We assume that the entire thickness h has the density of meteoric ice and therefore ignore the presence of firn. The star superscript on $n^\ast$
denotes the usual parameter n from Glen's flow law, adjusted to include the ratio of the principal horizontal strain rates $\alpha = \dot {\varepsilon }_{2} / \dot {\varepsilon }_{1}$
:
Equation (2) assumes that the perturbation to the ice shelf extends uniformly through the thickness such that there exists a relationship, α, between the principal strain rates (Bassis and Ma, Reference Bassis and Ma2015). The evolution equation is consistent with the damage formulation proposed by Albrecht and Levermann (Reference Albrecht and Levermann2012), with the addition of a dependence on the stress state and the basal melt rate. Although this model could be used in conjunction with, for example, the eigencalving model of Levermann and others (Reference Levermann2012) to calculate the average probability of calving at the terminus (Albrecht and Levermann, Reference Albrecht and Levermann2012), here we are primarily interested in the location where crevasses are predicted to penetrate the full ice thickness (i.e. where r = 1), which coincides with what we refer to below as the ‘fully-damaged terminus’.
In the linear stability analysis, when S 0 < 1, the tensile stress of the ice shelf pulls crevasses apart and allows them to deepen. By contrast, when S 0 > 1, the gravitational restoring force causes the ice to flow into the depressions associated with crevasses, healing them over time even in the absence of negative strain rates and accretion, which are the sole drivers of crevasse healing in other damage models (Pralong and Funk, Reference Pralong and Funk2005; Borstad and others, Reference Borstad2016). This viscous healing can nevertheless be overcome by a sufficiently large melt rate that deepens the crevasse faster than it can close. Ice shelves usually have S 0 > 2 across most of their horizontal extents, with S 0 < 1 occurring in small shear bands with high tensile stresses or where the ice is thin (Bassis and Ma, Reference Bassis and Ma2015). Below, we show that freely floating ice tongues have S 0 = 2. Not surprisingly, the gravitationally driven flow of ice tongues results in crevasses that close in the absence of thinning by melting.
Nothing in Eqn (2) explicitly prevents damage from exceeding r = 1, unlike other damage models (e.g. Albrecht and Levermann, Reference Albrecht and Levermann2012; Borstad and others, Reference Borstad2016). However, as we have defined r as the ratio of crevasse depth to ice thickness (Eqn (1)), we use r = 1 as a physical upper bound, consistent with the applications of damage mechanics. As damage approaches this upper limit, the perturbative approach used to derive Eqn (2) breaks down, similar to most other damage mechanics, which are technically valid only for small damage (Keller and Hutter, Reference Keller and Hutter2014).
Our model cannot form new crevasses to increase damage when r = 0 in Eqn (2) as the evolution results from the instability of already present perturbations. We therefore take the ice to always be densely crevassed along the basal surface to at least the depth where the tensile stress pulling them open balances the cryostatic pressure pushing them closed, also called the Nye zero-stress model (Nye, Reference Nye1957; Weertman, Reference Weertman1973; Jezek, Reference Jezek1984; Nick and others, Reference Nick, van der Veen, Vieli and Benn2010):
which assumes that ice has small tensile strength. With finite fracture strength, the initial value of, and lower-bound to, damage could be smaller (Weertman, Reference Weertman1977). For a freely floating ice tongue, r N ~ 0.44 everywhere, as we show below. For more complicated ice shelves r N is spatially varying: reduced relative to the freely floating ice tongue in the presence of buttressing and increased where ice is spreading. The damage mechanics of Sun and others (Reference Sun, Cornford, Moore, Gladstone and Zhao2017) advects this Nye zero-stress damage, which increases as it advects into thinner ice, but does not account for the growth of fractures in a manner akin to our shear-thinning necking to deepen an already present crevasse (the first term in the brackets on the right hand side of Eqn (2)).
In most scalar damage mechanics models, the damage affects the stress field of the ice by reducing the effective cross-sectional area over which forces act. Thus, damage is associated with a factor (1 − r) −1 multiplying either the deviatoric stress, resulting in a decreased viscosity (Borstad and others, Reference Borstad2016; Sun and others, Reference Sun, Cornford, Moore, Gladstone and Zhao2017), or the full stress tensor to affect buoyancy and allow hydrofracturing (Pralong and Funk, Reference Pralong and Funk2005; Duddu and Waisman, Reference Duddu and Waisman2012; Mobasher and others, Reference Mobasher, Duddu, Bassis and Waisman2016). Our model incorporates damage in an effective stress (including the hydrostatic component) in determining the growth rate of damage (Eqn (2)) and can be extended to affect the ice rheology in this way (e.g. as in the model of Sun and others, Reference Sun, Cornford, Moore, Gladstone and Zhao2017). In this work, however, we neglect this higher-order coupling and damage acts as a passive tracer. A future extension of this damage model will incorporate both the physical mechanisms for crevasse evolution presented here as well as the effects of damage on ice rheology.
Numerical method for damage
We implement damage evolving and advecting with the ice flow velocity u using
with dr/dt given by Eqn (2), in the BISICLES finite-volume, adaptive-mesh-refinement ice flow model, which solves a modified form of the L1L2 ice equations from Schoof and Hindmarsh (Reference Schoof and Hindmarsh2010). As BISICLES expresses the laws of motion in vertically integrated, conservative form, we write the damage evolution as
for some source σ and vertically integrated velocity $h\bar {{\bf u}} = \int _0^h{\bf u}{\rm d}z$
. To determine the form of σ, we integrate Eqn (6) vertically, noting that damage has no depth dependence (i.e. r = r(x, y) for horizontal coordinates x, y and vertical coordinate z)
We apply the product rule to obtain
The terms in the square bracket constitute the left-hand side horizontal transport equation (see, e.g. Cornford and others, Reference Cornford2013, Eqn (7))
The melt rate term is negative because we are defining positive melt rates to remove mass. Substituting Eqn (10) into Eqn (9) and changing sides gives:
Finally, comparing Eqns (8) and (11) we find that the vertically integrated, conservative force to be
which is now written in terms of the damage evolution expression of dr/dt in Eqn (2).
Incorporating the source term, we update the damage field using a piecewise parabolic method-based predictor-corrector method (see Cornford and others, Reference Cornford2013, Section 3.5):
for each cell in the domain, and the time step dt is calculated from a Courant–Friedrichs–Lewy condition. We constrain the damage to the range r ∈ [r N(x, y, t), 1] at all times, where r N is given by Eqn (5), so that crevasses always penetrate to at least the depth at which the tensile and compressive stresses balance.
Test cases
To assess model performance we first compute the steady-state damage of idealized representations of ice tongues and ice shelves. We specify the initial ice-shelf geometry, the sub-shelf melt rate and assume isothermal ice. The melt rate applied in all cases is spatially uniform over the whole ice shelf and constant in time.
Analytic solution for ice tongues
We begin by examining geometries that correspond to the Erebus Glacier Tongue and to the unembayed portion of the Drygalski Ice Tongue. We focus exclusively on the unembayed portion of the Drygalski Ice Tongue here to avoid the complicating effects of shear margins, which would make the system inherently 2-D, as we explore further in the ‘Idealized ice shelves’ Section.
Ice tongues have proven challenging to represent in numerical models because they range in length from tens to hundreds of kilometers. This variability makes them difficult to recreate in large-scale ice-sheet models, which often choose to ignore or remove ice tongues (see, e.g. Martin and others, Reference Martin2011), or instead predict them in places where they do not exist (see, e.g. Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012, their Fig. 8, which predicts a long ice tongue from Sermeq Kujalleq in Greenland). With our two examples, we demonstrate how damage evolution provides a plausible explanation for, and mechanism for modeling, their characteristics. Another reason for starting with ice tongues is that they lack the complications associated with buttressing and shear zones and can thus be approximated as 1-D by applying free-slip, lateral boundary conditions. This simplification allows us to generate an analytic, steady-state solution for the damage with which to verify and validate the damage evolution simulated by BISICLES.
With no bottom or lateral stresses, the depth-averaged driving stress for an ice tongue is, following Weertman (Reference Weertman1957):
where $h^\ast ( x)$
is the ice thickness at location x downstream of the grounding line and the star superscript represents profiles of the 1-D, floating ice tongue in steady state. Through Glen's flow law, Eqn (14) results in a strain rate of
where we define the parameter C as
that depends on the Glen's flow law parameter A, which is a function of temperature. Substituting the strain in Eqn (15) into the continuity Eqn (10) and integrating, van der Veen (Reference van der Veen2013, Eqn 5.80) derived an analytic solution for the thickness of an ice tongue experiencing a positive, uniform melt rate $\dot {m}$
of
The constant flux boundary thickness and velocity are given by h 0 and u 0, respectively, which can be specified from observations, as described below. To solve for the steady-state velocity $u^\ast$
, we use mass conservation: $h^\ast u^\ast = h_0 u_0 - \dot {m} x$
. The length of the ice tongue cannot exceed
where the thickness vanishes due to basal melting alone. Here, we refer to this distance as the ‘mass-balance terminus’. For non-negative melt rates (zero or positive accumulation), the extent is unbounded (van der Veen, Reference van der Veen2013).
We now simplify the damage evolution equation and solve for an analytic expression for the steady-state damage in a freely floating ice tongue confined to flow in only one direction. Given the forces in the ice tongue in Eqn (14), we find that S 0 = 2, α = 0, $n^\ast = n$
and r N = ρi/2ρw ≈ 0.44. Equation (2) becomes
We can quantify the location along the thickness profile where the viscous-healing and melt-driven terms balance in Eqn (19) by setting dr/dt = 0. Using the analytic expression for ice thickness in Eqn (17) and conservation of mass, we can derive an analytic expression for this critical position x cr (see Appendix A):
Downstream of x cr, damage production is dominated by melt, whereas crevasses would be closed upstream of x cr due to strain-driven healing except that they are held open by the Nye zero-stress criterion in Eqn (5). For non-negative melt rates (zero or positive accumulation), the melt term has the same sign as the compressive stress and damage does not increase beyond the Nye zero-stress damage.
In an advective steady state (∂r/∂t = 0), beyond the critical distance at which damage production exceeds healing, the damage satisfies the following differential equation in x alone:
Using conservation of mass and bringing some expressions into the derivative gives:
which we may readily integrate to find
Rearranging we arrive at an analytic form for the damage along a 1-D ice tongue with a uniform, positive melt rate:
In the low melt rate limit, which smoothly approaches the zero melt profiles (see van der Veen, Reference van der Veen2013, Eqn5.71), gradients in the ice thickness and velocity are small away from the grounding line and thus the damage increases minimally from the Nye zero-stress damage. We use Eqn (24) to compare against the damage evolution results predicted from BISICLES.
We may also find the distance beyond the grounding line at which the ice tongue becomes fully damaged, with crevasses penetrating the entire depth. This ‘fully-damaged terminus’ is located at L r such that r(L r) = 1, or
which we compute numerically by root finding. Note that while Eqn (24) is valid for any constant lower bound of damage, Eqn (25) does not assure that L r < L max so that crevasses terminate the shelf before melt does. As a reminder, we have chosen the Nye zero-stress damage r N to be the lower bound of damage. We may also use Eqn (25) to extract a quantity that is useful for comparison to observations: the thickness of the ice at the fully-damaged terminus. By substitution and a bit of algebra (see Appendix A) we find that
which depends only on the uniform melt rate ($\dot {m}$
) and the temperature and density of the ice and water (within the expression for C, see Eqn (16)).
We determine the model parameters for the ice tongues (grounding line thickness and velocity, h 0 and u 0, uniform melt rate, $\dot {m}$
, and the viscosity rate-parameter, A) by geodesic-accelerated least-squares minimization (Transtrum and others, Reference Transtrum, Machta and Sethna2011) between the analytic solutions for thickness and velocity, and observations along central flowlines (Holdsworth, Reference Holdsworth1974; Wuite and others, Reference Wuite, Jezek, Wu, Farness and Carande2009; Blankenship and others, Reference Blankenship2012). This fitting method corrects for the tendency of traditional algorithms to be sensitive to the starting guess and to spuriously prefer infinitely rigid ice (A → 0). The best-fit parameters, as well as those within 1 SD, are listed in the first two columns of Table 1. The melt rate we infer for Erebus Glacier Tongue is ~2 m a−1, which is consistent with the basal melt rates calculated by Holdsworth (Reference Holdsworth1982) that vary spatially between −0.4 and 4.0 m a−1 along the length of the ice tongue. Drygalski Ice Tongue's inferred melt rate of ~3.1 m a−1 is also consistent with the melt rates reported by Wuite and others (Reference Wuite, Jezek, Wu, Farness and Carande2009) across the full ice tongue, which vary spatially between −1 m a−1 (accumulation) and 21 m a−1 along the final 65 km of unembayed ice.
Table 1. Simulation parameters: grounding line thickness h 0 and velocity u 0, embayment length ℓem and width w em (for the ice shelves), rate constant A (found in Eqn (16)) and uniform basal melt rate $\dot {m}$
The ice tongue values result from fitting centerline thickness and velocity observations (Holdsworth, Reference Holdsworth1974; Wuite and others, Reference Wuite, Jezek, Wu, Farness and Carande2009; Blankenship and others, Reference Blankenship2012). Our simulations take Drygalski Ice Tongue's input flux as the approximate value at a position 52 km downstream from the grounding line, at the edge of the embayment. For the ice shelves, the tabulated values are our best estimates of observations at the grounding lines (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011a, Reference Rignot, Mouginot and Scheuchl2011b; Fretwell and others, Reference Fretwell2013).
Rectangular ice shelves
To allow for damage evolution simulations that include the effects of ice-shelf buttressing, we construct two idealized ice-shelf geometries with aspect ratios and grounding line fluxes similar to those of the Amery and Ross ice shelves, also in Antarctica. We idealize each ice shelf as a 2-D rectangular box. In each of the two geometries, we apply lateral no-slip boundary conditions along a rectangular embayment prescribed to be approximately the length and width of the average observed embayment: 505 km long and 100 km wide for the Amery-like geometry; 650 km long and 950 km wide for the Ross-like geometry, as recorded in Table 1. We use 1 and 6 km as our coarsest horizontal grid resolution for the Amery- and Ross-like geometries, respectively. We define the constant input fluxes as corresponding to the grounding lines of each ice shelf, and prescribe them along one edge of the box, perpendicular to the embayment walls, using a uniform thickness and a smooth velocity profile that goes to zero at the walls (see Supplementary material Section 1). For the grounding line ice thickness and centerline velocity, we take representative values from Bedmap2 (Fretwell and others, Reference Fretwell2013) and the MEaSUREs InSAR-Based Antarctica Ice Velocity Map (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011a, Reference Rignot, Mouginot and Scheuchl2011b). We take representative values for average melt rate from Depoorter and others (Reference Depoorter2013), Schodlok and others (Reference Schodlok, Menemenlis and Rignot2016) and Rignot and others (Reference Rignot, Jacobs, Mouginot and Scheuchl2013) to be $\dot {m} = 0.8$
m a−1 for the Amery Ice Shelf and $\dot {m} = 0.25$
m a−1 for the Ross Ice Shelf, but further explore melt rate dependence in the ‘Discussion’ Section. For both idealized ice shelves, we prescribe Glen's flow rate parameter to be 1.29 × 10−17 Pa−3 a−1, which corresponds to a temperature of ~− 9.5°C, the average of a linear vertical temperature profile from an air-surface temperature of −19°C to a pressure-melt temperature of 0°C.
Results
Ice tongues
Figures 2 and 3 show the simulated steady-state damage and thickness profiles for the Erebus Glacier Tongue and the unembayed portion of the Drygalski Ice Tongue, downstream to the fully-damaged terminus L r (where damage r = 1), alongside their analytic 1-D solutions. Values of damage across both ice tongues are shown directly within each ice column as the depths to which basal crevasses penetrate. These simulated profiles are overlain by observed cross-sectional thicknesses from Holdsworth (Reference Holdsworth1974) for Erebus Glacier Tongue and from Operation IceBridge data for Drygalski Ice Tongue (Blankenship and others, Reference Blankenship2012).
Fig. 2. (a) Damage profile along the Erebus Glacier Tongue, computed by BISICLES and from the analytic solution in Eqn (24). The Nye zero-stress damage is shown for reference. (b) Thickness profile of the Erebus Glacier Tongue, comparing the results from our model (solid blue line) to data reported by Holdsworth (Reference Holdsworth1974) (dashed red line, solid circles). The ice is shaded purple up to the depth to which basal crevasses penetrate. Axis marks are provided for the analytically derived x cr (Eqn (20), where damage begins to increase), L r (Eqn (25), where damage equals 1) and L max (Eqn (18), the maximum mass-balance length of the ice sheet). Error bars are shown for uncertainties in the observed terminus position (dashed) and uncertainties in L max and L r propagated from the uncertainties in the parameter fit.
Fig. 3. Same as Figure 2 for the unembayed portion of the Drygalski Ice Tongue, with data from Blankenship and others (Reference Blankenship2012). Estimated uncertainties in L r and L max are 0.32 and 0.27 km, respectively, and too small to visualize on this scale.
For Erebus Glacier Tongue, our model predicts a steady-state fully-damaged terminus at 15.0 km beyond the grounding line. This is almost within the 11.6 ± 3 km observed range of decadal terminus variability reported by Holdsworth (Reference Holdsworth1974), although the recent variability has been slightly smaller (± 2 km, Stevens and others, Reference Stevens, Sirguey, Leonard and Haskell2013). We also see good agreement between the shapes of the simulated and reported thickness profiles, with a terminal thickness of 67 m. Uncertainties in the fitted glaciological parameters contribute to the disagreement between the simulation output and the data. Sampling from the best-fit distribution, we find 1 − σ uncertainties of ± 1.5 km in the damage-terminated length, shown in Fig. 2, and ± 5 m in the terminus thickness. (See Supplementary material Section 2 for more information on fitting and uncertainties.)
Damage reaches unity on the unembayed portion of the Drygalski Ice Tongue at a steady-state length of 60.5 ± 0.5 km beyond the embayment walls, ~6 km shorter than the 67 km observations used from the November 2011 Operation IceBridge flyover (Blankenship and others, Reference Blankenship2012) but still well within the 40 km range of observed decadal terminus variability (Frezzotti and Mabin, Reference Frezzotti and Mabin1994). The terminus thickness of 82 ± 1 m is also consistent with the observations.
These 1-D results verify the numerical implementation of damage evolution within BISICLES using the analytic result of Eqn (24) and validate the model as reproducing, within observational error, the length and terminal thickness of these ice tongues. We now turn to applying the model to the idealized 2-D geometries for which no closed form solution exists.
Idealized ice shelves
Figures 4 and 5 show results for geometries with the approximate width, incoming flux and melt rate of the Amery and Ross ice shelves, as recorded in Table 1. When we introduce no-slip boundary conditions to simulate an embayed ice shelf, like the Amery and Ross ice shelves, we also introduce spatial variations in stress transverse to the direction of flow that affect how damage evolves.
Fig. 4. (a) Simulation results overlain on the observed grounding line and ice front marked in white (solid and dashed, respectively) from Bedmap2 (Fretwell and others, Reference Fretwell2013), for the Amery-like ice-shelf geometry. The purple shades show the damage within the simulated ice extent, with the hatching indicating fully damaged ice (r = 1), and the green contours show the ice thickness. The grounding line thickness, the centerline velocity at the grounding line, the temperature and the melt rate are as prescribed in Table 1. Analytic predictions from the flux-equivalent 1-D ice tongue (see Supplementary material Section 1) for x cr, L r and L max (Eqns (20), (25) and (18)) are shown along the top axis. (b) Evolved damage and the Nye damage from BISICLES along the center line, compared to the prediction from the 1-D flux-equivalent model. (c) Size and location of the idealized domain for Bedmap2 contours in panel (a). (d) Thickness profile along the center line, as in Fig. 2b, with thickness data (red, dotted) from Bedmap2.
Fig. 5. Same as Fig. 4, for the Ross-like ice shelf with the inset map in panel (c) rotated counterclockwise so the flux in the domain flows left-to-right.
Buttressing stresses from the walls extend the region where the stresses are compressive, reducing damage and allowing the distances to the fully-damaged terminus (identified as the start of hatching in Figs 4a and 5a) in these embayed ice shelves to extend further than that of their unembayed, or free-slip counterparts (shown as vertical lines labeled L r). This reduction in damage is more pronounced in the Amery-like ice shelf with its width of 100 km (Fig. 4a) than the Ross-like ice shelf (Fig. 5a) with the almost an order of magnitude larger width of its embayment (950 km).
In both simulations, the observed length between the grounding line and terminus position (solid and dashed white lines, respectively) is approximately reproduced by the centerline distance to the fully-damaged terminus, although this is not a significant improvement over the predicted length from the equivalent floating ice-tongue (L max). The thickness at the fully-damaged terminus is 74 and 52 m along the center line for the Amery- and Ross-like geometries, respectively. The fully-damaged terminus thickens near the embayment margins, where increased shear due to friction from the walls causes damage to reach the critical value of r = 1 in thicker ice.
Figures 4b, d and 5b, d display the damage (b) and thickness profiles (d) of the ice shelf along the center line in Figs 4a and 5a, respectively. We see that the damage follows the Nye zero-stress minimum damage (Eqn (5)) through most of the length of the ice shelf, where the flow is largely compressive. Furthermore, due to buttressing from the walls resisting the driving stress, the Nye zero-stress crevasse penetration depth initially decreases from the free-slip ice tongue value of r N = 0.44. Similarly, the buttressed shelf is thicker than the unbuttressed, free-slip tongue with the equivalent incoming flux, but is not as thick as the observed thicknesses.
The effect of buttressing on damage is more pronounced for the narrower confinement, with crevasses only penetrating 11% through the thickness at its minimum for the Amery-like ice shelf (compared to 30% for the Ross-like geometry and the constant 44% for the 1-D ice tongue). Damage that exceeds the Nye zero-stress damage is concentrated in a narrow band where melt rate and spreading dominate the damage evolution, and rapidly increases until the threshold is reached at 304 km for the Amery-like domain. This region points to the importance of the source term in the damage evolution (first term in Eqn (2)) compared to the Nye zero-stress specification of Sun and others (Reference Sun, Cornford, Moore, Gladstone and Zhao2017). In the Ross-like domain, due to its comparatively low melt rate (0.25 m a−1 compared to 0.8 m a−1), the departure from the Nye zero-stress damage along the center line is smaller, and the damage threshold is reached at 543 km along the center line – 12 km closer to the grounding line than the Nye zero-stress damage predicts.
Discussion
The fully damaged terminus reproduces the terminus positions of the Erebus and Drygalski ice tongues, within the range of observations, despite the two systems having lengths that differ by a factor of five. Figure 6 is a summary of the ice thickness at the fully damaged terminus from all the simulations described above, as compared to the analytic prediction for a confined ice shelf with free-slip walls (Eqn (26)). Our geometries with no-slip walls showed ice thicknesses at the fully damaged termini of 74 m along the center line for the Amery-like parameters and 52 m for the Ross-like parameters (solid points in Fig. 6). Although much smaller than the observed thicknesses of ~250 m for both ice shelves (Fretwell and others, Reference Fretwell2013), the fully damaged terminus is a significant improvement for simultaneous predictions of the location and thickness of the terminus relative to the mass-balance terminus – which predicts zero thickness at the terminus by definition – without heuristic specifications of how thick the ice should be at the terminus. Still, important discrepancies exist between the observations and our predictions that point to factors missing from our simulations that contribute to the evolution of crevasses and calving fronts.
Fig. 6. (a) Ice thickness at h(L r) and (b) distance from the grounding line to the fully damaged terminus as it depends on melt rate and average ice temperature. Lines in (a) show the analytic prediction for a 1-D ice tongue with temperatures corresponding to the best fits for Erebus Glacier Tongue (−5.9°C) and Drygalski Ice Tongue (−8.8°C), as well as for the prescribed ice-shelf temperature used in the Amery-like and Ross-like simulations (−9.4°C). Lines in (b) show the analytic prediction for each flux-equivalent tongue for each geometry. For the idealized Amery and Ross Ice Shelf geometries, which are 2-D by nature of the no-slip embayment walls, the thickness shown here is along the center line. Solid dots represent the experiments from Table 1, whereas open dots represent additional experiments discussed in text.
We have made several simplifying assumptions about the (1) mass-balance forcing, (2) boundary conditions, and (3) mechanical properties of damaged ice. In our simulations, we assumed uniform melt rates across the ice shelves to demonstrate the physics of the damage evolution equation (Eqn (2)) and facilitate the comparison to an analytic solution. We showed that for the ice tongue case, the thickness of the fully damaged terminus was entirely controlled by the uniform melt rate and the temperature of the ice (Eqn (26)), with no dependence on the flux across the grounding line, unlike the distance to the fully damaged terminus (Eqn (25)). The three curves in Fig. 6a show the thickness dependence on melt rate for the three temperatures used in the study and the four curves in Fig. 6b show the length dependence for all four geometries. Physically, higher melt rates increase damage production, pushing the location of the fully damaged terminus closer to the grounding line and into thicker ice. On the other hand, warmer ice deforms and thins more easily, leading to a smaller thickness, but a less significant effect on the length. The thicknesses we measure from the BISICLES simulations of Erebus and Drygalski ice tongues fall right on these curves, as expected from the agreement between the simulated and analytic profiles.
The confined ice shelves with no-slip walls show a generally similar melt rate dependency (blue, dashed lines in Fig. 6a). The Amery-like ice shelf, with its narrow confinements, produces a thicker terminus for smaller melt rates than predicted by the 1-D theory. The buttressing from the walls slows the ice, allowing it to thicken relative to an equivalent free-slip ice tongue. However, with a higher melt rate, the smaller ice thickness results in less contact with the wall, less buttressing, and thus the shelf thickness looks more like that of the analytic tongue.
The large width of the Ross-like shelf results in a more pronounced melt rate dependence along the center line, owing to a balance between the buttressing that causes both the damage to heal and the ice to thicken. At low values of the melt rate, the ice shelf is connected to the wall for longer, allowing the crevasses to heal more (the Nye zero-stress damage decreases with increased buttressing) so that the fully damaged terminus location is pushed farther from the grounding line (cf. the Ross points in Fig. 6b diverge from the theoretical 1-D prediction), into thinner ice, while the thickening effect is not as pronounced. For high melt rates, the ice shelf decouples from the walls almost immediately, and so the ice spreads out rapidly. Between these extremes, the healing is not as significant but the buttressing allows the ice to thicken, resulting in the non-monotonicity in Fig. 6a. These melt rate experiments suggest that a thicker fully damaged terminus could result from a more representative melt-rate profile (e.g. decreasing from the grounding line to the terminus at Amery as in Wen and others (Reference Wen2010) or increasing for Ross as in Stewart and others (Reference Stewart, Christoffersen, Nicholls, Williams and Dowdeswell2019)), but only in combination with a more complete representation of the boundaries of the ice shelves, as we describe next.
Our approximation of the domain as a rectangular embayment with no-flux, no-slip walls gives a general sense of how friction and lateral shearing affect the evolution of damage. In regions near the side walls, where strain rates are large, damage evolves more quickly, such that the length from grounding line to fully damaged terminus decreases (to ~160 km in the Amery-like geometry and 134 km in the Ross-like geometry), causing the ice shelf to terminate in thicker ice (221 and 91 m, respectively). However, comparing our boundaries to the observed grounding lines (white contours in Figs 4 and 5) highlights topographical features in the actual ice shelves that are absent in our simulations. This includes pinning points along the walls of the Amery Ice Shelf and islands in the Ross Sea, which would provide buttressing and allow the modeled ice shelves to thicken and extend farther beyond their grounding lines. These pinning points, and a more complex grounding-line flux, would also bring the modeled thicknesses closer to the observed thicknesses (compare solid-blue and dotted-red lines in Figs 4d and 5d).
We expect damage to be affected by the inclusion of these topographical features as well. As with the walls in our domains, ice flow around these features would experience increased shear rates, which would lead to faster accumulation of damage and a thicker fully damaged terminus. This is consistent with an empirical damage estimate along the Amery and Ross ice shelves by Bassis and Ma (Reference Bassis and Ma2015, their Figs 5a and c), who integrated the damage evolution Eqn (2) using satellite-derived thicknesses and velocity profiles. They found, as here, that damage is low throughout the majority of the floating ice shelf, with a rapid increase occurring near the terminus. What their estimates show that ours do not, however, are flowlines of high damage emanating from topographical features along the walls and from features along the grounding lines. Similarly, Indrigo and others (Reference Indrigo, Dow, Greenbaum and Morlighem2021) find that fracture propagation at Drygalski Ice Tongue is controlled by the thickness variations in the form of across- and along-flow basal channels that evolve from grounding line to calving front, suggesting an important interplay between grounding line morphology, melt and damage production. For the Drygalski Ice Tongue, although thicker ice between along-flow channels impedes fracture propagation, the large-scale calving behavior is determined by the thinning of the across-flow channels (Indrigo and others, Reference Indrigo, Dow, Greenbaum and Morlighem2021), in line with the necking instability we describe here, and captured by our 1-D model. For the Amery Ice Shelf, however, with buttressing forces and the accretion of marine ice (Fricker and others, Reference Fricker, Popov, Allison and Young2001), this interaction between grounding line and basal features will have a more pronounced effect on damage. We therefore believe that a more representative grounding line topography and the presence of pinning points would contribute to a more varied fully damaged terminus position: further from the grounding line in areas where buttressing lessens the shear stress and closer to the grounding line downstream of high shear zones.
Finally, we have assumed that damage has no effect on the dynamics of the ice and focused solely on the locations where damage reaches its maximum value r = 1. When a crevasse penetrates through the entire ice thickness, ice should detach, calve away and change the stress balance on the front. Calved ice does not evacuate to the open ocean immediately, however, but can persist in front of the ice as a mélange, pushing back as a buttressing force on the intact ice (e.g. Burton and others, Reference Burton, Amundson, Cassotto, Kuo and Dennin2018). Our results can be understood as one end-member on the spectrum between the presence of a mechanically strong mélange and ice that is immediately evacuated away from the terminus. Any decrease in the viscosity of the damaged ice as a result of coupling to the ice rheology will decrease the back-stress exerted on the still-intact ice and push the fully damaged terminus upstream, into thicker ice. Similarly, if calving occurs at a threshold lower than r = 1 (see, e.g. Pralong and Funk, Reference Pralong and Funk2005), the terminus would also fall closer to the grounding line, in thicker ice.
Challenges arise in implementing and interpreting simulations in which we remove damaged ice from the idealized ice-shelf geometries. Damage at the point where the moving ice, static-wall and ocean meet always grows due to a singularity in the stress, causing the ice shelf to decouple from the wall and thin considerably, similar to the rapid retreat observed by Sun and others (Reference Sun, Cornford, Moore, Gladstone and Zhao2017) in the MISMIP+ experiments with damage calving. Possible resolutions to this challenge will likely involve relaxing the no-slip boundary condition, coupling to a mélange model (e.g. Amundson and Burton, Reference Amundson and Burton2018; Schlemm and Levermann, Reference Schlemm and Levermann2021), or introducing a lubrication layer, such as a shear margin (e.g. Lhermitte and others, Reference Lhermitte2020), as we would get with geometries that depict more of the details of observed glacial boundaries.
Figures 4 and 5 show that that our model does not predict localized rifting or shear margins in these idealized geometries. This is particularly apparent in our idealized Amery-like simulations, as we do not resolve the tens of kilometers long rift system observed at the terminus of the real ice shelf (Bassis and others, Reference Bassis, Coleman, Fricker and Minster2005). Measurements made by Bassis and others (Reference Bassis, Fricker, Coleman and Minster2008) suggest that rift propagation depends primarily on the stress within the ice and not on short-term climate forcings or ice–ocean interactions. Our model (Eqn (2)) includes a dependence on the stress through the maximum principal strain rate and Glen's flow law, but this relationship also depends on the temperature through the effective ice viscosity. Coupling damage to the stress field via rheological feedbacks or incorporating more representative ice temperature evolution could therefore lead to the localization of damage in our model, similar to what was seen in Lhermitte and others (Reference Lhermitte2020). We have also omitted the process of brittle fracture that would directly cause such localized rifting on the scale of meters to tens of meters (Pralong and Funk, Reference Pralong and Funk2005; Aström and others, Reference Aström2013), which could possibly be resolved within BISICLES using its adaptive mesh refinement capability.
Conclusions
Our continuum damage mechanics model, which simulates the evolution of the ratio of crevasse depth to ice thickness according to a pseudo-plastic necking instability, provides a useful framework for modeling damage evolution and terminus characteristics without introducing additional parameters. When incorporated into an ice-sheet model, our fully damaged terminus model predicts broadly accurate steady-state extents for a suite of idealized, isothermal ice tongues and ice shelves forced by spatially uniform basal melt rates. We have provided an analytical expression for the damage in an ice tongue and predicted a fully damaged terminus thickness that increases with melt rate, decreases with ice temperature, and has no dependence on flux or system size. Whereas we are able to reasonably model observed lengths and terminus thicknesses for ice tongues, we identify more complex behavior in ice shelves with no-slip walls. Specifically, the fully damaged termini modeled in our idealized, 2-D ice-shelf geometries are thinner than those observed. The likely causes of these differences arise from certain idealizations in the present study (simplified grounding line, spatially uniform ocean forcing, flat, no-slip side walls and the omission of other pinning points) and to interactions between damage and the rheological properties of ice. Integration of these aspects of observed ice shelves with the necking instability we describe will be an effective way to model the future evolution of calving fronts in large-scale ice-sheet models.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2022.12.
Code and data availability
We used the damagedBISICLES branch of the publicly available version of the BISICLES ice-sheet model code, release version 1.0. Instructions for downloading and installing BISICLES after free registration with ANAG may be found in the ‘getting started’ section at http://bisicles.lbl.gov. The specific svn command for obtaining the relevant branch is: svn co https://anag-repo.lbl.gov/svn/BISICLES/public/branches/damagedBISICLES BISICLES.
BISICLES is written in a combination of C++ and FORTRAN and is built upon the Chombo AMR software framework. More information about Chombo may be found at http://Chombo.lbl.gov.
Static code, data, input, configuration files for the runs in this work are available at https://portal.nersc.gov/cfs/iceocean/iceshelfdamage or https://doi.org/10.5281/zenodo.5850262.
Acknowledgements
Support for this work was provided through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy (DOE), Office of Science, Biological and Environmental Research and Advanced Scientific Computing Research programs, as a part of the ProSPect SciDAC Partnership. This work is from the DOMINOS project, a component of the International Thwaites Glacier Collaboration (ITGC). Support was provided by National Science Foundation (NSF: Grant 1738896) and Natural Environment Research Council (NERC: Grant NE/S006605/1). Logistics provided by NSF-U.S. Antarctic Program and NERC-British Antarctic Survey. ITGC Contribution No. ITGC-063. Work at Berkeley Lab was supported by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center supported by the Office of Science of the U.S. DOE under contract no. DE-AC02-05CH11231. The authors also thank two anonymous reviewers for their thoughtful comments and Editor Matthew Siegfried for his guidance.
Appendix A
Detailed derivation of thickness at the fully damaged terminus of an ice tongue
In this appendix, we derive the thickness of the fully damaged terminus, by way of deriving x cr. We start by introducing a non-dimensionalization which is particularly useful for 1-D tongues with positive, uniform melt rates. As these have a finite extent where the total melt rate equals the flux in $L_{\max } \dot {m} = h_0 u_0$
, we scale distances with this length. Inspection of the analytic thickness suggests that the thickness should scale like $( \dot {m}/C) ^{1/( n + 1) }$
. We thus introduce non-dimensional variables $\tilde {x} = x/L_{\max }$
and $\tilde {h} = h/( \dot {m}/C) ^{1/( n + 1) }$
. The velocity scale turns out not to matter, so we leave it dimensional for the time being. In non-dimensional form, the thickness from grounding line ($\tilde {x} = 0$
) to mass-balance terminus ($\tilde {x} = 1)$
is
It is useful to note from this that
and
Remembering that the forces in the ice tongue give S 0 = 2, α = 0, $n^\ast = n$
and r N = ρ i/2ρ w ≈ 0.44, the location at which damage increases, which is also the location where the principal stress regime changes from compressive to tensile, can be found with
Simplifying Eqn (A4) for the thickness at this critical distance gives
We are now ready to simplify the damage evolution equation and solve for an analytic expression for the critical distance where damage increases in a freely floating ice tongue confined to flow in only one direction. Using Eqn (A1), we find that
which, after re-dimensionalizing, gives us Eqn (20).
We may now find the thickness of the fully damaged terminus at advective steady state, $h^\ast ( L_r)$
. We know from Eqn (25) that L r is defined such that
Upon substituting in Eqn (A2), non-dimensionalizing, and rearranging terms, we obtain
Using Eqns (A3) and (A5) and canceling some terms:
Thus, the non-dimensional thickness at the fully damaged terminus, $\tilde {h}^{\ast }( \tilde {L}_r) = \alpha$
, is the root to the polynomial
which has an analytic solution where n is an integer of 3 or less, but we find it numerically by root-finding. For the material parameters ρ i = 910 kg m−3, ρ w = 1028 kg m−3 and n = 3, we find that α = 0.25. Back in dimensional form, the thickness of the fully damaged terminus is
which is Eqn (26).
Open access
