A stress-based poro-damage phase field model for hydrofracturing of creeping glaciers and ice shelves

There is a need for computational models capable of predicting meltwater-assisted crevasse growth in glacial ice. Mass loss from glaciers and ice sheets is the largest contributor to sea-level rise and iceberg calving due to hydrofracture is one of the most prominent yet less understood glacial mass loss processes. To overcome the limitations of empirical and analytical approaches, we here propose a new phase field-based computational framework to simulate crevasse growth in both grounded ice sheets and floating ice shelves. The model incorporates the three elements needed to mechanistically simulate hydrofracture of surface and basal crevasses: (i) a constitutive description incorporating the non-linear viscous rheology of ice, (ii) a phase field formulation capable of capturing cracking phenomena of arbitrary complexity, such as 3D crevasse interaction, and (iii) a poro-damage representation to account for the role of meltwater pressure on crevasse growth. A stress-based phase field model is adopted to reduce the length-scale sensitivity, as needed to tackle the large scales of iceberg calving, and to adequately predict crevasse growth in tensile stress regions of incompressible solids. The potential of the computational framework presented is demonstrated by addressing a number of 2D and 3D case studies, involving single and multiple crevasses, and considering both grounded and floating conditions. The computational framework presented opens new horizons in the modelling of iceberg calving and, due to its ability to incorporate incompressible behaviour, can be readily incorporated into numerical ice sheet models for projecting sea-level rise.


Introduction
Ice sheets are large masses of glacial ice that inundate the surrounding landscape in Greenland and Antarctica today, and many other regions during ice ages [1]. These act as enormous stores of freshwater -containing approximately 70% of the planet's supply [2] -that assist in regulating a stable global climate, through maintaining global ocean-water levels and controlling surface temperatures by reflecting solar radiation due to its high albedo properties [3]. Ice sheets thin toward their margins, and if these are located in marine settings, they will form floating extensions known as ice shelves, which act to provide resistive buttressing to downslope flow and reduce the flux of grounded ice to the ocean. However, increasing global temperatures as a result of carbon emissions has lead to higher rates of ablation than accumulation, resulting in ice shelf and ice sheet thinning in some key areas where ice-sheet instability may follow [4]. Surface and basal crevasses can form within ice sheets as a consequence of ongoing deformations within the ice. These are deep crack-like defects that can propagate in an unstable manner and lead to large-scale iceberg calving events, and in extreme cases the catastrophic break up of ice shelves. The frequency of these events has grown in recent decades, beginning with the disintegration of Larsen A (1995) [5] and Larsen B (2002) [6] ice shelves, and more recently significant surface melting and iceberg calving on Larsen C (2017) [7], Pine Island and Thwaites (2018-2020) [8], and Conger (2022) ice shelves. Fracture within ice shelves can result in a loss of resistance to down slope glacial flow, leading to ice-sheet thinning, additional flotation of grounded ice and, thus, potentially irreversible grounding line retreats [9].
Deposition of grounded glacial ice into the ocean is one of the leading contributors to sea level rise [10], having direct implications within this Century on low-lying coastal regions through flooding, increased extreme environmental events, degradation of farmland and loss of habitat, among others. A key driving factor for their stability is the production of surface meltwater as a result of elevated surface temperatures [11]. When ice shelves and glaciers melt, meltwater flows down-slope into surface crevasses, causing additional tensile stresses to form within the crevasse.
This leads to crevasse instability, and with sufficient meltwater, the crevasse can propagate through the full thickness of the ice column. This process is generally referred to in the glaciological literature as hydrofracture [12]. A recent study by Lai et al. [13] found that approximately 60 ± 10 % of Antarctic ice shelves provide significant buttressing to downslope flow and are vulnerable to meltwater driven hydrofracture, highlighting the significance of studying the formation and propagation of crevasses in glaciers. An illustration of a grounded ice sheet, transitioning to a crevassed floating ice shelf is shown in Fig. 1. Ice sheet fracture and crevasse propagation have been mainly modelled previously using analytical methods. The estimation of crevasse penetration depths in an idealised glacier was first described by Nye in 1955 [14] based on the so-called 'zero stress' model. Nye assumed that ice has no tensile resistance to fracture and that a crevasse will stabilise at a depth where the longitudinal tensile stress is balanced by the lithostatic compressive stress [15]. This was later extended by Benn et al. to include the presence of meltwater within crevasses [16]. Linear elastic fracture mechanics models were introduced to provide a more accurate prediction of the depth of an isolated crevasse [17,18]. By exploiting the principle of superposition, stress intensity factors can be calculated by integrating over the crevasse depth for the normal tensile stress, the lithostatic compressive stress, and the meltwater pressure. In order for crevasses to stabilise, the net stress intensity factor K net must be equal to the material's fracture toughness K c . However, these analytical approaches have well-known limitations, such as (i) idealised scenarios and boundary conditions are assumed; (ii) creep effects, resulting from the continual movement of glaciers under their own weight, are neglected; and (iii) crevasse interaction cannot be captured.
Recently, computational methods have been used to predict crevasse growth and iceberg calving events. Local and non-local continuum damage mechanics formulations have been presented to predict ice sheet fracture [19][20][21][22]. These works have overcome some of the limitations intrinsic to analytical approaches, but often at the cost of using empirical parameters. Variational phase field fracture models offer an alternative approach, enabling the simulation of realistic conditions (3D geometries, multiple interacting crevasses, etc.) and providing a connection to fracture mechanics theory. Phase field fracture models have gained remarkable popularity in recent years due to their ability to predict complex cracking phenomena including crack bifurcation, coalescence and nucleation from arbitrary sites [23][24][25]. New phase field-based formulations have been presented for dynamic fracture [26,27], ductile damage [28,29], environmentally assisted cracking [30,31], fatigue crack growth [32,33], hydraulic fracture [34,35], and battery degradation [36,37]; among other (see Refs. [38,39] for an overview). In this work, we aim at extending the success of phase field fracture models to the area of glacier crevassing and iceberg calving. To this end, a new phase field formulation is presented capable of capturing the creep behaviour of glacial ice and the role of fluid pressure in driving crevasse growth. Also, for the first time, crevasse interaction is predicted in both 2D and 3D. Very recently, Sun et al. [40] used a phase field approach to predict hydrofracture in 2D linear elastic glaciers, assuming compressible behaviour and disregarding creep effects. Unlike them, we base our framework on a stress-based phase field fracture formulation, which offers several advantages in the context of hydrofracturing of glacier crevasses. First, strain energy-based approaches are unsuited for incompressible rheologies. This is not only important due to the incompressible nature of glacial ice, but also because it hinders its integration into largescale computational models for ice sheet evolution and sea level rise, which assume incompressible flow (see, e.g., the Community Ice Sheet Model (CISM) [41]). Second, ice-sheet fracture is driven by tensile stresses and not strains, with crevasses propagating solely in regions where the net longitudinal stress is positive [42]. This is naturally accounted for in a stress-based phase field model, while requiring a particular ad hoc split in strain energy-based formulations [40,43]. Third, a phase field length-scale insensitive driving force can be defined, enabling the use of coarser meshes, a key enabler given the large scales involved. These advantages provide further motivation for this work, presenting the first stress-based phase field computational framework for hydrofracturing of creeping glaciers and ice shelves.
The rest of the paper is outlined as follows. The theoretical and computational framework presented is described in Section 2. The model is then used in Section 3 to predict hydrofracturing in case studies of particular interest. First, the propagation of single crevasses in grounded ice considering both linear and non-linear rheologies is investigated. A parametric study is conducted to assess the role of relevant material parameters, seawater level and meltwater depth. Second, we simulate the growth of a field of densely spaced crevasses in a grounded glacier, comparing against the predictions of Nye's zero stress model. Third, the growth of basal and surface crevasses (and their interaction) is for the first time simulated for a floating ice shelf, using appropriate Robin boundary conditions. Fourth, the combined creep-phase field fracture model is used to predict the nucleation and growth of crevasses in a realistic geometry, corresponding to the Helheim glacier.
Finally, we provide the first 3D analysis of crevasse propagation in ice sheets. Concluding remarks end the manuscript in Section 4.

Numerical framework
In this section, we present our computational framework, which encompasses the three elements that are needed to resolve the hydrofracture process taking place in ice sheets; namely, the viscoplastic behaviour of ice, the propagation of meltwater-filled crevasses, and the role of meltwater pressure on crevasse propagation. These are modelled by means of Glen's flow law [44], a stress-based phase field description of fracture [45], and a meltwater-ice poro-damage model [21], 0 1 Fig. 2: Schematic diagram of a meltwater filled crevasse in glacial ice, illustrating the intact phase (φ = 0), fully cracked phase (φ = 1) and transition phase (0 < φ < 1). In the damaged and transition phases, there is a hydrostatic pressure contribution to damage arising from the meltwater. Relevant to the poro-damage part of the model, h s denotes the meltwater depth, and z s is the distance between the glacier base and the bottom of the crevasse, with z being the vertical height.
respectively. Fig. 2 illustrates upon a single crevasse the mechanistic and modelling assumptions of our framework. In the following, we present the kinematics of the problem (Section 2.1), formulate the energy functionals (Section 2.2), particularise the model upon suitable constitutive choices (Section 2.3), and briefly describe the finite element implementation (Section 2.4). Throughout, the formulation refers to a body occupying an arbitrary domain Ω ⊂ IR n (n ∈ [1, 2, 3]), with an external boundary ∂Ω ⊂ IR n−1 with outwards unit normal n.

Kinematics and general considerations
The primary variables are the displacement field vector u and the damage phase field φ. Restricting our attention to small strains and isothermal conditions, the strain tensor ε reads with the strain field itself being additively decomposed into its elastic and viscous parts, such that The growth of meltwater-filled crevasses is described by means of a smooth continuous scalar phase field, which takes a value of φ = 0 in intact ice and of φ = 1 in fully damaged regions (see Fig. 2). The aim is to overcome the need to track discrete crack surfaces, which is a cumbersome task. Thus, the intact ice-crack interface is not explicitly modelled but instead smeared over a finite domain, replacing interfacial boundary conditions by a differential equation that describes the evolution of the phase field φ. The smearing of the interface is controlled by a phase field length scale . Accordingly, a discontinuous surface Γ is regularised through the following crack surface functional [46]: where γ is the so-called crack surface density functional.

Energy functionals
A total potential energy can be defined by incorporating the contributions from the bulk strain energy density ψ s , which itself accounts for both viscous (ψ v s ) and elastic (ψ e s ) contributions, and the regularised fracture energy ψ f . Thus, considering the work done by external tractions T and body forces b, the total potential energy of the solid can be expressed as, As discussed below, it is important to consider as well the kinetic energy of the body, which is given by where ρ is the mass density of the material andu = ∂u/∂t. The Lagrangian for the coupled deformation-fracture problem can then be formulated by combining the kinetic and total potential energies, such that We shall now make constitutive assumptions and, building upon these, proceed to formulate the local force balances.

Constitutive theory
We proceed to particularise our choices with the aim of providing a suitable framework for predicting ice-sheet hydrofracture. To this end, the bulk strain energy density of the solid is given in terms of its elastic and viscous counterparts as, where σ 0 is the undamaged Cauchy stress tensor, λ and µ are the Lamé parameters, and g (φ) is a phase field degradation function, to be defined. Then, the homogenised (damaged) stress tensor can be estimated as σ = ∂ ε e ψ s . As described below, the viscous behaviour of the solid is described by Glen's flow law [44], a commonly used choice for glacial ice.
where A is the creep coefficient, σ 0 = σ 0 − tr(σ 0 )I/3 is the undamaged deviatoric stress tensor, n is the creep exponent, and σ e is an equivalent stress measure defined as σ e = 1 2 σ 0 : σ 0 . The creep coefficient A and the creep exponent n are typically calibrated with experiments, with the former exhibiting the following Arrhenius dependency with temperature, where T is the absolute temperature, Q is the activation energy, R is the universal gas constant, and A 0 is the creep coefficient at a reference temperature T 0 .

A stress-based phase field fracture model
The evolution of damage is driven by the phase field variable φ. A length-scale insensitive, stress-based approach is adopted, inspired by the work by Miehe et al. [45]. This choice enables us to capture purely stress-driven fractures in incompressible solids using relatively coarse meshes; as required to model hydrofractures in creeping glaciers. Accordingly, the fracture energy density is formulated as, Unlike conventional phase field fracture models, Eq. (10) shows that the present formulation introduces the phase field through a linear term. This naturally results in a damage threshold, below which φ = 0, preserving the elastic properties of uncracked regions. In (10), ψ c is a fracture energy density, which in a stress-based approach is defined as a function of a critical fracture stress or material strength σ c , such that [45]: Here, E denotes the material's Young's modulus. It remains to define the degradation function g(φ), which reduces the elastic stiffness of the solid -see Eq. (7). The choice of g(φ) must fulfill the following conditions, Here, we choose to adopt the following quadratic degradation function Finally, the phase field evolution law is given by [45], Where the left hand side is the geometric resistance and the right hand side corresponds to the driving force. Here, D d is the crack driving force state function, which is here defined based on the principal tensile stress criterion, such that Such a crack driving force state function is adequate for fractures resulting from the decohesion of surfaces perpendicular to the maximum principal stress and provides a quadratically increasing stress threshold for stress levels above a failure surface in the principal stress space, as determined by the material strength σ c . Also, Eq. (15) provides a criterion independent of the phase field length scale , which minimises the sensitivity of the results to this parameter. Given that the finite element mesh has to be sufficiently fine to resolve , typically requiring an element size seven times smaller [47], this facilitates tackling the large scales inherent to iceberg calving. For completeness, a non-dimensional parameter ζ has been introduced that, for ζ = 1 values, influences the slope of the stress-strain curve in the post-critical range. This is shown below by exploring the one-dimensional predictions of (14) and (15). Hence, the evolution of the phase field in a one-dimensional setting (∇φ = 0) is given by, and accordingly the damaged (homogenised) uniaxial stress is found by making use of the following relationship, where ε is the uniaxial strain. The responses obtained are shown in Fig. 3, for selected choices of the parameter ζ. A linear response is predicted until the critical fracture stress is reached, with the post-critical regime being sensitive to the value of ζ; higher values translate into a less dissipative damage process, with the response appearing to converge for ζ > 5. For simplicity, we will assume ζ = 1 but will also consider its influence in a parametric study.

A porodamage description of meltwater-driven crevasse growth
Meltwater plays a key role in crevasse propagation, introducing local tensile stresses that can become equal or larger than the lithostatic compressive stress. It is thus pivotal to incorporate the role of the water pressure p w in the damaged (φ = 1) and transition (0 < φ < 1) regions, as meltwater can accumulate in damaged zones and in the localised pore structure that arises in the transition region due to the nucleation, growth and coalescence of microvoids and microcracks. To this end, we follow Terzaghi's concept of an effective stress [48] and Biot's theory of poroelasticity [49]. Hence, the resulting stress tensor is defined as, where α is Biot's coefficient. In this work, α = 1. The use of degradation functions in Eq. (18) constrains the water pressure to damaged regions and removes the load carrying capacity of ice in fractured domains. Here, the water pressure p w is a hydrostatic term that is depth dependent.
For surface crevasses it is defined as, where ρ w is the density of freshwater, h s is the meltwater depth, z is the vertical height and z s is the distance between the glacier base and the bottom of the crevasse (see Fig. 2). The presence of the Macaulay brackets in Eq. (19) implies that the pressure is zero above the water surface. Also, it is important to note that z s is updated for every time increment, as defined by the minimum depth at which φ = 1. Consequently, the role of meltwater pressure extends beyond the initial damage zone and appropriately evolves with the propagating crevasse. On the other hand, for basal crevasses it is assumed that the crevasse is fully saturated with ocean-water at depths below the ocean-water level h w . The water pressure within basal crevasses is then given by In this context, the material density is interpolated as a function of the damage state, and the freshwater (ρ w ) and glacial ice (ρ i ) densities, reading

Finite Element implementation
Finally, we proceed to formulate the particularised coupled balance equations and briefly describe the finite element implementation. Considering the constitutive choices described in Section 2.3, the local force balances are given by, with the natural boundary conditionsσ Here, C 0 is the elastic stiffness tensor and the ansatz in the right hand side of Eq. (23) is introduced to ensure damage irreversibility. The discretised system resulting from the weak form of (22)- (23) is solved using a so-called multi-pass (alternate minimization) staggered scheme [50]. An implicit BDF time-stepping scheme is employed to solve, in a Backward Euler fashion, each set of equations.
The commercial finite element package COMSOL is used.

Results
In this section, we present a series of 2D and vertical coordinates. Gravitational load due to self-weight is applied as a uniform body force in the z-direction with a magnitude of −ρ i g. We also consider the surface meltwater pressure p w within a crevasse using the poro-damage approach presented in Eq. (19). A Neumann-type traction is applied normal to the ice-ocean interface at the terminus, with the hydrostatic ocean-water pressure varying linearly with depth and a magnitude of −ρ s g h w − z . Boundary conditions that are specific to the grounded glacier and floating ice shelf cases are discussed in sections 3.1 and 3.3, respectively. Our simulations deal with glacial ice, whose material properties are given in Table 1, along with the densities of seawater and meltwater. Poisson's ratio, ν [-] 0.35 [51] Density of glacial ice, ρ i [kg/m 3 ] 917 [52] Density of meltwater, ρ w [kg/m 3 ] 1000 [52] Density of seawater, ρ s [kg/m 3 ] 1020 [52] Fracture toughness, K c [MPa √ m] 0.10 [53] Critical fracture stress, σ c [MPa] 0.1185 [54] Creep exponent, n [-] 3 [22] Creep coefficient Tab. 1: Material properties assumed in this work (unless otherwise stated). The values are chosen to characterise the behaviour of glacial ice, with the subscript number denoting the relevant reference.
The strength σ c magnitude is chosen to be an intermediate magnitude within the experimentally reported values of the critical fracture stress in glacial ice, which are in the range 0.08-0.14 MPa [53,56,57]. An estimate of the phase field length scale, which plays a negligible role in this model, can be obtained through the Hillerborg et al. [58] relation, which for plane strain reads: Considering the toughness of glacial ice (K c = 0.1 MPa √ m), this gives a magnitude of = 0.625 m, which is the value adopted here (unless otherwise stated). To attain mesh-independent results, the characteristic element size along the crevasse propagation region is always chosen to be at least 5 times smaller than the phase field length scale .

Propagation of a single crevasse on a grounded glacier
We begin our numerical experiments by gaining insight into the behaviour of crevasses in   [52]. The effect of including the ocean-water pressure at the glacier terminus on the far field longitudinal stress can be observed by comparing Figs. 5a and 5b. Here, the ocean-water pressure provides a compressive stress that is constant with depth (in the far field region) and that decreases the extent of the tensile stress region near the top surface. If the ocean-water height is sufficiently large (≈ 90% of ice thickness), this can cause the glacier to become buoyant and form a floating ice shelf/tongue, resulting in an increased compressive stress regime. Vertical stress predictions (not shown) exhibit a behaviour that is also invariant with x-coordinate and that is compressive throughout the entire geometry, with the vertical stress being zero at the top surface and increasing linearly with depth.

Crevasse propagation
We next consider a grounded glacier with an isolated surface crevasse, represented by an initial rectangular notch of height d s = 2.5 m and width b = 10 m, which is located at mid-length of the top surface. This facilitates comparisons with LEFM. Following Ref. [40], we also consider a damage threshold F th , below which D d = 0. As discussed in Ref. [40] and shown below, this threshold has no influence on the final crevasse depth predicted but assists in localising damage.

Linear Elastic Rheology
We first consider a linear elastic rheology for the grounded glacier, so as to validate model predictions with those obtained using analytical LEFM methods. The computational predictions of normalised crevasse depth versus time are shown in Fig. 6a for an ocean-water height of h w = 0.5H and selected values of the meltwater depth. It can be seen that the crevasses propagate rapidly and stabilise at a constant depth. In agreement with expectations, larger meltwater depths lead to higher crevasse depths, with the crevasse propagating all the way to the base of the glacier for h s /d s > 0.5. A plot of the normalised stable crevasse depths for both the analytical LEFM and phase field models is given in Fig. 6b as a function of the meltwater and ocean-water depth ratios.
The stabilised crevasse depths estimated with the phase field model show a very good agreement with those predicted using LEFM for all values of meltwater depth ratio and ocean-water height.
It can be seen that land terminating glaciers (h w = 0) are susceptible to deeper fractures, even without the presence of meltwater, as there is no ocean-water compressive pressure at the terminus.
The crevasse depth reduces significantly when ocean-water is present. For example, a dry crevasse is predicted to propagate to 37.8% of the glacier height for an ocean-water depth of h w = 0.5H. The process of crevasse growth is shown in Fig. 7, through plots of phase field φ contours at

Parametric analysis
We shall now conduct sensitivity studies on relevant material, fracture and numerical parameters. The base model considered here is an isolated dry surface crevasse with an ocean-water level h w = 0.5H. We consider the individual effect on the stabilised crevasse depth of the mode I critical fracture stress or cohesive strength σ c , the crack driving force threshold F th , the post peak slope parameter ζ, and the phase field length scale , whilst keeping all other parameters constant. The results obtained are shown in Fig. 8.  Fig. 8c. A small influence is observed, with higher ζ values leading to larger crevasse depths, as they result in a higher D d magnitude for the same stress level.
This is also consistent with the sharper drop in the uniaxial stress-strain curve with increasing ζ shown in Fig. 3. Finally, the sensitivity to the phase field length scale is explored in Fig. 8d.
The results confirm the rather negligible sensitivity of the phase field formulation employed to the magnitude of .

Non-linear viscous rheology
We next investigate the influence of the rheology upon the final crevasse depth by considering with increasing values of meltwater depth ratio h s /d s and for an ocean-water height of h w = 0.5H.
A comparison between the stabilised crevasse depths from the phase field model and LEFM can be found in Fig. 9b. The influence of meltwater within the crevasse is qualitatively similar to the linear elastic case, with stabilised crevasse depths becoming progressively larger with increased meltwater. Full fracture is predicted at a meltwater depth ratio h s /d s = 0.5 or larger. Poisson's ratio was found to have a notable influence on the fracture driving force for elastic ice sheets.

Propagation of multiple surface crevasses in a grounded marine-terminating glacier
We next determine the penetration depths for a uniform field of densely spaced surface crevasses.
The same glacier geometry from the previous example is used, but we consider seven surface crevasses, each spaced at 50 m apart and located sufficiently far away from the glacier terminus, so that the edge effects do not influence crevasse growth (see Fig. 10). Here, we aim to study the effect of neighbouring crevasses, which are expected to provide crack shielding that reduces the final crevasse depth, and to compare the phase field model results with those predicted by the Nye zero stress model [14]. The results from the Nye zero stress model are found by computing the depth at which the far field longitudinal stress becomes zero, represented by the dashed purple line in Fig. 10. The model uses approximately 1.6 million linear triangular elements, with the mesh being refined ahead of each crevasse. Plots of the phase field damage variable can be found in Fig. 11, for an ocean-water height of h w = 0.5H and a meltwater depth ratio of h s /d s = 0.1. Qualitatively, the behaviour resembles that of the single crevasse model -crevasses propagate rapidly and subsequently arrest upon reaching the compressive region at the bottom. Each crevasse stabilises to a similar depth, although the outer crevasses penetrate slightly deeper because they experience shielding only from one side.

Interaction Between Surface and Basal Crevasses
In a floating ice shelf, iceberg calving can occur when the combined depth of surface and basal crevasses at a location reaches the full ice thickness [60]. Therefore, we consider the propagation of a surface and a basal crevasse within close proximity of each other and near the calving front.  The quantitative output of the calculations is shown in Fig. 17. Consider first Fig. 17a, where the predictions of crevasse depth are shown for the surface crevasse, as well as for the basal crevasse in isolation and at selected separation distances from the surface crevasse. First, a comparison with Fig. 15 (for h s /d s = 0.8) shows that the extent of surface crevasse penetration is the same with and without the presence of a basal crevasse. This is unlike the basal crevasse, which exhibits a stabilised crevasse depth that it is very sensitive to the proximity of a surface crevasse. As shown in Fig. 17a, the stabilised crevasse depth increases with the distance to the surface crevasse, with the limit case being given by the result obtained in the absence of a surface crevasse. The combined basal and surface crevasse depth is shown in Fig. 17b. It is interesting to note that the growth of the basal crevasse is hindered by the presence of the surface crevasse when they are aligned, and consequently calving is not observed. Also, since basal and surface crevasses do not coalescence, their combined depth exceeds the glacier height for sufficiently large separations. For basal crevasses directly beneath the surface crevasse, the crevasse propagates to 37.7% of the ice shelf depth, compared with 80.6% for the isolated basal crevasse.

Nucleation and growth of crevasses: application to the Helheim glacier
In this section, we simulate the initiation and propagation of crevasses from arbitrary sites in the Helheim glacier, one of the largest outlet glaciers in southeast Greenland. The aim is to show how the creep analysis can be used to determine the nucleation of crevasses, which are then predicted to grow in a coupled deformation-fracture simulation. To generate the glacier geometry, we take the surface elevation and basal topography data from field observations (see Refs. [54,62]).
A free slip boundary condition is applied normal to the base and the inlet flow velocity is restrained to zero at the left edge. Also, we apply an oceanwater pressure at the glacier terminus and assume an ocean water height of h w = 0.85H. The geometry is discretised using approximately 140,000 triangular quadratic plane strain elements.  Damage evolution is subsequently predicted using the phase field model with the updated geometry, assuming non-linear viscous ice rheology. As shown in Fig. 20c, we find that a field of densely spaced surface crevasses can initiate at sites both close to and away from the calving front. However, the depth at which they propagate to is shallow in comparison to the glacier geometry (approximately 40 m deep). This is in agreement with the field observations of Mottram and Benn [63], who measured crevasse depths close to the calving front of Breiðamerkurjökull in Iceland, finding that crevasses only penetrated tens of meters in depth. At the calving front, we also observe that damage can propagate to the full depth of the glacier, illustrating the possibility of ice cliff failure and retreat of the grounding line. This case study showcases the ability of the computational framework developed to combine creep and damage modelling to predict both the nucleation of crevasses and the subsequent propagation, for realistic geometries and conditions. propagation region is chosen to be at least 5 times smaller than and the model is discretised using 1.5M linear tetrahedral elements. The results obtained are shown in Fig. 22, through contours of the phase field variable in the fully damaged regime (φ = 1). Initially, the two crevasses propagate vertically (along the z-direction) and horizontally (along the y-direction). Subsequently, as the two crevasses approach each other, the crack tip stress state becomes mixed mode and this results in the two crevasses curving away from each other. This is followed by the development of a hook-shaped geometry before their coalescence. Similar fracture patterns have been observed in geological faults, with remote sections of the fault growing as purely tensile fractures, whilst in close proximity to each other the faults grow as mixed mode fractures [64,65]. This behavior has also been observed in laboratory experiments [61].

Conclusions
We have presented a new stress-based poro-damage phase field model for predicting hydrofrac- using Robin boundary conditions. Nucleation and growth of crevasses in a realistic geometry, that of the Helheim glacier, is predicted in the fourth case study, combining a sequential creep-damage analysis. Finally, the last case study provides the first simulation of interacting crevasses in 3D ice sheets. Several conclusions can be obtained from the model's insight into these case studies: • The model adequately predicts the propagation of crevasses in regions where the net longitudinal stress is tensile, without the need for ad hoc fracture driving force decompositions and exhibiting very little sensitivity to the choice of phase field length scale .
• Model predictions provide a good agreement with LEFM and Nye's zero stress model when particularised to the idealised conditions where these analytical approaches are relevant.
• Increasing amounts of meltwater, as a result of climate change, can significantly enhance crevasse propagation, with iceberg calving being predicted for meltwater depth ratios of 50% or larger.
• Predicted crevasse depths are greater when considering the incompressible stress state intrinsic to a non-linear viscous rheology. Thus, first-order estimates obtained from analytical LEFM approaches should consider a Poisson's ratio of ν = 0.5 to avoid underpredicting the impact of meltwater on ice-sheet stability.
• The model captures how the presence of neighbouring surface crevasses provides a shielding effect on the stress concentration and reduces the predicted crevasse depth.
• The model accurately predicts the growth of surface crevasses within floating ice shelves near the shelf front for large meltwater depth ratios. Also, if a surface crevasse is in close proximity to a basal crevasse then a reduction in basal crevasse penetration depth is observed.
• Crevasses are predicted to nucleate in areas with high surface gradients, highlighting the need for an adequate characterisation of the glacier's geometry.
• The large-scale 3D analyses conducted demonstrate the capabilities of the model of opening new horizons in the modelling of crevasse growth phenomena under the computationallydemanding conditions relevant to iceberg calving.
Potential future extensions of the present computational framework include incorporating basal melting, lateral and basal friction effects, and ice refreezing. We offer this novel approach as a means to capture the process of crevassing and calving within ice sheets and ice shelves, to better capture these processes in efforts to prognostically assess ice-sheet vulnerability to ice shelf stability and the resulting accelerated ice sheet flow to the ocean and sea level rise, and/or grounding line retreat (potentially driven by calving at a marine-terminating margin).

Acknowledgments
T. Clayton acknowledges financial support from the Natural Environment Research Coun- The equations of linear elasticity are then used (along with the plane strain assumption) to find the out of plane normal stress σ yy in relation to the in-plane normal stresses σ xx and σ zz : Setting ε yy = 0 and rearranging to find σ yy gives: σ yy = ν(σ xx + σ zz ) (A.12) Substituting this into the into the longitudinal strain equation gives: The membrane stress assumption is then adopted due to the thickness of glaciers being several orders of magnitude smaller than the length. The horizontal displacement is therefore invariant with depth, leading to the the following derivative: ∂ε xx ∂z = 0 (A.14) Applying this constraint to Eq. A.13 and rearranging in terms of the derivative of horizontal Since the longitudinal stress is invariant with x-coordinate and with the plane strain assumption, the longitudinal stress is only variant on the z-coordinate.
where C 1 is the indefinite integration constant that can be determined by considering the force equilibrium in the longitudinal direction for the lithostatic force of ice and the hydrostatic force of the ocean water F w = 1 2 ρ s gh 2 w as intensity factor is equal to the fracture toughness, the crack will propagate in an unstable manner; however, as the crack penetrates to greater depths (where the longitudinal stress reduces) the stress intensity factor decreases and the crack will arrest when K net I becomes less than K IC . The stress intensity factor is calculated using Eq. B.1 and is integrated over the entire crevasse depth due to the driving stress (far field longitudinal stress) varying linearly with depth. The use of σ net allows us to incorporate the contributions of the ice self weight, the ocean-water pressure and the meltwater pressure into the stress intensity factor. In Eq. B.1, M D is a weight function that is dependent upon the boundary conditions and specimen geometry. Owing to the boundary condition differences, the appropriate weight functions for the grounded glacier and the floating condition cases are different. For the grounded glacier condition, the 'double edge cracks' formulation gives good agreement with the stress intensity factors calculated using the displacement correlation method within FEM [52]. The weight function for the double edge cracks model is given by [67] M D = .

(B.5)
For the floating ice shelf condition, the stress intensity factors calculated using the weight function method in Krug et al. [54] and van der Veen [18] give better agreement with the stress