Supraglacial River Forcing of Subglacial Water Storage and Diurnal Ice Sheet Motion

Surface melting impacts ice sheet sliding by supplying water to the bed, but subglacial processes driving ice accelerations are complex. We examine linkages between surface runoff, transient subglacial water storage, and short-term ice motion from 168 consecutive hourly measurements of meltwater discharge (moulin input) and GPS-derived ice surface motion for Rio Behar, a moulin-terminating supraglacial river catchment on the southwest Greenland Ice Sheet. Short-term accelerations in ice speed correlate strongly with lag-corrected measures of supraglacial river discharge ( r = 0.9, τ = 0.7, p < 0.01). Though our 7 days record cannot address seasonal-scale forcing, diurnal ice accelerations align with normalized differenced supraglacial and proglacial discharge, a proxy for subglacial storage change, better than GPS-derived ice surface uplift. These basal in absolute drive accelerations. present of in situ in a large supraglacial river draining the ice sheet surface, just of it plummets into a GPS of ice surface motion record brief accelerations in ice sliding speed that follow daily cycles in meltwater entering the moulin. By comparing these measurements with proglacial river discharges leaving the ice sheet, we identify daily fluctuations in subglacial water storage that track short-term accelerations in ice motion. These findings affirm the importance of supraglacial rivers to subglacial water pressure and ice dynamics, even in relatively thick ice >40 km inland from the ice terminus.

Evidence for transient cavity evolution is drawn primarily from GPS-derived correlations of horizontal ice speed with vertical ice surface uplift (interpreted as a proxy for total subglacial water storage, S) or its first derivative (interpreted as subglacial water storage rate-of-change, ΔS). GrIS horizontal ice sliding speed broadly covaries with vertical surface uplift over seasonal time scales (e.g., Bartholomew et al., 2012Bartholomew et al., , 2010Hoffman et al., 2011), but variations at shorter timescales tend to correlate better with its derivative (Andrews et al., 2018;Cowton et al., 2016;Hoffman et al., 2011). Such correlations are typically weak and spatially variable due to a range of factors confounding estimation of basal uplift from ice surface elevation measurements (Andrews et al., 2018;Hoffman et al., 2011). Therefore, it is difficult to infer interactions between surface melting, subglacial water storage, cavity growth, and ice motion for the GrIS, despite previous success on mountain glaciers (e.g., Armstrong & Anderson, 2020;Bartholomaus et al., 2011Bartholomaus et al., , 2008 To study the links among supraglacial runoff, subglacial water storage fluctuations, and short-term ice motion, we present in situ measurements of moulin input (i.e., supraglacial discharge), ice surface speed, and ice surface uplift for Rio Behar, a large supraglacial river on the GrIS midelevation (>1,200 m a.s.l.) ablation zone ( Figure 1). We compare daily cycles in these variables with PROMICE automated weather station (AWS) measurements of surface energy balance and ablation (Fausto & van As, 2019), and with proglacial river discharges from three gauging stations downstream Tedstone et al., 2017;van As et al., 2019). We present GPS measurements of horizontal ice surface speed and vertical uplift, SMITH ET AL. 10.1029/2020GL091418 2 of 11 Figure 1. Study area in southwest Greenland. Black star shows location of our GPS measurements of ice surface motion and Acoustic Doppler Current Profiler (ADCP) measurements of moulin input (supraglacial discharge) in Rio Behar, a large supraglacial river penetrating the ice sheet >40 km from the ice edge. Field work was conducted ∼750 m upstream of the Rio Behar terminal moulin. Black outline delineates the surface catchment (60.02 km 2 in July 2016). White stars locate proglacial river gauging stations; black square locates PROMICE KAN_M automated weather station. Background is a July 26, 2016 true-color Landsat-8 satellite image. and use them to estimate subglacial storage S and rate-of-change ΔS, respectively. We also compute alternate proxies for S and ΔS by differencing normalized supraglacial and proglacial discharge hydrographs (adapted from Bartholomaus et al., 2008Bartholomaus et al., , 2011McGrath et al., 2011;Armstrong & Anderson, 2020). We conclude that diurnal cycles in supraglacial river discharge drive ice accelerations through ΔS, confirming that transient water storage and cavity growth are important influences on GrIS subglacial basal pressure and short-term ice motion.

Observational Data
In July 2016, the Rio Behar terminal moulin was located at 67.047°N, −49.033°W, with an upstream drainage catchment of ∼60.2 km 2 and mean surface elevation >1,200 m ( Figure 1). We established a field camp to measure moulin meltwater input and ice surface motion ∼750 m upstream (∼67.050°N, −49.018°W). During July 5-13, 2016, we collected 174 measurements of supraglacial river discharge using a SonTek Riv-erSurveyor M9 Acoustic Doppler Current Profiler (ADCP) and methods of Smith et al. (2017). A Tyrolean cableway was suspended over the river to safely and repeatedly tow the ADCP back and forth across the channel every hour (Figures S1-S4). In total 847 ADCP transects were acquired, of which 677 later passed rigorous quality-assurance screening and were used to compute 168 consecutive hourly discharge measurements from July 6-13 (Text S1-S3, Figure S5, Tables S1-S2, and Data sets S1-S2).

Proxies for S and ΔS
GPS-derived vertical positions and their first derivative were used to estimate subglacial storage S and rateof-change ΔS (e.g., Bartholomew et al., 2012;Cowton et al., 2016;Text S7). Proxies for S and ΔS were also computed by adapting a meltwater input-output approach (Armstrong & Anderson, 2020;Bartholomaus et al., 2008Bartholomaus et al., , 2011McGrath et al., 2011) comparing relative timings of supraglacial and proglacial river discharge hydrographs (Text S7). Hydrographs were normalized and differenced (supraglacial minus proglacial) to assess their relative timings and shapes at Rio Behar moulin versus at the ice edge. These "discharge-difference" proxies are unitless and do not satisfy mass conservation. They characterize instantaneous net water storage changes, not subglacial routing delays and/or storages known to retard proglacial discharges longer than 24 h (e.g., Chandler et al., 2013;Chu et al., 2016;Pitcher et al., 2020;Rennermalm, Smith, et al., 2013;Smith et al., 2015;van As et al., 2017). From dye tracing experiments, subglacial routing from ∼1,300 m elevation takes 1-3 days (Chandler et al., 2013, site L57), or ∼2-5 days from proglacial hydrograph analysis (van As et al., 2017). Such subglacial delays and storages are irrelevant to our purpose here, which is simply to characterize instantaneous subglacial conditions at our field site, not Lagrangian transport to the ice edge. Descriptions of all data, methods, and uncertainties are presented in SI (Texts S1-S9, Figures S7-S14, Tables S1-S5).

Correlations of Ice Speed with Other Variables
We find strong diurnal cycles in all variables except surface elevation, with daily accelerations in horizontal ice speed closely tracking moulin input and melt energy ( Figure 2 and Table S3). A consistent progression is observed in the timing of daily peaks, with melt energy and air temperatures peaking near local solar noon, followed by sequential peaks in ice ablation, proglacial discharge, moulin input, and horizontal ice speed ( Figure 3). The timing of daily peaks is most consistent for melt energy, proglacial discharge, moulin input, and ice speed, whereas peaks in air temperature, ablation, and especially ice surface elevation are more variable as indicated by their peak timing range (Table S3).
After lagging our GPS-derived horizontal ice speed time-series to correct for its mean timing offsets with the other variables, we compute correlations between potential forcing variables and ice speed using Pearson's r, which assumes a linear relationship, and Mann-Kendall's  , which does not assume a linear relationship between variables. We find that ice speed correlates strongly with moulin input (r = 0.90,  = 0.70), melt energy (r = 0.90,  = 0.67), and proglacial discharge (r = 0.88,  = 0.71) (Figure 4 and Table S4). We find moderately strong correlations with ice surface ablation (r = 0.74,  = 0.59) and air temperature (r = 0.70,  = 0.54), which are drivers of runoff and melt energy, respectively. Lowest correlation is found for SMITH ET AL.
10.1029/2020GL091418 4 of 11 Figure 2. In situ measurements of (a) melt energy, air temperature, and ice ablation from KAN_M; (b) Rio Behar moulin input (supraglacial discharge) and proglacial discharge measured at Kangerlussuaq (-5h timing correction applied); (c) ice surface speed and elevation. Colored envelopes (b), (c) represent measurement uncertainties. detrended ice surface elevation (i.e., uplift, r = 0.36,  = 0.24, Table S4 and Figure 4e). All correlations are statistically significant (p < 0.01). Unlike melt energy (which turns negative, implying nocturnal refreezing), moulin input persists throughout the night. Because (i) moulin input closely tracks (and derives from) melt energy; (ii) virtually all surface runoff generated within Rio Behar catchment flows to its terminal moulin; and (iii) an observed 6h time lag between peak melt energy and peak supraglacial discharge (Table S3) is similar to a previously calculated catchment routing delay for Rio Behar (i.e., estimated time-to-peak t p = 5.5 h, Smith et al., 2017) we infer that supraglacial river discharge, a product of catchment-integrated melt energy, is the dominant supraglacial forcing variable driving our locally recorded ice speed variations.

Comparison of Short-Term Ice Accelerations with S and ΔS
To further investigate drivers of short-term ice speed variations, we test proxies of subglacial water storage S and rate-of-change ΔS calculated from GPS-derived ice surface observations (following Anderson et al., 2004;Andrews et al., 2018;Cowton et al., 2016;Harper et al., 2007;Hoffman et al., 2011;Howat et al., 2008) and by differencing normalized hydrographs of supraglacial and proglacial river discharge (Texts S7-S9). Implicit in the latter "discharge-difference" calculations are assumptions that englacial storage is negligible; that en/subglacial melting is negligible; that subglacial routing delays are irrelevant to instantaneous net storage; and that proglacial discharge reflects overall regional basal water pressure, allowing Rio Behar moulin input to be compared with regional proglacial discharge despite its smaller spatial domain (60 km 2 versus ∼2,800 km 2 to 1,750 m a.s.l.) and absolute discharge magnitude (∼6-38 m s − 1 versus ∼800-1,300 m 3 s −1 ).
Comparison of our observed horizontal ice speeds with both proxies for S and ΔS suggests that ΔS drives short-term accelerations in ice speed (Figures 5, S11, S13, S14 and Table S5). This conclusion is clearest from the discharge-difference proxies, with ΔS tracking ice as well or better than S (see Figure 5c versus 5b; S13c versus S13b; Figures S11d, S14d vs. S11c, S14c). This same conclusion may be drawn, albeit less compellingly, from conventional GPS-derived S and ΔS proxies (Figure 5a versus Figure 2c; Figures S11b, S14b versus S11a, S14a). For both methods, ΔS generally correlates with ice accelerations better than S (Table S5) suggesting that changes in subglacial water storage force short-term ice speed accelerations at our field site.

Discussion and Conclusion
We find that diurnal cycles in moulin input are the primary driver of short-term accelerations in ice sliding velocity at our field site (Figure 4). This finding supports previous work  and the conclusion that over short time scales, diurnal variability in supraglacial river input imposes a first-order control on subglacial water pressure fluctuations. Furthermore, while short-term accelerations in ice speed closely follow moulin input, they also tend to track proxies of subglacial water storage change (ΔS) better than proxies of absolute storage (S) (Figures 5, S11, S13, and S14), suggesting that nightly peaks in subglacial water storage drive subglacial basal pressure and short-term ice motion.
This conclusion is more evident in discharge-difference proxies than conventional GPS-derived proxies (e.g., Figures 5c versus 5a; S13c versus S13a; S11d, S14d versus S11b, S14b). This is likely due to our inability to assess the impact of changes in the ice column due to variations in vertical strain, making our GPS-derived proxies susceptible to local and nonlocal forcings (e.g., Price et al., 2008;Ryser et al., 2014;see Text S7). Differencing supraglacial and proglacial hydrographs, therefore, may characterize subglacial water storage conditions more sensitively than small vertical ice surface elevation changes, which are inherently difficult to detect and have multiple sources of uncertainty (  . Mean daily peak timing (circles) and timing range (gray bars) of observed variables. Diurnal cycles in melt energy and air temperatures peak around solar noon, followed by peaks in ice ablation, proglacial discharge, moulin input, ice speed, and ice surface uplift. GPS-derived daily peaks in uplift are temporally variable (Table S3). Peak proglacial discharge has little timing variability and is shifted −5 h earlier to account for the mean timing offset between Kangerlussuaq and the ice edge. This figure presents only timing of daily peaks, not subglacial routing delays and/or storages of meltwater.
A discharge input-output approach (e.g., Armstrong & Anderson, 2020;Bartholomaus et al., 2011Bartholomaus et al., , 2008McGrath et al., 2011), comparing moulin inputs with proglacial outputs, offers an alternate strategy for characterizing subglacial water storage and its link to basal sliding laws. Future studies, for example, could develop discharge-difference proxies over longer time scales and larger study areas by pairing surface-routed climate model output (e.g., Smith et al., 2017;Yang et al., 2019) with proglacial discharge records van As et al., 2019), to relate net increases/decreases in ΔS to regional ice speed variations.
It is well-known that evolution of the subglacial system from inefficient to efficient states acts to modulate the ice dynamical response to supraglacial water inputs (e.g., Bartholomew et al., 2010;Hoffman et al., 2011). We find that peak or ascendant ΔS is associated with localized GrIS velocity accelerations (Figures 5c and S13c). This suggests that highest subglacial water pressure (and ice sliding speed) occurs when subglacial cavities are growing the fastest, not when their volume is largest (e.g., Cowton et al., 2016;Iken et al., 1983). As such, steady state theoretical basal sliding laws-which assume a relationship between cavity size and subglacial pressure -do not accurately represent transient behavior of the subglacial system.
It is important to note that the strong correlation between moulin input and ice velocity reported here (Figure 4d) is unlikely to hold over an entire melt season. Previous work has clearly established that Greenland ice sliding velocities are strongly influenced by long-term seasonal evolution of the subglacial hydrological system (e.g., Andrews et al., 2018;Bartholomew et al., 2010;Hoffman et al., 2011;Nienow et al., 2017). Our short 7 days record captures neither the early nor late melt season, when subglacial efficiency (and associated ice speeds) undergo extensive changes. Subglacial evolution makes melt-driven proxies inappropriate for estimating ice motion over the entire melt season Bartholomew et al., 2010) or multiple years (Davison et al., 2019;Tedstone et al., 2015). Over short time scales, however, we find that diurnal cycles in moulin input are the primary driver of fluctuating subglacial water pressures and associated ice accelerations-even in relatively thick ice (∼1 km) more than 40 km inland from the ice edge. Some slight SMITH ET AL.   Table S4. Moulin input displays strongest correlation with ice speed (r = 0.90, τ = 0.70, p < 0.01). differences in peak timings between ΔS and ice motion, as well as some nonlinear behavior on descending limbs (Figures 5c and 13c) are discussed further in Text S9.
This study adds to a small but growing collection of GrIS supraglacial streamflow measurements (Carver et al., 1994;Chandler et al., 2013;Echelmeyer & Harrison, 1990;Gleason et al., 2016;Holmes, 1955;McGrath et al., 2011;Smith et al., 2015Smith et al., , 2017. With peak daily discharges of 26.59-37.61 m 3 /s (Table S2), the discharges reported here are far larger than those collected in most supraglacial streams, but are typical for trunk supraglacial rivers in southwest Greenland . Nearly all of them terminate in moulins Yang & Smith, 2016), and the high diurnal variability we observe (ranging from 19.05 to 30.50 m 3 /s, Table S2) signifies that their subglacial channels are likely out of equilibrium with supraglacial inputs for large portions of the day, such that associated accelerations in ice speed are driven by addition or removal of water outside of the channelized system.
Based on satellite mapping (e.g., Lampkin & VanderBerg, 2014;Smith et al., 2015;Yang et al., 2015Yang et al., , 2016Yang & Smith, 2013 and topographic modeling (e.g., Banwell et al., 2016Banwell et al., , 2012Crozier et al., 2018;Karlstrom & Yang, 2016;King et al., 2016), we maintain that supraglacial rivers likely drive ice accelerations near hundreds of other terminal moulins as well. Process-level understanding and modeling of subglacial hydrology and associated ice dynamics should presume large, strongly diurnal inputs of meltwater entering hundreds of supraglacial river moulins distributed throughout Greenland's ablation zone. These inputs, countered by water output discharged beneath outlet glaciers, trigger short-term fluctuations in subglacial water storage that drive short-term accelerations in ice sheet motion.
SMITH ET AL.
10.1029/2020GL091418 7 of 11 , and discharge-difference (c) tracks ice speed better than GPS (a). Figure S13 presents another version of this figure using AK4 proglacial discharges.

Acknowledgments
We dedicate this paper to the memory of Dr. Konrad Steffen