Glacier algae accelerate melt rates on the western Greenland Ice Sheet

. Melting of the Greenland Ice Sheet (GrIS) is the largest single contributor to eustatic sea level and is amplified by the growth of pigmented algae on the ice surface that increase solar radiation absorption. This biological albedo reducing effect and its impact upon sea level rise has not previously been quantified. Here, we combine field spectroscopy with a novel radiative transfer model, supervised classification of UAV and satellite remote sensing data and runoff modelling to calculate biologically-driven ice surface ablation and compare it to the albedo reducing effects of local mineral dust. We demonstrate that algal growth led to an additional 5.5 – 8.0 Gt of runoff from the western sector of the GrIS in summer 2016, representing 6 - 9% of the total. Our analysis confirms the importance of the biological albedo feedback and that its omission from predictive models leads to the systematic underestimation of Greenland’s future sea level contribution, especially because both the bare ice zones available for algal colonization and the length of the active growth season are set to expand in the future.


Introduction
Mass loss from the Greenland Ice Sheet (GrIS) has increased over the past two decades (Shepherd et al., 2012;Hanna et al., 2013) and is the largest single contributor to cryospheric sea level rise, adding 37% or 0.69 mm yr -1 between 2012-2016 (Bamber et al., 2018). This is due to enhanced surface melting (Ngheim et al., 2012) that exceeds calving losses at the ice sheet's marine-terminating margins (Enderlin et al., 2014;van den Broeke et al., 2016). Surface melting is controlled by net solar radiation, which in turn depends upon the albedo of the ice surface, making albedo a critical factor for modulating ice sheet mass loss Ryan et al. 2018a). The most profound shift in albedo occurs when the winter snow retreats to expose bare glacier ice. However, there are several linked mechanisms that then modulate the albedo of the exposed ice and determine its rate of melting, including meltwater accumulation, ice surface weathering and the accumulation of light absorbing particles (LAPs), such as soot (Flanner et al., 2007) and mineral dust (Skiles et al., 2017). Photosynthetic algae and cyanobacteria also reduce the albedo of the GrIS (Uetake et al., 2010;Yallop et al., 2012;Stibal et al., 2017;Ryan et al., 2017Ryan et al., , 2018b, but their effects have not yet been quantified, mapped or incorporated into any predictive surface mass balance models Noel et al., 2017;. Hence, biological growth may play an important yet under-appreciated role in the melting of the Greenland Ice Sheet and its contributions to sea level rise (Benning et al., 2014).
The snow-free surface of the GrIS has a conspicuous dark stripe along the western margin which expands and contracts seasonally, covering 4 -10% of the ablating bare ice area (Shimada et al., 2016). The extent and darkness of this "Dark Zone" may be biologically and/or geologically controlled (Wientjes et al., 2011;Tedstone et al., 2017;Stibal et al., 2017).The algal community on the GrIS is dominated by Mesotaenium berggrenii, and Ancylonema nordenskioldii (Yallop et al., 2012;Stibal et al., 2017;Williamson et al., 2018;Lutz et al. 2018). The presence of these algae reduces the albedo of the ice surface, mostly due to a brown-purple purpurogallin-like pigment (Williamson et al., 2018;Stibal et al., 2017;Remias et al., 2012). Quantification of the biological albedo reduction, radiative forcing and melt acceleration has remained elusive due to the difficulty of separating biological from non-biological albedo reducing processes. Further, there is a need to determine diagnostic biosignatures for remote sensing, which would confirm the presence or absence of ice algal growth across the melting bare ice areas of the ice sheet. Upscaling of unmanned aerial vehicle (UAV) observations made in a small sector of the Dark Zone to satellite data was prohibited by a lack of spectral resolution and ground validation (Ryan et al., 2018). The Dark Zone is of the order 10 5 km 2 in extent and is undergoing long term expansion (Shimada et al., 2016;Tedstone et al., 2017). Quantifying the impact of algal colonization on the Dark Zone is therefore paramount.
Here, we directly address these issues, resolving a major knowledge gap limiting our ability to forecast ice sheet melt rates into the future. First, we use spectroscopy to quantify the effect of algae on albedo and radiative forcing in ice. We then use radiative transfer modelling to isolate the effects of individual light absorbing particles on the ice surface for the first time, enabling a comparison between local mineral dust and algae. To examine spatial coverage, we apply a supervisedclassification algorithm (random forest) to map ice algae in multispectral UAV and satellite data. Runoff modelling informed by our empirical measurements and remote sensing observations enables us to estimate the biological contribution to GrIS runoff for the first time.

Field Site
Experiments were carried out at the Black and Bloom Project field site (67.04 N, 49.07 W), near the IMAU Automatic Weather Station 'S6' on the western Greenland Ice Sheet between 10 -22 nd July 2017. The field site location is displayed in Fig 3C. Some ancillary directional reflectance spectra included in our remote sensing training set for the supervised classification algorithms were obtained between 15-25 th July 2016 at the same field site.

Field Spectroscopy
Albedo was measured using an ASD (Analytical Spectral Devices, Colorado) Field-Spec Pro 3 spectroradiometer with ASD cosine collector. The cosine collector was mounted horizontally on a 1.5 m crossbar levelled on a tripod with a constant height of 50 cm above the ice surface. The cosine collector was positioned over a sample site, connected to the spectroradiometer using an ASD fibre optic, then the spectroradiometer was controlled remotely from a laptop, meaning the operators could move away from the instrument to avoid shading it. Two upwards and two downwards looking measurements were made in a total time period of < 5 minutes at each site to account for any change in atmospheric conditions although the measurements presented were all made during constant conditions of clear skies at solar noon +/-1 hour. Each retrieval was the average of > 20 site replicates. The albedo spectra are provided in full along with all environmental and instrument metadata, and the codes used to process the raw data in our data repository.
3 Immediately after making the albedo measurements, the cosine collector was replaced with an 8 degree collimating lens, enabling a nadir-view hemispheric conical reflectance factor (HCRF) measurement and an albedo measurement to be obtained for each sample surface. The instrument setup and sample surface were identical between each pair of albedo and directional reflectance measurement except for the viewing footprint. For HCRF measurements the upwards looking measurements were replaced with HCRF measurements of a flat Spectralon panel with the spectroradiometer in reflectance mode. This protocol was followed for every sample surface with both albedo and directional reflectance measurements taking less than 5 minutes. The directional reflectance measurements were used for the surface classifications in this paper because they better approximate the measurements made by orbital remote sensing platforms. Albedo measurements were used to estimate radiative forcing because the hemispheric measurement is most appropriate for studying energy balance as it is less affected by scattering anisotropy. We closely followed the methodology described by Cook et al., (2017b).

Biological Measurements
The sample surfaces were destructively sampled into sterile whirlpak bags, melted in the dark and immediately fixed with 3% gluteraldehyde. The samples were then returned to the University of Bristol and University of Sheffield where microscopic analyses were undertaken. Samples were vortexed thoroughly before 20 µL was pipetted into a Fuchs-Rosenthal haemocytometer. The haemocytometer was divided into 4x4 image areas. These were used to count a minimum of 300 cells to ensure adequate representation of species diversity (where possible -low abundance samples had as few as one cell per haemocytometer). The volume of each image area was used to calculate cells per mL. Biovolume was determined by measuring the long and short axes of at least ten cells from each species in each sample using the "measure" tool in the GNU Image Manipulation Program (GIMP). These dimensions were then used to calculate the mean volume of each species in each sample, assuming the cells to be circular cylinders (after Hildebrand, 1999 andWilliamson et al., 2018). The average volume was multiplied by the number of cells for each species and then summed to provide the total biovolume for each sample.

Mineral dust sample preparation
High algal biomass ice samples were collected in sterile sample bags and melted at ambient temperatures (5-10° C). The thawed samples were filtered onto glass fiber filters (0.7 µm pore size), from which the solids were removed into a glass jar using a stainless steel spatula. In 50 mL centrifuge tubes, the samples were treated using 30% H2O2 (w/w) (Honeywell Fluka™) to remove the organic fraction. The samples (1-2 g) were sonicated (VWR ultrasonic cleaner) in 45 mL of the H2O2 treatment for 10 min to disaggregate the material. The samples were left in the H 2O2 treatment for 48 h, after which they were centrifuged for 10 min at 4000 rpm (Eppendorf centrifuge 5810). The supernatant was removed, and the H2O2 solution was replaced. This process was repeated up to ten times until no more organic oxidation was observed. The remaining mineral fraction was washed three times in water (Sartorius arium pro ultrapure water), with centrifugation after each wash.
A 5 mg of H2O2-treated sample was suspended in 10 mL of ultrapure water. The sample was sonicated to disaggregate the grains. The suspension was dispersed onto a 0.2 µm polycarbonate filter (Sartorius Track-Etch Membrane, 0.2 µm). Once dry, a section of each filter was adhered to a stainless steel SEM stub using an adhesive carbon tab. The sample was coated with 8 nm of Ir (Agar high resolution sputter coater). Samples were imaged using a Zeiss Ultra Plus field emission scanning electron microscope (FE-SEM) operated at 20kV. The PSD in each sample was determined by counting all particles in an area of ~1mm 2 that consisted of images of 100 x 100 and 250 x 250 µm2 areas that were stitched together using SmartStitch software. The particles on the individual and stitched images were counted using ImageJ2 (Rueden et al. 2017). Statistical estimations of the PSD were obtained using OriginPro 8 using a single peak fitting and a Gaussian function on the obtained histograms.

Optical properties of mineral dust and algae
The mineral dust sample was arranged into an optically thick layer on a microscope slide which was then pressed tightly against one aperture of a Thorlabs IS200-4 2" integrating sphere. The other apertures were covered with SM05CP2C caps and the sample reflectance was measured using the same ASD Field Spec Pro 3 as was used for field measurements. The PSD and reflectance of the sample were then used to determine the complex refractive index of the bulk mineral dust mixture by inverting the Discrete Ordinates radiative transfer model (DISORT) following the methodology of Skiles et al. (2017). Briefly, the measured reflectance is used as a target for repeated runs of the DISORT model with varying refractive indices. The refractive index that gives the lowest root mean square error across the solar spectrum is determined to be the real refractive index of the bulk mixture. This is then used to forward model the optical properties of the mineral dust using Mie scattering theory. This was undertaken for mineral dust of radii 0.01 -10 µm at a resolution of 0.01 µm. The measured PSD was then used to apply a weighted average to the estimated single scattering properties to provide a bulk refractive index for the measured mineral sample. This was then added to the lookup library for the radiative transfer model. The new radiative transfer package BioSNICAR_GO was used to predict the albedo of snow and ice surfaces with our field-measured mineral dust.
Geometrical optics was employed to determine the single scattering optical properties of the ice algae, since they are large (10 3 µm 3 : far outside the domain of Mie scattering) and best approximated as circular cylinders (Hildebrand, 1999). Our approach is adapted from the geometric optics parameterisation of van Diedenhoven (2014). The inputs to the geometric optics calculations are the cell dimensions and the complex refractive index. The imaginary part of the refractive index was calculated using a mixing model based upon Cook et al., (2017b) where the absolute mass of each pigment was measured and input into the model in mg. We updated Cook et al.'s (2017)  imaginary part of the refractive index of water and the algal pigments so that the cell looks like water at wavelengths where pigments are non-absorbing. We consider this to be more physically realistic than having cells that are completely nonabsorbing at wavelengths > 0.75 µm, especially since a water fraction (Xw) is used in the calculations to represent the nonpigmented cellular components of the total cell volume. This approach also prevents the refractive index from becoming infinite when the water fraction is zero, removing the constraint that 0 < Xw < 1 from the bio-optical scheme in the original BioSNICAR model. Based upon experimental evidence in Dauchet et al (2015) for the model species C. reinhardii, the real part of the refractive index has been altered from 1.5 (in Cook et al. 2017b) to 1.4. The absorption coefficients from which the imaginary refractive index is calculated are from Dauchet et al (2015) apart from the novel purpurogallin-type phenol whose optical properties were determined empirically.

Radiative transfer modelling
The ice optical properties were calculated using a parameterisation of geometric optics adapted from van Diedenhoven et al.
(2014) instead of the Mie scattering approach taken by the original SNICAR model. A geometrical optics approach to generating ice optical properties was chosen because it enables arbitrarily large ice grains with a hexagonal columnar shape to be simulated, in order to better estimate the albedo of glacier ice where grains are large and aspherical. The optical properties of the ice grains are modelled using refractive indices from Warren and Brandt (2008). BioSNICAR_GO is a major update to the the BioSNICAR model (Cook et al., 2017b;Flanner et al., 2007) that allows the user to select to predict single scattering optical properties for ice grains, mineral dusts and algal cells whose sizes exceed the upper limit of the Mie scattering domain using geometrical optics. For the entire workflow, documentation and annotated scripts are provided in our data repository. For the radiative transfer modelling presented in this study, the following model parameters were used: Diffuse illumination, ice crystal side-length and length per vertical layer = 3, 4, 5, 8, 10 mm, layer thicknesses = 0.001, 0.01, 0.01, 0.01, 0.01 m, underlying surface albedo = 0.15, layer densities = 300, 400, 500, 600, 700 kg m -3 .

Quantifying radiative forcing and biological melt acceleration
Incoming irradiance was measured at our field site using an ASD Field Spec Pro 3 spectroradiometer; however, these data were not available at sufficient temporal resolution for our hourly IRF calculations. Therefore, we simulated incoming irradiance at hourly temporal resolution and 1 nm spectral resolution using the PVSystems solar irradiance program (https://www2.pvlighthouse.com.au) following Dial et al., (2018). The spectral irradiance used to calculate IRF is available in our data repository.
The biological radiative forcing was calculated by first differencing the mean albedo for algal surfaces and the mean albedo for clean ice surfaces measured at our field site. The product of the difference and the incoming irradiance spectra provided the instantaneous power density absorbed by the algae. We assume that photosynthetic processes utilise 5% of this absorbed energy, and the remainder is conducted into the surrounding ice. Since these cells are coloured by the purple purpurogallin pigment, we assume the reflective radiative forcing to be negligible, as demonstrated by Dial et al. (2018). We approximate the daily radiative forcing by multiplying the hourly instantaneous forcing by 3600 s h -1 and summing. Melt was then determine by scaling by 10 4 to convert m 2 to cm 2 and dividing by the latent heat of fusion for melting (334 J cm 3 ). To calculate the percentage of daily melt caused by algal growth, we expressed the biological melting as a percentage of the mean ablative losses for the same day measured at a network of nine ablation stakes, converted into cm water equivalent using the ratio of density of ice (668 or 917 kg m -3 ) to density of water (1000 kg m -3 ). Two values were used for the density of ice, representing upper and lower bounds for local ice density resulting from weathering. These values were presented by Smith et al. (2017) based upon in-field measurements close to our field site during the same time of year. Using the upper value for ice density assumed the surface was composed of solid ice not weathered ice, making this a conservative estimate of biological contribution, whereas the lower value is more realistic for weathered summer ice in this area. Lower density ice represents a smaller water equivalent per unit surface lowering, such that the percentage accounted for by biological melting is greater.

UAV remote sensing
UAV mapping took place over a 200 x 200 m sample area that was kept pristine throughout the study period. Inside the sampling area we placed fifteen 10 x 10 cm Ground Control points (GCPs) whose precise location was measured using a Trimble differential GPS. At these markers we also made ground spectral measurements using an ASD-Field Spec Pro 3 immediately after each flight.
We integrated a MicaSense Red-Edge multispectral camera onto a Steadidrone Mavrik-M quadcopter. The camera was remotely triggered through the autopilot which was programmed along with the flight coordinates oin the open-source software Mission Planner. Multiband TIFF images were acquired at approximately 2 cm ground resolution with 60% overlap and 40% sidelap. The flights were less than 20 minutes long and at an altitude of 30 m above the ice surface.
We applied radiometric calibration and geometric distortion correction procedures to acquired imagery following MicaSense procedures (https://support.micasense.com/hc/en-us/articles/115000351194-RedEdge-Camera-Radiometric-Calibration-Model and https://github.com/atedstone/micasense_calibration). We then converted from radiance to reflectance using timedependent regression between images of the MicaSense Calibrated Reflectance Panel acquired before and after each flight.
Finally, the individual reflectance-corrected images were mosaiced using AgiSoft PhotoScan following procedures developed by the USGS (https://uas.usgs.gov/pdf/PhotoScanProcessingMicaSenseMar2017.pdf), yielding a multi-spectral ortho-mosaic with 5 cm ground resolution, georectified to our GCPs. While the agreement between the ground, UAV and 7 190 195 200 205 210 satellite derived albedo is generally good, there are some noticeable differences that we believe to be the result of different radiometric calibration techniques for satellite, UAV and ground measurements and the differing degrees of spatial integration.

Sentinel-2 Data Processing
Sentinel-2 L1C data were downloaded from SentinelHub (Sinergise, Slovenia). The ESA Sen2Cor algorithm was used to convert the images to L2A (surface reflectance).

Supervised Classification Algorithms
To map and quantify spatial coverage of algae over the ice sheet surface we employed a novel supervised classification scheme wherein a random forest classifier was trained on field spectra collected on the ice surface and then applied to spectral images gathered using a UAV and the Sentinel-2 satellite. We selected all tiles within the latitudinal range for the dark zone calculated by Tedstone et al. 2017 that were cloud free on 21 st July +/-3 days. Our directional reflectance measurements were first reduced to reflectance values at five key wavelengths coincident with the centre wavelengths measured by the MicaSense Red-Edge camera mounted to our UAV (blue: 0.475, green: 0.560, red: 0.668, red-edge: 0.717, NIR: 0.840 μm). A data table was produced with these five wavelengths as columns and a separate row for each individual sample surface measured during our field spectroscopy. This data table was then shuffled and split into a training set (80%) and a test set (20%). The training set was then used to train three individual supervised classification algorithms: Naive Bayes, K-nearest neighbours (KNN) and support vector machine (SVM). For the SVM, the parameters C and gamma were tuned using grid search cross validation. Two ensemble classifiers were also trained: a voting classifier that combined the predictions of each of the three individual classifiers, and a random forest (RF) algorithm. The performance of each classifier was measured using precision, accuracy, recall and F1 score and also by plotting the confusion matrix and normalised confusion matrix for each classifier. In all cases the RF outperformed the other classifiers according to all available metrics.
The performance was then measured on the test set, demonstrating the algorithm's ability to generalise to data outside of the training set. Overfitting is not usually associated with the RF classifier, and the strong performance on both our training and test sets confirms that the model generalizes well. The RF algorithm was then applied to the processed UAV image. This protocol was then repeated for Sentinel-2 imagery. In that case the directional reflectance data was reduced to eight bands coincident with the centre wavelengths measured by Sentinel-2 at 20m ground resolution (0.48, 0.56, 0.665, 0.705, 0.740, 0.788, 0.865, 1.610 μm). Sentinel-2 imagery was masked using the MeASUREs Greenland Ice Mapping Project ice mask (https://nsidc.org/data/nsidc-0714) to eliminate non-ice areas from our spatial analyses. Pixels with more than 50% probability of being obscured by cloud were masked using the Sentinel-2 L2A cloud product generated by the Sen2Cor processor. Training the classifiers using data from field spectroscopy ensures the quality of each labelled datapoint in the training set. The entire workflow was achieved using bespoke Python scripts that are available in our repository. For both the UAV and Sentinel classifiers, the final model are available in our data repository.

Comparing 2016 and 2017 using MODIS
Dark ice extent and duration, and snow depths and clearing dates, were calculated after Tedstone et al. (2017), where full methods may be found. Briefly, we used the MOD09GA Daily Land Surface Reflectance Collection 6 product to map bare and dark ice in 2016 and 2017. The MOD09GA uses reflectance from Terra and sensor degradation has been accounted for in the Collection 6 product. We detected bare ice and then dark ice within bare ice areas by applying thresholds to reflectance values. For bare ice we adopted R 0.841-0.871 μm < 0.60. For dark ice we used R 0.62 -0. 67 μm < 0.45. We used the common dark ice area defined by Tedstone et al. (2017) to define the spatial sampling area for comparing 2016 and 2017.
Annual dark ice extent corresponds to the extent (in km 2 ) covered by the pixels within the common area which were dark for at least 5 days each year. Annual duration was defined at each pixel in the common 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. The timing of bare ice appearance within the common area was calculated from MODIS data using a rolling-window approach on each pixel. Mean snow depths over the common area were extracted from outputs of the regional climate model MAR v3.8 , run at 7.5 km resolution forced by ECMWF ERA-Interim reanalysis data (Dee et al., 2011). To determine spatial coverage by each surface class in each year we identified those Sentinel-2 tiles that were cloud free and available on 21 st July +/-3 days in 2016 and 2017. For direct comparisons between the two years only tiles that were available and cloud-free in both years were used. The surface coverage counts were pooled from all available tiles in each year and the mean coverage assumed to be representative across the dark zone. In 2016, additional cloud-free tiles were available so for calculating the 2016 runoff an additional 2 tiles were included in the analysis.

Runoff Modelling
Runoff was calculated using a SMB model forced with local automatic weather station and MODIS albedo observations . The model interpolates meteorological and radiative measurements from three PROMICE automatic weather stations on the K-Transect (KAN_L, KAN_M and KAN_U) and bins them into twenty 100 m elevation bands (0 to 2,000 m a.s.l.). Surface albedo is adjusted from MODIS Terra MOD10A1 albedo and is averaged into the same 100 m elevation bins. For every one-hour time step, the model iteratively solves the surface energy balance for the surface temperature. If energy components cannot be balanced due to the 0 C surface temperature limit, a surplus energy sink for melting of snow or ice is included. If surface temperature is greater than the melting point, the surplus energy is used for melting of snow or ice. When calculating turbulent heat fluxes, aerodynamic surface roughness for momentum was set to 0.02 and 1 mm for snow and ice, respectively. We extrapolate modelled runoff across the western GrIS by deriving the areas of each elevation bin using the Greenland Ice Mapping Project (GIMP) DEM  constrained to the latitude range defined in Fig 4 (following Tedstone et al., 2017). Total summer runoff from bare ice was calculated by summing runoff in elevation bins that had mean daily albedo of less than 0.60. Total summer runoff from dark ice only was calculated in the same way but using a 0.39 threshold.
We multiplied the predicted runoff by the upper and lower estimates of melt attributed to algae from our radiative forcing experiments, weighted by the relative coverage determined in our remote sensing experiments. This provided an upper and lower estimate of the amount of melting attributed to biology.

Algae reduce ice albedo
Field measurements were made 38 km inland of the western margin of the GrIS (near Kangerlussuaq, Greenland) in July 2017 (Fig 1A, 3C). This site is within the GrIS Dark Zone, close to the IMAU Weather Station S6. The albedo of patches of the ice surface was measured using an ASD Field Spectrometer with a cosine collector, rotated to look upwards and downwards on the end of a 1.5 m horizontal tripod arm. These measurements were followed immediately by the physical removal of the upper 2 cm of the ice surface within the same patches. These ice samples were fixed in glutaraldehyde for mineral dust and algal cell identification and quantification via microscopy. The dominant species of algae were Mesotaenium bergrennii and Ancylonema nordenskioldii (Fig 1B), confirming observations made by Stibal et al. (2017) and Williamson et al. (2018) in the same region. Their long, thin and approximately cylindrical morphology has been shown to be near-optimal for light absorption (Kirk, 1976). Samples were divided into the following distinct classes, based upon qualitative observations of algal abundance: High algal abundance (Hbio), Low algal abundance (Lbio), Clean Ice (CI) and Snow (SN). The measured algal cell abundance (in cells mL -1 ) in each surface class was as follows: Hbio = 2.9 x 10 4 +/-2.01 x 10 4 , Lbio = 4.73 x 10 3 +/-2.57 x 10 3 , CI = 625 +/-381, SN = 0 +/-0. These cell abundances were significantly different between the classes (one-way ANOVA, F = 10.21, p = 3 x 10 -5 ) which Bonferroni-corrected t-tests indicated to be due to variance between all four groups. The albedo of the ice surface also varied significantly between the surface classes (oneway ANOVA for broadband albedo: F= 7.9, p = 2.8 x 10 -4 ; spectrally-resolved ANOVA shown in Supp Info 1F), again with Bonferroni-corrected t-tests showing variance between all four groups (Supp Info 1C,D). Greater algal abundance was associated with lower albedo. The albedo reduction was concentrated in the visible wavelengths ( Fig 1C) where both solar energy receipt and algal absorption peak (Cook et al., 2017b;Williamson et al., 2018), diminishing towards longer near infra-red (NIR: > 0.70 μm) wavelengths where ice absorption, represented by the effective grain size, is most likely to cause albedo differences (Warren, 1982). A strong inverse correlation (Pearson's R = 0.75, p = 2.74 x 10 -9 ) was observed between the natural logarithm of algal cell abundance (cells mL -1 ) in the surface ice samples and broadband albedo (Fig 1D). The linear regression coefficient of determination was 0.56, which is unsurprising since the physical structure of the ice also plays a primary role in controlling albedo independent of any light absorbing impurities (Warren, 1982). An inverse relationship was also observed between broadband albedo and biovolume (calculated as the sum of the products of the mean measured cell volumes and the cell counts for each algal species) but the coefficient of determination was lower (r 2 = 0.42).
This may well be the result of larger cells having a smaller effect on albedo than more numerous, smaller cells for a given total volume. The relationship between absorption and scattering coefficients and cell size may also not be straightforward for algal cells due to an increasingly important contribution to the cell optical properties from internal heterogeneity, organelles, cell walls and the packaging effect in larger cells (Morel and Bricaud, 1981;Haardt and Maske, 1987).
The albedo of Hbio and Lbio surfaces is depressed in the visible wavelengths (0.40 -0.70 μm, Fig 1A), creating a 'red-edge' spectrum commonly used in other environments as a marker for photosynthetic pigments (Seager et al., 2005). Chlorophyll has a specific absorption feature at 0.68 μm which is hard to discern in the raw spectra, but clear in the derivative spectra (Fig 2A) for Hbio and Lbio but not CI and SN. This feature has previously been described as 'uniquely biological' (Painter et al., 2001) and supports the hypothesis that the albedo reduction observed in these samples is primarily due to algae. Our measurements therefore strongly indicate a biological role in reducing the albedo of the GrIS surface; however to test that the lower broadband and spectral (Fig 1C) albedo observed on algal surfaces is primarily due to the presence of algal cells, it is also necessary to compare the albedo reducing effects of the algae to that of local mineral dust.

Algae have greater impact on albedo than mineral dust
We used radiative transfer modelling to compare the albedo reducing effects of local mineral particles and algal cells. The model BioSNICAR_GO was used, which uses Mie theory to model the optical properties of mineral dusts, geometric optics to model the optical properties of ice grains and algal cells (whose dimensions far exceed the upper limit of the Mie scattering domain) and a two stream radiative transfer code (Flanner et al., 2007;. To compare algal and mineral albedo reducing effects, the model was run with fixed irradiance and ice physical properties that were chosen to reduce the absolute error between the simulated albedo for ice without any impurities and our mean measured clean ice spectrum. To determine the optical properties of the mineral dust, organic matter was chemically removed from field samples and the spectral reflectance of the cleaned minerals was measured using a spectrometer and integrating sphere, providing a target spectrum for the DISORT radiative transfer model, thereby enabling the spectral complex refractive index of the minerals (Fig 2 B) to be estimated following Skiles et al., (2017). Given the measured mineral dust particle size distribution (PSD, Supp Info 1E), this enabled the optical properties of the bulk mixture of mineral dust to be calculated using Mie theory. The minerals were much more reflective than those measured by Skiles et al. (2017)  mineral mass was made up of very small particles. While the mean diameter was 1.86 μm, this was influenced by the presence of infrequent very large fragments. 75% of the particles had diameter < 1.68 μm, and 50% had diameter < 0.56 μm.
The optical properties of glacier algae were calculated using empirical measurements of pigment mass fractions from field samples, an empirically derived absorption spectrum for the purpurogallin-like phenolic pigment, a pigment mixing model (adapted from Cook et al. 2017a,b), a measured size distribution and a geometrical optics code that assumes chains of cells to be circular cylinders, after Lee and Pilon (2013).
The effect of adding local mineral dust to a simulated ice column was a small increase in broadband albedo (Fig 2C). In contrast, adding algal cells to the simulated ice caused an albedo reduction due to the enhanced absorption of incident light in the visible wavelengths. For example, adding a hypothetical 300 µg/g mineral dust increased the ice broadband albedo by 0.02, whereas the same mass mixing ratio of algal cells decreased the broadband albedo by 0.03. This mass mixing ratio was chosen to be within the range of cell masses measured in our field samples for algal ice with the cells concentrated into the upper 1 mm of ice to match our observations. The increase in albedo caused by the addition of mineral dust may seem like a surprising result, but the ice itself had a relatively low albedo due to having a long absorbing path length (our simulated ice column had ice grains with diameter 3 -10 mm) and the local mineral dust particles were sufficiently small, weakly absorbing and strongly scattering that their overall effect was to increase the light scattered skywards by the ice. In simulations with smaller ice grain radii (400 µm radii snow grains) and therefore higher initial albedo, the addition of 300 µg/g mineral dust had negligible effect (< 0.001) whereas 300 µg/g algae still reduced the albedo by > 0.02. These simulations indicate that mineral dust is not directly causing the albedo decline on the GrIS, although they may influence the ice albedo indirectly by acting as substrates for the formation of low albedo microbial-mineral aggregates known as cryoconite granules, which are often found in quasi-cylindrical melt holes or scattered over ice surfaces. This is seemingly in contrast to previous studies such as Wientjes et al. (2010;2011) and Bøggild (2010); however, we point out that neither of the Wientjes et al. (2010;2011) studies directly measured the surface albedo or any optical properties of the mineral dusts retrieved from their GrIS sampling sites and only inferred mineralogical darkening from low spectral resolution MODIS data that did not consider biological darkening as a potential explanation for the suppressed albedo in the visible wavelengths and the "wavy pattern" observed in the dark zone in MODIS imagery. We argue that while the "wavy pattern" may be indicative of geological outcropping onto the ablation zone, it does not necessarily follow that these minerals are responsible for surface darkening, but perhaps act as stimuli for biological growth in situ. In support of this, Wientjes et al. (2011) found strongly scattering and weakly absorbing quartz to be the dominant mineral in surface ice and speculated that biota may be having a darkening effect. Bøggild et al. (2010) found mineral dust to be an albedo reducer in Kronprinz Christian's Land (80N, 24W) but this area is geologically and climatologically distinct from our field site, and their transect only spanned ~8 km from the ice sheet margin, being an area prone to local dust deposition and not thought to mineral dust an enabler of biological albedo decline and exacerbating the so-called "wavy pattern" (Wientjes et al. (2011) created by outcropping dusts. In contrast, we have demonstrated that algae are potent albedo reducers. This is consistent with previous studies Yallop et al. 2012) that found mineral dust to be insignificant for explaining albedo variations in the same region.

Indirect effects of algae
Algae predominantly reduce the ice albedo in the visible wavelengths (0.40 -0.70 μm), whereas variations in the NIR result mainly from changes to ice grain radii and the presence of liquid water (Warren, 1982;Green et al., 2002). We compared the area of an absorption feature centered at 1.03 μm between the different surface types, finding significant differences between all four surface classes (one-way ANOVA, F = 12.8, p =7.16 x 10 -7 ) driven predominantly by variations between the two algal surfaces and the two clean surfaces. This absorption feature is linked to the optical properties of snow because it scales with grain size (Nolin and Dozier, 2000), so we interpret these variations as evidence that the optical properties of the ice surface differed between the surface classes, having an effect on the measured albedo. The feature area is smallest for H bio followed by Lbio, CI and largest for SN (Supp Info 1B, Supp Info 3). The features with the smaller areas also had lower albedo minima. The absorption features are also subtly, but systematically, left-asymmetric for the algal surfaces, consistent with the presence of liquid water in the fast-melting ice beneath algal blooms (Green et al., 1998;Cook et al., 2017b). These observations suggest that the lower albedo of algal surfaces is not explained entirely by enhanced absorption due to algae, but also due to the smoother, wetter ice surface with fewer opportunities for high-angle scattering of photons, compared to the well-drained and porous CI surfaces. Cause and effect is unclear because algae may cause this by enhancing melting of the weathered surface or may grow preferentially where there is already more melt. We expect the explanation to be a combination of these two interlinked processes, especially since melting liberates nutrients that stimulate algal growth. This supports the role of indirect feedbacks (Cook et al., 2017a,b) in biological darkening of ice sheets. This process is selfamplifying because algal growth is stimulated by melt, which can be enhanced by algal growth (Yallop et al., 2012;Ganey et al., 2017;Stibal et al., 2017;Cook et al., 2017a,b;Dial et al., 2018), enhancing the albedo lowering process.

Algae enhance radiative forcing and melt
Having determined that ice algae reduce the ice surface albedo, we took an empirical approach to quantifying their impact upon energy balance following Ganey et al. (2017), which includes both direct albedo effects (enhanced absorption of shortwave solar radiation by the algal cells) and the indirect effects explained above. We used the product of the difference in spectral albedo between algal and clean ice surfaces and the incoming spectral irradiance to calculate the hourly radiative forcing (RF) of algae, assuming mineral dust is not causing significant albedo reduction. Integrated over the entire day, this indicated a daily mean biological RF of 116 W m -2 and 65 W m -2 for Hbio and Lbio surfaces respectively, similar to RFs for Alaskan snow algae calculated by Ganey et al (2017). We used the biological radiative forcing integrated over the entire day and the latent heat of fusion for ice (334 J cm -3 ) to provide an estimate of melting due to algae under the assumption that the cold content of the ice is depleted (i.e. ice is at 0°C). We made two calculations with ice densities 917 kg m -3 (solid ice) and 668 kg m -3 (weathered ice) following Smith et al (2017), providing upper and lower bounds for melt attributed to algae depending upon local ice density. The melting on 21 st July was estimated to be 1.3 to 1.9 cm w.e. in Hbio areas, which represents between 21 and 29% of the total melting measured across a network of ablation stakes at our site (6.8 cm). For Lbio sites, biological melting on 21 st July 2017 was 0.76 to 1.07 cm w.e., which corresponds to between 12 and 18% of the observed ablation. These estimates are similar to those made by Ganey et al. (2017) for heavy algal blooms on snow (21% of total melt due to ablation).

Algae cover a large proportion of the ice sheet
Our analyses demonstrate that algae have a dramatic darkening effect on the ice surface, leading to increased melting.
However, the importance of this effect depends upon the spatial extent of the algal blooms over thousands of kilometers. We made 146 directional reflectance measurements immediately after albedo measurements for each sample surface, providing training data for supervised classification algorithms that were applied to multispectral imagery obtained by a UAV.
Directional reflectance is a more appropriate measurement than albedo for this purpose because it better approximates the measurements made by orbital remote sensing platforms and are less affected by surface heterogeneity, having smaller viewing footprints. Our ground spectra were reduced to five bands matching the centre-wavelengths measured by our UAV multispectral camera. The reflectance dataset was divided randomly into training (80%) and test (20%) sets. The this data was then used to train a random forest algorithm (chosen because of its high performance relative to other classifiers) to predict the surface class for each pixel in our UAV image, thereby enabling spatial upscaling of our field spectroscopic measurements. The accuracy, precision, recall and F1 score of the random forest classifier on the test set were all 95%.
The UAV image was obtained by flying a MicaSense Red-Edge camera integrated onto a quadcopter over a 0.04 km 2 area grid with a ground resolution of 5 cm at our field site on 21 st July 2017. The classified image indicated that 78.5% of the area was covered by algal blooms of which 61.1% was Lbio and 17.4% was Hbio (Table 2; Figure 3). The high ground resolution of the imagery enabled a qualitative assessment of the algorithm performance by visual comparison between the classifier and the raw imagery (following Ryan et al. 2018). The algorithm produced qualitatively realistic bloom shapes, correctly placed water in channels and individual cryoconite holes in their correct positions. The confusion matrix indicates that occasional misclassifications are generally between water and cryoconite (Supp Info 2). This is unsurprising since both cryoconite and water have relatively flat spectral shapes with few spectral features and cryoconite is often found beneath pools of surface water. We also point out that our cryoconite spectral reflectance measurements were made with cryoconite filling the entire field of view of the spectrometer, so best represent large cryoconite holes or dispersed cryoconite rather than surfaces peppered with many small holes. There was also some ambiguity between thin, wet snow and bare glacier ice, as these surfaces are spectrally similar. Nevertheless, these misclassifications affect a small area of the pixel and do not affect our estimate of algal bloom coverage.
A training set was then produced for application to Sentinel-2 satellite data by reducing the hyperspectral ground data to nine bands, coincident with those measured by Sentinel-2 at 20 m ground resolution. The confusion matrices (Supp Info 2) indicate similar misclassification types and frequencies to the UAV model. The RF model was applied to Sentinel-2 tiles that were within the Dark Zone, cloud-free and available on the date of our UAV flight +/-3 days. Non-ice areas were masked out prior to analyses. The classifier then predicted the surface type pixelwise across all Sentinel-2 tiles. The resulting data were then pooled into one large dataset. The mean algal coverage across all tiles was 18%, however there was significant spatial heterogeneity, with the tile covering our field site in the Kangerlussuaq region having 44% algal coverage. H bio surfaces were less common than Lbio ( Table 2). The spatial coverage by algae was different in the Sentinel and UAV datasets especially for Hbio, likely because a) the Sentinel-2 imagery includes ice that is outside of the Dark Zone, raising the overall reflectivity, and b) even in the UAV image, which was retrieved from within the Dark Zone, Hbio surfaces comprise just 15% of the ice surface and have a patchy distribution, meaning they may not be detected by the 20 m resolution Sentinel-2 data.
The lowest albedo surfaces -cryoconite and water -cover a small fraction (< 3%) of the total area in both UAV and Sentinel-2 images (Table 2), although we note that many individual cryoconite holes will not be detected due to being smaller than the spatial resolution of either Sentinel-2 or the UAV. The spatial coverages reported here from our multispectral UAV imagery are consistent with a k-nearest neighbours classification scheme applied to RGB (Red, Green, Blue) imagery from a fixed wing UAV flight over the Kangerlussuaq region by Ryan et al. (2018). They found up to 85% of the ice surface to be composed of 'ice containing uniformly distributed impurities' in the same region of the Dark Zone in July 2014, which our observations confirm were dominated by algae. They also found < 2% of the ice surface to be cryoconite covered and water coverage was < 5% (except for a supraglacial lake in their imaged area). This analysis demonstrates that algae are a major component of the ice surface.

Algae reduce ice sheet albedo at the landscape scale
Multispectral imagery acquired by the UAV was converted to albedo using the narrowband to broadband conversion of Knap et al. (1998). The additional bands for the Sentinel-2 data enabled the application of Liang et al.'s (2002) narrowband to broadband conversion. There was a significant difference between the albedos of each surface class in all four datasets, consistent with the findings from our ground spectroscopy (Figure 3; Table 1). The albedo of each surface class is approximately consistent between the datasets, despite the variation in spatial coverage (Table 2) Our satellite remote sensing data demonstrates that algal blooms are a major component of the ice surface in this area (Table   1,2, Figure 3, Supp Info 4), and where they are present the ice albedo is on average 0.13 lower for Lbio and 0.20 lower for Hbio compared to clean ice (Table 1). This, combined with our ground-based spectroscopy and radiative forcing calculations, provides robust evidence in support of algae having a significant melt-accelerating effect on the GrIS. We cannot yet explicitly separate mineral and biological effects, but our theoretical and empirical analyses indicate that: a) local mineral dust cannot explain the observed albedo reduction, b) low albedo areas had significantly elevated algal cell numbers relative to clean ice, c) uniquely biological features were detectable in the spectra and derivative spectra for the lower albedo sites, and d) radiative transfer models incorporating algal cells with realistic pigment profiles demonstrate the mechanism of albedo reduction. These reasons confirm that supervised classification of Hbio and Lbio surfaces is indeed detecting surfaces with high algal loading and can be used to make a conservative estimate of algal bloom extent. It is conservative because there is certain to be ice algae present in low numbers in some of the areas that are classified as clean.

Enhanced algal albedo reduction in a 'dark year'
Satellite remote sensing using MODIS (Figure 4) indicates that 2017 was a particularly high albedo year, where the Dark Zone was both smaller and brighter than in previous years. Figure 4 Figure 4C,D shows that summer 2016 was preceded by much thinner snow, which melted away to reveal bare ice over a month earlier than in 2017. Furthermore, there were several additional snowfall events (5-10 cm snow) during our field work period in 2017, which did not occur in the same period during 2016. Previous field evidence (Williamson et al., 2018) demonstrates that the ice was darkened by high concentrations of algae in 2016. We therefore applied our classification algorithms to Sentinel-2 data from the same date and precisely the same locations for both 2016 and for 2017 ( Figure 3, Table 1, Table 2 4). Interestingly, the Lbio coverage extends further inland in the darker year, whereas Hbio surface extend further towards the margin. This suggests that newly exposed ice at higher elevations was colonised by relatively low abundance blooms, whereas heavier blooms developed on ice at lower elevations that was exposed for longer. This suggests the intensity of the algal bloom is a function of exposure time (Tedstone et al., 2017) and corroborates the findings of space- interannual comparison demonstrates that more prolonged exposure of larger ablation areas under a warming climate (Stroeve et al. 2013;Shimada et al. 2016;Tedesco et al. 2016;Tedstone et al. 2017) are likely to be prone to more spatially expansive, darker algal blooms that enhance melt rates, leading to a potential positive feedback that is not currently accounted for in surface mass balance models. The existence of a bright area at the edge of the ablation zone may seem to undermine exposure time as a driver of algal growth; however, we do not discount the importance of outcropping mineral dusts from stratified ice layers which may well be critical nutrient stimuli for glacier algae which then darken the ice surface.
The emergence of these minerals may impose a lower boundary on the dark zone, or alternatively the melt runoff may be sufficient to wash light absorbing particles off the ice surface at these very low elevations. These observations therefore indicate that the omission of biological growth is leading current models to underestimate future GrIS contributions to sea level rise.

Algae cause enhanced GrIS runoff
We ran a surface mass balance (SMB) model forced with local automatic weather station and MODIS albedo observations  to estimate 94.1 Gt of runoff from bare ice on the western GrIS (albedo < 0.6) and 67.6 Gt runoff from dark ice (albedo < 0.39) within the latitudinal range defined in Fig 4. We used the mean spatial coverage determined using our remote sensing in each year and our radiative forcing calculations that attributed 12-18% of melting to algae in Lbio sites and 21-29% in Hbio sites to generate upper and lower estimates for the GrIS runoff caused by algal growth. We found that in 2017 between 0.08 -0.11 Gt of ice loss could be attributed to the growth of algae, representing 1 -2% of the total runoff from the western GrIS. In 2016 this contribution increased to 5.5 -8.0 Gt, representing 6 -9% of the total runoff from the western GrIS. These calculations confirm that algal growth is an important factor in contribution of the GrIS to global sea level rise. The interannual comparison strongly indicates that this contribution will increase if biologically-darkened areas expand or a greater proportion of the ice is covered by high biomass blooms under warmer climates.

Conclusions
Our measurements and modelling demonstrate that the growth of algae on the GrIS is accelerating the rate of melting and increasing the GrIS contribution to global sea level rise. Our field spectra show a dramatic depression of the surface albedo in the visible wavelengths for surfaces contaminated by algae. Derivative analysis of the same spectra show "uniquely biological" absorption features and an inverse relationship was observed between biomass and surface albedo. We employ a novel radiative transfer model to show that this albedo decline cannot be attributed to local mineral dusts. We demonstrate that the growth of algae occurs over a large proportion of the ablating area of the GrIS by training a random forest classifier to identify algal blooms in remote sensing data from a UAV and Sentinel 2. In the particularly dark 2016 melt season, coverage by algae was much higher than in the bright 2017 melt season. Our field measurements and remote sensing data inform a runoff model that estimates 5.5 -8.0 Gt or 6-9 % of the total runoff in summer 2016 could be attributed to the growth of algae. This study therefore unequivocally demonstrates that algae are important albedo-reducers and cause a meltenhancing feedback across the GrIS. The omission of these critical biological albedo feedbacks from predictive models of GrIS runoff is leading to underestimation of future ice mass loss and contribution to global sea level rise. This is particularly significant because larger ablation zones and longer growth seasons are expected in a future warmer climate.

Data Availability
Codes and datasets used in this study are available at the following doi's: The active repositories for ongoing development of BioSNICAR_GO and the ice surface classifiers codes are at www.github.com/jmcook1186/BioSNICAR_GO and www.github.com/jmcook1186/IceSurfClassifiers respectively.

Author Contribution Statement
JC developed the measurement protocol, gathered field measurements, analysed the data, wrote the main codes, curated the data repository, produced the figures and wrote the manuscript. OMcA was instrumental in building and testing the UAV. AT, CW, JMcC, SH, TG gathered crucial field data. CW provided essential advice regarding microscopy and biological sampling protocols and helped with experimental design, and also led the empirical measurements of glacier algae pigmentation and absorption coefficients. AT wrote the code for radiometric calibration of multispectral imagery from the UAV and post-   . Each year is labelled with dark ice extent. In each year, pixels that are dark for fewer than 5 days are not shown. (C,D) Average snow depth modelled by MAR (blue) and cumulative dark ice extent observed by MODIS (red) (Tedstone et al., 2017) during April to August. Vertical bars (grey) denote median date of snow clearing derived from MODIS; horizontal bars denote the interquartile range of the day of year of bare ice appearance.

815
The Cryosphere Discuss.  Table 1: A) summary of the albedo for each surface class as predicted from our classified UAV image. B) summary of the albedo for each surface class as predicted from our classified Sentinel-2 image for 2017. C) summary of the albedo for each surface class as predicted from our classified Sentinel-2 image for 2016. D) summary of the broadband albedo for each surface class as measured using field spectroscopy at our field site in 2017 (we do not have cosine-collector albedo measurements for water or cryoconite surfaces).