Simple models for the simulation of submarine melt for a Greenland glacial system model

Two hundred of marine-terminating Greenland outlet glaciers deliver more than half of the annually accumulated ice into the ocean and play an important role in the Greenland ice sheet mass loss observed since the mid 1990s. Submarine melt may play a crucial role in the mass balance and position of the grounding line of these outlet glaciers. As the ocean warms, it is expected that submarine melt will increase, potentially driving outlet glaciers retreat and contributing to sea level rise. Projections of the future contribution of outlet glaciers to sea level rise are hampered by the necessity to use models 5 with extremely high resolution of the order of a few hundred meters. That requirement in not only demanded when modelling outlet glaciers as a stand alone model but also when coupling them with high resolution 3D ocean models. In addition fjord bathymetry data are mostly missing or are inaccurate (errors of several 100s of meters), which questions the benefit of using computationally expensive 3D models for future predictions. Here we propose an alternative approach built on the use of a computationally efficient simple model of submarine melt based on turbulent plume theory. We show that such simple model 10 is in reasonable agreement with several available modeling studies. We performed a suite of experiments to analyze sensitivity of these simple models to model parameters and climate characteristics. We found that the computationally cheap plume model demonstrates qualitatively similar behaviour as 3D general circulation models. To match results of the 3D models in a quantitative manner, a scaling factor in the order of one is needed for the plume models. We applied this approach to model submarine melt for six representative Greenland glaciers and found that the application of a line plume can produce submarine 15 melt compatible with observational data. Our results show that the line plume model is more appropriate than the cone plume model for simulating the average submarine melting of real glaciers in Greenland.


Introduction
Since the 1990s the decadal loss of ice mass by the Greenland ice sheet (GrIS) has quadrupled (Straneo and Heimbach, 2013), with an average 1993-2010 contribution of 0.33 ± 0.08 mm yr −1 , which is about 10 % of the observed sea level rise during this period (Church and White, 2011;Church et al., 2013).This acceleration of the GrIS mass loss is attributed to increase of surface melt due to atmospheric warming (Khan et al., 2014) and speedup of the marine-terminating outlet glaciers (Rignot and Kanagaratnam, 2006).The latter has been connected, among other factors, to enhanced submarine melting, which in turn is caused by warming of the surrounding ocean (Straneo et al., 2012) and, probably, by increased subglacial water discharge (Straneo and Heimbach, 2013).While ice-ocean interaction potentially plays an important role in recent and future mass balance changes of the GrIS, the understanding of this interaction remains rather poor and represents one of the main source of the uncertainties in future sea level rise projection (Church et al., 2013).The ice sheet models used for the study of GrIS response to global warming and its contribution to sea level rise typically have resolution of 5 to 10 km (Bindschadler et al., 2013), which is too coarse to resolve most Greenland outlet glaciers.Instead, regional modeling at higher resolution is better suited to capture glacier dynamics.As an alternative to costly three-dimensional models, one-dimensional flow line models were convincingly applied to several major outlet glaciers (Nick et al., 2012(Nick et al., , 2013;; Published by Copernicus Publications on behalf of the European Geosciences Union.Lea et al., 2014;Carr et al., 2015).In particular, Nick et al. (2012) simulated with a flow line model the dynamical response of the Petermann Glacier to the abrupt breakup of its floating tongue in 2010 and investigated the influence of increased submarine melting on future stability of the glacier.They demonstrated the strong influence of increased submarine melt rate to the glacier's mass loss.In this study, submarine melt rate was prescribed and held constant.Using the same flow line model, Nick et al. (2013) implemented submarine melt proportional to the ocean temperature outside of the fjord.This study was performed for the four largest outlet glaciers.Under the assumption that the result of the four largest glaciers can be scaled up for the remaining glaciers, Nick et al. (2013) estimate a total contribution of the Greenland outlet glaciers to global sea level rise of up to 5 cm during the 21st century or about 50 % of the maximum expected GrIS contribution due to changes in surface mass balance.For the same period of time but using a three-dimensional ice sheet model, Fürst et al. (2015) estimated the contribution of enhanced ice discharge through outlet glaciers to be 20 to 40 % of the total mass loss.These large uncertainties are associated, among other factors, with the parameterization of the submarine melt rate.Note that in Fürst et al. (2015) the effect of ocean warming was parameterized through enhanced basal sliding rather than explicit treatment of submarine melt.Different approaches have been taken to estimate submarine melt rates of outlet glaciers by using empirical data (Motyka et al., 2013;Rignot et al., 2010): simplified one-dimensional models of line plumes (Jenkins, 1991(Jenkins, , 2011)), axisymmetric plume models (Cowton et al., 2015;Turner, 1973) and numerical two-and threedimensional ocean models (Holland et al., 2007;Little et al., 2009;Sciascia et al., 2013;Xu et al., 2013;Slater et al., 2015).Note that 2-D and 3-D modeling efforts also differ with respect to model formulation, in particular some authors use hydrostatic models (e.g., Holland et al., 2008;Little et al., 2009), while others use nonhydrostatic models (e.g., Sciascia et al., 2013).The experiments studied submarine melt with respect to subglacial discharge and its spatial pattern and vertical ocean temperature and salinity profiles.Additionally, the influence of subglacial discharge on the fjord circulation, which connects outlet glaciers with the surrounding ocean, was investigated with 3-D models (Cowton et al., 2015;Carroll et al., 2015).Different authors considered two main types of subglacial discharge.The first one is uniformly distributed along the grounding line (referred hereafter as "line plume", LP) (Jenkins, 1991(Jenkins, , 2011;;Sciascia et al., 2013;Slater et al., 2015;Xu et al., 2012) while the second one has localized subglacial discharge (the axisymmetric plume, referred hereafter as "cone plume", CP) (Cowton et al., 2015;Turner, 1973;Slater et al., 2015;Xu et al., 2013).The CP approach is motivated by the observation that a significant fraction of subglacial discharge during the melt season emerges through one or several channels underneath the glacier (Rignot et al., 2015;Stevens et al., 2016;Sole et al., 2011).All of the above-mentioned 2-D and 3-D model simulations show, in agreement with previous theoretical studies, that submarine melt strongly depends both on the ambient water temperature and the magnitude of subglacial discharge.However, different modeling studies revealed the complex dependence of submarine melting on temperature.Sciascia et al. (2013) investigated tidewater glaciers and found a linear dependence of the submarine melt rate on ambient water temperature above freezing point.In contrast, Holland et al. (2008) and Little et al. (2009) found a quadratic dependence on temperature under large ice shelves.Xu et al. (2013) detected that this relationship of melt rate to thermal forcing depends on the amount of subglacial discharge released through a single channel at a tidewater glacier: the melt rate dependence to temperature has a power of 1.76 for small discharges and is lower for higher discharge.Slater et al. (2016) found a power law dependence of melt rate on discharge, with the exponent 1/3 for both the CP and the LP models in a uniform stratification.For a linear stratification their study shows that the exponent enlarges to 3/4 for the CP model and to 2/3 for the LP model.A change in power law could also be detected by Xu et al. (2013).They determined an exponent of 0.5 at high discharge and 0.85 at low discharge for the CP.Simulations with 3-D models show the strong dependency of CP melt rate on stratification or other environmental factors, with maximum melt rate near the surface (e.g., Kimura et al., 2014, unstratified) or close to the bottom (e.g., Slater et al., 2015;Xu et al., 2013, stratified).While experiments with high-resolution (several to 10 m) nonhydrostatic 3-D ocean models demonstrate their potential to simulate rather realistically turbulent plumes and melt rates of marine-based glaciers, such models are too computationally expensive for modeling of the entire Greenland glacial system response to climate change on centennial timescale.An alternative is to use a method for submarine melt based on a simplified plume model (Jenkins, 2011;Cowton et al., 2015).Such a plume model can then be used to calculate submarine melt in, e.g., 1-D ice stream models.This would represent a step forward compared to a rather simplistic treatment of submarine melt used in previous works (e.g., Nick et al., 2013).The main purpose of this paper is to investigate the applicability of the simple plume models to simulation of melt rate of real glaciers in Greenland.To this end we first compared both cone and linear plume models with the available results of simulations from high-resolution 3-D ocean models.Then we compare results of plume models with the empirical estimates of submarine melt from several Greenland glaciers.The paper is organized as follows.The two versions of plume model are described in Sect. 2. We then study the plume model's sensitivity of simulated submarine melt rate to ocean temperature and salinity, the amount of subglacial discharge and the geometry of the calving front of the glacier itself.Results of simulations with the simple plume models are compared to results of numerical experiments with 3-D and 2-D ocean models in Sect. 4.
In Sect. 5 we compare our simulations to empirically estimated submarine melt rates for several selected Greenland glaciers.Finally, in Sect.6 we discuss the applicability of the plume model for the purpose of developing a comprehensive Greenland glacial system model.

The plume models
A plume model describes buoyancy-driven rise of subglacial meltwater after it exits subglacial channels, until it reaches neutral buoyancy near the surface.Two counteracting processes control its evolution: (a) submarine melting of the ice-ocean interface under the floating tongue (if any) and upwards along the calving front and (b) turbulent entrainment and mixing of surrounding fjord water.These processes act to maintain or reduce plume buoyancy, respectively.Subglacial meltwater discharge Q sg for a glacier can be estimated from surface runoff and basal melt over the catchment area of the glacier.How this discharge is distributed along the grounding line, however, is in general not known.It is believed that at least during the summer season, most of the subglacial discharge occurs through a network of channels (Chauche, 2016;Rignot and Steffen, 2008;Rignot et al., 2015;Schoof, 2010) but their precise number for different glaciers and relative importance is not known and can change throughout the season.We investigate two situations.The LP model corresponds to the simplest assumption that Q sg is uniformly distributed along the grounding line (Fig. 1a), while the CP assumes point-wise release of meltwater (Fig. 1b), i.e., from a channel whose dimensions are small compared to the plume diameter.Furthermore, the CP model assumes that a self-similar half-conical form is maintained.Note that there can be a number of discretely distributed plumes along the glacier, implying numerous CPs.

Model equations
Both models are formulated in one dimension, x, which is the distance from the grounding line upwards along the glacier front, or under the ice shelf, and depends on the glacier shape, described by its slope α (Fig. 1).The model equations are written under the assumption that the plume is in equilibrium and therefore do not explicitly account for time.All model parameters and their descriptions are listed in Table 1.

Line plume
The LP model after Jenkins (2011) accounts for a uniformly distributed subglacial discharge along the grounding line of a glacier (Fig. 1a).Far enough from the lateral boundaries, it assumes invariance by translation along the grounding line, so that the resulting equations only depend on x with d() dx = () :

Plume
Water line  Jenkins (2011).Subglacial freshwater flux q sg , which is uniformly distributed along the grounding line, enters the fjord, where it forms a plume that rises up due to buoyancy.The plume is described explicitly with its temperature T , salinity S, thickness D and velocity U .It rises along the ice shelf, following the shelf's slope α = 90 • − β, until it either reaches the water surface or has zero velocity due to the loss of buoyancy.The ambient water with salinity S a and temperature T a entrains into the plume with entrainment rate ė.Melting ( ṁ) occurs on the glacier front and adds to the plume buoyancy with water of the temperature T b and salinity S b .(b) Conceptual scheme of two-dimensional CP model modified after Jenkins (1991) and Cowton et al. (2015).Subglacial discharge enters the fjord, localized via a channel.The plume geometry is described as a half cone and the entrainment occurs around the arc.The subglacial

Bedrock
2 , where D 0 is the initial radius and U 0 is the initial velocity.
where the plume state variables D, U , T and S stand for its thickness, velocity in the x direction, temperature and salinity, respectively, all dependent on x.Equation (1) describes the conservation of volume flux q = DU (expressed per unit length in the lateral direction, i.e., m 2 s −1 ), which can increase by the entrainment of ambient seawater ė and by melting ṁ of ice from the glacier front.The momentum flux (Eq.2) is based on the balance between buoyancy flux and the drag C d U 2 of the glacier front.The buoyancy flux is proportional to the relative density contrast ρ ρ 0 between plume water and ambient water in the fjord (subscript a), parameterized in linear form as β S (S a − S) − β T (T a − T ), with coefficient β S and β T indicated in Table 1.The drag also results in a turbulent boundary layer (subscript b) at the icewater interface, where melting occurs, and heat and salt is exchanged by (turbulent) conduction-diffusion.The equations for T and S (Eqs.3, 4) account for the entrainment of ambient water and the addition of meltwater as well as for www.the-cryosphere.net/12/301/2018/The Cryosphere, 12, 301-323, 2018 where the subscript i for temperature and salinity refers to the inner ice, and c is the specific heat capacity.The system is closed by an expression of the freezing temperature T b , which can be linearly approximated as a function of depth Z (Z < 0) and salinity of the boundary layer S b : with coefficients λ i listed in Table 1.For a straight wall, Z = Z 0 + x • sin(α), where Z 0 is the negative depth at the grounding line (x = 0).Solving for Eqs.
(5-7) yields a second-order polynomial equation for the melt rate ṁ, as a function of plume state variables.Note that Jenkins (2011) also uses an approximation of the melt rate equations, which resolves in ṁ = M 0 U (T − T f ), where T − T f is the plume temperature above the plume freezing point, and M 0 is a slowly varying function of ice temperature below plume freezing point.After Jenkins (2011), with a simplified formulation for the heat balance at the ice-ocean interface, M 0 varies from 6 × 10 −6 to 7 × 10 −6 ( • C) −1 over a T i range from −20 to −2 • C and a plume freezing temperature calculated with Eq. ( 7) for the plume salinity varying from 0 to 35 psu (fast entrainment of ambient salinity), resulting roughly in 0 to −2 • C. Numerically, by calculating M 0 with M 0 = ṁ • (U (T − T f )) −1 from our experiments (fixed T i = −10 • C), M 0 is slightly higher and can vary from 7 × 10 −6 to 12 × 10 −6 ( • C) −1 (Fig. A1).
We do not use this approximation in our calculation, but this is nevertheless helpful to interpret some of the results presented in our paper, in particular in quantifying the amount of melt rate and simplifying the melt rate dependence on temperature and subglacial discharge (Appendix A).

Cone plume
The second plume model investigated in this paper is the CP model (Cowton et al., 2015).It differs from the LP model by the geometry of the plume, which resembles half of an upside-down cone (Fig. 1b).In that case, the plume has definite dimensions and fluxes are expressed in full units (m 3 s −1 ).A cross section of the plume is half a disk with area π 2 D 2 , where the length scale D is here the cone radius at a given x.Equations ( 1)-(4) now reform for the CP model by considering melting on the diameter 2D and entrainment around the arc π D: The Cryosphere, 12, 301-323, 2018 where variables, parameters and equations have the same meaning as for the LP model, and the volume flux 2 U is expressed in cubic meters per second.

Numerics
For the differential equation system of Eqs. ( 1)-( 4) and ( 8)-( 11) we choose a classical Runge-Kutta scheme with constant step size x.Note that T a , S a and sin(α) can vary as a function of X, so that the model can be applied on any glacier geometry.Initial conditions for D, U , T and S are needed at X = 0.Then, at every step, we first solve the boundary layer Eqs. ( 5)-( 7) and then compute the increments in the differential equation system of the LP (Eqs.1-4) or CP (Eqs.8-11) model.The procedure is continued until the plume reaches zero velocity or the water surface.The code is written in Python and Fortran for future coupling.

Initial conditions and balance velocity
In the rest of the paper, for simplicity, we refer to the boundary condition at x = 0 as "initial conditions" although the model equations are not time dependent.Since subglacial discharge consists of meltwater, the salinity and temperature of subglacial discharge water can be set to zero (S 0 = 0 and T 0 = 0).We choose T 0 = 0 since the temperature of subglacial water is unknown, but for obvious reasons it cannot deviate significantly from 0 • C. For conditions typical for the Greenlandic environment, we did not find any significant change in melt rate when using the pressure melting point instead of T 0 = 0 • C, since the plume temperature rapidly converges to a balance temperature close to ambient water temperature (Appendix, Fig. A3).
For both LP and CP models, initial dimensions (radius or thickness) D 0 and velocity U 0 are not known, but they are tied by subglacial discharge.In the CP case Q sg = π 2 D 2 0 U 0 (Fig. 1b), while for the LP case the subglacial discharge per glacier width W enters the model equations: 1a).
It turns out that for a given subglacial discharge, simulated velocity rapidly adjusts to a "balance" velocity, regardless of the initial velocity (Fig. 2a), as already noticed by Dallaston et al. (2015).Analytically, the balance velocity (noted U (x) below) is solution of the plume Eqs.(1)-( 2) and ( 8)-( 9) when the transient term U is neglected.The fast adjustment around x = 0 (where plume dimension is small) can be explained by some rearranging into a form analogous to a first-order lin- ear differential equation for U 2 (see Appendix A2.2).The balance velocity is not necessarily constant, but a simple expression for U 0 (at x = 0) can be derived if the plume dimension is expressed as a function of subglacial discharge, and the melt rate is neglected compared to entrainment in the volume flux Eqs. ( 1) and ( 8) (see Appendix A1).We obtain for the LP and for the CP Note that Eq. ( 12) is identical to the velocity derived by Jenkins ( 2011), and Eq. ( 13) is analogous to Eq. ( 5) in Slater et al. (2016), with the addition of the basal drag term.These balance solutions are only valid in the vicinity of the grounding line and velocity might then differ substantially as the plume develops, especially for small subglacial discharge (e.g., Magorrian and Wells, 2016).A more detailed discussion and a full, depth-dependent solution for the LP model are given in the Appendix.
Our sensitivity tests show that initial velocities higher than U 0 lead to maximum melting close above the grounding line of the glacier ("undercutting") while for lower velocities the melt rate increases with height and maximum melting is located further up the calving front (Fig. 2b).We checked that initial velocities smaller than the balance velocity yield very small difference in the cumulative melt rate (Fig. 3), although some differences occur for larger velocities.For the LP model an initial velocity 10 times larger than the balance velocity gives a 10 % higher melt rate while the CP model produces 25 % more melting (Fig. 3).Since the velocities of subglacial discharge are mostly unknown, these results prompted us to use the balance velocity of Eqs. ( 12) and ( 13) as initial condition in all experiments described below, unless stated otherwise.

Comparison between LP and CP models
A direct comparison between LP (defined per unit width of the grounding line) and CP (point-wise) models requires an assumption about a length scale W (for LP) and the number of sources (for CP) over which subglacial discharge is distributed.For the CP model we assumed that the entire subglacial discharge occurs through one channel in the center of the glacier (Q sg = 500 m 3 s −1 ).In the case of the LP model we assumed that the discharge is uniformly distributed over a fjord width W = 150 m, so that q sg = 3.6 m 3 s −1 .This width is about the maximum size of the plume in the CP model, near the surface, for this setting.Results in Fig. 4 show that simulated local melt rate can be higher in the CP model than in the LP model, but cumulative melt rate (i.e., the integral of the melt rate from the bottom and across entire surface area of the glacier front, of width W ) is much higher for the LP model because of the larger surface area over which melting occurs (roughly a factor of 2 in our chosen setting).
We shall see later in this paper (Sect.3.1) that the (local) melt rate in the LP model varies less than linearly with sub- glacial discharge parameter q sg , and thus for a given total discharge, Q sg , cumulative LP-induced melt increases with width.As a result, for a wide glacier (i.e., the glacier which is much wider than the maximum diameter of the CP), the LP model gives much higher cumulative melt rate compared to the CP model, when assuming the existence of a single subglacial channel.The situation when there is more than one channel is discussed in Sect.4.3.

Default experimental setting
In the next sections we perform a number of sensitivity studies with respect to key parameters.To that end we choose a default experimental setting as a benchmark.Unless otherwise stated, we consider a 500 m deep, well-mixed fjord with ambient temperature T a = 4 • C and salinity S a = 34.65 psu (maximal melting conditions for Greenlandic fjord), with total subglacial water discharge of q sg = 0.1 m 2 s −1 for the LP model or Q sg = 500 m 3 s −1 for the CP model (which corresponds approximately to the discharge in August 2010 of the 5 km wide Store Glacier; Xu et al., 2012).Since we apply our model to glaciers in Greenland fjords, most of which do not have a floating tongue (tidewater glaciers), we therefore generally perform experiments for a vertical wall (sin(α) = 1).Default model parameters, including entrainment rate E 0 , are indicated in Table 1.

Subglacial discharge
It is known that melt rate depends strongly on subglacial discharge.In agreement with previous studies (Jenkins, 2011;Slater et al., 2016) our model shows a cubic root dependence of the cumulative melt rate on discharge for the LP and for the CP (for the high discharge range) in a well-mixed environment.However, for small discharges (q sg → 0) the cumulative melt rate converges to a small but not insignificant value that does not obey the power law (Fig. 5).This value represents background melt rate, or "melt-driven convection", which does not depend on discharge and can be representative for winter melt rate when subglacial discharge is very small (Magorrian and Wells, 2016).
To explain this change of power law we undertook a dimensional analysis to obtain theoretical solutions for the LP in a well-mixed fjord (Appendix A).The velocity, which determines the melt rate (Eq.A1, Fig. A1), is initially controlled by subglacial discharge (see Sect. 2.3) and subsequently accelerates as the plume develops.Our analysis shows that the plume acceleration depends on ambient hydrographic conditions, entrainment and glacier front characteristics but is independent from subglacial discharge (Eq.A17-A19).As a consequence, for a given glacier, there is a high-end regime of purely discharge-driven convection (high discharge, high melt, near-constant velocity) and a lowend regime of purely melt-driven convection (zero discharge, low melt, pure acceleration) (Eqs.A21 and A22).We defined a critical discharge q c to characterize the transition between the two limiting regimes (Eq.A23).This critical discharge q c depends on glacier characteristics and especially on the presence and length of a floating tongue.For tidewater glaciers, it is very low (of the order of 10 −3 m 3 s −1 ), so that summer discharge is sufficient to trigger a purely discharge-driven plume.In a glacier with a long floating tongue, q c can be 2 orders of magnitude larger, and thus melt-driven convection may contribute all year long.Our approximate analytical solution for the cumulative melt rate (Eq.A20) is displayed in Fig. 5 and is in reasonable agreement with the line plume result over the full discharge range.
Note, however, that this analysis holds when the plume reaches a dynamic equilibrium.Slater et al. (2016) found deviation from the cubic root power law even for large discharge, for shallow fjord, when the adjustment length scale is large and a significant portion of the plume is not equilibrated.
In addition, as mentioned in the introduction, stratification can change this power law.We also performed experiments with stratification as in (Xu et al., 2013)   The red dashed and blue dashed lines indicate the limiting melt rate regimes M high (Eq.A21) and M low (Eq.A22), respectively.The transition of these regimes is defined by the critical discharge q c (Eq. A23), depicted here with black dashed line.(2013).Both models show an increasing exponent for lower discharge.

Entrainment rate
Entrainment is the mechanism through which the volume flux of the plume increases with distance from its source, as warmer, saltier fjord water mixes into the plume.This leads to more heat availability for melting, but also to decreased buoyancy -and velocity -as the plume gets saltier.Reduced velocity in turn reduces melting (Eq.5, 6) (Carroll et al., 2016) (note the plume also becomes thicker to accommodate for increased volume flux and decreased velocity).In this section we investigate the net effect of these processes on melting for typical plume configurations.
In both plume models, entrainment depends on an entrainment rate parameter and on glacier slope, as E 0 sin α (Sect.2.1).E 0 is not accurately known and can be regarded as a tunable parameter within a certain range of values known from previous work.Laboratory experiments for a pure vertical plume, model studies and theoretical considerations give a broad range for E 0 , from 0.036 to 0.16 (Jenkins, 2011;Kaye and Linden, 2004;Mugford and Dowdeswell, 2011;Kimura et al., 2014;Carroll et al., 2015;McConnochie and Kerr, 2016).Nevertheless, glacier slope sin α can vary by 2 orders of magnitude so that regardless of the value of E 0 within the reported range, tidewater glaciers (sin α ∼ 1) are characterized by a high-entrainment regime, while glaciers www.the-cryosphere.net/12/301/2018/The Cryosphere, 12, 301-323, 2018  Rearranging the LP equations (Eqs. 1 and 2) shows that entrainment acts as an effective drag C d + E 0 sin α (e.g., Eq.A15), with C d E 0 sin α for tidewater glaciers (i.e., controlled by entrainment), whereas C d E 0 sin α for glaciers with a long floating tongue (i.e., controlled by solid friction at the ice interface, insensitive to entrainment) (Fig. A3b).In contrast, plume temperature depends on the mixing ratio of meltwater to entrained water (Eqs.A9-A10), which is close to zero in tidewater glaciers, so that equilibrium plume temperature is nearly equal to ambient temperature in the full range of E 0 values, i.e., already at its maximum potential for melt and insensitive to E 0 (Fig. A3c).In the low-entrainment regime characterizing long floating tongue, the mixing ratio of melt to entrainment is significant, so that temperature is strongly controlled by E 0 (Fig. A3c).
As a result, for the LP model, an increase in E 0 leads to less melting in tidewater glaciers (Fig. 6a), where the plume slows down, but increases the melt rate for glaciers with long floating tongues, where enhanced mixing results in more heat available for melting (Fig. 6b).This effect is particularly strong for warm ambient temperature (Fig. 7) and strong discharge (Fig. 6).LP cumulative melt rate can decrease by over 42 % over the full E 0 range in a tidewater glacier (Fig. 7).In the CP model, the same physics applies and determines local melt rates, but entrainment also influences the plume radius, which grows faster with larger entrainment.As a result, for the CP model, even for tidewater glaciers, where local melt rates inside the CP ( ṁ) decreases with E 0 , the simultaneous effect of a increasing plume area in contact with the ice (2D) still dominates the calculation of the cumulative melt rate (M = 2 ṁ(x) • D(x)dx).
Thus, this leads to the statement that cumulative CP melt rate increases with increasing entrainment factor E 0 (Fig. 6c).

Ambient temperature and stratification
Different fjords are characterized by different temperature and salinity profiles.Since the temperature of the ocean is projected to increase with global warming, dependence of melt rate on ocean temperature is crucial to study glacier response to global warm- ing.Previous experiments with 2-D and 3-D ocean models, as well as analytical solutions (Jenkins, 2011;Sciascia et al., 2013;Carroll et al., 2015;Jenkins, 2011;Magorrian and Wells, 2016;Slater et al., 2016;Carroll et al., 2015), demonstrated the behavior of the cumulative melt rate as a function of the ambient temperature T a .Figure 8 shows for both plume models the dependence of cumulative melt rate on temperature in a well-mixed ambient environment for different values of subglacial discharge.Both models show for small discharge a nonlinear dependence of the melt rate on water temperature.If the discharge is very small, melt-driven Table 3. Numerical determination of the power law β for a melt rate dependence of m ∝ T β a with the thermal forcing T a = T a − T af for tidewater glaciers.The exponent was derived for the LP and CP model for lower and higher discharge ranges.A comparison to an analytically derived value of β for limiting discharge (low and high) ranges from literature and our study is listed for some cases additionally.convection dominates, and ambient properties determine the melting process (Slater et al., 2015).We assume a power law dependence of the cumulative melt rate per glacier area to the thermal forcing, i.e., m ∝ T β a , where T a = T a − T af and T af is the freezing temperature of the sea water at the fjord bottom.As listed in Table 3 we find that the exponent β increases with lower discharge: from 1.2 (high discharge q sg = 0.1 m 2 s −1 ) to 1.8 (q sg = 10 −6 m 2 s −1 ) for the LP and from 1.2 (high discharge 500 m 3 s −1 ) to 1.5 (low discharge Q sg = 0.005 m 3 s −1 ) for the CP.For the LP the range of this increase compares well to analytical solutions while the CP model can only be compared to the analytical solution of high discharge range (Table 3).The exponent also increases when using a realistic stratification (Fig. 9b).For the LP, we calculated an exponent of 1.2 for high discharge and 1.4 for low discharge, while the CP model shows a similar increase from 1.1 to 1.3.Carroll et al. (2015) showed that plume theory gives a good approximation of the outflow height for 3-D nonhydrostatic plume model but nevertheless the exponents differ slightly in the experiment by Xu et al. (2013).

Glacier front angle for the line plume
The glacier front angle sin(α) linearly impacts buoyancy (Eq.2) and entrainment.For glaciers with a floating tongue, and therefore a smaller angle (sin(α) 1), entrainment is reduced and so the temperature of the plume (Fig. 10b).The dependence of melt to the slope of the glacier for small discharges has been derived by Magorrian and Wells (2016).A glacier with a long floating tongue, and therefore a small angle, has a smaller average melt rate than a tidewater glacier but a higher cumulative melt rate (Fig. 6b).These high cumulative melt rates (Fig. 10) occur due to the longer distance under a floating tongue over which melting occurs.Furthermore, for long floating tongues, the plume velocity accelerates along the shelf (10b), with a square root dependency on Cumulative melt the distance (Eq.A19).This is consistent with our analysis for a LP model in a well-mixed environment, which demonstrates that glaciers with a long floating tongue have a high critical discharge (Eq.A23) and thus a larger contribution of melt-driven convection compared to tidewater glaciers.
However, for small α none of the plume models are applicable along the total shelf because they do not take into account Coriolis force and plume thickness is limited by the Ekman layer depth (Jenkins, 2011).Therefore, for the total shelf area, the plume models likely strongly overestimate plume velocity and melt rate (see more in Sect.5.1).
www.the-cryosphere.net/12/301/2018/The Cryosphere, 12, 301-323, 2018 The melt rate of the corresponding temperature profile is displayed in the same color as well as in the same style (dashed, dotted or solid) for the corresponding discharge.Note that a very high discharge (q sg = 0.5 m 2 s −1 ) is needed for the plume to reach the surface.
For each discharge value the corresponding cumulative melt rate is depicted (b) as a function of the thermal forcing ( T a = T a − T b , Eq. 7) at the grounding line.For ṁ ∼ T β a we found β values of 1.2 for (q sg = 0.5 m 2 s −1 ), 1.2 for (q sg = 0.1 m 2 s −1 ) and 1.4 for (q sg = 10 −6 m 2 s −1 ).

Background
Studies of turbulent plumes caused by subglacial discharge and their effect on submarine glacier melting have been performed using 2-D and 3-D nonhydrostatic general circulation ocean models (GCM) (Sciascia et al., 2013;Xu et al., 2012Xu et al., , 2013;;Kimura et al., 2014;Slater et al., 2015).These models are much more complex than our simplified 1-D equations, which enable them to simulate plume processes in greater details.However, they require multi-dimensional grids with high spatial resolution, which is computationally prohibitive for our purpose of simulating a large number of Greenland glaciers.
These models typically parameterize unresolved, subgridscale turbulence with a turbulent diffusivity.Kimura et al. (2014) and Slater et al. (2015) tuned the diffusivity in such a way that the axisymmetric simulated plume (without ice contact) showed the same characteristics as the analytical models of Turner (1973) and Morton et al. (1956).Xu et al. (2013) used a high spatial resolution in order to reduce the amount subgrid processes.These models were run for idealized fjord configuration with constant subglacial discharge and a verti- The fjord is well mixed with T a = 4 • C and S a = 34.2psu, and the discharge was set to q sg = 0.1 m 2 s −1 with E 0 = 0.1.Note that the profiles of α = 90 • and α = 10 • are very similar but the cumulative melt rate of the shelf glacier increased on 500 %.For the long floating tongue the cumulative melt rate is an order of magnitude higher.
cal ice front.In most LP experiments, where subglacial discharge was uniformly distributed along the glacier grounding line, 2-D settings were chosen.The melt rate in these experiments was computed using Eqs.( 5)-( 7).Since these models are more advanced compared to the simple plume model used in this study, it is informative to compare results of our plume models with these models.

Line plume simulations
We compare the melt rate profiles obtained in the experiments by Sciascia et al. ( 2013) with the LP model.Sciascia et al. ( 2013) used a 2-D GCM with a single 10 m wide grid cell for the width and a 600 m deep and 160 km long fjord with a resolution of 10 m × 10 m.For our simulation we used the same temperature and salinity profiles as in Sciascia et al. ( 2013) and the same subglacial discharge per unit of glacier front (q sg = 0.43 m 2 s −1 ).We used an entrainment factor of E 0 = 0.08 consistent with their experiments.The vertical melt rate profile of the simulated LP model resembled that of the melt rate simulated by the 2-D GCM model but is systematically overestimated by the LP model.If we apply a scaling factor of 0.48 to the results of the LP model, the two profiles are in reasonable agreement.Still, there are some differences.The melt rate simulated by Sciascia et al.
(2013) declines with height while the LP model simulates a constant melt rate over a broad depth interval.This is due to the fact that the plume model is not applicable in the vicinity of the fjord surface.A similar effect is seen in the 2-D experiment of Xu et al. (2012) in Fig. 11a.Again, the LP model overestimates the melt rate, but when scaled up by a factor of 0.75 it yields reasonable agreement with the GCM results of Xu et al. (2012).

Cone plume simulations
For the channelized subglacial discharge the most recent numerical experiments by Slater et al. (2015) and Xu et al. (2013) were compared with simulations of the CP model.We used the same experimental settings (discharge, salinity and temperature profiles) as in the experiments of the 3-D models, with an entrainment rate E 0 = 0.1.Xu et al. (2013) used results of a survey of Store Glacier (500 m deep and 5 km wide) performed in 2010, in particular the observed temperature and salinity profile.They performed simulations of plumes for different discharge values but the same diffusivity for a 150 m wide, 500 m deep fjord with a 1 m resolution near the glacier.Their sensitivity study showed that uncertainty in channel width yielded 15 % uncertainty in the cumulative melt rate.Figure 11b shows the dependence of the cumulative melt rate on the discharge for a single plume from Xu et al. (2013) and the CP model.Both models reveal a similar dependence of melt rate on discharge, but the CP model underestimates the melt rate compared to the 3-D GCM.To bring the two melt rates in better agreement a scaling factor of 3.4 was needed for the CP model.Slater et al. (2015) used a coarser-resolution GCM with parameterized turbulence.They calibrated the GCM (vertical plume, without ice) against pure plume theory for each applied discharge value by adjusting the diffusivity until plume properties (temperature, salinity, thickness and velocity) matched plume theory by Turner with E 0 = 0.1 (Donald Slater, personal communication, 2016).Turner's plume theory is similar to our CP model (Eqs.8-11) but omits the terms with melt rate ṁ and drag C d .After tuning, the GCM was applied to simulate the melt rate for the same discharge values and diffusivity for a vertical ice front.Furthermore a minimum velocity of U 0 = 0.04 m s −1 was introduced to create a background melting the is calculated with Eqs. ( 5)-( 7).In order to simulate the same cumulative melt rate as Slater et al. (2015), distributed by 1-10 channels, a scaling factor is needed of 2.46 without calculating background melting and 1.7 with calculating background melting is needed.In contrast, the result of the GCM for the same total subglacial discharge but uniformly distributed along the whole grounding line is rather close to the results of the LP model.Indeed, for this case Slater et al. (2015) received the cumulative melt rate of 3.69 m day −1 , while for the LP model we receive 2.42 m day −1 for E 0 = 0.1 and 3.71 m day −1 for E 0 = 0.036.

Outcome of comparison
From this comparison of simple models with physically based model it appears that the LP model needs to be scaled down (except for Slater et al., 2015) and the CP model scaled = 30 m 2 s −1 and E0 = 0.07, U 0 = 3 m s −1 and the same stratification as in Xu et al. (2012).A scaling factor of 0.74 is needed to match the two melt profiles (black, dashed line).(b) Average melt rate over a 150 m wide and 500 m deep glacier part as a function of discharge localized in one channel.Following Xu et al. (2013), for the x axis, the discharge Q sg was divided by the area of the ice face A ice = 150 × 500 m 2 so that q sg = 50 m day −1 corresponds to Q sg = 43.4m 3 s −1 .The numerical results of Xu et al. (2013) are displayed with the blue line.Taking the same conditions (T a S a , Q sg ) and an entrainment factor of E 0 = 0.1 the CP model gives the solid black line.To match the experiment a scaling factor of 3.40 is needed (black, dashed line).
up.The scaling factor is of the order of 1.Most importantly, CP and LP models reveal a similar qualitative behavior to much more complex and computationally demanding GCMs as shown in Xu et al. (2012), Slater et al. (2015) and Sciascia et al. (2013).

Comparison with empirical data
Few studies exist in which submarine melt has been calculated directly based on field measurements.We used here the available data to test the LP and CP models against observations.However, the results have to be observed with caution since a single temperature profile does not necessarily represent monthly or even annual temperature profile.As Jackson et al. (2014) shows, for Sermilik Fjord and Kangerdlugssuaq Fjord in the winter months the properties, including heat content, can undergo great variability within timescales of 3 to 10 days (Jackson et al., 2014).Furthermore, uncertainties in the estimation of melt rates from fjord flux gates have been analyzed in depth by Jackson and Straneo (2016) by considering the total heat and salt budget of the fjord.Thus, when considering the derived melt fluxes from field measurements we have to keep in mind the other possible melt rate flux contributors as sea ice or melting ice bergs.

Petermann Glacier
For the years 2002-2006Rignot and Steffen (2008) calculated the melt rate of the floating tongue of Petermann Glacier obtained from ice flux divergence.They detected four large channels incised into the underside of the floating tongue.Due to its long floating tongue, the estimated melt rate is reliable because it is less affected by errors in estimating the calving rate as it is the case for tidewater glaciers.For modeling the melt rate of Petermann Glacier we used temperature and salinity profile in the fjord in front of the floating tongue measured in the year 2003 by Johnson et al. (2011).We also use the data from Morlighem et al. (2014) to define the margins of the Petermann Glacier and to compute the average one-dimensional profile of the floating tongue.We then use a polynomial fit to smooth the profile of the floating tongue.Figure 12a shows the annual mean melt rate calculated with the LP model for E 0 = 0.16 and E 0 = 0.036.Our calculated melt rates were compared to the width-averaged melt rate derived by Rignot and Steffen (2008), which is mostly dominated by the four channels that have maximal melt rates of 30 m day −1 .Even for a minimum discharge of 10 −5 m 2 s −1 (as discussed in Sect.3.1) and with E 0 = 0.036, the LP model significantly overestimates the melt rate beyond a very narrow range (few kilometers) directly next to the grounding line.This is an expected result because for long floating tongues at a certain length L, Coriolis force becomes important for small to moderate Rossby numbers R and a (horizontal) plume velocity U if U f L < 1.On this length scale the plume flow gets deflected (to the right here) and will be dominated by geostrophic flow as modeled by Gladish et al. (2012).Yet this is not taken into account in our simple plume model, as discussed in Sect.3.4.However, when using the CP model and a large discharge given by the total runoff over the catchment area distributed over four identical subglacial channels, we receive very low melt rates (Fig. 12b).It is clear that the LP model is in better agreement than the CP model at simulating the melt rate near the grounding line of the Petermann glaciers but a correction for the Coriolis effect is required further from the grounding line.

West Greenland glaciers
In a small fjord in West Greenland the melt rate of four glaciers was determined by measuring the fjord salinity, temperature and velocity close to the glacier fronts (Rignot et al., 2010).In Torsukattak fjord (TOR) the average and cumulative melt represents the melt rates of both glacier fronts together (Sermeq Avangnardleq and Sermeq Kujalleq) since the fronts are situated in the same head of the fjord branch.The two other glaciers, Kangilerngata Sermia (KANGIL) and Eqip Sermia (EQIP), enter different fjords.The measured velocity in front of EQIP does not show an upwelling pattern but more a right-to-left circulation; nevertheless we also calculated the melt rate with our plume models for EQIP.For all glaciers we took the total width of the glacier to determine the subglacial discharge per unit of length for the LP model and determined the average depth of the grounding line as a starting point for the LP model.The total subglacial discharge Q sg was taken from the table of Rignot et al. (2010) We then compare our simulations to the average melt rate determined by Rignot et al. (2010).As shown in the experiment by Slater et al. (2015), a large number of channels act like an LP but we also computed cumulative melt, assuming the existence of one large single CP starting at the maximum depth of the grounding line.Table 4 shows the ratio between observed and simulated melt rate for two types of plume models and two values of entrainment rate factor E 0 .For KANGIL and EQIP results of the LP model are in reasonable agreement with measurements, especially for the smallest E 0 value (45-105 %).Although for EQIP the agreement is the best with the LP model, the lack of upwelling circulation indicates that the plume model may not be applicable to this glacier and therefore this agreement may be a pure coincidence.The melt rate ratio of one CP shows rather poor results (1-5 %).We also compared our model with the data from Fried et al. (2015) for Kangerlussuup Sermia Glacier, which is located in West Greenland, northward of the previously discussed glaciers.Realistic temperature stratifica- tion can lead to maximal melting at the bottom of tidewater glaciers near the grounding line (e.g., Fig. 9a).This maximal melting at the bottom may cause so-called undercutting, which may enhance mass loss by calving (Rignot et al., 2015).Fried et al. (2015) found that 80 % of the tidewater glacier is undercut by 45 m in average.The glacier releases subglacial discharge via two large channels, but their corresponding melting contributes only 15 % of the total melt of glacier front.Furthermore, Carroll et al. (2016) showed that the simulated melt rate of a single cone plume is about 2 orders of magnitudes lower than the spatially averaged melt rate by Fried et al. (2015).Thus we investigate whether the LP model can calculate the average melting by assuming that the 250 m deep glacier is undercut below 50 meters depth, with an angle of 77 • to achieve the observed undercutting (Fig. 13a).Bartholomaus et al. (2016) give the belonging CTD data and estimate an summer discharge.We use the CTD closest to the glacier front in summer 2013 and the mean summer discharge (208 m 3 s −1 ) per glacier width (3 km) as input data for the LP model.Fried et al. (2015) find a total melt rate of 2 m day −1 for the whole calving front.They assumed that the glacier is only undercut by submarine melting, such that the distance of grounding line to the overhang position subtracted by the glacier's velocity gives the submarine melt value.With these input data and an entrainment rate factor of E 0 = 0.036 we achieve an average melt rate of 1.5 m day −1 (Fig. 13).This value is close to the empirical data but this plume would not result in the mentioned undercutting depth, since it penetrates up to 10 m below the sea surface.The entrainment factor E 0 = 0.16-0.13lets the plume stop at 50 m depth but their melting corresponds only to 50 % of the empirical data for the total melt rate.If the LP model is correct, additional fjord circulation makes up 50 % of the melting.Furthermore, CTD profiles close to the glacier might be diluted by near local surface runoff or calving and thus cooling and freshening of the surface ocean waters.Thus deriving the melt rate from a such a CTD profile can lead to high uncertainty ranges.tures, salinities and velocities were measured at seven stations in the fjord to calculate the melt rate of Helheim Glacier.We applied the temperature and salinity profile closest to the glacier for the LP model to simulate the melt rate for comparison.We assume, following Sutherland and Straneo (2012), that Helheim Glacier is a tidewater glacier and has a depth of 700 m and a width of 6 km and the subglacial discharge of 5.1 ± 0.76 km 3 a −1 (summer 2008; Andersen et al., 2010).Figure 14 shows the simulated melt rates over for different E 0 with the average subglacial discharge.Our best fit computes an average melt rate of 1.7 m day −1 (1.8 m day −1 ; Sutherland and Straneo, 2012) with an entrainment factor E 0 = 0.036.

Store Glacier
Another well-documented glacier is Store Glacier.Xu et al. (2013) estimated an average submarine melt rate of 3.0 ± 1.0 m day −1 in summer (Sect.4) while new calculations, thanks to new bathymetry data, reveal a melt rate of 4.5 ± 1.5 m day −1 (Chauche, 2016).Additionally, Chauche (2016) conducted a survey to determine average melt rate and subwww.the-cryosphere.net/12/301/2018/The Cryosphere, 12, 301-323, 2018 Table 5.Comparison of the melt rate calculated with the LP model and the empirical data obtained with the Gade and Motyka model (Chauche, 2016).
Melt (m day −1 ) Melt (m day −1 ) LP (Chauche, 2016) (E 0 = 0.036) Gade 2.2 ± 0.5 0.6 ± 0.1 Motyka 1.6 ± 0.4 0.7 ± 0.3 Average 1.9 ± 0.5 0.7 ± 0.2 glacial discharge from November 2012 until May 2013.Two different techniques were used, which we reference as Gade (Gade, 1979) and Motyka (Motyka et al., 2003) (see Fig. 15a).The Motyka technique is based on conservation of heat, and volume.The Gade technique is based on the identification and quantification of different processes (i.e., submarine melting, runoff mixing, thermal cooling, local sea ice formation) that can be identified by their typical temperature-salinity (TS) properties in a TS plot (i.e., Straneo et al., 2011).To simulate melt rates, we used the LP model with E 0 = 0.036 and an input subglacial discharge determined by Motyka and Gade with the corresponding temperature and salinity profiles.Results from the LP model are biased low compared to the measurements (Fig. 15b), with melt rate underestimated by 75 % in average (Table 5).Note that the Motyka method comes with large error bars for both subglacial discharge and the corresponding melt rate, which accommodate for the LP model bias (Fig. 15).Stated uncertainties for the Gade method are smaller and are not consistent with the LP model results.6.Note that only one representative value for the Store Glacier in winter was used, as reported in Table 6.208-328 1.9-3.01.0-1.21.5-1.9Store (winter, Gade model) (ST) 8-73 1.7-2.7 0.5-0.7 0.9-1.1 Store (summer) (ST su) 201-291 3.0-6.01.4-1.7 2.4-2.9

Summary
We tested both line and cone plume models against available empirical data for melt rate, and the line plume was best suited to reproduce observations (Table 4).Table 6 and Fig. 16 provide for each glacier the measured discharge and melt rate, with error bars, and corresponding range in simulated melt rate (when errors in observed discharge are taken into account as input).When default drag, heat and salt transfer coefficients are used, the simulated melt rate tends to underestimate observed melt rate, and thus the best match was obtained with an entrainment rate, E 0 = 0.036, on the lower end of our range (e.g., see Fig. 7 for how melt rate varies with E 0 ).Nevertheless, three glaciers (Helheim, Eqip Sermia and Kangerlussuup Sermia) out of seven match observations within the error bars.Varying other model parameters can change the mean but not the spread of simulated melt rate across glaciers and discharge ranges.For instance, if the heat exchange coefficient is increased to T = 4.2×10 −2 (instead of the default T = 2.2×10 −2 ), the bias can be reduced and simulated melt rates are close with observations (Helheim and Eqip Sermia fall out, conversely).Figure 16 shows a comparison of measured and simulated melt rate with the modified heat transfer coefficient.Clearly, many fjord processes are not taken into account in this simplified approach.For example, the circulation in front of Eqip Sermia was mostly horizontal (Rignot et al., 2010) instead of the vertical upwelling represented in the model.There are also issues with the measurements themselves, such as time sampling or difficulties to retrieve discharge and melt rate, as seen for the Store Glacier (Fig. 15) or for Helheim, where CTD profiles for temperature and salinity were taken 1 year after discharge rates measurements.Nevertheless, the simple line plume model is in general agreement with the observations (Fig. 16) and shows a correlation coefficient of 0.7 (Fig. 16b) with the modified heat transfer coefficient.The theoretical background and similar dependency on discharge compared to more complex models (see previous sections) make it suitable for modeling studies over a larger number of Greenland glaciers, and to investigate melt rate response to future changes in subglacial discharge and fjord temperature.

Conclusions
We presented two simple models for simulation of the submarine melt rate of marine-terminating glaciers, the so-called cone plume and line plume models, and studied sensitivity of these two models to different forcings (fjord temperature, stratification, subglacial discharge) and model parameters (entrainment parameter E 0 ).We also compared these models with results of experiments performed with 2-D and 3-D ocean GCM by Slater et al. (2015) and Xu et al. (2013).Lastly we compared the results of simulations of the LP and CP models with empirical estimates of melt rate for several Greenland glaciers.Our analysis demonstrates that for small subglacial discharge, typical for winter conditions, cumulative melt does not depend on the discharge.For high discharge typical of summer conditions we found a power dependence of 1/3 of submarine melt on subglacial discharge for the LP models, which is consistent with previous studies.We found a theoretical explanation of this behavior, explained in Appendix A. Furthermore, we found that the power dependence to the ambient temperature in a well-mixed environment is 1.7-1.8 for lower discharges and is only 1.2 for the higher discharge for both models.
We investigated the sensitivity of the melt rate to the entrainment parameter E 0 that was used to parameterize turbulent entrainment into the plume.For a tidewater glacier the cumulative melt rate of the LP model increases with decreasing E 0 while it decreases for the CP model.This is explained by the fact that although in both cases higher entrainment rate www.the-cryosphere.net/12/301/2018/The Cryosphere, 12, 301-323, 2018 slows down the plume and reduces the melt rate per unit of area, for the CP, this effect is overcompensated by the widening of plume for the higher entrainment coefficient.In general, we found a notable effect of the entrainment parameter on the melt rate for the range of entrainment parameter given in the literature.The uncertainty range of E 0 can have the same effect as 1 • C change in ocean temperature.
Our comparison of the CP and LP models to results of 3-D GCM experiments showed qualitatively similar melt rate profiles.
In most the LP model overestimates the results of the GCM by approximately a factor of 2, while the CP model underestimates melt rate from GCMs.More importantly, we find the same power law dependence of melt rate on subglacial water discharge as in Slater et al. (2016), for given ambient hydrographic conditions.As a result, with a constant scaling factor of the order of 1, the simplified models can reproduce a wide range of melt rates spanning several orders of magnitude.
In the case of the long floating tongue, like the Petermann Glacier, the LP model significantly overestimates the melt rate outside of the narrow zone along the grounding line, which is probably due to the missing Coriolis force in the plume models.
Although it is known that in summer a part of the subglacial meltwater is delivered in the fjord through several channels, we found that the submarine melt rate associated with the discharge through the channels and better described by the CP model makes up only a small amount of the empirically estimated total melt rate of a glacier front.Furthermore the total number of channels for every summer is unknown for different glaciers.When we compare the LP model to empirical data, it is evident that the LP model is more appropriate than the CP model for simulation of both winter and summer melt of real Greenland glaciers.However, the model has to be adjusted for individual glaciers since the scaling parameter is not the same for different glaciers.Thus, for the future we will use the tuned LP model coupled to a 1-D ice flow model to determine the importance of submarine melt rate to glacier dynamics.
In this Appendix, we analyze the LP model equations in order to derive approximate analytical solutions.This in turn helps to interpret the results of the numerical experiments presented in this paper, performed with the more complete plume models from Jenkins (2011).First analytical solutions for the LP model were undertaken by Linden et al. (1990) and summarized in Straneo and Cenedese (2015).Slater et al. (2016) previously presented approximate analytical solutions for the LP model for higher discharge ranges.Jenkins (2011) noticed that for strong discharge, plume velocity in the LP model does not change much with depth and is thus similar to the initial balance velocity (our Eq. 12).Magorrian and Wells (2016) covered the case for small discharge.The reasoning in this Appendix provides a unifying solution for small and large discharge with the LP model applied at tidewater glacier and glacier with long floating tongues.
We restrict the analysis to the typical conditions of a 500 m deep Greenlandic fjord (T a ; 0-4 • C).

A1 Simplified melt rate equation
After Jenkins (2011), the melt rate can be approximated as where T = T − T f is the temperature above freezing and M 0 is a slowly varying function of ice temperature below freezing point, which can be considered constant for the purpose of this Appendix.Freezing point temperature is given by T f = λ 1 S +λ 2 +λ 3 Z.We run several experiments in a typical parameter range for tidewater and long floating tongue glaciers in Greenland's fjords and could confirm that the approximation is accurate for the LP model (Fig. A1a).The parameters were T a (0-4 • C), q sg (1 × 10 −5 -0.1), E 0 (0.036-0.16) and sin α (0.02-1) for constant depth of 500 m, T i = −10 • C and S a = 34.2.With linear regression we found an average value for M 0 = 8.8 × 10 −6 .Let T e = E 0 M 0 sin α, the entrainment-equivalent temperature ( • C), be a measure of the ratio of entrainment to melting (it corresponds to the temperature for which melting equates entrainment).We have ) in all experiments (Fig. A1b), consistent with the ranges for E 0 (0.036-0.16) and sin α (0.02-1), so that T e spans 2 orders of magnitude, roughly 10 2 -10 4 • C.

A2 Balance regime
In Fig. 2 we showed that CP velocity rapidly converges regardless of initial velocity.This has been also shown by Dallaston et al. (2015) for the LP and also holds for the plume temperature, salinity and melt rate.Here we derive analytical solutions for these convergence values (indicated with ) and associated length scales for the our approximation of the LP model (i.e., Eqs.A1 and A2), by using the equation for the volume flux (Eq. 1) so that where q = DU (the volume flux) and X can be either T , S or U .The convergence value X can be obtained by solving the corresponding equation (qX) = f (where f is the righthand-side term, e.g., Eqs. 2, 3 or 4) with X = 0.Moreover, when the right-hand side term is not or weakly dependent on X (i.e., for T and S, as will be detailed below), the equation is analogous to a first-order differential linear equation with convergence length scale L X = q q ≈ q ė = D E 0 sin α , i.e., with fast convergence near the grounding line, where plume thickness D is small.

A2.1 Balance temperature and salinity
Temperature and salinity Eqs. ( 3) and ( 4) can be rewritten as an intuitive mixing law by merging in Eqs. ( 5) and ( 6): where T m is an effective meltwater temperature, derived from Eq. ( 5): Variations of boundary layer temperature T b around 0 • C can be safely neglected compared to latent heat, so that we will treat T m as a constant.If T i = −10 • C, we have T m ≈ −89 • C. Nevertheless, for completeness note that T b can be expressed as a function of melt rate, plume and ice temperatures from Eq. ( 5).Using our simplified melt rate Eq. ( A1) and given that ṁ C 1 2 d T U by 2 orders of magnitude, an accurate approximation for T b is given by where we verify that boundary layer temperature is somewhat closer to freezing temperature than to plume temperature.In the case of plume salinity S b cancels out completely and S i = 0 (as can be verified straightforwardly using Eqs. 4 and 6), so no other term is needed.By decomposing Eq. (A4) as outlined in Eq. (A3), and searching for solutions when T = 0, with ṁ ė, we obtain an expression for balance temperature: which can be rearranged by using Eq.(A2) in "balance" regime, so that ṁ ė ≈ ṁ ė ≈ T T e , and neglecting the secondorder T a /T e , into The ratio −T m /T e spans about 10 −2 to 1 in our experiments, with a maximum of 0.02 for tidewater glaciers, i.e., T ≈ T a .Here the freezing temperature implied by should be taken for balance plume salinity, which is nearly the same as ambient salinity in first approximation (Eqs.A5, A2, A10): ≈ S a , (A11) so that T ≈ T − T fa and T a ≈ T a − T fa , where T fa is the freezing temperature for ambient salinity.

A2.2 Balance velocity
A similar reasoning as in the previous section (using Eq.A3 and q ≈ ė), Eqs. ( 1) and ( 2) can be rearranged into an equation for U 2 (note the identity (U 2 ) = 2U U ): where b = sin(α)g ρ ρ 0 and C e = E 0 sin α.This highlights in one equation basic plume dynamics, buoyancy-accelerated and balanced by drag and entrainment.Equation ( A12) is analogous to a first-order linear differential equation with equilibrium solution for x>>L u : and length scale Note that Eq. (A13) does not represent a strict equilibrium but a dynamic balance between velocity, plume thickness and buoyancy, which is maintained while the plume thickness and associated volume flux keeps increasing.Note that as the plume dimension increases due to entrainment (or for large discharge), so does the length scale L u , and U may lag behind U .At x = 0 for typical discharge and entrainment values L u is less than a centimeter.Our simulations show that velocity reaches dynamic balance U within the first few meters after the grounding line (not shown).This shows qualitative agreement with the above analysis but suggests that initial changes in plume dimension D and buoyancy b should be taken into account for more detailed analysis of the transient regime.In the present analysis we focus on the balance regime.The theoretical equilibration length scale for velocity is shorter than for temperature and salinity by a factor of 2 or more, since L T S /L U = 2(1 + C d E 0 sin α ), especially for long floating tongues.In the actual simulations the ratio is even larger, because the plume keeps growing with distance from its source.
The balance velocity is more conveniently expressed as a function of plume's volume flux q instead of thickness D. By taking the square of Eq. (A13) and multiplying by U , the identity q ≈ q = DU can be used, so that . (A15) This new expression shows that velocity can be written as a power law of the buoyancy flux qb.Its initial condition is known, since q = q sg and b = b 0 at x = 0.The evolution of qb as the plume develops can also conveniently be derived We used the fact that if U = (A + Bx) 1/2 , where A and B are constant (Eq.A19), a primitive U is 2 3B (A + Bx) 3/2 = 2 3B U 3 , provided that B = 0, i.e., that U = 0; additionally since B = (U 2 ) is proportional to 2 3 M 0 T (Eq.A17), these terms cancel out.
The error of Eq. (A20) compared to the cumulative melt rate of the LP model in the unstratified case for tidewater glaciers was 2 % for large discharge (q = 0.1 m 2 s −1 ) and 9 % for small discharge(q = 1 • 10 − 6 m 2 s −1 ).For the case of a long floating tongue and a discharge of q = 0.1 m 2 s −1 the error was in the range of 10 %.
We can identify two limiting cases for Eq.(A20), characterized by discharge-driven (U ≈ 0, i.e., large discharge M high ) or melt-driven (U 0 ≈ 0, i.e., small discharge M low ) convection (see also Eqs.A18 and A19).Let us define L as the full length of the plume (note that it is related to grounding line depth Z = L sin α).We have discharge-driven, q sg q c , (A21) , melt-driven, q sg q c , (A22) where the critical discharge q c is defined as the discharge for which M high = M low , equal to For a tidewater glacier q c is of the order of 10 −3 m 3 s −1 , and can reach 10 −1 m 3 s −1 for a few tens of kilometers of floating tongue.This indicates that in glaciers with a floating tongue, melt-driven convection tends to be more important than in tidewater glaciers.We also note from Eqs. (A21) and (A22) that the dependence of melt rate on entrainment coefficient C e (which contains both entrainment parameter E 0 and glacier slope sin α) is not monotonic.For moderately large values of C e (−M 0 T m /C e 1, e.g., tidewater glaciers), variations in the left-hand-side factor ("dynamics") dominate and the melt rate is a decreasing function of, but moderately sensitive to, C e .However, for small C e values (−M 0 T m /C e ∼ 1 and thus C e C d , i.e., long floating tongue), where variations in the right-hand-side factor ("thermodynamics") dominate, the melt rate is an increasing function of C e .For the lowest entrainment value E 0 = 0.036, "small" means a glacier front angle of the order of sin α = 0.02.

Figure 1 .
Figure1.(a) Conceptual scheme of the 1-D plume model afterJenkins (2011).Subglacial freshwater flux q sg , which is uniformly distributed along the grounding line, enters the fjord, where it forms a plume that rises up due to buoyancy.The plume is described explicitly with its temperature T , salinity S, thickness D and velocity U .It rises along the ice shelf, following the shelf's slope α = 90 • − β, until it either reaches the water surface or has zero velocity due to the loss of buoyancy.The ambient water with salinity S a and temperature T a entrains into the plume with entrainment rate ė.Melting ( ṁ) occurs on the glacier front and adds to the plume buoyancy with water of the temperature T b and salinity S b .(b) Conceptual scheme of two-dimensional CP model modified afterJenkins (1991) andCowton et al. (2015).Subglacial discharge enters the fjord, localized via a channel.The plume geometry is described as a half cone and the entrainment occurs around the arc.The subglacial

Figure 2 .
Figure 2. Different runs of the CP model for different initial velocities.Panel (a) depicts the velocity profile in the first 100 m.All initial velocities converge within 100 m to the trajectory of the balance velocity U 0 = 3.5 m s −1 (thick black line).The corresponding initial radii differ thus from 300 m (for U 0 = 3.5 × 10 −3 m s −1 ) to 3 m (for U 0 = 35 m s −1 ).Panel (b) shows the corresponding melt profile.Higher initial velocities give a maximal melt rates at deeper levels.All melt rate profiles converge to the same melt rate after a certain depth.

Figure 3 .
Figure3.Sensitivity of cumulative melt rate to different initial velocities for both plume models.Melt rate (black) is in percent of the cumulative melt achieved with initial balance velocity U 0 (red).Red dashed line shows 120 % mark.Only very high initial velocities can appreciably increase the cumulative melt rate for the CP.

Figure 4 .
Figure 4. Melt rate profiles in a well-mixed fjord simulated by the CP model (black) and LP model (blue) for a width W = 150 m and the total discharge of Q sg = 500 m 3 s −1 .In the case of the CP model the total discharge is delivered through one channel in the center of the glacier; in the case of LP model the discharge is uniformly distributed with the rate q sg = Q sg W = 3.6 m 2 s −1 .Both plumes start with a velocity of U 0 = 1 m s −1 .Solid lines show melt rate averaged across the plume in the case of CP model and across the full width W in the case of LP model.The dashed line shows the corresponding cumulative melt rate for the entire glacier.

Figure 5 .
Figure 5. Cumulative melt rate of the LP model for E 0 = 0.1 with different discharge values for a well-mixed environment with T a = 4 • C and S a = 34.65 psu.The gray dashed line is our analytical solution for the cumulative melt rate of the LP model (Eq.A20).The red dashed and blue dashed lines indicate the limiting melt rate regimes M high (Eq.A21) and M low (Eq.A22), respectively.The transition of these regimes is defined by the critical discharge q c (Eq. A23), depicted here with black dashed line.

Figure 6 .Figure 7 .
Figure 6.Cumulative melt rates of the different plume models as a function of entrainment parameter E 0 for four different discharge values.The cumulative melt rate is depicted for (a) LP of a tidewater glacier (sin(α) = 1), (b) LP of a long floating tongue and (c) CP of a tidewater glacier.For the LP model for sin(α) = 1 a higher E 0 leads to lower cumulative melting opposite to the other two cases.

Figure 8 .
Figure8.Cumulative melt rate per glacier width W for the LP model (black) and CP (blue) model as a function of the thermal forcing ( T a = T a − T af ) for high (solid lines) and low (dashed lines) discharge values.The experiment is for a well-mixed, 500 m deep and 5 km wide tidewater glacier (sin(α) = 1), with S a = 34 psu and E 0 = 0.1.

Figure 9 .
Figure9.Influence of stratification and discharge on the melt rate profile of the LP (a).The three different discharge values (q sg = 0.5, 0.1, 10 −6 m 2 s −1 , dashed, solid, dotted) in a stratified environment for a fixed salinity profile (d) and five different temperature profiles (c) result in 15 different melt rate profiles.The melt rate of the corresponding temperature profile is displayed in the same color as well as in the same style (dashed, dotted or solid) for the corresponding discharge.Note that a very high discharge (q sg = 0.5 m 2 s −1 ) is needed for the plume to reach the surface.For each discharge value the corresponding cumulative melt rate is depicted (b) as a function of the thermal forcing ( T a = T a − T b , Eq. 7) at the grounding line.For ṁ ∼ T

Figure 10 .
Figure 10.Melt profile (a) and corresponding plume velocity profiles (b) plume temperature (c) and salinity (d) for the LP model for different glacier types: a tidewater glacier (α = 90 • ), shelf glacier α = 10 • and a shelf glacier with a long floating tongue (α = 1.1 • ) of 25 km.The fjord is well mixed with T a = 4 • C and S a = 34.2psu, and the discharge was set to q sg = 0.1 m 2 s −1 with E 0 = 0.1.Note that the profiles of α = 90 • and α = 10 • are very similar but the cumulative melt rate of the shelf glacier increased on 500 %.For the long floating tongue the cumulative melt rate is an order of magnitude higher.

Figure 11 .
Figure 11.Comparison between LP, CP and GCM simulations.(a) Experimental results from Xu et al. (2012) (blue line) and LP model (black, solid line) for Q sg = 150 5= 30 m 2 s −1 and E0 = 0.07, U 0 = 3 m s −1 and the same stratification as inXu et al. (2012).A scaling factor of 0.74 is needed to match the two melt profiles (black, dashed line).(b) Average melt rate over a 150 m wide and 500 m deep glacier part as a function of discharge localized in one channel.FollowingXu et al. (2013), for the x axis, the discharge Q sg was divided by the area of the ice face A ice = 150 × 500 m 2 so that q sg = 50 m day −1 corresponds to Q sg = 43.4m 3 s −1 .The numerical results ofXu et al. (2013) are displayed with the blue line.Taking the same conditions (T a S a , Q sg ) and an entrainment factor of E 0 = 0.1 the CP model gives the solid black line.To match the experiment a scaling factor of 3.40 is needed (black, dashed line).

Figure 12 .
Figure12.Melt rate (a) LP model over the long floating tongue of Petermann Glacier with a minimal discharge of Q min = 10 −5 m 2 s −1 for the minimal (E 0 = 0.036) and maximal (E 0 = 0.16) value (black lines) of the entrainment parameter.In panel (b) we used the maximum discharge Q sg = 296 m 3 s −1 (total runoff assumed only in summer) distributed over four channels to compute the melt rate with the CP model.As forcing variables we used the fjord's temperature and salinity profile in front of the floating tongue for the year 2003 summarized byJohnson et al. (2011), and fromMorlighem et al. (2014) we determined the glacier thickness and depth of the floating tongue (see Sect. 5.1 for details).For both E 0 the melt rate is highly overestimated with the LP model and underestimated with the CP model.The empirical melt rate estimated byRignot and Steffen (2008) is displayed with the blue line.Note the different vertical scale on the panels.
Sutherland and Straneo (2012) used results of a field campaign in Sermilik fjord in summer 2009 in which tempera-

Figure 13 .
Figure 13.Sermia average undercut profile at the terminus (a) with the assumed temperature profile (b) give different melt rate profiles (c) simulated for different E 0 and same Q sg = 208 m 3 s −1 .All average melt rates are below the determined 2 m day −1 by Fried et al. (2015).

Figure 14 .
Figure 14.vertical melt rate profiles of the LP model (a) for three different entrainment coefficients E 0 for Helheim Glacier.With a discharge of 2.69 × 10 −2 m 2 s −1 and E 0 = 0.036 we obtain an average melt rate of m = 1.7 m a −1 close to Sutherland and Straneo (2012) (1.8 m day −1 ).

Figure 15 Figure 16 .
Figure 15.(a) Estimated discharge of Store Glacier for winter 2012/2013 from(Chauche, 2016).Red ranges give subglacial discharge estimates by the Motyka method and blue ranges by the Gade method and (b) the corresponding melt rate profiles.Simulated melt rates by the LP model with E 0 = 0.036 are depicted as the red dotted line (subglacial discharge from Motyka) and blue dotted line (subglacial discharge from Gade model).

Figure A1 .
Figure A1.Investigation of melting proportion in the plume equations for different LP experiments.The plume model was run in a well-mixed environment for different parameter settings:E 0 [0.036-0.16],sin(α)[0.02-1],q sg [10 −5 -0.1 m 2 s −1 ], T a ] [0-4 • C. Panel (a)shows the melt rate as a function of plume velocity U and plume temperature T and its freezing temperature T f ( ṁ = M 0 (T − T f )U ).The average slope of the run is M 0 ≈ 8.75 × 10 −6 (thick black line) while for the simplified equations of Jenkins (2011) a typical slope lies of the order of 6.5 × 10 −6 (black dashed line).The second panel illustrates that ṁ ė in this parameter range, but ṁ/ ė being largest for long floating tongues.

Table A1 .
Summary of Appendix variables.Illustrative value provided for T i = −10 • C, T a = 4 • C, S a = 34.65 psu, sin α = 1 (tidewater glacier) and range for E 0 = 0.036 − 0sin α(β S S a − β T T a ) buoyancy at x = 0 0.27 m s −2 b m g sin α(β S S a − β T (T a − T m )) buoyancy source term due to melting 0a − T f (S a ) ambient temperature above freezing ≈ 6 • C

Table 1 .
Model parameters of the LP and CP model with typical fjord default values.Note that the values may differ for specific experiments (as explicitly stated in the corresponding descriptions).
(Xu et al., 2013)charges with the CP model and LP model.By assuming a melt rate equation of ṁ = a • Q β + b as in(Xu et al., 2013), we derived β numerically for both models and listed them in(Table ).The CP model shows values close toXu et al.

Table 2 .
Determination of the power law β of melt rate in the equation

Table 4 .
Simulated cumulative melt rate (%) of empirical estimated cumulative melt rate for different entrainment rates for three West Greenland glaciers.

Table 6 .
Estimated subglacial discharge Q emp and melt rate m emp for a number of glaciers and corresponding melt rate m lp from LP model simulations.Values of Store winter are taken from the Gade model (Fig.15).For each glacier, local hydrography (temperature and salinity profiles) and measured subglacial discharge is used to drive the LP model.Ranges indicate measurement errors.Errors in subglacial discharge are propagated to errors in simulated melt rate via the LP model.For the Store Glacier in winter we only report a representative value in the table, where both mean and error were averaged.Simulated melt rate m LP is obtained with E 0 = 0.036, T = 2.2 × 10 −2 .Melt rate m LP with modified T = 4.2 × 10 −2 is also provided.