A model for tidewater glacier undercutting by submarine melting

Dynamic change at the marine‐terminating margins of the Greenland Ice Sheet may be initiated by the ocean, particularly where subglacial runoff drives vigorous ice‐marginal plumes and rapid submarine melting. Here we model submarine melt‐driven undercutting of tidewater glacier termini, simulating a process which is key to understanding ice‐ocean coupling. Where runoff emerges from broad subglacial channels we find that undercutting has only a weak impact on local submarine melt rate but increases total ablation by submarine melting due to the larger submerged ice surface area. Thus, the impact of melting is determined not only by the melt rate magnitude but also by the slope of the ice‐ocean interface. We suggest that the most severe undercutting occurs at the maximum height in the fjord reached by the plume, likely promoting calving of ice above. It remains unclear, however, whether undercutting proceeds sufficiently rapidly to influence calving at Greenland's fastest‐flowing glaciers.


Introduction
In recent decades ice loss from Greenland has accelerated due to both decreased surface mass balance and increased ice discharge into the ocean [van den Broeke et al., 2016]. Oceanic warming has been widely implicated as a driver of the latter [Holland et al., 2008;Straneo and Heimbach, 2013], yet process understanding of the interaction of the ocean with glaciers remains embryonic and must be improved to better forecast the contribution of the Greenland Ice Sheet to future sea level rise.
One such poorly understood process is submarine melting, which drives direct ice loss and may also impact glacier dynamics by undercutting calving fronts, amplifying calving, and initiating glacier retreat and acceleration . Submarine melting is thought to be amplified where subglacial discharge emerges at the grounding line, forming a buoyant plume rising up the terminus toward the ocean surface [Jenkins, 2011]. The consequent high ice-adjacent water velocities are expected to enhance heat transfer from ocean to ice [Holland and Jenkins, 1999], potentially driving melting at several meters per day [Slater et al., 2015].
Despite its potential importance, the interaction between submarine melting and terminus shape has not been investigated. Recent side-scan sonar surveys of glaciers in west Greenland have shown that calving fronts can be substantially undercut [Rignot et al., 2015;Fried et al., 2015]. It is thought that the most severely undercut sections arise from plume-driven enhanced melting. Existing models of submarine melting, however, assume vertical and static calving fronts and do not therefore capture the process of melt-driven undercutting.
In this study, we introduce a plume model with a dynamic ice-ocean boundary which evolves in response to calculated submarine melt rates. We investigate the effect of undercutting on plume dynamics and the effect of submarine melting on terminus morphology. In addressing these questions we aim to elucidate a key link between the ocean and tidewater glacier dynamics. along-fjord distance, depth and along-interface distance, respectively, taking z = l = 0 at the grounding line. The shape of the ice-ocean interface is defined by the angle (z), measured from the horizontal. The glacier flows from the left at velocity v. The plume induces submarine melting at rateṁ, which varies with depth, and is defined in the direction perpendicular to the ice-ocean interface. Thus, the position and shape of the terminus is determined by the competing influences of ice velocity and submarine melting. There is no calving in this model.
We use the line plume model of Jenkins [2011], which ensures conservation of volume, momentum, heat, and salt as the plume rises along the ice-ocean interface: Variables are illustrated in Figure 1a, with b, u, T, and S are respectively the plume width, velocity, temperature, and salinity, all assumed uniform in the plume cross section. The plume entrains fjord water of temperature T a and salinity S a at ratėe. The reduced gravity of the plume, g ′ , is defined through an equation of state [Fofonoff and Millard, 1983]. C d is a drag coefficient, and Γ T and Γ S are heat and salt transfer coefficients. Submarine melt rateṁ and interface temperature T b and salinity S b are defined using the standard three-equation formulation [Holland and Jenkins, 1999;Slater et al., 2016] in which melt rate is largely proportional to the product of plume velocity and thermal forcing. Finally, our entrainment parameterization follows Pedersen [1980] and Jenkins [2011], witḣe where = 0.1 is an experimentally determined coefficient [Turner, 1979]. Thus, the rate of entrainment of fjord water into the plume scales with plume velocity and the slope of the ice-ocean interface.
We depart from the model of Jenkins [2011] by allowing the ice-ocean interface to evolve in time. Ice flow moves the interface horizontally forward at velocity v while submarine melting acts in a perpendicular 10.1002/2016GL072374 direction, moving the interface backward and, if the terminus is undercut, upward (Figure 1a). At a fixed depth the horizontal movement of the interface may be expressed as (see supporting information) Some specific parameter choices are required. We assume a terminus fjord depth h = 500 m and idealized ambient fjord conditions (Figures 1b and 1c) which are generally representative of fjords in Greenland, with warm subtropical water overlain by cold and fresh polar water [Straneo and Cenedese, 2015]. We consider Greenland-relevant low (Q = 0.1 m 2 /s) and high (Q = 2 m 2 /s) subglacial discharge scenarios, which would, for example, be achieved with respective discharges of 25 m 3 /s and 500 m 3 /s emerging from a channel measuring 250 m in the across-fjord direction. For each discharge, we choose initial plume width and velocity ensuring that the plume has a balance of buoyancy and momentum [Slater et al., 2016]. Initial plume temperature and salinity are set to the in situ freezing point and zero, respectively, appropriate for the emerging subglacial discharge. Numerical details are given in the supporting information.
A few points should be noted regarding the idealized nature of this model. We consider a line plume, neglecting across-fjord variability in plume dynamics. As such this model is only strictly appropriate for "broad" subglacial channels (i.e., channels with significant size in the across-fjord direction), a fact which is tacit in many similar models [e.g., Jenkins, 2011;Sciascia et al., 2013]. A definition of broad is discussed in the supporting information; for h = 500 m, a channel breadth > 250 m might be considered broad, while for h = 250 m, a breadth > 125 m would suffice. The dimensions of real channels are not well known at present; however, observations of subsurface terminus morphology do show undercut caverns measuring a few hundred meters in breadth [Fried et al., 2015].
We assume that the plume follows the shape of the undercut terminus. In particular, our plume model does not apply if the terminus becomes overcut, due to the possibility of the plume detaching from the ice. We assume that ice velocity is invariant with depth, likely a good approximation for fast-flowing tidewater glaciers [Meier and Post, 1987], and neglect other ice-dynamic effects, such as adjustment due to buoyancy or acceleration in the along-flow direction.

Effect of Terminus Shape on Plume Dynamics
We first isolate the feedback of terminus shape on plume dynamics by simulating plumes at static calving fronts of various shapes ( Figure 2a) motivated by observations (for example, the ∼45 ∘ uniform undercutting identified in Fried et al. [2015]). It is immediately clear that, when plotted as a function of depth, the ice-perpendicular submarine melt rateṁ is similar for the various terminus shapes ( Figure 2b); the melt rates differ by at most 20% near the grounding line for the convex versus vertical terminus but generally by < 5%. This similarity arises due to a lack of feedback of terminus shape on plume dynamics when considered as a function of depth, which can be understood by rewriting equations (1a)-(1d) in terms of the vertical coordinate z (putting d∕dl → sin d∕dz) Thus, when considered as a function of depth, the terminus shape (sin ) influences plume dynamics only through the feedbacks of submarine melting (last terms in equations (4a), (4c) and (4d)) and plume-ice drag (last term in equation (4b)). In fact, for sufficient subglacial discharge (hence within plumes in Greenland), these feedbacks are negligible [Slater et al., 2016]. Plume dynamics are then independent of terminus shape, with the reduced buoyancy of the plume at an undercut terminus being compensated by reduced entrainment, as both buoyancy and entrainment scale with interface slope (equations (1b) and (2)). The neutral buoyancy level and maximum height reached by the plume are therefore little affected by terminus shape, exemplified by the low-discharge plumes all reaching a maximum height close to z = 330 m ( Figure 2b).
We note that use of an alternative Richardson number-dependent entrainment parameterization [Kochergin, 1987;Holland and Feltham, 2006] results in a slightly higher dependence of plume dynamics on shape (supporting information). We also stress that where either subglacial discharge is limited or sin is small, such as in winter or at ice shelves, terminus shape is expected to influence melt rate [Magorrian and Wells, 2016].
Returning to plumes at tidewater glaciers, it is crucial to note that more undercut calving fronts have a larger surface area ( Figure 2a). Therefore, a better measure of the effect of melting on the glacier, which takes into account the surface area, isṁ∕ sin . This might be thought of as a "horizontal melt rate" as it is the horizontal displacement of the ice due to melting and can therefore be compared directly to the ice velocity (equation (3)). Due to larger surface area, more undercut calving fronts have higher horizontal melt rate (Figure 2c). Total ablation due to submarine melting, defined as ∫ h 0ṁ ∕ sin dz, is therefore higher at more undercut calving fronts (Figure 2d). To emphasize once more, terminus shape has little effect on ice-perpendicular melt rateṁ (Figure 2b), but undercutting increases the impact of melting (Figure 2c) due to increased ice surface area.
For the remainder of this paper we neglect the feedback of terminus shape on ice-perpendicular melt rateṁ but do account for the increased surface area associated with undercutting. When considering the effect of melting on terminus shape we therefore use a fixed ice-perpendicular melt rate profile calculated for a vertical terminus.

Stationary Calving Fronts
It is instructive to consider whether a terminus may reach a state in which ice flux is balanced by submarine melting. At such a terminus, the horizontal melt rate would equal the ice velocity, with the terminus shape then given by Based on equation (5), if at any depth ice-perpendicular melt rateṁ exceeds ice velocity v, there is no stationary solution and the terminus will retreat. Conversely, if v exceedsṁ, a balance may still be achieved by the terminus becoming undercut; this gives a larger surface area of ice undergoing melting and thus higher total ablation by melting. Some example stationary calving fronts are plotted in Figure 3. Considering the low-discharge case (Figure 3a), the plume does not reach the fjord surface, and thus, melt rates near the surface are zero (Figure 2b). In this scenario, a balance of ice velocity and melting cannot be achieved in the upper portion of the terminus (plotted in dash-dotted lines in Figure 3a), and here removal of ice would have to occur by calving or submarine melting not driven by the plume. Where the plume reaches its maximum height the ice-ocean interface approaches horizontal (Figure 3a) as melt rates are low and a large surface area is required to balance the ice flux. The calving fronts are least severely undercut at the depth of maximum melting (z ∼ 50 m). In Figure 2b, melt rates near the grounding line are low because the plume is cold when emerging into the fjord. At a stationary terminus these low melt rates are compensated by a large surface area, and thus, the ice becomes more undercut close to the grounding line ( Figure 3a).
Similar considerations apply to the high-discharge cases (Figure 3b) where the plume reaches the fjord surface so that the ice flux may be balanced by melting over the full depth. As already motivated by Figure 2c, these stationary examples illustrate how the effect of melting depends on both the melt rate magnitude and on the slope of the ice-ocean interface. Under the same melt rate profile (Figure 2b), both calving fronts in Figure 3b are stationary despite differing ice flux. Finally, the stationary shapes obtained depend on the fjord stratification assumed (Figures 1b and 1c) with unstratified fjord conditions giving less concave stationary shapes (supporting information).

Effect of Submarine Melting on Terminus Shape and Position
We now model the time evolution of terminus position and shape under plume-driven submarine melting. We run each simulation for a typical melt season length of 100 days assuming for simplicity that subglacial discharge is constant in time. By assuming that terminus shape does not feedback on ice-perpendicular melt rateṁ, equation (3) becomes independent of equations (1)-(2). If the plume emerges from a subglacial channel, then near the grounding line x∕ z will be large and positive; thus, momentarily neglecting ice velocity, equation (3) gives which is a one-dimensional advection equation with advection velocityṁ. Solution of equations of this form requires an initial shape x(z, 0) and a boundary condition on the position x(0, t) or equivalently the slope x∕ z(0, t) at the bottom (but not top) of the domain. We choose the latter such that solution of equation (3) requires an initial terminus shape and specification of the slope of the ice-ocean interface near the grounding line through time.
It is not obvious how to specify these conditions. A vertical terminus would seem the default choice for initial shape, but low melt rates near the grounding line ( Figure 2b) will then result in an overcut toe rather than the severely undercut features that have been observed [Fried et al., 2015] and which we seek to model. One solution is to include a form of "subglacial channel" in the initial shape, which is physically reasonable given that we expect severely undercut regions to be associated with plumes emerging from subglacial channels. Thus, below the depth of maximum melting, our simulations begin with a subglacial channel shape taken from one of the stationary solutions previously described. Above maximum melt depth, the initial terminus shapes merge smoothly to vertical. The stationary solution also sets the slope boundary condition at the grounding line. Full details are in the supporting information, where we also show that the evolution of terminus shape is insensitive to the details of the imposed subglacial channel. Figures 4a-4c show evolution of calving fronts for the low-discharge scenario, with initial subglacial channel shape from the v = 2.5 m/d stationary solution (Figure 3a). Above the plume maximum height, melt rates are zero and the ice advances at the ice velocity. Below, the front becomes undercut because melt rates increase closer to the grounding line, until the depth of maximum melting. Below maximum melting, melt rates decrease, but this does not result in a protruding toe because the low melt rates are compensated by the increased surface area associated with the channel. After 100 days, terminus shape approaches that in Figure 3a, blue, which is stationary for an ice velocity of 2.5 m/d. Thus, for an ice velocity of 2.5 m/d, the grounding line is stationary, while for ice velocities of 3.5 and 10 m/d the grounding line advances at 1 and 7.5 m/d, respectively (Figures 4a-4c).
Figures 4d-4f show evolution for the same subglacial discharge but with an initial subglacial channel from the v = 3.5 m/d stationary solution (Figure 3a). The channel shape propagates upward such that the terminus becomes substantially more undercut and retreated than in Figures 4a-4c. Terminus shapes in Figures 4d-4f approach that in Figure 3a, red, which is stationary for an ice velocity of 3.5 m/d. Thus, the grounding line either retreats, is stationary or advances when the ice velocity is less than, equal to or exceeds 3.5 m/d (Figures 4d-4f ). Figures 4g-4l show similar results for the high-discharge case, with initial channels included from the v = 6.2 and 8 m/d stationary solutions (Figure 3b).

Controls on Submarine Melting
Previous studies investigating controls on submarine melting have considered fjord conditions, fjord depth, subglacial discharge, and subglacial hydrology [Xu et al., 2013;Sciascia et al., 2013;Carroll et al., 2015;Slater et al., 2015Slater et al., , 2016. This study adds terminus undercutting, demonstrating that for broad subglacial channels and under a commonly used entrainment parameterization, terminus undercutting does not significantly affect ice-perpendicular submarine melt rates. Total ablation by submarine melting, however, generally increases with the degree of undercutting due to the increased surface area undergoing melting. Thus, the impact of submarine melting depends not only on absolute melt rate but also on the slope of the interface undergoing melting.
Submarine melting may therefore account for a greater fraction of terminus mass balance where calving fronts are severely undercut. Equivalently, the same vertical profile of melting can balance a higher ice flux at a more undercut terminus. In principle, this provides a mechanism for melting to balance ice velocity even when ice-perpendicular melt rates are smaller than the ice velocity. The significance of this mechanism depends on the degree of undercutting and the profile of melting; undercutting increases total melt rate by between 3% and 36% in the idealized examples considered here (Figure 2d).
The importance of shape in determining total melt rate means that in our model, evolution of the terminus is sensitive to the slope boundary condition applied at the grounding line. Specifying a boundary condition with closer to 0 ∘ leads to a more undercut terminus which undergoes greater retreat (Figure 4). Therefore, if we are to simulate the effect of submarine melting on glacier dynamics using our model, we need a means of setting this condition. One might argue that this could be achieved by considering the dynamics of subglacial channels [Röthlisberger, 1972;Ng, 1998]. While channel dynamics are probably important near the grounding line, it seems unlikely that channel dynamics should exert such influence over the full height of the calving front.
It is perhaps more likely that the sensitivity to this condition arises from our application of a plume model arbitrarily close to the grounding line, a region where channel dynamics should dominate over plume dynamics. Additionally, the confinement of the plume when close to the grounding line may reduce entrainment, a factor which is not considered in current models. A more sophisticated approach which captures the transition from subglacial channel to plume may be needed to solve these issues.

Inferring Melt Rate From Shape
For a terminus with a steady state shape, equation (5) may provide a means of determining vertical variation in submarine melting, asṁ∕ sin must be constant with height (assuming also that ice velocity does not vary with depth). Given knowledge of terminus shape sin , potentially from side-scan sonar, the vertical variation in melt rate can be calculated. Thus, parts of a terminus with little curvature are indicative of depth-invariant melting (e.g., Figures 2b and 3b; 100-300 m) while high curvature is indicative of melt rates changing rapidly with depth (e.g., Figures 2b and 3b, near grounding line and fjord surface). Determining the absolute melt rate, rather than just vertical variation, requires additional information on terminus movement as a terminus could achieve a steady state shape but not position (e.g., Figure 4c). This argument relies on the assumption of steady state shape. In reality subglacial discharge varies in time, and calving may also impact terminus shape. At glaciers with stable terminus positions and infrequent subsurface calving [e.g., Luckman et al., 2015;Fried et al., 2015], there may, however, be sufficient time for the terminus to approach steady state. The described method of determining variation in melt rate is independent of entrainment or melt rate parameterizations and so could prove useful in constraining plume models.

Implications for Glacier Dynamics
The presence of subglacial channels appears to play a key role in determining terminus shape, allowing the calving front to become (or remain) undercut close to the grounding line despite low submarine melt rates, and guiding the calving front toward the stationary shapes. Modeled calving fronts become severely undercut where plumes reach their maximum height, whether at depth (Figures 4a-4f ) or the fjord surface (Figures 4g-4l). Since plumes often find neutral buoyancy at a strong pycnocline [Carroll et al., 2015], this is a likely depth for severe undercutting. Calving of ice above may occur when surface crevasses intersect the undercut terminus cavity [Fried et al., 2015], providing a coupling between calving dynamics and fjord stratification. Such a mechanism is perhaps most likely to occur at glaciers with infrequent subsurface calving, giving sufficient time for the calving front to approach a stationary shape. There is some evidence for this mechanism at glaciers in Alaska, Svalbard, and Greenland [Bartholomaus et al., 2013;Chauché et al., 2014;Luckman et al., 2015;Fried et al., 2015].
When ice velocity substantially exceeds submarine melting, as is likely at Greenland's fastest-flowing glaciers, the front may advance hundreds of meters over a period of weeks (e.g., Figure 4i). During this period, other calving processes may come into play [Benn et al., 2007], such as the full-depth buoyancy-driven calving identified at Helheim Glacier , which may occur before the terminus has had time to become substantially undercut or approach a stationary shape. Ice modeling does not yet provide consensus on the degree of undercutting required to influence full-depth calving, if indeed it does at all [O'Leary and Christoffersen, 2013;Todd and Christoffersen, 2014;Cook et al., 2014]. Alternatively, regularly resetting a terminus to vertical by full-depth calving could increase the likelihood of developing a protruding toe, which may promote buoyancy-driven calving [Wagner et al., 2016].

Limitations and Future Directions
We have here considered a line plume which is appropriate only for broad subglacial channels. Consideration of point sources is an obvious next step but would be considerably more complex due to difficulties in capturing the interaction between ice shape, entrainment, and plume structure in three dimensions. The sensitivity of our results to the prescribed grounding line boundary condition needs investigation, perhaps by detailed consideration of the subglacial channel to plume transition. Finally, as with all similar studies, we depend on parameterizations of submarine melting and entrainment, which remain untested in Greenland. Detailed field observations of plumes and terminus morphology are essential to constrain these parameterizations.

Conclusions
We have presented a model coupling buoyant plume theory to an evolving ice-ocean interface, facilitating study of feedbacks between terminus undercutting and submarine melting. We emphasize that the model and conclusions apply only where subglacial discharge emerges from a broad channel.
We find that undercutting has only a weak effect on plume dynamics and local submarine melt rate. All else being equal, local submarine melt rates will therefore not differ between vertical and undercut calving fronts, and undercutting will not affect the height reached by a plume. Undercut glaciers will, however, experience greater total ablation by submarine melting due to increased ice-ocean surface area. In consequence submarine melting may play a greater role in terminus mass balance and therefore glacier dynamics where calving fronts are severely undercut.
We present steady state terminus shapes in which melting and ice velocity are balanced and show how terminus shape may in some circumstances be used to infer vertical variation in melt rate. Time-dependent simulations suggest that the presence of subglacial channels plays a key role in driving severe undercutting. Given sufficient time, calving fronts become strongly undercut at the maximum height reached by a plume, which may promote calving of ice above and define calving style at slower-flowing glaciers. It remains unclear, however, whether undercutting proceeds sufficiently quickly to influence the full-depth calving characteristic of Greenland's fastest-flowing glaciers.