Moulin density controls drainage development beneath the Greenland ice sheet

Uncertainty remains about how the surface hydrology of the Greenland ice sheet influences its subglacial drainage system, affecting basal water pressures and ice velocities, particularly over intraseasonal and interseasonal timescales. Here we apply a high spatial (200 m) and temporal (1 h) resolution subglacial hydrological model to a marginal (extending ~25 km inland), land‐terminating, ~200 km2 domain in the Paakitsoq region, West Greenland. The model is based on that by Hewitt (2013) but adapted for use with both real topographic boundary conditions and calibrated modeled water inputs. The inputs consist of moulin hydrographs, calculated by a surface routing and lake‐filling/draining model, which is forced with distributed runoff from a surface energy‐balance model. Results suggest that the areal density of lake‐bottom moulins and their timing of opening during the melt season strongly affects subglacial drainage system development. A higher moulin density causes an earlier onset of subglacial channelization (i.e., water transport through channels rather than the distributed sheet), which becomes relatively widespread across the bed, whereas a lower moulin density results in a later onset of channelization that becomes less widespread across the bed. In turn, moulin density has a strong control on spatial and temporal variations in subglacial water pressures, which will influence basal sliding rates, and thus ice motion. The density of active surface‐to‐bed connections should be considered alongside surface melt intensity and extent in future predictions of the ice sheet's dynamics.


Introduction
The Greenland ice sheet (GrIS), which has lost mass at an accelerating rate over the past two decades, is predicted to be the largest cryospheric contributor to global sea level rise for the rest of this century [Graversen et al., 2011;Hanna et al., 2013]. To predict the effects of an increasingly warmer climate on the GrIS, it is crucial to better understand the physical processes that govern its surface and dynamic mass balance [Church et al., 2013]. Seasonal acceleration and deceleration of marginal areas is influenced by the dynamic response of the subglacial drainage system to spatial and temporal variability in surface meltwater delivery to the bed via crevasses and moulins [Thomsen and Olesen, 1989;Bartholomew et al., 2012;Cowton et al., 2013;Smith et al., 2015], which often open up as a result of surface lake drainage events [Das et al., 2008;Tedesco et al., 2013;Doyle et al., 2013;Stevens et al., 2015]. However, the precise manner in which surface meltwater is delivered to the subglacial drainage system remains poorly understood. Similarly, there is a lack of knowledge about how the subglacial drainage system responds to the spatially and temporally varying surface water inputs, and how it therefore evolves to affect basal water pressure and surface ice velocity [Vaughan et al., 2013].
Numerous fieldand model-based studies have been undertaken to address this lack of knowledge and understanding. Field-based studies have attempted to infer the behavior of the GrIS's subglacial drainage system from measurements of meteorology, borehole water pressures, dye tracing, and surface velocity and uplift [e.g., Zwally et al., 2002;Das et al., 2008;Chandler et al., 2013;Cowton et al., 2013;Doyle et al., 2013;Meierbachtol et al., 2013;Tedesco et al., 2013;Andrews et al., 2014;Ryser et al., 2014;Stevens et al., 2015]. Several catchment-scale modeling studies have been undertaken which have tended to fit into one of two approaches. The first approach has used real boundary conditions and calibrated, modeled water inputs to simulate subglacial drainage for specific catchments in particular years. However, a key limitation of this approach is the somewhat idealized representation of the subglacial drainage system, comprising either channels only (the locations of which needed to be prescribed) [e.g., Colgan et al., 2011;Banwell et al., 2013;Mayaud et al., 2014] or a weak-sediment layer/porous sheet only [e.g., Bougamont et al., 2014]. • Combining models of supraglacial and subglacial drainage is an effective way of assessing the impact of moulin density on subglacial drainage • A higher moulin density causes an earlier onset of subglacial channelization and results in a more widespread channel network across the bed • Moulin density has complex and contrasting effects on water pressure in the early and late melt season, with basal sliding implications Supporting Information: • Supporting Information S1 • Movie S1 Correspondence to: A. Banwell,afb39@cam.ac.uk a more complex representation of the subglacial drainage system, which can take the form of inefficient and efficient flow in a porous sheet [e.g., de Fleurian et al., 2014] or inefficient flow in a sheet and efficient flow in channels [e.g., Schoof, 2010;Bartholomaus et al., 2011;Hewitt, 2013;Werder et al., 2013;Hoffman and Price, 2014]. These configurations emerge within the model in response to variable meltwater inputs and pressure gradients, with interaction and switching between the inefficient and efficient components dependent upon water flux. To date, therefore, catchment-scale hydrological models for marginal areas of the GrIS have tended to sacrifice an element of reality having either (i) realistic topographic boundary conditions and meltwater inputs but with a simplified subglacial hydrology or (ii) a more realistic representation of the subglacial drainage system but with artificial topographic boundary conditions and water inputs.
Here, to benefit from the advantages that both modeling approaches have offered, we develop a high spatial (200 m) and temporal (1 h) resolution subglacial hydrological model and apply it to a drainage basin (extend-ing~25 km inland) within the Paakitsoq region, West Greenland. The model is based on that of Hewitt [2013] but is adapted for use with both the real topographic boundary conditions and the calibrated, modeled water inputs from Banwell et al. [2013] for the 2005 melt season. These water inputs are in the form of moulin hydrographs, which are calculated by a surface routing and lake-filling model [Banwell et al., 2012bArnold et al., 2014], combined with a surface lake drainage model [Clason et al., 2012;Banwell et al., 2013;Arnold et al., 2014]. The model is forced with distributed hourly runoff calculated by a surface mass balance (SMB) model, which incorporates the surface energy balance as well as subsurface conduction and melting and refreezing in the snowpack [Banwell et al., 2012a].
The primary aim of this study is to explore how the subglacial hydrological system of the GrIS evolves in space and time in response to varying modeled moulin inputs associated with fluctuating patterns of surface melt, refreezing, routing through snow and across ice, and lake drainage events. A key output is the spatial and temporal variation in subglacial water pressure, which is of interest in helping to explain patterns of surface velocity and uplift found in previous studies [e.g., Das et al., 2008;Tedesco et al., 2013;Doyle et al., 2015]. As a secondary aim, we compare the results of the current study to those from Banwell et al. [2013]. As the inputs and topographic boundary conditions for each model are identical, any difference between the results will be due solely to the different representations of the subglacial hydrology, particularly the inclusion of a distributed system in the current model. Although this is the secondary aim, for continuity with the previous work, we discuss this aspect first, before moving on to the more substantial primary aim of the paper.

Study Site and Available Data
The Paakitsoq region is defined as~2300 km 2 of predominantly land-terminating ice, located in western Greenland, north of Jakobshavn Isbrae [Banwell et al., 2013, Figure 1]. The region was initially chosen by Banwell et al. [2012aBanwell et al. [ , 2012bBanwell et al. [ , 2013Banwell et al. [ , 2014 and later by Arnold et al. [2014] and Mayaud et al. [2014] due to the availability of various data sets including (i) hourly meteorological data measured at three GC-Net stations, JAR 1, JAR 2, and Swiss Camp [Steffen and Box, 2001]; and (ii) coastal precipitation and temperature data from the Asiaq Greenland Survey Station 437 (190 m above sea level, 4 km west of the ice margin), all of which were used to drive the SMB model [Banwell et al., 2012a]; (iii) proglacial stream discharge data measured at the Asiaq station, which were used to validate the glacier hydrological model of Banwell et al. [2013]; and (iv) a 750 m resolution bed digital elevation model (DEM) [Plummer et al., 2008] for the subglacial routing model and a 30 m resolution surface DEM taken from the Advanced Spaceborne Thermal Emission and Reflection Radiometer global DEM for the surface melt and routing models. Both DEMs were resampled to 200 m using bilinear interpolation.
The Paakitsoq region includes the areas for which the SMB model and the surface routing and lake-filling model have both been calibrated [Banwell et al., 2012a[Banwell et al., , 2012b. Additionally, the lake drainage model has been calibrated for the entire region, through the comparison of modeled lake volumes and drainage dates to those observed from nine Landsat 7 ETM+ satellite images ]. In the current study, our model domain, as defined by Banwell et al. [2013, Figure 1], is a~200 km 2 subglacial catchment within the Paakitsoq region that feeds the proglacial Asiaq gauging station ( Figure 1). The proglacial discharge measured here was used by Banwell et al. [2013] to specifically calibrate the lake drainage model for this 200 km 2 catchment. The domain is entirely within the ablation area of the ice sheet, extends~25 km inland from the margin with a surface elevation range of~150 m to 1000 m, and includes ice thicknesses < 815 m.

The Subglacial Model
The subglacial routing model is based on that developed by Hewitt [2013]. Although Hewitt's [2013] model has two distinct components, a model for subglacial water flow and a model for ice flow, we use only the subglacial water flow component in this study. This is because we lack appropriate temporally and spatially resolved data to constrain the modeled surface velocities. We focus specifically on analyzing spatial and temporal variations in subglacial water pressure in response to changing surface-derived meltwater inputs through the 2005 melt season. A brief description of the model is given below; full details are provided in Hewitt [2013, section 2.3]. The model is almost identical to that of Werder et al. [2013], except for its numerical implementation.

Model Description
The model is implemented in a two-dimensional finite-difference framework. It is based on a continuum "sheet" connected to discrete conduits, which are arranged along the edges and diagonals of a rectangular mesh of nodes. For low discharge, almost all meltwater is accommodated in the distributed sheet, whereas for sufficiently high discharge, the model exhibits an instability that leads to the formation of a self-organized channel system, in which a small subset of the conduits carries most of the discharge. Water storage is accounted for in an englacial aquifer and in moulins, the latter also acting as point sources of water to the subglacial system (see section 3.2.2).
The distributed sheet is separated into a cavity sheet, with thickness h cav (x,y,t), and an elastic sheet, with thickness h el (x,y,t); both are described below. These two components are added together to obtain the overall sheet thickness h(x,y,t) = h cav + h el . This is one of two main differences to the model setup in Hewitt [2013], where only one or other of these two models was used. The cavity sheet evolves according to where t is time, ρ w is water density (1000 kg m À3 ), ρ i is ice density (910 kg m À3 ), h r and l r are bed roughness height (0.1 m) and length scales (10 m), respectively, U b is the basal sliding speed, N is the effective pressure (P i À P w , where P i is ice overburden pressure and P w is water pressure), and A and n are Glen's law parameters (6.8 × 10 À24 Pa À3 s À1 and 3, respectively). The second term in equation (1) represents the opening of cavities due to ice sliding over a rough bed, and the third term represents creep closure of the cavity roofs. As our model does not incorporate the feedback of hydrology causing variable basal sliding, we must choose a representative value for U b . Based on measured surface ice velocities in this area, which range betweeñ 50 and~200 m yr À1 [e.g., Colgan et al., 2012;Ryser et al., 2014], and acknowledging that some of the motion is due to internal ice deformation, we take a representative value of 100 m yr À1 for U b .
The elastic sheet is included to account simplistically for the instantaneous uplift of ice, often referred to as "hydraulic jacking" [e.g., Röthlisberger and Iken, 1981;Das et al., 2008;Tedesco et al., 2013], that we expect to occur at high water pressure. The functional form of the thickness h el in Hewitt [2013] was taken from Flowers and Clarke [2002]; here we prefer a different functional form, believing that the elastic uplift we are trying to capture is more likely controlled by the effective pressure. This is the second main difference to the model setup in Hewitt [2013]. Here h el depends directly on effective pressure according to the equation where N À = min(N,0), N + = max(N,0), C el is an elastic compliance, and N 0 is a small regularizing pressure to make this function smooth (10 3 Pa). This function is designed to be zero for positive N but to increase rapidly when N is negative. In other words, this component of the sheet is activated when water pressure reaches or exceeds ice overburden pressure. This method of accounting for jacking is only approximate, since it does not allow for elastic bending stresses in the ice. However, the method allows for lake drainage events to be accommodated by the subglacial drainage system, without the generation of unrealistically large pressures that would occur if the water was forced into cavities and conduits. The default value for C el , 1.02 × 10 5 m Pa À1 , allows 1 m of uplift for each 10 m of excess hydraulic head. Sensitivity tests show that the key effect of increasing C el is to reduce the amplitude of the short-term (<24 h) water-pressure spikes when lakes drain; the overall water-pressure trends and drainage system development over the melt season are relatively insensitive to the precise value of C el .
Discharge q(x,y,t) in the cavity sheet is given by where ϕ(x, y) = ρ w gz(x, y) + P w (x,y) is the hydraulic potential, K s h 3 represents an effective hydraulic transmissivity as a function of the sheet thickness (h), g is acceleration due to gravity (9.81 m s À2 ), and z is bed elevation (m).
Channel cross-sectional area (CSA), S(s,t), with s denoting distance along a channel, evolves according to where M is the melt rate of the channel walls (equation (6)). The final term represents creep closure of the channel walls due to the effective pressure. The discharge in the channels is given by where K c is a turbulent flow coefficient (0.1 m s À1 Pa À1/2 ).
Melting, M(s,t), is given by Journal of Geophysical Research: Earth Surface

10.1002/2015JF003801
where L is the latent heat of melting (freezing point dependence on pressure is neglected here), and ⋋ is the incipient channel width (i.e., the width scale over which basal-ice melting contributes to channel initiation; 10 m [from Hewitt, 2013]).
Mass conservation is expressed as where englacial storage ∑(x, y, t) is a function of water pressure In equations (7) and (8), σ is the englacial void fraction, A m is moulin (a vertical, cylindrical shaft) CSA (10 m 2 ), and the delta functions apply along the positions for the conduits x c (s) and moulins x m . The source term, E(x,y, t), is exclusively from moulin point sources (section 3.2.2). The englacial water storage volume in void space and moulins, Σ, has a strong influence on both the amplitude and timing of diurnal variations; the more storage, the more damped the amplitude of the pressure variations, and the greater the delay relative to the melt signal. As a final point, we note that although the physics of the channelized component in the current model is fundamentally the same as that of Banwell et al. [2013], the models have been parameterized in slightly different ways.

Parameter Values
Results described in Hewitt [2013] and initial sensitivity tests undertaken in the current study show that the changing structure and morphology of the subglacial drainage system (notably, whether it is predominantly distributed or channelized, and the timing of the transition from one morphology to the other) is particularly sensitive to two key parameters: the connected englacial void fraction, σ, and the sheet flux coefficient, K s . Therefore, we identify the most suitable values for σ and K s as follows. We perform multiple model runs with different parameter combinations. We choose values within the ranges 1 × 10 À4 ≤ σ ≤ 1 × 10 À2 and 1 × 10 À6 ≤ K s ≤ 1 × 10 À4 Pa À1 s À1 , varying by factors of 10. This includes the order of magnitude values used by Hewitt [2013] and gives us nine parameter combinations in total. For each model run, we compare the modeled proglacial discharge with that measured at the Asiaq station using two objective criteria. First, we compare the total discharge volumes over the time period 12 May to 31 August 2005. Second, we analyze the match between the daily mean discharge volumes through calculation of both the Nash-Sutcliffe model efficiency coefficient and the root-mean-square error (RMSE).
We note, however, that the measured and modeled proglacial discharge hydrographs are not directly comparable because a series of proglacial lakes introduce a temporally varying lag between the arrival of water at the ice margin and its arrival at the Asiaq station, which is not accounted for by the model . Thus, for parameterization purposes, we direct more attention to the results of the total discharge volume comparison.

Boundary Conditions and Forcings
The full model is run from  (Figure 1). In reality subglacial catchments are able to vary in size due to subglacial water-pressure perturbations [cf. Lindbäck et al., 2015;Chu et al., 2016], but the model domain needs to be defined and fixed for the current study.

10.1002/2015JF003801
Due to the way in which the catchment was defined , no water enters the domain across the supraglacial or subglacial boundaries. Atmospheric pressure is applied at the ice margin. Water is able to outflow from the subglacial system at any point along the ice margin; outflow points are not predefined, as they had to be in Banwell et al. [2013]. However, if there are multiple outflow points, all water is cumulated at each time step for calibration purposes.

Moulin Positions
Following Banwell et al. [2013], we assume that all depressions in the surface DEM potentially contain a moulin in their lowest cell, which can be activated if and when a lake in a depression drains through a simulated hydrofracture event (see section 3.3.2.). Although we realize that additional moulins, not associated with surface lakes, are likely to be present on the ice sheet surface [e.g., Catania et al., 2008;Colgan and Steffen, 2009;Smith et al., 2015;Yang et al., 2015], our assumption produces a maximum moulin density of 0.20 km À2 , which is comparable to estimates calculated by Colgan and Steffen [2009] (0-0.88 km À2 ) and Zwally et al. [2002] (0.20 km À2 ) for Paakitsoq.
The locations of the potential moulins (black dots), overlaid onto a subglacial upstream flow accumulation map, are shown in Figure 2. The moulins fall on or very close to the paths of highest subglacial flow accumulation [Shreve, 1972], along which conduit locations were prescribed by Banwell et al. [2013].

Input Hydrographs 3.3.1. Surface Routing and Lake-Filling Model
Input hydrographs for all depressions (and hence lakes) in the surface DEM are calculated by the surface routing and lake-filling (SRLF) model [Arnold et al., 1998Arnold, 2010;Banwell et al., 2012bBanwell et al., , 2013, which is driven with surface runoff, simulated by the SMB model [Banwell et al., 2012a]. Briefly, the SRLF model links a surface lake and catchment identification algorithm to a flow delay algorithm, to simulate water flow over both snow-covered (assuming Darcian flow) and bare ice (assuming open-channel flow). However, in order to calculate moulin input hydrographs (i.e., the meltwater discharge which exits the depression through a  Table 2. Journal of Geophysical Research: Earth Surface 10.1002/2015JF003801 moulin in its lowest cell once a lake drains), the SRLF model is linked to a surface lake drainage model [Clason et al., 2012;Banwell et al., 2013;Arnold et al., 2014].

Surface Lake Drainage Model
The surface lake drainage model uses a water-volume-based threshold to trigger lake drainage events. We assume that all depressions start to fill to become lakes, but that they can begin to drain directly into the subglacial drainage system through a simulated hydrofracture mechanism if the lake reaches or exceeds a critical volume threshold. This threshold assumes that drainage will occur through a fracture if V = F a · H, where V is lake volume, H is the local ice thickness beneath the lake, and F a is map-plane fracture area. This concept is based on the idea that once a fracture has been initiated, it is the availability of surface meltwater for filling the expanding fracture and offsetting freezing onto the walls which is crucial in controlling crevasse propagation [van der Veen, 2007;Krawczynski et al., 2009]. Although the F a is treated as a tunable model parameter, it is constant across the model domain for a given model run; thus, lakes over thicker ice have to reach a larger water volume in order to drain.
If a lake reaches the threshold volume for drainage, all water in the lake is assumed to drain rapidly through a fracture. If the critical volume threshold is not reached and the surface depression reaches capacity, the SRLF model allows lakes to overflow into the next downstream catchment. Slow lake drainage by overflow channel incision [cf. Tedesco et al., 2013] is not accounted for. Following Banwell et al. [2013], and based on observations of rapid lake drainage events on the GrIS [e.g., Das et al., 2008;Doyle et al., 2013;Tedesco et al., 2013], we add the total lake volume at the time of drainage to the subglacial drainage system over a 5 h period (a period that the model is relatively insensitive to). Subsequently, all surface meltwater that enters the topographic depression is able to reach the subglacial system directly through the open moulin, which is assumed to exist at the lowest part of the depression for the remainder of the melt season [e.g., Shepherd et al., 2009;Bartholomew et al., 2010;Catania and Neumann, 2010].
A large F a (and the resulting high critical volume for drainage) means that few lakes will drain during the melt season; most meltwater will be stored in lakes or will flow to the ice sheet margin supraglacially. Any lakes that drain will deliver large volumes of water to the bed and will do so relatively late in the melt season; thus, the duration over which subsequent surface runoff can enter the subglacial system via an open moulin will be relatively short. Conversely, a small F a (and the resulting low critical volume for drainage) means that most lakes will drain within a few weeks of the melt season commencing, the total volume of meltwater either stored in surface lakes or routed supraglacially will be low, the water volumes delivered to the bed by lake drainage will be small, and the duration over which subsequent surface runoff can enter the subglacial system will be long.
To identify the most appropriate F a to use together with the local ice thickness in the lake volume threshold relationship, Banwell et al. [2013] compared modeled and measured proglacial discharge for the same 200 km 2 model domain as that used here; the optimal value was 1000 m 2 . In the current study, we first use the value F a = 1000 m 2 to calibrate key subglacial model parameters through comparison of the modeled and measured proglacial discharge hydrographs, and then we compare the results of this model run to those from Banwell et al. [2013]. However, since the model used here is different, this value is no longer necessarily optimal, but there is in any case uncertainty associated with the use of a constant F a across the model domain (several recent studies have highlighted the complexity of lake drainage). For the main part of this study we therefore experiment with a range of fracture areas, which allows us to explore how varying the density and opening time of moulins affects the behavior of the subglacial drainage system.
We consider the following three separate scenarios: (i) a large F a of 2500 m 2 (hereafter called the low moulin density scenario), (ii) a small F a of 250 m 2 (hereafter called the high moulin density scenario), and (iii) a F a of zero (hereafter called the maximum moulin density scenario). The third scenario effectively prevents lakes from filling, effectively assuming that all surface depressions contain moulins, which are always open. Although allowing no lakes to fill is unrealistic, this scenario allows for the maximum possible volume of surface meltwater delivery to the bed, with no delay due to surface storage in lakes.

Calibration of Key Parameters
The difference between total cumulative volumes of modeled (ΣMo) and measured (ΣMe) proglacial discharge varies significantly (À15 to 4% of ΣMe) for different combinations of parameter values for the Journal of Geophysical Research: Earth Surface 10.1002/2015JF003801 connected englacial void fraction (σ) and the sheet flux coefficient (K s ) ( Table 1). The minimum difference between ΣMo and ΣMe discharge is produced with the values σ = 1 × 10 -4 and K s = 1 × 10 -5 Pa -1 s -1 , which also give a relatively low RMSE and a relatively high Nash-Sutcliffe coefficient. However, in terms of the lowest RMSE and highest Nash-Sutcliffe coefficient, the best match between the modeled and measured discharge hydrographs is produced with the values σ = 1 × 10 À3 and K s = 1 × 10 À4 Pa À1 s À1 , although for these values, the difference between ΣMo and ΣMe discharge is higher. As minimizing the difference between ΣMo and ΣMe discharge is considered to be the best parameterization method (section 3.1.2.), we use the parameter values σ = 1 × 10 À4 and K s = 1 × 10 À5 Pa À1 s À1 for the remainder of this study.

Comparison of the Two Models
Under the assumption that F a = 1000 m 2 , we analyze the modeled water pressure in the same eight moulins that were shown in Banwell et al. [2013, Figure 5] (Figure 3; see Figure 2 for moulin/lake locations and Table 2 for lake drainage dates). To enable water pressures to be compared to the results from Banwell et al. [2013] and also across the model domain, we express water pressure as a fraction of ice overburden pressure, i.e., P w /P i . The input hydrographs for the nine moulins that open up for this scenario are shown in Figure S1 in the supporting information.
In general, there are strong similarities between the time series of modeled pressures from the two models ( Figure 3). Both produce a general transition from relatively high average subglacial water pressure with low amplitude diurnal cycles to lower average pressure and higher-amplitude cycles, as the melt season progresses; this change is indicative of increasing hydraulic efficiency (i.e., channelization) within the subglacial system [e.g., Bartholomew et al., 2010;Colgan et al., 2012;Chandler et al., 2013;Cowton et al., 2013] (see also Figure S2, Movies S1 and S2 and section 5.2 for further discussion). Furthermore, lake drainage events in both studies cause short-term rapid increases in subglacial water pressure up to, or above, ice overburden pressure (i.e., P w /P i ≥ 1). However, there are five key differences between the model results: 1. Pressure spikes due to lake drainage events in the current study reach a maximum of P w /P i ≈ 1.3, the magnitude of which is controlled by the elastic compliance (C el ) of the elastic sheet ( Figure 3) (e.g., setting C el to half its current value generates spikes in P w /P i ≤ 1.5). In contrast, pressure spikes in the study by Banwell et al. [2013] were capped at P w /P i = 1.1, as water at this threshold pressure overflowed from moulins and onto the ice sheet surface; 2. In the current study, the early melt season water pressure across the domain is relatively high (P w /P i ≈ 0.7), and following pressure spikes due to lake drainage events, pressures tend to remain high (i.e., P w /P i ≈ 1) for up to 7 days. In contrast, Banwell et al. [2013] found much lower early melt season water pressures (P w /P i ≈ 0) and the return of pressures to P w /P i < 0.5 within 2 to 3 days of a rapid drainage event (e.g., compare plots for Moulins 582 and 444, where lakes drained on 21 May and 12 June, respectively; Figure 3a).  (RMSEs) and Nash-Sutcliffe model efficiency coefficients between the modeled and measured daily mean proglacial discharge time series are also given (note that mean seasonal measured proglacial daily discharge = 44.0 m 3 s À1 ). The row in bold indicates the chosen parameter combination.
Journal of Geophysical Research: Earth Surface 10.1002/2015JF003801 3. Lake drainage events in the current study substantially increase subglacial water pressures within a very localized area (<1 km). Conversely, Banwell et al. [2013] found that lake drainage events caused concurrent high-pressure spikes (P w /P i ≥ 1) up to~5 km from the location of the drained lake. For example, in that study, the drainage of Lake 551 on 1 June caused a pressure spike in Moulin 572 (despite Lake 572 ). (a) Highlights four moulins that experience short-term (<24 h) spikes in water pressure, while (b) highlights four moulins that experience periods of sustained high water pressure (see Figure 2 for lake/moulin locations). The timings of lake drainage events are shown along the x axes of each plot. The black horizontal dashed lines show where P w /P i = 1.  Figure 2.
Journal of Geophysical Research: Earth Surface 10.1002/2015JF003801 not draining until 10 June) (Figure 3a), and the drainage of Lake 444 caused a pressure spike in Moulin 468. This phenomenon is not seen in the current study (Figure 3b). 4. The timing of the drop to lower average pressures, and of the emergence of higher-amplitude diurnal cycles (i.e., indicative of a switch from an inefficient system to more efficient channels), generally occurs later in the current study than in that by Banwell et al. [2013], particularly for the moulins farthest inland. For example, Moulins 468 and 494 (the highest-elevation moulins to open) show diurnally varying pressure from early-mid August and mid-late July, respectively, whereas Banwell et al. [2013] observed diurnal variations in these moulins from mid-June and early July, respectively (Figure 3b). Meanwhile, Moulin 572 (at a lower elevation,~15 km from the ice margin) shows diurnally varying water pressure from early July, compared to from mid-June in Banwell et al. [2013]. However, diurnally varying water pressure in Moulin 582 (also~15 km from the ice margin) is observed from mid-June in both models. 5. Diurnal variations in subglacial water pressure (as seen in many moulins from 1 July) in the current study are generally lower in amplitude than those found by Banwell et al. [2013] (Figure 3b).
Excluding point (1), the likely reason for the differences between the two sets of modeled water pressures summarized in points (2)-(5) is that the current subglacial routing model includes both a distributed and channelized system, whereas that used by Banwell et al. [2013] was composed solely of channels. Regarding point (2) in particular, although the channels in Banwell et al. [2013] could enlarge through wall melting during the melt season, thus increasing their efficiency, their relatively low capacity in the early melt season (minimum CSA = 0.07 m 2 ) was greater than the capacity of the early-season distributed system in the current study. The current study also has a background basal melt rate, which did not exist in the Banwell et al. [2013] study.
Regarding point (3), in the current study, lake drainage events result in rapid rises in water pressure within a localized region (<1 km), followed by relatively sustained high pressure (P w /P i ≈ 1.0) due to the accommodation of water in the surrounding distributed system. In contrast, lake drainage events in the Banwell et al.
[2013] study caused immediate pressure spikes in the subglacial channels and moulins connected to them up to~5 km from the draining lake, but pressures quickly fell back to near-atmospheric levels as the channels rapidly routed the water away.
Regarding point (4), the later occurrence of higher-amplitude diurnal cycles in water pressure in the current study, which is indicative of a later onset of channelization, is also due to the presence of a distributed system and the lack of a predefined channel network. In the current study, water initially flows in many directions and its dissipative melting of the ice is not concentrated along any specific pathway, thus channels develop relatively slowly and high water pressure is maintained for a relatively long time. In contrast, water in the Banwell et al. [2013] model was forced to follow channels that had prescribed locations and a minimum CSA; dissipative heating from water flow was focused on widening those routes, allowing conduits to enlarge more quickly than in the current study, and thus enabling water pressures to drop more quickly. Although the first scenario seems a priori more realistic, the latter scenario may better represent reality if preferred subglacial pathways are maintained through the winter [e.g., Gulley et al., 2012], which might be possible under the relatively thin ice of the GrIS margins.
Regarding point (5), the lower amplitude diurnal water-pressure variations in the current study are because the pressure is moderated by the ability of water to exchange between the distributed and channelized systems. At times of high surface meltwater inflow (e.g., in the afternoon), the channels pressurize, forcing water into the distributed system, whereas when the surface water inflow reduces (e.g., at night), water will flow out of the distributed system and into the channels [e.g., Hubbard et al., 1995]. Conversely, in the channel-only model of Banwell et al. [2013], the channels were likely at high pressure during the day and at atmospheric pressure during the night.

Model Sensitivity to Moulin Density
For the low moulin density scenario (i.e., F a = 2500 m 2 ), five lakes drain rapidly during the melt season (and thus, five moulins open, giving a moulin density of 0.02 km À2 ), the mean lake volume at the time of drainage is 1.0 × 10 6 m 3 , and the mean drainage date is 12 June (Table 2). For the "optimum" scenario (i.e., F a = 1000 m 2 , from Banwell et al. The proportion of fast-draining lakes, as a percentage of the total number of surface lakes in the model domain (40), is 13% and 45% for the low and high moulin density scenarios, respectively. Thus, these two moulin density scenarios likely span the realistic range of lake-bottom moulin input distributions given that the lower of these values (13%) is comparable to the percentages of observed fast-draining lakes across the entire GrIS (13% [Sundal et al., 2011]) and for southwest Greenland (14% [Selmes et al., 2013]).
In the following three sections, we discuss further the results for the low, high, and maximum moulin density scenarios. For comparison, the results for the optimum moulin density scenario are also shown in Figure S2, Movies S1 and S2; they are intermediate between the low and high moulin density scenarios.

Low Moulin Density Scenario
Lake drainage events in the low moulin density scenario occur at higher elevations as the melt season progresses and produce localized zones of high discharge and water pressure, which subsequently move down-glacier (Figure 4a; Movies S3 and S4). For example, following the drainages of Lake 444 on 18 June The solid black dots indicate moulins that are closed (i.e., "potential moulins"), and the open black circles indicate moulins that are open (i.e., those that are receiving water), with their size proportional to the surface water discharge entering the moulin. In the left column, the intensity of blue shading represents the subglacial discharge, and the red contour lines indicate hydraulic equipotential. In the right column, the intensity of red shows P w /P i , with P w /P i = 0 indicative of water at atmospheric pressure, and P w /P i = 1 indicative of water at ice overburden pressure.

Journal of Geophysical
( Figure 4a, 28 June) and then Lake 468 on 4 July (Figure 4a, 14 July), high-pressure waves steadily move down-glacier over the remaining melt season (Figure 4a, 30 July to 31 August). Such pronounced pressure waves develop when a sudden pulse of meltwater is received by an area of the subglacial system that has previously experienced few (or no) large meltwater pulses associated with lake drainage events and little (or no) subsequent surface meltwater inputs via an open moulin.
Thus, for the low moulin density scenario, the lack of an already developed channelized system means that water from a drained lake cannot easily be accommodated, or routed efficiently, by the drainage system, which is still largely distributed. It is relatively hard for water from a lake drainage event to initiate channel formation, despite the relatively large water volumes involved [e.g., Dow et al., 2014]. Additionally, as the lakes drain relatively late in the season, there is little time remaining for surface meltwater to enter the subglacial system to help develop and sustain a channel once it has been initiated. Therefore, by the end of the season, the bulk of the meltwater is still flowing in the cavity sheet, while some is carried by a few poorly developed channels covering a small proportion of the domain (Figure 4a, 31 August).
For the melt season as a whole, the mean water pressure for the low moulin density scenario is highest up glacier and in the more southerly regions of the domain, and lowest in the ice marginal and more northerly regions ( Figure 5a). Additionally, there are three central areas that have a lower mean water pressure (P w /P i~0 .45) than their neighboring areas (P w / P i~0 .75), suggesting that channels have grown to a large size as water fluxes are concentrated there.
It is more informative, however, to look at patterns of water pressure and its variability through the melt season, as it is these that have important implications for basal sliding [e.g., Pimentel and  Flowers, 2010;Schoof, 2010;Colgan et al., 2012]. For the low density moulin scenario, large areas of the bed experience high, sustained pressures for many days (Figure 6a), due to the dominance of a distributed system during the melt season. For example, 18% of the subglacial domain experiences high water pressures (i.e., P w / P i ≥ 1) for five days or more, while 8% of the domain experiences P w /P i ≥ 1 for 10 days or more. The lack of channels means that diurnal variations in water pressure are generally low in amplitude over the majority of the bed (Figure 6b); only 6% of the domain experiences a diurnal range in P w /P i ≥ 0.2 for 30 days or more, while only 0.4% of the domain experiences a diurnal P w /P i range ≥ 0.2 for 50 days or more. The area where pronounced diurnal water-pressure fluctuations occur for the longest time is in the north of the domain, where channels first develop.

High Moulin Density Scenario
The high moulin density scenario means that lakes drain relatively early in the melt season such that when a particular lake drains, there is a high likelihood that at least one nearby lake (probably down-glacier) has already drained. As previous lake drainage events, followed by subsequent moulin water inputs, will have increased the efficiency of the major conduits (at least over some areas of the bed), additional pulses of meltwater to the subglacial system will be accommodated and routed more easily. Thus, large areas of high water pressures, with pressure waves traveling down-glacier, do not dominate (Figure 4b; Movie S6), in contrast to the results for a low moulin density (Figure 4a; Movie S4). Instead, for a high moulin density, we see initially highly fluctuating water pressures (corresponding to diurnal fluctuations in surface meltwater inputs) in localized areas around the numerous moulins that are receiving meltwater after the lakes have drained (e.g., as indicated by the rapidly shifting hydraulic potential contours in Movie S5 and the highly variable P w /P i in Movie S6). Such localized areas of highly fluctuating water pressures advance up-glacier through the melt season as lakes at higher elevations drain, and their amplitude decreases as conduits become more efficient and meltwater is more efficiently routed away from the moulins (Figure 4b, 28 June to 31 August).
By the end of the melt season, a denser network of well-developed channels exists with only a small proportion of water accommodated by the cavity sheet (Figure 4b, 31 August; Movie S5), in contrast to the low density moulin results (Figure 4a, 31 August; Movie S3). We note that the positions of these main drainage pathways, as well as the main marginal outlet, are very similar to those that were assumed (and fixed for that model run) by Banwell et al. [2013, Figure 2].
Over the entire melt season, the mean water pressure for the high moulin density scenario (Figure 5b; overall mean P w /P i = 0.613 (standard deviation = 0.255)) is very similar to that for the low density moulin scenario (Figure 5a; overall mean P w /P i = 0.625 (standard deviation = 0.254)), which can be expected, given that lowpressure channels take up little space relative to the higher-pressure distributed system. However, there are some key differences in the season mean water pressures between the two scenarios ( Figure 5c). Whereas the high moulin density scenario results in lower mean pressures in the ice marginal/down-glacier region of the domain (i.e., the red/orange areas, Figure 5c), this scenario also results in higher mean pressures in some of the more up-glacier regions of the domain (i.e., the blue areas, Figure 5c). Thus, opening additional moulins leads to reduced mean water pressures in some areas, as there is sufficient water flux and/or time for channelization to occur, which lowers water pressure. However, in other areas, opening additional moulins results in higher mean water pressures, as there is insufficient water flux and/or time for channelization to occur, and thus the distributed system pressurizes further.
In terms of intraseasonal water-pressure patterns for the high moulin density scenario, only small areas of the bed experience multiday periods of high water pressures. For example, only 9% of the bed experiences P w /P i ≥ 1 for 5 days or more, and only 5% of the bed experiences P w /P i ≥ 1 for 10 days or more ( Figure 6c). This is due to the relatively rapid growth of efficient channels across a large portion of the bed from early in the melt season (Figure 4b, Movies S5 and S6), which also explains why high-magnitude and widespread, diurnal fluctuations in pressure are seen for a significant proportion of the melt season. For example, 20% of the bed experiences a diurnal range in P w /P i ≥ 0.2 for 30 days or more, while 10% of the bed experiences a diurnal range in P w /P i ≥ 0.2 for 50 days or more ( Figure 6d). Additionally, 1% of the domain experiences a diurnal range in P w /P i ≥ 0.2 for 70 days or more, and these areas appear as "ribbons" (i.e., the red areas, Figure 6d), where channelization is most pronounced.

Maximum Moulin Density Scenario
The development of the subglacial drainage system for the maximum moulin density follows a similar pattern as for the high moulin density scenario, although channelization occurs from even earlier in the melt season and results in an even denser end-of-season channel network ( Figure S3 and Movies S7 and S8). This is a consequence of all 40 moulins being open at the start of the season; lakes do not fill and drain, and all surface meltwater is routed to the bed as soon as it reaches the lowest cell in a topographic depression, with minimal surface water storage.
For the maximum moulin density scenario, the mean seasonal water pressure across the domain ( Figure S4a; overall mean P w /P i = 0.616 (standard deviation = 0.261)) is again similar to that for the low and high moulin density scenarios (Figures 5a and 5b, respectively), although the key differences between the low and maximum moulin density scenarios (as outlined above) are even greater ( Figure S4b). Regarding intraseasonal water-pressure variations, more areas of the bed experience sustained periods of high water pressures for the maximum moulin density scenario ( Figure S3), but these areas are smaller than for the high moulin density scenario (Figure 4b). Thus, the proportions of the domain that experiences P w /P i ≥ 1 for 5 and 10 days (or more) are actually identical to those for the high moulin density scenario (9% and 5%, respectively) ( Figure 6e). However, compared to the high moulin density scenario, the maximum moulin density scenario experiences high-magnitude, diurnal water-pressure fluctuations that are even more widespread and last for longer. For example, 34% of the bed experiences a diurnal range in P w /P i ≥ 0.2 for 30 days or more, and 15% of the bed experiences a diurnal range in P w /P i ≥ 0.2 for 50 days or more ( Figure 6f).
However, for the maximum moulin density scenario, the proportion of the bed that experiences a diurnal range in P w /P i ≥ 0.2 for 70 days or more (i.e., the red areas, Figure 6f) is < 1%, which is actually slightly less than that for the high moulin density scenario (Figure 6d). This is indicative of a slightly different pattern of conduit development; although numerous channels are initiated as a result of the more spatially distributed surface meltwater inputs, individual channels do not develop as much as they did for the high moulin density scenario because meltwater inputs are less focused into specific areas of the bed. In contrast, the high density moulin run receives more focused inputs (i.e., not all moulins open, and thus some surface depressions fill to form lakes, and then overflow into other lakes down-glacier, contributing to their volume and/or moulin inputs).

Comparison of the Two Models With Field-Based Studies
In general, strong similarities are seen between the time series of modeled subglacial water pressures from this study and those from Banwell et al. simple model captures key aspects of the system behavior observed in the current, more complex, model. As both studies have identical inputs and boundary conditions, the five main differences that we observe (section 4.2) are due solely to the different representations of the subglacial hydrology, most importantly, the presence of a distributed system in the current model compared to the channel-only model of Banwell et al. [2013].
In comparison to the Asiaq measured proglacial discharge data, the current model produces a closer statistical match than that found by Banwell et al. [2013] (section 4.1). However, as we have no other field-based hydrology data for the precise region in 2005 (e.g., moulin/borehole water levels or englacial/subglacial water velocities from tracing experiments), and as our models do not predict basal sliding rates (which could be compared to fieldor satellite-derived ice velocity data), it is hard to validate the models to establish which one gives a better representation of reality. However, in the following two subsections, we compare our results with field-based measurements of (i) moulin/borehole water levels; (ii) water velocities from tracing experiments; and (iii) ice velocities, all from other marginal regions of the GrIS, and/or from different years.

Comparison With Moulin/Borehole Pressure Measurements
With reference to measured borehole pressure data from marginal regions of the GrIS [e.g., Lüthi et al., 2002;Andrews et al., 2014;Wright et al., 2016], we suggest that the higher early-season water pressure found by the current study (i.e., P w /P i ≈ 0.7 before any lake drainage events, cf. P w /P i ≈ 0 in Banwell et al. [2013]) is the more realistic of the two set of model results. Additionally, diurnally fluctuating borehole water pressures from marginal regions (<30 km inland) of the GrIS have been measured at 0.76 ≥ P w /P i < 1.17 van de Wal et al., 2015;Wright et al., 2016], suggesting that the middle-to late-season lower amplitude cycles found in the current model are more realistic than those found by Banwell et al. [2013]. Finally, Wright et al.
[2016] first observed diurnally varying, and a drop to lower average, borehole water pressures (measured 27 km inland from Isunnguata Sermia, West Greenland) in late June 2011 and 2012, which is more in line with the results of the current model. Cowton et al. [2013] and Chandler et al. [2013] used dye-trace data to infer that a degree of channelization had occurred in the lower 14 km of Leverett Glacier, West Greenland, by late May/early June 2010. This is 2 weeks earlier than when the channel-only model of Banwell et al. [2013] saw signs of a significant increase in channel efficiency~15 km from the ice margin and closer to 4 weeks earlier than when the current model starts to show signs of channelization (e.g., cf. Moulins 551 and 572, Figure 3a).

Comparison With Ice Velocity Measurements
In contrast to the results of Cowton et al. [2013] and Chandler et al. [2013], Hoffman et al. [2011] observed a sustained (>40 days) increase in surface ice motion (>50% above winter background) from middle to late June in 2007 in the area around JAR1 in the Paakitsoq region, suggesting that widespread channelization did not occur until middle to late July. Likewise, Ryser et al. [2014] did not observe any significant increases in ice velocity at FOXX (~20 km inland from the ice margin at Paakitsoq) until early to mid-June in 2012, again suggesting that channelization did not occur in this region for at least a few weeks after this. Thus, ice velocity measurements from Paakitsoq, but from other years, appear to be more in line with the results of the current study than the Banwell et al. [2013] study.

Comparison Summary
The current model, which includes both distributed and channelized drainage components, appears to better match previous studies of borehole water pressures and ice velocities, in terms of the magnitude and timing of water-pressure fluctuations and the onset of channelization to middle to late June, than the channelonly model of Banwell et al. [2013]. In contrast, the timing of the onset of channelization inferred from tracer experiments appears to be better matched in the channel-only model than the current model, although it still occurs about 2 weeks earlier than the channel-only model predicts.
Possible reasons for the delay of channel initialization in our current model compared to that inferred by Cowton et al. [2013] and Chandler et al. [2013] are as follows: (i) an overprediction of conduit creep closure rates during low discharge, (ii) the localized uplift mechanism during overpressured periods not fully accounting for the nonlocal effects of hydraulic jacking away from the moulins, and (iii) not accounting for the role of sediment erosion in conduit enlargement. However, the earlier channelization inferred by Cowton et al. [2013] and Chandler et al. [2013] is also likely to be because 2010 was a particularly intense melt year [e.g., van As et al., 2012].
Journal of Geophysical Research: Earth Surface 10.1002/2015JF003801

Effect of Moulin Density on Subglacial Drainage Development
Our results show that there is a gradual up-glacier transition from an inefficient (distributed) to an efficient (channelized) system through the melt season for all moulin density scenarios (Figures 3, 4, S2, and S3; Movies S1, S3, S5 and S7). This progression follows an increase in the magnitude and variability of meltwater inputs to the bed, which is due to increased melt rates and decreased refreezing and storage in the snowpack [cf. Bartholomew et al., 2010;Colgan et al., 2012;Chandler et al., 2013;Cowton et al., 2013]. In effect, the growth of a channelized system follows the up-glacier development of surface drainage and advances in a stepwise manner as more moulins open at progressively higher elevations over the melt season, similar to the situation that has been reported from valley glaciers [Nienow et al., 1998;Willis et al., 2002]. This general up-glacier progression of lake drainage events and moulin activation occurs because at lower elevations (i) surface melting commences earlier; (ii) melt rates are higher; and (iii) ice is generally thinner, meaning that lakes do not need to fill to such a large volume before they drain as they would for thicker ice at higher elevations [McMillan et al., 2007;Sundal et al., 2011;Banwell et al., 2013;Arnold et al., 2014].
However, as shown by the contrasting results of model runs with low, high, and maximum moulin densities (Figures 4, S3, and 5), the timing of channel initialization and the final extent of channelization across the bed are also strongly controlled by the opening time and overall density of moulins that develop during the melt season. For the low moulin density scenario, few lakes drain, they do so relatively late in the season, the onset of channelization is relatively late, and just a few well-developed channels exist by the end of the season (Figure 4a). In contrast, for the high moulin density scenario, many more lakes drain, they do so earlier in the season, thus, the onset of channelization is earlier, and by the end of the season, a denser network of channels exists over a larger proportion of the bed (Figure 4b).
Superimposed on the general trend of a seasonal increase in subglacial drainage efficiency, for all moulin density scenarios, are three distinctive phenomena: (i) water-pressure spikes (i.e., where P w /P i ≥ 1 for < 24 h), (ii) sustained high-pressure periods (i.e., where P w /P i ≥ 1 for > 5 days), and (iii) high-amplitude diurnal variations in water pressure (i.e., where P w /P i ≤ 0.2 and ≥ 0.8). Pressure spikes are generally caused by the rapid inflow of large water volumes at the time of lake drainage events (Figure 3a), which have also been observed in subglacial pressure borehole data [e.g., Lüthi et al., 2002]. Unsurprisingly, our results suggest that more spikes in water pressure are produced when the density of moulins is higher, because these moulins open due to rapid lake drainage events.
Periods of sustained high pressure are generally observed immediately after, and in areas down-glacier from, lake drainage events (Figures 3b, 4, S2, S3, 5a, 5c, and 5e), consistent with other fieldand modeling-based studies [Pimentel and Flowers, 2010;Hoffman et al., 2011;Palmer et al., 2011;Banwell et al., 2013;Dow et al., 2015]. This is because high volumes of surface meltwater continue to enter the distributed drainage system after the initial drainage event, causing the system to exhibit a positive relationship between water discharge and pressure until channels have had time to develop [Spring and Hutter, 1981;Pimentel and Flowers, 2010;Schoof, 2010]. However, if the discharge and pressure gradient are not sufficient to initiate channelization, waves of high sustained pressure will move down-glacier, as observed following the drainage of Lakes 444 and 468, for instance (section 4.3.1). This indicates the importance of a turbulent sheet for enabling water to move readily at the ice-bed interface in the absence of extensive channelization [e.g., Dow et al., 2014]. Compared to the higher moulin density scenario, which causes a greater rate of channelization across the bed (Figure 4b and Movie S6), such pressure waves are more prevalent in the low moulin density scenario due to the relatively slow development of a channelized system (Figure 4a and Movie S4).
Large diurnal variations in P w /P i occur middle to late melt season (Figures 3, 4, S2, S3, 5b, 5d, and 5f). These have also been observed in field-derived borehole data from western Greenland [Thomsen et al., 1991;Cowton et al., 2013;Meierbachtol et al., 2013;Andrews et al., 2014;Ryser et al., 2014;Doyle et al., 2015;van de Wal et al., 2015;Wright et al., 2016] and are due to the pressurization of already developed channels that are unable to enlarge fast enough to accommodate the surface-derived meltwater pulses. This is in support of the finding that the variability rather than magnitude of meltwater inflow has the greatest impact on water pressures [e.g., Schoof, 2010;Colgan et al., 2011;Bartholomew et al., 2012]. In the present study, the large diurnal P w /P i variations coincide with times of high diurnal variability in moulin inputs and occur in areas where channels have become most developed. Thus, a higher density of moulins results in longer time periods where large diurnal variations in P w /P i prevail, which persist over a larger area of the bed.
Journal of Geophysical Research: Earth Surface 10.1002/2015JF003801 Banwell et al. [2013] fixed their subglacial conduit locations under the assumption that P w /P i = 0.925, which was based on the theory that major drainage pathways are created when P w /P i is close to/at ice overburden pressure, but then remain in similar positions once water pressures drop, since they are unable to migrate laterally to areas of the bed with a lower hydraulic potential [e.g., Sharp et al., 1993;Hubbard et al., 1995]. The current study confirms that this was a reasonable assumption as our results (particularly for the higher moulin density scenarios) suggest that the overall drainage network structure is not only relatively static during the melt season (Figure 4b), but that the end-of-season drainage network (Figure 4b, 31 August) is also very similar to that prescribed by Banwell et al. [2013, Figure 2].

Implications of Moulin Density for Sliding
Although our subglacial routing model is not coupled to a sliding model, it is possible to infer where and when enhanced basal sliding is likely to occur by analyzing our modeled spatial and temporal patterns of P w /P i in the context of results from other model-and field-based studies. Such studies have indicated that the three distinct water pressure phenomena discussed earlier (section 5.2) have implications for basal sliding and thus ice velocity. First, water pressure spikes due to lake drainage events may cause or prolong short-term ice uplift and speedup in the early to middle melt season [e.g., Das et al., 2008;Van de Wal et al., 2008Pimentel and Flowers, 2010;Doyle et al., 2013;Tedesco et al., 2013;Ryser et al., 2014;Stevens et al., 2015]. Our results suggest that a high moulin density would lead to greater, and more widespread, short-term velocity increases during the early to middle melt season, because more water pressure spikes occur due to more lake drainage events. This result, however, is directly related to how we prescribe our moulins to open.
Second, sustained high water pressure during the early to middle melt season has been inferred from measurements of longer-term velocity increases of marginal areas of the GrIS during the so-called "summer regime" [Shepherd et al., 2009;Bartholomew et al., 2010;Hoffman et al., 2011;Palmer et al., 2011;Colgan et al., 2012;Ryser et al., 2014]. Compared to a high moulin density, our results suggest that a lower moulin density results in more localized channel development from relatively late in the season, causing more sustained high water pressures; this would provide a greater potential for a more sustained and widespread ice velocity increase.
Third, large, short-term (typically diurnal) variations in water pressure in the later part of the melt season have been shown theoretically to generate short-term high-velocity events, even once the drainage system is predominantly channelized [Pimentel and Flowers, 2010;Schoof, 2010]. Such high-velocity events have also been observed in field-based studies on the GrIS in the middle to late season [Shepherd et al., 2009;Hoffman et al., 2011;Bartholomew et al., 2012;Cowton et al., 2013;Andrews et al., 2014;Ryser et al., 2014]. As our results suggest that a higher moulin density leads to more widespread channelization from earlier in the season, which leads to larger areas of the bed exhibiting high diurnal water-pressure variations, a higher moulin density may result in a higher prevalence of short-term late-season high-velocity events, compared to a lower moulin density.
However, the picture that is emerging in terms of the overall impact of drainage system evolution on ice velocity is very complex. Early to midseason channelization, which we suggest may be correlated with a higher moulin density, leads to a gradual reduction in average water pressures and ice velocities over that time [e.g., Sundal et al., 2011;Bartholomew et al., 2012;Cowton et al., 2013;Tedstone et al., 2015]. However, early channelization may also cause greater fluctuations in water pressures later in the season, and this may lead to higher average late-season ice velocities. The extent to which short-term increases in late-season velocity may counteract the reduced early to middle season velocities requires further study.
Additionally, over annual and decadal timescales, the impact of spatially variable and temporally variable surface melting and water delivery to the bed on velocities requires further study. Recent work suggests that for marginal areas of the GrIS, increased surface melting produces faster ice motion during summer but slower movement over the following winter, as a result of the growth of an efficient drainage system that allows more water to drain from larger areas of the ice sheet bed with a high basal water pressure Tedstone et al., 2013]. Furthermore, there appears to be a cumulative effect, with successive years of high surface melting causing a year-on-year increase in the extent of channelization across the bed, and an overall reduction in annual velocities, at least for land-terminating regions with ice thicknesses <~1100 m [Tedstone et al., 2015].
The results of this study show is that it is not only the volume and areal extent of surface melt that is important in driving subglacial drainage development, water pressures, and therefore ice velocities but also areal Journal of Geophysical Research: Earth Surface 10.1002/2015JF003801 moulin density, how many moulins open, and when and where they open. If moulin density varies from yearto-year, so too will the timing of channel initialization and the peak extent of channelization across the bed. Therefore, the effects of this on water pressures and ice motion will be combined with any changes in water pressures or ice motion due to increased melt rates.

Model Limitations
This study has employed a surface lake drainage model, which uses a water-volume-based threshold equal to a prescribed fracture area (a tunable parameter) multiplied by the local ice thickness, to simulate lake drainage events (and thus moulin opening). The way in which the subglacial drainage system evolves, and the resultant spatial and temporal patterns of water flux and pressures, is very sensitive to the value of this parameter, which controls the pattern and timing of water delivery to the bed. It is for this reason that the primary aim of this study involves the use of a variety of fracture areas, to create low, high, and maximum moulin density scenarios.
To better understand the basal conditions beneath the GrIS, more research is required into the precise fracture mechanism and the factors that influence it [e.g., Stevens et al., 2015]. Similarly, there is a need to establish the importance of fractures outside of lake basins in capturing water and delivering it to the bed. Such research will enable a more accurate prediction of water delivery to the ice sheet bed, including lake drainage events, over space and through time.
Finally, to calculate and to predict changes in ice dynamics of marginal regions of the GrIS, this or similar models should be coupled to a sliding model [e.g., Hewitt, 2013] and forced with realistic boundary conditions and water inputs. Ideally, such a model should be constrained with field-derived data, such as borehole water pressures and surface ice velocities, to better establish how moulin density affects ice dynamics over intraseasonal and interseasonal timescales. Once this has been established for the past, such a model could be used to examine future ice sheet hydrology and dynamics. The need for such research is becoming increasingly important as future melt rates are predicted to rise [e.g., Hanna et al., 2013], which could lead to earlier lake drainage events (and thus earlier moulin opening times) at lower elevations [e.g., Mayaud et al., 2014], and possibly higher-elevation lake drainage events [cf., Leeson et al., 2014;Poinar et al., 2015] (thus higher-elevation moulin openings).

Conclusions
We have developed and applied a high spatial (200 m) and temporal (1 h) resolution subglacial hydrological model to a~200 km 2 land-terminating, ablation area of Paakitsoq, West Greenland, for the 2005 melt season. The subglacial model, which is based on that of Hewitt [2013], has both distributed and channelized components and has been adapted for use with both the real topographic boundary conditions and the calibrated, modeled moulin inputs from Banwell et al. [2013].
Using the optimum moulin density established by the earlier study of Banwell et al. [2013], we find strong similarities between the time series of modeled subglacial water pressures from this study and those from Banwell et al. [2013], particularly in the middle to late melt season when large areas of the bed have become channelized. Compared with the earlier study, we find a more complex set of behaviors and an improved fit with the observed proglacial discharge data, differences that are attributed to the presence of a distributed system in the current study.
Our results show that the spatial and temporal variability of surface meltwater reaching the ice sheet bed provides a fundamental control on subglacial drainage system development. In the current study, this process was investigated by varying the areal density of moulins, which was controlled by a simple ice-fracture mechanism affecting the location and timing of lake drainage events.
The different moulin densities cause the subglacial drainage system to evolve in contrasting ways, creating different patterns of water pressure, which might be expected to generate different ice dynamic responses. For a high moulin density, the change from a predominantly distributed to channelized system occurs earlier in the melt season and over a larger area of the bed, than for the low moulin density scenario. This leads to a more rapid and widespread reduction in water pressure from earlier in the melt season, which might be expected to cause a decrease in early-season basal sliding and thus the total annual displacement of ice. The Journal of Geophysical Research: Earth Surface 10.1002/2015JF003801 corollary of this, however, is that larger areas of the bed experience high diurnal water-pressure variations for longer throughout the late season for the high moulin density scenario compared to the low moulin density scenario, which might be expected to cause short-term increases in basal sliding at this time. The implications of all this for the overall summer and annual velocity regimes, both now and into the future, require further study.
In summary, our results suggest that although the intensity and areal extent of surface melt have a strong control on subglacial drainage development, water pressures, and therefore ice velocities, moulin density is also of crucial importance as this controls where and when this meltwater can access the bed during the melt season.