Modelling seasonal meltwater forcing of the velocity of the Greenland Ice Sheet

Surface runoff at the margin of the Greenland Ice Sheet 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 5 a multi-component model incorporating supraglacial, subglacial, and ice dynamic components applied to a land-terminating catchment in western Greenland produces modeled velocities which are in good agreement with those observed in GPS records for three melt seasons of varying melt intensities. This provides 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 development of efficient conduit drainage up-glacier from the ice sheet margin which develops 10 more extensively, and further inland, as melt intensity increases. This suggests current trends of decadal timescale slow-down in the ablation zone will continue in the near future, although the strong summer velocity scaling in our results could begin to offset potential future fall and winter velocity decreases for very high melt rates which are predicted for the end of the 21st century.


Introduction
Seasonal acceleration at the GrIS margin is driven by surface meltwater draining into the subglacial system.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 (Bartholomew et al., 2011;Fitzpatrick et al., 2013;Sundal et al., 2011).As the melt season progresses, a channelized system that efficiently drains water develops.This reduces water pressures and leads to a late summer deceleration (Bartholomew et al., 2010;Chandler et al., 2013;Cowton et al., 2013;Schoof, 2010).Understanding the impact of increased surface melting (Hanna et al., 2013;van den Broeke et al., 2009) 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, than a positive feedback between surface ing a subglacial hydrology model (Hewitt, 2013) with the ice flow model of Koziol and Arnold (2017).This model is driven by surface input from the surface hydrology-lake filling model ("SRLF" model) from Koziol et al. (2017), and initiated using the inversions from Koziol and Arnold (2017).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 (Morlighem et al., 2015), meteorological data constraining climatic input (Noël et al., 2015), and GPS data (Tedstone and Neinow, 2017) which provide a calibration and validation data set for model output.

Methods
The methods section begins with a description of the Russell Glacier study area and the data sets 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.

Study Area and Datasets
The Russell Glacier area is a land-terminating sector of the GrIS in Southwest Greenland.The study area 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 Yang and Smith (2016).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.
Two different topography datasets are used.The SRLF model is run with the 90 m resolution surface topography from the GIMP dataset (Howat et al., 2015).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 (Morlighem et al., 2014(Morlighem et al., , 2015) ) 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.This data is reinterpolated to 1000 m resolution.
Surface runoff and snow depth data for the SRLF model are provided by RACMO2.3 (Noël et al., 2015).Both runoff and snow depth are bilinearily interpolated from 11 km to 90 m.Three seasons with contrasting melt volumes are modelled: 2009, 2011, and 2012 (Fig .2) (Fig .2).Total melt over the SRLF study domain was 1.2 • 10 10 m 3 in 2009, 1.7 • 10 10 m 3 in 2011, and 2.1•10 10 m 3 in 2012.These three years serve as analogues for summers with with average, elevated, and extreme melt intensity respectively (following Koziol et al. (2017)).(Tedstone and Neinow, 2017).Purple diamonds show the locations of automatic weather stations (van de Wal et al., 2015).Cyan circles show the loactions of moulins used as tracer injections sites in Chandler et al. (2013).Inset shows the location in reference to Greenland.Ice Sheet Velocity Map dataset (Joughin et al., 2010a, b).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 Clason et al. (2015).
A crevassing threshold is selected by comparing the von Mises stress to observed patterns of crevassing in a Landsat 8 image.
A threshold value of 145 kPa is selected as optimal.
Moulin locations are specified as input data in the SRLF model.Moulin locations in the Russell Glacier area reported by Yang et al. (2015) 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 (Yang et al., 2015).As in Koziol et al. (2017), 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 (Tedstone and Neinow, 2017).A 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 in Fig. 1.

Supraglacial Hydrology (SRLF)
We use the supraglacial hydrology model (SRLF) from Koziol et al. (2017), run at 90 m resolution.A no-inflow boundary condition is imposed on all boundaries.Water is routed using a DEM of the surface of the ice sheet and a single flow direction algorithm (Arnold et al., 1998;?).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 (? Stevens et al., 2015), 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 Raymond and Nolan (2000).Overspill and channelized drainage can occur simultaneously if water enters a lake faster than can be evacuated by an existing channel alone.

Subglacial Hydrology
We use the subglacial hydrology model presented in Hewitt (2013) and Banwell et al. (2016).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 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.imposed, while the remaining boundaries have a no-flux condition.A concise model description is given here following Hewitt (2013) and Banwell et al. (2016) to provide context for the parameters used.However, for a detailed description the reader is referred to Hewitt (2013) and Banwell et al. (2016).
Discharge in the continuum sheet is modelled as: where h(x, y) is the thickness of the continuum sheet, K s is the sheet hydraulic conductivity, g is the acceleration due to gravity, and φ is the hydraulic potential.The hydraulic potential is defined as φ(x, y) = ρ w gb(x, y) + p w (x, y).
The distributed sheet thickness (h) is the sum of the thickness of the cavity sheet (h c ) and the elastic sheet (h el ).The cavity sheet evolves according to: where ρ i is the density of ice, m is the basal melting rate, U b is the basal sliding speed, h r is the bed roughness height scale, l r is the bed roughness length scale, A b 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 = ρ i gH − p w , where H is the ice thickness, p w is water pressure.
Basal melt rate is given by: The elastic sheet thickness is given by: where N − = min(N, 0), N + = max(N, 0), C el is an elastic compliance, and N 0 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 (Banwell et al., 2016).
Discharge in channels is modelled as: where K c is a turbulent flow coefficient, S is channel cross-section, and r is along channel distance.
The channel cross section evolves according to: where M is the melting rate along the channel wall.The melting rate along the channel walls is given by: where λ c is an incipient channel width.
The equation for mass conservation is: where Σ is englacial storage and R is the supraglacial input to moulins, The delta functions apply along channels (δ(x c )) and the positions of moulins (δ(x m )).
Englacial storage is represented as where σ is englacial void fraction and A m is moulin cross sectional area.
Model parameters held constant are shown in Table 1.Two parameters are assigned a spatially heterogeneous distribution in the calibration.The englacial storage parameter is assigned a background value of 10 −3 , with 50% of the cells then randomly set to 10 −4 .The effective sheet conductivity field is constructed using a background value of 10 −2 , with 15% of the cell nodes randomly assigned a value of 10 −7 .Since sheet conductivity is defined on the grid, neighboring nodes are averaged in the x and y directions to determine values on edges.At a sheet depth of 0.1 m, a sheet hydraulic conductivity of 10 −2 results in an effective hydraulic conductivity kh 2 of 10 −4 (Hewitt, 2013).This is at the upper end of values for till, which are inferred to be 10 −4 to 10 −9 (Fountain and Walder, 1998).The secondary value of 10 −7 for sheet conductivity assigned at nodes was selected to give an effective hydraulic conductivity at the opposite end of the spectrum.

Ice Flow/Inversion
The ice flow model implements the hybrid formulation of the ice sheet stress balance (Arthern et al., 2015;Goldberg, 2011), 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 (Hewitt, 2013), 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).
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, λ b is a bed roughness length (Hewitt, 2013).Following Hewitt (2013), negative effective pressures are eliminated by setting N + = max(N, 0), and regularized with a small regularization constant.
The linear sliding law (Eq.10) is used for the initial inversion of winter mean velocities, while the Weertman (Eq.11) (Budd et al., 1979;Hewitt, 2013) and Schoof (Eq.12) (Gagliardini et al., 2007;Schoof, 2005) 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 (Koziol and Arnold, 2017).It is based on automatic differentiation methods (Goldberg and Heimbach, 2013;Heimbach and Bugnion, 2009;Martin and Monnier, 2014), and uses the open source Matlab package AdiGator (Weinstein andRao, 2011-2016).The gradient of the cost function in this method is equivalent to one calculated using Lagrangian multiplier methods (MacAyeal, 1993;Morlighem et al., 2013) to generate the adjoint model (Heimbach and Bugnion, 2009).The cost function minimizes the weighted square of the difference of squares of measured and predicted velocities (Eq.13).A Tikanov regularization term is added for stability.
where γ 1 and γ 1 are scaling factors, Γ s is the surface domain, Γ b is the basal domain, w(x, y) is a weighting function, U obs (x, y) are observed surface ice speeds, U s 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 generalized Weertman sliding law, and µ b in the Schoof sliding law.The reported errors of surface velocities are used to calculate weights.

Model Integration
Englacial drainage is not considered in the SRLF model.For model integration, we assume that water drainage through the englacial system is strictly vertical, and that there is no horizontal transport in the englacial system.The SRLF model routes water into three different surface drainage pathways: moulins, lake hydrofracture, and crevasses.Moulins and surface to bed connections from lake hydrofracture are treated identically.All water entering these cells drains directly to the bed.However, drainage through crevasse field requires additional consideration.When water enters a cell determined as crevassed in the SRLF model, the water is removed from the model, and 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.Water drainage through crevasse fields is poorly understood, and the scheme implemented here (Figure 3) is motivated by simplicity.We assume all water in crevasse fields drain to bed, neglecting any refreezing.We also assume that contiguous areas of crevassed cells are hydrologically connected, perhaps by an internal water table.A crevassed cell in the SRLF model can accumulate water from two sources: 1) local ablation predicted by RACMO2; 2) a cell which is on the margin of a crevassed area may have water flowing into it from adjacent non-crevassed cells.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 • 10 5 m 3 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 veronoi 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 in it's 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 Hewitt (2013); the subglacial hydrology uses an implicit timestep using the current ice velocity distribution.After the state of the subglacial hydrology model in the next The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.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.

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.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; this corresponds to a run over a winter season (Sept 1 -April 30).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.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.As stated in Koziol and Arnold (2017), 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 (Colgan et al., 2012;van de Wal et al., 2015).
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 an inversion for each set of parameters tested.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.

Meltwater input partitioning
The majority of supraglacial meltwater drains into the englacial system (Fig. 5), consistent with observations (Zwally et al., 2002;Smith et al., 2015) and previous modelling (Koziol et al., 2017).Hydrofracture 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 10 concentrated in a small number of cells, with 50% of water routed over the ice sheet into crevasses entering in only 100 of the 7573 cells forming the perimeter of crevasse fields.Crevasse drainage is concentrated near the ice margin (Figure 6), while drainage into other pathways occurs throughout the study area.'Remaining Flow' refers to water still flowing over the ice sheet at the end of the model run.Water flowing into crevasses and moulins are 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.

Crevasse
Lake Hydrofracture    In summary, the early summer speedup and subsequent mid summer slow-down 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 (< 5 myr −1 ).
Ablation rates calculated from automatic weather stations near to sites S2, S4, and S6 are comparable to predicted RACMO2 surface runoff.The two data sets 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 correlate to predicted surface runoff, while they do show correlation at sites S4-S6.GPS measurements do not in general show correlation with the ablation rate at S2, and appear only weakly correlated at S4 and S6.

Model Sensitivity
Calibrating the integrated model is an underdetermined problem.Multiple parameters in each cell across the grid are constrained using only seven times 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 Werder et al. (2013) and Hewitt (2013).In this section, the focus is sensitivity of the model to the setup, and the parameters selected to have a spatially heterogeneous distribution.
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 • 10 5 resulted in crevasse input partitioning into 182 moulins and internal catchments.Changing the threshold value to 10 5 and 10 6 , 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 was tested by running a simulation where the water in the lake when hydrofracture occurs was not input to the base.The impact on ice velocities was found to be negligible.
A variable sheet conductivity 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 are overall very similar to the calibration runs, while decreasing the sheet conductivity value to 10 −3 leads to model output with prolonged periods of velocities exceeding 400 myr −1 .Increasing the initial coverage of lower conductivity nodes assigned a value of 10 −7 from 15% to 30% had a minor impact.Assigning 50% of the initial 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 conductivity nodes into 4x4 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.

Validation
The integrated model is validated against GPS velocity measurements from 2012 (Fig. 9), using the parameter values from the calibration run.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 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 is available at sites S5 and S7.Similarly to previous years, model output underpredicts GPS velocities at site S6.

Increased melt scenarios
The high melt scenarios accelerate the rate of the early summer speedup and result in higher peak velocities (Fig. 10).After the early summer speedup at sites S1-S3, simulations 2011, 2011x2, and 2011x4 all predict similar velocities.Modelled velocities can be observed to decrease slightly with increased melt season intensity during periods of slower flow (e.g.days 210-230 at S4).At sites S4-S6, increased melt input results in higher variability of ice flow in the first half of the summer season.Modelled velocity is observed to increase with melt input at these sites during the first half of summer.The relative increase of velocities between 2011x2 and 2011x4 is greater than between 2011 and 2011x2.At sites S4-S6 the model predicts similar velocities in all three simulations for the latter half of summer, from days 210 to 238.During this period at site S4, model velocities decrease with simulation melt input.At site S6, modelled velocities increase with greater melt input.Site S5 shows mixed behaviour, with model velocities from simulation 2011x4 higher than the simulation 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.Starting at day 238, a late season velocity spike is observed, with the magnitude of the velocity increase dependent on melt input.At site S7, the velocities from future simulations are slightly faster than 2011 velocities.The timing of events is similar in all three simulations.As melt input increases, the magnitude of short-term velocity spikes also increases, with a larger increase between 2011x2 and 2011x4 than between 2011 and 2011x2.

Channel Network Morphology/Extent
The development of the channelized system (see supplementary videos) is similar to that observed in previous modelling studies (e.g Banwell et al., 2016;Hewitt, 2013;Werder et al., 2013) and as inferred from observations (Bartholomew et al., 2011;Chandler et al., 2013).Channelization of the hydrological system begins at the margin and develops progressively up icesheet.As channelization develops up-ice, the system evolves to an arborescent morphology.The up-ice extent of channelization increases with summer melt intensity (Figure 13).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-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 2011x4 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 13 shows the locations of moulins used as tracer injections points in Chandler et al. (2013).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, nor that L41 and L57 would drain near L1.However, the model does predict that IS39 is on a different branch of the channelized system.Based on tracer measurements, Chandler et al. (2013)

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 (Figure 14).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-time scale variations in the distributed system but 10 with a dampened magnitude.At the end of the melt season, discharge in the distributed system decreases at a higher rate 22 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.than in the channelized system, so that there is a brief period in which discharge in channels is higher than in the distributed system.Under the future melt scenario 2011x4, 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 modelled melt seasons, discharge in the distributed system increases at a faster rate.However, peak discharge in the channelized system is nearly the same magnitude as peak discharge in the distributed system, and discharge in the channelized system exceeds that of the distributed system earlier in the year.Overall, model velocities are observed to be better at mid-elevation than either at 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 Chandler et al. (2013).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.Thin ice deviates from the assumption of a high aspect ratio in the hybrid formulation, while steep gradients are likely to lead to stresses assumed negligible in the 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 prescribed at the ice sheet margin, 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.Similarly to the study by Bougamont et al. (2014), modelled velocities do not capture this behavior.These velocity events may be the result of internal dynamics of water stored over winter, such as flooding events, that the subglacial hydrology model cannot capture, or early season melt which RACMO2 does not predict.
Modelled velocities at sites S6-S7 may be affected by excess capacity in the cavity system due to over prediction of basal ice velocities from the inversion process.The inversion process results in a sliding ratio of approximately 0.8 at the high elevations (Koziol and Arnold, 2017).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.73during the winter, increasing episodically to 0.9 during the summer (Ryser et al., 2014).The largest discrepancy between ablation at a weather stations and RACMO2 modelled surface runoff occurs at site S6, likely due to RACMO2 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.
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.

Model Sensitivity
Model velocities calculated with the two different sliding laws are comparable during much of the melt season.The timing of events are not effected 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 (Hoffman et al., 2016).Simulation results shows the Weertman 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 (Hewitt, 2013).
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 Hewitt (2013).Most model parameters of the integrated model are similar to previous studies applying the subglacial hydrology model (Hewitt, 2013;Banwell et al., 2016).The most significant parameter value difference is the sheet conductivity.The primary value of 10 −2 is between the order of magnitudes of 10 0 and 10 −5 used in Hewitt (2013) and Banwell et al. (2016) respectively.The parameter values for the model reported in Banwell et al. (2016), 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 (Fountain and Walder, 1998).
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 sheet height scale of 0.5 m is selected as a reasonable value in this paper.However, patterns of sheet thickness 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 lake hydrofracture events are not observed to have a pronounced effect in modelled velocities.This suggests that lake hydrofracture events are not a key process in the long term or large scale development of the subglacial hydrological system.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 exist.Rather, the primary impact of lake hydrofacture is in opening surface-to-bed connections, which then drain a significant proportion of the overall surface melt.
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.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 in 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.

Model Complexity
It is encouraging that the results provide a good match to 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 an important feedback (Hoffman and Price, 2014).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 sub-optimal occur throughout the study area.Such areas are characterized by: low 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 thermalmechanical 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 behavior is limited by the parameter space searched.However, a key behavior not replicated is the late summer/fall slowdown, and 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 inter-annual variability in the background state of the hydrological system.A model component simulating weakly connected regions of the hydrological system as incorporated in Hoffman et al. (2016) 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 Hoffman et al. (2011)).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 flow locally (Stevens et al., 2015;Tedesco et al., 2013).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 implies using observational records of lake drainages derived from satellite imagery (as in Bougamont et al. ( input to the ice-bed interface.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 may have an important role as the spatial density of water inputs are known to influence the development of the subglacial hydrological system (Banwell et al., 2016).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 need 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.

Implications
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, incorporating 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 on year changes, the practical value of the initialization process is implicitly validated through the subsequent fit to measured velocities.The model results also support the hypothesis that the margin of the GrIS is controlled by subglacial hydrology in a manner similar to alpine glaciers.
The timing of velocity variations are 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 impact of channelization on water pressures and water routing.
Modelling predicts that average 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 (van de Wal et al., 2015), but not in GPS records at North Lake, Western Greenland by Stevens et al. (2016).This implies that more intense melt seasons will result in a higher ice flux towards the margin during the summer.Whether this would be offset 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 (Tedstone et al., 2015;van de Wal et al., 2015), and is also postulated to result in slowdown in the subsequent winter and spring (Sole et al., 2013;Tedstone et al., 2015;van de Wal et al., 2015), model results suggest increasing summer melt intensity should 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.
However, model results show average summer velocities scaling with melt season intensity at mid-elevations.The modelled change in summer velocity at GPS stations 4 and 5 between the 2009 and 2011x4 melt scenarios is around 60 myr −1 .Assuming a winter velocity of 100 myr −1 and an 8 month-long winter, this increase would offset a reduction in winter velocity to around 60 myr −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 (Ryser et al., 2014).
Interpreting the model velocity output from the future melt scenarios is difficult.As melt season intensity increases, the validity of the initialization and calibration parameters 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 these will be offset 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 offsetting early summer speedup is less likely at higher elevations.

Conclusions
In this paper, multiple models are coordinated 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, Mean winter velocities are used for inversions of winter basal boundary conditions 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 3 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 1 .
Figure 1.Landsat 8 satellite image, 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

Figure 2 .
Figure 2. Daily surface runoff over the SRLF Russell Glacier study area for three contrasting summer melt seasons.
3) where τ b = (τ bx (x, y), τ by (x, y)) is the basal drag, u b = (u b , v b ) = (u(x, y, b), v(x, y, b)) 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.

Figure 3 .
Figure 3. Schematic drawing showing the conceptual model of the crevasse drainage implemented.Moulins are assumed to occur where high-flux supraglacial streams intersect the crevasse field.The crevasse field is than partitioned using Voronoi partitioning into internal catchments.All water in an internal catchments drains to the ice-bed interface at its corresponding moulin.

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

Figure 5 .
Figure 5. Pie chart of surface 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'.

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

3. 2
CalibrationModelled velocities are qualitatively calibrated against GPS measurements of horizontal surface velocities from 2009 (Fig7) and 2011(Fig.8).The calibration focused on matching the duration and magnitude of speedup events.The timing of the events is controlled by surface input, which was not varied.The model is calibrated using the Weertman sliding law, and the same parameter values are used in the simulations with the Schoof sliding law.Plots show model velocities output at noon for the Weertman sliding law, and daily averages calculated from output at 6 hr intervals for the Schoof sliding law.Subdaily variability in model output is relatively subdued, except during periods of high velocities in simulations applying the Schoof sliding law (seeKoziol (2017)).All 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 10 of ice velocities, surface runoff from RACMO2 and nearby surface ablation rates determined at weather stations(van de Wal et al., 2015) are plotted alongside velocities.RACMO2 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 accounts for both refreezing of meltwater and precipitation events.Refreezing, however, should The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.only be a small component (van deWal et al., 2015).An error of 5% is estimated for the calculated daily ablation rates (van deWal et al., 2015).

Figure 7 .
Figure 7. Modelled 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 Weertman sliding law are plotted in black and red respectively.Daily ablation from weather stations are shown in shaded blue, while RACMO2 surface runoff is shown in shaded red.Locations of GPS and weather station sites are shown Figure 1.Weather station ablation rates are plotted at the nearest GPS site.

Figure 8 .
Figure 8. Modelled 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 Weertman sliding law are plotted in black and red respectively.Daily ablation from weather stations are shown in shaded blue, while RACMO2 surface runoff is shown in shaded red.Locations of GPS and weather station sites are shown Figure 1.Weather station ablation rates are plotted at the nearest GPS site.

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 for 2009 and less englacial storage 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 2015, reducing the fit to GPS measurements.The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.
Melt season averaged modelled velocities at the GPS sites are shown in Figure11.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 ice margin.As melt season intensity increases, velocities in the upper ablation zone and at the equilibrium line (located at 1500 m elevation (van deWal et al., 2015), slightly above S6) are predicted to increase The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 9 .
Figure 9. Modelled 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 Weertman sliding law are plotted in black and red respectively.Daily ablation from weather stations are shown in shaded blue, while RACMO2 surface runoff is shown in shaded red.Locations of GPS and weather station sites are shown Figure 1.Weather station ablation rates are plotted at the nearest GPS site.

Figure 10 .
Figure 10.Modelled ice velocities using a Weertman Sliding law plotted for the 2011 melt season (blue), the 2x melt scenario (magenta), and for the 4x melt scenario (black).

Figure 11 .
Figure 11.Averaged modelled melt season velocities at each of the GPS sites.
report that channelization extends to at least L41, but not as far as L57.The modelled channelized system during 2009, 2011, and 2012 is inline with that result.The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 12 .
Figure 12. a) Map of melt season average velocities for 2009.b) Map of melt season average velocities for 2011x4.c) Change (%) between the melt season average velocities of 2009 and 2011x4.

Figure 13 .
Figure 13.Channelized system at maximum extent: a) 2009.b) 2011.c) 2012.d) 2011x4.Moulin locations used as tracer injections sites in Chandler et al. (2013) are shown in purple.Black lines correspond to 200 m surface topography contours at the same elevations as in Figure 1.

Figure 14 .
Figure 14.Time series of discharge in the distributed ("sheet") and channelized ("channels") system for the 2009, 2011 and 2012 summers, and the 2011x4 future melt scenario.
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 modeling 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 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.

Table 1 .
5194/tc-2017-225 Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.Constants used in the subglacial hydrology model during integrated runs in the Russell Glacier area.

Table 2 .
Constants used in the ice sheet/inversion model applied to the Russell Glacier Area.The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-225Manuscript under review for journal The Cryosphere Discussion started: 26 October 2017 c Author(s) 2017.CC BY 4.0 License.