The Impact of Bare Ice Duration and Geo‐Topographical Factors on the Darkening of the Greenland Ice Sheet

Dark (low albedo) surface ice on the Greenland Ice Sheet enhances melting and subsequent runoff, a major mass loss contributor during the ablation season. The accumulation of both biological (e.g., glacier ice algae) and abiotic (e.g., mineral dust) light‐absorbing particulates are important darkening factors, that are potentially influenced by the duration of snow‐free, bare ice (a phenological factor), and other geo‐topographical factors such as elevation, slope, aspect and the distance from the ice margin. Here, we present the first medium‐resolution (30 m) analysis of the phenological and geo‐topographical controls on the distribution of dark ice in SE and SW Greenland from statistical analysis of data derived from a harmonized satellite albedo product and ArcticDEM. The duration of bare ice primarily controls the distribution of dark surface ice, allowing for algae growth on inland ice surfaces in particular, whereas geo‐topographical factors are only secondary controls.

• The duration of snow-free, bare ice is the primary control on the distribution of dark ice • Longer durations of bare ice (median = 40 days) are prerequisites for ice to become dark • Surface slope and aspect influence darkening by modulating the duration of bare ice and glacier ice algae growth

Supporting Information:
Supporting Information may be found in the online version of this article.
and accumulation of meltwater (Greuell, 2000;Knap & Oerlemans, 1996) and light-absorbing particulates (LAPs) in the surface ice.These LAPs include dust derived from dry deposition of contemporary atmospheric aerosol, melt out of old dust from outcropping ice (McCutcheon et al., 2021;Wientjes & Oerlemans, 2010;Wientjes et al., 2011Wientjes et al., , 2012)), black carbon, and glacier ice algae (Anesio et al., 2017;Cook et al., 2017Cook et al., , 2020;;Lutz et al., 2018;Onuma et al., 2022;Perini et al., 2019;Ryan et al., 2018;Stibal et al., 2017;Williamson et al., 2018;Yallop et al., 2012), and cryoconites (Cook et al., 2018;Takeuchi et al., 2018).The accumulation of LAPs and the bio-albedo feedback amplify ablation and explain more than 70% of the spatial variability of albedo at regional scales (Ryan et al., 2018;Williamson et al., 2018).These studies focus mainly on the DZ, but dark ice is also present in other ablation zones around the GrIS, albeit in more restricted areas (Bøggild et al., 2010).Geo-topographical factors, such as slope and distance from the ice margin, have long been considered potent factors that influence the distribution of dark ice and melt.Terrain also has a strong influence on regulating the incoming solar radiation (J.He et al., 2019), the albedo (Colombo et al., 2023;Picard et al., 2020), and local wind-driven snow transport and melt (Cao & Barros, 2023).Field observations and remote sensing analysis show that sloping, south-facing snow/ice surfaces receive a higher amount of incoming solar radiation (Adhikary et al., 2000;Colombo et al., 2023;J. He et al., 2019) north of the Tropic of Cancer, which also increases surface melt.Finally, the GrIS is fairly flat inland, with steeper slopes near the ice margin, which helps promote katabatic and barrier winds, enhancing surface melt (Wang et al., 2021).The flatter surface ice at higher elevations enhances the accumulation of LAPs and meltwater (Greuell, 2000;Knap & Oerlemans, 1996;Wientjes & Oerlemans, 2010).
The duration of bare, snow-free surface ice is believed to be a key factor influencing the biological darkening through algae growth in summer (Tedstone et al., 2017), and is hereafter referred to as a "phenological" factor.This is strongly correlated with the mean air temperature in July (Shimada et al., 2016) since the average surface albedo is closely associated with the migration of the snowline to higher elevations, which exposes more bare surface ice (Ryan et al., 2019).Here, we aim to investigate whether the regional and local distribution of dark ice in the southeast and the southwest of GrIS can be explained by this simple phenological and other geo-topographical factors (elevation, slope, aspect, and distance from the ice margin) using a harmonized satellite albedo (HSA) product (Feng et al., 2023) and ArcticDEM (Porter et al., 2018).

Data Preprocessing
This study focuses on the southwest (SW) and southeast (SE) regions of the GrIS (Mouginot et al., 2017(Mouginot et al., , 2019;;Noël et al., 2019) (Figure 1a).The SW is less cloudy than other Arctic areas because of the persistent high atmospheric pressure (Noël et al., 2019).The increased frequency of anticyclones in central GrIS during the past decades induces more frequent advection of warm southerly winds over the western flank of the GrIS and enhances ablation (Fettweis et al., 2013;Hanna et al., 2014).Surface slopes are also shallower than in the SE (Figure 1).Consequently, the ablation zone in the SW is about 30% larger (Figures 1a and 1c).
ArcticDEM (Porter et al., 2018) provides digital elevation models (DEMs) from 16 August 2009 to 12 March 2017 with 2-m resolution, from which elevation, slope, and aspect were derived.Distance from the ice margin is the geodesic distance from a point on the ice surface to the ice sheet margin, defined by the shapefile of the study area (Mouginot & Rignot, 2019).The HSA (α HSA ) is a 30-m resolution albedo product derived from harmonized Landsat and Sentinel 2 data sets (Feng & Cook, 2023;Feng et al., 2023) developed in Google Earth Engine (GEE) (Gorelick et al., 2016).Sentinel 2 surface reflectance data is unavailable on GEE for the study period (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017), so the albedo is mainly derived from Landsat.
Both the ArcticDEMs and the HSA were masked by the GrIS basin data set from Mouginot and Rignot (2019) to remove non-ice areas.Ice surfaces with elevations above 2,000 m a.s.l (green line in Figure 1a) were masked out from our analysis as they are predominantly snow-covered (Fausto & the PROMICE team, 2018).The duration of snow-free, bare ice, hereafter the phenological factor, is derived from the HSA.Typically, ice albedo ranges from 0.20 to 0.65 (e.g., dark ice <0.45, 0.45 < clean ice < 0.55, 0.55 < superimposed ice < 0.65) (Cuffey & Paterson, 2010); and albedo values >0.65 were excluded to limit the analysis to (snow-free) ice surfaces only (Figure 1a).The pixel-wise snow-free duration at the peak of the melt season (July-August), when the annual spatial extent of the DZ is maximum (Tedstone et al., 2017), was calculated from the dates of the first and last day of albedo <0.65.Landsat has a revisiting time of 16 days, and the temporal resolution is 8 days by combining two Landsat satellites (Zhu et al., 2012).Limiting the study to July-August also helps to reduce uncertainties in estimating the phenological factor.

Random Sampling With Buffer and Filtering
The pre-processed albedo and ArcticDEM were matched by the timestamp (dt < 3 days), and then both data sets were sampled for analysis.The study area was divided into grid cells of 500 m (Figure 1b).Random sampling points (189,857 in the SW and 98,287 in the SE) were generated inside the grid cells with a strict buffer of 500 m to avoid oversampling and spatial autocorrelation, a known issue in geospatial data analysis, which leads to biased results (Gorelick, 2021;Ploton et al., 2020).Sampling points generated at the regional boundary were excluded.
Geo-topographical and phenological data were extracted using the generated random sample points at 30-m scale.

Regional Comparison and Statistical Analysis
The regional analysis tests the hypothesis that geo-topographical and phenological factors affect the spatial extent of dark ice by analyzing the statistical relationship between albedo and each factor.The Random Forest (RF) classifier (Breiman, 2001) has become increasingly popular in remote sensing because of both its reliability in classification and its capacity to accurately evaluate and rank importance of the predictor (Belgiu & Drăguţ, 2016).The averaged annual elevation, slope, aspect, distance, and duration data sets were normalized from 0 to 1 and split into training (70%) and test (30%) data sets for both regions.An RF model was trained to classify the ice types in the SE and SW using the training data set.The relative importance of the factors on ice types in both regions was analyzed by comparing the estimated out-of-bag permuted predictor importance (Breiman et al., 1984;Loh, 2002;Loh & Shih, 1997).The classification accuracy was evaluated using the test data set.A consequence of the random seed selections is that the final results differ even with the same input parameters.Hence, the RF model was run 10 times with the same input, and the resulting predictor importance values were averaged.The model uncertainty range is estimated using one standard deviation around the mean of the 10 runs.

Subregional Analysis in the Southwest GrIS
Dark ice in the SW extends further inland than in the SE (Figure 1a), enabling division of the sampled data into two groups determined by the median distance (Figures 2j and 2l) to the margin, namely margin (distance ≤9.57km) and inland (distance >9.57km).The number of sampling points for dark and bare ice in the subregional analysis are: (a) bare ice: 23,037 margin and 41,359 inland; (b) dark ice: 19,303 margin and 19,301 inland.
Comparison of ice at the margin and inland tests the hypothesis that topographic factors enable the dark ice to extend further inland in the SW than in the SE.

Regional Comparison and Analysis
Dark ice is more prevalent in the SW, with 37.5% of sampled ice pixels being dark ice, compared to only 25.2% in the SE (Figures 2a and 2b).The majority of the SW dark ice (interquartile range, IQR) falls in the elevation range of 673-1,214 m a.s.l (Figure 2c), and the median elevation is 1,013 m a.s.l.This is 500 m higher than the median elevation in the SE (612 m a.s.l.) and >100 m higher than the SE upper quartile (907 m a.s.l).The bare ice median elevation is higher than that of the dark ice in both regions.Generally, slopes in the SE are steeper than in the SW, and the slopes containing dark ice are steeper than those of bare ice (Figure 2f, right-tailed Wilcoxon rank sum test: p < 0.001).The distribution of aspects in the SE is smaller than in the SW (Figure 2i, left-tailed Wilcoxon rank sum test: p < 0.001).The median aspect of dark ice in the SE is 176°, quite close to the geographic south (180°), while the median aspect in the SW is 217°.The distribution of distance from the margin varies greatly in both regions.Half of the dark ice samples in the SE are limited to areas <2.7 km away from the ice margin, and most bare ice surfaces are found no further than 15.3 km (75th percentile) inland (Figures 1a, 2k, and 2l).Dark ice in the SW is also found near the margin (Figure 2j), but extends much further inland (upper whisker >50 km in Figure 2l).
The snow-free duration of bare ice in the SW is statistically longer than in the SE (Figure 2o, right-tailed Wilcoxon rank sum test: p < 0.001), although both the median values for dark ice are 40 days (Figure 2o).The estimated predictor importance of the RF model is shown in Figures 3a-3d.The overall classification accuracies for the SW and SE are 77% and 78%, respectively.The snow-free duration is the most important predictor in both regions, contributing 31 ± 5% and 33 ± 3% to the classification of dark and bare ice in the SW and SE, respectively.Elevation is the second most important predictor, with 29 ± 4% (SW) and 18 ± 3% (SE), followed by the distance from the ice margin scores (24 ± 4% (SW) and 22 ± 2% (SE)) and slope (6 ± 1% (SW) and 15 ± 2% (SE)).

Ice at the SW Margin and Inland
The median dark ice elevation is lower than bare ice for both the ice at the margin and inland (Figure 4a).Slopes of dark ice are steeper than bare ice at the margin (Figure 4b, right-tailed Wilcoxon rank sum test: p < 0.001).
However, the inland ice surfaces are mostly flat (IQR between 1 and 2°), regardless of the ice class.The distribution of aspect at the margin and inland differ (Figure 4c).Dark ice aspects are statistically smaller than for bare ice at the margin but are greater inland (one-sided Wilcoxon rank sum test: p < 0.001).The median of the aspect of inland dark ice is 235°, which is closer to the geographic west than the dark ice near the margin (median = 193°).
The snow-free duration of dark ice near the margin is typically longer than the inland ice (Figure 4d).

Geo-Topographical and Phenological Controls on Dark Ice in the SW and SE
The snow-free duration is the most critical factor (Figure 3) in the occurrence of dark ice in both regions (median = 40 days, Figures 2m-2o).Longer durations give the biological component of LAPs, glacial ice algae, more time to grow.Meltwater and the abiotic components of LAPs (dust, cryoconite, and black carbon) may also accumulate as the melt season progresses (Lloyd, 2021;Williamson et al., 2019Williamson et al., , 2020)).We hypothesize that the presence of algae helps to aggregate and retain abiotic impurities due to their production of cohesive compounds, such as extracellular polymeric substances (Holland et al., 2019;Stibal et al., 2012), which in turn reduce surface albedo and generate more melting.Snowlines are projected to retreat into the flatter ice sheet interior (Figures 1a  and 1c) under warmer climates, unless offset by enhanced winter accumulation (Ryan et al., 2019).Our analysis also suggests that future expansion of the bare ice area and prolonged snow-free duration may induce a rapid expansion of dark ice distribution, causing positive feedback resulting in even more surface melt.
The spatial distribution of dark ice is limited to the ice margin in the SE, while the dark ice can extend up to 50 km inland in the SW (Figures 1a and 2j-2l).The areas below the average snowline (Ryan et al., 2019) make up about 23.6% and 21.8% of the total area of the SW and SE, respectively (Figure 1c).However, the dark ice to bare ice ratio in these regions is disproportionally different, with 37.5% in the SW and 25.2% in the SE.The snow-free duration of bare ice in the SW is statistically longer than in the SE (Figure 2o), making it possible for microbial communities to colonize more ice surface.
Geo-topographical factors play a secondary role in regulating the distribution of dark ice, although elevation, surface slope and aspect may influence the snow-free duration.Generally, the snow-free duration decreases with elevation, but we note that dark ice in the SW occurs at higher elevations than in the SE.South-facing (median = 179°) slopes predominate in the dark ice of the SE.The influence of LAPs is more crucial on south-facing slopes than on flatter surfaces because the ice receives more intense solar radiation (Adhikary et al., 2000;Colombo et al., 2023); therefore, the snow-free duration is prolonged.

Geo-Topographical and Phenological Controls on Dark Ice at the Margin and Inland in SW
The bare ice and dark ice of the ice margin and inland in the SW have different relations with topographic factors (Figure 4).The slopes of dark ice are statistically steeper at the margin than inland (Figures 2d-2f and 4b).Steeper slopes near the ice margin expose bare ice at the onset of spring if winter accumulation is small (van den Broeke et al., 2009); therefore, the snow-free duration is longer (median duration = 47 days, Figure 4d).Melting downslope can be further enhanced by katabatic winds (Wang et al., 2021).However, the influence of slope is small across inland GrIS, since surface slopes are generally low (over 75% of the interior area <2°).
Dark ice is more dependent on the snow-free duration inland than at the margin, where snow-free periods are longer (Figure 4d).Aspect is a key parameter that adjusts the solar irradiance (Gutiérrez-Jurado & Vivoni, 2013;Rogeau, 2018), and dark ice aspects at the margin are closer to geographic south than bare ice (Figure 4c), as is also found in the SE (Figures 2h and 2i).However, inland dark ice favors aspects (median = 235°) that deviate further from the geographic south.This is because the darkening of the DZ is strongly associated with the growth of glacier ice algae (Cook et al., 2017(Cook et al., , 2020)), which can be suppressed under high light intensity (Jensen et al., 2023).The quantum efficiency of glacier algal photosystem II (PSII) is depressed under higher irradiances found on the GrIS (Williamson et al., 2020), and it could be that southwesterly aspects favor glacier algal growth over direct south-facing areas that receive the greatest irradiance.The more westerly facing flat ice surfaces in the SW may promote glacier algae colonization of areas with less intensive radiation.

Limitations
The definition of dark ice relies on an arbitrary albedo threshold (α < 0.45) and the duration that it remains dark.Future work should quantify the contributions of biological and abiotic darkening factors to dark ice formation.The snow-free duration was determined from the HSA.The accuracy is subject to the temporal resolution of clear (cloud-and saturation-free) observations (Zhang et al., 2021).The revisiting time of two Landsats is 8 days at best; therefore, it may not capture the exact timing of the onset and end of the melt season.The impact of this uncertainty is assumed to be limited as we restricted the analysis to the peak of the melt season (July-August).Saturated snow has an albedo <0.65 (Cuffey & Paterson, 2010), which may increase noise in the analysis.However, this higher threshold helps to capture the duration of melt season using low temporal resolution satellite data.A sensitivity experiment was run by changing the dark ice upper albedo threshold from 0.65 to 0.55 and the lower from 0.45 to 0.40.These changes do not have a great impact on the main statistical results (Table S1 in Supporting Information S1) or the rankings of the relative predictor importance (Figures S1 and S2 in Supporting Information S1).The difference between the estimated predictor importance of elevation (16.2 ± 2.4) and duration (16.3 ± 2.9) after changing the bare ice threshold to 0.55 (Figure S1 in Supporting Information S1) is very minor compared to Figure 3. Hence, the choice of dark ice thresholds in the ranges we have chosen does not change our main conclusions.

Conclusion
The darkening of ice on the GrIS results from a complex combination of abiotic and biological processes.Geo-topographical factors (elevation, slope, aspect, and the distance from the ice margin) and the duration of snow-free bare ice, a phenological factor were analyzed at regional and subregional scales.Results revealed that the snow-free duration is the primary factor governing the darkening of surface ice, enabling LAPs and meltwater to accumulate.Climate forcings regulate the extent of snow-free ablation areas, making room for new bare ice to darken.Geo-topographical factors are secondary controls on the darkening processes at an annual resolution.Surface slope and aspect influence the darkening by modulating the snow-free duration in the SE and the margin of the SW, where the dark ice predominates on steeper slopes and south-facing aspects.We hypothesize that aspect influences the biological darkening of the inland flat surface ice, where southwesterly slopes may be preferentially colonized over southerly slopes because of their lower solar irradiation receipts.Future upward migration of snowlines and longer snow-free duration may trigger a rapid expansion of dark ice on the GrIS, causing positive feedback that amplifies meltwater runoff production.innovation program (Grant agreement No 856416).Kathrin Naegeli is supported by the ESA PRODEX Trishna T-SEC project (PEA C4000133711).DEMs were provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691, and 1542736.The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.We would like to express our thanks to an anonymous reviewer, Chris Williamson, and the editor, Mathieu Morlighem, for their insightful commentary, which greatly improved the quality of this manuscript.Ian Stevens made valuable suggestions during the revision process.

Figure 1 .
Figure 1.Map of the study area (a) and exemplary random sampling points with strict 500 m buffers (b).The map shows the average summer (July-August) snow-free ice albedo (α mean < 0.65) in the study period (2009-2017).Contour lines are generated from the ArcticDEM mosaic (Porter et al., 2018).Subfigure c displays the ice sheet area at different elevations in the Southeast and Southwest Greenland Ice Sheet.The dashed lines denote the average snowline (SW:1,550 m a.s.l,SE: 1,453 m a.s.l, adapted from Ryan et al. (2019)).

Figure 2 .
Figure 2. Associations between (a and b) albedo and elevation, (d and e) slope, (g and h) aspect, (j and k) distance from ice margin, and (m and n) snow-free duration.The blue dashed lines indicate the threshold (albedo = 0.45) separating bare and dark ice.The distributions of each factor are shown in the notched boxplots in the column to the right (c, f, i, l, o), respectively.Ice classes are indicated by colors in the boxplots.

Figure 3 .
Figure 3. Out-of-bag permuted predictor importance estimated from the Random Forest model for geo-topographical and phenological factors in the (a and b) SW and (c and d) SE.The errorbars show the standard deviation of the predictor importance.