Submarine melt parameterization for a Greenland glacial system model

Two hundreds 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 plays 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 and outlet glaciers will retreat, contributing to sea level rise. Projections of the future contribution of outlet glaciers to sea level rise is hampered by the necessity to use extremely high resolution of the 5 order of a few hundred meters both for modelling of the outlet glaciers and as well as 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 based on using of computationally efficient parameterization of submarine melt based on turbulent plume theory. We show that such parameterization is in a reasonable agreement with several available modeling studies. We performed a suit 10 of experiments to analyse sensitivity of these parameterizations 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 parameterization of a line plume can produce submarine melt compatible with observational data. Our results show that the line plume model is 15 more appropriate than the cone plume model for simulating the 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 ± 8 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 related, 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 kilometers (Bindschadler et al., 2013), which is too coarse to resolve most of Greenland outlet glaciers.Instead, regional modelling at higher resolution is better suited to capture glacier dynamics.As an alternative to costly three-dimensional models, one-dimensional flowline models were convincingly applied to several major outlet glaciers (Nick et al., 2012(Nick et al., , 2013;;Lea et al., 2014;Carr et al., 2015).In particular, Nick et al. (2012) simulated with a flowline model the dynamical response of the Petermann glacier to the abrupt break up 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.Nick et al. (2013) using the same flowline model 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 with the parameterization of the rate of submarine melt.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 derived to calculate submarine melt rates of outlet glaciers by using empirical data (Motyka et al., 2013;Rignot et al., 2015); simplified one dimensional models of line plumes (Jenkins, 1991(Jenkins, , 2011)), axisymmetric plume models (Cowton et al., 2015;Turner, 1973) and numerical three-dimensional ocean models (3D models) (Holland et al., 2008a;Little et al., 2009;Sciascia et al., 2013;Xu et al., 2013;Slater et al., 2015).Note that 3D modelling efforts also differ with respect to model formulation, in particular some authors use non-hydrostatic models (e.g.Holland et al., 2008b;Little et al., 2009), while others use hydrostatic models (e.g.Sciascia et al., 2013).The experiments studied submarine melt with respect to subglacial discharge and its spatial pattern, vertical ocean temperature and salinity profiles.Additionally the influence of the fjord circulation, which connects outlet glaciers with the surrounding ocean, were investigated with the 3D models.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 is localized (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 observations 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).These simulations, in agreement with previous theoretical studies, show 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.On the other hand, Holland et al. (2008b) and Little et al. (2009) found a quadratic dependence on temperature under large ice shelves, where geostrophic flow becomes significant.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 and 0.85 at low discharge for the CP.Simulations with 3D models, which differ with respect to boundary conditions and turbulence parameters, show a variety of CP melt rate profiles: Kimura et al. (2014) showed a melt rate profile of the CP that reaches its maximum near the water surface while Slater et al. (2015) and Xu et al. (2013) found a CP melt rate profile with the maximum located near to the bottom.
While experiments with high-resolution (several to ten meters) nonhydrostatic 3D 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 at centennial time scale.An alternative is to use a parameterization of submarine melt based on a simplified plume model (Jenkins, 2011;Cowton et al., 2015).
Such parameterization can then be used to calculate submarime melt in a 1D 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 parameterizations to simulation of melt rate of real glaciers in Greenland.To this end we first compared both cone and linear plume parameterizations with the available results of simulations from high resolution 3D ocean models.Then we compare results of plume parameterizations 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 the section 2. There we study the plume models sensitivity of simulated submarine melt rate to ocean temperature and salinity, the amount of subglacial discharge and to the ice tongue geometry of the glacier itself.Results of simulations with the simple plume parameterizations are compared to results of numerical experiments with 3D ocean models in section 4. In section 5 we compare our simulations to empirically estimated submarine melt rates for several selected Greenland glaciers.Finally, in the section 6 we discuss the applicability of the plume parameterization 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, which are (a) additional melting under the floating tongue (if any) and along the glacier front, and (b) turbulent entrainment and mixing of surrounding fjord water.They 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 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 line plume (LP) model corresponds to the simplest assumption that Q sg is uniformly distributed along the grounding line (Fig. 1), while the cone plume (CP) assumes point-wise release of meltwater (Fig. 1b), i.e. from a channel whose dimensions are small compared to the plume diameter.Note there need not be only one CP: they can be a number of them discretely distributed along the glacier.

Model equations
Both models are formulated in one dimension, x, which is the distance along the glacier front from the grounding line and depends on the glacier shape, described by its slope α.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 description 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. 1).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 = () : where the plume state variables D, U , T and S stand for its thickness, velocity in the x-direction, temperature and salinity, all dependent on x.Equation (1) describes the conservation of volume flux q = DU (expressed per unit length in 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 density contrast ∆ρ 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 (subscribt b) at the ice-water interface, where melting occurs, and heat and salt is exchanged by (turbulent) conduction-diffusion.The Equations for T and S (Eq.3,4) account for the entrainment of ambient water and the addition of meltwater, as well as for conduction fluxes at the ice-water interface (i.e. between boundary layer and plume).The entrainment rate is calculated as ė = E 0 U sin(α), proportional to plume velocity and glacier slope, with coefficient E 0 .The melt rate is calculated by solving for heat and salt conservation at the ice-water boundary ( ṁ, T b and S b are unknown): 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 equations (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 , where T − T f is the plume temperature above freezing point, and M 0 is a slowly varying function of ice temperature below freezing point.Numerically, M 0 varies from 2.9 • 10 −6 to 0.910 −9 ( • C) −1 over a T i − T f range from −20 • C to 0 • C, respectively, and the freezing temperature is roughly −2 • C (Annex A).We do not use this approximation in our calculation, but this is nevertheless helpful to interpret some of the results presented in our manuscript, in particular the dependence of the melt rate on plume velocity (Annex 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 the 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.The equations ( 1)-(4) now reform for the CP model by considering melting on the diameter 2D and entrainment around the arc πD: where variables, parameters and equations have the same meaning as for the LP model, and the volume flux Q = πD 2 2 U is expressed in cubic meters per second.

Numerics
For the differential equation system of (1)-( 4) and ( 8)-( 11) we choose a classical Runge-Kutta-scheme in which we can regulate the regular grid size ∆x.Thus we have control of the numerical calculation time and can easily vary T a (Z) and S a (Z) for a stratified environment.Furthermore for glaciers with floating tongues sin(α) can vary as a function of X and therefore the model can adjust dynamically to the glacier in a coupled glacier-plume version.With these initial conditions for the plume T, S, U, D at x = 0 we solve the equations ( 5)-( 7)and firstly determine the melt rate ṁ| x0 and the boundary conditions S b | x0 and T b | x0 .These determined variables serve as the input parameters for the differential euqation system to determine the plume properties at the next step X i+1 = X i + ∆x.This routine is continued to determine the melt rate as a function of x 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 manuscript, 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 melt water, the salinity and temperature of subglacial discharge water can be set to zero (S 0 = 0 and T 0 = 0).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 the LP case, the subglacial discharge per glacier width W enters the model equations: q sg = Qsg W = U 0 D 0 (Fig. 1).
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 equations (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 linear differential equation for U 2 (see Appendix section 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 equations ( 1) and (8) (see Appendix section A1).We obtain for the LP: and for the CP: Note that equation ( 12) is identical to the velocity derived by Jenkins (2011), and equation ( 13) is analogous to equation ( 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).More detailed discussion and 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 near the bottom of the glacier ('undercutting') while for lower velocities the melt rate increases with height and would leave a so-called 'toe' at the glacier bottom (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 ten 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.

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 Greenlands 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 Greenland fjords, most of them do not have a floating tongue (tidewater glaciers) and we therrfore generally perform experiments for a vertical wall (sin(α) = 1).Default model parameters, including entrainment rate E 0 , are indicated in Table 1.
A direct comparison between LP (defined per unit length) 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 ).In the case of the LP model we assumed that the discharge is uniformly distributed over a fjord width W = 150m, so that q sg = 3.6m 2 /s.This width is about the maximum size of the plume in the CP model, near the surface.

Comparison between LP and CP models
Results in Figure 4 show that simulated, local melt rate is higher in the CP model than in the LP model practically for all depths, but cumulative melt rate (i.e.integral of the meltrate 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 two in our chosen setting).
We shall see later in this manuscript (sec.3.1) that the (local) melt rate in the LP model varies less than linearly with subglacial 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 are more than one channel is discussed in section 4.3.

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 (for the high discharge range) and for the CP in a well-mixed environment.Note that for the LP this dependence can already be determined by the look on the balance velocities U 0 (Eq.12 ).However for smaller discharge in a well-mixed environment, cumulative melt rate converges to a small but not insignificant value that does not obey the power law any more (Fig. 5.This value represents background melt rate which does not depend on discharge and can be representative for winter melt rate when subglacial discharge is very small.To explain this change of power law we undertook a dimensional analysis to obtain theoretical solutions for the plume model (A).Important is that, in a well mixed fjord, the melt rate is linear dependent on the velocity of the plume.This velocity is dependent on discharge and rather constant over the glacier front for big discharge or independent on sublgacial discharge and will accelerate along the x direction for small discharges (Annex A3, Fig. A2).Therefore one can separate the plume behaviour into regimes for high and low discharges.We also derived an analytical solution for the cumulative melt rate (Eg.A20) that is also displayed in Figure 5.
However, as mentioned in the introduction, stratification can change this power law.We also performed experiments with stratification as in (Xu et al., 2013) for different discharges with the CP model and LP model.The CP model shows values close to Xu et al. (2013).Both models show an increasing exponent for lower discharge (Tab.2).

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 available for melting, but on the other hand to decreased buoyancy -and velocity -as the plume gets saltier.Slower velocity as a result negatively affects melting (Eq.5, 6) (Carroll et al., 2016) (note the plume also becomes thicker to accomodate for increased volume flux and decreased velocity).In this section we investigate what is the net effect of these processes on melting for typical plume configurations.
In both plume models, entrainment depends on an entrainment rate parameter E 0 (sec.2.1), which is not accurately known and can be regarded as a tunable parameter within a certain range of values known from pevious work.Laboratory experiments for a pure vertical plume and model studies gives for E 0 a broad range from 0.036 to 0.16 (Jenkins, 2011;McConnochie and Kerr, 2016;Kaye and Linden, 2004;Mugford and Dowdeswell, 2011;Carroll et al., 2015b;Kimura et al., 2014).
It is known that the cumulative melt for tidewater glaciers decreases with increasing E 0 (Jenkins, 2011;Magorrian and Wells, 2016) (in the reported range of E 0 values).The melting is controlled in first order by the plume velocity and only to a lesser extent by availability of heat through mixing for stong subglacial discharge (Magorrian and Wells, 2016).Moreover the relative effect of entrainment can dominate the effect of the ambient temperature (Fig. ( 7)).E.g.A LP with low entrainment E 0 in a colder ambient temperature of T a = 3 • C will result in higher cumulative melt rate (due to its higher velocity) than for a plume with high E 0 in a warmer fjord of T a = 4 • C.
For the CP model (Fig. 6c) the dependence on E 0 is opposite: cumulative melt rate increases with the entrainment factor.This is due to the faster growing plume radius with higher entrainment, which leads to a larger area of the plume in contact with the ice, and thus to more melting overall, despite the lower local melt rate.

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 glaciers response to global warming.Previous experiments with 3D ocean models showed different 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 non-linear dependence of the melt rate on water temperature.If the discharge is very small, ambient properties dominate the melting process and one can speak of a 'melt driven convection' (Slater et al., 2015).
If we assume a power law dependence of the cumulative melt rate per glacier area to the thermal forcing, i.e. m ∝ T F β , where T F = T a − T af and T af is the freezing temperature of the sea water at the fjord bottom, we find that the exponent β increases with lower discharge.From 1.2 (high discharge q = 0.1 m 2 s −1 ) to 1.8 q = 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 increment compares well to analytical solutions while the CP model seems not to show this change in power law for analytical solutions (table 3).
An increment of the power law has also been detected when we use 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. (2015a) showed that plume thoery gives a good approximation of the outflow height for 3d nonhydrostatic plume model but nevertheless do the exponents differ slightly in the experiment by Xu et al. (2013) (table 3).

Glacier front angle
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 (10 c).The dependence of melt to the slope of the glacier has been derived by Magorrian and Wells (2016) but note that the choice of E 0 can have a similar effect as the choice of subglacial discharge as depicted in Figure 10.
A glacier with a long floating tongue, and therefore a smaller angle (i.e.sin(α) = 0.02), has a smaller average melt rate than a tidewater glacier.However in this case higher E 0 leads to higher cumulative melting (Fig. 6 b).These high cumulative melt rates (Fig. 10) occur due to the longer distance under a floating tongue in which the velocity accelerates (10 b).The theroetical explanation of the evolution of U , ṁ, T is explained in the Annex (A and summarized in A3).However, for small α both plume models are not applicable along the total shelf because they do not take into account Coriolis force and therefore likely strongly overestimate plume velocity and melt rate (see more in section 5.1).

Background
Studies of turbulent plumes caused by subglacial discharge and their effect on submarine glacier melting have been performed using 2D and 3D non hydrostatic general circulation ocean models (GCM) (Sciascia et al., 2013;Xu et al., 2012Xu et al., , 2013;;Kimura et al., 2014;Slater et al., 2015).Although these models contain the right physics to simulate plume dynamics, the problem is that it requires very high spatial resolution which is computationally too expensive for our purpose.In order not to resolve the small-scale turbulences, a parameterization for turbulent diffusivity is chosen to represent subgrid-scale mixing.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 resolve turbulence explicitly.These models were run for idealized fjord configuration with constant subglacial discharge and a vertical ice front.In most LP experiments, where subglacial discharge was uniformly distributed along the glacier grounding line, 2D settings were chosen.The melt rate in these experiments was computed using equations (5 -7).Since these models are more advanced compared to simple plume parameterization used in this study, it is informative to compare results of plume parameterization 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 2D GCM with a single 10 m wide grid cell for the width and a 600m deep and 160 km long with a resolution of 10 m × 10 m.For this 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 2D 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 2D experiment of Xu et al. (2012) in figure 11 a).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 (and most in agreement with plume theory) 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 3D models, with an entrainment rate E 0 = 0.1.Xu et al. (2013) used results of a survey to Store Glacier (500m deep and 5km wide) performed in 2010, in particular the observed temperature and salinity profile.They performed simulations of plumes for different discharge values but same diffusivity for a 150 m wide, 500 m deep fjord with a 1m resolution near the glacier.Their sensitivity study showed that uncertainty in channel width yielded 15% uncertainty in the cumulative melt rate.Fig. 11 b) 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 3D GCM.To bring the two melt rates in better agreement a scaling factor of 3.4 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 (personal communication from D. Slater).Turners plume theory is similar to our CP model (eq.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.04ms −1 was intoroduced to create a backround melting the is calclutated with Equation (5-7).In order to simulate the the same cumulative malt rate as Slater et al. (2015) distributed by 1 − 10 channels a scaling factor is needed of 2.46 without a backround melting and 1.7 with background melting is needed.On the other hand, 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 md −1 , while for the LP model we receive 2.42 md −1 for E 0 = 0.1 and 3.71 md −1 for E 0 = 0.036.

Conclusions
From these comparison of simple parameterizations 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 up.The scaling factor is in the order of one.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); Sciascia et al. (2013).

Comparision with empirical data
Few studies exist where 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 the average heat storage in a fjord since it can undergoe great variability (Jackson et al., 2014).

Petermann glacier
For the the years 2002-2006Rignot and Steffen (2008) calculated the melt rate of the floating tongue of Petermann glacier obtained from velocity measurements and mass balance.They detected 4 big channels of along 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. (2011b).We also use the data from Morlighem et al. (2014) to define the margins of the Petermann glacier and to compute average one-dimensional profile of the floating tongue.We then use a polynomial fit to smooth the profile of the floating tongue.Fig. 12 a) shows the annual mean melt rate calculated with the LP model for E 0 = 0.08 and E 0 = 0.036.Even for a minimum discharge of 10 −4 m 2 s −1 (as discussed in section 3.1) and with E 0 = 0.036, the LP model significantly overestimates the melt rate beyond a very narrow range (few km) 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.At this length scale the plume flow gets deflected (to the right here) and will be domintad by geostrophic flow as modelled by Gladish et al. (2012).Yet this is not taken into account in our simple plume model, as discussed in section 3.4.On the other hand, 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 correction for 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 salinity, temperature and velocity in the fjord near the glacier fronts (Rignot et al., 2010).In Torrsukatak fjord (TOR) the average and cumulative melt represents the melt rates of both glacier fronts together (Seermeq Avangnardleq and Sermeq Kujatdleq) since the fronts are situated in the same head of the fjord branch.The two other glaciers, Kangilerngata Sermia (KANGIL) and Equip Sermia (EQUIP) enter different fjords.Measured velocity in front of EQUIP 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 EQUIP.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.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 acts like a LP but we also computed cumulative melt assuming the existence of one big 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 EQUIP results of the LP model are in reasonable agreement with measurements, especially for the smallest E 0 value (45%-105%).Although for EQUIP the agreement is the best with the LP model, the lack of upwelling circulation indicates that the plume parameterization 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 previously discussed glaciers.We used data from Morlighem et al. (2014) for the glacier elevation and after averaging to a one dimensional profile we obtained a shelf of 3 km length.Note that caution is needed since the data set is averaged over 10 years and has a resolution of 300m.Realistic temperature stratification 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 meters in average.The glacier releases subglacial discharge via two big channels, but their corresponding melting contributes only 15% of the total melt of the glacier front.Furthermore Carroll et al. (2016) showed that the simulated melt rate of a single cone plume is about 2 magnitudes lower than the spatially averged melt rate by Fried et al. (2015).Thus we investigate whether the LP model can calculate the average melting by assuming that the 250 meter 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 closet to the glacier front in Summer 2013 and the mean summer discharge (208m 3 s −1 ) per glacier width (3km) as input data for the LP model.Fried et al. (2015) find a total melt rate of 2 md −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 this input data and an entrainment rate factor of E 0 = 0.036 we achieve an average melt rate of 1.5 md −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 meter below the sea surface.The entrainment factor E 0 = 0.16 − 0.13 lets 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 the LP model is correct it means, that additional fjord circulation make out 50% of the melting.Furthermore deriving the melt rate by one CTD profile close to the glacier might be diluted by near local surface runoff or calving an thus cooling and freshening of the surface ocean waters.

Helheim
Sutherland and Straneo (2012) used results of a field campaign in Semerlik fjord in summer 2009 where temperatures, 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 (section 7) 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.76km 3 a −1 (summer in 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.6 md −1 (Sutherland 1.8 md −1 ) 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 md −1 in summer (sec.4) while new calculations, due to new bathymetry data reveal a melt rate of 4.5 ± 1.5 md −1 (Chauche, 2016).
Additionally, Chauche (2016) conducted a survey to determine average melt rate and subglacial 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) in Figure 15 a).The Motyka technique is based on conservation of heat, salt 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 modify the water coloumn by it's tempreature and salinity in certain -for the process -typical gradients.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, to simulate melt rates.Results from the LP model are biased low compared to the measurements (Figure 15 b), 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 corresponding melt rate, which accomodate for the LP model bias (Figure 15).Stated uncertainties for the Gade method are smaller and are not consistent with the LP model results.

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 Figure 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, 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, Equip 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 the bias can be reduced and simulated melt rates are close with obervations (Helheim and Equip 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 Equip 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 one year after discharge rates measurements.
Nevertheless, the simple line plume model is in general agreement with the observations (Fig. 16) -and shows a a correlation coefficient of 0.7 ( 16 b) 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 modelling studies over a larger number of Greenland glaciers, and to investigate melt rate response to future changes in discharge and temperature.
6 Conclusions 1) We presented two parameterizations for simulation of the submarine melt rate of marine-terminated 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 parameter (entrainment parameter E 0 ).We also compared these models with results of experiments performed with 2D and 3D ocean GCM by Slater et al. (2016) and Xu et al. (2013).At last we compared the results of simulations of the LP and CP models with empirical estimates of melt rate for several Greenland glaciers.
2) We found that for small subglacial discharge, typical for winter conditions, cumulative melt does not depend on the discharge.For high discharge typical for summer conditions we found a power dependence of 1/3 of submarine melt on subglacial discharge for the LP models, and a power of 2/5 for the CP model, which is consistent with the previous studies.We found a theoretical explenation of this behaviour, explained in the Annex A. Furthermore we found that the power dependence to the ambient temperature in a well-mixed environment also is 1.7-1.8 for lower discharges and is only 1.2 for the higher discharge for both models.
3) We investigated the sensitivity of the melt rate to the entrainment parameter E 0 that was used parametrization of the turbulence of 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 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.4) We compared the CP and LP models to results of 3D GCM experiments, and find qualitatively similar melt rate profiles.
In most cases, the LP model overestimates the results of the GCM by approximately a factor two, while the CP model underestimates melt rate from GCMs.Such discrepancy is not surprising given the highly simplified parameterization of the LP and CP models compared to GCMs.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 one, the simplified models can reproduce a wide range of melt rates spanning several orders of magnitude.5) 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.6) 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 out 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 futue we will use the tuned LP model coupled to a 1d ice flow model to determine the importance of submarine melt rate to glacier dynamics.
By decomposing (A4) as outlined in (A3), and searching for solutions when T = 0, with ṁ ė, we obtain an expression for balance temperature: which can be rearranged by using (A2), and neglecting the second order T a /T e , into: The ratio −T m /T e spans about 10 −2 to 1 in our experiments.Here the the freezing temperature implied by ∆ should be taken for balance plume salinity, which is nearly the same as ambient salinity in first approximation (Eq.( A5), (A2), (A10)): where T f a is the freezing temperature for ambient salinity.

A2.2 Balance velocity
A similar reasoning as in the previous section (using (A3) and q ≈ ė), (1) and ( 2) can be rearranged into an equation for U 2 (note the identity (U 2 ) = 2U U ): where b = sin(α)g∆ρ 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 and length scale Note that equation (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 with entrainment, so does the length scale L u , and the feedback balance becomes looser.At x = 0 for typical discharge and entrainment values L u is less than a centimeter.Our simulations show that velocity reaches dynamic The error of (A20) compared to the cumulative melt rate of the LP model in the unstratified case for tidewater glaciers was 2 % for big discharge (q = 0.1m 2 s −1 ) and 9 % for small discharge(q = 1 • 10−6m 2 s −1 ).For the case of a long floating tongue and a discharge of q = 0.1m 2 s −1 the error was in the range of 10 %.
A3 The role of subglacial discharge and the shape of the glacier We investigated the plume properties and melt rate of a typical tidewater glacier and a glacier with a long floating tongue (order of Petermann glaicer).While for the tidewater glacier the plume temperature rapidly approaches the temperature of the ambient water (Fig. A3 c) the plume under a long floating tongue stays cooler since the melt-entrainment ratio becomes bigger (A8).A look on the velocity of the plume shows an acceleration under floating tongues.Equation (A19) reveals that for a tidewater glacier, a plume starting wit a velocity U 0 , which is dependent on the subglacial discharge, will accelerate with a slope independent of the subglacial discharge.Therefore plumes with small discharges will highly accelerate while the velocity of plumes with big discharge will remain almost constant along Z (Fig. A2) as reported by Dallaston et al. (2015).That explains the different exponents of melt rate as a function of subglacial discharge in the literature.In the case of a very small discharge q sg → 0 then U → 2 3 bm C d +Ce M0∆Ta 1−Tm/Te x and the melting becomes independent of the discharge and we speak of the background melting.Magorrian and Wells (2016) undertook a scaling analysis for the plume model for small discharges.For q sg → 0 comparision shows that in our analysis velocity accelerates with U (x) ∼ sin(α)X equal to the results of Magorrian and Wells (2016).The cumulative melt rate for small discharges in dependence to the thermal forcing ∆T a gives a power law of 3/2, also confirmed by the work of Magorrian and Wells (2016).For tidewater glaciers with very high discharge the acceleration term (A17) can be neglected and therefore the velocity -and thus the melt rate-depends on the initial conditions of subglacial discharge with the cubic root.Also the melt rate is then linear dependent on ∆T a .Our approximation of U , T and m = M 0 • ∆T • U are displayed along the LP models results in Figure A3 and show good agreements.The approximation of the cumulative melt rate (A20) shows the biggest deviation for the floating tongue with 10 %. 10 -10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 Q sg [m 2 /s]   high (q = 0.1 m 2 s ) 1.2 a 1 b,c,d LP, well mixed low (q = 1 • 10 −6 m 2 s ) 1.8 a 1.5 c,d LP, well mixed high (q = 0.1 m 2 s ) 1.2 a 2 (linear stratified) c LP, realistic stratified low (q = 1 • 10 −6 m 2 s ) 1.4 a 2(linear stratified) c LP, realistic stratified high (q = 500 m 3 s ) 1.2 a 1 e CP, well mixed low (q = 5 • 10 −3 m 2 s ) 1.5 a 1 e CP, well mixed high (q = 500 m 3 s ) 1.1 a -CP, realistic stratified low (q = 5 • 10 −3 m 2 s ) 1.3 a -CP, realistic stratified high (q > 5.8 m 3 s ) 1.2 f -3D CP, realistic stratified low (q < 5.8 m 2 s ) 1.6 f -3D CP, realistic stratified   q sg = 0.1 U 0 =0.64 q sg = 1e-06 U 0 =0.01

Figure 1 Figure 2 .Figure 3 .Figure 4 .
Figure 1.a) Conceptual scheme of 1D plume model afterJenkins (2011).Uniformly distributed along the grounding line subglacial freshwater flux qsg, enters the fjord and 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 .The plume rises along the ice shelf an slope α = 90 • − β until it either reaches the water surface or has zero velocity due to the loss its buoyancy.The ambient water with salinity Sa and temperature Ta entrains into the plume with an 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 and the entrainment occurs around the arc.The subglacial discharge isQsg = D 2 0 U 0 π 2where D0 is the initial radius and U0 is the initial velocity

Figure 5 .Figure 6 .Figure 7 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .Figure 14 .Figure 15
Figure 5. Cumulative melt rate of the LP model (a) for E0 = 0.1 with different discharge values for a well-mixed environment with Ta = 4 • C and Sa = 34.65 psu.For red line corresponds ṁ = 6.6 • 10 −5 • Q 1 3sg + 6.0 • 10 −6 for the higher discharge range.The grey dashed line (a) is our analytical solution for the cumulative melt rate of the LP model (Eq.A20) Determination of the power law β of melt rate on discharge Q or the equation ṁ = a(b•Q β +c).Seperation between high (Q > Qc) and low discharge (Q < Qc) Determination of the power law β of melt rate on thermal forcing T F = Ta − T af or the relation m ∝ T F β for low and high discharge with comparision to theoretical estimated values.Upper eight rows are from this study lower show results of other 3d simulaitons.

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: E0[0.036 − 0.16], sin(α)[0.02− 1], qsq[10 −5 − 0.1 m 2 s ], Ta[0 − 4] • C. Panel a) shows the melt rate as a function of plume velocity U and plume temperature T and it's freezing temperature T f ( ṁ = M0(T −T f )U ).The second panel illustrates that ṁ << ė in this parameter range,but beeing biggest for long floating tongues.

Figure A2 .Figure A3 .
Figure A2.Evolution of U for an initial velocity of U .The plume with a small discharge (qsg = 10 −6 m 2 s ) will accelerate strongly (red line,U = 0.14 m s ) while the plume velocity with larger discharge remains almost constant (black line, qsg = 0.1 m 2 s , ,U = 0.63 m s ).

Table 1 .
Model parameters of the LP and CP model with typical fjord default values.

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