A Transient Coupled Ice Flow‐Damage Model to Simulate Iceberg Calving From Tidewater Outlet Glaciers

Iceberg calving, the detachment of an ice block at the glacier front, is the main process responsible for the dynamic mass loss from the ice sheets to the ocean. Understanding this process is essential to accurately predict ice sheet response to the future climate. We present a transient multiphysics finite‐element model to simulate iceberg break‐off and geometry evolution of a marine‐terminating glacier. The model solves the coupled equations of ice flow, damage mechanics, oceanic melt, and geometry evolution on the same Lagrangian computational grid. A modeling sensitivity analysis shows that the choice of stress measure used for damage evolution strongly influences the resulting calving front geometries. Our analysis suggests that the von Mises stress measures produce the most realistic calving front geometry evolutions for tidewater glaciers. Submarine frontal melt is shown to have a strong impact on the calving front geometry. The presented multiphysics model includes all processes thus far shown to be relevant for the evolution of tidewater glaciers and can be readily adapted for 3‐D and arbitrary bedrock geometries.


Introduction
Iceberg calving, that is, the mechanical loss of ice from glaciers and ice shelves, is one of the main contributors to sea level rise from tidewater glaciers together with surface melt (van den Broeke et al., 2016) and is expected to further increase in the future (Nick et al., 2013). Consequently, it is crucial to understand the calving process to accurately predict the ice sheets' response to future climates (Benn, Warren, et al., 2007).
Iceberg calving is a dynamical process of material failure. When the local stress field near the the calving front reaches the strength threshold of ice, the material begins to weaken and ultimately fails. This drives the formation and propagation of microcracks that can coalesce and lead to the formation of crevasses at the surface and the bottom of the ice. Intense crevassing eventually leads to the detachment of blocks of ice from the glacier front, termed iceberg calving. The break-off of chunks of ice can occur above and below the water line, and, if basal and surface crevasses intersect, blocks of the entire thickness of the glacier can detach.
The local geometry of the terminus controls the stress field and therefore the fracture and failure process (Hanson & Hooke, 2000;Mercenier et al., 2018). Further processes can contribute to iceberg calving. When calving losses above the water line exceed those below, hydrostatic forces lead to a buoyancy-induced upward rotation of the ice below water, which subsequently detaches from the calving front. On the other hand, subaqueous melt erosion of the glacier terminus results in the formation of overhanging ice above the water line and therefore increased stresses and calving activity (Benn et al., 2017;Motyka et al., 2013;O'Leary & Christoffersen, 2013).
The majority of modeling studies on iceberg calving have used the zero stress crevasse-depth model (Benn, Hulton, et al., 2007, Benn, Warren, et al., 2007, which sets the terminus position at the location where crevasses penetrate below the water level. In this approach, the depth of a crevasse is determined based on the difference between the overburden ice pressure and the tensile stress. This dynamic approach for calving allowed for successful reproduction of calving front variations of ocean-terminating glaciers in Greenland and Antarctica (Cook et al., 2014;Otero et al., 2010Otero et al., , 2017Nick et al., 2010Nick et al., , 2013Todd et al., 2018). Although the crevasse-depth model can be calibrated to observations (Lea et al., 2014), the validation of the related processes with field observations is limited and is based on a snapshot of the stress balance, neglecting the preexistence of cracks and their effect on the stress state of the glacier (Krug et al., 2014;Mercenier et al., 2018). The earlier studies applying this model were limited to one or two dimensions and did not permit the model domain to evolve through time.
The application of fracture mechanics for crevasse penetration was introduced to represent the rapid propagation of preexisting fractures at very short time scales (Smith, 1976). In linear elastic fracture mechanics (LEFM) models, crevasses penetrate downward as long as a stress intensity factor exceeds the fracture toughness of ice (Rist et al., 1999;Smith, 1976). The LEFM approach can be expanded for two-dimensional plan-view crevasse modeling and permits the combination of different modes of crevassing (Colgan et al., 2016;Van der Veen, 1999). Due to the boundary-condition dependence of the stress intensity factor, it is important to note that calving models using LEFM approaches (Krug et al., 2014;Van der Veen, 1998b, 1998a) may be appropriate for floating ice tongues, but not for grounded glaciers (Jiménez & Duddu, 2018).
Recent more sophisticated approaches overcame some of the limitations of the zero stress and LEFM models. Benn et al. (2017) combined a flow model with a discrete element model to predict the calving front position, whereas Todd et al. (2018) developed a three-dimensional ice flow-calving model based on a modified crevasse-depth relation, including a rediscretization scheme. Another recent approach relies on the level set method in a two-dimensional map view calving model for which the calving front retreats when the von Mises tensile stress reaches a stress threshold (Morlighem et al., 2016, cf. section 2.2.2 for definition). This approach was used to model Store Gletscher's response to ocean thermal forcing. Further, Choi et al. (2018) compared four calving laws to determine calving front positions for several Greenland outlet glaciers. They found that the von Mises tensile stress calving law is best suited to simulate observed ice front positions.
As an alternative, a continuum description of microcrack accumulation is provided by isotropic damage mechanics. While other methods focus on individual crevasses in the ice, this method describes the presence of defects within the ice and their effect on the ice rheology . Pralong et al. (2003) first proposed an Eulerian damage mechanics description combining Stokes flow and continuum damage mechanics and were able to successfully predict lamella break-off from hanging glaciers (Pralong & Funk, 2005). Further glaciological studies used the combination of continuum damage mechanics with ice flow models (Borstad et al., 2012;Duddu & Waisman, 2012Jiménez et al., 2017;Jouvet et al., 2011;Krug et al., 2014Krug et al., , 2015Mercenier et al., 2018;Mobasher et al., 2016;Pralong & Funk, 2005). Applying this technique in an Eulerian reference frame facilitates the description of the physics of ice flow and fracture; however, the implementation of damage advection with flow is affected by artificial diffusion, and the numerical accuracy is mesh dependent (Jiménez et al., 2017). To overcome these issues, the use of a Lagrangian reference frame would be more appropriate to study fracture and damage mechanics (Duddu & Waisman, 2012;Jiménez et al., 2017).
In this paper, we develop a transient multiphysics finite-element model that simulates ice break-off and oceanic melt at the calving front and tracks the evolving geometry. The proposed calving model is based on the coupling of continuum damage mechanics with ice flow equations, and with a computational domain that evolves through time. A coupled Euler-Lagrange formulation is used to perform the model evolution without the need for rediscretization. Elements for which the damage variable or the accumulated melt exceeds a threshold value are removed from the mesh, which allows us to effectively model ice break-off at the calving front. We provide a complete description of the model and its implementation. We then investigate the sensitivity of the model to the choice of stress measure and to variations in model parameters such as damage rate, critical damage, ice thickness, and oceanic melt rate. These model experiments are performed for an initially idealized glacier geometry that evolves with ice flow, oceanic melt, and calving.

Methods
We use the finite-element library libMesh (Kirk et al., 2006) to implement a transient coupled ice flow-damage model. The model is formulated as a moving-mesh algorithm outlined in the flow chart of Figure 1.
The computational domain (the glacier) is discretized with a finite-element mesh. The continuum Stokes flow equations with strain-rate-dependent viscosity in an Eulerian reference frame are solved on this mesh at each time step. The resulting velocity field, defined on the mesh nodes, is used to calculate the evolution of a scalar damage variable, and for advection of the mesh nodes. An elimination algorithm deletes mesh elements that are considered destroyed and/or detached, representing glacier calving and oceanic melt. In all presented model experiments an idealized glacier geometry was used that initially consisted of a rectangular block of ice resting on a flat bed with an initial ice thickness H = 200 m and length L=10 H. The model domain was defined in a Cartesian coordinate system with horizontal axis x and vertical axis z with origin at sea level. The idealized glacier was partially submerged in oceanic water with a relative water level = H w H = 0.75 (H w is the water height) for all experiments.
The mesh size dependence of the model was tested and showed limited differences in total volume losses for the different mesh sizes ( Figure  S1 in the supporting information). It is also important to note that for the various mesh sizes, the modeled geometries displayed similar frontal shapes ( Figure S2). In order to reduce the computational time, the initial mesh was chosen with 20 elements in the vertical and 200 elements in the horizontal for the sensitivity analyses. The time step dependence of the model was tested and also showed limited differences in total volume losses ( Figure S3). All model runs were performed for 2,000 time steps of t=0.001 year corresponding to 2 years of geometry evolution.
First, we describe the ice flow model equations and its boundary conditions. Then we present the continuum damage mechanics model and its coupling to the ice flow model. Finally, we explain the geometry evolution and the calving and melting algorithms.

Ice Flow Model 2.1.1. Ice Flow and Rheology
Ice flow was computed at each time step on the model domain in an Eulerian reference frame. To model ice flow, we implemented the Stokes equations for incompressible fluid flow with continuum conservation of momentum and mass: where is the Cauchy stress tensor, i the ice density, g the gravitational force vector, and u the velocity vector (Table 1). For isotropic and incompressible ice the Cauchy stress tensor can be decomposed into an isotropic and a deviatoric part ′ : where m = 1 3 tr( ) = 1 3 ii is the mean stress and I the identity matrix. The undamaged ice rheology is described as a viscous power law fluid (Glen's flow law), linking the deviatoric stress tensor ′ to the strain rate tensor . : The effective shear viscosity is defined as where ) 1 2 is the effective strain rate, A the fluidity parameter, n=3 the power law exponent, and is a regularization parameter that imposes a linear rheology at small strain rates to avoid infinite viscosity at low stresses (Greve & Blatter, 2009).
The model domain was discretized with second-order square, isoparametric nine-node quadrilateral elements with Galerkin weighting. The model variables were approximated with a second-order approximation for the velocities u and v and a first-order approximation for the mean stress m (forming a LBB-stable set).
For each time step, the Stokes equations for ice flow with nonlinear rheology were assembled within libMesh and solved with the PETSc nonlinear solver SNES (Balay et al., 2017) to a relative accuracy of 10 −4 .

Boundary Conditions
In the Eulerian ice flow model the upper surface of the glacier was described as a traction free boundary. The upstream boundary was treated as a zero velocity Dirichlet boundary condition, independently of the inflow of ice prescribed for the geometry evolution (section 2.4), and without effect on the stress field at the terminus.
At the ice/ocean interface a normal stress boundary condition was imposed below the water level (z<0), while the surface above water was kept traction-free. The stress boundary condition thus reads nn = min( w gz, 0), where nn and nt i are the normal and tangential tractions applied on the calving front and w is water density ( Table 1). The traction nn is compressive (i.e., nn <0) since z<0 below water.
On the lower boundary of the domain the type of boundary condition depends on the normal stress nn =n·( ·n). For parts of the boundary where the normal stress exceeds the water pressure ( nn < w gz), a zero velocity Dirichlet boundary condition was imposed. Additional experiments with a zero velocity Dirichlet boundary condition normal to the bed and a prescribed velocity along the bed showed a negligible effect of nonzero velocity boundary conditions on the stress state. On the remaining parts of the bottom boundary the normal stress boundary condition (equation (6)) was imposed, which means that water pressure is applied. The latter condition assumes that there are no cohesive forces restraining the ice movement away from the bed (Durand et al., 2009).

Continuum Damage Mechanics Model 2.2.1. Damage Evolution
Continuum damage mechanics is based on the introduction of a damage variable D. In this study we assume isotropic damage using a scalar variable that represents degradation of mechanical properties at the mesoscale, which averages the accumulation of defects present in the material at the microscale (Lemaitre, 1992). Following the hypothesis of strain equivalence, an effective stress is introduced: where D is a state variable that takes the value D=0 in undamaged ice and D=1 when full damage is reached. Damage evolution depends on a stress measure and begins once a stress threshold th ( th = 0.11 MPa; Table 2) is crossed: Damage evolution is parametrized with the damage rate B, the stress measure chosen for damage evolution (section 2.2.2), and the power r=0.43 chosen according to Pralong and Funk (2005). A healing term h is easily added, but not considered in this study (i.e., h=0).
Damage evolution on all mesh nodes was calculated by evaluating the stress measure based on the solution of the ice flow model at each time step (Figure 1). In our model implementation damage is a state variable on the nodes of the moving mesh, therefore bypassing the use of an advection scheme that could introduce numerical diffusion of damage.

Stress Measure
The stress measure for damage evolution should be independent of the choice of coordinate system and can therefore be expressed as a function of the invariants and eigenvalues of the stress tensor. Our model permits the use of arbitrary linear combinations, and results from a selection are shown in the following.
Here we list the commonly used stress criteria in glaciological studies (see Table 2 for a summary of the stress measures and their notations): • The Hayhurst criterion H is the most widely used stress criterion in glaciological studies for damage evolution (Duddu & Waisman, 2012Jiménez et al., 2017;Mercenier et al., 2018;Mobasher et al., 2016;Pralong & Funk, 2005;Pralong et al., 2003). This criterion (Hayhurst, 1972) uses a linear combination of the maximum principal stress 1 , the first stress invariant I 1 = m = 1 3 ii , and the von Mises stress 1 n e as stress measure : The Hayhurst stress is generally used as a criterion for ductile and brittle materials that allows for damage under uniaxial compression (Krug et al., 2014;Pralong & Funk, 2005). • The maximum principal stress 1 has been used as a stress measure in combination with LEFM to simulate calving (Krug et al., 2014). As a stress criterion, the maximum principal stress only represents the purely tensile contribution and is therefore useful to locate the position of surface crevasses (Benn et al., 2017;Mercenier et al., 2018). However, without supplementary forcing such as water pressure in cracks, the stresses are compressive close to the bed due to the increasing overburden ice pressure (Krug et al., 2014;Ma et al., 2017) and are likely unable to induce behaviors such as basal crevassing. Combining the maximum principal stress distribution with the effect of water pressure in the submerged parts of the glacier front is suggested to lead to the formation of basal crevasses (Benn et al., 2017). This stress combination, namely, the effective principal stress 1,w = 1 +P w , is also considered in our stress measure analysis (P w = max(− w gz, 0) is the water pressure with z<0 below water). • The most widely used stress criterion in solid mechanics and structural analysis is the von Mises criterion J 2 = e , which is a natural quantity to investigate yielding of ductile materials. In glaciology, this criterion has been used to simulate crevasses and fractures on ice shelves (Albrecht & Levermann, 2014;Vaughan, 1993). • Since calving occurs mainly through stretching (Benn, Warren, et al., 2007), the von Mises tensile stress t is additionally considered as a stress measure for damage evolution (Morlighem et al., 2016). The von Mises tensile stress only takes into account the tensile parts of the von Mises stress and neglects compression. To formulate the von Mises tensile stress, an effective tensile strain ratẽ. e is defined: where .
i is the eigenvalues of the strain rate tensor. The von Mises tensile stress is then written as

Viscosity Feedback of Damage
The accumulation of microdefects weakens the ice through damage, which leads to increased deformation rates (Krug et al., 2014;Pralong & Funk, 2005;Pralong et al., 2003). Thus, a damage dependence in the material rheology is introduced, using the equivalence principle of equation (7) for isotropic damage, to describe the rheology of damaged ice (Pralong & Funk, 2005). Practically, the effective shear viscosity was adjusted such that This equation effectively replaces equation (5) in the ice flow model. The stress and flow regimes for undamaged ice (D=0) remain unaffected while damaged ice is softened, which leads to increased ice deformation rates. Introducing damage thus leads to a positive feedback on the stress and velocity fields. This feedback is represented in the model flowchart ( Figure 1) by a direct link between the damage output D box and the ice flow model box.

Submarine Melting
Melting of the ice in the submerged parts of the calving front, due to buoyancy-induced heat advection of oceanic water (Carr et al., 2013;Howat et al., 2010;Motyka et al., 2013;, is suggested to produce an overhanging ice face below the water line leading to increased calving activity (Benn et al., 2017;Krug et al., 2015;O'Leary & Christoffersen, 2013).
Submarine melting was implemented using a per element state variable m at the boundary between ice and ocean water. In every such element the melt from each side i with length L i exposed to the oceanic water was accumulated following where M(x,z,t) was the submarine melt rate imposed, which is an arbitrary function of space and time. If the total accumulated volumetric melt m exceeded the initial element volume V i , the element was removed from the mesh. The excess melt of this element was then evenly distributed among its neighbors. The choice of this implementation method for submarine melt is discussed in section 4.3. The updated geometry after the elements with D>D c were removed from the mesh. Note that elements disconnected from the main part of the mesh on the left side of the damaged area were also removed from the mesh.

Geometry Evolution and Calving
At each time step, the velocities u calculated by the ice flow model were used to update the model geometry.
Additionally, a constant horizontal velocity u in was prescribed that corresponds to the upstream inflow of ice into the computational domain and also served as a uniform basal motion. In the presented results, the ice moves from right to left toward negative x coordinates.
The mesh was modified with the following consecutive steps: 1. Elements in which the damage variable exceeded a threshold D c were removed from the mesh, as were elements no longer connected to the main mesh ( Figure 2). Details are given below. 2. Elements under the water line for which m>V i were removed from the mesh (section 2.3). 3. All mesh nodes were moved by the modeled velocity field u t. 4. All mesh nodes were additionally displaced by u in t with u in <0. 5. If the upstream end of the mesh had moved by more than one element length from its original position, a vertical column of elements was added at the upstream end of the domain. The damage variable on these elements was initialized to D=0. 6. The new mesh boundary was scanned, and appropriate boundary conditions for the flow model were set.
Iceberg calving was simulated by introducing a critical damage above which the ice was considered too weak. All elements that met the condition D>D c were removed from the mesh. This element deletion method is a heuristic approach that allows us to simulate the evolution of the frontal boundary in a simple and efficient way. After this step, groups of elements might exist that are no more connected to the main mesh, as illustrated in Figure 2. A simple search algorithm iteratively traversed the mesh and tagged elements that are connected to either the bed or to a tagged element. The search stopped when no more elements were tagged, and untagged elements were removed from the mesh. All elements removed due to this algorithm were considered as the ice volume calved during a time step. The calved and melted element volumes were stored for further processing. Note that accumulation/ablation fluxes at the upper surface of the glacier are not accounted for in the geometry evolution.

Results
In this section, we present the results of the different model sensitivity analyses and the initial adjustment of the calving front. The model formulation allows us to investigate the effect of different stress measures, as well as variations of damage rate B, critical damage D c , the initial ice thickness, and submarine melt. As a reference experiment, a model run with the von Mises stress e as stress measure, damage rate B=5 MPa −r a −1 , critical damage D c =0.5, initial ice thickness H=200 m, and a relative water depth =0.75 was chosen.
Note that for the analysis of the sensitivity experiments, our interest is to obtain calving front geometries that are representative of tidewater glacier geometries observed in nature. Calibration of the model parameters with observations would then likely allow us to reproduce the advance and retreat of the calving front positions.
For all experiments presented, the initially rectangular calving front adjusts within half a year to a quasi-stable shape. Depending on the model parameters, mainly u in , the front then keeps slightly advancing or retreating at a relatively constant rate (see Movies S1-S5). Various influx velocities u in were chosen for the different experiments to compensate frontal mass loss and therefore constraining the domain lengths for illustration purposes. It is important to note that the influx velocities do not influence the shape or mass loss evolution of the modeled idealized glaciers.

Choice of Stress Measure
To investigate the influence of different stress measures on glacier geometry, the coupled model was run for 1,000 time steps (1 year) starting from the same geometry the reference parameters (Table 3). Figure 3 and Movie S1 show the geometries for five different stress measures. For the experiments using the von Mises stress e and von Mises tensile stress t an influx velocity of u in = 2,500 m a −1 was chosen, whereas u in = 0 m a −1 was used for the other stress measures. For the Hayhurst stress, the linear combination parameters ( , , from equation (9)) were chosen according to Pralong and Funk (2005).
The results show distinct model geometry adjustments for the different stress measures (Figure 3 and Movie S1). Using the von Mises e and von Mises tensile t stress, similar calving front geometry evolutions and final shapes were obtained. For these stress measures the damage variable first evolves before a big crevasse forms that penetrates the whole thickness of the glacier at the calving front. The model geometry further evolves and forms an ice foot that regularly breaks off below water. This pattern is particularly pronounced for the von Mises tensile stress where the formation and breaking of the ice foot follow a repeating cycle. Although they have the same influx velocity and damage rate, these two stress measures produce different rates of calving, resulting in a small retreat for the von Mises stress geometry and an advance with the von Mises tensile stress geometry.
For the Hayhurst H and maximum principal stress 1 distributions (Figure 3), the damage is only concentrated at the surface and does not penetrate through the whole thickness of the glacier during the model experiment. Elements are only removed at the surface, and most of the ice below the water line remains intact. In reality, however, the damaged ice would likely remain at the surface, as it cannot leave the glacier unless it breaks off at the calving front. Further, the influx velocities u in were set to 0 to constrain the length of the domain and limit computation time, as the glacier front continuously advances without breaking off throughout the simulations.
For the Hayhurst stress (Figure 3), in contrast to using the maximum principal stress, the rapidly forming floating tongue is detached further away from the bed, and more ice is removed at the top surface. Further, if we take the water pressure in the submerged parts of the calving front into account using the effective principal stress measure 1,w , the elements at the front and below the floating tongue are progressively removed from the mesh during the simulations. However, damage does not penetrate through the whole thickness upstream of the glacier tongue. In this study we aim at providing a model with realistic tidewater glacier behavior such as regular calving and realistic terminus geometries. For the purpose of illustrating the model capabilities, we therefore restrict ourselves to the von Mises stress measure e . Only future analysis of the sparse field data with this model will allow us to identify the most suitable stress measures together with model parameters. Figure 4 and Movie S2 illustrate the quasi-stable geometries obtained from the coupled model after 1,000 time steps (1 year) with a critical damage D c ≃1 (with u in = 2,000 m a −1 ) and D c =0.5 (with u in = 2,500 m a −1 ). The model parameters for these experiments are shown in Table 4.

Influence of Critical Damage
Increasing the critical damage leads to a slower mass loss due to ice break-off at the calving front (Figures 4  and S4; Movie S2). Further, the resulting geometries differ slightly, with a more pronounced buoyancy driven up-rise for the fully damaged geometry (D c ≃1) and an ice foot closer to the bed for the lower critical damage (D c =0.5).   Figure 5 and Movie S3 show the geometries obtained from the coupled model after 1,000 time steps (1 year) with different damage rates B. Table 5 shows the model parameters that were chosen for all experiments.

Influence of Damage Rate
Varying the damage rate appears to rather change the position of the calving front, as a result of different calving rates, than the geometrical properties of the calving front (Figures 5 and S5; Movie S3). The model experiments follow the same geometrical evolution of the calving front, with a certain time lag due to the different damage rates. All geometries develop the distinctive inclined calving front with an ice foot developing below the water level. Thus, increasing the damage rate leads to a long-term higher rate of mass loss due to enhanced ice break-off at the calving front, but the effect on the calving front geometry is limited.

Influence of Ice Thickness
The stability of the glacier front strongly depends on the ice thickness and water depth (Bassis & Walker, 2012;Brown et al., 1982;Mercenier et al., 2018). Therefore, we performed sensitivity tests on the ice thickness of the initial mesh. The element size was kept identical for the different experiments, with initial square element dimensions of 10 m·10 m. For the ice thicknesses of H=250 m, H=200 m, and H=150 m, the influx velocities were set to 5,000 ma −1 , 2,500 ma −1 , and 1,000 ma −1 , respectively. The initial length L of the geometries was set to L=10 H. Table 6 shows the model parameters selected for the ice thickness sensitivity experiments.
In a similar way as for the damage rate sensitivity analysis, increasing the initial ice thickness leads to a higher rate of mass loss due to ice break-off at the calving front, with a limited effect on the calving front geometry ( Figure 6 and Movie S4). Though it appears in Figure 6 and Movie S4 that the thicker geometry loses mass at a lower rate, note that the constant horizontal velocities imposed strongly differ. The thicker glacier actually starts to lose mass earlier than the thinner glaciers and continues to lose mass at a higher rate ( Figure S6). Similar geometries are developed for all ice thicknesses, with a certain delay due to the different rates of mass loss.
For the purpose of this model description study, the thicknesses were limited to a maximum of 250 m in our modeling effort. The relative water depth was kept constant for the ice thickness experiments. However, further experiments using various relative water depths showed similar geometry evolutions as for the different ice thicknesses ( Figure S7), with a higher rate of frontal mass loss for smaller values of .   Figure 7 and Movie S5 show the geometries obtained from the coupled model after 1,000 time steps (1 year) for different melt scenarios. The reference parameters for all experiments are displayed in Table 7. A constant melt rate M, independent of time and depth, is prescribed below the water line where melt is accumulated in the frontal elements following equation (13). More realistic melt parametrizations could easily be prescribed.

Influence of Submarine Frontal Melt
The melt rate sensitivity analysis shows that submarine melt strongly affects the shape of the calving front ( Figure 7 and Movie S5). For the experiment with M = 1,000 m a −1 the front progressively develops a calving front that is oversteepened below the water line. During the first year of simulation, increasing the melt rate leads to a faster retreat of the calving front. However, once the calving front shape has adjusted itself after 1 year, mass loss occurs faster when submarine frontal melt is neglected ( Figure S8).

Discussion
The coupled ice flow-calving model presented in this study is capable of producing different geometries typically encountered in tidewater glaciers (Figures 3 to 7) and different calving styles by only varying a set of parameters and forcings. The dependence of model performance and glacier evolution on these parameters, and avenues for validating the model with measurements, are elaborated in the following sections.

Retreat and Advance
The calving front evolution is essentially controlled by four processes acting on different time scales: ice deformation, damage evolution, frontal melt, and basal sliding. Understanding these time scales is crucial, as they determine the rate of advance and retreat of the glacier front. Basal sliding and ice deformation determine the mass flux and thus affect the rate of advance of the glacier, while damage evolution and submarine melt control the calving rate and hence retreat.
While melt and inflow from upstream are external forcings, damage evolves at its own time scale, driven by mechanical stress and thus the glacier geometry. Damage evolution leads to two processes with different time scales: subaerial break-off and break-off below water. The latter process is an order of magnitude slower due to the effective gravity (buoyant force) under water of magnitude ( w − i ) ∕ i g ≃ 0.12 g and accordingly reduced deviatoric stresses.
The model results show that the shapes of the calving front obtained by varying the damage parameters remain very similar ( Figure 5). Thus, these parameters determine the time scales at which break-off occurs at the calving front in the coupled model rather than the shape of the calving front. The ice thickness variations exhibit a similar effect on the time scale of glacier geometry evolution ( Figure 6).  Note that for the ice thicknesses of H=250 m, H=200 m, and H=150 m, the influx velocities were set to 5,000 ma −1 , 2,500 ma −1 , and 1,000 ma −1 , respectively.
In order to reproduce observed rates of mass loss at the calving front the calibration of the critical damage and damage rate parameter are required, and the submarine melt rate would have to be prescribed.

Stress Measures and Calving Front Geometries
The choice of stress measure has a strong influence on the modeled calving front geometry. All our modeling experiments reproduce a vertical ice cliff above the water line (Figure 3), a feature that is commonly observed at calving glaciers (Benn, Warren, et al., 2007). In contrast, below the water line the modeled geometries show large differences for the different stress measures.
A noteworthy feature that is common to all model runs with von Mises stress measure is the presence of an ice foot below the water line. The development of subaqueous ice feet, which have been observed for a series of calving glaciers (Hunter & Powell, 1998;Motyka, 1997;O'Neel et al., 2007;Warren et al., 1995), can be caused by calving losses above the water line exceeding those below (Benn, Warren, et al., 2007). Buoyant forces then induce a rotation of the terminus either by ice creep or fracture propagation and calving (Benn, Warren, et al., 2007;Warren et al., 2001). This process can clearly be observed in our model experiments  involving the von Mises stress. However, this process is expected to be much slower than aerial calving since the buoyancy forces are ∼9 times smaller than gravitational forces, leading to an order-of-magnitude lower effective upward gravity and accordingly reduced stresses. Indeed, in our model subaqueous calving events are less frequent but of larger magnitude than above water (Movie S1), in agreement with observations (Hunter & Powell, 1998;Motyka, 1997;O'Neel et al., 2007;Warren et al., 1995).
Only for a few marine-terminating glaciers, the underwater geometry has been investigated. Recent studies using multibeam bathymetry in West Greenland revealed undercut ice faces below the water line as well as ice foot formations (Fried et al., 2015;Rignot et al., 2015). The observed glacier undercuts were larger at the sites of subglacial water discharge from the glacier. Buoyant plumes of subglacial water transport the warm oceanic water to the ice face and consequently enhance the local melt rate (Rignot et al., 2015;Slater et al., 2016). In the presented model runs, the persistent presence of an ice foot at the calving front could therefore be explained by the lack of undercutting by melt and the absence of preexisting cracks near the bed of the partly grounded glacier. As shown in section 3.5, undercutting by submarine melt has a significant impact on the underwater calving front geometry as the ice foot is melted away while it forms. Further, the preexistence of cracks at the bed, due to irregularities in the bed for example, would be expected to increase damaging and reduce the formation of an ice foot.
Model runs with the Hayhurst and maximum principal stress show damage concentration at the top surface rather than at the bed of the glacier. Therefore, cracks do not penetrate through the whole ice thickness, and elements are deleted only at the top surface. In reality such ice might not be removed from the glacier as it is in our model implementation (section 3.1). The lack of depth penetration is likely due to the importance of the overburden ice pressure at depth, which exceeds the tensile and shearing stress components. The no-slip basal boundary condition could further explain the lack of depth penetration, as free-slip boundary conditions led to full-thickness calving using the maximum principal stress (Ma et al., 2017).
Based on the above considerations, the von Mises and von Mises tensile stresses produce the most realistic calving front geometry evolutions for tidewater glaciers using a continuum damage mechanics model. This observation agrees with results from Choi et al. (2018), who found that the von Mises tensile stress provides more satisfactory results than the three other calving laws in their comparison effort.
Combining the maximum principal stress with the effect of water pressure in the submerged parts of the glacier, that is, using the 1,w stress measure, led to the formation of a floating ice tongue, but cracks at the surface and the bed were still not able to penetrate through the whole thickness of the glacier (Figure 3 and Movie S1). However, the presence of surface water would likely facilitate crevasse propagation through the whole thickness of the glacier (Benn, Warren, et al., 2007). Including the effect of such a process into the model with these stress measures would likely be effective in reproducing ice shelf disintegration.

Submarine Melt
The submarine melt parametrization used in this study is simple and independent of time and depth. However, any arbitrary melt parametrization as function of position and time can be prescribed in the presented model.
Implementing submarine melt at the calving front in a Lagrangian reference frame without a rediscretization scheme presents some challenges. While the implementation of oceanic melt used in this study is straightforward, it neglects certain effects caused by submarine melting. For instance, small geometrical changes of the calving front due to melt are neglected until the elements at the water boundary are removed. Conversely, a more accurate description of melt would first lead to a reduction of the size of the elements in contact with water before they are removed. Therefore, the length of the sides in contact with water are likely underestimated during the melt accumulation process, leading to an underestimation of melt with respect to the imposed melt rate. Further, the effect of these progressive geometrical changes on the stress state, during the melt accumulation process, at the calving front is also partially neglected. These shortcomings are mesh size dependent and negligible for small elements.
During model development, a potential improvement to the method described above was implemented but ultimately discarded. By moving the nodes in contact with the water according to the imposed melt at each time step, the size reduction of the elements in contact with water was tracked. This approach presents, however, some implementation issues. First, the solution of the damage distribution has to be remapped to the updated nodal positions of the shrinking elements. Second, the element shape often became very distorted when nodes were moved due to melt. These element distortions can cause convergence issues due to the presence of negative eigenvalues of the element Jacobian matrix. To remedy this issue, the distorted elements can be removed as soon as the element quality is too low. With such an approach elements would be removed that are not fully melted, which would lead to a nonnegligible overestimation of melt. Therefore, we decided to use the simple melt accumulation implementation presented before.
Undercutting by melt below and at the water line is generally suggested to lead to an increase in calving activity, resulting from the formation of an overhang in the submerged parts of the calving front (Benn et al., 2017;How et al., 2019;Motyka et al., 2013;O'Leary & Christoffersen, 2013;Petlicki et al., 2015). Indeed, such an overhang under the water line forms for the model experiment with a melt rate of 1,000 m a −1 (Figure 7). However, the long-term calving activity appears to be larger for the experiment without submarine melt. Therefore, the effect of melt undercutting on the rate of calving is not yet clear from our experiments. A detailed investigation will be necessary to understand the effect of modeled submarine melt on the calving rate.

Model Improvements
The coupled ice flow-calving model was applied to an idealized glacier geometry only. For arbitrary tidewater glacier geometries and comparison with data several adaptations are needed. Basal motion, now implemented as simple forward movement of mesh nodes, would be implemented as a frictional boundary with potential dependence on effective water pressure (e.g., Ryser et al., 2014). Ice flux into the model domain at the upstream boundary could be adapted in a straightforward manner to include variable vertical extent. Moreover, the evolution of the glacier over arbitrary bedrock geometries could induce local damage effects that might influence the calving front patterns.
The presented model was implemented in plane-strain formulation in a vertical section along the flow line.
Consequently, important effects such as lateral stress bridging are neglected (Choi et al., 2018;Todd et al., 2018). The current model implementation is generic and can be used in three dimensions by only changing the finite element types to three-dimensional hexahedral or tetrahedral elements. Computational costs would be much higher, which is alleviated by the fully parallel implementation of the libMesh finite element library and the employed PETSc solvers.

Conclusions
We developed a transient model of coupled damage evolution and ice flow to simulate ice break-off at the calving front of a tidewater glacier. A coupled Euler-Lagrange formulation was used for the model evolution, combined with the heuristic removal of elements at the calving front representing calving and oceanic melt, and addition of elements at the upstream influx boundary. The model therefore incorporates the main processes affecting tidewater glacier evolution and can be used to investigate the relative importance of individual processes, their time scales, and their intricate coupling.
The discussed model runs all started from an idealized geometry and evolved to a wide range of geometries and calving dynamics, which resemble those observed on tidewater glaciers. The emergent modeled calving front geometries and calving rates depend delicately on the choice of the stress measure and further model parameters, and oceanic melt. For a chosen stress measure the calving rate and hence the rate of advance and retreat of the glacier depend on the damage rate and basal sliding. Our analysis suggests that the von Mises stress measures produce the most realistic calving front geometry evolutions.