Contrasting regional variability of buried meltwater extent over two years across the Greenland Ice Sheet

. The Greenland Ice Sheet (GrIS) rapid mass loss is primarily driven by an increase in meltwater runoff, which highlights the importance of understanding the formation, evolution and impact of meltwater features on the ice sheet. Buried lakes are meltwater features that contain liquid water and exist under layers of snow, ﬁrn, and/or ice. These lakes are invisible in optical imagery, challenging the analysis of their evolution and implication for larger GrIS dynamics and mass change. Here, we present a method that uses a convolutional neural network, a deep learning method, to automatically detect buried lakes 5 across the GrIS. For the years 2018 and 2019 (which represent low and high melt years, respectively), we compare total areal extent of both buried and surface lakes across six regions, and use a regional climate model to explain the spatial and temporal differences. We ﬁnd that the total buried lake extent after the 2019 melt season is 56% larger than after the 2018 melt season across the entire ice sheet. Northern Greenland has the largest increase in buried lake extent after the 2019 melt season, which we attribute to late-summer surface melt and high autumn temperatures. We also provide evidence that different processes are 10 responsible for buried lake formation in different regions of the ice sheet. For example, in southwest Greenland, buried lakes often appear on the surface during the previous melt season, indicating that these meltwater features form when surface lakes partially freeze and become insulated as snowfall buries them. Conversely, in southeast Greenland, most buried lakes never appear on the surface, indicating that these features may form due to downward percolation of meltwater and/or subsurface penetration of shortwave radiation. We provide support for these processes via the use of a physics-based snow model. This 15 study provides additional perspective on the potential role of meltwater on GrIS dynamics and mass loss. during a relatively cold, low-melt year (2018) and a warmer, high-melt year (2019). We compare both buried lake distribution and surface lake


Introduction
The Greenland Ice Sheet (GrIS), which holds enough ice to raise sea level globally by more than 7 m (Church and Gregory, 2001;Smith et al., 2020) has experienced net mass loss every year since 1998 (Mouginot et al., 2019). Since 1972, the GrIS has contributed a total of 13.7 ± 1.1 mm of sea level rise. Prior to 2005, ice discharge was the primary driver of Greenland 20 mass loss (Enderlin and others, 2014). However, meltwater runoff, which has accelerated recently, is now the dominant factor in Greenland mass loss (Smith et al., 2020;Van Den Broeke et al., 2016;Enderlin and others, 2014). Thus, surface melt plays an increasingly important role in the mass balance of the GrIS.
Melting is exacerbated by the positive melt/albedo feedback, whereby melting acts to lower surface albedo, which in return allows for greater absorption of incoming solar radiation, thus further enhancing surface melt (Box et al., 2012;Lüthje et al., 25 2006; Tedesco et al., 2012). Meltwater often pools in surface lakes in the ablation zone during the summer from May to October (McMillan et al., 2007;Banwell et al., 2012). Water in these lakes then either runs off the ice sheet, drains via hydrofracture (Das et al., 2008;Tedesco et al., 2013;Williamson et al., 2018b), or refreezes in the firn (Bell et al., 2018). Firn takes up meltwater and buffers against mass loss (Harper et al., 2012). Water that refreezes in near-surface firn leads to the formation of ice lenses, which increases the density and decreases the porosity of the near-surface firn, thus reducing the capacity of 30 Greenland's firn to hold future meltwater . Ice lenses due to refrozen meltwater are rapidly increasing in areal extent across the GrIS, leading to increased runoff (MacFerrin et al., 2019;Culberg et al., 2021).
However, not all meltwater stored within the firn refreezes. In fact, some meltwater remains liquid, buried several meters below the surface (Koenig et al., 2015). These shallow buried lakes (or subsurface lakes) have been discovered later in the summer and at higher elevations on the GrIS than their surface counterparts (Miles et al., 2017;Lampkin et al., 2020). Liquid 35 water exists in buried lakes even throughout the winter (Schröder et al., 2020) and some have been shown to rapidly drain, delivering liquid water to the bedrock during the winter months (Benedek and Willis, 2021).
Because meltwater runoff is the dominant driver of mass loss from the GrIS, developing methods to detect surface water has been the focus of many recent studies (Yuan et al., 2020;Sundal and others, 2009;Williamson et al., 2017;Liang et al., 2012). Many of these methods, however, rely on optical imagery to detect meltwater. Buried lakes, in contrast, are invisible 40 in optical images, challenging the study of their extent, evolution, and interaction with drainage systems using conventional methods. More recently, studies have begun using Synthetic Aperture Radar (SAR) to detect both surface and buried lakes (Miles et al., 2017;Johansson and Brown, 2012;Dunmire et al., 2020;Schröder et al., 2020;Benedek and Willis, 2021). As an active sensor, SAR does not require a light source and thus is useful at night or during the polar winter. The European Space Agency's Sentinel-1 satellite uses C-band radiation, which is particularly useful because it is capable of penetrating clouds, 45 ice and snow up to several meters (Rignot et al., 2001). Water is a strong absorber of C-band SAR radiation. Because of this, Sentinel-1 microwave backscatter can be used as a diagnostic indicator of both the presence of melt and melt intensity, and can even detect liquid meltwater buried beneath several meters of ice and snow (Miles et al., 2017).
Fast, automatic detection of buried water features across large spatial and temporal scales is an important step in better understanding their impact on GrIS mass balance. In this study we develop a convolutional neural network (CNN), a deep 50 learning technique for automatic detection of features from images, to detect buried lakes across the GrIS. CNNs are beneficial for feature detection because of their ability to learn spatial relationships from two-dimensional images, and because they can easily be applied across large spatial and temporal scales. CNNs are becoming an increasingly popular choice for automatic feature detection in polar regions and have been used for detecting features such as glacier calving margins (Mohajerani et al., 2019;Zhang et al., 2019), icebergs (Rezvanbehbahani et al., 2019), sea ice concentration (Wang et al., 2016(Wang et al., , 2017Cooke 55 and Scott, 2019; Song et al., 2019), and surface water (Daneshgar et al., 2019;Yuan et al., 2020). Here, we develop and use a CNN to compare the distribution of buried lakes in Greenland's six subregions (Rignot and Mouginot, 2012) during a relatively cold, low-melt year (2018) and a warmer, high-melt year (2019). We compare both buried lake distribution and surface lake distribution between the two years, and subsequently use a regional climate model to help explain the spatial and temporal differences. Finally, by investigating optical imagery acquired prior to the date of the SAR imagery used here for buried-lake 60 detection, we hypothesize that different processes are responsible for the formation of these features in different regions of the GrIS.
2 Methods 2.1 Buried lake detection 2.1.1 Training/testing dataset preparation Our study region covers the entire GrIS. To create training and testing data sets we collected Sentinel-2 (S2) optical imagery and Sentinel-1 (S1) C-band SAR microwave backscatter imagery from all six GrIS subregions (Rignot and Mouginot, 2012). S2 images were collected during the late summer (September) of 2016 and 2017 and S1 images were collected during the winter between January 1 and January 7 following the 2018 and 2019 melt seasons. S1 images were chosen during the winter to minimize firn saturation than can obscure buried lakes in microwave imagery. Google Earth Engine (Gorelick et al., 2017) 70 (GEE) was used to collect all images used in this study. GEE preprocesses S1 images with the following steps: 1) thermal noise removal, 2) radiometric calibration, 3) terrain correction using ASTER DEM, and 4) values converted to decibels via log scaling. S2 images are available as top of atmosphere reflectance divided by 10000. All S1 and S2 images exported from GEE were rescaled to 30x30 m resolution (to increase processing speed) on the same grid using nearest neighbor resampling and reprojected to an Arctic stereographic grid (EPSG:3995). Images were broken up into 256x256 and 512x512 pixel tiles, 75 all resized to 256x256 pixels using a first order spline interpolation. For each tile, a false colour image was created using both the S1 and S2 imagery. The false colour images are made up of the following bands: -Band 1 = Greyscale (from S2) = 0.299 * Red + 0.5870 * Green + 0.1140 * Blue -Band 2 = Horizontally transmitted -vertically received (HV) band (from S1) normalized from 0 to 1 with respect to the individual image tile such that the minimum tile pixel is 0 and the maximum tile pixel is 1. The purpose of this band is 80 to emphasize microwave backscatter changes that result from the presence of liquid water.
-Band 3 = HV band (from S1) normalized between -30 and 0. This band functions to mute backscatter changes that are small with respect to the surrounding area, and thus not likely due to the presence of buried liquid water.
Image tiles were manually classified into seven different categories 1) uniform ice and snow (Fig. 1a, 1277 tiles total), 2) textured ice and snow (Fig. 1b, 1362 tiles total), 3) surface lake remnants (Fig. 1c, 1322 tiles total), 4) open ocean and sea water 85 (Fig. 1d, 99 tiles total), 5) buried lakes (Fig. 1e, 1399 tiles total), 6) land (Fig. 1f, 542 tiles total), and 7) mountains surrounded by snow (Fig. 1g, 227 tiles total). 70% of the image tiles for each class were used for model training, 15% for model validation during the training process, and 15% for a final model test. The training data set was augmented such that at each iteration of model training, image tiles were randomly flipped horizontally and/or rotated by a random increment of 30 • .

90
A CNN was the method of choice for detecting subsurface lakes across Greenland. Unsupervised learning has been shown to be unsuitable for similar analysis on Antarctica (Moussavi et al., 2020) and CNNs have proven superior over other supervised classification algorithms at detecting surface water on Greenland (Yuan et al., 2020). CNNs are being increasingly used for land and surface classification problems (Maggiori et al., 2017). Here, we use the pre-trained AlexNet architecture (Krizhevsky et al., 2012) because it outperformed other model architectures on our validation dataset (Appendix Table A1). To optimize the 95 model for our specific goal of detecting buried lakes, the fully connected layer was rebuilt using one which we designed. This layer used the Rectified Linear Unit (ReLU) activation function (Hara et al., 2015) to solve the vanishing gradient problem and implemented a dropout layer to prevent the model from overfitting (Srivastava et al., 2014). To optimize model parameters, we used the negative likelihood log loss function and Adam optimizer (Kingma and Ba, 2015). The learning rate began as 0.005, decreasing by a factor of 0.1 every 7 epochs using a scheduler. Our model was trained for 20 epochs with a batch size of 32 100 false colour training image tiles. Finally, our model was further optimized by feeding incorrectly classified tiles back into the training data set.

Model evaluation and testing
Two different methods were used to evaluate model performance: F1 score, and receiver operating characteristic (ROC) curves.  Our second model performance metric is a ROC curve. This metric compares true and false positive rates at different classification thresholds. For example, with a threshold of 0, any image tile that the model says has >0% chance of being a buried lake, is classified as such. Thus, at a threshold of 0, all image tiles will be classified as buried lakes and there will be high true and false positive rates. Similarly, for a threshold of 1, no image tiles will be classified as buried lakes and there will be low true and false positive rates. The area under the ROC curve (AUC) is used as an evaluation of model performance, with 115 larger AUC corresponding to better models.
Appendix Table A1 compares F1 score and AUC for all model architectures we trained to detect buried lakes. Our initial iteration of AlexNet had a relatively high recall and low precision. However, we decided that it was more important to mitigate false positives than to detect all possible subsurface lakes, thus prioritizing precision over recall. Therefore, to minimize false positives, we chose a relatively high threshold of 0.7 for classifying buried lakes. Using this threshold, our final model had a 120 buried lake precision of 0.929, recall of 0.880, and F1 score of 0.904. Precision for all classes can be seen in Appendix Figure   A1, which illustrates the confusion matrix for the test data set.

Classification
To detect buried lakes across the GrIS, we followed the workflow outlined in Figure 2. First, we collected S1 and S2 images across the entire periphery of the GrIS for the years 2018 and 2019, and followed the same data collection and pre-processing 125 steps as outlined above in Sect. 2.1.1. As described above, we used a threshold of 0.7 for buried lake classification such that if an image tile had > 70% probability of being a buried lake, according to the model, the tile was classified as such.
We also created a mountain/land mask to ensure that mountain areas were not erroneously classified as buried lakes. The mask was created by masking out areas where the grayscale band was <0.4 in image tiles that had >75% of being either a mountain or land tile, according to the model. The mask was dilated by 20 pixels to include regions of the ice that could have 130 been affected by mountain shadows.
We then used a combination of thresholding and morphological operations to outline individual buried lakes in the image tiles. First, we applied a threshold to band 1 of the false colour image (normalized HV backscatter from S1) such that all but the darkest 5% of pixels were masked out. We next performed a series of opening/closing operations using 3x3, 6x6, and finally 12x12 pixel kernels to fill holes in the lake and mask out darker pixels outside the lake. All binary image tiles were 135 then mosaicked and the individual buried lakes were contoured. We found that repeating this process of thresholding and using morphological operations improved the delineation of buried lakes so these steps were repeated for each contoured area.

False color image tiles (a)
-256x256 and 512x512 pixel tiles -Rescaled to 256x256 pixels -Band 1 = S2 grayscale -Band 2 = S1 HV tile-normalized -Band 3 = S1 HV normalized between 0 and 30 Finally, each contoured area had to meet a series of conditions to be considered a buried lake (Appendix Fig. A2). These conditions, which depend on the shape and size of the lake, were optimized by manual inspection to avoid false positives. For example, for smaller buried lakes to be considered, they had to have a lower (compared to larger lakes) microwave backscatter 140 compared to the surrounding background area.

Surface lake detection
We used all S2 images with <10% cloud cover, rescaled from the native resolution of 10x10 m to 30x30 m resolution (to match S1 image resolution for CNN work and for increased processing speed), from GEE to detect surface lakes across the GrIS during the 2018 and 2019 melt seasons. We followed Williamson et al. (2018a) to mask out clouds by removing pixels in 145 which the band 11 (SWIR) top-of-atmosphere (TOA) reflectance exceeded a threshold of 0.140. To mask ice-marginal areas, we followed Moussavi et al. (2020) by removing pixels with a Normalized Difference Snow Index (NDSI = Green−SW IR Green+SW IR ) < 0.85 and where band 2 (blue) < 0.4. Pixels above a Normalized Difference Water Index (NDWI = Blue−Red Blue+Red ) threshold of 0.5 (Miles et al., 2017) were considered to be a surface lake pixel. This threshold is higher than that used in other studies (Yang and Smith, 2013;Williamson et al., 2018a;Benedek and Willis, 2021), but we found that using a lower threshold resulted in the 150 erroneous inclusion of cloud and mountain shadows in the analysis. Additionally, use of this higher threshold was appropriate for our goal of comparing relative surface water ponding between the 2018 and 2019 melt seasons in each GrIS subregion. For lakes with small areas of floating surface ice, we filled in these areas by twice performing a morphological closing operation using a 5x5 pixel kernel. We only considered surface lakes that had an area > 0.05 km 2 , which is consistent with our area threshold for buried lake detection.

Buried lake analysis
We compared buried and surface lake elevations by calculating the mean elevation of each lake using the Greenland Ice Mapping Project (GIMP) elevation dataset . GIMP estimates of surface elevation have an ice-sheet-wide root-mean-square error relative to ICESat elevations of ± 10 m above sea level (asl) ranging from a minimum of ± 1 m asl in most ice areas to ± 30 m asl in high relief regions .

160
To investigate the buried lake formation process in different regions of the GrIS, we determined the percentage of detected buried lakes in each region that showed any evidence of surface meltwater during the previous summer. We define "evidence of surface meltwater" as the presence of > 5 pixels with an NDWI > 0.25 within the bounds of the buried lake. For several buried lakes, we also examined a combination of optical imagery from Sentinel-2, Landsat 8 and Landsat 7 from summers dating back to 2012. We also looked at surface topography profiles above select buried lakes using the Arctic DEM Mosaic at 165 2 m resolution (Porter et al., 2018).
Finally, we compared buried lake distribution with the locations of firn aquifers from 2015-2017 using the Miège (2018) firn aquifer dataset derived from Operation Ice Bridge radar.

Regional climate model analysis
To relate our buried lake detection results in the two consecutive years to the interannual variations in surface climate and 170 surface mass balance, we used output from the regional climate model RACMO2.3p2 ; RACMO2 hereafter).
RACMO2 provides output of Greenland climate and surface mass balance from 1958 to 2020 at 5.5 km horizontal resolution, which is then further statistically downscaled to 1 km horizontal resolution . We refer the reader to Noël et al. (2018) for a detailed description and evaluation of RACMO2 over the GrIS.
In this study, we analyzed RACMO2-derived monthly mean fields of 2m air temperature, surface melt, precipitation, and 175 surface mass balance for the GrIS from 1 June to 31 December in 2018 and 2019. We also compared these data to the long-term 1958-2017 climatology at elevations < 2500 m asl.

SNOWPACK modelling
To assess meltwater production, percolation and refreezing in the firn layer in areas with buried meltwater, we used the onedimensional, physics-based, multi-layer snow cover model SNOWPACK (Lehning et al., 2002b, a). We forced SNOWPACK 180 with RACMO2-derived precipitation, air temperature, relative humidity, incoming long and shortwave radiation, and wind speed for the closest grid point to a select buried lake in each of the CW, NW, and NE regions. Simulations were run from January 2000 to December 2019, starting with zero firn layer depth for sites with annual positive mass balance, and with 10 m deep ice for sites in the ablation zone.
Within SNOWPACK, the water flow in snow is solved using Richards Equation (Wever et al., 2014), which describes water 185 percolation based on capillarity and hydraulic conductivity. It enables the simulation of water accumulation on microstructural transitions in the snowpack. These transitions can be formed by contrasting microstructural properties (grain size and shape), as well as high variability in hydraulic conductivity (primarily governed by ice content) (Webb et al., 2018). We used the geometric mean to determine hydraulic conductivity at the interface nodes between two layers (Haverkamp and Vauclin, 1979). Using geometric mean maintains strong gradients in conductivity within SNOWPACK, allowing for water accumulation on capillary 190 and hydraulic barriers (Wever et al., 2018).
A zero flux bottom boundary condition was prescribed for temperature, and a free-drainage boundary condition was used for liquid water flow. Other model settings were chosen for optimal simulation of firn on ice sheets following Keenan et al. (2021). 195 We detect a total of 374 buried lakes (covering a total area of 241 km 2 ) following the 2018 melt season (Fig. 3a, Appendix Table B1) and a total of 599 buried lakes (covering 376 km 2 ) following the 2019 melt season (Fig. 3b, Appendix Table B1).

Surface and buried lake detection
In both years, buried lakes are concentrated along the western coast in the SW, CW, and NW regions of the GrIS. Detected buried lakes in SW Greenland have the largest mean area of all regions (0.77 km 2 following the 2018 melt season, and 0.90 ( Fig. B2). Buried lakes range in elevation from approximately 455 to 2450 m asl (Appendix Fig. B3), with the highest buried lakes located in the SW and SE (mean elevations of 1831 and 1710 m asl, respectively) and the lowest buried lakes located in the NE, NO, and NW (mean elevations of 1195 m asl, 1135 m asl, and 1144 m asl, respectively). These results are broadly consistent with Miles et al. (2017) and Koenig et al. (2015), who found that the majority of buried lakes ranged in elevation from 1000 to 2000 m asl. In general, buried lakes do not appear at substantially different elevations in 2019 than in 2018, with 205 the exception of the NO and NE regions. However, these exceptions may be the result of bias due to a relatively small number of buried lakes detected in 2018 in these regions (4 and 19 buried lakes in the NO and NE regions, respectively).
We also find that the total number of surface lakes and their total areal extent is much lower during the 2018 melt season than during the 2019 melt season. Using our method to detect surface lakes (NDWI > 0.5, and area > 0.05 km 2 ), we detect 3846 surface lakes (covering 1242 km 2 ) in 2018 (Appendix Fig. B1, Appendix Table B2) and 6146 surface lakes (covering 2569 210 km 2 ) in 2019 (Appendix Fig. B1, Appendix Table B2). Also similar to buried lake distribution, the mean elevation of detected surface lakes is highest in the SW and SE regions (mean 1365 m asl and 1335 m asl, respectively) and lowest in the NO, NW, and NE regions (mean 718 m asl, 860 m asl, and 962 m asl, respectively). In all regions for both years, the mean elevation of detected buried lakes is higher than for detected surface lakes (Appendix Fig. B3), which is consistent with previous work (Miles et al., 2017;Lampkin et al., 2020).

215
On an ice sheet wide scale, both the total surface and buried lake extents are much greater in 2019 than in 2018. Regionally, however, different patterns emerge between surface and buried lake distribution for the two years. Figure 4a shows that the ratio of surface lake area in 2018 to 2019 is much less than 1 across all GrIS subregions, ranging from 0.32 in NE Greenland to 0.67 in SE Greenland. The 2018:2019 buried lake area ratio, however, varies much more across the six subregions of the GrIS (Fig. 4a). In the SW and CW regions, we find approximately the same 2018:2019 buried lake area ratio (0.931 and 220 0.996, respectively), even though surface lake extent is much less in 2018. In contrast, in the NO and NE regions, the 2018 to 2019 buried lake area ratio is even smaller than the surface lake area ratio (  et al. (2020)). It makes sense that these features are co-located because they both require the similar climatic conditions of high 235 melt and high accumulation to exist Munneke et al., 2014;Dunmire et al., 2020). However, firn aquifers, unlike buried lakes, are buried too deep to be directly detected with S1 microwave imagery. The top surface of the perennial firn aquifer in SE Greenland is about 22±7 m below the ice surface (Miège et al., 2016) and can only be detected from S1 images by using the temporal change in microwave backscatter as a proxy to infer the locations of firn aquifers (Brangers et al., 2020). Although it is unclear what relationship the buried lakes and firn aquifers have, we suggest that the buried lakes may 240 feed firn aquifers by draining vertically.

Regional climate model and snow model analysis
To investigate the discrepancies between total surface and buried lake area across the six GrIS subregions, we analyzed RACMO2 data from 2018 and 2019 at elevations < 2500 m, comparing temperature and melt with the climatological mean (Appendix Fig. B4, B5). We find no spatial or temporal patterns associated with monthly precipitation or surface mass balance 245 anomalies that may explain surface and buried lake differences, so this analysis is not included here. Compared to 2018, 2019 is much warmer across elevations < 2500 m in every month from June to November (Fig. B4). During June and July 2019, high monthly mean temperatures are present across the entire ice sheet (+1.52 • C anomaly for June, +1.62 • C anomaly for July). It is likely that these anomalously high temperatures contribute to anomalously high melt during these months (Fig. B5, 1.51 and 1.02 standard deviations higher melt than the climatological mean for June and July, respectively). In contrast, in 2018, June 250 and July have mean temperature anomalies of −0.17 • C and −1.00 • C, respectively, and lower melt (0.24 and 0.33 standard deviations less than the climatological mean for June and July, respectively). Higher air temperatures in each region during June and July 2019 also contribute to higher ice-sheet-wide July 2019 subsurface temperatures (Fig. 7b). For sites X, Y, and Z respectively (Fig. 7a), the SNOWPACK-derived mean subsurface temperature in the top 7 m of the snow column is 2.06 • C,

Buried lake detection limitations
The number of buried lakes detected in this study is likely an underestimation of the actual number of buried lakes that existed across the GrIS in 2018 and 2019 for several reasons. Firstly, our detection of buried lakes was dependent on prescribed thresholds based on the microwave backscatter within the buried lake bounds relative to the surrounding area (Appendix Fig.   280 A2). We set these thresholds in a conservative manner that prioritized minimizing false positives, but as a result, our analysis likely missed some buried lakes that do not meet our thresholds. For example, if we increased the thresholds we use in this study by 5% (which makes the result less conservative), we saw 5.2% and 5.0% increases in the total detected buried lake volume in 2018 and 2019, respectively. However, this changed the 2018:2019 buried lake area ratio by only +0.18%. Conversely, a 5% decrease in the thresholds we use in this study (which makes the results even more conservative) resulted in 11.0% and 5.4% Another possible explanation is that buried lakes contained less water in 2018 than in 2019, making the lakes in 2018 appear 290 less different relative to their surrounding regions in the S1 images.
Another reason why our method may underestimate the number of buried lakes is because we do not consider lakes that are smaller than 0.05 km 2 (approximately 56 pixels), so our analysis may have missed smaller buried lakes. Additionally, as the penetration depth of C-band microwave radar in ice and snow is only several meters (Rignot et al., 2001), buried lakes located at depths greater than this will not be detected using our method. Finally, recent work has shown that buried lakes can drain, 295 and on very rare occasions, even during the winter months (Benedek and Willis, 2021); though this process will likely only contribute a very minor source of error.

Buried lake formation processes
Previous studies have generally assumed that buried lakes on the GrIS develop when a surface lake partially freezes through and then gets buried by snowfall, insulating any remaining liquid water under the surface (Koenig et al., 2015;Schröder et al., 300 2020). Modelling efforts have confirmed that this process is sufficient to allow liquid water to exist under the ice surface throughout the winter on both the Greenland and Antarctic ice sheets (Law et al., 2020;Dunmire et al., 2020;Lampkin et al., 2020). In SW Greenland, most of the buried lakes that we detected had some surface meltwater within their bounds during the previous melt season (Fig. 5), supporting the hypothesis that buried lakes form due to surface lake freeze-over followed by insulating snowfall, and indicating that this is the dominant process through which buried lakes form in SW Greenland. This However, in NW and SE Greenland, very few of the detected buried lakes had any evidence of previous summer surface melt within their bounds; suggesting that buried lakes in these two regions tend to form without prior surface ponding (Fig. 9). These results indicate that buried lakes in SE and NW Greenland may form via a different process than elsewhere on the GrIS. One 310 likely alternative explanation is that near-surface meltwater percolates downward and collects in depressions on impermeable ice layers, which are co-located with surface topological depressions (Drews et al., 2020). These subsurface ice layers could be the result of thermally controlled refreezing of downward percolating meltwater, providing impermeable layers for future meltwater to collect on top of. Evidence in support of this idea is that many of the buried lakes in the NW and SE are found to be located directly below surface depressions (Fig. 9e, j, k). Our SNOWPACK model simulations further support this idea. For 315 example, Figure 7c shows that meltwater accumulates on subsurface ice layers with low hydraulic conductivity. The simulated water content reaches up to 40-50% in some subsurface layers, indicating the formation of slush and subsequent pore space depletion. Substantial energy loss is required to refreeze slush layers, and autumn snowfall slows this refreezing process (e.g., compare autumn 2018 and 2019 in Fig. 7c). We also find that in autumn and winter, subsurface meltwater refreezing thickens the low permeability ice slabs, above which future meltwater may accumulate.

320
A second possible explanation for how buried lakes can form without initial surface meltwater is via penetration and absorption of shortwave solar radiation through the snowpack, which is capable of creating layers of slush and liquid water beneath the surface (Leppäranta et al., 2013;MacAyeal et al., 2019) without the presence of surface water. This is found to occur particularly when the surface is formed by pure ice. While downward percolation of meltwater is the more likely explanation for buried lake formation without the presence of surface meltwater in areas with a firn layer, subsurface penetration of solar 325 radiation may enhance melting in the vicinity of pre-existing buried lakes (Dunmire et al., 2020).
While Figure 5 indicates that the dominant buried lake formation process is different in different regions of the GrIS, it is likely that a combination of formation mechanisms exist in each region. Further, in some cases, it appears that a combination of formation processes even exist for individual lakes. For example, Figure 10 includes a series of images of a buried lake detected in CW Greenland. S1 and S2 imagery show that the areal extent of the buried lake is greater than the surface lake 330 detected during the previous melt season. These images therefore suggest that the formation of this buried lake resulted not only from the burial of a surface lake by snowfall, but also due subsurface melting and/or the downward percolation of surface meltwater.
We hypothesize that the dominant formation mechanism in different regions of the GrIS are driven by, in part, regional differences in annual precipitation. Generally, buried lakes that appear on the surface during the previous melt season are 335 located in areas with relatively lower total annual precipitation (Fig. 5b). Conversely, large concentrations of buried lakes that never appear on the surface during the previous melt season are located in areas with relatively higher total annual precipitation.
For example, in CW Greenland (Fig. 5c), the 1958-2017 climatological average annual precipitation that falls over the buried lakes detected in both 2018 and 2019 is 509 mm w.e/year for buried lakes which never appear on the surface during the previous melt season, and 451 mm w.e/year for buried lakes which do appear on the surface during the previous melt season. 340 The difference in these two means is statistically significant at the 99% confidence level. The observations described above therefore suggest that a mechanism by which annual precipitation may impact the buried lake formation process is through the availability of near-surface porous firn ( Figure B6). In areas with relatively low annual precipitation, there is less near-surface porous firn for meltwater to percolate, leading to increased surface ponding, and potentially a greater proportion of buried lakes that form following surface pond burial relative to other regions. In contrast, in 345 areas with relatively high precipitation, there is more near-surface firn available to store meltwater, which may initially pool on subsurface impermeable ice layers leading to less surface ponding prior to buried lake detection.

Climate drivers of buried lake distribution
Using a regional climate model to analyze temperature and melt anomalies from 2018 and 2019, we show that an increase in buried lake area in 2019 in northern Greenland (NW, NO and NE regions) is likely due to a combination of anomalously high 350 late-season melt and anomalously high near-surface autumn (i.e. August to November) air temperatures. Warmer late-summer and autumn air temperatures may prevent the water in surface lakes from freezing entirely and allow some liquid meltwater to remain buried. This analysis suggests that increasing air temperatures, which contribute to increased summer melt, will likely lead to an increase in the number and area of buried lakes, and therefore an increased volume of liquid water that is stored beneath the surface during the winter. Surface lake distribution across the GrIS has expanded inland to higher elevations 355 (Howat et al., 2013) and it is expected that surface lakes will continue to form at higher elevations as air temperatures continue to rise in the future (Leeson et al., 2015). Thus, we also expect that buried lakes will form further inland at higher elevations in a warming climate, with potential implications on changing surface hydrology and total ice loss.

Buried lake implications
Perennial buried lakes store meltwater that could otherwise runoff and eventually drain into the ocean. Thus, these features act to drain (Dunmire et al., 2020), even occasionally during the winter (Benedek and Willis, 2021). Because buried lakes exist perennially and can drain throughout the year, at lower elevations (< 1600 m, Poinar et al. (2015)), drainage of buried lakes could provide an influx of water to the bedrock at atypical times of year. However, as only a small fraction of buried lakes drain during the winter, the influx of meltwater to the bedrock during winter months is likely to be minimal. Additionally, as some 365 buried lakes in the SE and NW regions of the GrIS are located directly above observed firn aquifers, their potential drainage could provide an influx of meltwater to these firn aquifers. However, as the volume of water stored in buried lakes on an ice sheet wide scale is likely to be relatively small, these effects would mainly be relevant at local scales, particularly where large concentrations of buried lakes exist.

370
In this paper, we have developed a method that uses a convolutional neural network to detect buried lakes across the entire GrIS following each of the 2018 and 2019 melt seasons. We compared buried and surface lake spatial distributions across six GrIS subregions and used a regional climate model to explain the spatial and temporal differences. We find that while the total surface lake areal extent was less in 2018 than in 2019 across all six subregions, buried lake extent is roughly equal in 2018 and 2019 in the SW and CW regions, but much less in 2018 in the NO and NE regions. These regional differences in buried lake 375 extent between 2018 and 2019 can be explained by regional patterns of late-summer (August) melt and autumn (September -November) air temperatures. Anomalously high late-summer melt, coupled with anomalously high fall air temperatures in 2019 in northern Greenland likely contributed to increased buried lake extent in this region following the 2019 melt season.
Additionally, simulations using the one-dimensional SNOWPACK model confirm that subsurface meltwater in buried lakes remained unfrozen through autumn 2019. 380 We also examined the buried lake formation process by investigating S1 and S2 imagery prior to the detection of buried lake features. We find that different formation processes likely occur in different regions. In SW and CW Greenland, buried lakes likely form when surface lakes partially freeze over become insulated by following snowfall. In contrast, in SE and NW Greenland, no evidence of surface melt exists in optical imagery prior to the discovery of the majority of detected buried lakes; indicating that these features form via a different process. It is possible that subsurface penetration of shortwave radiation, 385 and/or downward percolation and subsurface pooling of near-surface melt, could be responsible for the formation of buried lakes in SE and NW Greenland. The simulation of the firn in NW Greenland showed liquid water accumulating on lowpermeability ice layers approximately 2 m below the surface.
Surface meltwater on the GrIS plays an important role in both ice sheet dynamics, as lake drainage events provide an influx of meltwater to the bed which may temporarily increase ice velocity, and direct mass loss, as meltwater runoff positively 390 contributes to sea level rise. The evolution of surface meltwater features has been thoroughly examined in prior studies using optical satellite imagery. In contrast, buried meltwater features are invisible in optical imagery, and as a result, these features are more poorly understood than their surface counterparts. Here, we provide a method for continent-wide mapping of buried lakes and a comprehensive look at the drivers of their formation and distribution across the GrIS.  Figure A1. CNN test data set confusion matrix. Values are normalized by the predicted conditions, and thus represent the precision for each class. For example, a value of 92.9% for predicted = buried lake, and actual = buried lake, means that 92.9% of the predicted buried lakes are actually buried lakes.

Buried lake contour
No Lake Is: BG diff < -0.2 Q4 < -0.17 Is: Area > 1.5 km 2 BG diff < -0.15 Q4 < -0.14 Lake No Lake Lake No Lake Lake Key Area: buried lake area (in km 2 ) Roundness: ratio of lake are to convex hull of lake area (1 = perfectly round) WL ratio: ratio of longest lake side to shortest lake side BG diff: average backscatter within lake bounds minus average backscatter outside lake bounds within 1 km Q4: average backscatter within lake bounds minus average backscatter outside lake bounds within 1 km on least-different lake side Figure A2. Thresholds used in the CNN for different buried lake shapes/sizes