How Oceanic Melt Controls Tidewater Glacier Evolution

The recent rapid retreat of many Arctic outlet glaciers has been attributed to increased oceanic melt, but the relationship between oceanic melt and iceberg calving remains poorly understood. Here, we employ a transient finite element model that simulates oceanic melt and ice break‐off at the terminus. The response of an idealized tidewater glacier to various submarine melt rates and seasonal variations is investigated. Our modeling shows that for zero to low oceanic melt, the rate of volume loss at the front is similar or higher than that for intermediate oceanic melt rates. Only very high melt rates lead to increasing volume losses. These results highlight the complex interplay between oceanic melt and calving and question the general assumption that increased submarine melt leads to higher calving fluxes and enhanced retreat. Models for tidewater glacier evolution should therefore consider calving and oceanic melt as tightly coupled processes rather than as simple, additive parameterizations.


Introduction
The current rapid retreat of ocean-terminating glaciers of the Arctic has been attributed to increased advection of warm ocean currents into the glacial fjords (e.g., Holland et al., 2008;Luckman et al., 2015;Slater et al., 2018;. Subaqueous melt erosion at the glacier terminus by warm water is generally assumed to result in the formation of an oversteepened calving front and therefore increased stresses and calving flux (Benn et al., 2017;Motyka et al., 2013;O'Leary & Christoffersen, 2013). This implies that increased oceanic melt drives enhanced volume loss through calving, which is also exploited in simple parameterizations of frontal ablation in models for tidewater glacier evolution (Amundson & Carroll, 2018;Bondzio et al., 2016;Morlighem et al., 2016). In a recent sophisticated modeling study, Todd et al. (2018) also demonstrated a direct influence of submarine melting on the evolution of Store Gletscher. However, their approach only allowed for a fully vertical calving front after ice break-off, and any buoyant submerged ice was removed as soon as it formed, which is not fully consistent with observations (Hunter & Powell, 1998;Fried et al., 2019;Motyka, 1997;O'Neel et al., 2007;Sugiyama et al., 2019;Sutherland et al., 2019;Warren et al., 1995). The presence of submerged ice induces buoyancy forces (Benn et al., 2007;Warren et al., 2001) that alter the stresses near the calving front and consequently the calving process. Other studies have demonstrated that melt undercutting has a limited effect on calving rates (Cook et al., 2014;Krug et al., 2015). Ma and Bassis (2019) found both an enhancing and a suppressing effect of melt on calving, depending on magnitude and vertical distribution of melt, but their simulation was limited to the onset phase of one calving event. Our own preliminary experiments with a transient Lagrangian ice flow and damage-based calving model showed similar effects and indicated that the relationship between submarine melt and iceberg calving may not be as straightforward as previously thought (Mercenier et al., 2019). However, that study lacked a systematic investigation of the influence of oceanic melt rates on calving flux and frontal volume loss.
In this paper, we aim to better understand the link between oceanic melt, iceberg calving and volume loss for a tidewater glacier geometry that evolves over several years. We use a transient Lagrangian finite element ice flow model that simulates ice break-off using a damage evolution law combined with the application of oceanic melt at the calving front. We investigate the response of an idealized glacier to variations of oceanic melt rate. Further, the effects of seasonality and vertical pattern in oceanic melt rate on the volume loss are evaluated.

Methods
We employ the transient Lagrangian multiphysics calving model developed in Mercenier et al. (2019) in two and three dimensions. This model is implemented using the libMesh finite element library (Kirk et al., 2006).

10.1029/2019GL086769
Key Points: • The effect of oceanic melt on tidewater glacier evolution is investigated using a transient calving model based on damage evolution • Oceanic melt has a complex influence on tidewater glacier evolution and increased melt rates may not necessarily lead to more volume loss • The calving and oceanic melt processes are not additive, which has implications on the forcing of models for tidewater glacier evolution Ice flow is calculated solving the Stokes equations for incompressible fluid flow with power law rheology (Glen's flow law). Material damage D is implemented as an evolving state variable (Pralong & Funk, 2005) that varies according to the damage rate B, a stress measure chosen for damage evolution, a stress threshold th , the power r, and a healing term h (for parameter settings, see Table S1 in the supporting information). Damage in the ice affects its viscosity with A the fluidity parameter, . e the effective strain rate, n = 3 the Glen's flow law exponent, and a regularization parameter. Therefore, damaged ice is softened while undamaged ice remains unaffected (D = 0) and a feedback between damage and the stress and velocity field is created. The model is fully Lagrangian with the state variables stored on the mesh nodes, which avoids issues with numerical diffusion as no advection problem needs to be solved. Details are given in Mercenier et al. (2019).
The model geometry is evolved over time in equal time steps calculating the velocity field and stresses. The state variable "damage" is updated according to equation (1), and elements in contact with the ocean accumulate melt according to their interface area and oceanic melt rate. Ice removal by break-off (calving) and ocean melt is simulated with element extinction once accumulated damage exceeds a critical value or melt accumulation exceeds the element volume. The mesh nodes are subsequently moved according to nodal velocities yielding a deformed geometry. The whole domain is then horizontally moved at a constant speed u in , which is chosen for each experiment to compensate average ice loss and to obtain a quasi-stationary calving front position. This horizontal movement does not influence the shape or volume loss of the modeled glaciers. The gap created at the upstream boundary by the horizontal movement is filled with new undeformed elements with their state variables set to zero. Details are given in Mercenier et al. (2019).
The damage evolution law (equation (1)) was tested with different types of stress measure (Mercenier et al., 2019, equation (8)). Damage evolves when the stress measure exceeds a threshold th , leading to ice weakening and ultimately failure. From all tested stress measures, only the "von Mises" and "von Mises tensile" stress measures produced realistic calving front geometries and significant calving activity (Mercenier et al., 2019). The von Mises stress is thus used throughout this study.
Starting from a rectangular block geometry, all model runs evolved for 5 years with 5,000 time steps of 0.001 year. The initial domain consisted of a simple geometry of 2,000 m length and 200 m thickness, discretized with quadratic isoparametric Q2Q1 elements at a spatial resolution of 10 m. The relative water depth ( = H w H ) for all experiments was set to 75% and the other model parameters are given in Table S1. The ice volume losses occurring through oceanic melt and calving were tracked separately to facilitate the analysis.
Different types of oceanic melt experiments were performed by varying the oceanic melt rate magnitude, its seasonality, and its vertical pattern. Combinations of melt-forcing parameters and their corresponding designations are listed in Table S2 and labeled with corresponding abbreviations. Oceanic melt was either continuously applied at a prescribed rate (denoted as C) or seasonally switched on and off (denoted as S). The maximum melt rate M was varied between 0 and 1,500 m a −1 in intervals of 250 m a −1 (denoted as Mx with x being the applied melt rate) which covers the main range of observed values from Greenland glaciers (Carroll et al., 2016;Rignot et al., 2010). The vertical pattern of oceanic melt was either set as constant with depth (C) or as linearly increasing from zero at the waterline to the prescribed melt rate at the bottom of the fjord (denoted as L), but any arbitrary oceanic melt parameterization as function of position and time can be prescribed in the coupled model. In the following, the term "melt" in general refers to submarine melt at the calving face.

Results
The whole suite of model experiments was performed on the same along-flow rectangular geometry that evolves by ice deformation, calving, and melt for 5 years. During the first 1.5 years the terminus geometries self-adjust from the initial condition and to the imposed melt rates. Therefore, the evolution of volume loss (with units m 3 ∕m) for all experiments is only shown from 1.5 years to the end of the simulation (Figure 2).

Melt Characteristics
In the first set of model experiments, labeled CMC (Table S2), a continuous oceanic melt that is constant with depth was applied. The results of this set are shown in Movie S1 in the supporting information and Figures 1 and 2a.
Different shapes of the calving front ( Figure 1) and volume loss rates (relative to the "no melt" scenario M0; Figure 2) are obtained for the imposed melt rates. The scenarios of high oceanic melt (M > 750 m a −1 ) show the development of an oversteepened calving face below the waterline (Movie S1 and Figure 1) and experience increasing volume losses with enhanced melt rates ( Figure 2a). In contrast, lower melt rate scenarios (M ≤ 750 m a −1 ) result in decreasing volume losses with increasing oceanic melt rates (Figure 2a), and an ice foot below the waterline develops near the base of the fjord that is more prominent for lower melt rates (Movie S1 and Figure 1). For intermediate melt rates, the terminus geometry is undulating, with an overcut in the upper part and an undercut shape toward the grounding line.
Movie S2 shows the evolution of the glacier geometry that experiences seasonal variations in depth-averaged melt (labeled SMC). For the same melt rates as before, smaller absolute and relative volume loss differences are obtained than for the CMC scenarios (Figures 2 and 3). During the melt season, the volume losses and geometries evolve similarly to the CMC scenarios, with the development of oversteepened calving faces ( Figure S1, between 2.6 and 3 years) and increased volume losses for the SMC1000 high melt scenario (Figure 2b). In this scenario, the damage accumulated in the ice below the waterline remains mostly below the threshold for break-off ( Figure S1, between 2.6 and 3 years), and thus oceanic melt almost directly translates into additional volume loss. After switching off oceanic melt, the volume loss briefly decreases ( Figure  S1, between 3 and 3.2 years and Movie S2), and the calving front adjusts to a similar shape and evolution as the "no melt" scenario M0, which exhibits a permanent large ice foot below the waterline. For lower melt rates the seasonal evolution of the front is very similar, although the undercutting during the melt phase is less pronounced and volume losses are slightly reduced.  Table S2). Similar geometries ( Figure S2 and Movie S3) and the same relationship between volume loss and oceanic melt are found as for the SMC scenarios. Compared to SMC the relative differences in volume loss are slightly subdued and a melt rate exceeding 1,250 m a −1 is necessary for a volume loss that is higher than the "no melt" scenario M0. Figure 3 summarizes the key characteristic of our sensitivity experiments, namely that moderate oceanic melt reduces volume loss rates. Volume losses are higher for the "no melt" scenario M0 than for most scenarios with melt. For our choice of model parameters the minimum average volume loss is obtained for scenario CMC750, which is ∼17% less than without melt. Only the scenarios with melt rates exceeding 1,000 m a −1 show a higher volume loss than scenario M0.
The relative contribution of calving and oceanic melt to volume loss is displayed in Figure 4 for the SMC scenarios. For high seasonal melt rates (>750 m a −1 ) the volume loss during the melt season is dominated by oceanic melt, with almost negligible damage-induced calving. During the period without melt, damage evolution quickly recovers and calving thereby rapidly reaches a constant rate again. For lower seasonal melt rates the effect of reduced calving due to melting is also clearly apparent.

Sensitivity to Damage Rate and Ice Thickness
The competition between damage evolution and ice removal through oceanic melt determines the total volume loss, and hence advance or retreat of the terminus. Additional experiments were performed with a doubled damage rate parameter, which showed the same contrasting relationship of oceanic melt with The ice thickness also has a strong control on the volume loss. We modeled the response of glaciers with different initial ice thicknesses (150, 200, and 250 m) to variations of oceanic melt. With increasing ice thickness, the minimum of the average volume loss occurs at higher melt rates (Figure 5b). For the thicker glacier, higher melt rates (>1,500 m a −1 ) lead to higher volume losses than the no melt scenario, similar to the other volume loss curves. Further, the geometries produced for the thicker glacier are similar to the glacier with an initial ice thickness of 200 m ( Figure S3).

The Effect of Melt on Tidewater Glacier Evolution
Melt undercutting is generally assumed to lead to higher calving activity due to the formation of an oversteepened calving face below the waterline (Benn et al., 2017;Hanson & Hooke, 2000;  O 'Leary & Christoffersen, 2013). Our model results show that the effect of oceanic melt on the calving front geometry could lead to a complex behavior of tidewater glacier termini. This complexity seems to stem from the different calving front geometries below the waterline that result from the competing processes of damage evolution and oceanic melt. While calving through damage leads to an overcut of the entire calving face, oceanic melt undercuts the submerged part of the terminus, as outlined in Figure 1 and Movies S1, S2, and S3.
For our choice of geometrical and model parameters the lowest volume loss was found for a melt rate of 750 m a −1 (Figure 3). For such a melt rate, the calving face undulates as a result of the competition between ice removal through damage and melt undercutting. At higher melt rates the volume losses are increasingly dominated by melt, and therefore calving through damage is effectively shut down, as the ice is melted away before critical damage is reached (Figure 4b). At lower melt rates, calving dominates the volume loss ( Figure 4a) and the calving front geometry is characterized by a subaqueous ice foot which is most pronounced for the "no melt" scenario M0. Such subaqueous feet are generated when calving loss above the waterline exceeds the volume loss below (Benn et al., 2007) and have been observed at several calving glaciers (Hunter & Powell, 1998;Motyka, 1997;O'Neel et al., 2007;Sugiyama et al., 2019;Warren et al., 1995).
The presence of an ice foot induces buoyancy forces (Benn et al., 2007;Warren et al., 2001) and enhances stresses at the calving face and around the grounding line (clearly visible in panels M0 to CMC500 of Figure 1). These stresses in turn lead to increased damaging and hence higher volume losses through calving (Figure 4). In contrast, increasing oceanic melt reduces the length of the ice foot faster than it can break off through damage accumulation, thus reducing volume loss rates. These model results demonstrate that moderate melt rates reduce volume loss from the glacier terminus and imply that for low melt rates an increase in fjord temperatures would not necessarily induce enhanced volume loss and terminus retreat.
Large ice feet have rarely been described at ocean-terminating tidewater glaciers (e.g., Motyka, 1997). Only recently, observations from tidewater and freshwater glaciers during summer conditions have been published. These studies illustrate that many calving fronts exhibit ice feet, which occasionally reach lengths exceeding 100 m (Bendtsen et al., 2017;Fried et al., 2015Fried et al., , 2019Slater et al., 2018;Sugiyama et al., 2019;Sutherland et al., 2019;Rignot et al., 2015). The shapes and extents of these documented terminus geometries show strong similarities with our model results. For oceanic melt rates typically estimated in Greenland fjords in summer, the modeled ice feet are reduced to 50 m or less (Figure 1 and Movie S2), and the modeled undulating front shapes are consistent with observations (e.g., Kangerlussuup Sermia, Fried et al., 2015).
Melt undercutting leads to the formation of oversteepened calving faces with large stresses (Figure 1), but calving through damage only occurs above the waterline and hence the loss through calving remains limited. Importantly, the contrasting relationships robustly map onto the two cases of (i) low and distributed oceanic melt rates and (ii) high melt rates from meltwater plumes (Fried et al., 2015;Sciascia et al., 2013;Slater et al., 2015Slater et al., , 2018. Under constant melt conditions, an increase in damage rate or ice thickness leads to an overall larger volume loss, with an increased contribution through ice break-off and consequently a reduced effect of oceanic melt. The transition from a volume loss reduction to an enhancement due to oceanic melt ( Figure 5), which, therefore, also depends on the ice thickness and the time scale of damage evolution (damage rate parameter).

Model Simplifications
The presented results apply for most tidewater glaciers with moderate thicknesses under 300 m but potentially also apply to the larger and thicker tidewater glaciers. Further experiments would be required to investigate the behavior of thicker glaciers.
In this study we aimed at distinguishing the effect of oceanic melt on volume loss alone, without the well-known strong control of bed geometry on tidewater glacier behavior. To use a more realistic geometry than our idealized block, and to validate the results with direct observations, several adaptations would be needed. Sliding of the glacier over the bedrock was implemented as a simple forward movement of the mesh nodes. Alternative implementations of basal motion use a frictional relationship that depends on the water pressure (e.g., Ryser et al., 2014). However, the sensitivity study of Mercenier et al. (2018) indicated that the effect of basal sliding on the stress regime at the calving front, and therefore the damage evolution, is limited, and the results presented here are robust in this regard.
The submarine melt profiles applied on the calving front were idealized to study a range of forcings. They therefore do not represent the exact melt rate distribution of any particular glacier (Ma & Bassis, 2019). The simple melt parameterizations are however sufficient to examine the effect of oceanic melt on tidewater glacier evolution.
Large submarine melt rates are often associated with narrow discharge outlets at the base of the terminus (Fried et al., 2015;Slater et al., 2018). While the runoff input drives the formation of localized plumes that only cover a portion of the terminus extent, fjord-wide circulations likely transport the warm water over most of the calving front (Slater et al., 2018). Our two-dimensional model results presented so far do not capture this spatial variability. To qualitatively assess three-dimensional effects of such localized high melt on the evolution of the calving front, we ran a preliminary simulation on an idealized three-dimensional glacier. The results of this preliminary experiment show the importance of localized high melt rates in plumes, which cut the front back and enhance the volume loss (for details, see Text S1).
Our three-dimensional modeling implies that increased melt in plumes may not only locally enhance volume loss but also trigger enhanced calving and retreat over the entire front, even if background melt rates remain relatively low. This three-dimensional effect (also found by Cowton et al., 2019) is particularly important as plume melt rates increase with subglacial discharge, which in turn is directly linked to surface melt, and has the intriguing consequence that atmospheric warming increases the frontal volume loss . This sensitivity to surface temperature is independent of any warming of ocean water.

Conclusions
Detailed model experiments have shown that oceanic melt has a more complex influence on tidewater glacier evolution than is commonly assumed. At very high oceanic melt rates, increased melt leads to increased volume losses. In contrast, at low to intermediate melt rates volume losses are almost constant or even decrease with increasing oceanic melt.
The complex interplay between ice fracturing processes and oceanic melt and its effect on the evolution of the calving front geometry is illustrated by our transient modeling results. In the scenarios with low or zero melt, the terminus geometry consists of a submerged ice foot that eventually breaks off due to buoyancy forces. Low to intermediate oceanic melt reduces the size of the ice foot, and consequently buoyancy 10.1029/2019GL086769 forces, and thus stabilizes the glacier geometry. Only at very high melt rates the glacier evolution is dominated by the removal of submerged ice with a negative feedback on ice break-off, even with the presence of over-steepened calving faces.
Our model results highlight the necessity to consider iceberg calving and oceanic melt as tightly coupled processes that both influence the terminus geometry, which in turn affects the calving process. Simple parameterizations of calving with oceanic melt or temperatures do not capture the complexity of tidewater glacier evolution and neglect the inverted relationship at low oceanic melt rates. The susceptibility of the terminus to changes in local external forcings from meltwater plumes highlights the need for further investigation of three-dimensional effects. Calibration of model parameters with detailed observations will also be necessary to reproduce the evolution of real tidewater glaciers. This will, together with additional processes, such as the buttressing effect from ice mélange, help to better understand the calving mechanism and its link to the climate system.