Dark ice dynamics of the south-west Greenland Ice Sheet

. Runoff from the Greenland Ice Sheet (GrIS) has increased in recent years due largely to changes in atmospheric circulation and atmospheric warming. Albedo reductions resulting from these changes have ampliﬁed surface melting. Some of the largest declines in GrIS albedo have occurred in the ablation zone of the south-west sector and are associated with the development of “dark” ice surfaces. Field observations at local scales reveal that a variety of light-absorbing impurities (LAIs) can be present on the surface, ranging from inorganic particulates to cryoconite materials and ice algae. Meanwhile, satellite observations show that the areal extent of dark ice has varied signiﬁcantly between recent successive melt seasons. However, the processes that drive such large interannual variability in dark ice extent remain essentially unconstrained. At present we are therefore unable to project how the albedo of bare ice sectors of the GrIS will evolve in the future, causing uncertainty in the pro-jected sea level contribution from the GrIS over the coming decades.Herewe use MODIS

Abstract. Runoff from the Greenland Ice Sheet (GrIS) has increased in recent years due largely to changes in atmospheric circulation and atmospheric warming. Albedo reductions resulting from these changes have amplified surface melting. Some of the largest declines in GrIS albedo have occurred in the ablation zone of the south-west sector and are associated with the development of "dark" ice surfaces. Field observations at local scales reveal that a variety of lightabsorbing impurities (LAIs) can be present on the surface, ranging from inorganic particulates to cryoconite materials and ice algae. Meanwhile, satellite observations show that the areal extent of dark ice has varied significantly between recent successive melt seasons. However, the processes that drive such large interannual variability in dark ice extent remain essentially unconstrained. At present we are therefore unable to project how the albedo of bare ice sectors of the GrIS will evolve in the future, causing uncertainty in the projected sea level contribution from the GrIS over the coming decades.
Here we use MODIS satellite imagery to examine dark ice dynamics on the south-west GrIS each year from 2000 to 2016. We quantify dark ice in terms of its annual extent, duration, intensity and timing of first appearance. Not only does dark ice extent vary significantly between years but so too does its duration (from 0 to > 80 % of June-July-August, JJA), intensity and the timing of its first appearance. Comparison of dark ice dynamics with potential meteorological drivers from the regional climate model MAR reveals that the JJA sensible heat flux, the number of positive minimum-air-temperature days and the timing of bare ice appearance are significant interannual synoptic controls.
We use these findings to identify the surface processes which are most likely to explain recent dark ice dynamics. We suggest that whilst the spatial distribution of dark ice is best explained by outcropping of particulates from ablating ice, these particulates alone do not drive dark ice dynamics. Instead, they may enable the growth of pigmented ice algal assemblages which cause visible surface darkening, but only when the climatological prerequisites of liquid meltwater presence and sufficient photosynthetically active radiation fluxes are met. Further field studies are required to fully constrain the processes by which ice algae growth proceeds and the apparent dependency of algae growth on melt-out particulates.

Introduction
Overall mass losses from the Greenland Ice Sheet (GrIS) have increased substantially since the early 1990s (Rignot and Kanagaratnam, 2006;Rignot et al., 2011;Shepherd et al., 2012). The average rate of mass loss increased from 34 Gt yr −1 during 1992-2001 to 215 Gt yr −1 during 2002-2011 (Sasgen et al., 2012). During 1991-2015 the GrIS lost mass at a rate equivalent to approximately 0.47 ± 0.23 mm yr −1 of sea level rise, with a peak contribution in 2012 of 1.2 mm (van den Broeke et al., 2016). Increases in mass losses since 2009 have been dominated by Published by Copernicus Publications on behalf of the European Geosciences Union.
increased surface runoff, with only 32 % of the total loss in this period attributable to solid ice discharge (Enderlin et al., 2014). It is therefore essential to understand the processes which control surface melting in order to be able to quantify the contribution of the GrIS to sea level rise over the coming century.
Surface albedo plays an important role in modulating the surface melt caused by incoming shortwave radiation. A lower albedo permits more absorption of shortwave radiation, which in turn leads to enhanced ice melting, and so albedo is the dominant factor governing surface melt variability in the ablation area (Box et al., 2012). The effective albedo of the GrIS is controlled by external factors including solar zenith angle, atmospheric composition and cloud cover, as well as the inherent optical properties of the surface. For both snow-covered and bare ice surfaces these inherent optical properties are modified by (a) ice grain metamorphism, (b) meltwater on the surface or in interstitial pores and (c) light-absorbing impurities (LAIs), including biological and mineralogical substances (Gardner and Sharp, 2010), each of which generally lead to reduced albedo.
Declines in GrIS bare ice albedo have an immediate impact on runoff from the GrIS and hence the surface mass balance (SMB). Decreases in the SMB since 1991 are predominantly due to enhanced runoff from low-lying (< 2000 m a.s.l.) bare ice parts of the ice sheet (van den Broeke et al., 2016). Since around 2000 the surface albedo of several sectors of the GrIS has often been significantly lower each summer than was observed during the 1990s (He et al., 2013). GrIS summer albedo showed a negative trend during 2000-2012, with the largest decreases observed in western Greenland (Stroeve et al., 2013). Some of the decline in albedo can be attributed to increases in bare ice extent. The GrIS-wide bare ice ablation zone extent increased by 7158 km 2 per year on average from 2000 to 2014, although with substantial interannual variability of between 5 % (89 975 km 2 ) and 16 % (279 075 km 2 ) of the ice sheet surface (Shimada et al., 2016).
In the ablation zone of the south-west GrIS, albedo lowered by as much as 18 % from 2000 to 2011 (Box et al., 2012). The south-west has seen the greatest increase in bare ice extent, by on average 5.8 % per year, with a mean extent of 56 603 km 2 during 2000-2014 (Shimada et al., 2016). However, increasing bare ice extent alone is insufficient to explain the declining albedo. Remotely sensed optical imagery for this sector shows a band of relatively darker ice within the bare ice ablation zone which recurred annually in the same location over the period 2001-2007, beginning 20-30 km inland from the ice sheet margin and extending up to ∼ 50 km wide, which has been postulated to be caused by LAIs (Wientjes and Oerlemans, 2010). LAIs on snow/ice surfaces reduce reflectance the most in the visible part of the solar spectrum (Warren, 1984;Painter et al., 2001;Bøggild et al., 2010), and this effect enabled Shimada et al. (2016) to quantify the interannual extent of dark ice -both GrIS-wide and for the south-west sector -by applying an empirically derived reflectance threshold to the 620-670 nm band of MODIS satellite imagery acquired in July each year. They found that dark ice extent varied substantially between years, both GrIS-wide (from 3575 to 26 975 km 2 ) and in the southwest (from 575 to 15 025 km 2 ).
There are a range of possible causes of dark ice on the GrIS. One consists of outcropping particulates in ablating ice. Wientjes et al. (2012) acquired shallow ice cores from the south-west sector in which they found dust that they dated to the Late Holocene. They therefore suggested that the dust was deposited in the accumulation zone and flowed with the ice down to the ablation zone, causing darkening of the surface as the ancient ice melts. However, they were not able to measure absolute concentrations of dust in their ice cores to compare to non-dark regions of the ice sheet. Meanwhile, Shimada et al. (2016) found a statistically significant correlation (r = 0.69) between July dark ice extent and air temperature (and hence surface melt rates, potentially depositing more particles from the ancient ice on the surface) in the south-west sector, but did not identify the responsible component of the surface energy balance.
Another potential source of darkening is the deposition of black carbon and other inorganic impurities by wet and dry atmospheric deposition, which has been investigated in ice and snow elsewhere (Warren and Wiscombe, 1980;Warren, 1984;Warren and Wiscombe, 1985;Gardner and Sharp, 2010). However, black carbon appears unlikely to explain variations in dark ice on the south-west GrIS. First, concentrations of black carbon in snowpack in the north-western snow sector are too low to cause any appreciable darkening and have been stable or even slightly declining over the past decade (Polashenski et al., 2015). Second, fire events in North America and Eurasia became rarer from 2002 to 2012 (Tedesco et al., 2016). Third, there is no recent statistically significant trend in aerosol flux deposition estimates along the south-west margin of the ice sheet (Tedesco et al., 2016).
In addition to inorganic impurities alone, the ice sheet can be darkened by ice surface habitats. Cryoconite is an aggregate of inorganic materials bound together by extracellular polymers produced by microorganisms, predominantly cyanobacteria (Wharton et al., 1985;Takeuchi et al., 2001;Hodson et al., 2008;Cook et al., 2016). Cryoconite absorbs more shortwave radiation than the surrounding ice and so, when the surface energy balance is dominated by shortwave radiation, ice overlain by cryoconite will melt more quickly than the surrounding ice. This produces water-filled cryoconite holes with a floor of biologically active sediment (Gribbon, 1979;Cook et al., 2016). These holes range from a few centimetres to several metres in diameter and depth (MacDonell and Fitzsimons, 2008), can cover a large part of the ablation zone (Hodson et al., 2008) and have been observed to occur in the south-west region of the GrIS (Stibal et al., , 2015Cook et al., 2012;Chandler et al., 2015;Cameron et al., 2016). Hole formation increases the albedo relative to dispersed cryoconite by sequestering the lowalbedo cryoconite from the ice surface at depth, resulting in a hemispheric albedo increase that will be further enhanced by specular reflection when covered by a reflective layer of meltwater (Bøggild et al., 2010). Occasional stripping events associated with high ambient air temperatures have been observed in the McMurdo Dry Valleys, Antarctica, resulting in cryoconite hole melt-out, the redistribution of cryoconite aggregates over the ice surface and subsequent formation of new holes (MacDonell and Fitzsimons, 2008). Similarly, in the south-west ablation zone of the GrIS, Chandler et al. (2015) observed that in warm, cloudy conditions some cryoconite holes melted out and release debris, but with little corresponding reduction in cryoconite hole coverage and set against an overall increasing trend in cryoconite hole coverage as the 2015 melt season progressed.
Distinct from the assemblages of microorganisms associated with cryoconite holes, ice algae can bloom in the upper few centimetres of bare melting ice. Abundant assemblages of ice algal communities have been reported on bare ice in both west (Uetake et al., 2010;Yallop et al., 2012) and east Greenland (Lutz et al., 2014). Ice algae produce specialist pigments which absorb UV and visible wavelengths, protecting the photosynthetic apparatus from excessive radiation (Dieser et al., 2010;Yallop et al., 2012;Remias et al., 2012). These pigments may be a significant source of darkening to GrIS surface ice Lutz et al., 2014).
The influence of cryoconite, cryoconite hole processes and/or ice algal assemblages on the substantial interannual variability apparent in dark ice extent of the GrIS is currently unknown. Whilst Shimada et al. (2016) proposed cryoconite sequestration into cryoconite holes as the mechanism underlying the negative correlation (r = −0.52) between ice-sheetwide July dark ice extent and shortwave radiation, this relationship did not hold when examined for the south-west sector alone. Additionally, although opposed pyranometer measurements (300-1100 nm) demonstrated that local algal bloom patches on snow had lower albedo at these wavelengths than snow without visible blooms (Lutz et al., 2014), broadband albedo measurements relevant for energy balance have not been isolated from grain evolution, meltwater ponding and abiotic impurities.
In this study we aim to identify the "top-down" controls of significant variability in dark ice extent between successive melt seasons in the south-west of the GrIS. We first characterise the interannual dark ice dynamics of the southwest GrIS using visible satellite imagery to quantify dark ice in terms of its extent, duration, intensity and the timing of its appearance each year. We then examine the extent to which interannual variations in dark ice dynamics are controlled by prevailing seasonal meteorological and climatological conditions and how they could drive surface darkening through three potential processes: (1) inorganic particulate deposition or redistribution, (2) cryoconite hole processes and (3) growth of ice algal assemblages.

Identification of dark ice
We used the MOD09GA Daily Land Surface Reflectance Collection 6 product, which is derived from data acquired by the MODIS sensor on board NASA's Terra satellite, to map bare and dark ice. Collection 6 products include improved calibration algorithms to correct for MODIS sensor degradation on Terra (Lyapustin et al., 2014) which was responsible for an apparent decline in GrIS dry snow albedo over the last decade or so (Polashenski et al., 2015;Casey et al., 2017). We used the clouds discrimination layer of the MOD10A1 Daily Snow Albedo Collection 6 product to identify and discard pixels covered by cloud. Our full time series encompasses daily observations between May and September from 2000 to 2016 but here we concentrate mainly on observations made during JJA.
Both MODIS Level-2 products are delivered on a sinusoidal grid at 500 m nominal resolution which exhibits significant distortion over the GrIS and prevents simple comparison with meteorological fields output by the regional climate model Modèle Atmosphérique Régional (MAR) (Sect. 2.5). We therefore first re-projected the MODIS data to the polar stereographic projection using nearest-neighbour resampling, yielding a spatial resolution of ∼ 600 m × 600 m. When undertaking comparisons with MAR outputs, cloudfree MODIS data were binned into 7.5 km 2 pixels to match MAR's resolution.
We detected bare ice and then dark ice within bare ice areas by applying thresholds to reflectance values (R) (Shimada et al., 2016). For bare ice we adopted R 841−876 nm < 0.6. To detect dark bare ice we used R 620−670 nm < 0.45. Pilot field spectra acquired in July 2016 indicate that this slightly higher threshold -compared with R 620−670 nm < 0.4 used by Shimada et al. (2016) -captures dark ice more accurately (Appendix A). However, we note that we do not know precisely what this dark ice threshold represents physically. The red band (620-670 nm) sits within the visible wavelengths and is therefore affected mostly by LAIs rather than grain evolution or water ponding, which mostly affect the near-infrared wavelengths (700-1100 nm). However, we caveat that other mechanisms can also reduce the reflectance across the entire solar spectrum, including in the red waveband. These include reduction in volumescattering due to wind or water "polishing" the ice surface, infilling of interstitial air spaces with meltwater and "trapping" by roughness features such as crevasses. Nevertheless, by combining these thresholds we are able to distinguish at first order between clean and dark (LAI-laden) ice surfaces.

Selection of common area
We defined a common area of maximum dark ice extent. This enabled the spatial sampling area to be held constant when www.the-cryosphere.net/11/2491/2017/ The Cryosphere, 11, 2491-2506, 2017 calculating interannual statistics. We chose this approach over defining different areas of dark ice for each year because then the spatial sampling area would have changed dramatically from one year to the next, whereas we are mainly interested in the primary drivers of interannual variability in dark ice dynamics.
To define the common area, first, in each year, we identified the pixels which were flagged as dark on at least 10 days during June-July-August (JJA). Then, we retained only those pixels which went dark in at least 4 years of our time series. Finally, we removed all dark pixels which occurred within ∼ 1 km of the ice sheet margin as defined by the Greenland Ice Mapping Project (GIMP; Howat et al., 2014) in order to remove errant pixels consisting of mixed land and ice cover which remained after applying the GIMP ice area mask. The common area is depicted in Fig. 1 and covers ∼ 10 400 km 2 .

Metrics of dark ice dynamics
We derived four metrics to characterise spatiotemporal variations in dark ice. First, annual extent (D E ) corresponds to the extent (in km 2 ) covered by the pixels within the common area which were dark for at least 5 days in each year. Second, annual duration (D D ) was defined at each pixel in the com-mon area as the percentage of daily cloud-free observations made in each JJA period which were classified as dark and is thereby normalised for cloud cover. Third, intensity (D I ) was defined as the mean daily reflectance over 620-670 nm of all cloud-free pixels in the common area and annual intensity (D I ) as the JJA mean of D I . A lower value of D I or D I therefore means that dark ice intensity was greater. D I is on a continuous scale and so is independent of the stringent dark ice presence threshold defined in Sect. 2.1. Only days in which at least 50 % of the common area was cloud-free were included in the calculation, and only the cloud-free pixels within the common area were used. Fourth, normalised darkness (D N ) was expressed as and therefore provides a combined indicator of both the duration and intensity of dark ice presence. We note that cloud cover was present to some degree over the common area in almost every day of our time series, which prevented us from quantifying daily dark ice extent. At each pixel and for each year we identified the date on which (a) bare ice emerged from underneath the melted snowpack (t B ) and (b) dark ice appeared (t D ), if at all. In both cases we used a 7-day rolling window on the relevant time series of reflectance at each pixel. Each year, we identified the first rolling window at each pixel that contained at least 3 days of bare or dark ice (not necessarily consecutive) and 0 days of non-bare or non-dark ice, which therefore permitted up to 4 days of cloud cover in the window. We then selected the first day of bare or dark ice appearance from within the chosen window. This windowing strategy enabled us to minimise the likelihood of false-positive identification of bare and dark ice appearance dates which would have occurred if only looking at daily observations in isolation and also allowed us to ameliorate for cloud cover.
Finally, we calculated the median and interquartile range (25th and 75th percentile) of the day of year of bare ice appearance for the common area each year ( t B ) from the pixellevel data.

Meteorological and climatological data
We performed simulations of meteorological conditions over the GrIS using version 3.6.2 of MAR, a regional climate model (Fettweis et al., 2017). The model was run on an equal-area 7.5 km × 7.5 km resolution grid for the whole of Greenland and was forced at its boundaries every 6 h by ECMWF ERA-Interim reanalysis data (Dee et al., 2011). For comparison with dark ice dynamics we downsampled the MODIS-defined common area to MAR's resolution. We calculated mean shortwave-down (SW↓ ), longwavedown (LW↓ ) and sensible heat flux (SHF ) anomalies (positive when into the ice sheet surface) in the common area for each JJA relative to 1981-2000. We also calculated the mean daily snow depth in the common area from April to August each year, total snowfall (from t B to 31 August) and total rainfall (during JJA).
We characterised near-surface air temperatures in two ways. First, we defined the mean air temperature during JJA as T . Second, we defined the number of days in each JJA period in which the common area's daily minimum near-surface air temperature exceeded 0 • C as (T > 0).
As introduced previously, cryoconite holes form and tend to be sustained under SW↓ dominant conditions. Conversely, this suggests that they are likely to melt out, depositing cryoconite onto the ice surface, if the surface energy balance shifts to LW↓ or SHF dominant conditions. We therefore characterised the conditions which could cause melt-out of cryoconite holes as the "melt-out flux" (MOF) using where α was the daily mean MOD10A1 albedo over the common area (only on days with <50 % cloud cover) and LW↑ was 315.6 W m −2 for melting ice surfaces as defined by Cuffey and Paterson (2010). SHF is positive when into the ice sheet surface. MOF corresponds to the mean JJA MOF. We used the monthly Greenland Blocking Index (GBI) (Hanna et al., 2016) to consider the role of the synoptic atmospheric circulation in dark ice dynamics. The GBI is the mean 500 hPa geopotential height for the 60-80 • N, 20-80 • W region and therefore provides a measure of the extent of high-pressure blocking over Greenland. We calculated the mean GBI for each JJA period.

Characteristics of dark ice dynamics
Shimada et al. (2016) identified a general trend of increasing D E over time but also saw that D E on the south-west GrIS varied dramatically between years. We found similar characteristics in our expanded time series (Fig. 1). D E ranged from almost no dark ice identified (2000, 2001 and 2015) to wide, contiguous areas of dark ice stretching from 65.5 to 69 • N (2007, 2010, 2011, 2012, 2014 and 2016).
There was substantial interannual variability in D D during JJA. Generally, when D E was high, D D was also high, especially in 2010, 2012 and 2016. Moreover, the extension of our time series relative to Shimada et al. (2016) which only examined July to encompass June through August revealed relatively large D E and D D in 2005, 2007, 2008, 2009, 2011 and 2014. Examination of D I (Fig. 2) shows that most dark ice presence was concentrated into the months of July and August. In some years (2010 and 2016) more significant darkening of the ice sheet surface was observed as early as mid-June. In years when substantial darkening occurred there was a time lag between t B and the first widespread occurrence of dark ice of ∼ 10-15 days, although in some cases (e.g. 2010) this may be attributable to the large interquartile range in date of bare ice appearance (Fig. 2). D I tended to increase over the season. Variability in D I at daily to weekly timescales was minimal compared to the magnitude of variability over interannual timescales. Dark ice usually persisted until the onset of anti-cyclonic, cloudy conditions (Fig. 2, days shaded grey) and snowfall during late August and September, which buried the bare ice surface under snowpack for the winter period.
We did not find any evidence that the dynamics of dark ice in one year controlled dark ice dynamics the following year. There were years of higher D N recently (2012,2014,2016)  Moreover, D I values at end of one melt season were generally significantly different to those in the period after t B the following year (Fig. 2). We used t D to calculate cumulative D E in the common area through the summer (Fig. 3, red lines). In several years D E was very small (2000,2001,2003,2015). In years when medium D E occurred (e.g. 2005, 2006, 2008, 2009, 2013), dark ice appeared step-wise through July and into August. This step-wise appearance also occurred in the high D E years of 2007 and 2010. In contrast, the widespread expansion of D E in 2011, 2012, 2014 and 2016 occurred rapidly over just a few days in July. In particular, we found 28 large, single-day expansions in dark ice extent in our time series (defined as > 520 km 2 , equivalent to ∼ 5 % of the common area). These large changes cannot be explained by gaps in our time series owing to cloud cover: the median number of preceding days when cloud cover was > 50 % was 0, and the mean common area covered by cloud in the preceding 7 days was 34 %. There tended to be minimal further dark ice expansion in August. As shown by D E (Fig. 1) and D I (Fig. 2), the extent and intensity of dark ice tended to persist for the rest of the summer season (Figs. 1 and 2).

Controls on dark ice presence
In years when D E was higher then D D (Fig. 1) and t D (Fig. 3) tended to be spatially invariant across the common area. This suggests that the driver(s) of dark ice presence is (are) synoptic, governing dark ice dynamics over the whole common area. We therefore used meteorological and climatological variables representative of the common area to examine their potential impact upon dark ice dynamics.

Snow
Snow can control dark ice dynamics in at least two major ways: (a) the thickness of the snowpack from the preceding winter will, in combination with air temperatures, control t B ; and (b) snowfall that occurs during the melt season will at least temporarily obscure the bare ice surface. Figure 3 shows the mean snowpack depth and t B from April through to August each year. In the years of longest D D and greatest D I (2007,2010,2012,2016) bare ice appeared by roughly mid-June according to both MODIS and MAR, compared to other years when bare ice did not appear until early to mid-July.
Earlier t B is not just a function of total snowfall during the preceding winter but is also strongly dependent on the progression of melting during spring which, in extreme cases such as 2016, began as early as April, introducing liquid water to the snowpack and accelerating its warming despite additional snowfall in May. In contrast, in years such as 2015, significant melting did not occur for the first time until mid-June. Not only was 2014-2015 winter snowfall relatively large compared to other years in our study, but more snowfall occurred around the start of June just before the melt season started. Nevertheless, in general a thinner winter snowpack favoured earlier t B , and earlier t B in turn favoured increased D N (R 2 = 0.51, p < 0.01; Fig. 5f). Last, when further snowfall occurred during summer (Fig. 4b) then D N tended to be lower (R 2 = 0.36, p < 0.05).

Atmospheric energy fluxes
SW↓ was consistently positive from 2007 onwards but continued to show substantial interannual variability (Fig. 4a).
There was no statistically significant relationship between JJA SW↓ and D N . Unlike SW↓ , from 2000 to 2007 LW↓ anomalies were consistently positive and then after 2007 the sign became more variable, with both positive and negative anomalies occurring (Fig. 4a). Like SW↓ there was no statistically significant relationship with D N . SHF was consistently positive throughout the time series. There was a significant positive correlation between SHF and D N (R 2 = 0.41, p < 0.01; Fig. 5e). We examined the likelihood for cryoconite hole melt-out (causing redistribution of cryoconite materials onto the ice sheet surface) using MOF (Fig. 4c), which is suggestive of the energy balance conditions that are needed to melt cryoconite holes out of their weathering crust. Positive MOF signifies that longwave and sensible heat fluxes dominate the energy balance, which could cause spatially "even" surface melting as opposed to spatially heterogeneous melting permitted by stronger absorption of SW↓ where cryoconite material is present. MOF was negative in all years, and all days within 3σ of the mean were also negative, which suggests that the positive MOF conditions required for the widespread melt-out of cryoconite holes were seldom (if ever) met.
As liquid meltwater constitutes a prerequisite for algal growth, we assessed the likelihood of continuous liquid meltwater presence on the ice surface over each 24 h cycle using (T > 0). We found a positive correlation between (T > 0) and D I (R 2 = 0.37, p < 0.05; Fig. 5c) and to a lesser extent with D N (R 2 = 0.27, p < 0.05; Fig. 5g). Greater p < 0.01). Moreover, we found that single days of large dark ice area expansion were associated with a median of 3 days of continuous (24 h) melting, compared to 0 days for the rest of the time series. These sudden increases in dark ice extent were associated with higher absolute sensible heat fluxes, with a mean of 57 ± 24 W m −2 , equivalent to 96 % more sensible heat than the SHF mean in the period from the start of dark ice expansion until a maximum D E of 90 % of the common area. Over daily timescales, higher SHF was associated with warmer near-surface air temperatures (R 2 = 0.54, p < 0.01) but more strongly with higher near-surface wind speeds (R 2 = 0.67, p < 0.01). Days on which the minimum air temperature was greater than 0 • C had mean wind speeds of 6.5 ± 1.8 m s −1 ±1σ compared to 4.9 ± 1.3 m s −1 ±1σ on days when the minimum air temperature was 0 • C or less.
Last, we examined whether dark ice dynamics have any relationship with the GBI. We found a positive correlation between the JJA GBI and D N (R 2 = 0.46, p < 0.05).

Rainfall
Rainfall can occur in the ablation zone of the GrIS during summer. Limited observations from elsewhere in the cryosphere indicate that whilst the direct melt impact of rainfall upon melt rates is generally limited, rain can affect melt indirectly by increasing the liquid water content of the ice surface, reducing its albedo (Hock, 2005). Total JJA rainfall in this sector ranged from 30 mm w.e. to as much as 140 mm w.e. (Fig. 4b). However, there was no statistically significant relationship between JJA total rainfall and D N . Eyewitnesses on the ice sheet surface have observed that the generally widespread porous surface weathering crust, typically on the order of 20-30 cm deep, can be stripped down towards the underlying high-density bare ice dur-ing rainfall events (potentially dispersing cryoconite material from cryoconite holes). However, they also observed that it tends to thicken again within days due primarily to renewed subsurface melt by incoming shortwave radiation (Cook et al., 2017a). We therefore also examined the impact of each JJA rainfall event upon D I . We selected all rainfall (and snowfall) events of > 1 mm w.e. day −1 across our common area in our time series. Then, we calculated the change in D I using the closest observations immediately before and after the rainfall event. We found no systematic change in D I caused by rainfall events: in some cases D I increased while in other cases it decreased. This was the case whether or not mixed rainfall and snowfall events were excluded from analysis, although we note that MAR may not adequately discriminate between rainfall and snowfall over ice surfaces.

Discussion
At outlined in Sect. 1, a number of processes have been proposed to explain the dark ice dynamics on the south-west of the GrIS. Our characterisation of dark ice in terms of D E , D D , D I and t D , when combined with analysis of the prevailing meteorological conditions estimated by MAR, allows us to consider the extent to which each proposed process fits with our observations of dark ice dynamics.

Variability driven by inorganic particulate deposition or redistribution
There are two primary ways in which inorganic particulate matter can arrive on the ice-sheet surface: (1) by wet and dry atmospheric deposition and (2) deposition of material previously trapped in the ablating ice. Previous research indicates that there is no relationship between albedo reductions and the number of fires occurring over North America and Eurasia or with modelled atmospheric aerosol fluxes (Tedesco et al., 2016). Measurements of black carbon in snow on the present-day surface of the GrIS (as opposed to in ice cores) have been made in the north-west (Aoki et al., 2014;Polashenski et al., 2015) and high in the accumulation zone (Hegg et al., 2010;Doherty et al., 2013). However, at only a few ppb, these measurements of black carbon are insufficient to explain the substantial reduction in reflectance in the south-west (see Warren, 2013). Atmospheric deposition events would presumably have to occur in only years of high D N in order to explain the spatiotemporal patterns in dark ice that we observed. In dark years, D I increased gradually over the summer and so we postulate that deposition would need to occur over at least several days after t B . More problematic is that maximum D E is relatively invariant and begins ca. 20 km inland rather than at the margin, whereas atmospheric deposition would presumably occur over a more dispersed area.
The well-defined geometry of the dark-ice area between years lends support to the hypothesis that the dark-ice area is caused by surface deposition of particulates previously trapped in the ablating ice (Wientjes and Oerlemans, 2010;Wientjes et al., 2012). Warmer air temperatures in darker years (R 2 = 0.43, p < 0.01; Fig. 5d and h) support the idea that more material is deposited on the surface from melted ancient ice in these years, contributing to darkening and possibly acting as a positive feedback mechanism. However, our dynamics observations suggest that particulates from ancient ice are not the primary source of darkening. First, there is more variability in D E and D D than we would expect if total summer ablation alone determined darkening by controlling the corresponding quantity of particulates being released from ancient ice. D E was negligible in several years (Fig. 1), yet melting in this area occurred in all years (e.g. van As et al., 2016). Second, particulates are unlikely to be dispersed homogeneously through the ice column and so the concentration of particulates emerging at the surface will change non-linearly with respect to the icemelt rate. This could explain why D E is negligible in several high-melt years. However, the wavy patterns of surface darkening at decametre scales observed by Wientjes and Oerlemans (2010) are indicative of dispersion of previously welldefined particulate horizons by vertical shear due to ice flow, which suggests that particulates from ancient ice were deposited in each summer of our time series. Third and most critically, in order to explain non-dark years such as 2013 and 2015, the particulate material responsible for substantial darkening during 2010-2012 and 2014 would have to be evacuated from the ice-sheet surface at the end of summer to explain both dark ice that summer and the lack of dark ice the year after. More broadly, we did not find a year-on-year increase in D I that we would expect if particulates from ancient ice were accumulating at the ice-sheet surface (Fig. 2).
We also found that in years of high D N the onset of high D I was delayed by ∼ 10-15 days after t B (Fig. 2). This may be attributable to remaining snow patches or superimposed ice formation, although the prevalence of the latter process on the GrIS remains poorly understood (Larose et al., 2013;Chandler et al., 2015). Nevertheless, if particulate materials were still present on the surface from the previous melt season we would expect high D I almost immediately after bare ice appearance.
Overall, our observations provide little support to the hypothesis that particulate deposition causes surface darkening. We are unable to identify a mechanism by which the icesurface mass flux of particulate material could change over timescales commensurate with dark ice dynamics. There are several noteworthy limitations to investigating cryoconite hole processes from satellite observations. First, cryoconite hole processes occur over decimetre scales and so we may not be able to capture their variability using 500 m MODIS imagery. Second, the reflectance of ice surfaces with cryoconite holes varies strongly as a function of viewing angle (Bøggild et al., 2010), and so observations made by MODIS -which has a "push-broom" scanning assemblywill vary depending on how near to nadir the angular instantaneous field of view (IFOV) is. This is likely to impact the broadband albedo value from MOD10A1 used to calculate SW net as part of MOF in a way that is not currently known. Shimada et al. (2016) suggested that low D E in July 2011 followed by widespread D E in July 2012 could be attributable to cryoconite hole wash-out during anti-cyclonic conditions in late August and September 2011. However, our results reveal a different spatiotemporal pattern of darkening: we found that the common area did go dark in 2011 but that this did not begin until late in July. Maximum D I was reached during August, which Shimada et al. (2016) omitted from their analysis. In 2012 we observed an early onset of high D E , with rising D I as the season continued, whereas under a return to SW↓ dominant conditions we would expect D I to gradually decrease as cryoconite hole surfaces deepen, albeit non-linearly because cryoconite holes cease deepening once they are in equilibrium with their surroundings (Gribbon, 1979). This also makes it difficult to explain low D E in 2013, as cryoconite holes would have needed to form over a short period at the end of summer 2012 in order to sequester cryoconite particles at depth, unless the presence of snow patches and/or superimposed ice at the surface was so prolonged that only in a few pixels did enough melting take place to expose bare/dark ice.

Variability driven by cryoconite hole processes
We also looked at intra-annual variations in dark ice dynamics when considering the potential for variability to be driven by cryoconite hole processes. Field observations of cryoconite hole morphology show that cryoconite holes form within the ice weathering crust over timescales of a few days (Cook et al., 2016). However, in dark years, t D was relatively synchronous across the common area. Importantly, D I then increased as the season continued, suggesting that episodic cryoconite hole flushing and reformation is unlikely. This is supported by field observations made in the south-west ablation zone of the GrIS by Chandler et al. (2015), who found that cryoconite hole coverage increased over the course of the 2015 melt season despite transiently warm, cloudy conditions experienced during July that caused some holes to melt out and release their debris. Moreover, the lack of energy available for cryoconite hole melt-out over seasonal timescales (Fig. 4c) suggests that variability in the areal extent of cryoconite holes forced by changes in the dominant component(s) of the surface energy balance cannot explain dark ice dynamics. Our only evidence in support of cryoconite hole processes is that large single-day increases in dark ice extent were associated with higher absolute SHF, which could cause transient hole-flushing events during the melt season. However, the rest of our evidence strongly suggests that cryoconite hole processes are not responsible for interannual dark ice dynamics.

Variability driven by ice algal assemblages
Last, we examined evidence for the role of ice algal assemblages as the principal driver of dark ice dynamics. In addition to typical light-harvesting pigments characteristic of green microalgae (Remias et al., 2009), ice algae produce a unique UV-VIS absorbing purpurogallin pigment that presumably affords protection from the significant radiation experienced in GrIS surface habitats (Remias et al., 2012). Given this pigmentation, ice algal blooms are known to impact visible reflectance at local (metre) scales Lutz et al., 2014). Knowledge of the regulation of temporal and spatial patterns in ice algal biomass (and thus pigmentation) in surface habitats is limited Chandler et al., 2015), but the fundamental prerequisites for algal life are known, including liquid water, nutrient resources and photosynthetically active radiation (PAR; 400-700 nm).
The significant positive relationship identified between (T > 0) and D I (R 2 = 0.37, p < 0.01; Fig. 5g) supports the role of ice algae in ice sheet darkening, as do the singleday increases in D E of > 5 % of the common area, which were generally preceded by several days of continuous positive air temperatures; both of these observations are indicative of liquid meltwater presence. Ice algae require liquid water in order to grow, and ice surfaces are reservoirs of potentially viable propagules that can become active when they encounter sufficient liquid water of appropriate chem-istry (Webster-Brown et al., 2015). Minimum air temperatures above 0 • C required for the presence of liquid water will facilitate growth of ice algae. As blooming progresses, the relationship between liquid water availability and algal proliferation may be strengthened by the establishment of a positive feedback loop via albedo reduction. For example, blooms of snow algae have been shown to result in surface albedo reduction, increased heat retention at the snow surface and thus enhanced melting and liquid water availability for continued algal growth (Lutz et al., 2016). Climatically, enhanced liquid meltwater presence in dark years -especially continuing through the night when SW↓ tends to zero -is also partially attributable to increased SHF , with a positive correlation observed between SHF and D N (R 2 = 0.41, p < 0.01; Fig. 5e) in our study. Thus a combination of greater (T > 0) and higher SHF may regulate interannual liquid water availability in ways critical to ice algae growth, thereby governing whether or not dark ice appears.
t B has a first-order impact on whether the common area is dark in any given year, with later appearance associated with lower D E (Fig. 5b and f). t B will significantly impact PAR availability at the ice surface. If bare ice appears in early to mid-June, it will receive PAR over several complete diurnal cycles, unlike in years when bare ice does not appear until July. Although ice algae likely experience excessive irradiance over the ablation season, as evidenced by their production of "sunscreen" pigments (Remias et al., 2012), a minimum threshold of PAR (or photo-period duration) may be required to allow bloom initiation, which would be favoured by earlier t B . Alternatively, variability in t B may impact algal blooms (and thus darkening) via the timing of nutrient inputs to surface ice or due to the formation of superimposed ice. With delayed snow line retreat, percolating snow melt in spring/early summer may release snowpack nutrients to surface ice (Larose et al., 2013) before PAR is available to allow algal utilisation, stalling bloom formation. It may also result in sustained presence of superimposed ice (Larose et al., 2013;Chandler et al., 2015), preventing PAR penetration to the previous year's ice algal cells and initiation of growth. However, we found no significant relationship between SW (which corresponds approximately to PAR) and D N , so the role of seasonal PAR fluxes in algal growth remains unclear. More broadly, field studies are required in order to identify precisely the transition from snow to a bare ice surface could impact ice algal assemblages.
If prerequisites for the initiation of an algal bloom are achieved then an increase in algal biomass is likely, with a concomitant increase in D I . This is consistent with increases in D I after the first appearance of dark ice, as opposed to "flickering" between less-and more-dark states. Increases in D I could be driven by an increase in the spatial extent of ice algal assemblages and/or an increase in algal concentrations per unit area, although we note that increases in D I which occur immediately after snow retreat could also be due to the melting away of superimposed ice (e.g. Larose et al., 2013).
Previously, Chandler et al. (2015) recorded an increase in "dirty ice" extent within our common area over an ablation period, though they did not assess algal cell numbers within dirty ice. Whilst algal concentrations likely increase until a limiting factor becomes apparent, analogous to algal blooms in aquatic systems (Teeling et al., 2016), progressive colonisation of clean ice at the sub-MODIS pixel scale would still result in continued increases in D I at the regional scale. Indeed, we do not know how much of a given MODIS pixel must be covered in a light algal bloom before the pixel reflectance dips below the dark ice threshold or how much of an impact processes such as overnight refreezing of the ice surface will have on diurnal variability in reflectance. Such considerations are particularly important when considering variations in area-wide D I above and/or within a few percent of the field-derived reflectance threshold of 0.45, as physical changes in the weathering crust (e.g. liquid water content) will also force variations in D I independently of variability in algal growth. Lastly, given the confounding impacts of cloud cover on MODIS observations, assessing the relative contribution of increases in the extent (D E ) versus concentration of algae (D I ) on regional variability in dark ice dynamics is not possible. We suggest, however, that intra-annual patterns in D I over ablation periods are more consistent with the progression of ice algal blooms than with dynamics in other darkening agents previously discussed.
We interpret our observations as support for the role of ice algae in controlling interannual dynamics in darkening, but note that there is currently not sufficient evidence to formally test this assertion. In particular, one major aspect of dark ice variability which ice algae cannot explain is the well-defined maximal spatial extent of dark ice, both in the south-west and GrIS-wide. However, if algal growth were the only factor causing interannual variability in dark ice presence, we would expect to see dark ice present wherever the climatological prerequisites for algal growth are met. These climatological prerequisites can be found elsewhere, most notably in the 20 km wide zone from the ice-sheet margin to the start of the common area examined in this study. This suggests that algal growth controlled by climatology alone cannot fully explain dark ice dynamics in the south-west sector of the GrIS. In light of our findings, we hypothesise that interannual variability in dark ice presence -both in the south-west sector and GrIS-wide -requires (1) particulates outcropping from ancient ablating ice and (2) blooming of ice algal assemblages. Specifically, we suggest that the outcropping particulates defines the spatial extent of dark ice. Algal blooms are likely to exert a control on dark ice intensity in areas where outcropping particulates are present but only when the climatological prerequisites for growth are also met. For our hypothesis to be correct, outcropping particulates must enable ice algae growth of sufficient magnitude to cause appreciable darkening, for instance as a source of nutrients, but they do not need to be present in high enough concentrations to cause darkening by themselves, nor do they need to have mineralogical characteristics that appear dark in the visible spectrum. In addition our observations imply that over interannual timescales in the south-west GrIS, outcropping dust is always present due either to continuous delivery or to long residence times. This is unlikely to hold over longer timescales if darkening by algal assemblages is contingent on the delivery of outcropping dust from ancient ablating ice.

Conclusions
We detected hitherto overlooked dynamics of the dark ice zone of south-west of the GrIS using remotely sensed imagery. Our results show that GrIS dark ice dynamics must be examined across the full duration of the melt season in order to understand the processes most likely to be reducing the albedo of bare ice surfaces. We found that in years when the south-west sector of the GrIS darkens, this usually occurs within several days and then remains widespread for the rest of the melt season, indicating that the darkening occurs in response to a common synoptic forcing. The seasons of longest dark ice duration (D D ) tend to be associated with earlier retreat of the winter snowpack. Once the ice goes dark, the dark ice intensity (D I ) tends to increase gradually through the melt season. Daily variations in D I are fairly small.
In our analysis, the JJA sensible heat flux anomaly and the date of bare ice appearance represent the most important climatic controls on dark ice extent (D E ), D D and D I , with higher sensible heat fluxes (associated with higher wind speeds) and earlier bare ice appearance favouring more dark ice. Higher JJA air temperatures and a greater number of days during JJA on which continuous surface melting occurs are also associated with darker years. There is a positive correlation between D N and the JJA GBI (R 2 = 0.46, p < 0.05), which indicates that the climatic conditions which drive darker years can be attributed at least partly to the summer presence of high-pressure blocking systems over the ice sheet.
Our observations suggest that neither deposition of particulates nor cryoconite hole processes can independently explain interannual variability in dark ice presence. Our observations tentatively support the proposal that algal blooming is the primary cause of albedo reductions in dark years, likely driven by earlier winter snowpack retreat and positive sensible heat flux anomalies. However, climatological controls on algal growth alone cannot explain the spatial distribution of interannual dark ice presence. We therefore suggest that interannual variability in dark ice in the south-west sector of the GrIS is enabled first by the deposition of particulates from melting ancient ice. These particulates play an as-yet unknown role in facilitating the growth of ice algal assemblages, which is also controlled by physical and/or climatic prerequisites that remain to be identified conclusively.
Future research has several key challenges. First, the spatial distribution, mineralogy and ice-darkening potential of all LAIs outcropping from ancient ice needs to be quantified. Second, the spatial distribution and hence ice-darkening potential of ice algae needs to be examined not just at plot scales but also at scales of hundreds of metres and more. Third, if algal cells are found to be abundant and to be the primary driver of dark ice, then the physical/climatic and nutrient controls on the growth of ice algae need to be established. Last, all these findings should be assimilated into a model of ice surface albedo that can be embedded within a regional climate model in order to project the impact of dark ice upon runoff from the GrIS during the 21st century.

Appendix A: Choice of spectral thresholds
We validated the spectral thresholds used in this study through comparison to hemispherical-conical reflectance factor (HCRF) measurements made in the field on 19 July 2016 in the vicinity of S6 (67 • 4 28.6 N, 49 • 21 32.4 W). We made HCRF measurements for three qualitatively identified surface types: (1) white ice, (2) light algal bloom (characterised by a light brown colouration to the ice surface) and (3) heavy algal bloom (characterised by a dark brown colouration to the ice surface). Subsequent microscopic examination of surface samples confirmed the presence of algal cells, along with a very low concentration of mostly clear quartz particles. The HCRF measurements were made following the HCRF measurement protocol described by Cook et al. (2017b). Briefly, an ASD Field Spec Pro spectral radiometer with an 8 • fore optic was positioned 30 cm above the sample surface with a nadir-viewing angle. This device measures reflected radiance in the wavelength range 350-2500 nm and therefore senses reflected radiance over about 95 % of the solar spectrum. The sample surface was qualitatively homogenous in a buffer zone of at least 30 cm around the viewing footprint of the sensor. We calculated the mean of at least 20 sample replicates, all of which were made within 1 min without changing the sensor position. All measurements were acquired within a 2 h sampling window around solar noon, thereby minimising error due to changing solar zenith. The sky was cloud-free throughout the measurement window. Naturally illuminated nadirview HCRF is reported for consistency with the reported MOD09GA data.
The field spectra (Fig. A1) show that the bare ice threshold used by Shimada et al. (2016) adequately captures whiteice surfaces. Their threshold of R 620−670 nm < 0.4 to define dark ice is conservative and prevents positive identification of light algal blooms. In this study we used a threshold of R 620−670 nm < 0.45, which is set to just below our field observations of light algal bloom reflectance in order to reduce the likelihood of false positives.  Figure A1. Field HCRF spectra acquired on the GrIS (see Appendix A). For each surface type, solid lines denote mean reflectance and the shaded bounds are delimited by the minimum and maximum reflectances. The grey shaded box corresponds to MODIS Band 2 (841-876 nm) and the red shaded box to MODIS Band 2 (620-670 nm). White divisions in each box correspond to the spectral thresholds utilised in this study to define bare and dark ice areas.