A numerical model for meltwater channel evolution in glaciers

Meltwater channels form an integral part of the hydrological system of a glacier. Better understanding of how meltwater channels develop and evolve is required to fully comprehend supraglacial and englacial meltwater drainage. Incision of supraglacial stream channels and subsequent roof closure by ice deformation has been proposed in recent literature as a possible englacial conduit formation process. Field evidence for supraglacial stream incision has been found in Svalbard and Nepal. In Iceland, where volcanic activity provides meltwater with temperatures above 0 °C, rapid enlargement of supraglacial channels has been observed. Supraglacial channels provide meltwater through englacial passages to the subglacial hydrological systems of big ice sheets, which in turn affects ice sheet motion and their contribution to eustatic sea level change. By coupling, for the first time, a numerical ice dynamic model to a hydraulic model which includes heat transfer, we investigate the evolution of meltwater channels and their incision behaviour. We present results for different, constant meltwater fluxes, different channel slopes, different meltwater temperatures, different melt rate distributions in the channel as well as temporal variations in meltwater flux. The key parameters governing incision rate and depth are channel slope, meltwater temperature loss to the ice and meltwater flux. Channel width and geometry are controlled by melt rate distribution along the channel wall. Calculated Nusselt numbers suggest that turbulent mixing is the main heat transfer mechanism in the meltwater channels studied.


Introduction
Flow of water through glaciers has received considerable attention from the scientific community since theoretical treatment of the phenomena began with two publications in 1972 (Röthlisberger, 1972;Shreve, 1972).Recently, two detailed review articles have summarized the current state of knowledge, one focusing on jökulhlaups (Björnsson, 2010), also known as glacial outburst floods, and the other focusing on Röthlisberger channels (Walder, 2010), which are waterfilled, pressurized englacial channels.The formation and evolution of surface meltwater channels has been studied in the field (e.g.Knighton, 1981;Marston, 1983) as well as in the laboratory (e.g.Isenko et al., 2005) and treated analytically, as evolving from surface crevasses (e.g.Fountain and Walder, 1998) or forming during drainage of surface lakes (e.g.Raymond and Nolan, 2000).Supraglacial channels evolving into englacial conduits have been considered as a mechanism for the formation of englacial passages (e.g.Fountain and Walder, 1998;Benn et al., 2009;Gulley et al., 2009).The role of supraglacial drainage systems during englacial tuya eruptions is still not fully understood, but they are thought to play an important role as a controlling mechanism for eruption site lake levels (e.g.Smellie, 2006).
Simulating temporal and spatial evolution of meltwater channels in ice, which is key to understand the processes involved, requires adequate numerical models capable of resolving the underlying physics in great detail.Subglacial drainage systems have been studied extensively, both numerically and theoretically (e.g.Röthlisberger, 1972;Shreve, 1972;Nye, 1976;Cutler, 1998;Ng, 1998;Flowers and Clarke, 2002;Clarke, 2003;Creyts and Schoof, 2009;Schoof, 2010) bed has been identified to cause acceleration of ice-sheet flow in Greenland (e.g.Zwally et al., 2002;Bartholomew et al., 2010) which in turn contributes significantly to eustatic sea level rise (Rignot and Kanagaratnam, 2006).New models of supraglacial channel evolution are required to quantify and further understand surface meltwater transport to the bed of large ice sheets.In combination with advanced models of subglacial hydrology (e.g.Schoof, 2010) they hold great promise for detailed modelling of the process chain from surface meltwater availability to eustatic sea level contribution of ice sheets.
In this contribution we present a new model that for the first time provides explicit numerical simulation of supraglacial meltwater channels with a focus on their evolution as well as their incision behaviour.We introduce the physical basis of our model in Sect.2, followed by notes on the numerical implementation in Sect.3. Subsequently, we present model results for different key model parameters in Sect. 4 and close with conclusions and an outlook on future research.

Model physics
Our model consists of three components: (1) ice dynamics, (2) turbulent meltwater flow in open channels, and (3) thermal transfer between meltwater and ice.Figure 1 displays the principal geometry of our model setup and introduces several model parameters.
Along the contact area between meltwater in an open channel and ice, which forms the channel geometry, heat exchange plays a key role in the evolution of the system.Within the scope of this paper, we focus on a specific interaction along this boundary: outward growth of the channel geometry driven by thermal melting of the ice walls, counteracted by creep closure of the same channel.
To successfully model the evolution of such a system, all three aforementioned components have to be combined, which mathematically form a Stefan problem (Lamé and Clapeyron, 1831;Stefan, 1891).We will first introduce the physics we use for each component and subsequently describe the numerical details in Sect.3.

Ice dynamics
Ice is simulated as a Stokes fluid with a non-linear viscous behaviour.Starting with the Stokes equations for which u denotes the velocity vector, p the pressure field, ρ ice the ice density, and g the gravity vector, we further include the standard rheology of ice (Glen, 1955;Nye, 1957) in our physical representation, which leads to a non-linear ice viscosity such that: (3) Here we write the Glen rate factor as A, the Glen exponent as n and the effective strain rate as ˙ = 1/2 ˙ ij ˙ j i .This physical description of ice is commonly referred to as a "Full Stokes" model.

Open channel hydraulics
To simulate the turbulent flow of water in an open channel, we use the Gauckler-Manning formula (Gauckler, 1867;Manning, 1891), which relates cross-sectional, average velocity (V ) in an open channel to its slope (β) and hydraulic radius (R h ).Although this is an empirical formula, it can be derived analytically from the phenomenological theory of turbulence (Gioia and Bombardelli, 2002).Based on the Gauckler-Manning formula, the water flux (Q) in an open channel can be expressed as Here n c is the Gauckler-Manning coefficient (cf.Table 1) and the hydraulic radius can be written as R h = A c /P , with A c the cross-sectional area of the flow and P the wetted perimeter.

Water -ice thermal transfer
Following the approach of Raymond and Nolan (2000), we can relate the melting at the channel walls to the energy loss in the water flow:  (2010, p. 75).n c for ice channels is taken from Fountain and Walder (1998).
Symbol Value with L the latent heat of fusion per unit mass, ρ w the water density, and C w the heat capacity of water per unit mass.dθ/dl is the change in water temperature dθ per unit length dl along the channel.To simplify Eq. ( 5), we introduce in accordance with Raymond and Nolan ( 2000) and assume dθ/dl to be the input from the water temperature above freezing θ divided by the distance l along the channel over which the temperature drops to freezing.Based on Eq. ( 5), a mean melt rate m along the channel perimeter is computed as: We assume that the local melt velocity at any point along the wetted perimeter depends on depth below the water surface in channel, where H max is the instantaneous elevation of the water surface in the channel, and s is arc length measured along the ice surface boundary.For a linear scaling of local melt rate with water depth, ν = 1, but it is not a priori clear if this is the 'correct' choice.Other values for ν are plausible.A uniform m loc along the channel perimeter leads to ν = 0. Assuming a 'shallow channel' of Ng (1998) would result in non-linear scaling with ν = 3/2.
3 Numerical implementation

Power law fluid
The power law fluid described by Eqs.(1-3) can be solved for complex geometries using the finite element method (FEM, e.g., Zienkiewicz et al., 2005).We use icetools (Jarosch, 2008) in its revised version1 to solve the ice dynamics for our transient complex geometries.To achieve convergence of the non-linear viscosity problem we substitute successively u and η in a classical Picard iteration.Further details on the numerical implementation, e.g.how we handle Eq. (3) in case of vanishing strain rate, can be found in Jarosch (2008).

Ice geometry and boundaries
We embed the evolving meltwater channel into a rectangular, two dimensional cross-section (along the x and z coordinate) of an idealized glacier.The point of origin for our model coordinate system is at the base in the centre of the glacier cross-section (cf.Figs. 1 and 2).Any forming meltwater channel is located at x = 0 m and at the initial glacier surface, where a predefined small surface depression facilitates the initial location of the meltwater stream (cf.Fig. 3a).
In the model it is possible for the glacier surface to slope perpendicularly towards the channel by setting the slope angle α > 0 (cf.Fig. 1).We set u = 0 at the lateral ice boundaries as well as at the basal boundary.To avoid any influence from these static boundaries on the evolving channel, we place the lateral boundaries at x = ±1900 m and set the initial ice surface at z = 500 m.The ice surface, including the meltwater channel wall, is modelled as a stress free boundary (for which η ∇u + (∇u) T • n + pn = 0, with the unit normal n) and we use a forward Euler scheme on the same grid to evolve the surface for each time step (cf. Jarosch, 2008).For the sake of simplicity and to purely focus on the meltwater channel evolution, we disregard surface mass balance for all model results in this paper.Considering the placement of our coordinate system, we could impose a symmetry boundary at x = 0 m to simplify the model and only simulate one half-space of the domain.As a test for the numerical stability of our model, we choose to do otherwise and model the full domain for which we only impose an initially symmetric geometry.As we will later demonstrate (cf.Fig. 4), the transient simulations remain almost perfectly symmetric over long time spans.
The interface between meltwater and ice (red line in Figs. 1 and 2a) demands closer attention, as there are several physical processes (cf.Sects.2.2 and 2.3) involved in the temporal evolution of this boundary.In the next section we describe in detail how the interaction between ice dynamics, channel hydraulics and wall melting is implemented numerically.

Moving channel boundary
At each time step of the evolving model, we have to account for all three aforementioned processes (cf.Sects.2.1-2.3).If we were to model the turbulent water flow in the channel  explicitly within the FEM model, we could use an iterative scheme commonly used in fluid structure interaction models (e.g.Bungartz, 2010) to simulate the moving boundary.Here we choose a much simpler approach as we model ice dynamics with Eqs.(1-3), meltwater hydraulics with Eq. ( 4), and heat transfer from meltwater to ice with Eqs. ( 7) and ( 8).The free surface above and below the water line (H max ) evolves as where n is the outward-pointing unit normal to the surface, and x is a (non-material) point on the boundary.In Fig. 2 we illustrate individual sub-steps of moving the boundary forward in time by using Eq. ( 9).Discretizing Eq. ( 9) with a forward Euler scheme requires a time step, t, which can be chosen to be either fixed at a value small enough to avoid numerical instabilities (i.e.t < 3 days) or by utilizing a CFL condition (e.g.Courant et al., 1928) on the forward moving grid, which calculates an optimal value for t.As the surface grid points move at each time step, we also regenerate the numerical mesh for the FEM model using the updated geometry.Figure 3 displays the evolution of a typical channel in our model.In this example, a predefined flux of Q = 1 m 3 s −1 is used to fill an initial, small surface depression and subsequently the model evolves over 20 days ( t = 2 days) to form a meltwater channel.Note that the channel width changes over time (Fig. 3a-c) and the channel width after 20 days (Fig. 3c) is a result of the model physics and not prescribed by us.

Channel pinch-off and flow regime transition
After a certain time (t p ) of model evolution, the channel pinches off from the surface due to ice deformation (cf.Table 2 for values of t p and Fig. 3d-i).At this stage, the numerical mesh along the ice walls in the channel merges above the channel bottom (cf.Fig. 4).The closure velocity of the channel walls at that location decreases to zero as the mesh re-merges, simulating a direct contact and instant refreezing of the channel walls.Numerically, a control algorithm monitors the grid points at the channel walls and detects overlaps, for which the numerical mesh is merged at these locations.
As time progresses, ice flow continues to act on the channel geometry and the channel's cross-sectional area decreases until it is sufficiently small to fill the whole channel with meltwater (cf.Fig 6c-f).At this timestep (t f ), Eq. ( 4) is not valid anymore, channel hydraulics switch to a pressurized channel flow, and the model stops.The treatment of pressurized channel flow is not readily possible with our current model and we discuss extensions to include pressurized channel flow in Sect. 5. Generally, in pressurized channels with circular cross-sections, radial melting is balanced by radial inward creep of ice, thus no significant downward motion is to be expected (e.g.Röthlisberger, 1972).Including gravity and density differences within the channel water in a model would lead to a presumably small downward motion of the channel.

Results
In this section we investigate the behavior of our model for different model parameters.Obviously there exists a large model parameter space, which could be explored, but we will focus here on a small sub-set, which is intended to demonstrate the general behaviour of the system.For our sensitivity analysis of final channel incision depth (d f ) and pinch-off depth (d p ) we vary key model parameters over typical ranges (cf.Sect.4.2).Data from a meltwater channel formed after the 1996 Gjálp volcanic eruption in Iceland (Gudmundsson et al., 2004) is included for comparison.
For all results presented here, we set several model parameters to fixed values, which are listed in Table 1.Moreover, we keep the initial ice geometry with its boundary conditions predefined (cf.Sect.3.2 and Fig. 3a) except the ice Table 2. Tabular form of results displayed in Figs. 5 and 7. Final channel incision depth, d f , is given for model time t f , at which the channel switches from an open channel flow regime to a pressurized channel, while at t p the channel pinch-off occurs.Error between melt-only and full model, d f , is given in meters and relative percentage.Mean channel width w c is averaged between 5 m below the surface and depth at t p .N u = 0.0078P r 1/3 Re 0.927 are Nusselt numbers with the respective Prandtl number P r (Lunardini et al., 1986)  surface slope angle α.The meltwater temperature within the channel is assumed to be zero for most results, which leads to dθ/dl = 0 • C m −1 and therefore γ = 0.An increase of meltwater temperature provides more energy to melt ice, but does not change the general behaviour of the model.The ability to model meltwater temperatures significantly above zero is important for glaciovolcanological applications.Note that the contribution of either channel slope (β) or meltwater temperature gradient factor (γ ) to cross-sectional area change (dA c /dt) is linear in Eq. ( 5), and thus it is sufficient to vary one for a demonstration of model behaviour.For example, increasing the channel slope (β) by 0.03 is equivalent to T ≈ 0.35 • C over a channel distance l = 5000 m, or dθ/dl = 7 × 10 −5 • C m −1 and thus γ = 0.0301.We define our reference parameter set as Q = 1 m 3 s −1 , α = 0.0, β = 0.03, γ = 0.0, and ν = 1.0 which is marked as a bold red line in Figs. 5 and 7 as well as with bold font in Table 2.

Constant meltwater flux
In Fig. 5, we present results for Q = 0.1 and 1 m 3 s −1 , glacier surface slopes α = 0.0, 1.0, and 2.0, channel slopes β = 0.01, 0.03, 0.06, and 0.09, melt rate distribution factor ν = 0.01, 0.5, 1.0, and 1.5 and no meltwater temperature loss along the channel (γ = 0).Most simulations are computed until t f is reached.Results marked with (*) in Fig. 5 are stopped after 2000 days, as their incision rate is very slow.In one case, in which the channel cross-section had a large width to depth ratio, the channel closes by roof deformation before the channel gets pressurized (marked with (**) in Fig. 5).Numerical values for all runs are listed in Table 2.
Obviously different combination of model parameters lead to different incision depths as well as rates.A common feature of all runs is the decrease in incision rate towards the end of each run, shortly before t f is reached.Final channel incision depth, d f , varies more between different runs than pinch-off depth d p (cf.Table 2 and Fig. 8).As a first sensitivity test, we vary ν between 0.01 and 1.5 (Q = 1 m 3 s −1 , α = 0.0, β = 0.03, and γ = 0.0).The possible range of incision results is marked as a grey area in Fig. 5. ν = 0 would represent a uniform distribution of melt rate along the channel wall.This case is numerically challenging, as it creates an outward melting of the channel at the water line, which results in an undercutting of the channel wall.Numerically more stable is the case of ν = 0.01, which leads to a melt rate distribution very close to the case of ν = 0, but with the distinct advantage of m loc = 0 exactly at the water line.
One can anticipate, based on Eq. ( 8), that ν does not only control final incision depth, d f , but also channel shape.In Fig. 6a we demonstrate the influence of ν on channel shape after 30 days of evolution.Lower values of ν create wider channels with slower incision rate.The three channels with ν = 0.5, 1.0, and 1.5 do not vary much in width (1.12 m, 0.9 m, and 0.74 m) but the channel with ν close to zero is significantly wider (3.52 m).Sensitivity of d f on ν is shown later (cf.Sect.4.2).Differences in evolution of meltwater channels due to different values of ν are visualized in an Fig. 5. Results of incision depth evolution for meltwater fluxes Q = 0.1 and 1.0 m 3 s −1 , different glacier surface slopes α, channel slopes β, water temperature loss parameter γ , and melt rate distribution factor ν are plotted with different colors.Potential incision behaviour due to variations of ν between 0.01 and 1.5 is indicated as a grey area.All results are based on our standard parameter set (Q = 1 m 3 s −1 , α = 0.0, β = 0.03, γ = 0.0, and ν = 1.0) with modifications indicated in the legend.The timestep (t p ) at which the channel pinches off is marked with a colored dot for each model.At the end of each plotted line, the respective channel becomes pressurized (t f ), except for runs marked with (*), which are stopped after 2000 days due to their slow incision rate and the run marked with (**), where the channel closes by roof deformation before the channel gets pressurized.animation (cf.supplementary material).Channel width averaged over the whole channel evolution, w c , are listed in Table 2 for all model runs.In case of fixed meltwater flux, ν is controlling channel width, but similar widths can be achieved by significantly increasing meltwater flux.A channel with Q = 1 m 3 s −1 and ν = 0.01 creates the same width as a channel with Q = 100 m 3 s −1 and ν = 1.0.Note that channel shape, flat tip vs. round tip, is solely controlled by ν (cf.Fig. 6a), which should allow for a conclusive identification of parameters causing a distinct channel width and geometry combination.Supraglacial channel width to depth ratios measured in the field range from 3.4 to 12.0 (Knighton, 1981), which correspond to values for ν between 0.5 and  7) and (8) as a dashed line.Data from the post 1996 Gjálp eruption meltwater channel (Gudmundsson et al., 2004) is plotted as a dashed, black line with corresponding error estimates as grey areas (detailed view in insert).At the end of each plotted line, the respective channel becomes pressurized (t f ), except for the run marked with (*), which was stopped after 50 days, the length of available Gjálp data.
0.01 in our simulation (cf.Fig. 6a).Thus smaller values of ν (≤ 0.5) appear to be appropriate to model supraglacial channels as they are consistent with measured channel crosssections.
We initiate each meltwater channel with a given surface depression (cf.Fig. 3a) and so the question arises if the channel shape and evolution is dependent on this initial geometry.In Fig. 6b, different initial surface depression shapes (half ellipses) are displayed and their corresponding channels after 30 days of evolution.All presented cases converge towards the same channel shape, which is a result of Eqs. ( 7) and ( 8) and not predefined in the model.Note that initial differences in incision rate do not play a significant role on the long timescales of t p and t f .Different channel geometries at t f are displayed in Figs.6c-f, which are all close to circular.This supports our assumption that the channels get pressurized at t f and no significant downward motion is to be expected2 .
Figure 7 presents model results for higher meltwater fluxes (Q = 10.0 and 100.0 m 3 s −1 ), includes meltwater temperature loss along the channel (γ > 0), and compares our results with meltwater channel data from the 1996 Gjálp eruption, Iceland (based on Fig. 8 and Table 3 in Gudmundsson et al., 2004).Again we observe incision rate decrease as the channels reach t f and less variation of d p in comparison with d f over all results presented.
To demonstrate that transient modelling is important to correctly simulate the system, i.e. accounting for ice flow and channel melt at each time step, we compare incision behaviour of the fully coupled model with a melt-only simulation.For Q = 10 m 3 s −1 (α = 0.0, β = 0.03 γ = 0.0, and ν = 1.0), the result of our model is plotted as a solid green line in Fig. 7 and the incision behaviour based solely on Eqs. ( 7) and (8) as a dashed green line.After 514 days of evolution (t f reached) the melt-only model overestimated d f by d f ∼ 50 m, which is an error of ∼ 29 %.Overestima-tion values for all of our simulations are given in Table 2 and range between 18 to 47 %.A melt-only model has many drawbacks besides this overestimation.It can not estimate d p , t p , d f , and t f , nor can it predict correct channel shapes at any time as it lacks the influence of ice dynamics.
During the 1996 Gjálp subglacial eruption, meltwater draining from the eruption site created a large channel in the surrounding ice cap.We use information presented in Gudmundsson et al. (2004) to infer dθ/dl ∼ 0.001 • C m −1 based on the observed heat loss of 40 % from the initially 20 • C meltwater over an 8 km long subglacial path.Meltwater fluxes in the supraglacial channel ranged from Q = 10 to 100 m 3 s −1 , β was about 0.01 and γ ∼ 0.4304 (dθ/dl ∼ 0.001).Observed incision rates are plotted as black, dashed lines in Fig. 7 together with their respective error estimates as grey shaded areas.It is not the scope of this paper to simulate and reconstruct channel formation at Gjálp.However, to demonstrated the potential of our model we have included one run which reproduces the Gjálp channel incision within error range of the observed data.With Q = 10 m 3 s −1 , α = 0.0, β = 0.01, γ = 0.2152, and ν = 1.0, our model computes a similar incision behaviour as observed at Gjálp (cf.cyan line in Fig. 7).

Sensitivity analysis
Results from our sensitivity analysis on pinch-off depth d p and final incision depth d f are plotted in Fig. 8.We use our reference run (Q = 1 m 3 s −1 , α = 0.0, β = 0.03, γ = 0.0, and ν = 1.0) to infer sensitivities of d p (dashed lines in Fig. 8) and d f (solid lines in Fig. 8) by varying one parameter at a time.
An increase in α reduces d p due to an increase in cross channel ice flow.Increasing Q results in an non-linear increase of d p as more meltwater flux generally leads to faster incision (cf.Fig. 7 and Table 2).β and ν do not influence d p significantly.
In the case of d f , α displays the same, almost linear decreasing effect as found for d p .The non-linear Q sensitivity of d f is similar as seen for d p but more pronounced, because more time is available for the interplay between melt and ice flow.We find a strong, linearly increasing effect of β on d f .Note that increasing Q from 1 to 100 m 3 s −1 leads to the same increase in d f as increasing β from 0.03 to 0.09.ν is found to have the strongest control on d f .Small values of ν result in shallow d f .The maximum of d f is reached for ν = 0.5, followed by a decrease of d f with increasing ν and a subsequent small increase of d f with further increase of ν (cf.Fig. 8).It is important to realize that this complex behaviour of the model is partly due to the fact that, in the case of ν = 0.01, the channel closes by roof deformation before the channel gets pressurized.

Variable meltwater flux
To simulate an idealized meltwater flux cycle throughout a year, we vary Q over time such that which results in a period of half a year in which Q increases to Q max and decreases again to zero as well as half a year period with no meltwater flux.This is intended to represent more realistic meltwater runoff conditions on glaciers.We use a similar set of model parameters as in Fig. 5 for Q = 1 m 3 s −1 , but now we set Q max = 1 m 3 s −1 and use Eq. ( 10) to vary Q with time.The influence of different values for α on incision depth is also investigated.Results of these simulations are displayed in Fig. 9. Similar results as for constant meltwater fluxes can be observed.Again a twofold increase in β leads to a deeper incision (∼17 m more at t = 182.5 days) and increased ice flow towards the channel (α > 0) counteracts the channel evolution.The pronounced difference to the constant meltwater flux case is the period of no incision during which Q = 0 m 3 s −1 .In the case of increased ice flow towards the channel, the channel bottom even rises again during this period and in the case of α = 1.0, leads to channel closure by roof deformation due to ice flow (cf.Fig. 9).

Conclusions and outlook
In this paper we have presented a new model that for the first time provides explicit numerical simulation of meltwater channel evolution in glaciers, based on the combination of ice dynamics, open channel hydraulics, and ice-water thermal transfer.The model is capable of simulating channel incision over time for a given meltwater flux, meltwater temperature, channel slope, melt rate distribution in the channel, and initial ice geometry.Shape and evolution of the channel are purely driven by model physics and are not pre-defined.
To demonstrate the principal model behaviour, we have computed solutions for a set of model parameters.In case of constant meltwater fluxes, we have identified channel slope β, meltwater flux Q, and meltwater temperature loss to ice γ to be the main controlling parameters for channel incision depth and rate.Increased ice flow towards the channel (α > 0) counteracts channel incision.Melt rate distribution along the channel walls (ν) has a complex influence on incision behaviour and controls channel geometry, mainly channel tip shape.Comparison with melt-only model results of maximum incision depth clearly demonstrates the importance of resolving the transient nature of the system.Calculated Nusselt numbers of order 10 4 suggest that turbulent mixing is the main heat transfer mechanism in the studied meltwater channels.All solutions displayed here use γ = 0. Model time t = 182.5 days, after which Q(t) = 0 m 3 s −1 for half a year, is marked with the vertical dashed line.In case of α = 1.0 (cyan line), the the channel closes by roof deformation at the end of the plotted line due to ice flow.In the case of β = 0.06 (magenta line), we plot results to t = 220 days, highlighting the different model behaviour for the first half of a year.All other results are manually stopped after 1.5 years.
We have also computed results for Q varying over a synthetic annual cycle, creating a half year period of increase and subsequent decrease followed by a half year period of no water flux.Similar model behaviour can be observed.Again β is an important model parameter with respect to incision depth.The period of Q = 0 m 3 s −1 shows no downward motion of the channel bottom and can be considered to be representative of wintertime conditions.In the case of α > 0, ice deformation not only counteracts the incision process, but causes an uplift of the channel bottom during the winter regime.In reality, net accumulation within the channel in combination with dormant incision during wintertime adds complexity to channel closure processes.Channels can completely fill with snow over wintertime and disappear until spring.It is also possible that snow filled channels again carry meltwater streams at their base in the next melt season.As the stream starts in spring, snow in the lower part of the channel is washed away, which can lead to a snow bridge across the top of the channel, effectively resulting in a form of 'pinching off' the channel without the aid of ice creep.A complete removal of the filled-in snow by spring meltwater can also occurs.Thus purely ice creep driven channel closure is compounded in reality by surface mass balance processes and will be hard to observe in the field.
Currently the model is limited to simulate cases of open channel flow.A useful future expansion of the model would be to simulate the meltwater flow within the channel explicitly in a FEM simulation and couple it to the existing model.This would enable an explicit simulation of heat transfer processes at the water-ice boundary, which would allow for an evaluation of our melt rate distribution scheme and would be crucial for the treatment of cold ice conditions.A welcomed benefit would be to explicitly model the transition from open to pressurized channel flow.The simulation of fully developed turbulent water flow in channels and pipes is a numerically complex and difficult task, especially in the case of free surface flow, so this model extension remains a challenge for future research.
Because the model presented in this contribution is capable of simulating meltwater channel evolution dynamically, based on a state of the art ice dynamics model, we foresee that this approach holds great promise for glacierhydrological modelling applied to jökulhlaup evolution, moulin formation and evolution, surface lake drainage, quantifying eustatic sea level contribution of ice sheets, and even englacial ice-volcano interaction.

Fig. 1 .
Fig. 1.Sketch of the model geometry which indicates several model parameters as defined in Sects. 2 and 3.

Fig. 2 .
Fig.2.Forward stepping scheme for the moving melt boundary at each time step.In (a), melt processes expand the channel below the water level along the wetted perimeter (red line).Subsequently ice deformation closes the channel (b), which leads to the final geometry at this time step (c).The dotted line marks the channel shape from the prior time step, the dashed line the intermediate geometry after melt, and H max denotes the maximal water height inside the channel.Please note that ice deformation is greatly enhanced in (b) and neglected above the water line in this sketch to facilitate clarity.

Fig. 4 .
Fig. 4. Magnitude of velocity field for a simulation with Q = 1 m 3 s −1 , α = 0.0, β = 0.03, γ = 0.0, and ν = 1.0 at t p = 320 days (a) and t = 322 days (b).Between the two presented time steps, the channel pinches off from the surface and in (b), the numerical mesh already merged above the channel bottom.

Fig. 6 .Fig. 7 .
Fig. 6.Different channel shapes based on different values of ν are displayed in (a) and (b) demonstrates the independence of channel shape from the chosen initial surface depression (half ellipses with axes a and b).All results are based on our standard parameter set (Q = 1 m 3 s −1 , α = 0.0, β = 0.03, γ = 0.0, and ν = 1.0) with modifications indicated in the legend.Final channel geometries at t f for the standard set in (c), for ν = 1.5 in (d), for α = 1.0 in (e), and for Q = 10.0 m 3 s −1 in (f).

Fig. 9 .
Fig. 9. Evolution of channel incision for a time varying flux Q, different glacier surface slopes α, and different channel slopes β are plotted with different colors.All results are based on our standard parameter set (Q max = 1 m 3 s −1 , α = 0.0, β = 0.03, γ = 0.0, and ν = 1.0) with modifications indicated in the legend.The timestep (t p ) at which each channel pinches off is marked with a colored dot.All solutions displayed here use γ = 0. Model time t = 182.5 days, after which Q(t) = 0 m 3 s −1 for half a year, is marked with the vertical dashed line.In case of α = 1.0 (cyan line), the the channel closes by roof deformation at the end of the plotted line due to ice flow.In the case of β = 0.06 (magenta line), we plot results to t = 220 days, highlighting the different model behaviour for the first half of a year.All other results are manually stopped after 1.5 years.

The Cryosphere, 6, 493-503, 2012 www.the-cryosphere.net/6/493/2012/Table 1 .
Constants used in this manuscript.We assume temperate ice and take the corresponding value for A from Cuffey and Walder and Re the Reynolds numbers.As N u and Re change over time, mean values are reported.Values for our reference run are bold.† In this case the channel closes by roof deformation before the channel gets pressurizes.