Dynamic response of the Greenland ice sheet to recent cooling

The subglacial hydrological system critically controls ice motion at the margins of the Greenland Ice Sheet. However, over multi-annual timescales, the net impact of hydro-dynamic coupling on ice motion remains poorly understood. Here, we present annual ice velocities from 1992–2019 across a ~10,600 km2 land-terminating area of southwest Greenland. From the early-2000s through to ~2012, we observe a slowdown in ice motion in response to increased surface melt, consistent with previous research. From 2013 to 2019 however, we observe an acceleration in ice motion coincident with atmospheric cooling and a ~15% reduction in mean surface melt production relative to 2003–2012. We find that ice velocity speed-up is greater in marginal areas, and is strongly correlated with ice thickness. We hypothesise that under thinner ice, increases in basal water pressure offset a larger proportion of the ice overburden pressure, leading to reduced effective pressure and thus greater acceleration when compared to thicker ice further inland. Our findings indicate that hydro-dynamic coupling provides the major control on changes in ice motion across the ablation zone of land terminating margins of the Greenland Ice Sheet over multi-annual timescales.


Results
Spatial pattern of velocity change. Ice velocity displays a clear slowdown across the study site ( Fig. 1A) between the periods 1992-2003 and 2003-2012, with 93.1% (9895 km 2 ) of the region exhibiting reduced velocity in the latter period, and a mean regional slowdown of ~15.3% ( Fig. 2A). The slowdown is greatest at lower ice thicknesses and decreases in magnitude inland as surface elevation and ice thickness increase (Fig. 2C), although deceleration characterises all ice thicknesses assessed and extends further inland than in previous work 40 . We would expect this to be the case as at higher elevations further inland, lower surface melt rates, thicker ice and shallower surface slopes slow channel growth, allowing subglacial water pressures to remain higher in smaller subglacial channels 26,49,50 .
Subsequently, we observe a mean region-wide acceleration of ~7.9% between the periods 2010-2012 and 2017-2019 (Fig. 1B), with 89.6% (8218 km 2 ) of this second change map exhibiting increased velocity during the latter period (Fig. 2B). This acceleration in ice motion is proportionally largest at lower ice thicknesses and decreases as ice thickness increases but is observed across all ice thicknesses studied (Fig. 2D).

Ice velocity and surface melt/runoff time series.
To investigate changes in ice motion (Fig. 3A), we follow recent work by assessing anomalies in annual velocity (see SI) whereby only pixels common to all the velocity fields presented in the time series are used in the computation of the median velocity anomalies 50 . Most of the common points fall between 600 m and 1000 m elevation (a.s.l.) due to the hypsometry of our study region (Fig. 4), although we observe >50% coverage in each 100 m elevation band between 100 m and 900 m. In total, we compute velocity anomalies across 71703 common measurements (4130 km 2 ), an order of magnitude greater than previous work 40 (their Fig. 2C).
Over the time period from 1992 to 2019, the ice velocity anomaly time series reveals a mean slowdown of −0.72 ± 0.08 m yr −2 (R 2 = 0.70, p < 0.01). The quality of the fit is further improved via break point analysis (see SI) which recognises three statistically distinct periods based on the pattern of velocity changes (R 2 = 0.85 (Fig. 3A)). Between 1992 and 2003, we observe a period exhibiting no significant trend (−0.03 ± 0.23 m yr −2 , p = 0.48). A substantial slowdown of −1.64 ± 0.37 m yr −2 (p < 0.05), starting around 2003 (but which could range from www.nature.com/scientificreports www.nature.com/scientificreports/ 2002-2004 due to a gap in our data), occurs until 2012, after which the most recent period of ice motion exhibits a significant increasing trend (0.58 ± 0.37 m yr −2 , p < 0.05).
We use the regional climate model MAR v3. 10 7,9 to calculate the mean surface melt production anomaly during the three distinct periods identified in our ice velocity anomaly time series; 1992-2003, 2003-2012, and 2012-2019 (Fig. 3B). Mean surface melt production anomalies across our study region rose by ~30% from −0.18 m yr −1 to 0.45 m yr −1 (w.e.) (~2.15 m yr −1 to ~2.77 m yr −1 (w.e.)) between 1992-2003 and 2003-2012. Following this period of sustained higher surface melt production, a ~15% decrease to 0.06 m yr −1 (~2.39 m yr −1 ) is observed for 2012-2019, when compared to the 2003-2012 mean. The period of ice velocity decrease therefore coincides with a period of enhanced ice surface melt, while rates of constant and accelerating ice velocity occur during periods of lower surface melt.
Previous work has argued that there is a significant relationship between antecedent surface melt production and ice velocity 40 . More recent work however suggests that this may be a statistical construct 43 , whereby as more antecedent melt years are included in the average value, R 2 tends further towards 1 ( Supplementary Fig. 12) as the data becomes increasingly smoothed and the individual points become less independent. With only one period of sustained slowdown, we are unable to test whether ice velocity slowdown is in response to a gradual increase in surface melt or to the passing of some surface melt threshold. Moreover, the response times of ice velocity to increases and decreases in surface melt forcing appear to differ -surface melt displays a long-term increasing trend from the early-mid 1990s before ice velocities begin to decrease in ~2003, whereas ice velocities stabilise and begin accelerating almost instantly in response to the large reduction in surface melt forcing from 2013 onwards. To investigate the impact of year-to-year variability in surface melt production on year-to-year velocity, we calculate a linear regression through detrended velocity and melt production anomaly time series (Fig. S13), which gives an R 2 of 0.08 (p = 0.11), indicating that there is no significant relationship between annual ice velocity and annual runoff, consistent with earlier work 41,43,48 .

Discussion
Changes in ice thickness and surface gradient can contribute to changes in ice motion through their associated impacts on driving stress. Thinning at the margins of the Greenland Ice Sheet has been observed since the early 1990s 2 , and from 1992-1998, the South West Greenland land-terminating sector thinned by ~0.02-0.4 m yr −1 51,52 . A period of stability characterised the late 1990s and early to mid-2000s 52 , followed by rapid thinning (~  40 . Ice surface elevation contours (grey lines) are from Howat et al. 88 , with the 1000 m contour bold to make clear the area across which the velocity time series was calculated (Fig. 3). The two tidewater glaciers to the north of the study region are masked out as they are undergoing different dynamic processes to the rest of the region. The base image is a MODIS (Terra) corrected reflectance image from EOSDIS NASA Worldview (https://worldview. earthdata.nasa.gov/). 1-1.5 m yr −1 ). Whilst some datasets assume a constant thinning of a similar magnitude through to 2014 53 , recent work points to a maximum thinning rate in 2012, after which thinning has continued at a more moderate rate 54 , albeit with differences in thinning trends between drainage basins 52,55 . Prior modelling work in our study region 40 indicates that changes in driving stress due to thinning could only explain 17-33% of the slowdown signal beyond 10 km from the ice-margin, and none of the signal beyond 50 km. This modelling work applied 20 m of thinning over 30 years, equating to a thinning rate of −0.66 m yr −1 . Outside of the period 2007-2011, there is no evidence for thinning rates greater than this during our study period [51][52][53]56 , and we thus do not believe that thinning rates could be responsible for more slowdown than already reported in Tedstone et al. 40 . Moreover, continued thinning from 2012 would be expected to continue to reduce ice motion, which contrasts with the acceleration in ice velocity observed during this period. As such, these findings imply that most of the ice velocity signal is controlled by processes operating at the ice-bed interface. The subglacial hydrological system exerts a critical control on ice motion at land-terminating margins of the Greenland Ice Sheet 38,57 . Our results provide support for a previously postulated mechanism 40 , whereby under a sustained multi-annual increase in surface meltwater production, the subglacial drainage system is characterised by both a gradual increase in the extent of channelisation, as well as the time during which these channels remain open and are thus able to evacuate water from the surrounding distributed drainage system 58 . These conditions promote the drainage of waters from weakly-connected 39 regions of the distributed drainage system, thereby reducing basal water pressure and associated basal lubrication and hence ice velocity over longer (decadal) timescales as effective pressures increase. Whilst it is hypothesised that increasing runoff to the bed will increase sediment pore water pressure, resulting in reduced sediment shear strength, increased sediment deformation and thus enhanced ice flow 59 , such behaviour contrasts with the observations reported here and elsewhere [40][41][42]44 . Consequently, our results suggest either that extensive layers of subglacial till are not ubiquitous across our study region or if they are present, the tills have not deformed pervasively in response to the sustained period of enhanced meltwater input.
Surface melt production in our study region peaks in 2012, decreasing thereafter (although surface mass balance remains negative), consistent with trends observed across the ice-sheet 46 . Concurrent with decreasing meltwater runoff, we observe an acceleration in ice velocity. During the period of sustained high melt, the drainage of basal waters from a larger area of the distributed drainage system causes a gradual depressurisation of the background water pressure with an associated increase in effective pressure. Reduced surface meltwater production would result in seasonally less extensive and smaller subglacial channels and consequently these channels will undergo faster creep closure 32 . This results in a reduced time-period over which the main subglacial channels   www.nature.com/scientificreports www.nature.com/scientificreports/ can evacuate water from the surrounding weakly connected, hydraulically distributed regions. As water pressure is no longer being systematically reduced, the background water pressure can increase through a gradual re-pressurisation of the subglacial hydrological system, as meltwater can recharge from basal melting without being evacuated by efficient subglacial channels.
Effective pressure (N) is calculated as N = Pi − Pw, where Pi refers to the ice overburden pressure and Pw to the basal water pressure 32,60 . Repressurising the subglacial hydrological system causes basal water pressure (Pw) to increase, thus reducing effective pressure and causing ice to accelerate. As effective pressure is a function of ice thickness, for the same increase in basal water pressure we would anticipate acceleration to be greater where ice is thinner. Although we expect channel closure to occur earlier and faster under thicker ice (due to higher ice overburden pressure) 26,49,50 , modelling studies suggest that once channel pressure falls to 90% of overburden pressure, there is less than 2 days difference in closure time for channels with a cross-sectional area of 10 m 2 compared to <1 m 2 50 . Moreover, it has been suggested that maximum channel closure rates are similar at both marginal and interior sites 49 . Consequently, we argue that any variation in repressurisation due to channel closure rates will only have a small impact upon effective pressures when compared to variability in ice thickness (which ranges from 0-1400 m across our study area). As a result, we hypothesise that the observed acceleration is greater for marginal regions of thinner ice due to the impact of basal meltwater recharge on water pressure in these areas offsetting a larger proportion of the ice overburden pressure (Fig. 2D).
The inland limit and spatial extent of efficient channel formation is subject to considerable debate 38,61 . Borehole and tracer studies and ice velocity records have been used to infer channels extending 40-80 km inl and 19,25,26,34,[62][63][64] , with high flow velocities of traced waters and the rapid transmission of pulses of meltwater from the ice sheet surface to the margin indicating efficient drainage 38,65 . These observations are supported by modelling studies, whereby efficient channels have been modelled up to 50 km inland, under 900 m thick ice 66,67 . The spatial extent of channelisation is influenced by the distribution of surface-to-bed connections, with high moulin density conducive to widespread and rapid channel development 68 . Further inland, ice thickness increases, surface slopes become shallower and runoff production decreases, leading to enhanced creep closure due to higher effective pressures and reductions in hydropotential gradients and wall-melting 32,61,69,70 . Under these conditions, modelling studies suggest that at a certain elevation and ice-thickness, efficient subglacial channels will not develop 49,50 in which case melt and ice-velocity are expected to scale positively 71 . Velocity measurements 140 km from the margin 71 , within the accumulation area of the Leverett catchment, indicate a 2.2% increase in annual velocity between 2009 and 2012. In our study, we are unable to determine whether a positive change in ice velocity has occurred at elevations above 1300 m (a.s.l.) due to increasing noise in our dataset further inland. However, GPS observations reveal reduced velocities to at least 80 km inland (at ~1500 m a.s.l) in the years of record surface melt in 2010 48 44 . At a regional scale, this radar data suggests a weaker slowdown trend than that reported either here, or by Tedstone et al. 40 44 . Methodological differences between the studies may affect the conclusions drawn as this study and that of Tedstone et al. 40 applied feature tracking to image pairs either side of an annual baseline (352-400 days in Tedstone et al. 40 , 352-384 days here), whereas Joughin et al. 44 composed mosaics of winter velocities, comprised of data sampled between September and May. Regardless, this raises the possibility that the previously observed slowdown 40 may be related to some local process, rather than regional-scale surface melt forcing 44 , for example water piracy as a result of dynamic thinning on Jakobshavn Isbrae has been suggested as a driver of slowdown in the north of the West Greenland land-terminating sector 44 . However, by extending the study area ~120 km to the south, and by increasing the number of common pixels used in our time series by an order of magnitude compared to that of Tedstone et al. 40 , we extend confidence in the observed slowdown of ice motion from the early-mid 2000s. Moreover, by calculating an ice velocity anomaly 72 , we remove the influence of any biases between sensors from the observed trends, further increasing confidence that this sector of the ice sheet underwent significant deceleration during this period.
A variety of models exist to simulate the subglacial drainage system and its impact on ice dynamics in response to increasing meltwater fluxes in a warming climate 73 . However, whilst current models can reproduce observed dynamics of Greenland's land-terminating margins over days/weeks 32,50 , on a seasonal timescale 39,67,74 , and across several years 66,75 , there is a need to consider the longer, decadal and multi-decadal response of ice motion to surface meltwater forcing in order to better reproduce observed dynamics and project future change. On a decadal scale, recent modelling predicts increased ice-motion and thus ice flux into the ablation zone under enhanced summer melt within 44 km of the ice-margin 76 (their Figure 5C). These predictions contrast with the results presented both here and from other GPS and satellite-derived observations 37,[40][41][42][43][44] . Further work is therefore required to investigate what aspects of the model set-up or framework cause the modelled ice-motion to accelerate on a multi-annual timescale in response to increasing surface runoff, in contrast to numerous observations reported in this and other studies. Whilst we do observe reduced slowdown as ice thickness increases, our results and those of others show that ice motion within the ablation area has slowed as a result of increased surface runoff, and we do not observe any appreciable or spatially extensive increases in ice velocity as a result of an extended period of high surface melt. We therefore argue that in order to better represent ice dynamical processes and project future changes in coupled hydro-dynamic models, the models must better utilise the available data from field and satellite observations both in model set-up and for assessing model performance, and there is a need for changes in the numeric representation of processes, including ice-hydrology coupling.

Scientific RepoRtS |
(2020) 10:1647 | https://doi.org/10.1038/s41598-020-58355-2 www.nature.com/scientificreports www.nature.com/scientificreports/ Our findings demonstrate that ice at surface elevations below 1300 m (a.s.l.) in the South West land-terminating region of the GrIS underwent deceleration during the mid to late-2000s and early-2010s with subsequent acceleration following sustained reductions in surface melt production, with this acceleration proportionally greater in areas of thinner ice. We hypothesise that this recent acceleration is driven by increases in basal water pressure offsetting a larger proportion of ice overburden pressure, and thus causing a greater reduction in effective pressure under thinner ice. This behaviour has not been observed previously on a multi-annual timescale. Whilst we have demonstrated that changes affecting the subglacial hydrological system are the most likely driver of ice velocity change in our study region, we are unable to conclusively state whether antecedent melt production is the predominant factor controlling ice motion, or whether some threshold in the amount of melt is required in order to drive a change in velocity. However, it is clear from our results that the period of ice deceleration occurred during a period of sustained high melt production, and ice motion accelerated subsequently as melt production decreased to consistently lower values from 2013 onwards. Moreover, we observed no indication of an ongoing slowdown once surface melt production exhibited a sustained decrease. Thus, the behaviour observed supports a process-based understanding of the links between hydrology and ice dynamics in this land-terminating sector of the Greenland ice sheet. The results displayed in this study show that the observed trends in ice velocity are not statistical artefacts or the result of biases between sensors 72 . Further, we see no evidence of speed-up at the elevations studied under a warming climate, irrespective of the bed conditions. Future work is required to improve observations of change at higher elevations, and at other land-terminating sectors of the ice sheet, in order to assess the extent to which these processes operate on an ice sheet scale.
Furthermore, while the results presented here relate specifically to land-terminating margins, it has been argued that the dynamics of tidewater glaciers are a product of both oceanic and atmospheric forcing 77,78 , with the latter potentially driving a positive feedback as a result of its influence upon fjord water circulation and thus submarine melt 77,79,80 . This positive feedback contrasts with our findings at land-terminating margins, indicating that at an ice-sheet scale, the relationship between hydrology and ice dynamics is complex and spatially variable. Since the extent to which surface melt processes impact on tidewater glacier dynamics remains uncertain, analysis of existing large-scale, multi-decadal datasets 5 at marine-terminating margins should be undertaken with respect to surface melt change, in order to better project the future response of the Greenland ice sheet to changing atmospheric and oceanic conditions.

Methods
Remote sensing of ice velocity. We apply optical feature tracking 40 Figures S1-S3). A pair length of 368 ± 16 days is used to minimise the impact of any seasonal variability in interannual ice velocity (Fig. S4). We use oriented correlation, matching the feature of gradient orientation for each pixel 82 , and use a combination of spectral bands at wavelengths ~0.52-0.69 μm (bands 2 + 3 for Landsat 5 and 7, 3 + 4 for Landsat 8). We enhance the images by applying a principal component analysis to these bands, and subsequently, a highpass gaussian filter is used to enhance surface features such as crevasses and reduce the impact of temporally stable features relating to the basal topography 40,81 . We use a reference window of 80 pixels and set the search window based upon prior velocity estimates 83 . Following Tedstone et al. 40 , we apply a median coregistration to the output velocity fields in order to remove errors associated with georeferencing.
The coregistered velocity fields are then fused via a spatio-temporal median over annual or multi-year periods. Low-quality velocity estimates are removed through filtering by a threshold signal-to-noise ratio value (snr threshold = 6), identified by calculating the value beyond which the median absolute deviation of velocities over stable ground becomes asymptotic 81 . The final velocity field for each time period is composed of the median value of all velocity fields within the time period at each pixel. To calculate uncertainty at each pixel, we compute a 1σ confidence interval for each component of the velocity field in the form: Where MAD is the median absolute deviation over stable ground, N is the number of velocity fields used in the computation of the median velocity, σ is the 1σ confidence interval, and k and α are parameters determined for each time period from the stable ground velocity which is known to equal zero. This relationship is extrapolated on-ice using the appropriate values of MAD and N at each pixel. Following Tedstone et al. 40 , we discard pixels with σ > 60 myr −1 from the subsequent analyses. In addition, we mask out the tidewater glaciers to the very north of the study region, limit our analysis to pixels below 1000 m above sea-level, and retain only the velocity fields with an area coverage above 9250 km 2 . The median velocity is subsequently calculated across the 71703 pixels common to all the remaining merged velocity fields (Fig. S5) in order to avoid spatial bias influencing the change signal.
Anomaly-based time series. Recent work indicates that velocity magnitude has a biased mean, with this bias increasing with the standard deviation of the velocity components (and so with noise), causing an artificial negative velocity trend 72 . To mitigate this effect, we follow the velocity anomaly approach of Dehecq et al. 72 . First, we calculate a mean of all velocity pairs covering the period 1992-2019 described as V 0 . The velocity anomaly is defined as the value of the difference vector V t − V 0 projected on the mean velocity vector: www.nature.com/scientificreports www.nature.com/scientificreports/ The result of this approach is to centre the noise distribution symmetrically around zero (Fig. S6) such that there is no bias in the mean value, removing any artificial slowdown trends due to variability in noise magnitude between sensors 72 . The resultant velocity anomaly fields are displayed in Fig. S7.
To assess the long-term trends in our time series, we first compute a simple linear regression through the data, and then test whether the data can be divided into three statistically different segments (Figures S8-S10) using the non-parametric Mann-Whitney Wilcoxon Test. We select the pair of breakpoints with the lowest RMSE, and that are statistically significant with 99% certainty, which gives breaks at 2003 and 2012. Spatial trends in ice velocity. We  where c 1 is the first period (i.e. 1992-2003) and c 2 is the second period (i.e. 2003-2012). We remove pixels with uncertainty greater than 60 m yr −1 . We also filter by Velocity Vector Coherence (VVC) 81 , which follows the form: where T is the set of N velocity estimates V(i,j,t) merged to obtain the median velocity V (i,j) at pixel(i,j). We filter out all pixels with VVC < 0.45. In addition, we erode the edge of our ice mask by 3 pixels in order to limit the influence of noise at the ice margin.
Surface mass balance. We obtain surface mass balance (SMB) data from the MAR v3.10 regional climate model 7 , forced by NCEP-NCARv1 from 1992 to 2019. We limit our analyses of SMB below an ice surface elevation upper limit of 1600 m (a.s.l.). Little lake drainage occurs above this elevation [84][85][86] and it has been argued that this is an approximate maximum elevation where crevasses, and thus moulins, are likely to form 87 . Whilst surface meltwater can runoff from elevations above this, surficial drainage is less likely to occur in high elevation regions due to the shallower surface slope 86 . Consequently, below this threshold elevation, we can be confident that surface meltwater drains to the ice-bed interface, and so influences ice motion. Regardless, melt at all elevations has increased from the 1958-1987 average in the period 1988-2013 87 , with the difference between the means of the period being positive at p < 0.05 from 400-2600 m (a.s.l.).

Data availability
The Landsat imagery used in this study was provided by the United States Geological Survey and the European Space Agency third party missions program and are freely available.