Impact of runoff temporal distribution on ice dynamics.

. Records of meltwater production at the surface of the Greenland ice sheet have been recorded with a surprisingly high recurrence over the last decades. Those longer and/or more intense melt seasons have a direct impact on the surface mass balance of the ice sheet and on its contribution to sea level rise. Moreover, the surface melt also affects the ice dynamics through the meltwater lubrication feedback. It is still not clear how the meltwater lubrication feedback impacts the long term ice velocities on the Greenland ice sheet. Here we take a modelling approach with simpliﬁed ice sheet geometry and climate 5 forcings to investigate in more detail the impacts of the changing characteristics of the melt season on ice dynamics. We model the ice dynamics through the coupling of the Double Continuum (DoCo) subglacial hydrology model with a shallow shelf approximation for the ice dynamics in the Ice-sheet and Sea-level System Model (ISSM). The climate forcing is generated from the ERA5 dataset to allow the length and intensity of the melt season to be varied in a comparable range of values. Our simulations present different behaviours between the lower and higher part of the glacier but overall, a longer melt season will 10 yield a faster glacier for a given runoff value. Furthermore, an increase in the intensity of the melt season, even under increasing runoff, tends to reduce glacier velocities. Those results emphasise the complexity of the meltwater lubrication feedback and urge us to use subglacial drainage models with efﬁcient drainage components to give an accurate assessment of its impact on the overall dynamics of the Greenland Ice sheet. speciﬁc combination of time-steps gave consistent results while keeping the computational cost at a manageable level. The management of the different time-steps for the coupling is performed through the averaging of the effective pressure over the length of the ice ﬂow time-step and the averaged effective pressure is then used to compute the ﬂow. The geometry of the ice model needed as an input to the hydrology model is then kept ﬁxed for the four sub-steps of the subglacial hydrology 100 model. originally leads to a higher maximum velocity at low elevations but then translates to a lower mean summer velocity over the whole domain and an overall slowdown of the glacier. This behaviour is mostly due to the development of a very transmissive efﬁcient drainage system at low elevations which keeps the effective pressures at a high level throughout the melt season, and so leads to relatively low summer velocities when averaged over the whole domain. This study can not give a deﬁnite answer on the impact of an increase in runoff on the long term evolution of glacier velocities. Indeed, the impact of the melt season amplitude is radically different at different altitudes and these regions are delimited by the existence or not of an efﬁcient 450 drainage system. Our experiments show that an increase in the amplitude of the melt season leads to a larger extent of the region where the drainage is controlled by the efﬁcient subglacial drainage system. This change of regime on the highest regions of the glaciers might signiﬁcantly alter the velocity proﬁle of the glacier with regions at higher altitudes then experiencing spring speed-up events. This ﬁnding emphasises the fact that subglacial drainage models with efﬁcient drainage components should be used if one wants to give an accurate assessment of the effect of meltwater on the overall dynamics of the Greenland Ice sheet.


15
Since the 2000's a large number of studies have pointed towards large increases in the amount of melt recorded at the surface of the Greenland ice sheet (e.g. Steffen et al., 2004;Mote, 2007;Hanna et al., 2008) This is reflected by the fact that during the first two decades of the twenty-first century we have seen the melt record over Greenland beaten four times already (Nghiem et al., 2005;Mote, 2007;Mernild et al., 2009;Tedesco et al., 2011Tedesco et al., , 2013aTedesco and Fettweis, 2020). Ahlstrøm et al. (2017) identified in southwest Greenland that a shift in the runoff regime took place in 2003 with an 80% increase in runoff for the 20 following decade compare to the period 1976-2002. The changes in the melt season are clearly observed in the distribution of the melt (Zwally et al., 2011;Sasgen et al., 2012) culminating in the 2012 season when the whole surface of the Greenland ice sheet experienced melt at some point during the year (Nghiem et al., 2012;Tedesco et al., 2013a). But even if it is less visible, the length of the melt season has been increasing since the late 70's (Colosio et al., 2021) and that lengthening has a large and Werder (2018) showed that it could account for a volume loss significantly higher than what was estimated from the simple parameterisation of Shannon et al. (2013); Fürst et al. (2015). 60 We will first give an overview of the component of the model which are specific to this study before describing in more details the set-up and forcing that are used in the study. We first focus on the results from the reference simulation which give us the opportunity to assess the behaviour of the system before the impact of the different forcings are presented. Finally we give a broader interpretation of the results of our experiments in light of recent findings on the meltwater lubrication feedback.

Model
In order to investigate the impact of meltwater availability on ice dynamics we perform coupled subglacial hydrology and ice dynamics simulations within the Ice-sheet and Sea-level System Model (ISSM Larour et al., 2012). Within the ISSM, the ice flow is treated following the Shallow Shelf Approximation (SSA) (Morland, 1987;MacAyeal, 1989). The choice of this approximation is motivated by a need to run on relatively short time-steps, and so it is necessary to have a relatively 70 computationally-cheap ice flow model. Since our interest is in the sliding of the glacier (rather than its deformation) SSA is well suited. The lubrication feedback and the impact of subglacial water on ice flow dynamics depends on the choice of the sliding law linking the ice dynamics to the subglacial hydrology. Here we use a non-linear friction law as described by Schoof (2005) and Gagliardini et al. (2007) which links sliding velocities (u b ) and basal shear stress (τ b ): (1) 75 Where the parameters A s and C are the sliding parameter in the absence of cavities and Iken's bound parameter respectively.
Iken's bound (Iken, 1981) represent the maximum value that can be taken by τ b /N and is only determined by the maximum up-slope of the bed. N is the effective pressure which is produced by the subglacial hydrology model in response to the surface melt forcing. The subglacial hydrology model used is the Double Continuum (DoCo) approach as described in de Fleurian et al. (2016). This model is based on a double sediment layer system. The first layer, with a low conductivity (K s ) represents 80 the inefficient drainage system (IDS). The second layer, with a much higher conductivity (K e ) is only activated when the local effective pressure is equal to 0. Once activated, the thickness (e e ) of this efficient drainage system (EDS) evolves from its initial thickness following equations describing the size of subglacial channels (Röthlisberger, 1972;Nye, 1976) and scaled to take into account the specific geometry of the EDS.
∂e e ∂t = gρ w e e K e ρ i L (∇h e ) 2 − 2An −n N n e e ,  ice. The other term represents the closing of the system through ice creep where A and n are Glen's parameter and exponent.
As the pressure in the EDS decreases, it will get thinner until the point where its transmitivity (T e = K e × e e ) is lower than that of the IDS at which point the EDS is deactivated.

90
The parameters of the model are given in Table 1, and detailed information about the model can be found in de Fleurian et al. (2016) andde Fleurian et al. (2014).
One challenge when running a coupled ice-flow subglacial hydrology model is that the two systems are responding on different timescales. The subglacial hydrology model needs short time-steps to achieve stability and provide reliable results while the ice flow model can run on a longer time-step, which saves on computing time. In this study we used a 15 minutes time-

95
step for the hydrology model and a one hour stepping for the ice flow model. Testing of a number of different options showed that this specific combination of time-steps gave consistent results while keeping the computational cost at a manageable level.
The management of the different time-steps for the coupling is performed through the averaging of the effective pressure over the length of the ice flow time-step and the averaged effective pressure is then used to compute the flow. The geometry of the ice model needed as an input to the hydrology model is then kept fixed for the four sub-steps of the subglacial hydrology 100 model.

Geometry and spin-up
The geometry of the system presented on Fig. 1 is kept as simple as possible in order to obtain a consistent response to applied perturbations. The domain on which the simulations are performed is a synthetic representation of a land terminating ice sheet margin. The glacier is 150 km long and 20 km wide with a flat bedrock of elevation z b = 465 m and an initial surface elevation 105 z s defined by a parabolic function: The parabolic function is chosen such that it resembles the topography of West Greenland as used in the ERA5 dataset. From this geometry, a long spin-up is run to achieve stability of the different components of the system. This involved an offline coupling of the ice sheet model where the surface mass balance was given by the reference forcing described in Sect. 2.3 and 110 the subglacial hydrology model. A relatively stable state with only a small volume loss was attained after roughly 5000 years of simulation. This final state is then used for all the perturbation simulations.

Forcing
Here we take an idealised forcing using the method first described by Hewitt (2013). In this formulation the temperature is first defined at a reference height using the equation: where -T ref is the temperature at the reference elevation of 465 m at the front of the glacier [°C]. r m is the positive degree day at the reference elevation [°Cdays].

120
t spr is the beginning of the melt season (day 100) -∆ t is the length of the initiation of the melt season [days].
The runoff itself is then computed using a given lapse rate where r s is the lapse rate in°Km −1 , ddf is the degree day factor or conversion rate from temperature to runoff.

125
The three parameters of this model that we chose to test the sensitivity are the length of the melt season (∆ m ), the PDD at the reference elevation (r m ), and the length of the initiation (∆ t ) as shown on the inset of Fig. 1. We use the ERA5 reanalysis data to derive a realistic positive degree day (PDD), season length and lapse rate for our model (Hersbach et al., 2020). First we extracted the daily surface air temperature from 1979 to 2018 for south-western Greenland at a fixed latitude and for a longitude band going from close to the coast (465 m above sea level) to near the highest point of the land at 2256 m above 130 sea level (67°N, 45-50°W). To calculate the length of melt season for each year we smoothed the daily temperatures using a 15-day midpoint running mean. This smoothing was applied to avoid creating anomalously long melt seasons due to a single day of high temperatures occurring far outside of the rest of the melt season of a given year. The length of the melt season was then calculated using the first and the last day of the year with daily-mean temperatures greater than 0°C. The lapse rate for each day was calculated by taking the gradient of a least-squares linear regression of temperature against elevation across the 135 domain. The daily lapse rate was then averaged for the every day where the temperature was above 0°C at at least one grid point, to obtain a typical lapse rate for the melt season. To calculate the positive degree days (PDD) for each year we simply summed the positive temperatures in Celsius over all n levels of the ERA5 reanalysis data within our domain: This procedure gave us a single representative value for the PDD, length of the melt season, and the typical lapse rate for each 140 year from 1979-2018. We then fit a normal distribution to the values for all years to determine the mean and standard deviation (SD) of the inter-annual variability in each of these variables. These were used to define low (mean-SD), medium (mean), and high (mean+SD) values for our sensitivity studies. We also tested the steepness of the temperature variation at the beginning and end of the melt season, which is controlled by the parameter ∆ t and called initiation further on in the manuscript (Fig. 1).
Since this cannot be derived from ERA5 data we instead chose these values pragmatically to create a spread of reasonable 145 seasonal cycles. The different initiation and melt season length are given in Table 2. The value for r m is presented there with the ratio r m /∆ m which represent the maximum temperature at the reference elevation.
Regardless of the length of the melt season we define the summer as the period between days 100 (t spr , April 10 th ) and 241 (t spr + ∆ m ) for the reference simulation, August 29 th of the simulation.  To overcome this distribution issue, we decided to use the Wilcoxon signed-rank test to investigate the difference of every parameter set to the reference simulation. We use a non parametric Wilcoxon signed-rank test to assess the differences between the reference simulation and the perturbed simulation 160 and test if the perturbation lead to significant differences in the response of the model. Every time that the differences from one simulation with respect to the reference is said to be significant it means that the null hypothesis of the Wilcoxon signed-rank test (the median of the differences is zero) is rejected with a 1% confidence level. These higher velocities are due to the continuous lowering of the effective pressure throughout the melt season. These spikes 180 are linked to a secondary reduction in the effective pressure which happens when the EDS is fully developed and so the water going into the system tends to overload this system. These accelerations have been observed in Greenland where it is usually linked to more important melt or rainfall events after the development of an efficient drainage system (e.g. Cowton et al., 2013;Doyle et al., 2015;van de Wal et al., 2015). We will further refer to this late acceleration phase as the autumn acceleration. At the end of the melt season, the drop in runoff leads to a fast increase in the effective pressure which goes back to its winter 185 level. At the lower elevations ( Fig. 3g-h), we see an overshoot of those winter values as the subglacial drainage capacity is higher than the recharge at this time. is driven by downstream activity as there is no runoff at these elevations. As a result, the effective pressure shows very small variations which get more and more out of phase with the melt season as we go higher up the glacier (Fig. 3j). We will in the following use the acceleration of the glacier as a reference metric to define the elevation at which the glacier shifts from the spring speed-up driven velocity pattern to the more gradual one. This shift in behaviour (henceforth called SSU max ) is set We can also observe the two different behaviours of the system in Fig. 4. Here we see that the fast spring speed-up is confined bellow the SSU max elevation (white line on all panels) as expected from its definition (Fig. 4b). The EDS transmitivity (Fig. 4d) 205 is a proxy for the capacity of the subglacial drainage system to drain water. The white region corresponds to periods where the EDS is not active, while the higher values indicate a highly developed and efficient system. The lowering of the transmitivity at  The runoff is given in panel (a) for reference. The white line across panels (b) to (d) represents SSUmax as defined in Sect. 3.1 the end of the melt season shows how a drop in water pressure leads to the contraction of the EDS and ultimately to its collapse and deactivation. From the EDS transmitivity we can define the maximum altitude at which the EDS is activated (EDS max ). EDS max is slightly higher than SSU max and is the effective boundary between the lower part of the glacier where the effective 210 pressure is mostly controlled by the efficient drainage system, and the upper part of the glacier where the effective pressure evolution shows more gradual evolution in line with the weaker conductivity of the inefficient drainage system.

Melt Season Length Forcing
We first investigate the effect that a shorter or longer melt season has on the glacier's velocities. Since we chose to keep the runoff constant for this set of simulations, changes in the length of the melt season simultaneously impact the melt intensity.

215
This choice allows to compare both length and intensity of the melt season and investigate which parameter has the stronger effect on the glacier's velocity. In Sec. 3.3 we will present more details about the specific impact of a change in melt intensity or melt season length when they are applied separately. Figure 5 presents the comparison between the reference simulation (grey), the long and low intensity melt season (red) and the short and high intensity melt season (blue) is shown for two different altitudes.

220
Starting with the velocities we see quite different evolution for the long and short melt seasons with a distinct evolution above and bellow the SSU max elevation. We must note here that SSU max also varies with the different forcing ranging from (c) (d)  996 m for the longest melt season to 1100 m for the shortest one as reported in Table 3. At low elevation, the shorter and more intense melt season leads to a sharper, shorter lived, but also more intense spring speed-up (Fig. 5a). This differs with the evolution modelled for the longer melt season where the spring speed-up is less marked initially and followed by a second At higher elevations, the changes in the length and intensity of the melt season leads to contrasted evolution for the effective pressure. For the short melt season, the response is similar to the one that was observed at lower elevation. Indeed, at this elevation the runoff is quite substantial leading to the rapid decrease of the effective pressure at the beginning of the melt season. The effective pressure then reaches values that are low enough to activate the EDS as seen in Fig. 5f-g, which in turn  Table 3 summarises the impact of the length of the melt season on the dynamics of the glacier. The shorter and more intense melt season shows a significant reduction in the overall glacier velocity which is mostly driven by a reduction of the summer velocities below the SSU max altitude (1071 m in this case). In the case of the long melt season, the acceleration of the glacier is mostly driven by the summer acceleration of the lowest parts of the glacier. There is a slight deceleration at higher elevations but that is not sufficient to offset the doubling of the mean summer velocities that is observed during summer at the lowest 260 elevations. The winter mean velocities of the longer melt season experiments are somewhat deceiving as the longer melt season means that part of the melt is now happening during the so called winter. Coming back to Fig. 5a it is clear that the winter  where the shorter and more intense melt season leads to a slower glacier at low elevation and faster at higher elevations where 265 the reversed response is observed for a longer, less intense melt season. In both cases however it is the summer velocities at low elevations which are driving the changes in the annual velocities of the glacier.

Intensity vs length
To discriminate between intensity and length of the melt season the requirement that the runoff is to be equal in all simulations has been released. First we run a series of experiments with the same intensity but in which the length of the melt season varies 270 (Fig. 6). Second, we keep the melt season length fixed, but vary the intensity (Fig. 7).
As expected the evolution at the beginning of the melt season of those experiments is very similar as they are all experiencing the same forcing. The simulations start to differ at the point when the short melt season ends. At this time, the effective pressure for this simulation goes back to its winter level with a slight overshoot (Fig. 6c). This timing also coincides with a period during which the effective pressure of the two longer melt seasons are dropping again. This decrease seems to be due to the efficiency 275 of the drainage system reaching a maximum before slowly decreasing and driving the effective pressure down (Fig. 6e). The effect of this effective pressure decline is transferred to the ice velocities where we see an accelerating trend at the end of summer for the simulations with the longer melt seasons (Fig. 6a). At higher elevations the overall behaviour of the glacier (c) (d)  shortest melt season finishes the melt season the effective pressure goes back to its winter levels, while the simulations with 280 longer melt seasons have a continued decrease of the effective pressure (Fig. 6d). In terms of velocity, only the duration of the summer acceleration is significantly different under the different forcings, with a mean summer velocity that is comparable for all simulations (Fig. 6b). However, the extreme values for the longer melt seasons tend to show more important acceleration events happening at the end of the melt season. The overall changes in velocities are presented in Table 4, contrasting with the experiments performed with a constant PDD, 285 we see here that the changes in the melt season length act in the same direction (either acceleration or slow down) at all altitudes. We also note that the symmetry in the forcing difference leads to a linear response in term of velocities where the annual mean velocities for the shorter melt season forcing are decreased by a similar amount with respect to the reference simulation and vice versa.
Comparing the simulations with different intensities yields more significant differences between simulations (Fig.7). We find 290 here a similar pattern to the simulations with varying intensity and duration (Sec. 3.2). At low elevation the initiation of the melt season yields the same pattern with a strong spring speed-up which is followed by a secondary acceleration in the case of the low intensity forcing (Fig.7a). The higher intensity melt season also shows a marked overshoot of the winter effective pressure value at the end of the melt season (Fig.7c). Above SSU max the difference in the melt intensity means that the experiments with the lower intensity do not experience runoff at this elevation. This leads to the effective pressure in this case being driven by the 295 downstream evolution of the effective pressure which results in a slow decrease of the effective pressure throughout the melt season (Fig.7d). The velocity response follows the variations in effective pressure with a slow increase towards a rather low maximum summer velocity towards the middle of the melt season (Fig.7b). Despite showing some recharge at this elevation, the reference simulation shows a similar velocity pattern as the meltwater availability at this altitude does not lead to very low water pressures. The response to the more intense (warmer) melt season is quite different. Here the water input is sufficient to 300 drive a quick drop in effective pressure and the following rebound once the efficient drainage system is activated (Fig.7f). In (c) (d)  The general impact of the melt season intensity on the glacier's velocity shows counter-intuitive results ( Table 5). The mean velocities over the whole domain indicate here that the stronger intensity forcing does not lead to any significant difference in 305 velocity compared to the reference simulation. However, the colder (lower intensity) melt season shows a substantial increase in the mean annual velocities. This acceleration is mostly due to a two fold increase of the summer velocities at low elevation, which in our experiments can not be offset by the slowdown at higher elevation or during winter. The increase of the melt area that is driven by the larger intensity of the melt season drives an acceleration of the upper part of the glacier, but this acceleration is not sufficient to impact the overall velocity of the glacier.
310 Figure 8 summarises the effect of intensity and length of the melt season on both the subglacial hydrology system and the glacier dynamics. In Fig. 8a we see a clear linear relation between the increase in the runoff volume and the transmitivity of the efficient drainage system, that shows that on the whole glacier, an increase in meltwater volume leads to a more developed subglacial drainage system either in areal extent (increase of EDS max ) or transmitivity of the system. This relationship is not related to the distribution of the melt throughout the year, and we see that an increase in runoff driven by either a longer melt 315 season (growing marker sizes) or a intensifying melt season (darkening markers) lead to a similar increase in transmitivity.
The response is more contrasted for the areal development of the efficient drainage system. This characteristic is described by EDS max (Fig. 8b), the maximum altitude of the glacier surface under which the efficient drainage system is active. As for the transmitivity, Fig. 8b shows that the efficient drainage system develops further upstream if the increase in runoff is due to an increase in the intensity of the melt season. However, EDS max shows a more gradual increase if the runoff is only 320 increased by a lengthening of the melt season. A result of these two different behaviour is that for a given runoff, the spread of the efficient drainage system will be greater for a short and intense melt season (small black dot of Fig. 8b) than for a long and colder melt season (big light grey dot of Fig. 8b). The interplay between these two relations leads to a complex relationship for the mean velocity of the glacier as seen on increase in runoff and so a longer melt season will lead to an increase in the glaciers velocity. However, for a fixed season 325 length an increase in the melt intensity first drives a sharp decrease in velocity followed by a slow velocity increase. Focusing again on a constant runoff, the velocities in our model are increasing if the length of the melt season increases. That is the reverse scenario of the one of EDS max which was expected. Indeed, a more widespread efficient drainage system allows to drain away a larger amount of water and in the end allows the velocities during summer to settle at a lower point than on a glacier with a lower effective pressure.

Initiation Length Forcing
The definition of our forcing means that a change in amplitude would lead to small variations in the steepness of the initiation of the melt season. In order to evaluate the effect of this change in the recharge increase, here we present a set of experiments where the length of the initiation has been changed which leads to initiations with different rates. Figure 9 shows the evolution for simulations in which the initiation rate is altered from the common reference. These simulations have a common runoff 335 which means that the simulation set with the shorter initiation has the lowest melt intensity.
At low elevations, the simulations with the steeper initiation show a very sharp and short-lived spring speed-up Fig. 9a. That contrasts with the simulation with a more gentle initiation where the spring speed-up reaches lower peak velocities and has a longer duration with a secondary peak later in spring. Figure 9a also shows that with a steeper initiation the velocities after the initial speed-up are relatively slow which contrasts with the more gentle speed up where a number of the ensemble members 340 show quite strong velocities throughout the melt season. These two contrasting responses can be explained by the development speed of the efficient drainage system as seen in Fig. 9e. The simulations with steep initiation periods show the fast development of a very efficient drainage system, whereas the simulations with longer initiation periods show a later development of a less efficient drainage system. This difference explains why the effective pressure in the case of the steep initiation simulation shows a very steep rebound to a rather high summer value for the effective pressure, while the effective pressure for a more gentle 345 melt season initiation stays rather low throughout the summer, which drives the observed fast velocities. The development of a rather inefficient efficient draining system in the case of the gentle initiation also leads to a large migration upstream of this system and its rather low efficiency produces rather low effective pressures on the higher part of the glacier.
We see on Table 6 that the effect of a change in the initiation length only impacts significantly the lower part of the domain where the efficient drainage system controls the effective pressure response. The domain means show that the gentler initiation 350 has a larger impact on the velocities with an acceleration of roughly 15% where the steeper initiation only drives a 6% decrease in ice velocity.
It is interesting here to compare these results with the ones from the preceding experiment Fig (7c. In Fig. 9c the simulation with the lowest intensity (sharper initiation) shows a behaviour that is closer to the one of the highest intensity (which is also the one with the sharper initiation) of Fig. 7c, where the drop in effective pressure is quickly followed by a fast rebound to 355 rather high summer velocities. (c) (d)  This reinforces the hypothesis that the rate of recharge of the subglacial drainage system is might have a similar or even a larger impact than the volume of water that is actually injected into the system (Bartholomaus et al., 2008;Hoffman et al., 2011;Bartholomew et al., 2012;Cowton et al., 2016). As with any modelling study, the results presented here might be impacted by the model design and the experiment set-up. It is important to note that subglacial hydrological models have not converged on a standard way to treat the subglacial drainage system yet and the SHMIP exercise (de Fleurian et al., 2018) has shown that there is some discrepancies between the different approaches. However, we think that the results presented here are robust in this model and have a physical explanation which is coherent with existing theories and observations of subglacial drainage.

365
One specific shortcoming in our experiments is that the recharge of the subglacial drainage system has the same distribution and timing as the surface runoff. This is not what is expected in a natural setting where water produced at the surface will transit at the surface of the ice for a given time before going through the ice through moulins and entering the subglacial drainage system through these localised injection points. Scholzen et al. (2021) showed a limited impact on the subglacial drainage system between simulations lead with a spatially homogeneous or discrete basal recharge. Their results show a difference in 370 the timing of the effective pressure response to the runoff but only slight variations in amplitude or length of the pressure pulse.
This is coherent with the results of the SHMIP intercomparison (de Fleurian et al., 2018), it seems that the moulin location does not have a large impact on the overall evolution of the subglacial water pressure. It seems however that the highest injection point on the domain could alter the model results, in this study this highest injection point is defined as the highest point where melt exists on the glacier. However, if the supra and intraglacial drainage was consider it is most likely that the highest 375 injection point would not change on a yearly basis to adapt to the changes in the melt season. Gagliardini and Werder (2018) has modelled moulin migration rates from roughly one to ten meters of elevation per year depending on the setting.
The results of our model suggest that the relationship between ice velocity and meltwater runoff is also strongly influenced by the distribution of the melt during the year. Hence if an increase in runoff is driven by a longer melt season the mean annual 380 velocities of the glaciers will show a strong increase at all elevations. However, an increase in the intensity of the melt season first drives a reduction of the ice velocities before they start to increase again with what seems to be a smaller rate. The mean annual velocities can however not explain the complexity of the lubrication feedback. The effect of a change in the length of the melt season has a similar effect on all regions of the glacier, and the velocity increases linearly with the length of the melt season (Table 4). This differs from the response to an increase in the intensity of the melt season. There the response is different at low 385 elevations, where the subglacial drainage is controlled by the efficient component, and at high elevations, where the drainage is controlled by the inefficient components. At higher elevations, our results compare well with those of Gagliardini and Werder (2018) where we see that an increase of the recharge in the regions controlled by the inefficient drainage system leads to an acceleration of the glacier. However, at lower elevations our results diverge from preceding studies. In the lower parts of the glacier, where the subglacial drainage is controlled by an efficient drainage system, an increase in the intensity of the melt 390 season leads to a decrease in the ice-flow velocities. We explain this result through the evolution of the efficient drainage system in those simulations ( Fig. 7e-g). In these figures we can see a faster development of the efficient drainage system to higher elevations for the more intense melt season. The development of this system leads to an increase in the effective pressure on the lower part of the glacier, and so keeps the velocities at levels that are comparable to the ones of our reference simulations. In the case of the low intensity melt season however, the water recharge is not sufficient to trigger the development 395 of a well-developed efficient drainage system, that leads to low effective pressures throughout the melt season which in turn induce a large increase in velocities compared to our reference simulation. It must be noted however that these conclusions only hold on seasonal velocities and that a higher intensity runoff will lead to higher maximum velocities. This result can be compared to the observations made on lake drainage by Tedesco et al. (2013b) where a fast draining lake triggered a large and short lived speed-up while a slower draining lake only generated a mild speed-up after which the velocities stabilised at a 400 higher level than what was recorded before the lake drainage. We can see the effect of the rate of recharge on the experiments where the length of the initiation of the melt season was changed (Fig. 9). In this case, the variation in the duration of the initiation of the melt season leads to quite large variations in the overall velocities of the glacier. Here, the sharper initiation leads to slower velocities, but this simulation also has a slightly smaller amplitude to keep the overall runoff identical for all simulations. This is in contradiction with the the preceding experiments where the smaller amplitudes were driving a faster 405 glacier. But the similarities here lie in the sharpness of the runoff curve at the begining of the melt season. Hence a sharp rise of the temperatures at the begining of the melt season has the potential to produce a large amount of meltwater which in turn could trigger an early activation of the subglacial drainage system and lead to gentler velocities throughout the summer (after a large spring speed-up event). This large impact of the slope of the temperature rise at the begining of the melt season is problematic to provide estimates of the impact of the lubrication feedback as this parameter is highly variable and complex to characterise 410 in the existing dataset.
In our model, the observed mean velocities are mainly driven by the lower regions of the glaciers where the velocities are significantly higher. That is clear when comparing the velocity evolution at different altitudes to the domain mean velocity on The large difference in behaviour between the lower and higher part of the glacier make it challenging to extrapolate the results to a longer term evolution of the ice dynamics without performing the actual simulations. The two different behaviours are tightly linked to the region in which the efficient drainage system develops which itself is quite sensitive to the intensity of 420 the melt season as seen on Fig. 8b. An increase in the length of the melt season should not alter in a large way the spread of the efficient drainage. Moreover, the variation in melt season length is affecting the whole glacier in the same way so that the response here is quite clear and longer melt season should lead to overall faster glaciers. However, it is not expected that the current evolution in climate would only alter the length of the melt season in Greenland and our model shows that the impact of lengthening the melt season is actually only one third of the acceleration that we observe when we reduce the temperature 425 by a comparable amount. This shows that at least in our model, the effect of the intensity of the melt season is more marked than its length. However, this second effect as larger implications on the migration of the efficient drainage system upstream.
In our experiments, an increase in runoff amplitude leads to a slight slow down of the glacier which is linked to the larger spread of the efficient drainage system. The large runoff at higher altitudes also drives a faster glacier there and in the long term that could induce a larger region of faster flow which would lead to faster velocities overall. On the other hand, a longer melt 430 season of lower intensity drives a faster glacier in the marginal region, but also shows a decrease in the upstream velocities.
This should however be taken with caution as a lowering of the glacier surface would lead to an intensification of the melt which we have shown to be a potential negative feedback on ice velocities. To better understand the impact of these processes on future Greenland Ice Sheet velocities, studies could focus on finding out the main expected trends in length and intensity of future melt seasons.

435
These different scenarios, with counter-intuitive results for the seasonal velocities, show that the full subglacial drainage system, and particularly its efficient component, must be included in studies that aim to quantify the effect of meltwater lubrication feedback.

Conclusions
We developed a set of experiments that allows the comparison of the effect of different parameters impacting the distribution 440 of runoff throughout the melt season. The use of forcing scenarios based on the ERA5 reanalysis dataset gives us confidence that the variations in intensity and length of the melt season that we tested here are representative of the range of existing melt seasons. Our results show that under a constant runoff, an increase in the length of the melt season drives an overall faster glacier. With simulations spanning a range of different runoff intensities we show that an increase in melt intensity originally leads to a higher maximum velocity at low elevations but then translates to a lower mean summer velocity over the 445 whole domain and an overall slowdown of the glacier. This behaviour is mostly due to the development of a very transmissive efficient drainage system at low elevations which keeps the effective pressures at a high level throughout the melt season, and so leads to relatively low summer velocities when averaged over the whole domain. This study can not give a definite answer on the impact of an increase in runoff on the long term evolution of glacier velocities. Indeed, the impact of the melt season amplitude is radically different at different altitudes and these regions are delimited by the existence or not of an efficient sheet.
Code and data availability. The Ice-sheet and Sea-level System Model is freely available at https://issm.jpl.nasa.gov/. The model outputs coresponding to this study are available on zenodo https://doi.org/10.5281/zenodo.5959181 Author contributions. Experiment design was collectively done by all co-authors. RD designed the forcing. BdF ran the experiments and wrote the manuscript with input from all co-authors.