Explorer Effect of near-terminus subglacial hydrology on tidewater glacier submarine melt rates

Submarine melting of Greenlandic tidewater glacier termini is proposed as a possible mechanism driving their recent thinning and retreat. We use a general circulation model, MITgcm, to simulate water circulation driven by subglacial discharge at the terminus of an idealized tidewater glacier. We vary the spatial distribution of subglacial discharge emerging at the grounding line of the glacier and examine the eﬀect on submarine melt volume and distribution. We ﬁnd that subglacial hydrology exerts an important control on submarine melting; under certain conditions a distributed system can induce a factor 5 more melt than a channelized system, with plumes from a single channel inducing melt over only a localized area. Subglacial hydrology also controls the spatial distribution of melt, which has the potential to control terminus morphology and calving style. Our results highlight the need to constrain near-terminus subglacial hydrology at tidewater glaciers if we are to represent ocean forcing accurately.


Introduction
Observations of the mass balance of the Greenland ice sheet in recent decades have shown significant losses at the coastal margins [van den Broeke et al., 2009;Pritchard et al., 2009], much of which has been attributed to the thinning [Pritchard et al., 2009], speedup [Rignot and Kanagaratnam, 2006], and retreat [Jiskoot et al., 2012] of tidewater glaciers. During this period, coastal waters were observed to warm [Holland et al., 2008;Christoffersen et al., 2011] with water from the Irminger current , raising the possibility that Greenlandic tidewater glaciers reacted sensitively to ocean forcing.
One way ocean forcing would manifest itself is by melting the calving fronts of tidewater glaciers, which might provide a significant direct contribution to glacier mass balance [Inall et al., 2014], and could also amplify calving rates [O'Leary and Christoffersen, 2013]. At tidewater glaciers in Greenland, which typically lie at the end of deep and narrow fjords, melt water from the glacier catchment is expected to emerge into the fjord at the glacier grounding line, thereafter rising as a buoyant plume . In combination with the presence of warm subsurface ocean water in Greenlandic fjords [Straneo et al., 2010] these plumes lead to submarine melting of the calving front.
However, the way in which subglacial discharge emerges at the glacier grounding line remains poorly understood. Subglacial drainage is commonly classified as either channelized or distributed; the former characterized by a small number of large channels, routing water quickly beneath a glacier, while the latter is often conceptualized as a network of linked subglacial cavities, providing a tortuous and slow pathway for transit of meltwater [Fountain and Walder, 1998]. At the terminus of a tidewater glacier, a channelized system would input subglacial discharge to the fjord from one or a few large outlets, while a distributed system would spread the input more uniformly across the grounding line. The effect of this distinction on tidewater glacier submarine melt rates remains largely unstudied and is the subject of this paper.
Submarine melt rates at tidewater glaciers in Greenland have been estimated from hydrographic data, and by theoretical and numerical modeling. Rates calculated using hydrographic surveys range from 0.7 to 10 m/d [Rignot et al., 2010;Sutherland and Straneo, 2012;Xu et al., 2013;Inall et al., 2014]. Jenkins [2011] proposed a one-dimensional model which coupled the theory of buoyant plumes with a melt rate parameterization following McPhee [1992]. Xu et al. [2012] and Sciascia et al. [2013Sciascia et al. [ , 2014 used two-dimensional numerical models to investigate fjord circulation and submarine melt rate, while Cowton et al. [2015] recently combined buoyant plume theory and numerical modeling in a three-dimensional fjord model. Finally, Xu et al. [2013] and Kimura et al. [2014] presented high-resolution models of subglacial discharge plumes in three dimensions.
When comparing their hydrographic data and modeled submarine melt rates, Xu et al. [2013] implicitly assumed a fairly distributed system consisting of channels spaced at ∼150 m intervals across the grounding line. Kimura et al. [2014] suggested that the submarine melt rate per unit subglacial discharge decreases as the subglacial system becomes more channelized, but the system on which this conclusion was based was a single channel of variable shape or two channels in close proximity. In addition, the modeled fjord was unstratified, so that the resulting plumes could never reach neutral buoyancy, which will affect submarine melt [e.g., Xu et al., 2013].
This study investigates the effect of variation in subglacial hydrology on tidewater glacier submarine melt rates in greater detail and aims to capture the potential complexity and significance of hydraulic structure at the glacier terminus. The fjord is stratified in a fashion which is typical of conditions in Greenland. For three fixed values of total subglacial discharge, we vary the subglacial hydrological configuration between the conceptual end-members of "channelized" and "distributed," model the induced near-ice water circulation in three dimensions and study the resulting submarine melt.

Model Setup
We use the MITgcm [Marshall et al., 1997a[Marshall et al., , 1997b in nonhydrostatic configuration to model water circulation in an idealized fjord. The domain is 2 km in length and width, and 500 m in depth. Resolution is 5 m across-fjord and vertically, while along-fjord resolution is 5 m for the 250 m closest to the glacier, thereafter increasing to ∼40 m and remaining constant for the 1 km closest to the ocean end of the domain. In all experiments the vertical motion of the plumes is contained within the high-resolution section of the domain.
Initial temperature (T a ) and salinity (S a ) (Figure 1a) are set throughout the domain with data from Chauché et al. [2014], who report conductivity-temperature-depth casts close to Store and Rink Glaciers. The fjord sides are closed boundaries, while at the ocean boundary we impose a sponge layer [e.g., Sciascia et al., 2013] which restores properties toward the T a ∕S a profiles shown in Figure 1a. The glacier end of the domain consists of a vertical calving front 2 km wide by 500 m high. This boundary is closed, except for subglacial channels at the base of the calving front, through which fresh (0 psu) subglacial discharge is injected at the pressure melting point (−0.29 ∘ C). Melting of the calving front is treated using the "icefront" package [Losch, 2008;Xu et al., 2012] with a slightly modified melt rate parameterization described in section 2.3 below. In common with Xu et al. [2013] and Kimura et al. [2014] we implement a free slip condition on the calving front. Quoted melt rates are temporal averages over a 7000 s period after melt reaches a steady state, which provides a period of sufficient length to smooth out turbulence-induced fluctuations (which are no more than 8% of the mean melt rate).

Subgrid-Scale Mixing
The turbulent nature of subglacial discharge plumes leads to entrainment of surrounding fjord water, with the rate of entrainment strongly affecting plume dynamics, and therefore the submarine melt and fjord circulation induced. At 5 m resolution, turbulence is not fully resolved, so some parameterization of the entrainment process is required. This is achieved using Laplacian eddy diffusion of momentum, heat, and salt. With isotropic model resolution, our diffusivities are isotropic [e.g., Kimura et al., 2014], and further, we set the Prandtl number to 1 [e.g., Sciascia et al., 2013] leaving 1 degree of freedom (K) to set the magnitude of subgrid-scale mixing. In common with previous studies, we calibrate this degree of freedom using buoyant plume theory. Specifically, we model a vertically issuing plume of initial discharge q = 1, 3, 10, 30, 100, or 500 m 3 /s in a cubic domain of side 500 m and compare the resulting plume width, velocity, temperature, and salinity to the buoyant plume theory of Morton et al. [1956]. By varying diffusivity K in increments of 0.01 m 2 /s, we identify the value of K which gives the best fit to buoyant plume theory for each value of q. We find that this value of K is approximately proportional to q 1∕4 and use a power law fit ( Figure 1b  The relationship between model diffusivity and per plume discharge, obtained by comparison of plume theory and MITgcm simulations in a simplified setup and applied throughout our main simulations. Vertical black bars represent ± 0.01 m 2 /s on the best fit diffusivity, which is the increment by which we varied K in the comparison procedure. (c) The relationship between channel discharge, size, and velocity obtained by balancing wall melt and creep closure in a Röthlisberger channel.

Melt Rate Parameterization
The melt rate parameterization at the calving front is the three-equation model [Holland and Jenkins, 1999] which has frequently been used in this setting [e.g., Xu et al., 2013;Sciascia et al., 2013;Kimura et al., 2014]. Free stream properties for input to the melt rate parameterization are taken from the grid cells adjacent to the ice. We use velocity-dependent turbulent transfer coefficients T, where v∕w are the tangential water velocities in the cells adjacent to the ice. Use of this parameterization introduces a minimum velocity U 0 into the melt calculation, motivated as follows. Laboratory experiments [Huppert and Turner, 1980] and theoretical modeling [Wells and Worster, 2008;Jenkins, 2011] suggest that in the absence of subglacial discharge or away from discharge-driven plumes, the calving front should still melt leading to the formation of weak convection cells. This is the regime of melt-driven convection rather than the convection-driven melt which is the main focus of this paper. Sciascia et al. [2013] modeled melt-driven convection in a two-layer stratification, achieving vertical velocities averaging 0.04 m/s. We find that with a model resolution of 5 m, and even at the low values of K used in this study, we do not resolve the delicate plumes formed in the absence of subglacial discharge. Thus, in order to represent this "zero-discharge" or background melt, we impose a minimum velocity U 0 = 0.04 m/s in the melt calculation.

Subglacial Hydrology
Input of a total subglacial discharge Q into the fjord beneath the vertical calving front of the glacier requires the choice of the number of discharging channels n, channel size X, velocity V, and channel shape. In this study we vary n and assume that the channels are approximately semicircular in shape. In order to relate X and V to a discharge q = Q∕n in each channel, we balance wall melt and creep closure in a Röthlisberger channel [Röthlisberger, 1972;Schoof , 2010], giving relationships X ∝ q 6∕7 and V ∝ q 1∕7 (Figure 1c, see also supporting information). Note that this calculation is only applicable to grounded termini; indeed, the assumption of a vertical calving front with discrete subglacial channels may only be relevant to grounded termini.

Description of Experiments
This study aims to determine the effect of variation in near-terminus subglacial hydrology on tidewater glacier submarine melt rates. We use three values of total subglacial discharge Q = 125, 250, and 500 m 3 /s, and vary the subglacial hydrology between the end-members of "channelized" and "distributed" drainage. We split the discharge Q over a number n = 1, 2, 3, 5, 10 or 50 identical channels spaced evenly along the glacier grounding line. These runs are indexed by the parameter = 1∕n. To refer to a certain simulation we introduce some notation; Q 250 0.1 refers to the simulation with a total discharge of 250 m 3 /s emerging through 10 identical channels spaced evenly. We also include experiments in which the discharge is uniformly distributed across the grounding line (identified with = 0) and a number of sensitivity experiments, discussed in section 3.2. A full list of model parameters is given in Table S1.

Results and Discussion
Velocity, temperature, and submarine melt rate distributions for various experiments are shown in Figure 2, with rates spatially averaged over the 2 km wide by 500 m high calving front displayed in Figure 3. When Q = 0 m 3 /s, and at the lowest diffusivity used in any of the experiments (0.025 m 2 /s), water velocities in the cells adjacent to the ice do not exceed U 0 ; thus, U 0 is the velocity used in the melt parameterization. Melt rate is then determined by water temperature and is therefore greater at depth. Averaged over the calving front, we obtain 0.12 m/d of background melt in this case (Figure 3).

Varying Channel Number ( )
Model snapshots of ice-tangential water velocity at the calving front are shown in Figures 2a-2d. After emerging from a subglacial channel, the plumes rise turbulently and spread to form a conical shape. For each plume, maximum velocities are achieved at depth. Weak plumes (q < 2.5 m 3 /s) reach neutral buoyancy before the surface, giving low near-surface velocities (Figure 2d). Water velocities increase with increasing discharge per plume.
For > 0.02, high water velocities are contained within the conical plumes. Thus, in the Q 500 1 simulation, water velocities exceed U 0 over just 29% of the calving front. Furthermore, the plumes remain visibly independent as their spacing is decreased until = 0.2 (Figures 2b and 2c). In contrast, for Q 125 0.02 (Figure 2d) SLATER ET AL.
©2015. The Authors. the plumes coalesce within 200 m above the grounding line, and water velocities exceed U 0 over 92% of the calving front.
Snapshots of ice-adjacent water temperature are shown in Figures 2e-2h. After emerging from a channel, plume temperature increases due to mixing with ambient water. Consistent with plume theory, plumes with a higher initial discharge stay colder until closer to the surface, for example, the Q 500 1 plume (Figure 2e) remains colder than 1 ∘ C for ∼200 m, while the Q 500 0.1 plumes are warmer than 1 ∘ C within 30 m above the grounding line ( Figure 2g).
All of the plumes simulated advect the warm and deep ambient water upward, resulting in warming where the plumes reach the surface. Most of this warm plume water subsequently flows down fjord and away from the ice, though there is also some lateral spreading across the calving front resulting in slight surface warming away from the plumes (e.g., Figure 2e). With many plumes ( = 0.33 or below, Figure 2f ), the cold surface water is replaced by warm plume water across most of the calving front, except when the plumes do not reach the surface (e.g., Figure 2h). In the Q 125 0.02 simulation (Figure 2h), water warmer than 2 ∘ C covers almost the entire calving front and, in general, an increased number of plumes results in a larger fraction of the calving front experiencing warm water.
Submarine melt resulting from the distribution of velocity and temperature is shown in Figures 2i-2l. Elevated melt rates are contained within a conical region above each subglacial channel. Maximum melt rates are generally achieved between 350 m and 450 m depth, coinciding with maximum velocities. In the Q 500 1 experiment (Figure 2i), maximum velocities are obtained at ∼375 m depth, but the plume remains colder until closer to the surface and this shifts the maximum melt location up the calving front to ∼300 m depth. When plumes reach neutral buoyancy before the surface (Figure 2l), melt rates are low near the surface. The melt distributions emphasize the independence of the plumes for ≥ 0.2 (Figures 2j and 2k). In general, the distribution of melt follows that of velocity; areas of the glacier in contact with warm water but low velocities undergo little melt (e.g., Figure 2e and 2i). This highlights that the melt rate parameterization used in this study is strongly dependent on water velocity, a dependence which has yet to be tested with observations from a real system. Spatial averaging of the submarine melt rate gives the results displayed in Figure 3. For all three total discharges considered, the spatially averaged submarine melt rate increases rapidly as the subglacial drainage system becomes more distributed. For Q = 500 m 3 /s, the presence of 50 subglacial channels increases total submarine melt by a factor 4.5 over the single channel case. Such a trend can be motivated theoretically: Cowton et al. [2015] suggest melt rate for a single plume scales with discharge asṁ ∝ Q 2∕5 . If the plumes remain independent, splitting Q over n channels gives a per plume melt rate (Q∕n) 2∕5 and a total melt rate n(Q∕n) 2∕5 = −3∕5 Q 2∕5 . Then for a fixed total discharge but varying number of channels, melt should scale as −3∕5 . Application to our Q = 500 m 3 /s results provides a good fit from = 1 to 0.1, the latter of which is the point in our simulations when the plumes first merge. Finally, the = 0 experiments are effectively two-dimensional ( Figures S3-S5) and give the highest total melt of all, reaching 3.6 m/d for Q 500 0 . We argue that it is unlikely that water would emerge in a sufficiently uniform fashion that this case is realistic, and thus our sensitivity experiments use = 0.02 as a distributed end-member.
The increase in total melt as discharge becomes more distributed can be viewed as arising from the sublinear dependence of melt on discharge for a single plume. In terms of our modeled distributions of SLATER ET AL.
Channel width (m) 20 10 60 n/a n/a −6.5 +8.8 velocity and temperature, it arises because a distributed drainage system leads to both water motion and warmer water over a larger proportion of the calving front. We also note that the spatially averaged submarine melt rate is sensitive to glacier width when there is a single plume, while more distributed cases will be less sensitive to variation in glacier width. The theoretical argument above suggests that melt rate will increase with the number of plumes provided these plumes remain independent. As such, our results apply to "wide" glaciers, where the glacier width is much larger than the width of a single plume. We expect the vast majority of tidewater glaciers in Greenland to satisfy this condition.

Sensitivity Experiments
Numerical modeling of tidewater glacier submarine melt rates involves the choice of several inputs which are poorly constrained. We therefore present a number of experiments designed to test the sensitivity of our results to variation in these inputs (Table 1). As base cases we take the Q 250 0.02 and Q 250 1 simulations. A brief discussion of the results is presented here with more detail in the supporting information.
The sensitivity of our results to the melt calculation parameters T i and C 1∕2 d Γ T,S , and to a uniform shift in ambient temperature by ΔT a at all depths can be understood largely by consideration of the melt calculation alone [Holland and Jenkins, 1999]. Sensitivity to channel width and the slip condition is in line with results from Kimura et al. [2014].
Of relevance specifically to this study are the sensitivities to U 0 , V, K, and resolution. Sensitivity to U 0 is related to the proportion of the calving front which is significantly affected by a plume; thus, the Q 250 0.02 case is insensitive to a doubling of U 0 , while Q 250 1 experiences a 17% increase in total melt (Table 1). Total melt proves relatively insensitive to channel velocity V; change in velocity by a factor 2 affects melt by at most 5% (Table 1).
Change in the diffusivity (K) affects the rate at which the plumes entrain ambient fjord water. We have fitted K to our experiments with a precision of 0.01 m 2 /s, and variation in K by 0.02 m 2 /s affects total melt by at most 5%. We also investigate resolutions of 2.5 and 10 m with diffusivities again chosen to fit plume theory. A decrease in resolution to 10 m results in decreases in total melt reaching ∼15%, while increasing resolution to 2.5 m gives modest increases in melt. In general, none of the sensitivity experiments performed suggest that the main conclusions of this paper would be affected by variation in our choice of model inputs.

Melt Rates
The results of this study suggest that the configuration of the near-terminus subglacial hydrological system is an important factor controlling tidewater glacier submarine melt rates. For a fixed total discharge, we predict that spatially averaged submarine melt rates increase quickly as discharge becomes more distributed. In particular, the modeling suggests that, for the temperature and salinity profiles used in this study, a plume resulting from a single, large subglacial channel is unable to induce spatially averaged submarine melt rates exceeding ∼0.65 m/d, or lower if a wider calving front was considered. However, numerous small plumes of subglacial discharge can induce significant submarine melt but may not reach the fjord surface. Therefore, plumes visible at the fjord surface may not be the dominant contributors to total submarine melt.
The melt rates achieved in this study display order of magnitude agreement with estimates from hydrographic data. However, it is important to stress that significant uncertainty exists in the submarine melt calculation outlined in section 2.3, notably within the turbulent transfer coefficients C 1∕2 d Γ T,S U [e.g., Jenkins, 2011] and regarding the relative importance of temperature and water velocity. It is also important to acknowledge that hydrographic data provide a snapshot of conditions in fjords which display significant short-term variability [Straneo et al., 2010;Jackson et al., 2014], such that calculated heat transport may not be representative of the mean . Finally, we emphasize that this study focuses in high resolution on a small section of fjord adjacent to the glacier terminus, using fixed T a ∕S a profiles. As such, we have neglected the potential effect of variations in subglacial hydrology on wider fjord circulation which might in turn affect melt rates [Jackson et al., 2014;Sciascia et al., 2014].

Calving
Aside from influencing total melt, variations in subglacial hydrology might also influence calving rate and style due to spatially varying melt. In general, within the modeled plumes, melt rates are greater at depth, which could lead to undercutting of the terminus and amplification of calving [O'Leary and Christoffersen, 2013]. Considering across-ice variation in melt rate, isolated plumes (e.g., Q 500 1 ) might lead to the formation of calving bays and unstable headlands, consistent with the calving style observed at Store Glacier [Chauché et al., 2014]. It is also possible that calving front morphology might impact plume dynamics, with the potential for important feedbacks between plume dynamics, submarine melt, and calving. Indeed, the melt rates modeled in this study are not sufficient to directly explain the observed retreat of tidewater glaciers in Greenland; thus, if ocean forcing is the key driver of this retreat, it is likely a result of a close coupling between submarine melt and calving mechanics.

Conclusion
A general circulation model, MITgcm, has been used to model near-ice water circulation and submarine melt rates driven by buoyant subglacial discharge at the terminus of an idealized tidewater glacier. The emergence of discharge at the grounding line of the glacier is varied between the end-members of a "distributed" and "channelized" subglacial hydrological system, focusing on the effect of this variation on submarine melt.
The results suggest that variation in subglacial hydrology results in large changes in both the distribution and total volume of submarine melt. In particular, we find that (i) total melt volume is greater when discharge emerges in a distributed rather than channelized fashion, with enhancement by a factor 5 possible under certain conditions, (ii) strong plumes emerging from large subglacial channels are not, in isolation, able to induce large total melt volumes, but (iii) numerous and distributed small inputs of subglacial discharge are able to drive significant submarine melt. The distribution of melt rate has the potential to influence calving front morphology and calving style, a coupling which remains poorly understood but is likely important for tidewater glacier dynamics. Our results identify a need to constrain near-terminus subglacial hydrology at tidewater glaciers if we are to represent ocean forcing accurately and ultimately to understand and even predict the behavior of Greenland's outlet glaciers. The ocean model and relevant packages are available for download from mitgcm.org. Data used in the production of Figures 2 and 3, and Table 1 are available free of charge by emailing the corresponding author. D. A. S is supported by a NERC PhD studentship. T. R. C is supported by NERC grant NE/K014609/1 awarded to P. W. N and A. J. S.
The Editor thanks Fiamma Straneo and an anonymous reviewer for their assistance in evaluating this paper.