Modelling the impact of submarine frontal melting and ice mélange on glacier dynamics

. Submarine melting of the calving face of tidewater glaciers and the mechanical back force applied by the ice mélange layer are two mechanisms generally proposed to explain seasonal variations at the calving front of tidewater glaciers. However, the way these processes affect the calving rate and glacier dynamics remains uncertain. In this study, we used a ﬁnite element-based model that solves the full Stokes equations to simulate the impact of these forcings on two-dimensional theoretical ﬂow line glacier conﬁgura-tions. The model, which includes calving processes, suggests that frontal melting affects the position of the terminus only slightly (less than a few hundred metres) and does not affect the multiannual glacier mass balance at all. However, the ice mélange has a greater impact on the advance and retreat cycles of the glacier front (more than several kilometres) and its consequences for the mass balance are not completely negligible, stressing the need for better characterization of forcing properties. We also show that ice mélange forcing against the calving face can mechanically prevent crevasse propagation at sea level and hence prevent calving. Results also reveal different behaviours in grounded and ﬂoating glaciers: in the


Introduction
In the context of global warming, the cryosphere's contribution to sea level rise is a major concern. Depending on the four RCP scenarios (representative concentration pathways) considered in the IPCC fifth assessment (Church et al., 2013), the sea level is predicted to rise between 0.26 and 0.82 m in 2081-2100 relative to 1986-2005. The Greenland Ice Sheet (GIS) mass loss, which was 142 ± 49 Gt a −1 on average over the past 2 decades, has increased in recent years to reach an estimated value of 263 ± 30 Gt a −1 between 2005 and 2010 (Schrama and Wouters, 2011;Shepherd et al., 2012) and 359.8 ± 28.9 Gt a −1 between April 2009 and 2012 (Khan et al., 2014). This mass loss extended over a large part of the GIS (Khan et al., 2010;Schrama and Wouters, 2011;Khan et al., 2014), which is thus becoming a major contributor to sea level rise (Cazenave, 2006;Rignot et al., 2011).
Increasing ice loss highlights the need for accurate estimations of the future mass balance, but the large discrepancies in the behaviour of Greenland's outlet glaciers make a simple mass balance extrapolation unreliable unless we understand the processes that control their dynamics (Howat et al., 2011;Seale et al., 2011). The mass loss from the GIS is the consequence of two main mechanisms: the dynamic ice discharge (through calving and frontal melting) and the negative surface mass balance. Ice discharge was estimated to represent 40 to 60 % of the total mass loss (Rignot et al., 2008;van den Broeke et al., 2009;Khan et al., 2014), corresponding to −156.3 ± 40.9 Gt a −1 . Ice discharge is therefore a significant mechanism, and the two related processes (melting and calving) not only directly affect the position of the front but also affect the forces at the front; feedback between calving processes and ice dynamics are therefore to be expected. Holland et al. (2008) hypothesized that the increased discharge in Greenland may have been triggered by an increase in the subsurface ocean temperature. This claim was sup-

Published by Copernicus Publications on behalf of the European Geosciences Union.
ported by Straneo et al. (2010), who stated that a rapid advective pathway exists between North Atlantic's oceanic variability and the margin of the ice shelf in the vicinity of Sermilik Fjord, east Greenland. The underlying process suggests that submarine frontal melting promotes the emergence of an ice block overhanging the water line, which calves rapidly due to an undercutting effect. Remote observations of east GIS glaciers revealed a correlation between variations in the position of the terminus and variations in the temperature of the ocean (Seale et al., 2011). However, melting intensity is hard to measure accurately and is usually inferred from hydrographic measurements (water velocity, temperature, and salinity) of the heat transport within the water layers. Summer melt rates vary between 1 and 17 m day −1 depending on the glacier and the associated fjord system (Motyka et al., 2003;Rignot et al., 2010;Sutherland and Straneo, 2012;Bartholomaus et al., 2013;Inall et al., 2014).
Another quantity, that of ice mélange, a heterogeneous mixture of sea ice, marine ice, blown snow, and fragments of icebergs, is suspected to play an important role in the seasonal cycles of the glacier front (Higgins, 1991;Sohn et al., 1998;Reeh et al., 2001;Khazendar and Jenkins, 2003;Fricker et al., 2005;Copland et al., 2007;Joughin et al., 2008a, b, c). Observations showed that winter freezing of the ice mélange is correlated with a decrease in the calving rate, an advance of the glacier front, and a slowing down of the ice flow (Sohn et al., 1998;Joughin et al., 2008c) and that summer decay is followed by an increase in the calving rate, a retreat of the front, and accelerated ice flow (Higgins, 1991;Copland et al., 2007). Some authors argue that ice mélange may directly resist the ice flow (Walter et al., 2012), while others suggest that it only maintains the integrity of the terminal part of the glacier (Sohn et al., 1998;Amundson et al., 2010). Thus, although variations in the position of the terminus and the existence of a layer of ice mélange layer are clearly correlated, the underlying processes that control this behaviour are still poorly understood.
Several attempts have been made to incorporate frontal melting and the ice mélange back force in ice flow models. In particular, the relation between calving and undercutting was investigated by Vieli et al. (2002), who applied a seasonal calving pattern on a simplified geometry of Hansbreen Glacier in Svalbard, assuming that calving was controlled by melting at the water line. These authors concluded that melting-driven calving only had a minor impact on glacier dynamics. On the contrary, O'Leary and Christoffersen (2013) used a fixed geometry to investigate the effect of different melting patterns on the stress field in the ice. These authors showed that undercutting can be a strong driver of calving due to the concentration of stress that occurs at the upper surface. However, recent studies using calving parameterization based on an instantaneous stress balance (Benn et al., 2007a, b;Nick et al., 2009Nick et al., , 2010) applied on two-dimensional flow line geometries of Helheim (Cook et al., 2014) and Store glaciers (Todd and Christoffersen, 2014) tempered these conclusions. According to Cook et al. (2014), when undercut, the upper surface of the glacier drops, thereby reducing tensile stress. The effect of ice mélange on glacier dynamics was analyzed by Nick et al. (2010) and Vieli and Nick (2011) using a depth-and width-integrated model combined with the calving parameterization of Nick et al. (2010). Both studies managed to reproduce the cycles of advance and retreat of the glacier fronts with realistic amplitudes. However, the authors did not undertake further investigation of the underlying processes. Cook et al. (2014) stated that only unrealistically high back pressure would be able to change the position of the front, highlighting the need for further modelling focused on processes, whereas Todd and Christoffersen (2014) applied a back force similar to the one evaluated by Walter et al. (2012) to Store's calving front and showed that this realistic forcing could have a significant impact on the advance and retreat cycles of the glacier front.
In this article we examine the consequence of submarine frontal melting and the ice mélange on glacier dynamics and on the behaviour of the glacier front using a full Stokes iceflow finite element model combined with calving parameterization based on damage and fracture mechanics. This enables a complete representation of the stress field in the vicinity of the front and provides a reliable tool to study front dynamics. To be sure our conclusions are robust for a number of glacier geometries and flow specifications, we ran more than 200 simulations combining a wide range of glacier sizes, flow and damage parameters, and forcing constraints. We provide a brief description of the model in Sect. 2, and in Sect. 3 we describe the setup and list the parameters. In Sect. 4 we describe glacier responses to seasonally variable forcings, and in Sect. 5 we provide a deeper analysis of the processes and mechanisms and compare the behaviour of grounded and floating glaciers.

Ice-flow model
We considered an incompressible, isothermal, and gravitydriven ice flow. The ice exhibits non-linear viscosity, and the flow is ruled by the Stokes equations, which reads where σ represents the Cauchy stress tensor, g the gravity force vector, ρ i the density of ice, and u the velocity vector. The Cauchy stress tensor can be expressed as a function of the deviatoric stress tensor S and the cryostatic pressure p with σ = S − pI and p = − tr(σ )/3. Ice rheology is represented by a non-linear Norton-Hoff-type flow law called Glen's flow law, which can be expressed as This equation links the deviatoric stress tensor S to the strain rate tensor˙ . The effective viscosity η is written as where I 2 2 represents the square of the second invariant of the strain rate tensor, A is the fluidity parameter, and E is an enhancement factor. A complete description of the model can be found in Gagliardini et al. (2013).

Damage and calving model
The ice-flow model described above was coupled with a calving model based on damage and fracture mechanics. Damage mechanics were used to describe the slow degradation of the mechanical properties of ice under the stress field, averaged at a mesoscale. Linear elastic fracture mechanics (LEFM) were used to describe the brittle initiation and propagation of crevasses.
Our damage model is inspired from the work of Pralong and Funk (2005) and relies on damage mechanics (Lemaitre et al., 1988). The level of isotropic damage in the ice is quantified by a scalar variable D called the damage variable, which equals 0 for undamaged ice and tends to 1 for fully damaged ice. In order to avoid singularity when D = 1, an upper bound is set such that D cannot exceed 0.7, accounting for the fact that D = 0.6 usually refers to fully damaged ice (see Pralong and Funk, 2005;Borstad et al., 2013). Damage is advected with the ice flow and it follows where B is a numerical parameter called damage enhancement factor. Damage increase depends on the stress state: Here σ I is the maximum principal stress and σ th represents a stress threshold for damage initiation. χ is called the damage criterion and quantifies the damage source term. Then, damage alters the deviatoric part of the Cauchy stress tensor S by introducing an effective deviatoric stress: This new effective stress reduces the effective viscosity of the ice through the expression of the enhancement factor that enters Eq. (4): The depth of crevasse fields can be represented using a damage contour and a stress history of the ice can be recorded.
This damage model was coupled with an LEFM model (inspired from van der Veen, 1998a, b), which was used to represent the rapid propagation of crevasses that characterizes calving events. In LEFM, the initiation of crevasse propagation depends on the stress intensity factor K I . To trigger propagation, the stress intensity factor, which depends on the size of the initial flaw and on the stress field, must be higher than the ice toughness K Ic . An initial crevasse depth provided by the previously computed damage field was used to compute the stress intensity factor. Once propagation is initiated, the stress intensity factor is computed at sea level. If the stress intensity factor is higher than an arrest criterion K Ia , the crevasse continues to propagate until it reaches the bottom of the glacier and triggers calving. Crevasse propagation may be facilitated by surface meltwater entering the crevasse. However, our model currently does not incorporate this process especially because of the lack of field observation that would be required to constrain it (see Krug et al., 2014, for details).
Among the numerical parameters required to run the model, three have to be calibrated, and are discussed below: the damage critical value D c , the stress threshold σ th , and the damage enhancement factor B. The damage contour given by D = D c corresponds to the depth of pre-existing flaws (from which LEFM is applied), σ th is the load that has to be applied to trigger ice damaging, and B quantifies the rate at which damage increases. The criterion for calving was initially proposed by Benn et al. (2007b) for the calculation of penetration depth of surface crevasses. It was then expanded by Nick et al. (2010) in order to incorporate the growth of basal crevasses. The main difference with our model is that the crevasse propagation, in the formulation of Benn et al. (2007b), does not rely on linear elastic fracture mechanics.
The model summarized here is described in detail in Krug et al. (2014) (along with all sensitivity tests) and implemented in the finite element open-source model Elmer/Ice (see Gagliardini et al., 2013).

Geometries and boundary conditions
We wanted to generalize our conclusions to a wide range of two-dimensional synthetic flow line glacier geometries of time-varying length (L) and thickness (H ). To this end, we built 60 geometries that depend on five parameters: the inlet ice flux (F inlet ), water depth (H w ), and damage parameters (σ th , B, D c ). These parameters were sampled using a Latin hypercube sampling (LHS) method, which ensures that each probability distribution in the model is evenly sampled, using a given number of simulations and a given number of parameters to sample. The glaciers were built up from ice slabs initially grounded on a linear analytical prograde slope (1 %). The choice of prograde slopes is motivated by the fact that dynamical instabilities arising from retrograde sloped glaciers make intercomparison of a large set of geometries difficult. The meshes comprise ∼ 7000 quadrilateral elements. Horizontally, the refinement is higher at the front (10 to 15 m) and in the vicinity of the upper surface (2 to 3 m) to account for the processes that occur at the calving front as well as the production of damage and advection at the upper surface. The setup is illustrated in Fig. 1. Boundary conditions are the same as those given in Krug et al. (2014). In addition, some specific conditions are given below.
-At the bed, the glacier can be either grounded or floating. The grounding line position is obtained through the resolution of a contact problem (Durand et al., 2009). The basal friction is linearly decreasing along the flow from 1.5 × 10 −2 MPa m −1/3 a 1/3 at the inlet boundary to 1.0 × 10 −4 MPa m −1/3 a 1/3 at x = 10 km. Feedbacks on basal friction arising from changes in glacier geometry are not studied here.
-As glacier thickness can vary with time, the total depthintegrated flux through the inlet boundary is kept constant (F inlet ).
-A lateral friction is prescribed to account for a constant fjord width of 10 km. This parameterization follows Gagliardini et al. (2010).
For the purpose of comparison, we chose to apply submarine melting and ice mélange forcing on glaciers in a quasi-steady-state (QSS) mode, i.e. their front has to stabilize within a given range lower than the length of one calving event. Among the 60 simulations generated by the LHS sampling, 20 had this feature and are listed in Table 1. Other geometries either advanced too far without calving or collapsed because of prolific calving. The sets of damage parameters with which a QSS was reached generally differed slightly from those calibrated in Krug et al. (2014). B ranged from 1.5 to 3 MPa −1 and σ th from 0.01 to 0.11 MPa (compared with 0.5 to 2 MPa −1 and 0.01 to 0.2 MPa respectively in Krug et al., 2014). The explanation for these differences is straightforward: the geometries studied here flow on a linear bedrock with no bumps or roughness. Consequently, except near the front, no high-velocity gradients appeared in the upper surface. Damage is consequently more difficult to initiate than in cases of rough bedrock and consequently has to be promoted. In addition, since thinner glaciers are subject to less internal stress, they require parameters that promote damaging, unlike thicker glaciers (see Table 1).
For the sake of clarity, out of the 20 representative geometries, unless otherwise specified, in Sects. 4 and 5, we only use one to illustrate the model's response. This terminusfloating geometry (hereafter referred to as Geo 9) is shown in bold in Table 1. However, the conclusions obtained in this study are robust against all the geometries considered, as discussed in Sects. 4.1 and 4.2.

Model experiments
The 20 setups summarized in Table 1 were run for 7 years to reach the QSS discussed above (spin-up). After this spinup, for each setup we imposed eight perturbations in melting or ice mélange as well as a control run (CR), in which the glacier continues its QSS evolution (i.e. without any melting or ice mélange forcing). These forcings were maintained for 5 years, after which they were removed, and we let the geometries evolve freely for 5 more years (relaxation period). In total, we performed 180 simulations of 17 years each.
The perturbations described in Sect. 3.2.1 and 3.2.2 are listed in Table 2. Perturbations in submarine frontal melting are named U1 to U4 (for "undercutting"). Ice mélange perturbations are named S1 to S4 (for "Sikussak", the Greenlandic word for ice mélange).

Submarine melting parameterization
Glacier frontal melting usually results from warm saline ocean water entering the fjord and mixing with fresh and cold subglacial freshwater flow. The resulting current melts the ice it meets as it rises along the calving face (Motyka H T refers to the QSS mean ice thickness of the terminus, F inlet represents the ice flux at the inlet boundary, and H w T is the QSS mean water depth at the terminus. H AB is the mean height above buoyancy at the front, where the terminus is grounded. The letter F is used instead of H AB when the glacier is afloat. σ th is the stress threshold that starts damage, B is the damage enhancement factor, and D c the damage contour. The line in bold is the representative simulation.  Table 1. Runs U1 to U4 refer to the undercutting experiments. The maximal melt rate (MMR) is indicated byṁ. Runs S1 to S4 refer to the ice mélange experiments: σ b max is the maximum ice mélange back pressure applied over a depth h and results in a maximum back force of max(σ b (t) · h). The control run did not undergo either ice mélange or melting. For each forcing, the "realistic" cases are in bold.
Run Following these studies and measurements, we tested different MMR summer values ranging from 0.41 to 12 m day −1 , following a 4-month sinusoidal peak and decay. In winter, following Cook et al. (2014) we prescribed a constant MMR of 0.41 m day −1 (see Fig. 2a). The melt rate was imposed in the front-normal direction, and its value is listed in Table 2. We deliberately chose to ignore melting at the bottom surface in the cases when the glacier started to float. We do not deny that this choice is a limitation, but we made it because we had no information on how the measured melt rate is distributed below the floating tongue and along the calving face. Considering that prescribing the same MMR under the glacier tongue would lead to its rapid collapse and that we wanted to compare the different behaviours of grounded and floating glaciers, taking into account melting under the tongue would require a more complex melting parameterization which is beyond the scope of this study.
In addition, some modelling work suggest that "the melting increases with height above the freshwater subglacial discharge", leading to an overcutting effect rather than the classical undercutting effect (Kimura et al., 2014). We do not consider this distribution, because the subsequent calving process would rely on basal crevasses, which are currently not incorporated into the model.

Ice mélange parameterization
Although ice mélange and its effect on glacier dynamics have been studied for a few decades (Rignot and MacAyeal, 1998;Reeh et al., 2001;Joughin et al., 2008c), two major unknowns remain: (i) the speed at which it becomes rigid and collapse and (ii) the force it applies against the glacier front.
i. Seale et al. (2011) studied the correlation between ice mélange disintegration and front retreat in fjords in Greenland using MODIS imagery and showed that disintegration can occur in a very short time, from a few days to a couple of weeks, whereas sea ice stiffening can take much longer. However, due to the lack of solar illumination, they did not obtain reliable information regarding ice mélange formation. We thus chose to simulate a growing period of 5 months (150 days) followed by a 20-day decay period.
ii. The question of the force transmitted up glacier by the ice mélange is more complex: using a two-dimensional flow line model, the back force must be represented as pressure applied over a thickness. Measurements of the back pressure σ b max are difficult to obtain. The main studies inferred and used a broad range of values between 0.02 to 3 MPa (Nick et al., 2009;Vieli and Nick, 2011;Walter et al., 2012;Cook et al., 2014;Todd and Christoffersen, 2014). To investigate glacier response to maximum back pressure, we tested values up to 1 MPa. Most of studies consider the mélange strength to be much lower than 1 MPa. Considering the fact that sea ice strength depends on many parameters (temperature, salinity) which are poorly constrained for mélange, we think that this value is a reliable upper bound for σ b max . The mélange thickness h was broadly estimated from 70 to 130 m in several studies (Fricker et al., 2005;Seneca Lindsey and Dupont, 2012;Cook et al., 2014;Todd and Christoffersen, 2014).
Considering that sea ice binds fragments of icebergs together, we can reasonably assume its stiffness is the firstorder control on mélange strength. Anderson (1961) linked the increase in sea ice thickness to the square root of time and to the gradient between oceanic and atmospheric temperatures. Thus, considering sea ice strength to be closely correlated with its thickness and keeping the same kind of kinetics, we expressed the back pressure applied by the mélange on the glacier as √ t(mod 365 days) at the end of winter 0(mod 365 days) in summer .
Ice mélange growth and decay are depicted in Fig. 2b. Finally, the ice mélange was prescribed through a timevarying back-stress assumed to be homogeneous over its thickness and resulting in a total back force equal to the product σ b (t) · h. We tested several combinations of mélange thickness and back pressure parameters, but for the rest of this study we only illustrate the most representative (S1 to S4, see Table 2).

Melting impact
The simulation starts on 1 January (time = 0), when the prescribed melt rate is set to its minimum value. The distribution of maximum melt rate for U1 to U4 and for the control run is given in Fig. 3a. The control run was not subject to any melt rate and its front position never moved by more than a few tens of metres (see Fig. 3b). Figure   summer season (day 173). Simulations U1 to U4 produced slight oscillations, resulting in a slight advance of the terminus compared to the control run, but it never moved more than 400 m downstream. These advances may seem counterintuitive as most research suggests that submarine melting causes the front to retreat. However, our model revealed that when an advance takes place, it is not triggered by the same mechanism in all the setups. It is related to (i) a decrease in the frequency of calving events (this process is described in Sect. 5.2 below) and/or (ii) a torque effect caused by the retreat of the lowest point of the front (due to melting) and the advance of the highest point. In the case presented here, the advance of the front is due to a decrease in the frequency of calving events. Its geometry is illustrated in Fig. 4. As soon as the forcing was removed, all the fronts reached their QSS position within a few months, except in simulation U4 (whose specific behaviour is discussed in Sect. 5.2). The ice velocity at the front varied within a range of 200 m yr −1 , which is very similar to the natural variation in the control run after a calving event (see Fig. 3c).
Considering the contribution of melting and calving to ice loss, the volume of calved ice always appears to be larger than the melted volume. During summer, melting accounted for up to ∼ 23 % of the total mass loss (Fig. 5). Comparing winter and summer suggests that an increase in the intensity of undercutting does not significantly alter the total loss: more ice is melted but less ice is calved, meaning the cumulated volume does not vary significantly over the seasons.
To account for the different geometries, we summarized them as a function of their QSS mean thickness and velocity at the terminus for a given realistic forcing in melting (U2) (Fig. 6). The conclusions drawn above were qualitatively confirmed for all the other geometries: mean ice loss during the melting season was comparable to the ice loss during the rest of the year, whatever the size and velocity of the glacier. In summer, the calved volume was reduced by the increasing melt rate, but the cumulative loss remained unchanged.

Ice mélange impact
To measure the impact of ice mélange, we ran simulations S1 to S4 (Table 2), as well as the control run. Figure 7a shows the ice mélange intensity. Figure 7b shows changes in the position of the front as a function of time in the four corresponding experiments and in the control run (dashed black line). Two types of behaviour were observed during the first 5 years. In winter, the ice mélange strengthened, and calving frequency decreased or stopped. As a consequence, the front advanced. In summer, the decay of the ice mélange back force was immediately followed by a rapid sequence of largescale iceberg calving events. In all the perturbation simulations, the glacier front was always located further downstream than in the control run. Each winter, the gaps between the positions of the S1-S4 fronts and that of the control run increased to reach a value of 500 m in S1 and 3 km in S4. Moreover, in perturbations S2, S3, and S4 after each year, the front did not retreat back to its QSS position, suggesting a consequence for interannual mass loss. These behaviours are consistent with observations, confirming the hypothesis that a strong mélange reduces calving discharge (Sohn et al., 1998;Joughin et al., 2008c). Figure 7c shows the ice velocity at the front using the same colour scale. The position of the terminus was inversely correlated with the velocity of the ice. The advance of the front in winter led to a decrease in ice velocity due to the increasing buttressing effect of the glacier sliding against the fjord walls and, to a lesser extent, to increasing back pressure with increasing ice mélange strength (e.g. S3 and S4). As can be seen for the three highest back forces (S2, S3, S4), when the mélange collapsed, the ice flow at the front accelerated to a faster speed than the maximum in the control run. The increase in speed can be explained by the following chronology: first, the release of the mélange back force accelerated the ice flow; second, after the first calving event was triggered, the resulting geometry was a high vertical ice cliff. Velocity vectors were no longer parallel to sea level and a torque appeared, leading to a force imbalance that further increased the ice velocity at the front. The red inset in Fig. 7 shows the first year of the forcing and underlines this phenomenon: at stage 1, mélange strength was maximum and its decay accelerated ice flow. At stage 2, the major front retreat in simulation S4 occurred, further accelerating the flow.
Such rapid acceleration of the flow was observed during calving events at the front of the Jakobshavn glacier in May 2007 by Amundson et al. (2008). Near the front, GPS stations recorded an increase in ice velocity from 11 315 m a −1 to 12 775 m a −1 in 20 days, during which the glacier underwent three calving events that were attributed to an adjustment in the stress field. The freshly calved glacier was shorter, the buttressing effect was reduced, and glacier flow accelerated (Benn et al., 2007b). Walter et al. (2012) monitored the magnitude of the speedup at 550 m a −1 during the days following the break up of the ice mélange at the terminus of Store Glacier. Our results are in agreement with these variations in measured velocity. Figure 8 shows the mass loss for all geometries for perturbation S2. Whatever the glacier geometry, the loss in summer was greater than in winter, as the winter ice mélange layer reduced calving activity. However, the same back force does not have the same effect on small and large glaciers. In the case of smaller glaciers, it completely prevents calving; in the case of larger geometries, it only decreases the iceberg discharge. This explains why in the winter inset in Fig. 8, the thickest glaciers show the smallest contrast between winter and summer. This feature is also visible in Fig. 7b: large ice mélange intensities prevent calving (light blue, yellow, and red curves), while smaller intensities simply reduce calving frequency (dark blue). This consideration reinforces the need for better knowledge of the properties of ice mélange.

Mechanical impact of the ice mélange on the glacier front
According to Amundson et al. (2010), to prevent the rotation of a calved iceberg away from the terminus, the required back force is between ∼ 1.0 × 10 7 N m −1 and 10.0 × 10 7 N m −1 , depending on the glacier flotation and the inclination angle of the iceberg. Concerning perturbations caused by the ice mélange, our model suggests that calving ceases as soon as the applied back force reaches a given value. We investigated model sensitivity to the applied back force by evaluating the value of  Figure 6. Daily averaged ice loss over the winter and summer seasons for the 5 years of the simulation for the setups listed in Table 1 (experiment U2). The area of the disks represent the volume of ice lost by the glacier associated with a mean ice front velocity V T and ice front thickness H T . The fraction of ice loss due to melting is in green (also indicated by the percentage) and the fraction due to calving is in blue. The melted fraction accounts for 5 to 39 % of the summer glacier loss (highest for smaller glacier) and agrees with the lower bounds of the calculations of Rignot et al. (2010), especially for thinnest glaciers. this threshold. To this end, we isolated the pairs of parameters (σ b max , h) for which the winter season was characterized by the absence of calving events. For each of the five winter seasons, we then calculated the back pressure applied when calving ceased using Eq. (5). Multiplying this back-stress by the thickness of the ice mélange gives a back force per metre of lateral width. The corresponding distribution for the 45 values is given in Fig. 9, of which the mean value is around 1.1 × 10 7 N m −1 . Thus, the value of the back force that prevents the iceberg rotation (as calculated by Amund-  Table 1 (experiment M2). The area of the disks represents the volume of ice lost by the glacier associated with a mean ice front velocity V T and ice front thickness H T . Our coupled ice flow and calving model enabled us to distinguish between the different processes that culminate in iceberg calving that could be affected by the ice mélange. In the first stage, development of the crevasse field is determined by the damage criterion χ, which quantifies the damage increase in the ice. Figure 10a shows changes in the position of the terminus over a period of 150 days, i.e. a full winter season. Three key events are highlighted by diamond symbols. The red diamond corresponds to a situation in which the glacier is about to calve, the yellow diamond illustrates a case where the glacier is subject to ice mélange, and the blue diamond corresponds to a situation in which the glacier has just calved. Figure 10b shows the value of χ along the upper surface where the tensile stress is the highest. During the ice mélange season (yellow curve), damage production is slightly lower than that in the pre-and post-calving situations but remains positive. This means that ice mélange reduces the production of damage but that the effect is too weak to completely halt damage to the ice.
Following damage to the ice, three criteria have to be fulfilled to trigger calving: the condition on damage contour D = D c , the initiation of fracture propagation at a depth given by the D c contour, and propagation to sea level. These criteria are shown in Fig. 10c, where the stress intensity factor is represented on the vertical axis. The horizontal extent of the solid coloured lines shows the length of the damage envelope (the front is located on the right side of the figure): a longer one illustrates a more extended crevasse field. As expected, the solid blue curve is almost nonexistent: as calving has just occurred, the crevasse field is not deep enough to apply LEFM model (D < D c ) or the stress intensity factor is lower than the propagation threshold (K I < K Ic ). In contrast, the solid red and yellow lines show that the surface is sufficiently damaged to reach the criterion D = D c over a larger surface area. Wherever the stress intensity factor becomes higher than the ice toughness (K I > K Ic ), crevasse propagation begins. This condition was satisfied in the case of the red and yellow lines near to the upper surface, so propagation can begin. The criterion K I > K Ia must then be validated at sea level. The stress intensity at sea level is represented by the coloured dashed curves. The red curve satisfies this criterion at some points, meaning that the calving event can be- gin. However, when the ice mélange layer is present (yellow dashed line), this criterion is not fulfilled, and as a consequence the crevasse cannot propagate down to sea level. Regarding model sensitivity to the thickness of the ice mélange, we observed that a thinner layer associated with a stronger back force reduced the calving rate more than a thicker layer associated with a weaker back force, with the same total back force (data not shown). This is because the thinnest layer of ice mélange is concentrated at sea level, which significantly reduces the stress intensity factor at this depth, thereby preventing crevasse propagation.
These mechanisms could explain the behaviour observed in Fig. 7b. For simulations S2, S3, and S4, the figure shows that the decay of the mélange layer was followed by a "cascade" of calving events. However, the glacier did not immediately retreat to its QSS position. This rate of retreat depends on the degree of damage at the surface, which depends on the driving force of the ice mélange. This suggests that a stronger or longer winter season could alter the position of the front over a period of more than 1 year, when one considers the "stress history" of the glacier, and does not only rely on a one-off record of the stress balance.
The results presented in this section are in direct contrast with those of Cook et al. (2014). These authors observed no remarkable changes in the position of the front unless they applied a back force of 50.0 × 10 7 N m −1 , although they simulated a difference of 25 m in the longitudinal extent of the crevasse field near the front for smaller values of ice mélange (5.0 × 10 7 N m −1 ). Conversely, the results of Todd and Christoffersen (2014) clearly agree with ours, as they simulated a comparable advance of the front (∼ 1.5 km) with a back force close to ours in simulation S2. Finally, for an applied force of the same order of magnitude, our model shows that the ice mélange acts sooner than suggested by Amundson et al. (2010) by preventing the propagation of the fracture down to the glacier base. Figure 11a shows the maximum difference in the surface along-flow component of the deviatoric stress tensor S xx between the U2 and the CR simulations, in the vicinity of the front (< 600 m) during the middle of the first summer period. Undercutting grounded glaciers slightly increased the tensile stress at the upper surface (red dots). It increased the frequency of calving events but reduced the duration of each event (see red dots and crosses in Fig. 11b). Conversely, in the case of floating glaciers, the surface adjustment of the tongue decreased the tensile stress compared with the control run (Fig. 11a, blue dots). Consequently, the frequency of calving events decreased slightly, but the distance the front retreated at each event increased slightly (blue dots and crosses in Fig. 11b).

Differences between floating and grounded termini
Concerning the behaviour of the grounded geometry of Helheim Glacier, Cook et al. (2014) Table 1 between the U2 experiment and the control run (CR). The difference is computed within the last 600 m before the calving front in the middle of the first summer (day 182). Red dots indicate grounded fronts, while blue dots highlight floating termini. (b) Ratio between the summer and the winter frequency of calving event (dots) and ratio between the summer and the winter length of calving front retreat (crosses) for all the geometries listed in Table 1 forced by the U2 experiment (maximum melt rate: 6 m day −1 ), as a function of the inlet flux. Grounded glaciers are in red, floating glaciers in blue. Circles and crosses refer to the mean number of calving events and to the mean length of the front retreat ratios respectively. front, unless the prescribed melt rates are extremely high (up to 20 m day −1 ). On Store Glacier, Todd and Christoffersen (2014) modelled a slight increase in frequency with undercutting, as well as a decrease in the amplitude of the retreat of the calving front, and they attributed this interannual stability to the glacier's topographic setting. As the geometry of Todd and Christoffersen (2014) was grounded for most of the melt season, their modelling results are in agreement with ours. Here, it is worth to be mentioned that we only managed to observe a front retreat with especially high melt rates (≥ 12 m day −1 ). Incorporating melting from underneath the floating ice tongue would probably make the glacier front to retreat. In our simulations, another difference appeared between grounded and floating glaciers. In Sect. 4, we showed that for most perturbations, after the relaxation period, glacier fronts usually reached their QSS position. However, this was not the case for the glacier undergoing perturbation U4 illustrated in Fig. 3b. Indeed, its front rapidly advanced further downstream than the others and appeared to stabilize at an extent of 6.2 km, compared with 6.0 km for the other fronts. When extended to other geometries, the same result also was obtained in some of the ice mélange experiments. However, it only concerned glaciers with a floating tongue and only occurred under the strongest forcings.
Concerning these processes, we propose an explanation for this phenomenon (see Fig. 12). The melting perturba-tion applied on the glacier front affects the shape of the floating tongue (Fig. 12b). It reduces its area along with the subsequent buttressing effect. As a consequence, the whole glacier accelerates and thins, and the grounding line retreats (Fig. 12c). The ice flux at the grounding line is therefore modified, and a new equilibrium is established that relies on interactions between the ice flow, damage production, and the calving law. Considering the ice mélange, the concept is similar but the process is reversed (Fig. 12d). As the ice mélange prevents the floating tongue from calving, the area of the tongue increases. Consequently, the glacier slows down and thickens, and the grounding line advances (Fig. 12e). Again, a new equilibrium may be established with an associated quasi-steady state front dynamics.

Conclusions
Ice mélange and melting of the glacier front have been reported by many authors to influence the behaviour of tidewater glaciers. In particular, they have been cited as a possible explanation for the seasonal advance and retreat cycles of glacier fronts, among other external forcings. However, although some correlations between these mechanisms and the advance/retreat of the front have been established on many outlet glaciers in Greenland, little is known about the exact role of these forcings.  In this study, we combined a full Stokes ice flow model with a calving framework using damage and fracture mechanics to investigate the impact of these forcings on glacier dynamics. This allowed us to represent the slow degradation of the mechanical properties of the ice and the initiation and propagation of pre-existing fractures, which are essential to describe the processes occurring at the front. We performed experiments on a large set of synthetic geometries using different values for melting and ice mélange back-stress and thickness, and the conclusions we have drawn are robust in all these experiments. However, it is important to note that the model used here considers surface crevasses only. A deeper analysis would require the modelling of the development and propagation of basal crevasses, and their feedback with the stress field and the glacier dynamics.
Our modelling showed that frontal melting has an impact on the calving rate and on the position of the front (less than a few hundred metres), but no effect on inter-/multiannual mass loss. On the contrary, applying an ice mélange layer against the front affects its position to a larger extent (up to several kilometres) compared to melting. In addition, its consequences for the inter-/multiannual mass loss, when slight, may not be completely negligible and thus support Joughin et al. (2008c), according to whom "It is likely that the processes that control the seasonal calving cycle may also influence the interannual variability". By investigating the processes occurring during calving events, we have shown firstly that the ice mélange reduces the rate of surface damage by reducing the tensile stress in the glacier upper surface and secondly prevents fracture propagation at sea level and hence calving. Better field characterization of undercutting and ice mélange properties should increase the accuracy of further modelling.
Finally, our results also reveal a feature that is specific to glaciers with floating termini, i.e. that strong perturbations (either in melting or in ice mélange) may affect their multi-annual behaviour. By affecting the buttressing effect of the tongue, the perturbation may modify the subsequent glacier equilibrium and lead to a new stable geometry for the same model parameters. This new stable position then depends on feedback between glacier flow and calving law parameters.