Journal topic
The Cryosphere, 12, 971–991, 2018
https://doi.org/10.5194/tc-12-971-2018
The Cryosphere, 12, 971–991, 2018
https://doi.org/10.5194/tc-12-971-2018

Research article 22 Mar 2018

Research article | 22 Mar 2018

# Modelling seasonal meltwater forcing of the velocity of land-terminating margins of the Greenland Ice Sheet

Modelling seasonal meltwater forcing of the velocity of land-terminating margins of the Greenland Ice Sheet
Conrad P. Koziol1,a and Neil Arnold1 Conrad P. Koziol and Neil Arnold
• 1Scott Polar Research Institute, Cambridge, UK
• anow at: School of Geosciences, University of Edinburgh, Edinburgh, UK

Abstract

Surface runoff at the margin of the Greenland Ice Sheet (GrIS) drains to the ice-sheet bed, leading to enhanced summer ice flow. Ice velocities show a pattern of early summer acceleration followed by mid-summer deceleration due to evolution of the subglacial hydrology system in response to meltwater forcing. Modelling the integrated hydrological–ice dynamics system to reproduce measured velocities at the ice margin remains a key challenge for validating the present understanding of the system and constraining the impact of increasing surface runoff rates on dynamic ice mass loss from the GrIS. Here we show that a multi-component model incorporating supraglacial, subglacial, and ice dynamic components applied to a land-terminating catchment in western Greenland produces modelled velocities which are in reasonable agreement with those observed in GPS records for three melt seasons of varying melt intensities. This provides numerical support for the hypothesis that the subglacial system develops analogously to alpine glaciers and supports recent model formulations capturing the transition between distributed and channelized states. The model shows the growth of efficient conduit-based drainage up-glacier from the ice sheet margin, which develops more extensively, and further inland, as melt intensity increases. This suggests current trends of decadal-timescale slowdown of ice velocities in the ablation zone may continue in the near future. The model results also show a strong scaling between average summer velocities and melt season intensity, particularly in the upper ablation area. Assuming winter velocities are not impacted by channelization, our model suggests an upper bound of a 25 % increase in annual surface velocities as surface melt increases to present levels.

1 Introduction

Surface meltwater draining into the subglacial system drives seasonal acceleration of ice velocities at land-terminating sectors of the Greenland Ice Sheet (GrIS) margin. It may also be an important factor for seasonal acceleration of marine-terminating sectors . Increased water pressures reduce basal drag by decreasing ice–bed coupling, leading to faster ice flow. Early in the summer, surface runoff drains into an inefficient hydrological system, elevating water pressures, and accelerating ice flow . As the melt season progresses, a channelized system that efficiently drains water develops. This reduces water pressures and leads to a late summer deceleration . Understanding the impact of increased surface melting on the spatial and temporal evolution of basal hydrology is important for constraining the GrIS's future evolution. If increased summer melt intensity drives faster mean annual velocities, then a positive feedback between surface melt and ice flow would contribute to mass loss from the GrIS in a warming climate . Faster ice flow would draw ice down to lower elevations, where the melting is greater, which in turn drives faster ice flow.

Observations do not, however, show a simple relationship between surface runoff and ice velocities. Decadal-timescale observations in southwest Greenland of land-terminating sectors show mean annual velocities decreasing in the ablation zone . However, reported correlations between summer melt intensity and mean annual ice velocities from these studies are either slightly negative or nonexistent. In the accumulation zone, decadal-timescale measurements are sparse and the data inconclusive about velocity trends . Measurements on a daily timescale in the ablation zone show that increased melt intensity can lead to faster ice flow early in the summer. However, the impact of increased ice motion early in the summer on the average annual velocity can be offset by an earlier onset of channelization and corresponding deceleration as the melt season progresses . Increases in channelization extent may also lead to slower mean winter flow due to more extensive drainage of the subglacial system leading to lower water pressures during winter . As melt season intensity continues to increase, it remains unclear how ice velocities will be forced by water input at higher elevations where ice thickness is greater and whether patterns of water input at higher elevations will change .

Numerical models can provide insight into the hydrological processes driving faster summer flow. Recent subglacial hydrology models have progressed to simultaneously incorporating both distributed and efficient systems, explicitly treating the interaction between the two . Current models can reproduce the observed up-glacier development of the efficient system through the melt season. When coupled to an ice sheet model, the results broadly reproduce the observed velocity patterns of the GrIS margin . However, recent hydrological models coupled to ice flow models have not been applied to real domains of the GrIS to model large-scale behaviour of the ice margin during the summer melt season. Rather, applications to real domains of the GrIS for the summer melt season have either omitted ice flow , used a simplified hydrological model coupled to ice flow , or focused on a small domain . Coupling recent hydrological models with ice sheet models allows for important feedback between the distributed system and ice velocities and allows explicit comparison between GPS velocities and model output. Comparisons to surface ice velocity measurements are an important means for validating subglacial hydrological models and provide a method for constraining poorly understood aspects of subglacial hydrology (Flowers2015). Present challenges in applying coupled ice dynamics–hydrology models to the GrIS margin for modelling seasonal evolution include the values of parameters, the form of the sliding law which relates water pressures to basal drag, and whether the models presently include the necessary elements. Additionally, modelling surface hydrological input to drive the subglacial hydrology model is in itself a challenge. A variety of methods have been employed, incorporating different drainage elements (Banwell et al.2016; Bougamont et al.2014; de Fleurian et al.2016). However, no model including drainage via all of crevasses, moulins, and lake hydrofracture has been used to force a subglacial hydrology model to date.

This paper aims to model summer ice flow in the land-terminating Russell Glacier area of western Greenland for three contrasting melt seasons using a multicomponent model approach similar to the previous work of and on alpine glaciers. The model is then used to test the ice sheet response to higher melt input. A coupled hydrology–ice flow model is produced by integrating a subglacial hydrology model (Hewitt2013) with the ice flow model of . This coupled model is driven by surface input from the surface hydrology lake filling (SRLF) model from and initiated using the inversions from . The Russell Glacier area is selected as a study site to take advantage of the numerous observations available. These observations include radar flight lines constraining bed topography , meteorological data constraining climatic input , and GPS data , which provide a calibration and validation dataset for model output.

Figure 1Landsat 8 satellite image acquired on 19 August 2013, band 2, showing the Russell Glacier area. Black solid rectangle outlines the study domain for the integrated model, while the black dashed rectangle outlines the SRLF study domain. The blue triangles show the locations of GPS stations . Purple diamonds show the locations of automatic weather stations . Cyan circles show the locations of moulins used as tracer injections sites in . Inset shows the location in reference to Greenland.

2 Methods

The methods section begins with a description of the Russell Glacier study area and the datasets used. The study site is presented first so that the domain can be referred to when describing the boundary conditions applied in the models. Each individual model is then briefly described, before detailing how the models are linked. The coupled ice flow–subglacial hydrology model is referred to as the “integrated model” for simplicity. Finally, the modelling workflow is described.

## 2.1 Study area and datasets

The Russell Glacier area is a land-terminating sector of the GrIS in southwest Greenland. The study area's boundaries for the SRLF model and the integrated model are shown in Fig. 1. The domain of the SRLF runs is selected to be larger than the integrated model domain to minimize the impact of boundary conditions. A 6 km buffer is used at the northern and southern boundaries of the SRLF domain, based on the reported internally drained catchments by . The SRLF domain extends 8.5 km to east of the integrated model study site to capture as much higher-elevation melting as possible. The domain of the SRLF model is discretized at a 90 m resolution, while the domain of the integrated model is discretized at a 1000 m resolution.

Figure 2Daily surface runoff over the SRLF Russell Glacier study area for three contrasting summer melt seasons.

Two different topography datasets are used. The SRLF model is run using surface topography from the GIMP dataset . The high-resolution surface topography is necessary for accurate water routing and so that lake basin topography is accurately preserved. The integrated model is run with surface and bed topography from BedMachine2 to take advantage of the mass-conservation methods used to determine basal topography. BedMachine2 provides both topographic datasets at 150 m, although the true resolution is reported as 400 m. These data are reinterpolated to 1000 m resolution.

Surface runoff and snow depth data for the SRLF model are provided by RACMO2.3 . Both runoff and snow depth are bilinearly interpolated from 11 km. Three seasons with contrasting melt volumes are modelled: 2009, 2011, and 2012 (Fig. 2). Total melt over the SRLF study domain was 1.2×1010m3 in 2009, 1.7×1010m3 in 2011, and 2.1×1010m3 in 2012. Following , we use these 3 years as representative of summers with average, elevated, and extreme melt intensity, respectively.

Mean winter velocities are used for inversions of winter basal boundary conditions (see Sect. 2.6) and to determine crevasse locations as an input to the SRLF model. Mean winter velocities for 2008–2009 are provided at 500 m resolution by the MEaSUREs Greenland Ice Sheet Velocity Map dataset . For the inversion procedure, the winter velocities, along with their associated errors, are reinterpolated to 1000 m. Velocities at 500 m resolution are used to determine surface stresses, assuming an ice temperature of 5 C. Crevassed areas are then calculated using a von Mises stress criterion following . A crevassing threshold is selected by comparing the von Mises stress to observed patterns of crevassing in a Landsat 8 image, acquired on 19 August 2013. A threshold value of 145 kPa gave the best visual match.

Moulin locations are specified as input data in the SRLF model. Moulin locations in the Russell Glacier area reported by are used. These were derived automatically from a Landsat 8 image acquired on 19 August 2013, using an algorithm which determines where streams are observed to abruptly disappear . As in , moulin locations which do not coincide with a stream location calculated by the surface routing algorithm are slightly adjusted, such that they are located on a stream. A small number of moulins from the dataset are deleted, as they were not near a calculated stream and hence would drain negligible water.

A key validation dataset in the Russell Glacier area is GPS surface velocity measurements for 2009–2012 . Time series of hourly and daily averaged surface speeds are provided in the dataset. Here, the daily averaged speeds are used for comparison with model results. The locations of GPS stations are shown in Fig. 1.

## 2.2 Supraglacial hydrology (SRLF)

We use the SRLF model from , run at 90 m resolution. A no-inflow boundary condition is imposed on all boundaries. Water is routed using a digital elevation model (DEM) of the surface of the ice sheet and a single flow direction algorithm . Water can enter the subglacial system via drainage into pre-existing moulins or crevasses and is also allowed to drain off the edge of the domain or over the western ice margin. Water also collects in depressions in the surface DEM forming lakes. Lakes which are predicted to hydrofracture , using a fracture area criterion, drain to the ice–bed interface and create a surface-to-bed connection (treated as a moulin) for the remainder of the melt season. Lakes can also drain over the surface of the ice sheet via “overspill” drainage and “channelized” drainage. Overspill drainage refers to when water exceeding the capacity of the lake is routed downstream, with no incision of a channel at the lake edge. Channelized drainage refers to when water is routed downstream but incises a channel at the lake edge, which allows slow lake drainage. Channel incision is modelled following . Overspill and channelized drainage can occur simultaneously when water enters a lake faster than can be evacuated by an existing channel alone.

Table 1Constants used in the subglacial hydrology model during integrated runs in the Russell Glacier area.

## 2.3 Subglacial hydrology

We use the subglacial hydrology model presented in and . Distributed flow occurs through a continuum “sheet”, composed of a cavity sheet component and an elastic sheet component. The latter is included so that during lake hydrofracture events “hydraulic jacking” is simulated. Channels can form along the edges and diagonals of the rectangular finite difference mesh. Dissipative heating over an incipient channel-width length scale provides the initial perturbation for channel initialization. Water input occurs at moulins located at cell nodes, which, along with an englacial aquifer, allow for water storage. The model is run at 1000 m resolution. At the ice margin edge an atmospheric pressure boundary condition is imposed, while the remaining boundaries have a no-flux condition. A concise model description is given here following and to provide context for the parameters used. However, for a detailed description the reader is referred to and .

Table 2Constants used in the ice sheet/inversion model applied to the Russell Glacier area.

Discharge in the continuum sheet is modelled as

$\begin{array}{}\text{(1)}& \mathbit{q}=-\frac{{K}_{\mathrm{s}}{h}^{\mathrm{3}}}{{\mathit{\rho }}_{\mathrm{w}}g}\mathrm{\nabla }\mathit{\varphi },\end{array}$

where h(x,y) is the thickness of the continuum sheet, Ks is the sheet flux coefficient, g is the acceleration due to gravity, ρw is the density of water, and ϕ is the hydraulic potential. The hydraulic potential is defined as $\mathit{\varphi }\left(x,y\right)={\mathit{\rho }}_{\mathrm{w}}gb\left(x,y\right)+{p}_{\mathrm{w}}\left(x,y\right)$, where pw is water pressure and b(x,y) is the bed elevation.

The distributed sheet thickness (h) is the sum of the thickness of the cavity sheet (hcav) and the elastic sheet (hel). The cavity sheet evolves according to

$\begin{array}{}\text{(2)}& \frac{\partial {h}_{\text{cav}}}{\partial t}=\frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}m+{U}_{\mathrm{b}}\frac{{h}_{\mathrm{r}}-{h}_{\text{cav}}}{{l}_{\mathrm{r}}}-\frac{\mathrm{2}{A}_{\mathrm{b}}}{{n}^{n}}{h}_{\text{cav}}|N{|}^{n-\mathrm{1}}N,\end{array}$

where ρi is the density of ice, m is the basal melting rate, Ub is the basal sliding speed, hr is the bed roughness height scale, lr is the bed roughness length scale, Ab is the ice creep parameter, n is the exponent from Glen's flow law, and N(x,y) is the effective pressure. The effective pressure is defined as $N={\mathit{\rho }}_{\mathrm{i}}gH-{p}_{\mathrm{w}}$, where H is the ice thickness.

Basal melt rate is given by

$\begin{array}{}\text{(3)}& m=\frac{G+{\mathbit{\tau }}_{\mathbf{b}}\cdot {u}_{\mathrm{b}}}{{\mathit{\rho }}_{\mathrm{w}}L},\end{array}$

where ${\mathbit{\tau }}_{\mathbf{b}}=\left({\mathit{\tau }}_{\text{b}x}\left(x,y\right),{\mathit{\tau }}_{\text{b}y}\left(x,y\right)\right)$ is the basal drag, ${u}_{\mathrm{b}}=\left({u}_{\text{b}},{v}_{\text{b}}\right)=\left(u\left(x,y,b\right),v\left(x,y,b\right)\right)$ is the basal velocity, G is the net conductive flux, defined as the geothermal heat flux minus conductive loss into the ice, and L is latent heat.

The elastic sheet thickness is given by

$\begin{array}{}\text{(4)}& {h}_{\text{el}}={C}_{\text{el}}\left[-{N}_{-}+\frac{\mathrm{1}}{\mathrm{2}}{N}_{\mathrm{0}}max{\left(\mathrm{0},\mathrm{1}-\frac{{N}_{+}}{{N}_{\mathrm{0}}}\right)}^{\mathrm{2}}\right],\end{array}$

where ${N}_{-}=min\left(N,\mathrm{0}\right)$, ${N}_{+}=max\left(N,\mathrm{0}\right)$, Cel is an elastic compliance, and N0 is a regularization parameter. When effective pressure is positive, this layer is designed to be zero. As effective pressure approaches zero or is negative, the thickness is determined by the product of the elastic compliance and effective pressure .

Discharge in channels is modelled as

$\begin{array}{}\text{(5)}& Q=-{K}_{\mathrm{c}}{S}^{\mathrm{5}/\mathrm{4}}{\left|\frac{\partial \mathit{\varphi }}{\partial r}\right|}^{-\frac{\mathrm{1}}{\mathrm{2}}}\frac{\partial \mathit{\varphi }}{\partial r},\end{array}$

where Kc is a turbulent flow coefficient, S is channel cross section, and r is along channel distance.

The channel cross section evolves according to

$\begin{array}{}\text{(6)}& \frac{\partial S}{\partial t}=\frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}M-\frac{\mathrm{2}{A}_{\mathrm{b}}}{{n}^{n}}S|N{|}^{n-\mathrm{1}}N,\end{array}$

where M is the melting rate along the channel wall.

The melting rate along the channel walls is given by

$\begin{array}{}\text{(7)}& M=\frac{|Q\frac{\partial \mathit{\varphi }}{\partial r}|+{\mathit{\lambda }}_{\mathrm{c}}|\mathbit{q}\cdot \mathrm{\nabla }\mathit{\varphi }|}{{\mathit{\rho }}_{\mathrm{w}}L},\end{array}$

where λc is an incipient channel-width length scale (melting of basal ice over this scale contributes to channel initialization).

The equation for mass conservation is

$\begin{array}{ll}\frac{\partial h}{\partial t}& +\mathrm{\nabla }\cdot \mathbit{q}+\left[\frac{\partial S}{\partial t}+\frac{\partial Q}{\partial r}\right]\mathit{\delta }\left({\mathbit{x}}_{\mathrm{c}}\right)+\frac{\partial \mathrm{\Sigma }}{\partial t}\\ \text{(8)}& & =m+M\mathit{\delta }\left({\mathbit{x}}_{\mathrm{c}}\right)+R\mathit{\delta }\left({\mathbit{x}}_{\mathrm{m}}\right),\end{array}$

where Σ is englacial storage and R is the supraglacial input rate to moulins. The delta functions apply along channels (δ(xc)) and the positions of moulins (δ(xm)).

Englacial storage is represented as

$\begin{array}{}\text{(9)}& \mathrm{\Sigma }=\mathit{\sigma }\frac{{p}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}g}+{A}_{\mathrm{m}}\frac{{p}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}g}\mathit{\delta }\left({\mathbit{x}}_{\mathrm{m}}\right),\end{array}$

where σ is englacial void fraction and Am is moulin cross-sectional area.

Model parameters held constant are shown in Table 1. Two parameters of the subglacial hydrology model, Ks and σ are the focus of calibration experiments. They were identified in and as key parameters determining the morphology of the subglacial hydrology system.

## 2.4 Ice flow/inversion

The ice flow model implements the hybrid formulation of the ice sheet stress balance , which can be considered a combination of shallow ice approximation and shallow shelf approximation. The model implicitly accounts for depth-varying ice flow, and surface velocities can be explicitly calculated when comparing model output to GPS measurements. This model is similar to the one used in , except the conservation of momentum equations are a function of depth-integrated velocities rather than basal velocities. Parameters for the model are listed in Table 2. A Dirichlet boundary condition is imposed on all lateral domain margins except the ice margin, where the standard boundary condition based on the continuity of stress is used. A no-penetration boundary condition is applied at the edge of the nunatak (Fig. 1). Three sliding laws are implemented:

$\begin{array}{}\text{(10)}& & {\mathbit{\tau }}_{\mathbf{b}}={\mathit{\beta }}^{\mathrm{2}}{u}_{\mathrm{b}},\text{(11)}& & {\mathbit{\tau }}_{\mathbf{b}}={\mathit{\mu }}_{\mathrm{a}}\phantom{\rule{0.25em}{0ex}}{N}_{\text{sl}}^{p}{{U}_{\mathrm{b}}}^{q}\frac{{u}_{\mathrm{b}}}{{U}_{\mathrm{b}}},\text{(12)}& & {\mathbit{\tau }}_{\mathbf{b}}={\mathit{\mu }}_{\mathrm{b}}{N}_{\text{sl}}{\left(\frac{{U}_{\mathrm{b}}}{{U}_{\mathrm{b}}+{\mathit{\lambda }}_{\mathrm{b}}{A}_{\mathrm{b}}{N}_{\text{sl}}^{n}}\right)}^{\frac{\mathrm{1}}{n}}\frac{{u}_{\mathrm{b}}}{{U}_{\mathrm{b}}},\end{array}$

where β(x,y) is a basal drag coefficient, μa(x,y) is a drag coefficient, p and q are positive exponents, μb(x,y) is a limiting roughness slope, and λb is a bed roughness length (Hewitt2013). Following , negative effective pressures are eliminated by setting ${N}_{\text{sl}}=max\left(N,\mathrm{0}\right)$ and regularized with a small constant (102Pa).

The linear sliding law (Eq. 10) is used for the initial inversion of winter mean velocities, while the Budd (Eq. 11) and Schoof (Eq. 12) sliding laws are used subsequently. The linear sliding law uses a single parameter to represent all the processes at the ice–bed interface, while the non-linear sliding laws attempt to explicitly incorporate the impact of effective pressure and have a more complex dependence on velocity.

The inversion code used in this paper is described in . It is based on automatic differentiation methods and uses the open-source MATLAB package ADiGator . The gradient of the cost function in this method is equivalent to one calculated using Lagrangian multiplier methods to generate the adjoint model . The cost function (Eq. 13) has two terms. The first is the weighted square of the differences of measured and predicted velocities. The second is a Tikanov regularization term added for stability.

$\begin{array}{}\text{(13)}& J={\mathit{\gamma }}_{\mathrm{1}}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }w\cdot \left({U}_{\text{obs}}-{U}_{\mathrm{s}}\left(\mathit{\alpha }\right){\right)}^{\mathrm{2}}\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{s}}+{\mathit{\gamma }}_{\mathrm{2}}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }\left(\mathrm{\nabla }\mathit{\alpha }\cdot \mathrm{\nabla }\mathit{\alpha }\right)\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{b}},\end{array}$

where γ1 and γ2 are scaling factors, Γs is the surface domain, Γb is the basal domain, w(x,y) is a weighting function, Uobs(x,y) are observed surface ice speeds, ${U}_{\mathrm{s}}\left(\mathit{\alpha },x,y\right)$ are modelled surface speeds, and α(x,y) is the control parameter. The control parameter depends on the sliding law, and represents β in the linear sliding law, μa in the Budd sliding law, and μb in the Schoof sliding law. The inverse of reported errors of surface velocities are used as weights. Modelled surface velocities depend on the control parameter via the sliding law. The inversion procedure minimizes the cost function with respect to the control parameter.

## 2.5 Model integration

The SRLF model is used to determine supraglacial input rates to the subglacial system. For model integration, we assume that water drainage through surface-to-bed connections is strictly vertical, with no horizontal component. The SRLF model routes water into three different surface-to-bed pathways: moulins, lake hydrofracture, and crevasses. Moulins and surface-to-bed connections from lake hydrofracture are treated identically. All water entering these cells drains into the subglacial hydrology system at that location. However, drainage through crevasse fields requires additional consideration. When water enters a crevassed cell in the SRLF model, no further routing occurs. Since it is unlikely that every crevassed grid cell drains water locally to the ice–bed interface, postprocessing of SRLF output is necessary.

Figure 3Schematic drawing showing the conceptual model of the crevasse drainage implemented. Shaded area indicates a crevasse field. Moulins are assumed to occur where high-flux supraglacial streams intersect the crevasse field. The crevasse field is then partitioned using Voronoi partitioning into internal catchments. All melt which occurs within an internal catchment is assumed to drain into the subglacial system at the corresponding moulin.

Water drainage through crevasse fields is poorly understood, and the scheme implemented here (Fig. 3) is motivated by simplicity. We assume all water in crevasse fields drains to bed, neglecting any refreezing. We also assume that contiguous areas of crevassed cells are hydrologically connected. A crevassed cell in the SRLF model can accumulate water from two sources: (1) local ablation predicted by RACMO2.3 and (2) water flow from adjacent non-crevassed cells if the cell is on the margin of a crevassed area. Modelling predicts approximately 70 % of the water drained by crevasses is intercepted water flow over the ice sheet surface (source 2). This water is concentrated at the points where supraglacial streams intersect the crevasse fields. The model assumes that moulins exist at these points, as high water input would be favourable to nucleating and sustaining moulins. Moulins are only placed in cells with sufficient drainage, determined by a volume threshold. A value of 5×105m3 is selected, corresponding to approximately the median volume drained by moulins outside of lake basins. A lower threshold results in a rapidly increasing number of moulins draining smaller amounts of water. A Voronoi partitioning is then used around the inferred moulins to create internal catchments within the crevasse field. All water in a catchment is assumed to drain into its corresponding moulin. As stated, the SRLF model does not route water within crevasse fields; there is no travel time associated with melt in the internal catchments of crevasses and the moulin.

The supraglacial model is run independently to determine a time series and location of water inputs to the base. These are then used as input to the coupled subglacial hydrology–ice flow model. A potential feedback omitted by this approach is the influence of surface velocity on lake hydrofracture. However, the current design and computational requirements of the SRLF model make it impossible to run in a fully coupled manner with the integrated model.

The integration of the subglacial hydrology and ice flow models mirrors that of ; the subglacial hydrology uses an implicit timestep using the current ice velocity distribution. After the state of the subglacial hydrology model in the next timestep is calculated, the ice model is called to update ice velocities. At each timestep, the basal melting rate is updated. The geometry of the domain is kept constant for the whole run.

Figure 4Flow chart showing the work flow for initializing and running the integrated model.

## 2.6 Workflow

Figure 4 shows the workflow for initializing and running the integrated model. The initial step is to perform an inversion using the linear sliding law over the study area (see , for details of inversions, which are done over the same study area but with a resolution of 500 m). All linear inversions are run using mean winter velocities from 2009, the most recent year for which data were available. This inversion provides an initial distribution of basal drag and basal velocities to calculate the basal melt rate (Eq. 3). The subglacial hydrology model is then run for 240 days holding basal velocities fixed, corresponding to a run over a winter season (1 September–30 April). By the end of the run, effective pressures reach an approximate steady state . The effective pressures at the end of the subglacial hydrology simulation are then incorporated into an inversion with a non-linear sliding law to determine the background values of the coefficients, μa and μb. These sliding law coefficients, the basal water pressures, and the surface runoff input from the SRLF model form the inputs to the integrated model. The integrated model is then run for the summer melt season using the end of winter effective pressures as an initial condition. As stated in , a key assumption of this procedure is that the mean winter velocities are valid both at the beginning and end of the winter season. Although winter velocities are not constant, published GPS records in southwest Greenland of winter velocities show limited variability .

The inversions are run with constant parameters. Both the winter subglacial hydrology run and the subsequent integrated model runs use the same parameters. A parameter search therefore requires performing the inversion using an effective pressure-dependent sliding law for each set of parameters tested.

## 2.7 Simulations

We run five main simulations along with those used in the sensitivity analysis (not shown). Two simulations are run calibrating the model using data and inputs for 2009 and 2011. Another simulation is then run validating the model with data and inputs for 2012. Two potential future melt scenarios are simulated by using and the modelled supraglacial input to the subglacial system for 2011. These are referred to as “2011×2” and “2011×4”, respectively. The aim of these scenarios is to investigate potential changes in the behaviour of the subglacial system rather than to model a melt season or reliably predict future ice velocities. Accurate predictions of ice velocities would not only require predicted surface runoff but also depend on predicting changes in ice sheet topography and predicting the future distribution of supraglacial drainage pathways. Addressing these issues requires careful consideration and is beyond the scope of this paper.

3 Results

## 3.1 Meltwater input partitioning

The majority of supraglacial meltwater drains into the englacial system (Table 3), consistent with observations and previous modelling . Approximately 12 % of supraglacial lakes are predicted to hydrofracture. These events drain only a small percentage of surface runoff (1.3 %). Most drainage (86.1 %) occurs through features modelled as moulins: crevasses, surface-to-bed connections subsequent to lake hydrofracture, and moulins outside of lake basins. Of the water drained by crevasses, approximately 30 % is generated locally via ablation in crevassed cells, while 70 % is routed into crevasses. Water routing into crevasses is concentrated in a small number of cells, with 50 % of the water routed into crevasses entering in only 100 of the 7573 cells forming the perimeter of crevasse fields. Crevasse drainage is concentrated near the ice margin (Fig. 5), while drainage into other pathways occurs throughout the study area.

Table 3Surface runoff partitioning into different meltwater pathways for the 2009 melt season in the SRLF domain. Water flowing over the western boundary is categorized as “ice margin”, while water flow over the lateral boundaries is labelled as “lateral outflow”. “Remaining flow” refers to water still flowing over the ice sheet at the end of the model run. Water flowing into crevasses and moulins is in categories “crevasses” and “moulins”, respectively. “Lake storage” refers to water in lakes at the end of the simulation. “Lake hydrofracture lake” refers to the water in lakes that is drained by hydrofracture events themselves. “Lake hydrofracture moulin” refers to water drainage into the subsequent surface-to-bed connections from hydrofracture events.

Figure 5Modelled supraglacial input in the Russell Glacier area for the integrated model domain in 2009. Meltwater pathways are denoted by circles of different colours, with red, green, and blue corresponding to moulins, crevasses, and lakes, respectively. Circle areas are scaled by volume. GPS stations are shown as black triangles. Hatch marks show grid cells calculated as crevassed. Crevasse inputs appear within hatched areas due to resampling from 90 to 1000 m resolution. Background is basal topography from BedMachine2 reinterpolated at 1000 m. Light grey contours correspond to 100 m basal contours. Black lines correspond to 200 m surface topography contours at the same elevations as in Fig. 1.

## 3.2 Calibration

Model parameters are calibrated by qualitatively comparing modelled velocities to GPS measurements of horizontal surface velocities from 2009 (Fig. 6) and 2011 (Fig. 7). The focus of the calibration is on parameters identified as key to determining the morphology of the subglacial system, Ks and σ . The calibration resulted in the two parameters being assigned spatially heterogeneous distributions. The σ field is assigned a background value of 10−3, with 50 % of the cells then randomly set to 10−4. The Ks field is constructed using a background value of 10−2${\mathrm{Pa}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, with 15 % of the cell nodes randomly assigned a value of 10−7${\mathrm{Pa}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$. Since Ks is defined on the grid, neighbouring nodes are averaged in the x and y directions to determine values on edges. At a sheet depth of 0.1 m, a Ks value of 10−2 results in an effective hydraulic conductivity Ksh2 of 10−4m s−1 (Hewitt2013). This is at the upper end of values for till, which are inferred to be 10−4 to 10−9m s−1 . The secondary value of 10−7${\mathrm{Pa}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ assigned to Ks was selected to give an effective hydraulic conductivity at the opposite end of the spectrum.

The calibration focused on matching the duration and magnitude of speedup events, as the timing of the events is controlled by surface input, which was not sensitivity tested, given we used the same time series of surface runoff for each year from RACMO2.3 in all SRLF runs. The model is calibrated using the Budd sliding law, and the same parameter values are used in the simulations with the Schoof sliding law. Figures 6 and 7 show model velocities output at noon for the Budd sliding law as representative of the daily average and daily averages calculated from output at 6 h intervals for the Schoof sliding law. Subdaily variability in model output is subdued, except during periods of high velocities in simulations using the Schoof sliding law (see ). Model output values shown are from the summer immediately following the winter initialization. There are only minor differences between this model output and from running the model for an additional year and using the output from the second summer. Since surface water input to the subglacial hydrological system is a key driver of ice velocities, surface runoff from RACMO2.3 and nearby surface ablation rates determined at weather stations are plotted alongside velocities. RACMO2.3 surface runoff forces modelled ice flow, while the weather station ablation rate is taken as representative of the water input driving measured ice velocities. Some caution is necessary comparing the datasets, since RACMO2.3 accounts for both refreezing of meltwater and precipitation events. Refreezing, however, should only be a small component . An error of 5 % is estimated for the calculated daily ablation rates .

Figure 6Modelled ice velocities plotted against GPS measurements for the 2009 melt season. Daily average horizontal velocity from GPS measurements are plotted in blue. Modelled velocities using the Schoof sliding law and Budd sliding law are plotted in black and red, respectively. Daily ablation from weather stations are shown in shaded blue, while RACMO2.3 surface runoff is shown in shaded red. Locations of GPS and weather station sites are shown Fig. 1. Weather station ablation rates are plotted at the nearest GPS site.

Figure 7Modelled ice velocities plotted against GPS measurements for the 2011 melt season. Daily average horizontal velocity from GPS measurements are plotted in blue. Modelled velocities using the Schoof sliding law and Budd sliding law are plotted in black and red, respectively. Daily ablation from weather stations are shown in shaded blue, while RACMO2.3 surface runoff is shown in shaded red. Locations of GPS and weather station sites are shown Fig. 1. Weather station ablation rates are plotted at the nearest GPS site.

The Schoof and Budd sliding laws result in model output of comparable fit to the measured velocities for large segments of the velocity time series. However, during periods of high velocities, the Schoof law can overpredict the magnitude of the velocity by a factor of 3. Model output with the Schoof sliding law is also observed to have a sharper and higher magnitude summer speedup, as well as a slight increase in velocity variability. Since the Budd sliding law results in an overall better match to the measured velocities, the analysis of the velocity time series in the remainder of the paper focuses on those results.

The model predicts low ice velocities throughout most of the summer melt season at the three GPS sites closest to the margin. In general, measured GPS velocities are also relatively low, except for early season high magnitude variability observed at sites S1–S3 (Fig. 6a–c) in 2009 and at sites S1–S2 in 2011 (Fig. 7a–b). This observed variability in the GPS velocities precedes melt predicted by RACMO2.3 and is not reproduced by the model. At sites S1 and S2, modelled velocities show some limited acceleration during the early summer in both 2009 and 2011, but to a lesser extent than in the observed velocities. The fit improves at site S3 for both years, as modelled velocities in 2009 approximate observed velocities (after the early onset of increased summer velocity in the observed data), and in 2011 the model predicts the early summer speedup more effectively.

Modelled velocities at sites S4–S5 (Figs. 6d–e and 7d–e) capture the seasonal trend of ice flow and mirror some of the observed short-term speedup events. Modelled velocities at S4 match the general flow while diverging from GPS measurements during periods of observed and modelled enhanced flow. In 2011 modelled ice flow shows similar short-term speedup events as the GPS measurements, such as those beginning on days 198 and 237. At site S5, the model does not predict the gradual speedup observed in the GPS velocities. Similar to the GPS measurements, there is a brief period of enhanced flow in mid-summer, followed by a slowdown. In 2011, however, the model captures the early velocity speedup and the general trend through the remainder of the summer, including the same speedup events observed at site S4.

Model velocities underpredict the measured velocities at the highest sites. In 2009, measured velocities at site S6 (Fig. 6f) show a gradual increase in the first half of the melt season, followed by a gradual decline in the second half. Neither the increase nor decrease in velocity mirrors the weather station ablation rate. In contrast, model velocities are observed to be enhanced in the middle of summer, mirroring modelled melt. Site S6 (Fig. 7f) measurements show faster flow in 2011 than in 2009. The model velocities match the initial velocity increase observed in GPS velocities but do not reach the same magnitude. A late summer slowdown is observed in both the modelled and measured velocities, as are short-term increases in velocities at days 200 and 240. At site S7, modelled velocities depart from the winter mean by a few metres per year in both 2009 and 2011 (Figs. 6g and 7g). Measurements show an increase on the order of 10–20 m yr−1 in both 2009 and 2011.

In summary, the early summer speedup and subsequent mid-summer slowdown at mid-elevations are captured. The model is also able to reproduce the pattern of synchronous speedups observed at multiple adjacent GPS stations. Consistent features not captured are early summer variability at low sites, short-term variability, and late summer deceleration below the winter mean. Modelled velocities are only observed to flow slower than the winter velocity mean for a period of a few days and by a small magnitude (<5m yr−1).

Ablation rates calculated from automatic weather stations near to sites S2, S4, and S6 are comparable to predicted RACMO2.3 surface runoff. The two datasets show similar magnitude at sites S2 and S4, with higher variability in ablation than runoff. At S6, both ablation and predicted runoff are similar in 2009, while in 2011 ablation is approximately twice the magnitude of surface runoff and has a much higher variability. Qualitatively, model velocities at sites S1–S3 do not correspond with predicted surface runoff, while they do correspond at sites S4–S6. GPS measurements do not in general correspond with the ablation rate at S2 and only weakly correspond at S4 and S6.

## 3.3 Model sensitivity

Calibrating the integrated model is an underdetermined problem. Multiple parameters in each cell across the grid are constrained using only seven time series of point GPS data. The parameters selected for the subglacial hydrological model are not unique in giving a qualitatively good fit. Within the parameter space searched, different sets of parameters either enhanced or dampened the magnitude of the velocity output or resulted in a velocity signal that significantly diverged from GPS measurements. Extensive sensitivity analysis of the subglacial hydrology component of the integrated model to parameters are conducted in and . In this section we focus on the sensitivity of the model to the setup, to Ks, and to σ.

Drainage through crevasses is poorly constrained, and hence the impact of varying crevasse drainage is tested. Velocities at the GPS stations are not found to be sensitive to variations in crevasse drainage. The standard value of the moulin volume threshold of 5×105m3 resulted in crevasse input partitioning into 182 moulins and internal catchments. Changing the threshold value to 105m3 and 106m3 resulted in 337 and 122 internal catchments, respectively. Model output in both scenarios showed negligible changes. Similarly, neglecting water generated over crevasse fields and only using water flowing into the crevasse fields from external streams had little impact on modelled velocities at the GPS stations.

Lake hydrofracture events result in a large volume of water rapidly draining to the base during the event itself and a surface-to-bed connection which drains water for the remainder of the melt season. The impact of the initial rapid delivery of water on the model behaviour is tested by running a simulation where the water in the lake, when hydrofracture occurs, was not input to the base, but removed from the system. The impact on modelled ice velocities at GPS stations (none of which are near lakes which undergo hydrofracture) was found to be negligible. The season-long average velocity across the catchment was also not affected.

A heterogeneous sheet flux coefficient field is found to benefit the fit of modelled velocities by increasing the magnitude of the early summer speedup. The results using a constant value of 10−2${\mathrm{Pa}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ are overall very similar to the calibrated runs, while decreasing the value to 10−3${\mathrm{Pa}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ leads to model output with prolonged periods of velocities exceeding 400 m yr−1. Increasing the number of grid nodes with the sheet flux coefficient assigned a value of 10−7${\mathrm{Pa}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ from 15 to 30 % had a minor impact. Assigning 50 % of the grid nodes resulted in a worsening fit early in the summer at site S4 but had little impact at other sites or beyond the initial speedup. Patterning low value nodes into 4×4 patches, randomly seeded at 125 points, was also tested. The number of patches was selected so that if there was no overlap of the patches, 20 % of the nodes would be assigned a lower conductivity. Two simulations were conducted with different random locations of patches. One simulation strongly impacted the early summer speedup at site S3 and S4, while the other had a similar effect but on sites S4 and S5.

The GPS records in 2009 and 2011 show differing characteristics, with ice velocities in 2009 showing much less variability and more gradual changes than 2011. The choice of the parameter value setting for englacial storage (σ) attempts to balance the fit in both years. A better fit was observed with increased englacial storage ($\mathit{\sigma }={\mathrm{10}}^{-\mathrm{3}}$) for 2009 and decreased englacial storage ($\mathit{\sigma }={\mathrm{10}}^{-\mathrm{4}}$) in 2011. Increased capacity of englacial storage had the effect of dampening the velocity output. In 2009, this increased the fit of the model predictions by reducing the high velocities observed at sites S3–S5 between days 180 and 200. However, increased englacial storage also reduced the velocity speedups observed in 2011, particularly around day 205, reducing the fit to GPS measurements.

Figure 8Modelled ice velocities plotted against GPS measurements for the 2012 melt season. Daily average horizontal velocity from GPS measurements are plotted in blue. Modelled velocities using the Schoof sliding law and Budd sliding law are plotted in black and red, respectively. Daily ablation from weather stations are shown in shaded blue, while RACMO2.3 surface runoff is shown in shaded red. Locations of GPS and weather station sites are shown Fig. 1. Weather station ablation rates are plotted at the nearest GPS site.

## 3.4 Validation

The integrated model is validated against GPS velocity measurements from 2012 (Fig. 8), using the parameter values determined in the calibration. The pattern of modelled velocities at sites S1 and S2 are similar to those in 2011, with a moderate early velocity speedup followed by a gradual slowdown for the remainder of the summer. Unlike previous years, GPS velocities at site S1 do not exhibit high magnitude velocity variations, improving the match of the modelled velocities. Although the integrated model does not respond strongly to melt input for most of the summer at site S1, it does predict slightly elevated velocities driven by late season input around days 255 and 265, in line with GPS measurements. At site S2, the general pattern of speedup observed in the GPS velocities is mirrored by the modelled velocities. However, the magnitudes are consistently under predicted, particularly those of the short-term high magnitude speedups. The magnitude of modelled velocities improves at site S3 and especially S4, matching the timing but overpredicting the magnitude at S3, and matching both magnitude and timing of events at site S4. Little GPS data are available at sites S5 and S7. Similarly to previous years, model output underpredicts GPS velocities at site S6.

## 3.5 Increased melt scenarios

The high melt scenarios show faster flow early in the summer and higher peak velocities (Fig. 9). After the early summer speedup at sites S1–S3, simulations 2011, $\mathrm{2011}×\mathrm{2},$ and 2011×4 predict broadly similar velocities, although the 2011×2 and 2011×4 runs show slightly lower values overall; this effect can be seen most clearly at site S3. As elevation increases, modelled velocities in the high melt intensity runs decrease more during slow-flow/low-melt periods (e.g. days 210–230 at S4). At sites S4–S6, increased melt intensity results in higher variability of ice flow, with higher peak velocities in the first half of the summer season, but with similar low velocities during low-melt periods. The relative increase of peak velocities between 2011×2 and 2011×4 is greater than between 2011 and $\mathrm{2011}×\mathrm{2}.$ At these sites, the model predicts broadly similar velocities in all three simulations for the latter half of summer, from days 210 to 238, but at site S4, model velocities are lowest for the 2011×4 scenario, whereas at site S6, modelled velocities increase slightly with greater melt input. Site S5 shows mixed behaviour, with model velocities from simulation 2011×4 higher than the other simulations between days 210 and 222 but lower between days 223 and 236. Between days 210 and 238 at sites S4–S6, model velocities are low and only slightly elevated above their winter values. At site S7, the velocities in the increased melt intensity simulations are faster than 2011 velocities, and with some periods where the 2011×4 scenario generates the slowest velocities (e.g. days 198 and 210), when short-term melt rates drop after a period of higher melt and velocity. Starting at day 238, a late season velocity spike is observed, most strongly at sites S4 and S5, though apparent at other sites. The melt input during this period decreases with elevation, but the impact of this event increases strongly with elevation to sites S4–5 and then decreases at sites S6 and S7.

Figure 9Modelled ice velocities using a Budd Sliding law plotted for the 2011 melt season (blue), the melt scenario (magenta), and for the melt scenario (black).

Figure 10Averaged modelled melt season velocities at each of the GPS sites for the different years and future melt scenarios. Ice surface elevation also shown.

## 3.6 Average melt season velocities

Melt season averaged modelled velocities at the GPS sites are shown in Fig. 10. Average velocities are highest at GPS site S4 and decrease towards the ice margin and at high elevations. Average velocities increase with melt season intensity at all GPS sites, with a pattern skewed away from the sites closest to the ice margin which show the least sensitivity. As melt season intensity increases, velocities in the upper ablation zone and at the equilibrium line (located at 1500 m elevation, slightly above S6; ) are predicted to increase the most, scaling with melt (comparing 2011, $\mathrm{2011}×\mathrm{2},$ and 2011×4); site S7, with the lowest melt, shows more limited sensitivity. The pattern observed at the GPS stations is generally reflective of that across the study domain (Fig. 11); areas between around 800 and 1400 m show the largest increase, but with areas of slower flow predicted to accelerate more than areas of faster flow. Average velocities between 2009 and 2011×4 increase by up to 70 %.

Figure 11(a) Map of melt season average velocities for 2009. (b) Map of melt season average velocities for $\mathrm{2011}\phantom{\rule{0.125em}{0ex}}×\phantom{\rule{0.125em}{0ex}}\mathrm{4}.$ (c) Change (%) between the melt season average velocities of 2009 and $\mathrm{2011}\phantom{\rule{0.125em}{0ex}}×\phantom{\rule{0.125em}{0ex}}\mathrm{4}.$

## 3.7 Channel network morphology and extent

The development of the channelized system (see videos in Supplement) is similar to that observed in previous modelling studies (Banwell et al.2016; Hewitt2013; Werder et al.2013) and as inferred from observations . Channelization of the hydrological system begins at the margin and develops progressively up ice sheet. As channelization develops up-ice, the system evolves to an arborescent morphology. The up-ice extent of channelization increases with summer melt intensity (Fig. 12). In 2009, channels occur primarily below the 1000 m surface elevation contour. The extent increases past 1100 m, and approaches 1200 m, in 2012. As melt season intensity increases from 2009 to 2012, pockets of channelization at higher elevations are seen. The maximum extent of channelization occurs at approximately the same time in each modelled melt season and was qualitatively identified to occur between days 220 and 225 in all three melt seasons. Although the extent of channelization varies between 2009, 2011, and 2012, there are no significant differences in the organization of the channelized system. In the future scenario 2011×4 the morphology of the channelized system is similar to that in the modelled melt seasons. However, the extent increases further upstream past 1300 m and approaches 1400 m.

Figure 12Channelized system at maximum extent: (a) 2009; (b) 2011; (c) 2012; (d) 2011×4. Moulin locations used as tracer injections sites in are shown in purple. Two moulin locations are unlabelled for clarity. Black lines correspond to 200 m surface topography contours at the same elevations as in Fig. 1.

Figure 13Time series of discharge in the distributed (“sheet”) and channelized (“channels”) system for the 2009, 2011, and 2012 summers and the 2011×4 future melt scenario.

Figure 12 shows the locations of moulins used as tracer injection points in . Dye-tracing experiments by were performed in the summers of 2009 to 2011. Except for moulin IS39, tracers injected into the moulins drained from the subglacial system at an outlet located near moulin L1. Tracers injected into IS39 are reported to drain from an outlet of an adjacent catchment. The channel morphology in the modelled melt season output does not predict a major outlet located near L1 or that L41 and L57 will drain near L1. However, the model does predict that IS39 is on a different branch of the channelized system. Based on tracer measurements in 2011, report that channelization extends to at least L41, but not as far as L57. The modelled channelized system during 2009, 2011, and 2012 is in line with that result.

## 3.8 Distributed and channelized discharge

Water flow beneath the ice sheet is modelled to occur in interacting distributed and channelized systems. The discharge in each system follows similar trends for all three modelled melt seasons (Fig. 13). In 2009, 2011, and 2012, integrated discharge over the summer melt season in the channelized system is slightly less than half (43–48 %) of the integrated discharge in the distributed system. Modelled discharge begins to increase simultaneously in both systems at the start of the melt season. In 2011 and 2012, discharge in the distributed system rapidly increases in the early melt season. This is followed by a long period with overall high flow but with strong variations. At the end of the melt season, discharge in the distributed system rapidly decreases. In 2009, the early season increase in discharge is less rapid and more prolonged, and discharge peaks before decreasing to a plateau, after which it rapidly decreases. Discharge in the channelized system increases at a much slower rate and tends to increase until mid–late summer. It mirrors many of the short-timescale variations in the distributed system but with a dampened magnitude. At the end of the melt season, discharge in the distributed system decreases at a higher rate than in the channelized system. In 2011 and 2012, this results in a brief period (after around day 240) in which discharge in channels is higher than in the distributed system. Under the future melt scenario $\mathrm{2011}×\mathrm{4},$ the integrated discharge in the channelized system increases to 77 % of the integrated discharge in the distributed system. Early in the melt season, discharge increases in both the channelized system and distributed system simultaneously. Similar to the other melt seasons, discharge in the distributed system increases at a faster rate. However, drainage in the channelized system equals or exceeds that of the distributed system much earlier in the year (day 200).

4 Discussion

## 4.1 Model fit

The modelled velocities are the combined result of five models (RACMO2.3, SRLF, an ice sheet model, the associated adjoint model, and a subglacial hydrology model) and several datasets. Most parameters used in the models are assigned standard values, with calibrated parameter values for the subglacial hydrology model. The validation simulation affirms the calibrated model and theory, as measured velocities are reproduced to the same qualitative level of fit. Although each model has biases and is limited by assumptions, their combined result reproduces measured ice velocities to a first order. Many of the features observed in the GPS time series are captured in the modelled velocities. This gives confidence that the models and datasets are representative of their respective component. The complexity of the models, and the process, makes assigning a model uncertainty infeasible. It is unclear how to partition the cause of model mismatch between errors in inputs such as topography, model–theoretical uncertainty such as the form of the sliding law, or computationally imposed limitations such as grid resolution or choice of ice sheet model.

Overall, model velocities are observed to be better at mid-elevation than at either the lowest or highest sites. Model velocities at sites S1–S2 are likely affected by the model not recreating the subglacial water routing inferred by . A number of factors could contribute to differences in water routing, including errors in topographic data, the spatial distribution of inputs from crevasse fields, and model assumptions and boundary conditions. In general, thin ice and steep gradients in topography make ice flow and hydrology modelling near the margin difficult. Steep surface gradients may lead to stresses assumed negligible by the hybrid formulation in the ice stress balance. Drainage components in the subglacial hydrology model are formulated in terms of effective pressure, on the implicit assumption that they remain full. Underneath thin ice, or when there are steep gradients, both channels and cavities could be expected to exist while partially full or empty. The atmospheric pressure boundary condition prescribed at the ice sheet margin in the subglacial hydrology model may, in reality, extend inland for periods in the summer. Additionally, the high-velocity spring–early summer events observed in the GPS records occur before any melt is predicted by RACMO2.3. Similar to the study by , modelled velocities do not capture this behaviour. These velocity events may be the result of internal dynamics of water stored over winter , such as flooding events, that the subglacial hydrology model does not capture or early season melt which RACMO2.3 does not predict.

Modelled velocities at sites S6–S7 may be affected by excess capacity in the cavity system due to overprediction of basal ice velocities from the inversion process. The inversion process results in a sliding ratio of approximately 0.8 at the high elevations . However, internal deformation can be expected to be dominate over basal sliding so far inland, suggesting a much lower sliding ratio. Measurements at boreholes in the Paakitsoq region at lower elevations show a sliding ratio of 0.44–0.73 during the winter, increasing episodically to 0.9 during the summer . The largest discrepancy between ablation at a weather stations and RACMO2.3 modelled surface runoff occurs at site S6, likely due to RACMO2.3, allowing for refreezing of surface melt. This additional complexity increases the uncertainty in runoff predictions, and surface input to the base may be underestimated at sites S6 and S7.

The development of the subglacial hydrology system is driven by surface runoff input to the bed and is a key control on velocities across the domain. Measurements in the Russell Glacier area of a single 63.1 km2 moulin terminating catchment at approximately 1250 m elevation over a 72 h period found that RACMO2.3 overestimated runoff by approximately 60 % . In general, an average of multiple regional climate models was reported to overpredict surface runoff by +21 to 58 % for the catchment . The impact of any discrepancies between modelled and actual surface runoff is not necessarily limited to a local temporal and/or local spatial scale. However, to what extent the results of the study generalize across the model domain is unresolved.

Spatial maps of modelled velocities show some numerical artifacts. Although these do not appear to have a strong direct impact on the velocities at the GPS stations, numerical artifacts are a cause for concern and should be mitigated in future work. One likely cause is high-velocity gradients near the lateral margins due to the Dirichlet boundary conditions. A second is strong variations in basal drag due to subglacial hydrology likely results in non-negligible horizontal gradients in vertical velocities, contrary to the assumptions of the hybrid formulation. Alternative boundary conditions or numerical schemes to improve convergence may mitigate these effects but were not pursued further in this study.

## 4.2 Model sensitivity

Model velocities calculated with the two different sliding laws are comparable during much of the melt season. The timing of events is not affected by the choice of sliding law, and the primary difference observed is the magnitude of velocities during short-term speedup events. The overprediction of speedup during events with the Schoof sliding law suggests adding a regularization constant, such that a minimum basal drag exists. Such a term could reflect the fact that the subglacial hydrological system may not extend throughout a gridcell or that part of the cell has a weakly connected system with a different water pressure . Simulation results show the Budd sliding law with standard exponent values has practical value in simulations. However, the form and parameters of the sliding law remain uncertain, and the Schoof law has greater theoretical support (Hewitt2013).

Calibrating the integrated model is an underdetermined problem, as the number of observations is not sufficient to constrain the parameters in all the models. The calibration therefore focuses on the key parameters of the subglacial hydrology model, while keeping parameters of the ice sheet model and surface hydrology model constant. The calibration was achieved mainly by trial and error, starting with values used in . Most model parameters of the integrated model are similar to previous studies applying the subglacial hydrology model . The most significant parameter value difference is the sheet flux coefficient (Ks). The primary value of 10−2${\mathrm{Pa}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ in our study is greater than the value of 10−5${\mathrm{Pa}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ used in . The parameter values for the model reported in , which are calibrated against observed water discharge at an outlet in the Paakitsoq region, were found not to reproduce GPS velocity records, as water at mid–high elevations was not effectively evacuated. The difference in parameters suggests that care needs to be taken transferring parameter values between study sites in different areas and at different scales.

The calibrated value for sheet conductivity is at the higher end of inferred values for till . Although model results are no longer comparable when sheet conductivity decreases by an order of magnitude, model results are resilient to heterogeneity. The simple tests conducted suggest that random heterogeneity in sheet conductivity has a lower impact than larger-scale spatial patterns. Heterogeneity in sheet conductivity could arise from local topography, variable till coverage, and till properties (including deformational history). A constant bed roughness height scale of 0.5 m is used in this paper. However, patterns of bed roughness would also provide a strong control on discharge at the base. Overall, model results suggest it is necessary for the distributed system to be able to sustain a high discharge.

The initial rapid delivery of a large volume of water to the bed during individual lake hydrofracture events is not observed to have a lasting or widespread effect on modelled velocities (in line with observations by ). This suggests that lake hydrofracture events themselves are not a key process in the long-term or large-scale development of the subglacial hydrological system, as at lower elevations the numerous conduits and high water input drive channelization, while at higher elevations a combination of insufficient input and conditions unfavourable for channelization exists. Rather, the primary impact of lake hydrofracture is in opening surface-to-bed connections. The spatial density of such events has been shown to affect the rate of development of channelized drainage and to act as a key mechanism for the creation of moulins away from crevasse fields or current lake basins , which then drain a significant proportion of the overall surface melt as we find in this study and previous work .

The configuration of internal catchments and moulins which drain crevasses was not found to have a strong impact; neither was eliminating drainage of water generated from ablation in internal catchments within crevasse fields. However, the GPS sites at which model velocities are compared do not capture spatial heterogeneity of crevasse drainage, which occurs along the length of the ice margin. Hence, the impact may be much stronger at other locations within the study area. However, since model velocities at higher GPS sites were not observed to vary with changes in crevasse drainage, the impact of crevasse drainage should be limited to the margin.

## 4.3 Model complexity

It is encouraging that the results provide a clear match to many features seen in velocity observations, particularly at the relatively coarse resolution used. However, the models and workflow applied in this paper are characterized by a high degree of complexity. An important consideration is where simplifications can be applied and where further complexity may be justified.

The use of a higher-order ice sheet model/inversion code should be explored due to increased accuracy in basal velocity calculations. Basal velocities are a key control of the subglacial hydrological system since they determine cavity spacing and provide important feedback . A higher-order model may perform more robustly throughout the study area. Areas where the performance of the hybrid model may be expected to be suboptimal occur throughout the study area. Such areas are characterized by high aspect ratio, high variability in basal topography, or low sliding ratio. The ice flow model is also constrained by the assumption of a uniform temperature distribution throughout the ice. Calculating a thermal–mechanical steady state, or alternatively inverting for the structure, would increase accuracy of calculated basal velocities. Either of these options could be incorporated in the step with the linear inversion at limited cost since this step is only executed once. Importantly, both the use of a higher-order ice sheet model and determination of the thermal state can be implemented without adding further assumptions or unconstrained parameters.

The subglacial hydrology model is the least constrained model in the workflow. Many parameters remain unknown and the exploration of its behaviour is limited by the parameter space searched. However, a key behaviour not replicated is the late summer–fall slowdown and the subsequent gradual winter acceleration. The integrated model returns to its initial state at the end of summer. This indicates a need for a component of the model which operates on a longer timescale than is currently included. The difficulty in recreating both the smoother 2009 velocity record and the more variable 2011 record also suggests interannual variability in the background state of the hydrological system. A model component simulating weakly connected regions of the hydrological system as incorporated in may be key to reproducing these observations. These regions are conceptualized as parts of the distributed system with a much lower hydraulic connectivity. The connectivity of these regions may be temporally variable.

The SRLF model offers the best opportunity for simplification. To at least a first order, lakes which hydrofracture can be modelled as moulins (in line with observations by ). This suggests using the locations of moulins derived from satellite imagery acquired at the end of the melt season as representative of both moulins outside of lake basins and hydrofractured lakes. Lake hydrofracture events are observed to result in temporarily faster local flow . In order for a model to capture these events, the specific location, timing, and volume of lakes will need to incorporated into the model. Given the ongoing uncertainties around the processes controlling hydrofracture, this suggests that using observational records of lake drainages derived from satellite imagery (as in ) to derive hydrofracture input to the ice–bed interface forms a valid strategy for present-day studies, though such an approach would not work for prognostic tests. Crevasses also drain a significant proportion of water, most of which travels over the ice surface into crevasse fields from upstream rather than being generated locally. The controls on water drainage through crevasses to the ice sheet bed are poorly understood, but they may have an important role as the spatial density of water inputs are known to influence the development of the subglacial hydrological system . Since moulins and crevasses drain water in a continuous manner with a high spatial density, a simpler surface hydrology scheme approximating input into each drainage pathway from its local catchment may be effective. The output of each catchment into the corresponding drainage pathway may be simplified to two output hydrographs, one for snow-covered and the other for bare-ice conditions. For internal catchments of crevasse fields routing can likely be neglected. This calculation needs only be done once; moulin input at each time step could then be calculated at little computational cost based on total surface runoff and the dominant surface cover in the catchment.

## 4.4 Implications

The existence of channelized and distributed systems beneath the GrIS is inferred indirectly through borehole observations, dye tracing experiments, and patterns of GPS velocities, building on extensive observations and theoretical developments derived from studies on alpine glaciers. The key result of this paper is to provide numerical support for the understanding of subglacial hydrology of the GrIS, based on theories derived from studies of alpine glaciers, as well as support for the explicit description of the model components we include (i.e. the equations used). We show that these theories can quantitatively reproduce measurements to a first order and, in the sense of our validation, predict ice velocities. This builds on previous work which shows that this understanding can be used to reproduce idealized seasonal patterns of ice velocity as well as effective pressures in line with ice velocities .

The timing of velocity variations is controlled by surface input and modulated by subglacial hydrology. At high elevations where channelization is not observed, variations in model velocities track modelled surface runoff closely. GPS velocities, however, do not show the same fidelity to the time series of ablation from automatic weather stations, which are qualitatively more variable than modelled runoff. This suggests dampening of the variability of surface input by the supraglacial and subglacial hydrology and that variability in daily ablation rates are not simply correlated to faster flow. A quantitative analysis of the two time series may provide better insight into the relationship between surface melt and ice velocities. However, ice velocities are driven by the cumulative melt over a larger upstream area from the point of measurement, which may not be well represented by the variability of melt at a single point. At lower elevations, channelization is important in modulating the impact of surface water on ice velocities. The low modelled and observed velocities closer to the ice margin imply a consistently high effective pressure at the GPS sites due to the impact of channelization on water pressures and water routing.

Modelling predicts that average summer ice velocities over the melt season will increase with melt season intensity. A similar correlation was observed in GPS records over the upper ablation zone of the Russell Glacier region by , but not in GPS records at North Lake, western Greenland, by . This implies that more intense melt seasons will result in a higher ice flux towards the margin during the summer. Whether this would be compensated for in terms of the average annual ice velocity by decreased ice flux during the winter is unresolved by the model.

Channelization is observed to develop more extensively and further inland as melt intensity increases. This trend is observed in the three modelled melt seasons and continues into the two future melt scenarios. This suggest that the subglacial hydrological system will continue to drain surface meltwater input in a similar manner as melt intensity increases beyond 2012 levels. Since channelization is thought to result in the observed slowdown in mid–late summer and is also postulated to result in slowdown in the subsequent winter and spring , model results suggest increasing summer melt intensity could lead to a more spatially extensive annual velocity slowdown. The slowdown may also become more pronounced in the future as the channelized system is predicted to drain an increased proportion of water and accesses a larger proportion of the model domain.

However, as channelization increases up-ice in our model, we do not see a marked impact on model velocities. Model velocities at the higher GPS stations in model runs 2011×2 and 2011×4 both show a similar pattern to 2011, with a higher magnitude of ice flow during speedup events. We do not observe a shift in velocity patterns towards that of lower GPS stations, with acceleration early in the melt season transitioning to deceleration in the latter part of the melt season. This suggests that channelization may have a more limited impact on annual velocities in the accumulation zone. The magnitude of any impact is unresolved by our model. In particular, although there are periods when velocities in the 2011×4 run are lower than 2011×2 and 2011, the magnitude of this decrease is bounded. This is a limitation of the model, which is only able to decrease velocities nominally below the initial winter values. This also implies that we are unable to model the winter season accurately.

The key question for the longer-term response of the ice sheet to increased melt is whether the potential summer increase in velocity due to increased melt will outweigh any late summer and winter decrease due to the evolution of a more efficient system under higher melt conditions. While observations show long-term decreases in ice velocities in the lower ablation zone , this question remains unresolved at higher elevations. Although we cannot directly predict annual velocities with the model presented in this study, we can investigate how annual velocities may change at the limits of winter behaviour.

One limit of winter velocities is that integrated ice flow over the winter decreases faster than integrated ice flow over the summer. This is observed in GPS measurements in the upper ablation zone in the Russell Glacier region by . Under this limit, channelization has a similar impact at high elevations as in the ablation zone. Velocity measurements near the vicinity of a lake hydrofracture at approximately 1450 m elevation suggest that channelization occurs even at high elevations . Winter velocities in this limit will decrease until they hit a lower bound where flow is purely deformational, with no contribution from basal sliding. The maximum increase in mean summer velocities is approximately 60 m yr−1, at GPS stations 4 and 5 between the 2009 and 2011×4 melt scenarios. Assuming a winter velocity of 100 m yr−1 and an 7-month winter, the summer increase predicted by the model would compensate for a possible reduction in winter velocity to around 60 m yr−1. This approaches the lower bound for winter velocity suggested by borehole measurements showing that internal deformation accounts for 25–50 % of the total ice velocity in the Paakitsoq region of western Greenland . Climate model predictions suggest surface runoff rates quadrupling from present levels by circa 2100 .

The second limit occurs if winter velocities at higher elevations are not impacted by channelization and summer velocities dominate the annual signal. The argument that thick ice and shallow surface slopes inhibit channel growth at high elevations favour this limit of behaviour , as do observations suggesting limited changes in the efficiency of the channelized system . Under this scenario, a change of mean summer velocity of 110 to 170 m yr−1 with winter velocities remaining constant at 100 m yr−1 would result in mean annual velocities increasing from 104 to 129 m yr−1 between the 2009 and 2011×4 simulations. Under this scenario, annual velocities would increase by approximately 25 % by circa 2100, when surface runoff is predicted to quadruple.

Interpreting the model velocity output from the future melt scenarios is difficult, however, and our bounds should be interpreted cautiously. We do not evolve our model geometry or evolve the distribution of surface drainage locations, nor do we use climate model predictions of surface runoff for the future. As melt season intensity increases, the validity of the initialization and calibration parameters also becomes more uncertain. Further, the model has bias towards capturing short-term speedup events rather than prolonged slowdowns due to model velocities remaining near or above their winter values. The modelled velocities show higher variability and a significant increase in the magnitude of short-term speedup events. However, quantifying whether the impact of these events on annual velocity will be compensated for by a corresponding late summer slowdown or by a winter slowdown is beyond the capability of the current model. Model output can be interpreted to suggest that a late summer velocity slowdown compensating for an early summer speedup is less likely at higher elevations. It is not evident, however, whether the suggested upper bound of a 25 % increase in annual velocities in this limit would have an impact on the overall mass budget of the ice sheet as great as that from a increase in surface runoff in itself.

The success of the model in recreating features in the measured velocities provides validation for each model component, as well as their integration. The work supports integrating models of high complexity that incorporate a range of processes. Further model refinement and data acquisition should continue to improve the fit between modelled and measured velocities. A key uncertainty in the initialization process was the subglacial hydrology model run during winter and the subsequent inversion for background basal parameters. Although the process used cannot capture year-to-year changes, the practical value of the initialization process is implicitly validated through the subsequent fit to measured velocities.

5 Conclusions

In this paper, we couple multiple models in order to predict summer ice velocities at the southwest margin of GrIS from topographic and climatic input data. These models represent the main components of the ice sheet system: supraglacial hydrology, subglacial hydrology, and ice flow. The key component of the simulations presented in this paper is a coupled hydrology–ice flow model. This integrated model is initialized using a workflow incorporating the adjoint ice flow model and is forced during the simulations using surface input from a surface hydrology model. Calibration of the integrated model takes advantage of GPS velocities from two summer melt seasons: 2009 and 2011. The model validation on 2012 GPS data reproduces measured ice velocities to a similar degree as in 2009 and 2011. To a first order, the magnitude and timing of the measured velocities are replicated in modelled velocities at multiple sites.

The success of the multicomponent modelling to recreate summer velocities reflects on the integrity of each individual model and dataset. This work should encourage further model coupling as it suggests that individual components and datasets are robust. However, limitations of the multicomponent model are evident in the model output, particularly that the model velocity does not significantly drop below its initialized winter value. Additional data and theory will be necessary to address these issues. Together, the models also form a quantitative test of the hypothesis proposed by numerous authors (Chandler et al.2013; Colgan et al.2011; Cowton et al.2013; Hewitt2013; Hoffman et al.2011; Nienow et al.2017; Schoof2010; van de Wal et al.2015) that the summer acceleration of the GrIS margin is controlled by the evolution of the subglacial hydrological system in a manner analogous to the seasonal speedup of alpine glaciers. The key result of this paper is quantitative support in favour of this hypothesis.

The observed decadal-timescale slowdown at the margin of the GrIS is attributed to increased channelization reducing late summer and winter water pressures , and hence velocities. Our results suggest that the decadal slowdown will continue in the near future, particularly close to the ice margin. However, the model predicts a strong scaling of the average summer velocity with melt season intensity. We investigate the impact of this under two limits. If integrated ice flow over the winter decreases faster than integrated ice flow over the summer at higher elevations, our modelling suggests that annual velocities in the upper ablation zone would begin to increase by around 2100 (when surface runoff is predicted to quadruple from present levels), as predicted summer velocity increases offset likely winter velocity decreases. In the second limit, in which winter velocities remain at present levels while summer velocities increase, our model suggests an upper bound of a 25 % increase in annual velocities by around 2100.

Data availability
Data availability.

All datasets used are publicly available. BedMachine data are available from Morlighem et al. (2015). MEaSUREs Greenland Ice Sheet Velocity Map data are available from Joughin et al. (2010b). Moulin location data are available from Yang and Smith (2016). RACMO2.3 data are available on request from Michiel R. van den Broeke and Brice Noël. GPS data are available via the NERC Polar Data Centre (Tedstone and Nienow, 2017). Code is currently not available.

Supplement
Supplement.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

We would like to thank Mathieu Morlighem and Ian Joughin for the BedMachine2 and MEaSUREs datasets and Michiel R. van den Broeke and Brice Noel for providing RACMO2.3 data. CPK would like to acknowledge Robert Arthern for guidance on writing an ice sheet model and inversion code and Brent Minchew, Poul Christoffersen, and Daniel Goldberg for thoughtful discussions. Both authors thank Ian Hewitt for sharing model code, the two referees for their careful and constructive reviews, as well as the editor Andreas Vieli. Conrad P. Koziol was funded through St. John's College, Cambridge, and in part by UK Natural Environment Research Council Grant NE/M003590/1.

Edited by: Andreas Vieli
Reviewed by: two anonymous referees

References

Andrews, L. C., Catania, G. A., Hoffman, M. J., Gulley, J. D., Lüthi, M. P., Ryser, C., Hawley, R. L., and Neumann, T. A.: Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet, Nature, 514, 80–83, https://doi.org/10.1038/nature13796, 2014. a

Arnold, N., Richards, K., Willis, I., and Sharp, M.: Initial results from a distributed, physically based model of glacier hydrology, Hydrol. Process., 12, 191–219, https://doi.org/10.1002/(SICI)1099-1085(199802)12:2<191::AID-HYP571>3.0.CO;2-C, 1998. a, b

Arthern, R. J., Hindmarsh, R. C. A., and Williams, C. R.: Flow speed within the Antarctic ice sheet and its controls inferred from satellite observations, J. Geophys. Res.-Earth, 120, 1171–1188, https://doi.org/10.1002/2014JF003239, 2015. a

Banwell, A., Hewitt, I., Willis, I., and Arnold, N.: Moulin density controls drainage development beneath the Greenland Ice Sheet, J. Geophys. Res.-Earth, 121, 2248–2269, https://doi.org/10.1002/2015JF003801, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Bartholomaus, T., Anderson, R., and Anderson, S.: Growth and collapse of the distributed subglacial hydrologic system of Kennicott Glacier, Alaska, USA, and its effects on basal motion, J. Glaciol., 57, 985–1002, https://doi.org/10.3189/002214311798843269, 2011. a

Bartholomew, I., Nienow, P., Mair, D., Hubbard, A., King, M. A., and Sole, A.: Seasonal evolution of subglacial drainage and acceleration in a Greenland outlet glacier, Nat. Geosci., 3, 408–411, https://doi.org/10.1038/ngeo863, 2010. a

Bartholomew, I., Nienow, P., Sole, A., Mair, D., Cowton, T., Palmer, S., and Wadham, J.: Supraglacial forcing of subglacial drainage in the ablation zone of the Greenland ice sheet, Geophys. Res. Lett., 38, L08502, https://doi.org/10.1029/2011GL047063, 2011. a, b

Bartholomew, I., Nienow, P., Sole, A., Mair, D., Cowton, T., and King, M. A.: Short-term variability in Greenland Ice Sheet motion forced by time-varying meltwater drainage: Implications for the relationship between subglacial drainage system behavior and ice velocity, J. Geophys. Res., 117, F03002, https://doi.org/10.1029/2011JF002220, 2012. a

Bougamont, M., Christoffersen, P., Hubbard, A. L., Fitzpatrick, A., Doyle, S., and Carter, S.: Sensitive response of the Greenland Ice Sheet to surface melt drainage over a soft bed, Nat. Commun., 5, 5052, https://doi.org/10.1038/ncomms6052, 2014. a, b, c, d

Budd, W., Keage, P., and Blundy, N.: Empirical studies of ice sliding, J. Glaciol., 23, 157–170, 1979. a

Chandler, D. M., Wadham, J. L., Lis, G. P., Cowton, T., Sole, A., Bartholomew, I., Telling, J., Nienow, P., Bagshaw, E. B., Mair, D., Vinen, S., and Hubbard, A.: Evolution of the subglacial drainage system beneath the Greenland Ice Sheet revealed by tracers, Nat. Geosci., 6, 195–198, https://doi.org/10.1038/ngeo1737, 2013. a, b, c, d, e, f, g, h, i

Clason, C. C., Mair, D. W. F., Nienow, P. W., Bartholomew, I. D., Sole, A., Palmer, S., and Schwanghart, W.: Modelling the transfer of supraglacial meltwater to the bed of Leverett Glacier, Southwest Greenland, The Cryosphere, 9, 123–138, https://doi.org/10.5194/tc-9-123-2015, 2015. a

Colgan, W., Rajaram, H., Anderson, R., Steffen, K., Phillips, T., Joughin, I., Zwally, H. J., and Abdalati, W.: The annual glaciohydrology cycle in the ablation zone of the Greenland ice sheet: Part 1. Hydrology model, J. Glaciol., 57, 697–709, 2011. a

Colgan, W., Rajaram, H., Anderson, R., Steffen, K., Zwally, J., Phillips, T., and Abdalati, W.: The annual glaciohydrology cycle in the ablation zone of the Greenland ice sheet: Part 2. Observed and modeled ice flow, J. Glaciol., 58, 51–64, https://doi.org/10.3189/2012JoG11J081, 2012. a, b, c

Cooley, S. W. and Christoffersen, P.: Observation bias correction reveals more rapidly draining lakes on the Greenland Ice Sheet, J. Geophys. Res.-Earth, 122, 1867–1881, 2017. a

Cowton, T., Nienow, P., Sole, A., Wadham, J., Lis, G., Bartholomew, I., Mair, D., and Chandler, D.: Evolution of drainage system morphology at a land-terminating Greenlandic outlet glacier, J. Geophys. Res.-Earth, 118, 29–41, https://doi.org/10.1029/2012JF002540, 2013. a, b

Das, S. B., Joughin, I., Behn, M. D., Howat, I. M., King, M. A., Lizarralde, D., and Bhatia, M. P.: Fracture propagation to the base of the Greenland Ice Sheet during supraglacial lake drainage, Science, 320, 778–781, https://doi.org/10.1126/science.1153360, 2008. a

de Fleurian, B., Morlighem, M., Seroussi, H., Rignot, E., van den Broeke, M. R., Munneke, P. K., Mouginot, J., Smeets, P. C. J. P., and Tedstone, A. J.: A modeling study of the effect of runoff variability on the effective pressure beneath Russell Glacier, West Greenland, J. Geophys. Res.-Earth, 121, 1834–1848, https://doi.org/10.1002/2016JF003842, 2016. a, b, c, d

Dow, C. F., Kulessa, B., Rutt, I., Doyle, S., and Hubbard, A.: Upper bounds on subglacial channel development for interior regions of the Greenland ice sheet, J. Glaciol., 60, 1044–1052, https://doi.org/10.3189/2014JoG14J093, 2014. a

Doyle, S. H., Hubbard, A., Fitzpatrick, A. A. W., van As, D., Mikkelsen, A. B., Pettersson, R., and Hubbard, B.: Persistent flow acceleration within the interior of the Greenland ice sheet, Geophys. Res. Lett., 41, 899–905, https://doi.org/10.1002/2013GL058933, 2014. a

Fitzpatrick, A. A. W., Hubbard, A., Joughin, I., Quincey, D. J., Van As, D., Mikkelsen, A. P., Doyle, S. H., Hasholt, B., and Jones, G. A.: Ice flow dynamics and surface meltwater flux at a land-terminating sector of the Greenland ice sheet, J. Glaciol., 59, 687–696, https://doi.org/10.3189/2013JoG12J143, 2013. a

Flowers, G. E.: Modelling water flow under glaciers and ice sheets, P. Roy. Soc. A-Math. Phy., 471, 20140907, https://doi.org/10.1098/rspa.2014.0907, 2015. a

Flowers, G. E. and Clarke, G. K. C.: A multicomponent coupled model of glacier hydrology. 2. Application to Trapridge Glacier, Yukon, Canada, J. Geophys. Res., 107, 2288, https://doi.org/10.1029/2001JB001124, 2002. a

Fountain, A. G. and Walder, J. S.: Water flow through temperate glaciers, Rev. Geophys., 36, 299–328, https://doi.org/10.1029/97rg03579, 1998. a, b

Gagliardini, O., Cohen, D., Raback, P., and Zwinger, T.: Finite-element modeling of subglacial cavities and related friction law, J. Geophys. Res.-Earth, 112, f02027, https://doi.org/10.1029/2006JF000576, 2007. a

Goldberg, D. N.: A variationally derived, depth-integrated approximation to a higher-order glaciological flow model, J. Glaciol., 57, 157–170, https://doi.org/10.3189/002214311795306763, 2011. a

Goldberg, D. N. and Heimbach, P.: Parameter and state estimation with a time-dependent adjoint marine ice sheet model, The Cryosphere, 7, 1659–1678, https://doi.org/10.5194/tc-7-1659-2013, 2013. a

Hanna, E., Navarro, F. J., Pattyn, F., Domingues, C. M., Fettweis, X., Ivins, E. R., Nicholls, R. J., Ritz, C., Smith, B., Tulaczyk, S., Whitehouse, P. L., and Zwally, H. J.: Ice-sheet mass balance and climate change, Nature, 498, 51–59, https://doi.org/10.1038/nature12238, 2013. a

Heimbach, P. and Bugnion, V.: Greenland ice-sheet volume sensitivity to basal, surface and initial conditions derived from an adjoint model, Ann. Glaciol., 50, 67–80, https://doi.org/10.3189/172756409789624256, 2009. a, b

Hewitt, I.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sc. Lett., 371–372, 16–25, https://doi.org/10.1016/j.epsl.2013.04.022, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u

Hoffman, M. and Price, S.: Feedbacks between coupled subglacial hydrology and glacier dynamics, J. Geophys. Res.-Earth, 119, 414–436, https://doi.org/10.1002/2013JF002943, 2014. a, b, c

Hoffman, M. J., Catania, G. A., Neumann, T. A., Andrews, L. C., and Rumrill, J. A.: Links between acceleration, melting, and supraglacial lake drainage of the western Greenland Ice Sheet, J. Geophys. Res., 116, F04035, https://doi.org/10.1029/2010JF001934, 2011. a, b, c

Hoffman, M. J., Andrews, L. C., Price, S. A., Catania, G. A., Neumann, T. A., Luethi, M. P., Gulley, J., Ryser, C., Hawley, R. L., and Morriss, B. F.: Greenland subglacial drainage evolution regulated by weakly-connected regions of the bed, Nat. Commun., 7, 13903, https://doi.org/10.1038/ncomms13903, 2016. a, b, c

Hoffman, M. J., Perego, M., Andrews, L. C., Price, S. F., Neumann, T. A., Johnson, J. V., Catania, G., and Lüthi, M. P.: Widespread Moulin Formation During Supraglacial Lake Drainages in Greenland, Geophys. Res. Lett., 45, 778–788, https://doi.org/10.1002/2017GL075659, 2018. a

Howat, I., Negrete, A., and Smith, B.: MEaSURES Greenland Ice Mapping Project (GIMP) Digital Elevation Model, Version 1. 90 m resolution, NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado USA, https://doi.org/10.5067/NV34YUIXLP9W, 2015. a

Howat, I. M., Box, J. E., Ahn, Y., Herrington, A., and McFadden, E. M.: Seasonal variability in the dynamics of marine-terminating outlet glaciers in Greenland, J. Glaciol., 56, 601–613, https://doi.org/10.3189/002214310793146232, 2010. a

Joughin, I., Smith, B. E., Howat, I. M., Scambos, T., and Moon, T.: Greenland flow variability from ice-sheet-wide velocity mapping, J. Glaciol., 56, 415–430, https://doi.org/10.3189/002214310792447734, 2010a.  a

Joughin, I., Smith, B. E., Howat, I. M., Scambos, T., and Moon, T.: MEaSUREs Greenland Ice Velocity Map from InSAR Data, NASA DAAC at the National Snow and Ice Data Center, Boulder, Colorado, USA, https://doi.org/10.5067/MEASURES/CRYOSPHERE/nsidc-0478.001, 2010b. a

Koziol, C., Arnold, N., Pope, A., and Colgan, W.: Quantifying supraglacial meltwater pathways in the Paakitsoq region, West Greenland, J. Glaciol., 63, 464–476, 2017. a, b, c, d, e, f

Koziol, C. P.: Modelling the impact of surface melt on the hydrology and dynamics of the Greenland Ice Sheet, PhD thesis, University of Cambridge, Cambridge, UK, 2017. a

Koziol, C. P. and Arnold, N.: Incorporating modelled subglacial hydrology into inversions for basal drag, The Cryosphere, 11, 2783–2797, https://doi.org/10.5194/tc-11-2783-2017, 2017. a, b, c, d, e, f, g

Leeson, A., Shepherd, A., Briggs, K., Howat, I., Fettweis, X., Morlighem, M., and Rignot, E.: Supraglacial lakes on the Greenland ice sheet advance inland under warming climate, Nat. Clim. Change, 5, 51–55, 2015. a

MacAyeal, D.: A tutorial on the use of control methods in ice sheet modeling, J. Glaciol, 39, 91–98, 1993. a

Martin, N. and Monnier, J.: Adjoint accuracy for the full Stokes ice flow model: limits to the transmission of basal friction variability to the surface, The Cryosphere, 8, 721–741, https://doi.org/10.5194/tc-8-721-2014, 2014. a

Meierbachtol, T., Harper, J., and Humphrey, N.: Basal drainage system response to increasing surface melt on the Greenland Ice Sheet, Science, 341, 777–779, https://doi.org/10.1126/science.1235905, 2013. a

Moon, T., Joughin, I., Smith, B., van den Broeke, M. R., van de Berg, W. J., Noel, B., and Usher, M.: Distinct patterns of seasonal Greenland glacier velocity, Geophys. Res. Lett., 41, 7209–7216, https://doi.org/10.1002/2014GL061836, 2014. a

Morlighem, M., Seroussi, H., Larour, E., and Rignot, E.: Inversion of basal friction in Antarctica using exact and incomplete adjoints of a higher-order model, J. Geophys. Res.-Earth, 118, 1746–1753, https://doi.org/10.1002/jgrf.20125, 2013. a

Morlighem, M., Rignot, E., Mouginot, J., Seroussi, H., and Larour, E.: Deeply incised submarine glacial valleys beneath the Greenland ice sheet, Nat. Geosci., 7, 418–422, https://doi.org/10.1038/ngeo2167, 2014. a

Morlighem, M., Rignot, E., Mouginot, J., Seroussi, H., and Larour, E.: IceBridge BedMachine Greenland, Version 2, Thickness, NASA DAAC at the National Snow and Ice Data Center, Boulder, Colorado USA, https://doi.org/10.5067/AD7B0HQNSJ29, 2015. a, b

Nienow, P., Sole, A., Slater, D., and Cowton, T.: Recent Advances in Our Understanding of the Role of Meltwater in the Greenland Ice Sheet System, Current Climate Change Reports, 3, 330–344, 2017. a, b

Noël, B., van de Berg, W. J., van Meijgaard, E., Kuipers Munneke, P., van de Wal, R. S. W., and van den Broeke, M. R.: Evaluation of the updated regional climate model RACMO2.3: summer snowfall impact on the Greenland Ice Sheet, The Cryosphere, 9, 1831–1844, https://doi.org/10.5194/tc-9-1831-2015, 2015. a, b

Pimentel, S. and Flowers, G. E.: A numerical study of hydrologically driven glacier dynamics and subglacial flooding, P. Roy. Soc. A-Math. Phy., 467, 537–558, https://doi.org/10.1098/rspa.2010.0211, 2010.  a, b, c

Poinar, K., Joughin, I., Das, S. B., Behn, M. D., Lenaerts, J. T. M., and Broeke, M. R.: Limits to future expansion of surface-melt-enhanced ice flow into the interior of western Greenland, Geophys. Res. Lett., 42, 1800–1807, https://doi.org/10.1002/2015GL063192, 2015. a

Raymond, C. and Nolan, M.: Drainage of a glacial lake through an ice spillway, Intl. Assoc. Hydrol. Sci. Publ., 264, 199–210, available at: http://iahs.info/uploads/dms/iahs_264_0199.pdf (last access: May 2017), 2000. a

Ryser, C., Luthi, M. P., Andrews, L. C., Hoffman, M. J., Catania, G. A., Hawley, R. L., Neumann, T. A., and Kristensen, S. S.: Sustained high basal motion of the Greenland ice sheet revealed by borehole deformation, J. Glaciol., 60, 647–660, https://doi.org/10.3189/2014JoG13J196, 2014. a, b

Schoof, C.: The effect of cavitation on glacier sliding, P. Roy. Soc. A-Math. Phy., 461, 609–627, https://doi.org/10.1098/rspa.2004.1350, 2005. a

Schoof, C.: Ice-sheet acceleration driven by melt supply variability, Nature, 468, 803–806, https://doi.org/10.1038/nature09618, 2010. a, b, c

Schoof, C., Rada, C. A., Wilson, N. J., Flowers, G. E., and Haseloff, M.: Oscillatory subglacial drainage in the absence of surface melt, The Cryosphere, 8, 959–976, https://doi.org/10.5194/tc-8-959-2014, 2014. a

Shannon, S. R., Payne, A. J., Bartholomew, I. D., van den Broeke, M. R., Edwards, T. L., Fettweis, X., Gagliardini, O., Gillet-Chaulet, F., Goelzer, H., Hoffman, M. J., Huybrechts, P., Mair, D. W. F., Nienow, P. W., Perego, M., Price, S. F., Smeets, C. J. P. P., Sole, A. J., van de Wal, R. S. W., and Zwinger, T.: Enhanced basal lubrication and the contribution of the Greenland ice sheet to future sea-level rise, P. Natl. Acad. Sci. USA, 110, 14156–14161, https://doi.org/10.1073/pnas.1212647110, 2013. a

Smith, L. C., Chu, V. W., Yang, K., Gleason, C. J., Pitcher, L. H., Rennermalm, A. K., Legleiter, C. J., Behar, A. E., Overstreet, B. T., Moustafa, S. E., Tedesco, M., Forster, R. R., LeWinter, A. L., Finnegan, D. C., Sheng, Y., and Balog, J.: Efficient meltwater drainage through supraglacial streams and rivers on the southwest Greenland ice sheet, P. Natl. Acad. Sci. USA, 2000, 201413024, https://doi.org/10.1073/pnas.1413024112, 2015. a

Smith, L. C., Yang, K., Pitcher, L. H., Overstreet, B. T., Chu, V. W., Rennermalm, Å. K., Ryan, J. C., Cooper, M. G., Gleason, C. J., Tedesco, M., Jeyaratnam, J., van As, D., van den Broeke, M. R., van de Berg, W. J., Noël, B., Langen, P. L., Cullather, R. I., Zhao, B., Willis, M. J., Hubbard, A., Box,, J. E., Jenner, B. A., and Behar, A. E.: Direct measurements of meltwater runoff on the Greenland ice sheet surface, P. Natl. Acad. Sci. USA, 114, E10622–E10631, 2017. a, b, c

Sole, A., Nienow, P., Bartholomew, I., Mair, D., Cowton, T., Tedstone, A., and King, M. A.: Winter motion mediates dynamic response of the Greenland Ice Sheet to warmer summers, Geophys. Res. Lett., 40, 3940–3944, https://doi.org/10.1002/grl.50764, 2013. a, b

Sole, A. J., Mair, D. W. F., Nienow, P. W., Bartholomew, I. D., King, M. A., Burke, M. J., and Joughin, I.: Seasonal speedup of a Greenland marine-terminating outlet glacier forced by surface melt-induced changes in subglacial hydrology, J. Geophys. Res., 116, F03014, https://doi.org/doi:10.1029/2010JF001948, 2011. a

Stevens, L. A., Behn, M. D., McGuire, J. J., Das, S. B., Joughin, I., Herring, T., Shean, D. E., and King, M. A.: Greenland supraglacial lake drainages triggered by hydrologically induced basal slip, Nature, 522, 73–76, https://doi.org/10.1038/nature14480, 2015. a, b, c

Stevens, L. A., Behn, M. D., Das, S. B., Joughin, I., Noel, B. P. Y., van den Broeke, M. R., and Herring, T.: Greenland Ice Sheet flow response to runoff variability, Geophys. Res. Lett., 43, 11295–11303, https://doi.org/10.1002/2016GL070414, 2016. a, b, c

Sundal, A. V., Shepherd, A., Nienow, P., Hanna, E., Palmer, S., and Huybrechts, P.: Melt-induced speed-up of Greenland ice sheet offset by efficient subglacial drainage, Nature, 469, 521–524, https://doi.org/10.1038/nature09740, 2011. a, b

Tarboton, G.: A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resour. Res., 33, 309–319, 1997. a

Tedesco, M., Willis, I. C., Hoffman, M. J., Banwell, A. F., Alexander, P., and Arnold, N. S.: Ice dynamic response to two modes of surface lake drainage on the Greenland ice sheet, Environ. Res. Lett., 8, 034007, https://doi.org/10.1088/1748-9326/8/3/034007, 2013. a

Tedstone, A. and Neinow, P.: Ice motion measurements, southwest Greenland Ice Sheet, NERC, Polar Data Centre, http://doi.org/b2f3, 2017. a, b, c

Tedstone, A. J., Nienow, P. W., Gourmelen, N., Dehecq, A., Goldberg, D., and Hanna, E.: Decadal slowdown of a land-terminating sector of the Greenland Ice Sheet despite warming, Nature, 526, 692–695, https://doi.org/10.1038/nature15722, 2015. a, b, c, d, e

Van Angelen, J. H., Lenaerts, J. T. M., Van Den Broeke, M. R., Fettweis, X., and Van Meijgaard, E.: Rapid loss of firn pore space accelerates 21st century Greenland mass loss, Geophys. Res. Lett., 40, 2109–2113, https://doi.org/10.1002/grl.50490, 2013.  a

van de Wal, R. S. W., Smeets, C. J. P. P., Boot, W., Stoffelen, M., van Kampen, R., Doyle, S. H., Wilhelms, F., van den Broeke, M. R., Reijmer, C. H., Oerlemans, J., and Hubbard, A.: Self-regulation of ice flow varies across the ablation area in south-west Greenland, The Cryosphere, 9, 603–611, https://doi.org/10.5194/tc-9-603-2015, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

van den Broeke, M., Bamber, J., Ettema, J., Rignot, E., Schrama, E., van de Berg, W. J., van Meijgaard, E., Velicogna, I., and Wouters, B.: Partitioning recent Greenland mass loss, Science, 326, 984–986, https://doi.org/10.1126/science.1178176, 2009. a

Weinstein, M. and Rao, A. V.: ADiGator: A MATLAB Automatic Differentiation Tool, available at: https://sourceforge.net/projects/adigator/ (last access: 1 February 2015), 2011–2016. a

Werder, M. A., Hewitt, I. J., Schoof, C. G., and Flowers, G. E.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res.-Earth, 118, 2140–2158, https://doi.org/10.1002/jgrf.20146, 2013. a, b, c, d

Yang, K. and Smith, L. C.: Internally drained catchments dominate supraglacial hydrology of the southwest Greenland Ice Sheet, J. Geophys. Res. Earth, 121, 1891–1910, https://doi.org/10.1002/2016JF003927, 2016. a

Yang, K., Smith, L. C., Chu, V. W., Gleason, C. J., and Li, M.: A caution on the use of surface digital elevation models to simulate supraglacial hydrology of the Greenland ice sheet, IEEE J. Sel. Top. Appl., 8, 5212–5224, https://doi.org/10.1109/JSTARS.2015.2483483, 2015. a, b

Zwally, H. J., Abdalati, W., Herring, T., Larson, K., Saba, J., and Steffen, K.: Surface melt-induced acceleration of Greenland ice-sheet flow, Science, 297, 218–222, https://doi.org/10.1126/science.1072708, 2002. a, b