Controls on Water Storage and Drainage in Crevasses on the Greenland Ice Sheet

Surface crevasses on the Greenland Ice Sheet (GrIS) capture nearly half of the seasonal runoff, yet their role in transferring meltwater to the bed has received little attention relative to that of supraglacial lakes and moulins. Here, we present observations of crevasse ponding and investigate controls on their hydrological behavior at a fast‐moving, marine‐terminating sector of the GrIS. We map surface meltwater, crevasses, and surface‐parallel stress across a ∼2,700 km2 region using satellite data and contemporaneous uncrewed aerial vehicle (UAV) surveys. From 2017 to 2019 an average of 26% of the crevassed area exhibited ponding at locations that remained persistent between years despite rapid advection. We find that the spatial distribution of ponded crevasses does not relate to previously proposed controls on the distribution of supraglacial lakes (elevation and topography) or crevasses (von Mises stress thresholds), suggesting the operation of some other physical control(s). Ponded crevasse fields were preferentially located in regions of compressive surface‐parallel mean stress, which we interpret to result from the hydraulic isolation of these systems. This contrasts with unponded crevasse fields, which we suggest are readily able to transport meltwater into the wider supraglacial and englacial network. UAV observations show that ponded crevasses can drain episodically and rapidly, likely through hydrofracture. We therefore propose that the surface stress regime influences a spatially heterogeneous transfer of meltwater through crevasses to the bed of ice sheets, with consequences for processes, such as subglacial drainage and the heating of ice via latent heat release by refreezing meltwater.

Crevasses are an important pathway for the transfer of water to the subglacial environment of glaciers and ice sheets, while water itself can drive the propagation of crevasses via hydrofracture (Alley et al., 2005;Krawczynski et al., 2009;van der Veen, 2007;Weertman, 1973). Once full-depth hydrofracture has occurred, water flow forms an efficient route for continued meltwater delivery to the bed in the form of moulins. Many studies of the Greenland Ice Sheet (GrIS) have largely focused on supraglacial lake drainage as the primary method of routing surface meltwater to the ice sheet bed (Banwell et al., 2016;Christoffersen et al., 2018;Hoffman et al., 2018). However, supraglacial lakes deliver less total meltwater volume to the ice sheet bed than do crevasse fields, which may capture as much as half of seasonal surface runoff (Koziol et al., 2017;McGrath et al., 2011). However, the limited studies available of crevasse field hydrology describe variable, sometimes mutually contradictory, drainage processes. Some studies observe discrete crevasse drainage (Cavanagh et al., 2017;Lampkin et al., 2013), that appears to result from episodic full-depth hydrofracture, displaying parallels to supraglacial lake drainage events. In contrast, other studies suggest crevasse fields continuously, but inefficiently, transmit a low water flux to the subglacial system without the need for fulldepth hydrofracture McGrath et al., 2011). So far, no studies account for the full spectrum of observations and assumptions surrounding the routing of water through crevasse fields.
As well as a paucity of observational studies, crevasse field hydrology has largely been omitted from numerical models of GrIS surface hydrology (e.g., Arnold et al., 2014;Banwell et al., 2013Banwell et al., , 2016, to be included in only the most recent studies (e.g., Clason et al., 2015;Koziol & Arnold, 2018;Koziol et al., 2017). Where crevasse hydrology is included, its presence is predicted using spatially homogeneous stress thresholds. Indeed, the use of simple stress thresholds to predict hydrological behavior is common across studies of GrIS supraglacial hydrology (Clason et al., 2015;Everett et al., 2016;Koziol et al., 2017;Poinar et al., 2015;Williamson, Banwell, et al., 2018;Williamson, Willis, et al., 2018). However, these thresholds are identified from observational studies performed with the aim of identifying suitable predictors of crevasse presence, not of crevasse hydrology (Hambrey & Müller, 1978;Harper et al., 1998;van der Veen, 1998;Vaughan, 1993); they may not even be suitable for that purpose (Mottram & Benn, 2009), as ice fracture is increasingly understood to be complex and multi-dimensional (Colgan et al., 2016;Hubbard et al., 2021;van der Veen, 1999). To date, no observational studies exist to support the use of any such controls, stress or otherwise, on the hydrological behavior of crevasses.
Here, we use remote sensing to identify the diversity in, and controls on, the hydrological behavior of crevasses on a fast-moving sector of the GrIS. We make use of large-scale, satellite-derived data to investigate the spatial variation of crevasse hydrological state across a ∼2,700 km 2 area between 2017 and 2019. Using this data, we identify the spatial heterogeneity in crevasse hydrological behavior; the interannual variability of crevasse fields; and we test potential topographic and dynamic controls on such behavior. We supplement these regional observations using repeat surveys of a ∼7 km 2 sector of a fast-flowing crevasse field from an uncrewed aerial vehicle (UAV). These high spatial and temporal resolution observations allow for the identification of filling and drainage processes occurring at the scale of individual crevasses, as well as the local validation of satellite observations.

Study Area
Our study area is a ∼2,700 km 2 sector of the western GrIS (Figure 1), which comprises six marine-terminating outlets extending from Sermeq Kujalleq (also known as Store Glacier; 70.4°N 50.6°W) in the south to Perlerfiup Sermia (71.0°N, −50.9°W) in the north. Within this large-scale region of interest (hereafter the "satellite ROI"), we used UAV surveys and Structure-from-Motion with Multi-View Stereo (SfM-MVS) photogrammetry to assess, at high resolution, a crevasse field in the Store Glacier drainage basin, 25 km from the calving front (hereafter the "UAV ROI"). The UAV ROI is 1.5 km wide and 5 km long (Figure 1), and was chosen based on its coverage of an initiating crevasse field, thereby including areas that ranged from hosting no visible crevasses through to areas where crevasse width exceeds 50 m.

Crevasse Classification
Crevasse identification from digital elevation models (DEM) can be approached in a variety of ways (Florinsky & Bliakharskii, 2019), but we use a simple method that identifies crevasses from the residuals between the original and a smoothed DEM (Kodde et al., 2007;Ryan et al., 2015). We output a binary crevasse mask of the satellite ROI using ArcticDEM v3 mosaic data at 2 m resolution (Porter et al., 2018), processing the data in Google Earth Engine ( Figure S1; Gorelick et al., 2017) to allow for efficient computation over the ∼2,700 km 2 study area. We first cropped the ArcticDEM to the GIMP ice mask (Howat et al., 2014), before smoothing the elevation model by performing an image convolution with a circular kernel of 50 m radius (manually selected after inspecting 10-100 m outputs). Residuals greater than 1 m between the smoothed and raw elevation values were identified as crevasses. To compare with glaciological stress estimates, the 2 m dataset was aggregated into grid cells to match the resolution (200 m) and projection (NSIDC sea ice polar stereographic north) of the Making Earth System Data Records for Use in Research Environments (MEaSUREs) 2018 ice sheet surface velocity grid. Aggregated values represent the fraction of grid cell area classified as crevasses. When partitioning into crevassed and non-crevassed regions of the ice sheet, we define a "crevassed" pixel as that with >1% crevasse coverage. This threshold value was chosen to ensure that grid cells were still classified even where crevasses were small (e.g., upstream crevasse fields), whilst still -based upon visual inspection -being sufficient to filter out grid cells with false positive pixels.
We also classified crevasse initiation zones by manually locating the upstream boundary between crevasse fields and bare ice from the 2 m crevasse dataset. The flow direction was determined from 2018 MEaSUREs velocity data and a 200 m buffer around the linear boundary was used to identify pixels in the dataset that should be classified as being in crevasse initiation zones.

Water Classification
We produced a single binary map of water presence across the satellite ROI through the 2017-2019 melt seasons using Sentinel-2 imagery in Google Earth Engine ( Figure S2). These seasons were selected based on the availability of Sentinel-2 data, with 2018 specifically selected to match the timing of the UAV surveys on Store Glacier. We first identified all Sentinel-2 scenes with <40% cloud cover and <70° solar zenith angle from May to October of each year. We clipped the images to the GIMP ice mask (Howat et al., 2014) and converted digital number values to top of atmosphere reflectance. The latter have been shown to be suitable for identifying surface water on the GrIS from medium-resolution optical imagery (Pope, 2016), and have been used for surface water classification in Sentinel-2 data (Williamson, Banwell, et al., 2018). We then calculated the normalised difference water index (NDWI) from bands 2 (blue) and 4 (red) for all images. Following Williamson, Banwell, et al.'s (2018) parameterization for the Store Glacier region (falling within the satellite ROI), we used an NDWI threshold of 0.25 to create binary water classification maps for each Sentinel-2 image. In order to avoid false positive identification of shaded regions, we masked areas in topographic shadow, using the ArcticDEM (resampled to Sentinel-2 resolution and projection) for topography and the solar zenith angle from Sentinel-2 image metadata. Finally, we summed the image stack to count the number of times through the 2018 melt season that a pixel was identified as water. To reduce false positive classifications (e.g., cloud shadow, ephemeral slush zones at the beginning of the melt season) we classify as water any pixel that was identified as water in ≥2 images through the melt season. As for crevasse maps, we aggregate this data onto the velocity grid with a unit of fractional coverage of water within each grid cell. When partitioning the ice sheet surface into "ponded" and "unponded" regions (i.e., where water was or was not, respectively, observed at the surface of crevassed regions), we define a ponded pixel as that with >1% water coverage (following our 1% threshold for crevassed pixels), and furthermore a "lake-filled" pixel as that with >95% water coverage. The latter value was selected as an appropriate classification threshold for lakes by comparing a range of thresholds against the high-resolution annual water-presence maps.

Topographic Analysis
In order to explore the extent of topographic controls on crevasse surface hydrology, we identified topographic sinks that would capture the surface flow of water by filling closed depressions in the ice sheet surface, similar to previous studies of potential lake sites elsewhere on the GrIS (Ignéczi et al., 2016). Before filling, the ArcticDEM data were resampled to the resolution and projection of the velocity data. This process removed false depressions at high spatial scales and allowed intercomparison with other data.

Stress Analysis
Strain rate and stress thresholds are commonly used to predict crevasse formation and supraglacial lake drainage from estimates of surface velocity (e.g., Christoffersen et al., 2018;Poinar et al., 2015;Stevens et al., 2015). Inferring stress from velocity-derived strain rate requires an additional estimate or assumption of ice temperature, but allows for comparability of critical failure thresholds between ice masses of varying temperature (Colgan et al., 2016;Vaughan, 1993). We calculate simple stress thresholds previously proposed to control surface water routing in lake drainage studies, and assess their applicability to crevasse ponding. Specifically, we estimated the stress in the first and second principal directions (as applied by Poinar et al., 2015;Williamson, Banwell, et al., 2018;Williamson, Willis, et al., 2018), as well as the von Mises yield criterion (as applied by Clason et al., 2015;Everett et al., 2016;Koziol et al., 2017;Williamson, Willis, et al., 2018), using Glen's flow law as the constitutive equation linking ice stress and strain rate. As a proxy for whether the dominant stress regime is extensional or compressive, we further calculated the mean surface-parallel stress from the first and second principal stresses.
Surface strain rates were derived from MEaSUREs gridded GrIS annual velocity data for 2018 (Joughin et al., 2010), with extensional strain rates defined as positive. The surface strain rate tensor ε ij is calculated from the surface-parallel components of velocity, u and v (in NSIDC Polar Stereographic North grid directions x and y), as Derivatives were approximated using the finite difference of the MEaSUREs velocity field. The stress tensor, σ ij , was calculated following the Nye-Glen isotropic flow law (Glen, 1955;Nye, 1957) as where n is the flow law exponent with value 3 and B is a viscosity parameter, which we assign a value of 324 kPa a 1/3 (Cuffey & Paterson, 2010) based on an assumed 10 m ice temperature of −5°C. This uniform temperature assumption follows that made in other regional studies examining similar spatial scales (Clason et al., 2015;Koziol et al., 2017;Williamson, Willis, et al., 2018), and matches observations made using distributed temperature sensing measurements at Store Glacier . The effective strain rate, ε e , was then calculated following Cuffey and Paterson (2010) as Because only surface-parallel stresses are considered, σ ij can be expressed by two principal stresses. The first surface-parallel principal stress, σ 1 , was calculated as the highest (most extensional) eigenvalue of the stress tensor, and second surface-parallel principal stress, σ 2 , as the lowest (most compressive) eigenvalue.
The von Mises stress, σ V , was calculated from the surface-parallel principal stresses following Vaughan (1993) as A von Mises failure envelope was prescribed as the 95th percentile of the von Mises stress across the non-crevassed area, allowing up to 5% of the data to be misclassified (Vaughan, 1993).
The mean surface-parallel principal stress, σ m , was calculated as the arithmetic mean of the first and second surface-parallel principal stresses as To explore the interaction between stress state and crevasse behavior, we followed Vaughan (1993) in presenting data in the form of failure maps, presented in terms of the two surface-parallel principal stresses.

UAV Photogrammetry
We acquired aerial imagery over a 13-day period in July 2018 (Table S1) using a custom-built, fixed-wing UAV with 2.1 m wing span. Imagery was collected using a Sony α6000 24 MP camera with a fixed 16-mm lens, processed using Structure from-Motion with Multi-View Stereo (SfM-MVS) photogrammetry, and used to derive velocity fields within the UAV ROI as described by Chudley, Christoffersen, Doyle, Abellan, and Snooke (2019). In brief, photogrammetry was performed using AgiSoft Metashape v.1.4.3 software, and geolocated by using an on-board L1 carrier-phase GPS unit (post-processed against a local on-ice ground station) to locate the position of aerial photos (Chudley, Christoffersen, Doyle, Abellan, & Snooke, 2019). Outputs from the photogrammetric process were 0.15 m resolution orthophotos and 0.2 m resolution DEMs.

Surface Classification
To date, UAV-based crevasse detection has been based on topographic analysis of DEMs (e.g., Florinsky & Bliakharskii, 2019;Ryan et al., 2015). Although useful from a hazard assessment perspective (Florinsky & Bliakharskii, 2019), DEM-based methods alone cannot be used to identify features, such as ponded or healed crevasses, while crevasse detection is also sensitive to threshold choice and ultimately DEM resolution (Florinsky & Bliakharskii, 2019;Jones et al., 2018). To take advantage of the high spatial resolution and multi-dimensional outputs of UAV surveys, we used a combination of object-based image analysis (OBIA) and supervised classification to identify crevasses and their hydrological state in a survey of the UAV ROI flown on 2018-07-08. OBIA is based not on the numerical characteristics of individual pixels but on objects, that is, groups of meaningfully similar pixels segmented according to spectral homogeneity (Blaschke, 2010). This has been used successfully in a glaciological context by Kraaijenbrink et al. (2016Kraaijenbrink et al. ( , 2018 for mapping cliff/ pond systems and emissivity on a debris-covered glacier. We again used Google Earth Engine to perform the full segmentation and supervised classification workflow ( Figure S2).
We identified a number of variables that could be used as inputs for a supervised classification algorithm to identify crevasse field surface features. This included: the red, green, and blue values of the orthophoto ( (Figures 2d and 2j); the DEM slope, which effectively highlighted small crevasses with widths on the order of a few meters (Figures 2e and 2k). Following Kodde et al. (2007), we also used DEM values that were black-top-hat filtered with a 30 m structuring element that was useful in identifying large crevasses with widths on the order of tens of meters (Figures 2f and 2l). A black top-hat filter morphologically closes the glacier surface at scales smaller than the structuring element, before subtracting the closed surface from the original data. This process was performed in MatLab prior to input into Google Earth Engine.
We performed image segmentation using Simple Non-Iterative Clustering (SNIC; Achanta & Süsstrunk, 2017), a computationally efficient implementation of superpixel-based clustering. Rather than segmenting an image into semantically meaningful objects, superpixel-based segmentation aims to simplify the image into small, uniform, and compact clusters of similar pixels ("superpixels"), with a focus on boundary adherence. The variables described above were used as the input to the segmentation algorithm. We manually selected a seed spacing of 15 pixels (2.25 m) and a (relatively high) compactness factor of 10.1029/2021JF006287 6 of 18 200. This resulted in superpixels small enough to display strong boundary adherence to small and healed crevasses at the scale of meters, whilst still clearly delineating the margins of larger features such as water bodies. As an input to the supervised classification, we calculated the average and standard deviation of values in each superpixel from the variables described above, as well as the perimeter-to-area ratio of the superpixel, and normalised the results. We adopted a supervised approach to surface classification (Kraaijenbrink et al., 2016(Kraaijenbrink et al., , 2018Ryan et al., 2018) by training a random forest classifier. Random forests make use of an ensemble of decision trees that classify objects by applying a series of if-then-else logical conditions determined by training data. By utilizing an ensemble of independent decision trees, random forests aim to avoid overfitting that may occur when using a single tree. Each tree utilizes a randomised subset of training data, and the final result is gathered from a majority vote. To reduce the amount of redundant information used to train the classifier, we performed a non-parametric mutual information (MI) test on our training data as a proxy for the predictive power of each input variable. Rejecting input variables beneath the median MI value ( Figure S3) did not noticeably reduce the accuracy of the output data ( Figure S4). Therefore, we used only the nine most significant variables as inputs to the classifier. We manually constructed training datasets of 90 points each, distributed across the ROI, for six distinct surface types: bare ice, snow, healed crevasses, "small" crevasses, "large" crevasses, and water. We separated "small" and "large" crevasses (those with widths of meters vs. tens of meters) into two training datasets as they displayed distinctly different values for properties such as brightness, slope, and the top hat filtered DEM (Figure 2). We trained the random forest classifier with 128 trees on two-thirds of the dataset (60 points per classification) and retained one-third (30 points per classification) for validation. Output classification performed well visually ( Figure S5) and validation data showed that a >95% F 1 accuracy score was observed for all surface types ( Figure S4 and Table S2), apart from for snow and bare ice, which for our purposes were not important to distinguish. Although we identified six surface types, for this analysis we were only interested in three distinctions: crevasses (combining "small" and "large" crevasses), ice (combining bare ice, snow, and healed crevasses), and water.

Ponded Crevasse Distribution
From the ArcticDEM elevation model and Sentinel-2 optical imagery, we mapped the distribution of crevasses ( Figure 3a) and surface water (Figure 3b) across the study region. Of the total ice area assessed, ∼34% (∼960 km 2 out of ∼2,695 km 2 ) was classified as being crevassed (i.e., where a 200 m 2 grid cell has a crevasse fraction >1%). On average, 26% of this total crevassed region was observed to exhibit surface ponding (i.e., where a 200 m 2 grid cell has a water fraction >1%) across the 2017-2019 ablation seasons, ranging from 20% to 30% (Table 1). Although total ponded area varies from year-to-year, the spatial pattern of ponded crevasses persisted between years. This can be observed qualitatively (Figures 3i-3iii), and is supported by Pearson's correlation coefficient tests, which return statistically significant (   some areas, suggesting the ability of crevasses to pond is externally controlled rather than a property of individual crevasses.

Topographic Analysis
If the ability of a crevasse to pond was depended solely on meltwater availability, ponded crevasses would be more prevalent at lower elevations, where air temperature is higher. We tested this hypothesis by comparing the spatial distribution of satellite-derived water fraction (Figure 3b) with elevation. Within crevassed grid cells, there was no significant relation between elevation and water fraction (R 2 value = 0.0014, p < 0.01). This analysis indicates that (air-temperature driven) meltwater availability does not exert a major control on ponded crevasse formation up to the limit of the elevation range of the satellite ROI (∼1,500 m above sea level).
If crevasse systems pond due to receiving meltwater that has been transported laterally across supraglacial drainage networks, ponding should occur at the bottom of surface depressions. We tested this hypothesis by assessing the prevalence of meltwater ponding within surface depressions as identified from ArcticDEM data ( Table 2). This analysis indicated that 78.9% of grid cells classified as lakes occurred in topographic depressions, whilst only 5.0% of crevassed grid cells and 9.7% of ponded crevasse grid cells were similarly located. Thus, while supraglacial lakes were predominantly located within surface basins (as expected), ponded crevasse fields were predominantly located outside such basins (Figure 4). Infilling of surface basins by lateral supraglacial water transport therefore appears to explain only a minority of crevasse ponding locations.

Stress Analysis
We present strain-rate-derived stresses as failure maps, plotted in the form of density plots ( Figure 5). To aid visualisation, each data point is plotted twice with assignments of σ 1 and σ 2 reversed, giving symmetry across the line σ 1 = σ 2 . Based on the stress distribution of the non-crevassed area, the von Mises failure envelope was prescribed at 76 kPa (marked with ellipses in Figure 5). However, this threshold does not differentiate either crevasse incidence or hydrological status. Crevasses plot both inside and outside the von Mises failure envelope (Figure 5b), while initiating crevasses (Figure 5c) plot predominantly within the envelope. The von Mises failure envelope is also not useful for differentiating crevasse ponding (Figure 5d). Inspection of Figure 5 does, however, reveal that initiating crevasses ( Figure 5c) and ponding crevasses (Figure 5d) are separated by the line defined by σ 1 = −σ 2 (dashed line in Figure 5). This line marks the transition in mean surface-parallel stress state (σ m ) from negative (compressional; below the line) to positive (extensional; above the line). Thus, this analysis reveals that 89% of initiating crevasses are located in areas of extensional stress (cf. 54% of all crevasses). In contrast, 68% of ponded crevasses are located within areas of compressive stress (cf. 46% of all crevasses). An example of the contrasting relations between Surface type Grid cell classification thresholds

Proportion of surface type within depressions [%]
Crevasse >1% crevasse fraction 5.0 Ponded crevasses (excluding lakes) >1% crevasse fraction and 1%-95% water fraction 9.7 Lake >95% water fraction 78.9  crevasse status and mean surface-parallel stress are plotted spatially in Figure 6. While not all crevasses that are located in compressive regimes pond, the transition between extensional and compressive σ m regimes represents a convincing boundary between regions of crevasse initiation and ponding. This suggests that a compressive mean stress regime is a necessary, but not sufficient, condition of ponding.

UAV-Based Image Analysis
UAV-derived observations of crevasse initiation and ponding (Figure 7a) follow similar patterns to those revealed by the regional-scale satellite data (Figure 7b). Crevasses initiated-or at least become identifiable in the decimeter-resolution data-in the upstream section of the study zone. As they are advected down-glacier, crevasses width increased from decimeters to a maximum of ∼10-60 m (Figures 7a, 7i, and 7ii). The region of crevasse initiation was coincident with a zone of extensional σ m in the satellite-derived data (Figure 7b). In the downstream sector of the UAV ROI, crevasse size remained relatively stable, but displayed a higher propensity to pond in the down-glacier direction (Figure 7a). This region of crevasse ponding occurs where satellite-derived σ m is observed to be compressive (Figure 7b).
Repeat UAV surveys provide insight into processes occurring at the scale of individual crevasses. Over the 13-day period in July 2018 over which surveys were undertaken (Table S1), three crevasse systems in the UAV ROI were observed to drain, and six underwent significant filling. Crevasse drainages appear to be rapid. Of the three drainages identified, two represent crevasses that had a constant or rising water level in sequential imagery acquired prior to drainage. All three exhibited significant water loss between subsequent adjacent surveys (e.g., Figures 8a-8d). One crevasse system lost a substantial volume of water in less than 24 h (Figures 8a and 8b), and water levels continued to drop for the rest of the 12 days survey period ( Figure S6b and S6c). This suggests that either a moulin had formed, and that water therefore continued to drain into the subglacial system, or that small open fractures continued to transfer water inefficiently into the englacial system. The filling crevasses were clustered tightly at the upstream side of the ponded crevasse system, in a location where crevasses are advecting from an extensional to compressive σ m regime. Of the three crevasse drainage events identified, two occurred within the larger ponded system and compressive mean stress regime, while one occurred in a smaller crevasse at the periphery of the system, in a weakly extensional regime. These observations are consistent with the satellite-based observations in that ponded crevasses are observed to occur within the compressional mean stress regime (Figure 7c), but also show that drainages of these ponded crevasses occur discreetly and rapidly.
Within the UAV-derived data, our observations do not indicate significant lateral meltwater routing. Where supraglacial streams routed water between ponded crevasses, they were easily identified in imagery (Figure S6a), but this was not common across the UAV ROI. When individual crevasses drained, observations of any consequential effects on the surrounding system was limited. For the two most prominent crevasse drainages (Figure 8), we identify the adjacent crevasses that also drained, either through visible supraglacial networks (marked "S" in Figure 8), or without visible supraglacial connections, which we thus interpret to be connected englacially (marked "E: in Figure 8). In the first case (Figure 8a), an overflowing crevasse system formed local supraglacial networks, and after one crevasse had drained, water levels across the entire network dropped (Figure 8b). However, only one adjacent crevasse drained without visible surface routing (Figure 8b). In the days following the drainage event, incised supraglacial channels formed between the previously overflowing system ( Figures S6b and S6c). In the second drainage case, an individual crevasse drained without affecting water levels in the adjacent crevasses at all (Figures 8c and 8d). As such, in the ponded crevasse system within the UAV ROI, we do not observe hydrological responses to drainage events that extend more than one to two crevasses (∼100-200 m) from the initiating crevasse.

Performance of Satellite-and UAV-Based Classification Methods
The object-based random forest classification of UAV data (Figure 7a) enabled the identification of crevasses and surface water with >95% accuracy ( Figure S4 and Table S2). Comparison of the UAV and satellite-derived surface classifications (cf. Figures 7a and 7b) shows clear agreement in the distribution of surface features. The cutoff width below which crevasses are unable to be identified from ArcticDEM v3 data is ∼10 m (∼5 pixels). Although this cutoff means the satellite data do not identify the smallest crevasses (such as those in Figure 7ii), the resolvable size of a crevasse is approximately equal to the resolution of the Sentinel-2 bands used for NDWI calculation (10 m), allowing the two satellite-derived datasets to be compared directly. Despite the ArcticDEM mosaic being derived from multitemporal data (individual tiles across the satellite ROI range from 2009 to 2017), the distribution of crevasses is consistent with the 2018 UAV dataset (cf. Figures 7a and 7b). This indicates that, regardless of the advection of individual crevasses (at a rate of ∼650 m a −1 in this area), interannual variation in crevasse field extent is minor across the time period in this study. We can therefore be confident that the ArcticDEM 2009-2017 crevasse distribution can be meaningfully compared to the Sentinel-2 2017-2019 surface water distribution. Sentinel-2-and UAV-derived surface water locations also agree consistently (cf. Figures 7a and 7b), with individual ponded crevasses able to be co-located between the Sentinel-2 and UAV datasets. The Sentinel-2 data additionally identifies ponded crevasses that were not water-filled on the date of the UAV survey. In summary, the UAV data shows that we can be largely confident in our satellite-derived crevasse and water mapping at spatial scales ≥10 m. It is likely that the ArcticDEM-based crevasse mapping underestimated crevasse extent at higher elevations, as the optically derived dataset will not be sensitive to snow-filled crevasses; however, false negative classifications do not affect the conclusions drawn from this dataset.

Spatial Variability in Crevasse Hydrology
Across the ∼2,700 km 2 study area, an average of 26% of the crevassed region (containing crevasses >10 m in width) exhibited visible ponding between 2017 and 2019-although we note that hydrological activity may still be occurring in regions without visible ponding at the surface (see Section 4.4.2). The inter-annual variation in ponded crevasse coverage was lowest in 2017 and highest in 2019, consistent with the lowest and highest reported ice-sheet-wide melt data (e.g., Tedesco & Fettweis, 2020) which identified 2019 as an exceptional melt year. This suggests that the extent of crevasse ponding is in part controlled by melt intensity, consistent with previous models conceptualizing crevasses as linear reservoirs McGrath et al., 2011), that is,, as melt intensity increases beyond the capacity of crevasses to discharge water into the englacial system, more ponding will be observed at the surface. However, inter-annual melt intensity cannot explain the spatial distribution of observed ponding. Only a minority of crevasses are observed to pond, and the spatial distribution of these is consistent from year-to-year (Figures 3i-3iii). These patterns occur on scales <1 km, smaller than any reasonable spatial boundary in melt intensity resulting from surface mass-balance drivers, such as the vertical gradient in air temperature. This inference is supported by statistical tests that rejected the hypothesis that crevasse ponding was more prevalent at lower elevations within the ablation zone, where surface melt intensity is generally higher. These local-scale patterns of crevasse ponding are also stable in space regardless of ice velocity, suggesting that ponding incidence is not advected with individual crevasses. We therefore conclude that likely controls on the incidence of crevasse ponding are: (i) distinct from melt intensity; (ii) not associated with the properties of individual crevasses; and (iii) spatially variable on the order of 10 2 -10 3 m.

Controls on Crevasse Hydrology
Previous studies have predicted the location of current and future supraglacial lakes on the GrIS by identifying depressions in the surface topography that would capture supraglacial meltwater (e.g., Ignéczi et al., 2016). This could be reasonably applied to crevasse ponding if meltwater could be routed, without obstruction, for long distances laterally along hydrological gradients across crevasse fields. This would result in localized crevasse ponding within surface depressions, forming a single surface lake if water supply exceeded open crevasse volume. However, our satellite observations suggest that-at least in the wide, deep and active crevasses examined here-only a small fraction of crevassed area behaves in this way. Only 9.7% of the ponded crevasse region is located within surface depressions (cf. 5.0% of all crevasses and 78.9% of the supraglacial lake area). This suggests that lateral supraglacial transport of water into topographic basins is not the principal cause of crevasse ponding. Indeed, because ponding occurs across topographic highs (Figure 4), we infer that in ponded crevasse systems drainage into the wider supraglacial and/or englacial drainage system is being restricted. This inference is supported by our UAV repeat surveys, which show that hydrological connections between ponded crevasses-whether supraglacial or englacial-are rare and have limited spatial extent, often between only a few crevasses on the scale of 100-200 m (Figures 8b and 8d). Even where hydrological connections exist, they appear to form as a consequence, rather than a cause, of crevasse drainage (Figures S6b and S6c). Our UAV surveys cover only a small temporal window, and hence we could be missing supraglacial pathways that exist at other points in the melt season. However, comparing the UAV-derived ponding extent to the satellite-derived data that covers the whole melt season (cf. Figures 7a and 7b) shows that the crevasse field was close to its maximum capacity when the UAV surveys were undertaken, suggesting a majority of supraglacial pathways were filled. If lateral meltwater flow is restricted, the drainage of water from ponded crevasse systems cannot, for the most part, be caused by water being routed into wider supraglacial networks (and from there to moulins, lakes, or the englacial system). This contrasts with previous studies, which often assume that lateral supraglacial drainage can occur unrestricted across crevasse fields (e.g., Clason et al., 2015;Poinar, 2015).
In our satellite-and UAV-derived datasets, ponded crevasses are largely restricted to regions of compressive mean stress. This association may be explained by supraglacial and englacial drainage remaining open and well-connected across crevasse fields in areas characterized by extensional mean stress regimes. This is consistent with the view of crevasse systems on temperate valley glaciers as hydraulically connected, albeit inefficiently, to englacial and/or subglacial drainage systems through a linked network of fractures (Fountain et al., 2005), and supported by radar observations at Sermeq Kujalleq (Store Glacier), where englacial meltwater storage in crevasse-damaged ice has been inferred down to a depth of 48 m (Kendrick et al., 2018). In these systems, meltwater availability rarely exceeds drainage rate, explaining the limited crevasse ponding observed in such extensional regions. In contrast, we suggest that in compressive regimes, englacial connections undergo what Irvine-Fynn et al. (2011) described as "pinch-off," whereby crevasse closure by ice creep hydraulically isolates crevasse systems from the wider supraglacial and englacial drainage system, resulting in subsequent ponding. This hypothesis is supported by our UAV data, which show that englacial connections between crevasses in ponded regions are limited in prevalence and extent. While a compressive mean stress regime appears to be a necessary condition for crevasse ponding, not all crevasses located within such regimes are water-filled ( Figure 6). This may be because the englacial system does not always close entirely: we note that, at the UAV ROI, the upstream ∼400 m of the compressive stress field does not exhibit ponding (Figure 7b), perhaps indicating a delay before the crevasse system is fully isolated. As the velocity at the UAV ROI is ∼650 m a −1 , this would indicate a closure time of ∼7 months to isolate the system.

Rapid Drainage of Ponded Crevasses
Given that we found little direct evidence for hydrological connections between crevasses in ponded regimes (Section 4.2), we consider full-depth hydrofracture and drainage to the subglacial environment to be a likely mechanism by which ponded crevasses drain (Boon & Sharp, 2003;Krawczynski et al., 2009;van der Veen, 2007;Weertman, 1973). This would be consistent with the rapid and heterogenous drainages observed in UAV data. Analogous to supraglacial lake drainages via rapid hydrofracture (Chudley, Christoffersen, Doyle, Bougamont, et al., 2019;Das et al., 2008;Doyle et al., 2013;Stevens et al., 2015), rapid crevasse drainage events are expected to deliver distinct, isolated pulses of meltwater to the bed. The full dynamic consequences of such events are explored in detail elsewhere (e.g., Nienow et al., 2017), but it is apparent that meltwater inputs to the bed that are rapid (Schoof, 2010) and spatially discrete (Banwell et al., 2016) can influence ice dynamics. For example, rapidly draining crevasse systems at the shear margin of Jakobshavn Isbrae have been shown to deliver meltwater to the bed at sufficient rates and volumes to overwhelm the capacity of the subglacial system (Lampkin et al., 2013), increasing ice mass flux across the shear margin and enhancing glacier discharge (Cavanagh et al., 2017;Lampkin et al., 2018).
There are, however, several features of rapid crevasse drainages that are distinct from more widely studied lake drainage events. After hydrofracture, ongoing meltwater delivery via the newly open moulin represents an important hydrological component of lake drainages (Hoffman et al., 2018;Koziol et al., 2017). In contrast, the smaller catchments of individual crevasses means that this effect is likely less important following crevasse drainage (although crevasses are more numerous than lakes). Furthermore, unlike lakes, it appears to be relatively common that crevasse systems can drain multiple times through a single ablation season (Cavanagh et al., 2017). However, the net effect of this on ice dynamics has yet to be identified. Additionally, crevasses that remain ponded and then refreeze at the end of the season will release latent heat and facilitate ice warming (Lüthi et al., 2015) at depths of up to hundreds of meters .

Slow Drainage of Unponded Crevasses
Approximately 74% of crevasse fields display no evidence of ponded meltwater at the surface, despite there being no difference in local meltwater availability compared to adjacent ponded regions. Since meltwater is inevitably also routed into these unponded crevasses, we suggest that this water is accommodated by the wider supraglacial and englacial hydrological system, stored and/or routed through pre-existing englacial pathways (e.g., Figure S7) or linked fracture networks (Fountain et al., 2005;Kendrick et al., 2018), all of which are maintained by an extensional stress regime. Since no ponding is observed in these regions, the condition for hydrofracture is restricted, meaning that these crevasses are unlikely to route meltwater directly to the bed. Instead, we suggest that this laterally routed meltwater must eventually intersect pre-existing moulins (Catania & Neumann, 2010), terminate at supraglacial lakes, remain as a liquid reservoir (Kendrick et al., 2018), or freeze during the winter season.
We conceptualize that, due to low rates of lateral meltwater transport, drainage in unponded crevasse systems has a long total transit time to the bed. Such slow, continuous crevasse drainage has previously been applied in crevasse hydrological models at the GrIS by Colgan et al. (2011) and McGrath et al. (2011 suggested crevasse surface-to-bed delivery rates may be 200-fold slower than moulins (∼12 h for a 0.1 m wide crevasse cf. ∼1 h for a 1 m 2 moulin), whilst McGrath et al. (2011) suggested that crevasses may slow englacial drainage to such an extent that a diurnal cycle of meltwater input can be damped to a quasi-steady state discharge on the timescale of hours-days. This slow and sustained delivery of meltwater through crevasses to the glacier bed would be less likely to overwhelm temporarily the transmission capacity of the subglacial drainage system. Therefore, both studies argue that regions of the bed subject to this style of meltwater delivery are less likely to exhibit ephemerally enhanced basal sliding compared to regions experiencing episodic, efficient meltwater pulses (as in Section 4.4.1). Additionally, this slower englacial drainage style associated with crevasses may have distinct thermal consequences. It has been argued that slow meltwater delivery through crevasses would deliver latent heat that results in more cryo-hydrologic warming relative to regions fed by discrete moulins  because densely packed and slow hydrological pathways increase the volume of ice warmed by the latent heat release of englacial freezing relative to efficient drainage pathways (Lüthi et al., 2015;Phillips et al., 2010). As full-depth hydrofracture is likely restricted in unponded crevasse systems, drainage in these regions may deliver less latent heat into depths on the order of hundreds of meters compared to ponded crevasses, where propagation is facilitated by hydrofracture (e.g., Poinar, 2015).

Implications for Hydrological Routing Models
In the past, regional models of ice sheet hydrology and dynamics have often failed to include crevasse drainage, instead focusing exclusively on supraglacial lake drainage (e.g., Arnold et al., 2014;Banwell et al., 2013Banwell et al., , 2016). Recent regional hydrological models have begun to include crevasse drainage in simple ways. For example, Clason et al. (2015) incorporated crevasse drainage, but considered it similarly to supraglacial lake drainage. These authors identified crevassed regions based on a σ v threshold, which were then allowed to fill and hydrofracture according to a LEFM model (van der Veen, 2007). After full-depth hydrofracture and drainage, a moulin was formed that continued to drain any further meltwater continuously to the bed. More recently, Koziol et al. (2017) allowed crevasses to continuously drain, with meltwater produced at the surface of crevasse fields (again identified according to a σ v threshold) drained immediately to the bed without requiring hydrofracture. These behaviors, reflecting a paucity of observations available at the time, were assumed to be spatially uniform. The observational results we present herein highlight ways, summarized below, in which future studies may be able to account further for a wide diversity of crevasse hydrology while keeping inputs and classifications as simple as possible.
Our first recommendation is to avoid using simple stress thresholds or zero stress models to predict crevasse presence in surface routing models. Several models that incorporate crevasse drainage have used a von Mises yield criterion (following Vaughan, 1993) to estimate the location of crevasses for water routing (Clason et al., 2015;Koziol et al., 2017), as well as to approximate tensile stress in crevasse hydrofracture (Clason et al., 2015;Morlighem et al., 2016). However, our data indicate that von Mises yield stress, σ v , is not the most appropriate predictor of crevasse incidence and hydrology. Indeed, our analysis indicates that crevasses exist across a range of σ v values, both above and below the yield threshold prescribed following the method of Vaughan (1993; Figure 5b). Furthermore, our data indicate that even broad regions of ice failure are not predicted accurately by a von Mises yield criterion (Figure 5c), with initiating crevasses existing predominantly: (a) below the 76 kPa threshold prescribed following Vaughan (1993) and (b) in regions of positive mean stress. This is unsurprising considering that the compressive strength of ice greatly exceeds its tensile strength, a factor that von Mises stress is insensitive to. As such, we do not recommend a von Mises yield criterion as a suitable threshold for identifying regions of crevasse incidence. Additionally, using the alternative measure used by Vaughan (1993), the strain energy dissipation criterion, does not significantly improve the fit to our data (see Text S1 and Figure S8). While alternative measures, such as the first principal (Benn et al., 2007) or longitudinal (Harper et al., 1998) stress, may be more appropriate for predicting brittle ice failure, they are unlikely to also be a useful threshold for predicting crevasse field distribution (Colgan et al., 2016), as crevasses advecting into unviable stress regimes take time to adjust to the new equilibrium (Mottram & Benn, 2009). For studies of contemporary crevasse fields, we instead recommend the use of direct observations as the simplest and most practical way to map crevasse locations. Herein, we present a simple method to achieve this based on the analysis of high-resolution DEMs, but other studies have adopted alternative methods such as using convolutional neural networks to classify optical imagery (Lai et al., 2020). However, using satellite-derived data comes with weaknesses, such as being unable to identify crevasses under snow cover or below the spatial resolution of the input data. The clustering of incipient crevasses in extensional mean stress regimes (Figure 5c) suggests that models could properly map crevasses from surface stresses if a more suitable stress threshold were to be identified and subsequent advection and closure could be accounted for (e.g., Albrecht & Levermann, 2014).
Our second recommendation is to begin accounting for the diverse representations of crevasse drainage. Our data suggests that the delivery of supraglacial meltwater to the glacier bed through crevasses falls into a spectrum of behavior, ranging from episodic rapid drainage via full-depth hydrofracture (Section 4.4.1) to slow and continuous englacial drainage (Section 4.4.2). At the coarsest scale, it may be desirable to implement a binary system to account for the end-members of this spectrum. Our data suggest that the mean surface-parallel stress is a useful way of segregating crevasse hydrological behavior, with crevasses in compressive regimes being hydrologically isolated and exhibiting episodic rapid drainage. In contrast, crevasses in extensional regimes are hydraulically connected, exhibiting continuous drainage into the wider supraglacial and englacal system. Thus, in compressive regimes, drainage could be modeled as episodic rapid hydrofracture following Clason et al. (2015), while also restricting the lateral flow of meltwater between grid cells. In extensional regimes, meltwater could be routed laterally to the nearest moulin or supraglacial lake. Simple thresholding such as this could be used in regional hydrological models to investigate the seasonal and long-term effects of spatial heterogeneity in crevasse hydrology on subglacial hydrology and ice sheet dynamics (see, e.g., Poinar et al., 2019).

Conclusions
Previous work on Greenland's surface hydrology has assumed that water flow through crevasses, if considered at all, is spatially homogeneous and can be predicted using simple physical thresholds. Our analysis of regional satellite data and local UAV surveys has demonstrated that crevasses instead exhibit spatially variable but inter-annually persistent hydrology across a ∼2,700 km 2 marine-terminating sector of the western GrIS. Only 26% of crevasses are observed to pond, which we infer to result from the hydraulic isolation of crevasse systems in areas of compressive mean surface-parallel stress. Through UAV surveys, these ponded crevasses were shown to drain rapidly and episodically, likely through full-depth hydrofracture. The remaining 74% of crevasses remain unponded throughout the observational period, which we infer to be due to water draining into the wider englacial and supraglacial system, connected through linked fracture networks that are actively maintained in extensional stress regimes. Our findings indicate that controls on crevasse ponding are distinct from better-studied processes, such as supraglacial lake formation and crevasse opening, that are often used to represent the mechanics of crevasse hydrology in surface routing models. This highlights the need for a better implementation of crevasse hydrology in surface routing models, particularly as early implementations show that crevasses can act as pathways for nearly half of all meltwater (Koziol et al., 2017). Our observations indicate that some form of simple stress threshold may still be suitable to drive parameterisations in regional-scale models of meltwater drainage. However, further observations are necessary to improve our process-based understanding of crevasse hydrology, including in situ observations of crevasse-scale mechanics and ice-sheet-scale satellite observations of spatio-temporal variability. Understanding the full spectrum of variability in crevasse hydrology is not only essential to be able to model the GrIS's thermodynamic response to increasing surface runoff in the 21st century, but has wider implications across other cryospheric contexts, such as ice shelf hydrology and breakup (e.g., Lai et al., 2020).