Challenges in predicting Greenland supraglacial lake drainages at the regional scale

A leading hypothesis for the mechanism of fast supraglacial lake drainages is that transient extensional stresses briefly allow crevassing in otherwise-compressive ice flow regimes. Lake water can then hydrofracture the crevasse to the base of the ice sheet, and river inputs can maintain this connection as a moulin. If future ice-sheet models are to accurately represent moulins, we must understand their formation processes, timescales, and locations. Here, we use remote-sensing velocity products to constrain the relationship between strain rates and lake drainages across ∼1600 km in Pâkitsoq, western 5 Greenland, between 2002–2019. We find significantly more-extensional background strain rates at moulins associated with fastdraining lakes than at slow-draining or non-draining lake moulins. We test whether moulins in more-extensional background settings drain their lakes earlier, but find insignificant correlation. To investigate the frequency that strain-rate transients are associated with fast lake drainage, we examined Landsat-derived strain rates over 16and 32-day periods at moulins associated with 240 fast lake drainage events over 18 years. A low signal-to-noise ratio, the presence of water, and the multi-week 10 repeat cycle obscured any resolution of the hypothesized transient strain rates. Our results support the hypothesis that transient strain rates drive fast lake drainages. However, the current generation of ice-sheet velocity products, even when stacked across hundreds of fast lake drainages, cannot resolve these transients. Thus, observational progress in understanding lake drainage initiation will rely on field-based tools such as GPS networks and photogrammetry.

average uncertainties in our study region are ∼140 m/yr, or ∼70% of the average velocity. For additional quality control, we discarded observations at all pixels whose flow directions deviated more than 90 • from the wintertime flow direction, which we defined from the twenty-year average (Joughin et al., 2016a, b), which amounted to 5% of all observations. Overall, we used a subset of the TUD dataset that contained 181 velocity maps associated with 97 fast-draining lake events at 18 low-elevation 125 lakes over the 14-year period 2002-2015. To study strain rates in melt seasons 2013-2019, we used the GoLIVE product (Scambos et al., 2016;Fahnestock et al., 2016) (Fig. 1). We selected velocity observations constructed from Landsat scenes spaced by 16 days. This differs from our 16-32 day choice for the TUD velocities for two reasons: (1) There were more GoLIVE observations available overall, allowing tighter constraints in image spacing; and (2) including 32-day observations significantly smoothed any strain rate signal. This 130 dataset covers our entire study area. The GoLIVE pixel size ∆x = ∆y is 300 m (20×20 Landsat OLI Band 8 pixels) and average uncertainties in our study region are ∼55 m/yr, or ∼30% of the average velocity. As with the TUD velocities, we discarded observations with flow directions >90 • from background (Joughin et al., 2016a, b), which here amounted to 3% of all observations. Thus, we used 97 velocity maps associated with 152 fast-draining lake events at 59 distinct lakes over the 7-year period 2013-2019. 135

Calculation of principal strain rates
We analyzed each velocity dataset described above separately. For each dataset, we smoothed the velocities with a 1 km × 1 km boxcar filter, which carries through many high-strain-rate features (King, 2018). We propagated errors in the velocity observations through this filtering using the following formula: (1) 140 for an n-element filter (in the case of our 3×3 filter, n = 9), for both components of the surface velocity v (eastward and northward). We then calculated the two horizontal principal strain rates˙ 1 and˙ 3 by taking centered spatial derivatives of these smoothed velocities, as follows: These definitions follow Harper and Humphrey (Harper et al., 1998). We used twice the native pixel size, 2∆x and 2∆y, as the differential spatial length, effectively averaging strain rates across three adjacent pixels in each direction. We assumed plane

Moulin -lake distance
We located the moulins that drained the supraglacial lakes analyzed by Morriss et al. (2013Morriss et al. ( ) each year, using 2009Morriss et al. ( -2019 WorldView images. Due to inconsistent coverage, we were not able to identify the moulin location for every lake in every year; for 78 lakes over 7 years, we identified 262 moulins associated with 66 lake locations. The moulin history for Lake #31 is 275 shown in Fig. 5. We chose this lake as an example due to the exceptionally good WorldView coverage: 7 out of 11 of the years have cloud-free scenes where we were able to visually locate the lake-draining moulin. We see multi-year reuse of the Lake #31 moulin as it advects downstream; for example, the lake drained into this moulin in 2009, 2010, and 2012, when it was respectively located 340, 350, and 520 m from the lake center ( Fig. 5b- (2015) and 570 m (2016) from the lake center ( Fig. 5e-f). In 2017, the lake drained rapidly through the same drainage channel and likely to the same moulin, but snow drifts obscured the view of the full path of the water (Fig. 5g).
We found that 210 (80%) of the 262 lake-draining moulins were located within 1 km of the centers of their respective lakes.

285
The greatest observed distance between a lake and its moulin was 3.8 km, where in 2016 a slow-draining lake at ∼1010 m elevation (Lake #36) drained into a moulin at ∼980 m elevation. That moulin was associated with the rapid drainage of a closer lake (Lake #42), 400 m from the moulin, as well as a second slow-draining lake (Lake #41) located 2.1 km upstream of the moulin. Slow-draining lakes are typically associated with water overtopping the basin and flowing downstream into a moulin (Tedesco et al., 2013;Kingslake et al., 2015). Thus, the timing of slow lake drainage reveals more about local topography and 290 runoff rates than local stress or strain rate changes, and for that reason, we did not analyze the timing of slow-draining lakes.
Our 2009-2019 dataset of 262 lake-draining moulins contained moulins associated with 95 fast-draining lake events. The mean and median distance between fast-draining lake centers and their moulins were 463 m and 298 m, respectively (Fig. 5i).
The maximum observed distance was 2.9 km, in 2010 at an elevation of ∼1200 m (Lake #18).
3.3 Lake drainage type 295 We analyzed our dataset of 78 lakes over the 11 summers (2006-2010 and 2013-2018) that followed winters with Greenland Ice Sheet velocity mosaics (Fig. 1). We sorted the lakes into four drainage types each year: Fast, slow, non-draining, and unknown. We classified all lakes as one of these four types in 2010-2019, but for 2006-2009, we use data from Morriss et al.
(2013), who only identified fast-draining lakes. We continued with analysis for only lakes whose type we identified with high confidence and all fast-draining lakes identified by Morriss et al. (2013), for a total of 319 lakes. This is a subset of our full  3.3.1 Effect of elevation on lake drainage type Figure 6 shows the distributions of high-confidence fast-, slow-, and non-draining lakes across elevation and strain rate regime.
Because it shows only lakes whose drainage types we identified with high confidence, Fig. 6 does not accurately represent population sizes (Sect. 3.3); however, it does provide a means for comparison of elevation and strain rates across lake-drainage type. We find that fast-draining lakes (N = 212) are located at lower elevations than slow-or non-draining lakes (p < 10 −11 ), 320 and non-draining lakes (N = 32) are located at higher elevations than fast-or slow-draining lakes (p < 10 −12 ), as shown in

Effect of wintertime strain rate on lake drainage type
The lake drainage type also varies with winter principal strain rates at the lake-draining moulin. The variation of the lessextensive winter principal strain rate,˙ 1 , with lake drainage type is shown in Fig. 6c-d. The˙ 1 at fast-draining and at non-325 draining lakes was higher (more extensional) than at lake-draining moulins (p = 0.04 and p = 0.007, respectively). Though significant, we do not further analyze this difference, as lake drainage fundamentally requires extension. Thus, we concentrate on the more-extensional winter principal strain rate,˙ 3 . The variation of˙ 3 with lake drainage type is shown in Fig. 5e-f. We find that moulins associated with fast-draining lakes (N = 212) are sited in more extensional areas compared to slow or nondraining lakes (p = 0.01), but with a substantially weaker relationship than that with moulin elevation, as seen in the staircase 330 plots ( Fig. 6b and f). Slow-draining and non-draining lakes did not have significantly different background˙ 3 than fast-draining lakes.

Lake drainage completeness
We classified each high-confidence lake drainage event as "complete" or "partial", following Chudley et al. (2019). Over the 11 summers (2006-2010 and 2013-2018) for which we have wintertime strain rate data, we found 73 complete fast-draining 335 events at 35 unique lakes, 16 partial fast-draining events at 13 unique lakes, 31 complete slow-draining events at 25 unique lakes, and 42 partial slow-draining events at 27 unique lakes. This subset of our drainage data spans 60 unique lakes. Figure 7 shows the distributions of complete and partial high-confidence fast-and slow-draining lakes across elevation and strain rate regime. We find that complete fast-draining events (N = 73) drain into moulins at significantly lower elevations 340 than the other types (p < 10 −4 ), and that partial slow-draining events (N = 42) drain into moulins at significantly higher elevations (p < 10 −9 ), as shown in Fig. 7a-b. Complete slow-drainage events and partial fast-drainage events did not occur at significantly different moulin elevations.

Effect of wintertime strain rate on lake drainage completeness
We found no significant variation of more-compressive wintertime principal strain rate˙ 1 across lake drainage types (p > 0.1,

Lake drainage dates
We examined our dataset of high-confidence fast-draining lakes (N = 212) to find relationships between the date of fast-lake drainage and wintertime strain rates at the lake-draining moulin, moulin elevation, or the date of melt onset that year.
3.5.1 Effect of wintertime strain rate on lake drainage date  Figure 8 shows the comparison of local wintertime strain rates to lake drainage dates the following summer. We tested for correlations between the drainage date of high-confidence, fast-draining lakes and the principal strain rates at the site of the moulin connected to each lake in each year (N = 212). We compared the drainage date to both˙ 1 (more-compressive principal strain rate) and˙ 3 (moreextensive principal strain rate) at each lake-draining moulin. Within each melt season, there is a weak trend for lakes connected 360 to moulins sited in more-extensional regimes to drain earlier. Both principal strain rates,˙ 1 and˙ 3 , correlate inversely with drainage date, although neither is statistically significant in any of the 11 melt seasons we studied.

Effect of elevation and runoff on lake drainage date
We also compared the drainage dates of high-confidence, fast-draining lakes (N = 212) to the surface elevation of the moulin and the date of melt onset for each year. We defined the date of first melt as the first day when the 2-meter air temperature at 365 the Swiss Camp meteorological station exceeded zero, averaged over the four available sensors there (Steffen et al., 1996). We simultaneously regressed moulin elevation, melt onset date, winter˙ 1 , and winter˙ 3 onto the dates of high-confidence fast lake drainage over the 11-year study interval. We found that only moulin elevation (p < 10 −9 ) and date of melt onset (p < 10 −8 ) correlated significantly. Winter strain rates˙ 1 (p = 0.08) and˙ 3 (p = 0.62) showed insignificant correlation with lake-drainage date. These higher p-values indicate that much of the univariate correlation of˙ 1 and˙ 3 with drainage date, shown in Fig. 9 370 and presented in Sect. 3.5.1, is due to relationships between˙ 1 and˙ 3 and moulin elevation, most likely due to the inverse correlation of the magnitude of principal strain rates with elevation on the ice sheet (Poinar et al., 2015). Overall, our findings suggest that wintertime strain rate is not a useful predictor of lake drainage date the following summer.

Evolution of strain rates at GPS station-pair midpoints
We next calculate principal strain rates from the 2011-2012 GPS network data from Andrews et al. (2014,2018). These data 375 were recorded at 15-second intervals at eleven GPS stations separated by 2 to 20 km within a local network. Thus, these data are temporally precise, but spatially coarse (Table 1) finding that 25% of lakes in Pâkitsoq drained rapidly, suggesting that the true ratio is higher than our result.
For insight, we look to slow-draining lakes in our dataset. We incorporated the cloud state during the lake drainage into our confidence ratings (high, medium, and low; Sect. 2.2.2). For instance, when clouds obscured observations for >6 days during 530 a lake drainage, we classified the lake drainage as slow but assigned it medium confidence to denote the possibility that it may have been a fast drainage. Our dataset contains N = 49 medium-confidence slow-draining lakes. These average 98 m higher in elevation than the high-confidence fast-draining lakes. If fast-draining and slow-draining lakes in our study area were to be sited at indistinguishable elevations (p > 0.05), to match the results of Cooley and Christoffersen (2017), this would require that the 26 (53%) highest-elevation medium-confidence slow-draining lakes (1110 m < s < 1410 m) were actually fast-draining. While 535 possible, this seems unrealistic, as cloudiness should affect identification of slow-draining lakes all elevations (530 m < s < 1410 m) equally. Instead, we randomly selected a subset of the medium-confidence slow-draining lakes and reassigned them as fast-draining. We found that the elevations of fast-draining lakes were still significantly different from slow-and non-draining lakes, with average p = 0.008. In fact, we were not able to match the indistinguishibility found by Cooley and Christoffersen (2017) to p > 0.05 in ten million sets of random selections.

540
Although we represented uncertainty due to cloudiness, we were unable to reproduce the findings of Cooley and Christoffersen (2017) that lake elevation is not correlated to drainage type. Our results agree with previous findings that fast-draining lakes occupy lower regions of the ice sheet (Selmes et al., 2013;Morriss et al., 2013;Miles et al., 2017).

Moulin background strain rate and lake drainage timing
There is suggestive observational evidence for lake-drainage cascades, whereby water supplied to the subglacial system by We investigated whether lakes in more-extensional regimes drained first and triggered regional cascades but found that this was not likely. In some years, higher-elevation lakes were the first to drain: notably 2007, 2009, 2015 (Fig. 8), but the winter 550 strain rates at these lakes were not different from the population average of lake-draining moulins. Background principal strain rates can influence the orientation of moulin-forming fractures in the bottoms of fast-draining lakes (Chudley et al., 2019), but we found no significant relationship to lake drainage date (Fig. 8).
Our results are consistent with those of Williamson et al. (2018b), who tested whether fast-draining lakes were statistically different from slow-draining and non-draining lakes with respect to a full suite of parameters, including winter strain rates Indeed, the earliest fast-draining lake in Pâktisoq each year is not generally sited in a background extensional area, but rather most often a lower-elevation lake (Fig. 8). Earlier access of meltwater to the bed at lower elevation is well documented 565 (McGrath et al., 2011;Bartholomew et al., 2011;Doyle et al., 2013), but such access is often through crevasse fields or streamfed moulins, not necessarily through moulins fed by fast-draining lakes (Catania et al., 2008). Indeed, the first fast-draining lake of each year occurred 10-60 days after the day of first runoff recorded at Swiss Camp (Fig. 8), whereas melt-induced ice flow has been observed 3-10 days after melt onset in our study area (Andrews et al., 2018). This contrast suggests that water reaches the bed via other means before fast lake drainages begin.

.1 Progress in mapping current features versus identifying predictive metrics
A multitude of intensive, independent, and international efforts over the past decade has yielded multiple reliable algorithms and workflows that date and detect lake drainages from satellite imagery (e.g., McMillan et al., 2007;Sundal et al., 2009;Selmes et al., 2011;Liang et al., 2012;Morriss et al., 2013;Johansson et al., 2013;Fitzpatrick et al., 2014;Williamson et al., 575 2017). A related remote-sensing analysis push has yielded techniques that produce detailed maps of lakes and streams on the ice-sheet surface, which are readily applied to locate moulins at stream terminations (Yang and Smith, 2013;Smith et al., 2015;Yang and Smith, 2016;King et al., 2016;Andrews et al., 2018;Yuan et al., 2020). These sophisticated remote-sensing observations of supraglacial lakes, streams, and moulins cover the entirety of Southwest Greenland and extend more than a decade in time. The high data volume should have utility for identifying patterns in lake-drainage date and moulin location, as 580 attempted in this study and by Williamson et al. (2018b). However, beyond coarse trends in lake elevation, which are disputed (Cooley and Christoffersen, 2017), runoff timing ( Fig. 9), and lake size (Williamson et al., 2018b), we find no lake or local ice-sheet properties that are predictive of drainage date. Thus, despite the community's state-of-the-art near-real-time analysis capabilities of ever-increasing image data, we regard a "more data" approach to identifying predictive lake-drainage metrics with pessimism. 585 4.2.2 Current remote-sensing products for regional-scale strain rates Precise measurements from GPS networks show short-lived, high-magnitude strain rate events that precede fast lake drainages (e.g., Fig. 10; Stevens et al. (2015)). Recent research has provided a mechanistic understanding that these events induce lake drainage (Tedesco et al., 2013;Stevens et al., 2015;Hoffman et al., 2018;Christoffersen et al., 2018). The GPS network data shown in Fig. 10 are typical, however, in that their spatial resolution is limited (∼4 km in our example). Notably different is the km, as used by Stevens et al. (2015), is likely necessary. This would increase the number of stations required by an order of magnitude. Design of any such GPS network will require careful consideration of the trade-offs in spatial resolution, spatial coverage, and the cost and feasibility to install and maintain stations in the fidgety conditions of ice-sheet ablation zones.

Parameterizing moulins in ice-sheet models
Ice-sheet models, which are based around detailed physical representation of ice flow, are central tools for state-of-the-art 665 projections of ice-sheet mass balance (Nowicki et al., 2016(Nowicki et al., , 2020. Societal ability to accurately predict global sea levels in future climates depends, in part, on the fidelity of ice flow patterns predicted by these models. This, in turn, depends on sliding speeds (Rückamp et al., 2020) via basal hydrologic conditions, which evolve according to the spatial distribution and timing of melt water delivery to the subglacial system (e.g., Banwell et al., 2016). State-of-the-art ice-sheet projections, such as those from ISMIP6, currently do not explicitly include meltwater inputs to the bed (Nowicki et al., 2020), and most current meltwater accumulation data. They forecasted the date of meltwater input to the bed as the day the crevasse depth reached the full ice thickness. Overall, existing methods to locate moulins in space are largely reliant on present-day observations, and some current methods to construct time-varying basal input hydrographs ignore the episodic nature of moulin activation. Thus, we currently lack a method to accurately represent moulin locations and hydrographs in models of ice sheets in future climates.
We propose a potential path forward for the spatial and temporal representation of moulins in ice-sheet models. First, we 685 look at spatial forecasting. At low elevations, surface-to-bed connections are ubiquitous via stream-fed moulins (Fig. 3) and crevasse fields (McGrath et al., 2011). At high elevations, most moulins are associated with lakes (Smith et al., 2015;Yang and Smith, 2016). Lake locations are static (Lampkin and VanderBerg, 2011), and near-future lake locations, in topographic lows that are currently above typical elevations with runoff, are predictable from current DEMs (Leeson et al., 2014;Ignéczi et al., 2016). Finally, locations of lakes in the distant future, or in any scenario in which ice-sheet size and geometry are 690 substantially different from today, are predictable in ice-sheet models run with sufficient resolution (<∼1 km) to resolve local topographic depressions. If lake locations are forecastable, what remains is predicting whether the lakes will drain to local moulins or to moulins far downglacier, via commonly observed long outflow streams (Poinar et al., 2015). For this, we look to Stock (2020), who used an efficient transfer-function model (Gudmundsson, 2003(Gudmundsson, , 2008 to quantify the probability of moulin is applicable to arbitrary ice-sheet configuration, so it could be used to generate moulin locations in future ice-sheet conditions. Next, we look toward forecasting the date of activation of an arbitrary moulin. Important metrics are elevation on the icesheet surface and cumulative runoff, although together these explain less than half of the variance in the date of fast lake drainage (Sect. 3.5.2). The remaining variance in moulin activation date may lie in strain-rate precursors (Christoffersen et al., 2018) instigated by variations in basal hydrology. Resolving these precursors is within reach of future coupled subglacial 700 hydrology -ice sheet models, although they will require fine spatial (<∼1 km) and temporal (<∼1 day) resolution. Elevation, runoff, and time varying strain rates are potential predictor variables that are readily available within ice-sheet models and their climate forcings; thus, there is potential to construct accurate basal water input hydrographs even in future ice-sheet states.
The factors that determine the spatio-temporal distribution of moulins are finer than spatial and temporal resolution of most ice-sheet-scale models. Thus, we anticipate that some parameterization will be needed to achieve computational feasibility; 705 we envision a probabilistic or ensemble approach to representing meltwater delivery to the bed. For instance, in the absence of discrete moulin locations, moulin density may suffice (Banwell et al., 2016). The spatial density of moulins could be approximated using bedrock geometry, ice-sheet thickness, and slip ratio (Stock, 2020), all of which are common ice-sheetmodel variables, as well as less well represented parameters related to local basal water perturbations throughout a melt season.
These hydrological quantities could potentially be calibrated using currently available image-or field-based surveys of moulin 710 locations, then applied to future climate scenarios and ice-sheet geometries. Ultimately, a probabilistic map of moulin location and activation date within the melt season could be supplied to subglacial hydrology models. In time, coupled ice sheetsubglacial hydrology models, driven with the parameterized basal input hydrographs we envision, could supply better-informed projections of ice flow variability in future climates.

715
In this study of 78 lakes in the Pâkitsoq area of western Greenland over 18 melt seasons, we identified 527 lake-drainage events with high confidence. We found that lake elevation and date of melt onset, but not the annual-average strain rate, were meaningful predictors of the timing of fast lake drainage (N = 212). We also investigated the hypothesis that transient, highmagnitude strain rate events regularly precede fast lake drainages, but fundamental issues with data coverage and temporal resolution limited our ability to make significant conclusions.

720
Our study is motivated by ongoing rapid development and application of ice-sheet models that forecast likely states of the Greenland and Antarctic Ice Sheets in future climate scenarios. Especially for Greenland, the coming generation of ice-sheet models will need to parameterize the timing and locations of melt delivery to the bed. State variables like ice thickness and strain rates are attractive potential parameters because these are available at every model time step. However, our finding that background strain rates do not significantly correlate with lake drainage date implies that this is not likely to provide a path 725 forward. Better metrics for predicting lake-drainage dates are elevation on the ice-sheet surface (lakes associated with lower moulins drain first), and meteorological data such as the date of surface melt onset (lakes drain earlier in years with early melt). WorldView images (light purple), when available, to locate the moulins that drain the lakes in our study area each year. Ideal sensor for regional-scale lake-drainage-> 30 km < 1 km >2 months/yr ≤ 24 hours inducing strain rate perturbations