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

Abstract

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 finite element-based model that solves the full Stokes equations to simulate the impact of these forcings on two-dimensional theoretical flow line glacier configurations. 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 floating glaciers: in the case of a floating extension, the strongest forcings can disrupt the glacier equilibrium by modifying its buttressing and ice flux at the grounding line.

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 , 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 $\pm$ 49 Gt a${}^{-\mathrm{1}}$ on average over the past 2 decades, has increased in recent years to reach an estimated value of 263 $\pm$ 30 Gt a${}^{-\mathrm{1}}$ between 2005 and 2010 and 359.8 $\pm$ 28.9 Gt a${}^{-\mathrm{1}}$ between April 2009 and 2012 . This mass loss extended over a large part of the GIS , which is thus becoming a major contributor to sea level rise .

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 . 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 , corresponding to $-$156.3 $\pm$ 40.9 Gt a${}^{-\mathrm{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.

hypothesized that the increased discharge in Greenland may have been triggered by an increase in the subsurface ocean temperature. This claim was supported by , 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 . 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${}^{-\mathrm{1}}$ depending on the glacier and the associated fjord system .

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 . 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 and that summer decay is followed by an increase in the calving rate, a retreat of the front, and accelerated ice flow . Some authors argue that ice mélange may directly resist the ice flow , while others suggest that it only maintains the integrity of the terminal part of the glacier . 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 , 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, 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 applied on two-dimensional flow line geometries of Helheim and Store glaciers tempered these conclusions. According to , 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 and using a depth- and width-integrated model combined with the calving parameterization of . 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. 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 applied a back force similar to the one evaluated by 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 ice-flow 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. , and in Sect.  we describe the setup and list the parameters. In Sect.  we describe glacier responses to seasonally variable forcings, and in Sect.  we provide a deeper analysis of the processes and mechanisms and compare the behaviour of grounded and floating glaciers.

Model presentation

Ice-flow model

We considered an incompressible, isothermal, and gravity-driven ice flow. The ice exhibits non-linear viscosity, and the flow is ruled by the Stokes equations, which reads $$\begin{array}{}& & \text{div}\left(\mathit{\sigma }\right)+{\mathit{\rho }}_{\text{i}}\mathit{g}=\mathrm{0}& & \text{div}\left(\mathit{u}\right)=\mathrm{0},\end{array}$$ where $\mathit{\sigma }$ represents the Cauchy stress tensor, $\mathit{g}$ the gravity force vector, ${\mathit{\rho }}_{\text{i}}$ the density of ice, and $\mathit{u}$ the velocity vector. The Cauchy stress tensor can be expressed as a function of the deviatoric stress tensor $\mathbf{S}$ and the cryostatic pressure $p$ with $\mathit{\sigma }$ $=$ $\mathbf{S}$ $-$ $p\mathbf{I}$ and $p$ $=$ $-$ $tr\left(\mathit{\sigma }\right)/\mathrm{3}$. Ice rheology is represented by a non-linear Norton–Hoff-type flow law called Glen's flow law, which can be expressed as $$\mathbf{S}=\mathrm{2}\mathit{\eta }\dot{\mathit{ϵ}}.$$ This equation links the deviatoric stress tensor $\mathbf{S}$ to the strain rate tensor $\dot{\mathit{ϵ}}$. The effective viscosity $\mathit{\eta }$ is written as $$\mathit{\eta }=\frac{\mathrm{1}}{\mathrm{2}}(EA{)}^{-\mathrm{1}/n}{\mathbf{I}}_{{\dot{\mathit{ϵ}}}_{\mathrm{2}}}^{(\mathrm{1}-n)/n},$$ where ${\mathbf{I}}_{{\dot{\mathit{ϵ}}}_{\mathrm{2}}}^{\mathrm{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 .

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 and relies on damage mechanics . 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. Damage is advected with the ice flow and it follows $$\frac{\partial D}{\partial t}+\mathit{u}\mathrm{\nabla }D=\left\{\begin{array}{ll}B\cdot \mathit{\chi }& \text{if}B\cdot \mathit{\chi }>\mathrm{0}\\ \mathrm{0}& \text{otherwise},\end{array}\right.$$ where $B$ is a numerical parameter called damage enhancement factor. Damage increase depends on the stress state: $$\mathit{\chi }\left({\mathit{\sigma }}_{\text{I}},,,{\mathit{\sigma }}_{\text{th}},,,D\right)=\text{max}\left\{\mathrm{0},,,\frac{{\mathit{\sigma }}_{\text{I}}}{(\mathrm{1}-D)},-,{\mathit{\sigma }}_{\text{th}}\right\}.$$ Here ${\mathit{\sigma }}_{\text{I}}$ is the maximum principal stress and ${\mathit{\sigma }}_{\text{th}}$ represents a stress threshold for damage initiation. $\mathit{\chi }$ is called the damage criterion and quantifies the damage source term. Then, damage alters the deviatoric part of the Cauchy stress tensor $\mathbf{S}$ by introducing an effective deviatoric stress: $$\stackrel{\mathrm{̃}}{\mathbf{S}}=\frac{\mathbf{S}}{(\mathrm{1}-D)}.$$ This new effective stress reduces the effective viscosity of the ice through the expression of the enhancement factor that enters Eq. (): $$E=\frac{\mathrm{1}}{(\mathrm{1}-D{)}^n}.$$ 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, 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_{\text{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_{\text{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_{\text{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 seefor details.

Among the numerical parameters required to run the model, three have to be calibrated, and are discussed below: the damage critical value $D_{\text{c}}$, the stress threshold ${\mathit{\sigma }}_{\text{th}}$, and the damage enhancement factor $B$. The damage contour given by $D$ $=$ $D_{\text{c}}$ corresponds to the depth of pre-existing flaws (from which LEFM is applied), ${\mathit{\sigma }}_{\text{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 for the calculation of penetration depth of surface crevasses. It was then expanded by in order to incorporate the growth of basal crevasses. The main difference with our model is that the crevasse propagation, in the formulation of , does not rely on linear elastic fracture mechanics.

The model summarized here is described in detail in (along with all sensitivity tests) and implemented in the finite element open-source model Elmer/Ice see.

Setup of the experiment. $H_{\text{T}}\left(x\right)$, $H_{w_{\text{T}}}\left(x\right)$, and $h$ represent glacier thickness, water depth, and ice mélange thickness respectively. The glacier is grounded on a solid bedrock with a slightly positive slope (exaggerated here) represented by the thick brown line. The grounding line is indicated by the red dot, at the abscissa $x$ $=$ $x_{\text{G}}$. The blue arrow shows the direction of the ice flow.

[Figure omitted. See PDF]

Setup and forcing parameterization

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_{\text{inlet}}$), water depth ($H_{\text{w}}$), and damage parameters (${\mathit{\sigma }}_{\text{th}}$, $B$, $D_{\text{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 $\sim$ 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. .

Boundary conditions are the same as those given in . 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 . The basal friction is linearly decreasing along the flow from 1.5 $\times$ 10${}^{-\mathrm{2}}$ MPa m${}^{-\mathrm{1}/\mathrm{3}}$ a${}^{\mathrm{1}/\mathrm{3}}$ at the inlet boundary to 1.0 $\times$ 10${}^{-\mathrm{4}}$ MPa m${}^{-\mathrm{1}/\mathrm{3}}$ a${}^{\mathrm{1}/\mathrm{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 depth-integrated flux through the inlet boundary is kept constant ($F_{\text{inlet}}$).

• A lateral friction is prescribed to account for a constant fjord width of 10 km. This parameterization follows .

List of geometries and their associated parameters used for the model experiments. ${\stackrel{\mathrm{‾}}{H}}_{\text{T}}$ refers to the QSS mean ice thickness of the terminus, $F_{\text{inlet}}$ represents the ice flux at the inlet boundary, and ${\stackrel{\mathrm{‾}}{H}}_{w_{\text{T}}}$ is the QSS mean water depth at the terminus. ${\stackrel{\mathrm{‾}}{H}}_{\text{AB}}$ is the mean height above buoyancy at the front, where the terminus is grounded. The letter $F$ is used instead of ${\stackrel{\mathrm{‾}}{H}}_{\text{AB}}$ when the glacier is afloat. ${\mathit{\sigma }}_{\text{th}}$ is the stress threshold that starts damage, $B$ is the damage enhancement factor, and $D_{\text{c}}$ the damage contour. The line in bold is the representative simulation.

Geometry	${\stackrel{\mathrm{‾}}{H}}_{\text{T}}$	$F_{\text{inlet}}$	${\stackrel{\mathrm{‾}}{H}}_{w_{\text{T}}}$	$H_{\text{AB}}$	${\mathit{\sigma }}_{\text{th}}$	$B$	$D_{\text{c}}$
	(m)	($\times$ 10${}^{\mathrm{3}}$ m${}^{\mathrm{2}}$ a${}^{-\mathrm{1}}$)		(m)	(MPa)	(MPa${}^{-\mathrm{1}}$)
Geo 1	358	710	308	16	0.017	2.9	0.47
Geo 2	356	679	307	14	0.014	1.6	0.42
Geo 3	362	488	319	7	0.025	1.8	0.44
Geo 4	362	572	317	10	0.041	2.2	0.46
Geo 5	456	1210	405	6	0.026	2.9	0.41
Geo 6	474	1133	427	F	0.021	1.5	0.44
Geo 7	465	1528	412	8	0.013	1.8	0.43
Geo 8	461	1432	409	7	0.037	2.2	0.41
Geo 9	631	3940	623	F	0.059	2.4	0.54
Geo 10	638	4406	632	F	0.024	2.1	0.6
Geo 11	597	2273	609	F	0.068	2.1	0.48
Geo 12	627	3535	632	F	0.021	1.7	0.53
Geo 13	824	7719	908	F	0.047	2.3	0.5
Geo 14	810	7143	909	F	0.050	1.6	0.41
Geo 15	860	10 203	975	F	0.032	1.5	0.48
Geo 16	842	8953	970	F	0.079	2.0	0.56
Geo 17	863	11078	923	F	0.068	2.0	0.41
Geo 18	866	11 389	942	F	0.016	2.2	0.59
Geo 19	854	9979	942	F	0.115	2.4	0.43
Geo 20	810	7233	924	F	0.100	1.9	0.46

Complete list of the experiments performed for each setup listed in Table . Runs U1 to U4 refer to the undercutting experiments. The maximal melt rate (MMR) is indicated by $\dot{m}$. Runs S1 to S4 refer to the ice mélange experiments: ${\mathit{\sigma }}_{b_{\text{max}}}$ is the maximum ice mélange back pressure applied over a depth $h$ and results in a maximum back force of max(${\mathit{\sigma }}_{\text{b}}\left(t\right)$ $\cdot$ $h$). The control run did not undergo either ice mélange or melting. For each forcing, the “realistic” cases are in bold.

Run	${\mathit{\sigma }}_{b_{\text{max}}}$	$h$	max(${\mathit{\sigma }}_{\text{b}}\left(t\right)$ $\cdot$ $h$)	MMR
name	(kPa)	(m)	($\times$ 10${}^{\mathrm{7}}$ N m${}^{-\mathrm{1}}$)	(m day${}^{-\mathrm{1}}$)
U1	0	0	0.0	3
U2	0	0	0.0	6
U3	0	0	0.0	9
U4	0	0	0.0	12
S1	170	80	1.36	0
S2	200	100	2.0	0
S3	350	120	4.2	0
S4	750	80	6.0	0
CR	0	0	0.0	0

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 . 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 . $B$ ranged from 1.5 to 3 MPa${}^{-\mathrm{1}}$ and ${\mathit{\sigma }}_{\text{th}}$ from 0.01 to 0.11 MPa compared with 0.5 to 2 MPa${}^{-\mathrm{1}}$ and 0.01 to 0.2 MPa respectively in. 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 ).

For the sake of clarity, out of the 20 representative geometries, unless otherwise specified, in Sects.  and , we only use one to illustrate the model's response. This terminus-floating geometry (hereafter referred to as Geo 9) is shown in bold in Table . However, the conclusions obtained in this study are robust against all the geometries considered, as discussed in Sects.  and .

Model experiments

The 20 setups summarized in Table  were run for 7 years to reach the QSS discussed above (spin-up). After this spin-up, 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.  and  are listed in Table . 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 . The melting intensity appears to be tightly linked with ocean water circulation, water stratification and its variability, as well as the specificities of the fjord, topography, size, or runoff seasonality and could partly explain the wide range of different measurements from one glacier to another. measured summer melt rates ranging from 0.6 to 3.8 m day${}^{-\mathrm{1}}$ at the face of four calving glaciers in west Greenland. calculated an annual mean melt rate of around 2 m day${}^{-\mathrm{1}}$ in the area of the Sermilik Fjord–Helheim Glacier, whereas measured summer values of around 10 m day${}^{-\mathrm{1}}$ along the face of Kangerdlugssuaq Glacier, southeast Greenland. Using a similar technique, in Alaska, obtained a range between 9 and 17 m day${}^{-\mathrm{1}}$ for Yahtse Glacier, and measured 12 m day${}^{-\mathrm{1}}$ at the calving front of LeConte Glacier.

Most parameterizations of frontal melt in ice flow models published so far assume a linear variation of melt from 0 at sea level to a maximum value at the lowest point of the front. Following these parameterizations, the maximum value is used to characterize the intensity of the melt and is referred to as maximal melt rate (MMR). applied a MMR of 8 m day${}^{-\mathrm{1}}$ during 3 summer months and 0 m day${}^{-\mathrm{1}}$ for the rest of the year. tested different values of MMR ranging from 2.7 to 13 m day${}^{-\mathrm{1}}$ during the 5-month summer period and a winter constant MMR of 0.41 m day${}^{-\mathrm{1}}$ (actually 150 m yr${}^{-\mathrm{1}}$).

Following these studies and measurements, we tested different MMR summer values ranging from 0.41 to 12 m day${}^{-\mathrm{1}}$, following a 4-month sinusoidal peak and decay. In winter, following we prescribed a constant MMR of 0.41 m day${}^{-\mathrm{1}}$ (see Fig. a). The melt rate was imposed in the front-normal direction, and its value is listed in Table . 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 . 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 , 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.

  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 ${\mathit{\sigma }}_{b_{\text{max}}}$ are difficult to obtain. The main studies inferred and used a broad range of values between 0.02 to 3 MPa . 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 ${\mathit{\sigma }}_{b_{\text{max}}}$. The mélange thickness $h$ was broadly estimated from 70 to 130 m in several studies .

Considering that sea ice binds fragments of icebergs together, we can reasonably assume its stiffness is the first-order control on mélange strength. 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 $${\mathit{\sigma }}_{\text{b}}\left(t\right)=\left\{\begin{array}{ll}\frac{{\mathit{\sigma }}_{b_{\text{max}}}}{\sqrt{150}}\sqrt{t}\left(\text{mod}365\text{days}\right)& \text{in}\text{winter}\\ {\mathit{\sigma }}_{b_{\text{max}}}-\frac{{\mathit{\sigma }}_{b_{\text{max}}}}{\sqrt{20}}\sqrt{t}\left(\text{mod}365\text{days}\right)& \text{at}\text{the}\text{end}\text{of}\text{winter}\\ \mathrm{0}\left(\text{mod}365\text{days}\right)& \text{in}\text{summer}\end{array}\right..$$ Ice mélange growth and decay are depicted in Fig. b.

Shape of perturbations over a period of 1 year. (a) Melting parameterization at glacier bottom. (b) Ice mélange parameterization.

[Figure omitted. See PDF]

Finally, the ice mélange was prescribed through a time-varying back-stress assumed to be homogeneous over its thickness and resulting in a total back force equal to the product ${\mathit{\sigma }}_{\text{b}}\left(t\right)$ $\cdot$ $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 ).

Glacier response in the undercutting experiments (Geo 9, Table ) U1 to U4 (coloured lines) and CR (dashed black line). (a) Variation in the maximum melt rate imposed at the base of the calving front. During the 5-month summer season, melting follows a sinusoidal pattern. Otherwise, a constant melt rate of 0.41 m day${}^{-\mathrm{1}}$ is prescribed. (b) Variation in the position of the front as a function of time and (c) ice velocity at the terminus. For the sake of clarity, the velocity was smoothed with a 10-day moving average.

[Figure omitted. See PDF]

Illustrative example of the glacier shape. (a) Velocity field for the representative geometry (Geo 9, Table ), undergoing perturbation U3 at day 173 (first summer season). The red dot shows the position of the grounding line. The glacier is about 6000 m long and 600 m thick. (b) Zoom in the black rectangle, glacier front geometry and mesh.

[Figure omitted. See PDF]

Results

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. a. 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. b). Figure  shows the shape of the glacier under perturbation U3 in the middle of the 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.  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. .

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. ). The ice velocity at the front varied within a range of 200 m yr${}^{-\mathrm{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 $\sim$ 23 % of the total mass loss (Fig. ). 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.

Mean daily winter and summer ice loss due to calving (blue) and melting (green) for five values of melt rate, over 5 years. Computed from the control run (CR) and Geo 9, perturbations U1 to U4.

[Figure omitted. See PDF]

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. ). 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.

Daily averaged ice loss over the winter and summer seasons for the 5 years of the simulation for the setups listed in Table  (experiment U2). The area of the disks represent the volume of ice lost by the glacier associated with a mean ice front velocity ${\stackrel{\mathrm{‾}}{V}}_{\text{T}}$ and ice front thickness ${\stackrel{\mathrm{‾}}{H}}_{\text{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 , especially for thinnest glaciers.

[Figure omitted. See PDF]

Ice mélange impact

To measure the impact of ice mélange, we ran simulations S1 to S4 (Table ), as well as the control run. Figure a shows the ice mélange intensity. Figure b 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 large-scale 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 .

Glacier response in ice mélange experiments (Geo 9, Table ) S1 to S4 (coloured lines) and in the control run (CR) (dashed black line). (a) Variation in mélange back force ${\mathit{\sigma }}_{\text{b}}\left(t\right)$ $\cdot$ $h$ as a function of time (per metre of lateral width). The value of the maximal back force depends on the pairs of parameters (${\mathit{\sigma }}_{b_{\text{max}}}$, $h$). (b) Variation in the position of the front as a function of time and (c) ice velocity at the terminus. For the sake of clarity, velocity is smoothed with a 10-day moving average. The red inset shows the ice velocity at the terminus without smoothing and the precise chronology of variations in velocity: the dashed vertical line 1 corresponds to the maximum mélange strength and line 2 corresponds to the major calving event in experiment S4.

[Figure omitted. See PDF]

Figure c 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.  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 . Near the front, GPS stations recorded an increase in ice velocity from 11 315 m a${}^{-\mathrm{1}}$ to 12 775 m a${}^{-\mathrm{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 . monitored the magnitude of the speedup at 550 m a${}^{-\mathrm{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  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. , the thickest glaciers show the smallest contrast between winter and summer. This feature is also visible in Fig. b: 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.

Daily averaged ice loss over the winter and summer seasons for the 5 years of the simulation for the setups listed in Table  (experiment M2). The area of the disks represents the volume of ice lost by the glacier associated with a mean ice front velocity ${\stackrel{\mathrm{‾}}{V}}_{\text{T}}$ and ice front thickness ${\stackrel{\mathrm{‾}}{H}}_{\text{T}}$.

[Figure omitted. See PDF]

Discussion

Mechanical impact of the ice mélange on the glacier front

According to , to prevent the rotation of a calved iceberg away from the terminus, the required back force is between $\sim$ 1.0 $\times$ 10${}^{\mathrm{7}}$ N m${}^{-\mathrm{1}}$ and 10.0 $\times$ 10${}^{\mathrm{7}}$ N m${}^{-\mathrm{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 this threshold. To this end, we isolated the pairs of parameters (${\mathit{\sigma }}_{b_{\text{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. (). 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. , of which the mean value is around 1.1 $\times$ 10${}^{\mathrm{7}}$ N m${}^{-\mathrm{1}}$. Thus, the value of the back force that prevents the iceberg rotation as calculated by and the one which prevents fracture propagation are on the same order of magnitude.

Histogram of the back force from an ice mélange required to prevent calving. The mean of the distribution is 1.1 $\times$ 10${}^{\mathrm{7}}$ N m${}^{-\mathrm{1}}$ and its standard deviation is 1.3 $\times$ 10${}^{\mathrm{6}}$ N m${}^{-\mathrm{1}}$ per metre of lateral width.

[Figure omitted. See PDF]

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 $\mathit{\chi }$, which quantifies the damage increase in the ice. Figure a 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 b shows the value of $\mathit{\chi }$ 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.

(a) Changes in the position of the front over a winter season for a glacier forced by an ice mélange characterized by a back-stress ${\mathit{\sigma }}_{b_{\text{max}}}$ $=$ 0.2 MPa and a thickness $h$ $=$ 100 m. The solid black line refers to the front position. Red diamond corresponds to day 46 (before winter season, before a calving event), yellow diamond to day 140 (during the winter season), blue diamond to day 161 (after the winter season, after a calving event). (b) Variation in the damage criterion $\mathit{\chi }$ at the upper surface along the flow line using the same colour code. (c) Variation in the stress intensity factor. Coloured curves show the stress intensity factor computed at the depth at which the damage threshold $D$ $=$ $D_{\text{c}}$ is reached (solid curves) and at sea level (dashed curves). The horizontal dashed black lines represent ice toughness $K_{\text{Ic}}$ $=$ 0.2 MPa m${}^{\mathrm{1}/\mathrm{2}}$ and arrest criterion $K_{\text{Ia}}$ $=$ 0.1 MPa m${}^{\mathrm{1}/\mathrm{2}}$.

[Figure omitted. See PDF]

Following damage to the ice, three criteria have to be fulfilled to trigger calving: the condition on damage contour $D$ $=$ $D_{\text{c}}$, the initiation of fracture propagation at a depth given by the $D_{\text{c}}$ contour, and propagation to sea level. These criteria are shown in Fig. c, 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_{\text{c}}$) or the stress intensity factor is lower than the propagation threshold ($K_{\text{I}}$ $<$ $K_{\text{Ic}}$). In contrast, the solid red and yellow lines show that the surface is sufficiently damaged to reach the criterion $D$ $=$ $D_{\text{c}}$ over a larger surface area. Wherever the stress intensity factor becomes higher than the ice toughness ($K_{\text{I}}$ $>$ $K_{\text{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_{\text{I}}$ $>$ $K_{\text{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 begin. 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. b. 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 . These authors observed no remarkable changes in the position of the front unless they applied a back force of 50.0 $\times$ 10${}^{\mathrm{7}}$ N m${}^{-\mathrm{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 $\times$ 10${}^{\mathrm{7}}$ N m${}^{-\mathrm{1}}$). Conversely, the results of clearly agree with ours, as they simulated a comparable advance of the front ($\sim$ 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 by preventing the propagation of the fracture down to the glacier base.

(a) Maximum difference in the along-flow component of the deviatoric stress tensor at the glacier surface for each glacier listed in Table  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  forced by the U2 experiment (maximum melt rate: 6 m day${}^{-\mathrm{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.

[Figure omitted. See PDF]

Differences between floating and grounded termini

Figure a 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. b). Conversely, in the case of floating glaciers, the surface adjustment of the tongue decreased the tensile stress compared with the control run (Fig. a, 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. b).

Concerning the behaviour of the grounded geometry of Helheim Glacier, stated that the melting of the front has relatively little effect on the position of the front, unless the prescribed melt rates are extremely high (up to 20 m day${}^{-\mathrm{1}}$). On Store Glacier, 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 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 ($\ge$ 12 m day${}^{-\mathrm{1}}$). Incorporating melting from underneath the floating ice tongue would probably make the glacier front to retreat.

Sketch of the process suggested for glacier equilibrium destabilization. (a) Control run disturbed by (b) increased melting and (c) the resulting stable geometry. (d) Ice mélange perturbation and (e) resulting stable geometry.

[Figure omitted. See PDF]

In our simulations, another difference appeared between grounded and floating glaciers. In Sect. , 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. b. 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. ). The melting perturbation applied on the glacier front affects the shape of the floating tongue (Fig. b). 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. c). 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. d). 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. e). 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 , 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 multiannual 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.

Acknowledgements

This study was funded by the Agence Nationale pour la Recherche (ANR) through the SUMER, Blanc SIMI 6-2012. All the computations presented in this paper were performed using the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône–Alpes region (GRANT CPER07_13 CIRA: http://www.ci-ra.org). We also thank the CSC-IT Center for Science Ltd (Finland) for their support in Elmer/Ice development. We also thank Andreas Vieli, Joe Todd and Martin Truffer, whose comments and advice greatly improved the quality of the paper. Edited by: A. Vieli