Map of Russell Glacier, in West Greenland, with the model outline (red rectangle), the location of moulins (red stars), and supraglacial lakes (blue dots). The supraglacial watersheds (light blue) are delineated from the supraglacial rivers (dark blue) using a Landsat 8 image from August 2013 (Landsat8_2013_231_B8) which is used as an overlay. The green squares are the locations of in situ measurements of GPS velocity presented in] and used for comparison with the model results. The yellow triangle is the location of the SHR station where water pressure was recorded []. The black lines are the limits of the different regions used to describe the results: the frontal region (F) below 400 m, the lower region (L) between 400 and 800 m, the median region (M) between 800 and 1200 m, the upper region (U) between 1200 and 1600 m, and finally, the highest region (H) above 1600 m.

In this study, we focus on the modeling of subglacial water pressure without coupling the hydrological model with a numerical ice flow model. Our objective is to investigate the effect of interannual runoff volume variability on water pressure at the glacier base. By doing so, we avoid the difficulty of introducing a friction law in our modeling scheme that adds complexity to the analysis. We base our study on the Russell Glacier region in West Greenland, where surface mass balance models [ van de Wal et al. , 2012 ] are well constrained by in situ data. The computed mass balance then allows to reconstruct runoff needed for the forcing of the subglacial hydrology model. This region also provides reliable displacement data [ Tedstone et al. , 2013 ] and water pressure measurements [ Smeets et al. , 2012 ], which are used to assess the reliability of the hydrological model output (Figure 1 ). The subglacial hydrological model is based on the double‐continuum approach described in de Fleurian et al. [ 2014 ] with improvements on the treatment of the efficient component of the drainage system. The effective pressure computed by the model is compared to GPS measurements of vertical displacements to evaluate the model and its parameters. Finally, we discuss on the role of runoff in controlling basal water pressure and explain how the glacier responds dynamically depending on the state of the subglacial hydrological system.

Subglacial hydrological modeling may provide useful insights into the future evolution of the Greenland Ice Sheet because the models can be employed at the spatial and temporal resolutions needed to study the response of water pressure to changes in runoff. A number of subglacial hydrological models [e.g., Werder et al. , 2013 ; de Fleurian et al. , 2014 ; Hoffman and Price , 2014 ] now simulate water pressure based on our current understanding of subglacial drainage, which helps improve our understanding of the feedbacks between basal water pressure and runoff availability. These numerical models are generally designed to be coupled to higher‐order ice flow models in order to include the feedback between meltwater availability and glacier dynamics. Previous studies on both synthetic [e.g., Pimentel and Flowers , 2010 ; Schoof , 2010b ] and real cases [e.g., Bougamont et al. , 2014 ] have already shown the capability of these models to improve our knowledge of the interactions between surface runoff and glacier dynamics.

On Russell Glacier, observations show that large changes in meltwater availability (e.g., in year 2012 runoff was twice as large as in year 2009) seem to lead to a stagnation or even a slowdown in annual ice velocity at low elevation (below the Equilibrium Line Altitude, ELA) [ Sundal et al. , 2011 ; Sole et al. , 2013 ; Tedstone et al. , 2013 ; van de Wal et al. , 2015 ]. The seasonal variability shows an acceleration at the onset of the melt season, linked to the availability of surface meltwater, while the subglacial drainage system is still weakly developed [e.g., Bartholomew et al. , 2010 ; Palmer et al. , 2011 ], which is followed by a significant slowdown when the volume of runoff decreases, typically at the end of the summer [ Hewitt , 2013 ; Sundal et al. , 2011 ; Schoof , 2010a ]. However, sparse GPS velocities from higher‐elevation areas show a persistent interannual acceleration of the ice sheet [ Doyle et al. , 2014 ]. Recent studies based on remote sensing data allow to obtain the spatiotemporal resolution needed to determine the effect of subglacial hydrology on the interannual velocity of glaciers [ Tedstone et al. , 2015 ], but the spatial coverage is still mostly limited to the lower part of the glacier.

Contrary to internal deformation, which is well known and constant in time in the absence of major glacier thinning, basal slip varies seasonally and interannually [e.g., Hewitt and Fowler , 2008 ]. It is commonly agreed that water pressure at the glacier base is the main driver for fluctuations in basal velocity of land‐terminating glaciers [e.g., Iken et al. , 1993 ], but a simple relationship between runoff volume and basal water pressure does not exist and the modeling of subglacial hydrology remains an open question [ Flowers , 2176 ]. Moreover, the difficulty of accessing the glacier base makes it challenging to obtain observations required to evaluate theories of basal slip.

Ice surface velocity is the result of both internal deformation of the ice body and basal slip defined as the combination of glacier sliding on bedrock and, if present, the deformation of a sediment layer between the glacier and bedrock [ Cuffey and Paterson , 2010 ]. Basal slip has been measured on a number of mountain glaciers to determine that this component accounts for about half of the observed surface displacement [ Cuffey and Paterson , 2010 ], but it can also account as much as 95% of the overall surface velocity in some cases [e.g., surging glaciers Kamb et al. , 1985 ]. Observations on ice sheets are more scarce, but the available data suggest that basal slip plays a major role in the observed surface velocity of fast‐flow regions of the ice sheets [e.g., Morlighem et al. , 2013 ; Luthi et al. , 2002 ; Engelhardt and Kamb , 1998 ].

The material and physical parameters of the model are given in Table 2 . For the hydrological parameters, the selection is difficult. This is a result of a lack of data to constrain the physical parameters and also of our particular modeling approach where the use of equivalent layers is based on parameters that are not measurable. Hydrological parameters were selected by fitting the model results of the 2008 simulation to the velocity recorded at three locations throughout the 2008 melt season [ Bartholomew et al. , 2010 ]. Our main goal while setting the parameters is to get the correct timing for the evolution of water pressure with respect to the velocity variations observed in 2008. The hydrological parameters chosen for the simulation are listed in Table 1 . This fitting procedure is one of the limitations of our model; however, sensitivity studies conducted in a prior study [ de Fleurian et al. , 2014 ] showed a reasonable robustness of the model results to significant changes in model parameters.

Boundary conditions for the subglacial hydrological system are straightforward to define at the front of the glacier where the meltwater outlets define the altitude of the water table. There are a number of outlets in our region of interest, which are visible on satellite imagery and were mapped by Lewis and Smith [ 2009 ]. At these locations, the water heads of both the inefficient and efficient systems are fixed at the flotation limit. All the other boundaries of the modeled region are treated as no‐flow boundaries. On the upstream boundary (where the surface elevation reaches 1850 m), this hypothesis is motivated by the fact that our domain extends to the altitude where the bed is assumed to be frozen and is thus not producing meltwater [ Poinar et al. , 2015 ].

The lack of data to describe the subglacial hydrological system makes it difficult to set both the initial and boundary conditions of the model. While we focus our study on the Russell Glacier, we do not limit the model domain to the glacier catchment region, we define a larger zone of interest to make sure that the runoff of neighboring regions, which is routed subglacially to Russell Glacier, is also taken into account (Figure 1 ). The geometry of the glacier surface is taken from the Greenland Ice Mapping Project digital elevation model (GIMP DEM) [ Howat et al. , 2014 ], and the bedrock topography is from the mass conservation method by Morlighem et al. [ 2014 ]. We rely on an unstructured, nonuniform mesh with characteristic element size ranging from 250 to 2500 m with an hourly time step. Ice dynamics and changes in ice geometry are ignored. We are aware that changes in geometry and the feedbacks between ice dynamics and the subglacial drainage system could impact the results of the hydrological model, but we choose to keep the model as simple as possible in order to isolate the effect of the runoff variation on the subglacial drainage system.

As it is less critical for the hydrological system due to its small contribution to the overall water input, the basal water source is treated in a simpler way. The geothermal heat flux is assumed to be homogeneous and constant in time over our study region with a value of 56 mW m −2 , which is consistent with prior studies [ Fox Maule et al. , 2009 ; Dahl‐Jensen et al. , 1998 ]. The other basal source, which is due to frictional heat, is ignored in this study as we are not running a coupled model, which would provide the required values of basal velocities. This second term would increase water production in the fast‐flowing region of the glacier. The induced feedback is still to be studied in detail, but the small amount of water that would be generated in a region where the drainage system is already well developed should not have a significant impact on the overall effective pressure.

We are interested in modeling the effect of varying runoff volume on water pressure at the base of the glacier. Considering this objective, we need a forcing that captures both the intraannual and interannual variability of the runoff. The runoff itself is well reconstructed by regional atmospheric climate models that give both reliable and extensive (in terms of space and time) time series of water production. These inputs, however, need to be preprocessed before being used as source terms in the hydrological model. The water input is here derived from runoff reconstructed by RACMO2.3 (Regional Atmospheric Climate MOdel) [ Noel et al. , 2015 ]. The locations of the input points are determined from the detection of supraglacial moulins and rivers on a Landsat 8 image. More information on the preparation of the forcing is provided in Appendix A .

Water input for the subglacial hydrological system comes from two different sources: a basal source and a surface source. The basal input is due to the geothermal and frictional heat fluxes which, by melting the ice, generate a relatively small amount of water. The surface input is driven by surface runoff, which is routed supraglacially and then englacially to the bed of the glacier. In Greenland, the surface source is the main contributor to the subglacial drainage system on the lower part of the ice sheet. However, in the interior, the basal source dominates. The surface source is more critical for this study as it presents both an intraannual and interannual variability due to the variability of surface runoff, the effects of which are the main interest of this study.

Computation of the size of a subglacial channel involves a widening term due to the melting of the channel walls in contact with water and a closing term due to the creeping of ice into the channel cavity. Using the equation developed by] for the computation of the cross section of a channel and scaling it to our specific geometry and fluxes, we obtain the following equation to describe the evolution of the EPL thickness:whereandare the density and latent heat of fusion for the ice, respectively,is the effective pressure, andandare the parameters of Glen's flow law. A collapsing thickness () is introduced in the model to define a threshold below which the EPL will be deactivated. Introducing this new feature in the model makes it possible to model multiannual cycles of water head variations, which, coupled to simulated water input, enables us to study the interactions between surface runoff and water pressure at the base of glaciers.

More details about this model are provided in de Fleurian et al. [ 2014 ]. For this study, the double‐continuum approach has been implemented in the Ice Sheet System Model [ Larour et al. , 2012 ], along with a new scheme to control the evolution of the thickness of the EPL. This improvement was required in order to allow for the collapse of the EPL at the end of the melt season, which is necessary to perform multiannual simulations. The thickness of the EPL is defined at its activation ( e j value for the EPL given in Table 1 ) and then left free to evolve. The thickness evolution is based on the equations developed for the computation of the size of a subglacial channel [ Röthlisberger , 1972 ; Nye , 1976 ].

The major difference between IDS and EPL, apart from the parameter values and the water input source, is that the IDS is always active, whereas the EPL has a specific activation procedure. The EPL is activated when the local effective pressure (computed as the difference between water pressure and ice overburden pressure) drops to zero. Upon activation, the boundary condition of the EPL region is defined as a no‐flow boundary which, considering mass conservation in the system, tends to propagate the opening of the EPL downstream (i.e., following the IDS gradient). When the opening of the EPL reaches the glacier front, the same Dirichlet boundary condition as the one used for the IDS applies.

The source termis treated in different ways depending on the system considered. Water forcing is applied only to the inefficient component of the drainage system, whereas the source for the efficient drainage system is a transfer term in between the two systems.whereis the forcing input discussed in section 2.2 andis the transfer flux between the two layers, which is driven by the water head difference of the two systems aswhereis the leakage time defining the efficiency of the coupling of the two systems.

Our subglacial hydrological model is based on a double‐continuum approach developed by] to model karstic hydrology and adapted to glaciological applications in]. This approach uses two porous layers to model both the inefficient (or distributed) and efficient (or channelized) components of the subglacial hydrological system (Figure 2 ). The two components of the system are referred to as the Inefficient Drainage System (IDS) to model the distributed part of the system and the Equivalent Porous Layer (EPL) to model the efficient drainage subsystem (Figure 2 ). The computation of the water heads () in the two systems is achieved by solving a vertically integrated diffusion equation:where the subscriptrepresents either the IDS () or the EPL (). The physical characteristics of the two layers, transmitivitiesand storing coefficients, are defined aswhereandare the thickness and conductivity of IDS or EPL, respectively. The other parameters are material characteristics of the porous media (porosityand compressibility) and of water (densityand compressibility), andis the acceleration due to gravity.

To complement these velocity observations, we use the pressure data recorded in a borehole on the lower part of Russell Glacier. The position of this borehole (SHR station on the K transect) is indicated by the yellow triangle in Figure 1 . The details on the sensors and installation are found in Smeets et al. [ 2012 ] and van de Wal et al. [ 2015 ]. The measurements started in July 2010, and we compare the daily mean to the model results.

(a–f) Evolution of the Flotation Fraction (FF, blue left axis) and the detrended GPS measurement of vertical displacement (Vert. Disp., red right axis) for stations S2 to S7 shown in Figure 1 . (g) The comparison between the measured FF (green curve right axis) at station SHR (yellow triangle in Figure 1 ) and the modeled FF at the same point (magenta, left axis) and a location at the same elevation where the EPL is developing (cyan, left axis). (h) Evolution of the integrated runoff (gray left axis) and cumulative integrated runoff (black right axis). The gray zones represent the summer period.

To evaluate the model results, we use an extensive data set of GPS measurements at the surface of the glacier. This data set consists of seven stations (S1 to S7 in Figure 1 ) which were deployed in spring 2009 and recorded the movement of the glacier until spring 2013. Sole et al. [ 2013 ] and Tedstone et al. [ 2013 ] describe the acquisition method. From the full GPS displacement, we use only the vertical displacement recorded at each station. The measurement itself is averaged as a daily mean and linearly detrended to subtract the evolution of the station elevation as it moves downstream. The vertical displacements of stations S2 to S7 are shown in Figure 3 with the reference set at an arbitrary elevation.

Our model suggests that the record melt year of 2012 does not lead to extreme FF values on the lower part of our domain. The summer means are in the range of values modeled for the other years of the simulation on the lower median and upper regions. In contrast, the highest region displays a FF that increases twice as fast as during the entire modeling period.

Year 2011 marks the return of a more average melt year with meltwater levels closer to those of year 2008. This period shows different patterns in the evolution of the FF across the region. In the lower and median regions, the summer FF is back to the level observed in 2008 but the winter mean value is lower than the previous years leading to a smaller mean annual FF. The evolution is rather different in the upper region: the FF mean values are more comparable to the ones modeled during 2010, with only a slight decrease in annual mean FF driven by a decrease of the summer value of the water pressure. Finally, the highest zone shows a continuous increase in FF even during this rather low input melt season.

As stated before, there is no significant trend in interannual or intraannual evolution of the FF on the frontal region. The evolution is more significant on the rest of the model domain. For the period starting at the beginning of the 2008 melt season and finishing at the end of the 2009 winter, the evolution of the FF is the same for all regions with a decrease in annual FF. This decrease is triggered by lower summer FF in the lower and median regions but is mainly driven by the reduction in winter FF in the highest region. The upper region reduction of the annual FF is due to a lowering of both the summer FF and winter FF.

The values presented in Figure 6 are the percentages difference of the FF seasonal means with respect to the MWV averaged over the given regions. The seasons are defined as typical for this region: the summer is from 1 May to 31 August, winter from 1 September to 30 April, and the annual mean is from 1 May to 30 April. The figure only shows the summer mean for 2012 as the simulation does not span the entire winter of 2012.

An examination of the seasonal mean of the FF provides useful insights about the long‐term evolution of water pressure. From the previous results, we define regions with different responses to the runoff forcing: (1) the frontal region (from the ice margin to 400 m of elevation), (2) the lower region (from 400 to 800 m), (3) the median region (from 800 to 1200 m), (4) the upper region (from 1200 to 1600 m), and (5) the highest region (from 1600 m to the top of the domain). The five regions are delineated in Figure 1 .

Performing long simulations makes it possible to investigate the interannual variations of the basal water pressure. To get a more synoptic view of the water pressure evolution throughout the years, we compute the mean value of the FF in given surface elevation bins. Figure 5 presents the FF variations for 50 m surface elevation bins throughout the simulation. The reference for each elevation bin is the mean value of the FF for all the simulated winters (mean winter value, MWV). The pattern of the FF evolution is similar from year to year with differences in the duration and amplitude of the signal. The water pressure remains at or just below the MWV for the beginning of the summer and then increases as meltwater becomes available. The region below 400 m elevation does not show any significant variation in water pressure. Above 400 m, there is a strong seasonal cycle in water pressure. The pressure first increases around 800 to 1200 m elevation, and then the increase in FF spreads downstream as the pressure builds up in the inefficient drainage system and upstream as the melt season goes on and a larger number of moulins feed the subglacial hydrological system. At the end of the summer, the water pressure drops to its winter level. The first place where the FF returns to its winter level is again between 800 and 1200 m elevation, where the development of an efficient drainage system allows to move the water downstream while lowering the water pressure at the base of the glacier. The winter season is marked by a lower FF, but the dynamics of the variability in FF varies across the region. Between the front of the glaciers and 800 m, the FF drops at the end of summer and then stays at the same level for the whole winter. Higher up on the domain, between 800 and 1200 m, the FF quickly drops at the end of summer to get to a low and then increases to get to the MWV at the end of winter. The behavior is reversed at higher elevation (from 1200 to 1600 m) where the FF decreases during all winter. Finally, the highest part of the glacier, above 1600 m, exhibits a different behavior with a moderate decrease of the FF during winter, which is particularly visible after the 2010 melt season.

During the melt season (4 June and 3 August in Figure 4 ) the moulins are feeding the subglacial drainage system, leading to an increase in FF on most of the domain. The FF increase starts on the lower part of the glacier with the emergence of the first moulins at lower elevation. This increase in pressure triggers the development of the EPL at these elevations, where a few well‐developed channel‐like structures are formed to evacuate the water. Later in the melt season, the FF increase extends over most of the model domain with only a small region at the upper limit of the domain, which is not impacted by the meltwater input. Note that the higher‐pressure front is reaching a point higher than the highest moulin. On 3 August, the EPL is at its maximum extent. In this configuration, we observe a few well‐developed channel‐like structures on the lower part of the domain. Higher up, the drainage system is more widespread and less efficient (smaller EPL thickness) over almost all the glacier bed.

Maps of the Flotation Fraction (FF), its difference to the mean winter value (MWV) of the FF, and the thickness of the Equivalent Porous Layer (EPL) at different dates during the 2010 melt season. The mean winter value is defined as the mean winter FF average over all the simulated winters. The FF presents the evolution of the water pressure, whereas the thickness of the EPL represents the development of the efficient drainage system at the base of the glacier. Inactive EPL points are transparent on the third column. The results are overlaid on the same Landsat 8 image as in Figure 1 . See supporting information for an animation of the complete simulation (Movie S1 ).

The FF at different times in the 2010 melt season, along with the differences to the mean winter FF and the thickness of the EPL, is presented in Figure 4 . An animation of the full simulation is available in the supporting information (Movie S1 ). The state of the model at the end of winter (10 April in Figure 4 ) is representative of the mean winter values as shown by the small difference between the FF at this date and the mean winter FF. The FF computed by the model is high at the glacier front, which is forced either by the boundary condition fixing the water head at the flotation limit or small ice thickness. Away from these outlets, the FF drops to more plausible values and then increases with the distance from the front. The FF is higher in large troughs in the glacier bed, which indicates a preferential drainage in these areas. On the highest parts of the domain, the FF decreases as the ice thickness increases, whereas the water pressure is stable due to the lack of water influx at these altitudes. The EPL thickness, which represents the development of an efficient drainage system, shows remnants of the draining system from the preceding melt season. These remnants of an efficient drainage system are likely comparable to an unconnected water system as they have a very small thickness and are not connected to the glacier front.

The evolution of the FF makes it possible to investigate the effect of higher meltwater quantities on the interannual evolution of the water pressure. The cumulative runoff integrated over the entire domain (Figure 3 h, black curve) shows the significant amount of meltwater that is available during the melt season of 2010 and 2012. The runoff integrated over the entire domain (Figure 3 h, gray curve) shows the differences between these melt seasons, 2010 being an average year in terms of amplitude but sustained over a long period, while 2012 had an amplitude higher than normal. These two situations lead to different evolutions for the FF: in 2010, the water pressure reaches a value in the range of the previous years that lasts a long period, whereas in 2012, the FF is higher than the one that was modeled for the previous melt seasons.

To complement these observations, we compare the FF to the measured daily mean FF at the SHR station on the K transect (yellow triangle in Figure 1 ) [ Smeets et al. , 2012 ; van de Wal et al. , 2015 ]. For this comparison, we use the IDS pressure when the EPL is inactive and the EPL pressure when it is active. Comparison between the measured pressure (green curve in Figure 3 g) and the modeled pressure at the same location (magenta) shows a poor agreement due to the fact that no EPL develops in this region according to the model. The comparison with a second point located at the same surface elevation but at a location where the EPL develops (cyan) gives a better agreement with the observations. The steps observed in the modeled FF are due to the pressure difference in the two systems, which appears at the activation and closure of the EPL. The modeled FF where the EPL occurs shows a good agreement with the measured FF in terms of timing of the response. There is, however, a significant difference in amplitude, with a larger amplitude of the signal for the modeled FF.

We first compare the model results (in terms of FF) and the observed vertical displacements at the surface of the glacier. We chose to compare the modeled FF to the detrended daily mean vertical displacement. The vertical displacement is a good indicator of the hydraulic ice bed separation which is directly related to an increase in subglacial water pressure [e.g., Iken and Bindschadler , 1986 ]. This comparison (Figure 3 ) shows a good qualitative agreement between the model and the observed displacements on a multiannual simulation. We note that there is a lag between the model and the measurements at station S7 and a very low response of the model at station S2.

5 Discussion

The application of a new generation, multicomponent, subglacial hydrological model provides new insight into the varying responses of the subglacial hydrological system to the evolution of the meltwater volume in Greenland. The analysis of these distinct subglacial hydrological modes is key to understanding the causes of the different velocity evolutions that are observed on Russell Glacier.

The comparison between the modeled water pressure and the measured vertical surface motion of the glacier shows a good qualitative agreement (Figure 3). The results are consistent with the data for both the seasonal and shorter time scale variations of the signal. The results at stations S2 and S7 are less convincing, which is probably due to various reasons. At station S2, the lack of evolution in the modeled water pressure is probably due to the fact that this station is close to the glacier front and still impacted by the boundary condition that is fixed at this point (Dirichlet). At the uppermost station (station S7), we are more confident in the model results, but, at this point, the ratio between surface uplift and vertical displacement due to the downstream motion of the ice is small, leading to uncertainties in the reliability of the detrended data at this location. The water pressure modeled at the SHR station is not consistent with the observations [Smeets et al., 2012], which was anticipated given the fact that the EPL does not develop at this location in the model. However, the comparison with the data from a grid point where the EPL develops shows a good agreement in terms of timing of the pressure variations. The general shape of the modeled pressure also agrees reasonably well to the measured pressure, even for the late season water pulse observed during the winter of 2010. The amplitudes of the modeled pressure are higher than the measured ones, which is probably due to the rather low pressure that is modeled on the lower part of the glacier. We note that there is still a lack of velocity data during the winter season. Acquiring more data during winter as done by van de Wal et al. [2015] will enable a better evaluation of the models.

The modeling of a large region around the Russell Glacier catchment helps us investigate the characteristics of the drainage system in the region. The snapshots in Figure 4 and particularly the EPL thickness on 3 August give us a good idea about the characteristics of the drainage system when its efficient component is at its maximum extent. On this date, the EPL thickness shows large channel‐like structures extending up to 1200 m in elevation, which is comparable to 41 to 57 km from the glacier front which was defined as the limit of the development of the efficient drainage system by Chandler et al. [2013]. The modeled effective system also shows a different localization from year to year as proposed by Chu et al. [2016]. The displacement of these major drainage pathways can be consulted in the supporting information (Movie S1). Farther upglacier, the EPL is still active but shows a smaller thickness and no channel‐like structure, which is consistent with a less efficient system that would be formed by a larger network of interconnected cavities. Finally, on the highest part of the ice sheet, the EPL is nonexistent, even at the peak of the melt season, which in our model is mainly due to the fact that the water input is too low (no moulin) to allow water pressure to reach flotation level and activate the EPL. Contrary to the results from Chandler et al. [2013], our model shows a developed efficient system, which exists year round near the front of the glacier. This is probably caused by the imposed boundary conditions at the ice margin.

The evolution of water pressure in Figure 5 sheds light on the observed surface velocity and its link to the subglacial water pressure. Observations made in the region have shown an increase in summer velocity during the summers of 2010 and 2012 [Sole et al., 2013; van de Wal et al., 2015], which had a high meltwater input to the glacier. These unusually high velocities can be explained by the broader and higher water pressure peak that are modeled during these years. The model results, however, do not show a lower winter water pressure during these years as hypothesized by Pimentel and Flowers [2010] and observed by Sundal et al. [2011] in the lower part of the ice sheet. The FF evolution during 2010, which is the only high melt season for which we also modeled winter values of the water pressure, yields interesting results because even with a larger melt event, the winter FF is not impacted and is similar to the FF modeled for the winters of 2008 and 2009. This result indicates that the observed lower winter velocities subsequent to extreme summer melt [Tedstone et al., 2015] are probably not entirely due to the response of the subglacial drainage system. Our results, however, are only addressing the evolution of water pressure from a nonfully coupled system. A fully coupled model could yield different results due to feedbacks [Hoffman and Price, 2014; Hewitt, 2013]. Another hypothesis would be the one presented by Andrews et al. [2014] and Meierbachtol et al. [2016] about the importance of hydraulically isolated region of the beds on the variation in late season water pressure.

The results of the model in the highest part of the domain give insight into the processes leading to the distinct behavior observed when investigating the evolution of the ice velocity on this glacier. Above 1600 m, the model shows a constant increase in subglacial water pressure after 2010. This increase cannot be due to a migration of the moulins upstream as these injection points have a fixed location in space throughout the simulation. Our model results indicate that following the high meltwater forcing of the 2010 melt season, the water heads in the moulin region reached a height sufficient to invert the gradient of the hydraulic potential in the highest part of the domain. This inversion allowed the inefficient system to fill‐up without reaching a volume of water necessary to build the water pressure at a level consistent with the development of an efficient drainage system. This phenomenon shows that even if there is no way for the water to reach the bed at these high elevations, the upstream propagation of a pressure wave in the inefficient drainage system has the ability to increase the subglacial water pressure in the interior of the ice sheet. These results provide an explanation for the distinct velocity regimes that are observed on the glacier. On the lower part of the glacier, the pressure builds up during the beginning of the melt season, leading to a speedup. The speedup is hindered by the development of an efficient drainage system that lowers the water pressure to values close to the mean winter value at the end of summer. This pressure drop leads to the return of the glacier velocities to their winter values. On the other hand, the buildup of pressure on the highest part of the glacier is not followed by a drop in pressure as the efficient drainage system does not develop at these altitudes.

This study has allowed to differentiate two regions in the drainage dynamics of the Russell region. In the lower part of the glacier, the water pressure is driven by an efficient drainage system which develops early in the melt season, decreases subglacial water pressure, and hence produces the surface velocities observed by Sole et al. [2013] and Tedstone et al. [2013]. The subglacial water pressure on the upper part of the glacier, however, is mostly driven by an inefficient drainage system in which the pressure builds up from one melt season to the next one. This increase in subglacial water pressure could be the source of the acceleration in surface velocity observed by Doyle et al. [2014] on the higher part of Russell Glacier.

The mechanism explained above still leaves uncertainties about the long‐term effect of an increase in meltwater production in Greenland. Our results suggest that the water pressure on the highest part of the glacier will continue to build up if the meltwater input is sustained at its current level. This would have the effect of maintaining the interannual speedup, which is currently observed. If this process is sustained for a certain amount of time, it could build up a subglacial pressure wave with the potential to destabilize the highest part of the ice sheet as was proposed by Zwally et al. [2002]. Another possibility is that the increase in pressure at higher elevation leads to the development of an efficient channel system at high elevation. This development would lead to a mitigation of the interannual acceleration, and the upper part of the glacier would exhibit a more marked seasonal evolution similar to the one presently occurring below the ELA. This would necessitate the development of an efficient drainage system under thick ice, which is unlikely [Dow et al., 2014] if the input of water at higher elevation is not sustained by large and sustained meltwater production.

The main limitation of our model is in the choice of its parameters that we cannot measure in the field and hence need to be fitted. Sensitivity studies performed in earlier experiments [de Fleurian et al., 2014] show that the model gives consistent results in terms of variations in the timing and amplitude of water pressure when different parameter sets are used. Another problem that arose in this study is the impact of the ice margin conditions on water pressure. However, this impact seems to be limited to the lowest part of the glacier tongues, which are not very extensive and do not seem to impact the results on the remainder of the domain. The coupling of this model to an ice dynamics model should be an important objective for future studies in order to fully ascertain if the hypothesis presented herein remains valid in the case of a fully coupled system.