The role of microbes in snowmelt and radiative forcing on an Alaskan icefield

A lack of liquid water limits life on glaciers worldwide but specialized microbes still colonize these environments. These microbes reduce surface albedo, which, in turn, could lead to warming and enhanced glacier melt. Here we present results from a replicated, controlled field experiment to quantify the impact of microbes on snowmelt in red-snow communities. Addition of nitrogen–phosphorous–potassium fertilizer increased alga cell counts nearly fourfold, to levels similar to nitrogen–phosphorusenriched lakes; water alone increased counts by half. The manipulated alga abundance explained a third of the observed variability in snowmelt. Using a normalized-di erence spectral index we estimated alga abundance from satellite imagery and calculated microbial contribution to snowmelt on an icefield of 1,900 km. The red-snow area extended over about 700 km, and in this area we determined that microbial communities were responsible for 17% of the total snowmelt there. Our results support hypotheses that snow-dwelling microbes increase glacier melt directly in a bio-geophysical feedback by lowering albedo and indirectly by exposing low-albedo glacier ice. Radiative forcing due to perennial populations of microbes may match that of non-living particulates at high latitudes. Their contribution to climate warming is likely to grow with increased melt and nutrient input.

Glacier microbiomes are water-limited 21,22 , because ice is generally not metabolically available, and oligotrophic, because their nutrient content equals that of precipitation plus deposition by airborne dust, pollen, and so on, with only limited N-fixation by local cyanobacteria 21,22 . Moreover, rapidly percolating water through large-grained snow may exacerbate both water-and nutrient limitation for algae in supraglacial snow-covered habitat 7 . Surface meltwater on glaciers derives principally from radiative flux into the snowpack 2 , with the spectral distribution of albedo determining absorbed insolation 1,3 . Thus, snow algae increase liquid water availability through 'bioalbedo' 18 , determined by the absorptive properties of their photosynthetic and secondary pigments. The chlorophylls and carotenoids in red-snow algae-particularly the heat-dissipating tetraterpenoid astaxanthin that colours cells red as protection from excess radiation 23,24 -reduce snow-surface albedo in the visible spectrum. Nutrient depletion in the snowpack triggers production of carotenoids 20 in non-motile cells that ultimately overwinter as buried cysts. Spring meltwater then stimulates germination when propagules actively disperse through the snowpack, absorbing nutrients, photosynthesizing, and reproducing there and at the surface 20 .
This study experimentally demonstrates the potent role of redsnow algae in glacier melt. It presents novel results that establish these algae as both nutrient-and water-limited, and shows redsnow algae satisfy two assumptions of a bio-geophysical feedback loop for snowmelt. By manipulating alga abundance and simultaneously measuring melt on an icefield's surface, a causal, quantitative 'alga-melt model' is provided. Simultaneous field measures of spectral reflectance and alga abundance motivate a normalizeddifference snow-darkening index (NDI) to estimate alga abundance through remote sensing. Instantaneous radiative forcing (IRF) is calculated as a function of alga abundance using abundance-specific reflectance spectra integrated with instantaneous irradiance. Spatially mapping red-snow algae with satellite-derived NDI, together with the alga-IRF and alga-melt models, allows estimation of the red-snow ecosystem's contribution to snowmelt on Alaska's Harding Icefield (60 • N, 150 • W; Supplementary Fig. 1). Adjacent to the Gulf of Alaska, this rapidly thinning 25 network of >30 glaciers covers 1,890 km 2 and drains a snow-and ice-covered plateau, providing a potential model for ice sheets in Greenland and Antarctica.
Interval photography (Supplementary Movie) and an automated weather station on the Eklutna Glacier (61.21 • N 148.96 • W; Supplementary Fig. 1) demonstrated previous observations [16][17][18][19][20][21][22][23][24] that the appearance of red snow coincides with increased melt rate and temporal trends in melt predictors (Fig. 1a). A PCA rotation of melt predictors showed that the first principal component (PC1) explained 89% of overall variance in melt predictors and 31% of the variability in daily melt rate (Fig. 1b,  red snow, melt predictors, and melt rate are consistent with red-snow algae acting as both a cause and an effect of snowmelt. We focus here on the red-snow sub-season, when snow algae are most likely involved in snowmelt.

Nutrients and water limit snow-alga abundance
Experimental additions of an aqueous solution of nitrogenphosphorus-potassium (NPK) fertilizer, and of water alone, showed strong nutrient-and moderate water limitation in snow algae on the Harding Icefield because their abundance increased with the addition of each (Fig. 2a, Table 1). Treatment effects also manifested as timing of visible red snow, first detected from repeat photography (22 June) and in experimental plots (12 July). By 19 July, >90% of nutrient-enriched plots, 50% of water and control plots, but only 10% of alga-reduction plots showed red snow. A single pulse of N-enrichment in-frame at the start of interval photography on the Eklutna Glacier also suggested a response of red-snow algae to nutrient enrichment, with earlier, more intensely red-coloured snow at and near the enrichment point (Supplementary Movie). Previous field observations 16,26 and in vitro microcosm nutrient additions 17,27 suggest that nutrient limitation in snow algae is study-dependent. For example, Lutz et al. 16 in a wide-ranging Arctic study found no correlation between alga abundance and nutrient concentrations in filtered red snow. In contrast, Fujii et al. 26 attributed Antarctic snow-alga distribution and abundance to allochthonous bird deposition of limiting nutrients. A field microcosm study on North American volcanoes 27 found no NP-enrichment effect on snow-alga abundance, while a laboratory microcosm of Greenland cryoconite debris 17 showed NPenrichment doubled both organic carbon and chlorophyll-a. Here, the response to NPK-enrichment in 2-m 2 plots exceeds the mean response of approximately 3.4× reported for pelagic phytoplankton in NP-enriched freshwater lakes 28 .
This experimental study tests for liquid-water limitation in snow algae, a critical first assumption of the feedback hypothesis; however, snow-alga's elevation range and life-cycle anticipate the observed increase in cell counts with water-only treatment. Drier than lowelevation snow, cold, high-elevation snow lacks alga colonization 29 . Classic studies of snow algae 22 suggest that the appearance of liquid water and light in the spring initiates germination.

Manipulating alga abundance manipulates melt
The second assumption of the feedback hypothesis, that algae cause melt, was also supported (Supplementary Movie). Experimentally varying alga abundance correspondingly varied melt rate at the plot level following algae's appearance (22 June; Supplementary Fig. 4). Alga abundance predicted melt rate well as a square-root function of cell surface area (alga-melt model), suggesting a possible asymptote ( Fig. 2b and Supplementary Table 1). Six weeks and three treatments into the 100-day experiment, melt rates were highest on average in enriched plots and lowest in reduced plots ( Supplementary  Fig. 4). Melt rate appeared decoupled from degree-day sum after 12 July, when red snow was visible. Including alga abundance or treatment as an interaction term significantly improved the  Table 2). c, Treatment e ects on cell counts. Line segments link plot counts within experimental blocks; slope equal to response ratio (treatment:control count). Font size (Control = 1) of text equals geometric-mean response ratio averaged across blocks. Black arrow points at zero count. degree-day model after 12 July when comparing traditional degreeday melt models using a Likelihood Ratio Test (LRT: cell area: p < 0.001, AIC = 21.7; treatment: p = 0.003, AIC = 7.9); however, before 12 July, including treatment did not improve model performance (p = 0.51). Degree-day models, commonly applied in spatial calculations of snowmelt 30 , have not previously accounted for bioalbedo on snow. Plots with experimentally increased algae also exposed bare ice earlier than reduction and control plots (Supplementary Movie). Twice as many enriched plots (28% total), 50% more water plots (20%), and two-thirds as many reduced plots (10%) as control plots (14%) melted to ice or saturated firn by 10 August. Enriched plots were almost three times more likely than reduced plots to melt to bare ice or slush. By directly ablating clean snow (albedo >0.75) from glacial surfaces to reveal bare ice (albedo <0.35 16 ) that may not otherwise be exposed during the melt season, and by exposing ice sooner, algae indirectly amplify their melt effects 16 , and do so repeatedly due to their perennial life history.
We present three lines of evidence that experimental treatments manipulated melt through alga abundance and bioalbedo, rather than through physio-chemical processes. First, observed effects on melt occurred only after the appearance of red snow. There were no significant treatment effects on melt ( Supplementary  Fig. 4) until four treatment periods into the experiment. If physical processes related to treatments caused melt, then treatment effects would probably have occurred earlier. Second, alga-melt models (Supplementary Table 2) with and without treatment as a factor variable did not differ (LRT interaction p = 0.69; additive p = 0.83).
This is a statistical result consistent with a negligible role for physical processes. Third, an enriched plot (1.4 × 10 4 cell ml −1 ) and a reduced plot (no cells detected) differed in melt rates as ∆H = 10 mm w.e. d −1 -where w.e. is water equivalent-well within the range predicted by a published physical model 18 Supplementary Fig. 4). The generally smaller differences than theoretically predicted 18 could be due to non-zero alga counts in reduction plots (Fig. 2c), an increase in snow-grain size, and so melt by HOCl, or both. Overall the results that suggest increased melt with increasing microbes occurred primarily through biologically driven radiative forcing, not physio-chemical change in grain size.

Remote sensing of microbial impact
The relationship between alga abundance and melt (alga-melt model; Fig. 2b) allowed estimation of snow algae's contribution to snowmelt across the Harding Icefield using Landsat-8 satellite imagery. We sampled alga abundance and measured snow-surface reflectance using a field spectrometer ( Fig. 3a) for comparison to modelled reflectance spectra 3 for dust and black carbon (BC) at commonly observed concentrations (Fig. 3b,c). The field-measured spectrum of the lowest density alga sample (7 ppm) was similar to the modelled spectrum of moderate density dust (172 ppm) with absorption <650 nm. Modelled BC reflectance spectra were uniform at these concentrations (7-333 ppb). The alga reflectance spectra here generally agreed with previously published spectra for red-snow alga. The spectrum at 333 ppm was similar in shape, but with lower values of reflectance, to a theoretical spectral albedo of snow with very high densities of algae 18 . The spectrum at 172 ppm closely matched the field reflectance from a California sample 9 of 2.1 × 10 4 cells ml −1 and less closely a red-snow sample from the Harding Icefield 10 , two studies that observed strong linear relationships between alga abundance and indices based on reflectance features in the 500-700 nm range. We, too, observed a strong linear relationship (NDI-alga model, Fig. 4a, R 2 = 0.93) between field measures of alga abundance and a normalized snowdarkening index (NDI), where reflectance, R i , is measured in the spectral regions corresponding to Landsat-8 bands 4 (red: 636-673 nm) and 3 (green: 533-590 nm). | Spectral reflectance varies by particle abundance and source. a, Mean reflectance spectra (mean ± s.e.m., n = 10) for snow samples of varying alga abundance. b,c, Modelled reflectance spectra for dust (b) and black carbon (uncoated in ppb) (c) using SNICAR. NDI = (R red − R green ) (R red + R green ) −1 for each curve. d, Instantaneous radiative forcing by alga abundance. IRF = [(R A (λ) − R o (λ))]IR(λ)∆λ over 350 ≤ λ ≤ 850, with R A (λ) = reflection spectrum for alga abundance A, R o (λ) for alga-free snow and ∆λ = 10 nm. IR(λ) = direct plus indirect solar irradiance normal to insolation at 60 • N 150 • W on 29 July 2013 at 11:09 local sidereal time. Red line is best fit through origin, black line is given by equation with R 2 .
Because R red ≤ R green on clean ice and snow surfaces 1 , we interpreted NDI ≤ 0 on cloud-free snow as a sufficient indicator for an absence of red-snow algae at the pixel scale. The NDI index developed here for algae (Fig. 3a) has been applied to measure mineral dust (Fig. 3b) on alpine snow 31 , and is insensitive to BC (Fig. 3c), which lacks algae's strong and dust's weak absorption peak <600 nm (refs 18,23,31,32) (see Fig. 3a,c). NDI is more sensitive to alga abundance than equivalently common particle concentrations: at 172-333 ppm NDI > 0.10 for algae but NDI < 0.05 for dust and NDI = 0 for BC at 172-333 ppb. We did not attempt mixture modelling of algae, dust, and BC using NDI, but suggest this as an avenue for further study, particularly regarding dust particles that vary by colour. Here, four lines of evidence suggest NDI measured primarily algae on the Landsat-8 images of the Harding Icefield during the 2013 melt season used in this study: no large fires or volcanic eruptions occurred nearby; the maritime icefield is distant from dust and soot sources; the maritime icefield is subject primarily to coastal, windward onshore winds and continental, leeward downglacier winds; and we found good visual agreement between a GPS track and satellite-predicted alga abundance on a Landsat-8 scene captured three days before we surveyed an isolated 1.25 ha dark red-snow patch on the icefield's surface with a sub-metre GPS (Fig. 4b). Nevertheless, it is most likely that the spectra in Fig. 3a represent a mixture of algae, dust, BC, and other substances; however, at this point we have no estimates for those values.
We present a landscape estimate of snowmelt due directly to snow algae. Pixels with NDI ≤ 0 were omitted and those with NDI > 0 (red-snow pixels) used to algebraically compose the alga-melt model ( Fig. 2b) with the NDI-alga model (Fig. 4a), yielding snowmelt due to algae. On a cloud-free, unsaturated Landsat-8 scene (29 July 2013) we spatially mapped total snowmelt in red-snow pixel-i as  Expressed as a percentage increase of snowmelt in the absence of snow algae, bioalbedo increased melt 21%, an empirical estimate within the lower range of modelled predictions 18 (20-40%) for bioalbedo on a snow-covered Greenland glacier. Additional indirect effects (for example, unveiling ice) are not included in this single day estimate. Musilova et al. 17 calculated that 5% of the Greenland ice sheet's annual meltwater can be attributed to microbes dwelling in cryoconite. In general, microbial contribution to discharge will depend on the porosity of the glacier and its hypsometry, scaling with the length of the glacier's equilibrium-line altitude contour (snow algae's primary habitat 21,29 ) and the cosecant of terrain slope there.
High albedo of snow-covered glaciers and ice sheets diminishes solar forcing 2-4 . Thus, red-snow algae could affect climate by reducing surface albedo. How great an effect can be assessed as a first approximation by calculating instantaneous radiative forcing (IRF). We calculated IRF using modelled irradiance and fieldmeasured reflectance values (Fig. 3a) sampled simultaneously with alga abundance. IRF was modelled as in refs 4,31 for each of seven alga densities (172-333 ppm) by subtracting each spectrum's wavelength-specific reflectance value from that observed for algae at 7 ppm (indistinguishable from approximately 170 ppm of dust; Fig. 3). Regressing IRF on alga abundance (alga-IRF model, Fig. 3d, R 2 = 0.92) predicted IRF across the Harding Icefield ( Supplementary Fig. 5) with the Landsat-8 scene used for snowmelt by composing the alga-IRF model (Fig. 3d) with the NDI-ppm abundance model (Supplementary Table 1). The area-averaged IRF was 21.6 W m −2 (sd = 5.9 W m −2 ), with a maximum (87.9 W m −2 ), similar to instantaneous radiative forcing values for large-scale dust events 4,31 . Multiplying melt rate (PI 95% 0.38-0.52 cm d −1 ) due to algae by pixel area and latent heat of fusion (334 J g −1 ) suggested a daily-averaged radiative forcing PI 95% of 14.7-20.1 W m −2 , near to but less than the area-averaged value of IRF (21.6 ± 5.9 W m −2 ).

Implications for high-latitude ice sheets
This study highlights the substantial impact of red-snow communities on glacier melt at high elevations and latitudes. Experimental results presented here, together with previous correlative observations [8][9][10][11][12][13][14][15][16] , laboratory experiments 17 , and theoretical calculations 18 , provide a compelling case for the magnitude of the glacier microbiome's effect on hydrology and climate. Snow algae amplify their albedo reduction through life history, population growth, dispersal, and physiology. They probably evolved under natural selection to actively melt snow, thereby generating needed liquid in a frozen environment. Snow algae contrast with episodic deposition of non-living particulates, which sink into the snowpack if hydrophilic or non-polar, and remain atop a single season's melting snowpack if hydrophobic 33,34 . Because algae remain on the surface over much of the melt season and perennially resurface across melt seasons, they compound their effects over time. Like arctic shrubs 35 , high-latitude microbes potentially drive a warming climate feedback if their radiative forcing at the regional scale increases the extent of near-freezing snowpack, their primary habitat. Given an upward-elevation shift with warming, algae will increase most rapidly across flat, snowcovered topography, such as Greenlandic and Antarctic ice sheets, regions with critical albedo effects on global climate 2,3,36 . Worldwide ash from biomass burning 37 and dust from agricultural regions 38 are increasingly deposited on these high-latitude ice sheets. This airborne nutrient input, together with growing meltwater availability 39 , will certainly increase the glacier microbiome's impact on polar albedo. Climate and melt models that ignore the ecology of microbial radiative forcing risk underestimating rates of warming and consequent sea level rise.

Methods
Methods, including statements of data availability and any associated accession codes and references, are available in the online version of this paper.

Methods
Melt predictors and the emergence of red snow. Standard variables used for mass-balance studies 40 on glaciers (for example, temperature, snow depth, depth-specific density, incoming and outgoing solar radiation) and interval photography of snow surfaces were collected from south-central Alaska's Harding Icefield on the Kenai Peninsula and the Eklutna Glacier in the Chugach Mountains ( Supplementary Fig. 1). The Harding Icefield has been the site of both previous mass-balance 25 and glacial biota studies 10,41 . We supplemented data from an automated weather station (AWS; RAWS Harding Icefield station http://www.raws.dri.edu/cgi-bin/rawMAIN.pl?akAHAR) 210 m above and <3 km distant from our experimental study site (60.156029 • N 149.794006 • W) with three temperature sensors (Onset, model HOBO Pendant datalogger in radiation shield) maintained 2 m above the snow surface at the vertices of the experimental study site (Supplementary Fig. 1) from 2 May to 10 August 2014. Positive degree-day (>0 • C) sum was spatially modelled across the study site through a polynomial interpolation (ArcGIS, vers. 10.2) using temperature data collected from the vertices. Interval photography documented phenology of red snow. The alpine Eklutna Glacier has undergone mass-balance study since 2008 42  Snowmelt was calculated as the observed change in snow thickness converted to water equivalent (cm w.e.) based on density observations from snow profiles. Snowmelt (cm w.e. d −1 ) over a given time interval was estimated by integrating depth-specific density ρ (x) over depth, where z and z − ∆ were the snow depths at time 1 and time 2, and ∆ was the difference in snow depth at time 1 and time 2. On the Eklutna Glacier, difference in depth, ∆, was calculated as the difference between sonic ranger measurements; on the Harding Icefield, as the difference between ablation cable length measurements. We measured depth-specific snow density (ref.   Fig. 1) were recorded at sub-metre accuracy (Trimble, model GeoXH 2008 series GPS) on 20 August (1,074-1,117 m asl; mean elevation 1,100 m asl). We sampled snow-alga abundance in all plots on 19 July. We investigated the predictive ability of melt using three metrics of alga abundance: density of cell counts, density of cell surface area, and density of cell volume. Cell volume has been used in previous studies 10,16,18 . Four 10 ml surface snow sub-samples were combined as a single sample per plot, melted and preserved in 3-6% formalin under ambient field conditions. Abundance metrics were estimated from 10-mega-pixel micrographs (14-23 images plot −1 ) of algae within a 1 ml Sedgewick-Rafter counting chamber viewed through an optical microscope.
Counts and diameters, x, of round cells (5 µm < diameter < 45 µm; x = 11.3 µm, sd = 2.8 µm, n = 8,049 cells) were recorded from 1.09 mm 3 samples. We used ImageJ (version 1.49n) to record cell count, n, mean diameter,x, and standard deviation, sd, in µm from each micrograph. These statistics provided bootstrap estimates of individual cell diameters, d, circular surface area, A = 0.25πd 2 , and spherical volume, V = 0.167πd 3 . For each of the 1,670 micrographs, we generated 1,000 bootstrapped samples of n diameters from a gamma distribution with scale = sd 2x −1 and shape =(xsd) 2 . For each bootstrap sample, the cell area and volume of each of the n cells were summed for the sample, then the 1,000 bootstrap samples were averaged and aggregated by plot. Some samples included cryoconite that masked uncounted cells. Oval-shaped cells were rare and excluded from counts. Abundance metrics were scaled as count density (cells ml −1 ), cell surface area density (cm 2 l −1 ), and density of cell volume (ppm).
Particle abundance and spectral reflectance. We compared spectral reflectance of red-snow patches to the alga abundance of each snow patch by measuring reflectance with a portable spectrometer and sampling alga abundance as described Immediately following spectral measurements, four 10 ml surface snow samples were collected and combined, then cells counted (n = 20 counts sample −1 ) as above for the manipulation experiment. We used the date, time, and location of the reflectance measurements to find solar zenith angle from a solar insolation (ref. 43) calculator at PV Lighthouse website (https://www.pvlighthouse.com.au). The solar angle was low enough to use diffuse incident radiation in SNICAR models 3 of particulates as black carbon (BC) and dust reflectance spectra (http://snow.engin.umich.edu). Details of parameters for SNICAR are indirect incident radiation, 1,000 µm snow-grain radius, 3 m snowpack, 500 kg m −3 density, 0.3 visible and 0.7 NIR underlying albedo, uncoated BC with MAC scaling factor = 1.0, 5-10 µm dust particles with 'Surface spectral distribution' set to 'Mid-latitude winter, clear-sky'; 'Greenland summit, clear-sky' gave similar results.
Instantaneous radiative forcing (IRF) for spectrometer measures of snow with algae was calculated as a finite sum approximating the integral where the difference in wavelength-specific reflectance between clean snow, R s (λ), and snow with particulates, R p (λ), is weighted by global irradiance, I (λ), and multiplied by the wavelength interval in the sum (here 10 nm, from 350 to 850 nm). R 7 (λ) was the reflectance spectrum from the field spectrometer for algae at 7 ppm and R i (λ) was the reflectance spectra for algae at i = ppm (Fig. 3a). We downloaded I (λ) for N 60 • E-150 • on 29 July 2014 at 13:10 AK DST (converted to local sidereal time) using https://www.pvlighthouse.com.au. This instantaneous irradiance spectrum was used to calculate IRF for each of the seven samples of alga abundance. A linear function fits the relationship of IRF regressed on alga abundance well (Supplementary Fig. 6) and is referred to as the 'alga-IRF model' .
Statistical analyses. Statistical and spatial modelling was undertaken with R (v. 3.3.2) and included packages lmer, lme4, nlme, lmtest, raster, and rgdal. Analysis of treatment effects used linear (LMM) and generalized linear mixed-effects (GLMM) models with control plots as reference in an effects parameterization of a balanced design of three treatment plots and a control (four levels of a fixed factor) per random block (m = 21 blocks; n = 84 plots total). We tested treatment effect on raw cell counts, cell density, cell area density, and cell volume density (ppm), the latter frequently used in the literature. We also used GLMMs with plot (variable numbers of micrographs per plot) and block as random effects, treatment as a fixed factor, and micrograph cell count as a negative binomial and as a Poisson response variable. We chose the negative binomial because Poisson-regression residuals were overdispersed. Scaled alga abundance (counts ml −1 , area ml −1 and volume ml −1 = ppm) was right-skewed with variance increasing with mean, so abundance was log or square root transformed for use in LMMs, establishing the effect of treatments on scaled abundance measures. Treatment effect on cell counts is reported as the geometric mean of the response ratio of treatment to control counts as in a review of enrichment experiments 28 . P values compare treatment effects to control as one-tailed Wald tests, with α = 0.05 the critical value for significant deviation from statistical null hypotheses.
We investigated treatment effect on melt rate with three analyses. One analysis derived a function that could predict melt rate given alga abundance ('alga-melt model'). We used linear regression of melt rate from 12 to 27 July on square-root-transformed alga abundance (cell surface area in cm 2 l −1 ) sampled on 19 July with all n = 84 plots and treatment as an additional independent variable.
Second, because we anticipated that treatment effects on melt rate would differ by melt sub-season depending on algal presence (see Fig. 1a), we applied an LMM (block as random, treatment as fixed effect) to each of six time periods, with a Tukey HSD post hoc test when Pr(F 3,79 ) < 0.05. In our third analysis, we parameterized a degree-day sum (DS) model for melt (H ) passing through the origin with alga abundance (A) included as an interaction term, DS. Phenomenological models apply thawing degree-day sum to predict glacial melt using temperature 30,40 . Our assumption of interaction imposes the following two conditions. First, if no days are above 0 • C (DS = 0) yet alga abundance is positive, then melt remains zero. Second, the effect of algae on snowmelt is to increase the degree-day factor, b = b o + b 1 A. We compared degree-day sum models with alga abundance (or treatment before June 22) to models without using Likelihood Ratio Tests (LRT). After we experimentally established an alga-melt effect, we determined the best form and variables of a predictive model among the scaled alga abundance measures and their various transformations, f , using simple linear regression to estimate melt, H (cm w.e. d −1 ), as a function, f , of abundance, Spatial modelling. A previous study 10 of snow algae on the Harding Icefield applied satellite remote sensing (SPOT) after sampling alga cell abundance and regressing it against field measures of spectral reflectance in the 600-700 nm range. In our study, after simultaneously measuring snow-alga abundance and reflectance values, R, we constructed two normalized-difference indices, NDI, where 'red' is Band 4 of Landsat-8 (636-673 nm) and 'green' is Band 3 (533-590 nm), and NDI blue , where the 'blue' Band 2 (452-512 nm) replaced green, Astaxanthin absorbs 44 in the 350-570 nm region (green and blue) and reflects >600 nm, supporting these band choices. Reflectance spectra, R i (λ) , at wavelength λ and alga abundance i were recorded at ∆λ = 0.5 nm intervals using a field spectrometer (see Particle abundance and spectral reflectance above). Reflectance values were averaged across each band and across samples for use in NDI. Alga abundance, A, regressed against NDI and NDI blue in simple linear models provided the 'NDI-alga model' , A = c o + c 1 NDI, used to predict alga abundance given Landsat-8 scenes. We chose NDI over NDI blue because the R 2 values were higher and the intercept was closer to zero with NDI. Among the three metrics of alga abundance, NDI predicted cell surface area best (Supplementary Table 1).
Note that if R red − R green ≤ 0, when R red < R green , then NDI ≤ 0. This occurs when algae are absent, and so pixel-i with NDI i ≤ 0 was removed from further analysis, such as melt or IRF estimates across the Harding Icefield. Only pixels with NDI > 0 ('red-snow' pixels) were used.
As a measure of validation of NDI for identifying algae from Landsat-8 scenes, with a sub-metre GPS (Trimble, model GeoXH 2008) we outlined the boundaries of an isolated 1.25 ha, dark red-snow patch on the Harding Icefield at the south-west corner of the study area ( Supplementary Fig. 1) on 17 August 2013. We compared the outline of this red-snow patch to red snow predicted from a Landsat-8 scene (LC80680182013226LGN00, collected 14 August 2013). The Landsat scene was transformed to NDI values, with NDI ≤ 0 removed. Then we applied the NDI-alga (ppm) model to predict alga abundance (Supplementary Table 1) to generate alga abundance pixel by pixel and overlayed the GPS path.
We estimated snowmelt due to algae on the Harding Icefield for a single day as follows. A cloud-free, unsaturated Landsat-8 scene (LC80680182013210LGN01) of the Harding Icefield acquired at 21:09 GMT (11:09 local sidereal time) on 29 July 2013 was transformed to NDI values with NDI ≤ 0 removed, leaving only red-snow pixels. Pixels with NDI < 0 were removed because a negative NDI occurs if and only if reflectance in the green band surpasses reflectance in the red band (R red < R green ), which suggests red-snow algae are absent based on reflectance spectra of algae, dust, black carbon (Fig. 3) and clean snow. Each pixel i with a positive NDI received a snow-alga abundance value, A i , estimated with the NDI-alga model as A i = c 1 NDI i , where c 1 was estimated as the regression coefficient using least-squares linear regression through the origin. The predicted alga abundance in pixel i was then algebraically composed with the alga-melt model as per-pixel snowmelt rate, Per-pixel i snowmelt rate multiplied by 30 m pixel area (9 × 10 6 cm 2 ) gave daily volume flux for each pixel i in ml d −1 , as v i = 9H i × 10 6 . The volume flux was then summed across all pixels on the Harding Icefield with non-zero NDI (that is, 'red-snow pixels' where algae were present) as an estimate of total icefield snowmelt, T, where algae were detected using NDI. We then set alga abundance equal to zero in the alga-melt model for each red-snow pixel as H o = s o + s 1 f (0) = s o , converted to volume flux = 9s o × 10 6 , and summed across red-snow pixels to get the non-alga-derived snowmelt, G. The difference, T − G, gave alga-derived snowmelt on the day of imagery scaled to Gl d −1 by multiplying by 10 −12 . The 95% prediction error was found as T ± ε. If N = the number of red-snow pixels, then total melt across the red-snow area was estimated as the sum over i = 1, . . . , N pixels. i + σ 2 as the total prediction error, where for pixel i the value se i was pixel i 's standard error of predicted melt and σ was the standard deviation of residuals from the regression of melt on square root of algal area. Pixel values of expected melt were found using the predict.lm function in R applied to extracted pixel values of algal area, themselves found using NDI values.
We used a similar method of function composition to calculate IRF across the red-snow pixels (NDI > 0) of the Harding Icefield at the day and time of the Landsat-8 scene capture used for melt estimation. In this case, however, both functions were linear least squares. First we used the regression of alga abundance in ppm against NDI as an NDI-alga model, ppm = p o + p 1 NDI, (Supplementary Table 1), then we used the regression of IRF against alga abundance through the origin as IRF = g 1 ppm ('alga-IRF model' , Supplementary Fig. 5). In summary, given red-snow pixel i with Landsat-scene value NDI i , the pixel's instantaneous radiative forcing was IRF i = g 1 ppm = g 1 p o + g 1 p 1 NDI i .
In the case of IRF we did not calculate the total contribution to radiative forcing over time or over the icefield, but instead determined the statistical distribution of pixel values of IRF to compare to the heat of enthalpy needed to generate predicted melt on the Harding Icefield.
We checked the melt rate against instantaneous radiative forcing using the latent heat of fusion, Ω = 334 J g −1 . near to but less than the area-averaged value of IRF (21.6 ± 5.9 W m −2 ) found above by composing the alga-IRF model (Fig. 3d) with the NDI-ppm model (Supplementary Table 1 Code availability. The code that supports the findings of this study is available from the corresponding author upon request. Data availability. Landsat-8 images LC80680182013226LGN00 (Product identifier on Earth Explorer LC08_L1TP_068018_20130814_20170309_01_T1) and LC80680182013210LGN00 (Product identifier LC08_L1TP_068018_20130729_20170309_01_T1) were downloaded from USGS Earth Explorer (https://earthexplorer.usgs.gov). Reflectance spectra for dust and black carbon were downloaded from SNICAR online (http://snow.engin.umich.edu). Insolation data for satellite scene date, time, and location were downloaded from the PV Lighthouse website (https://www.pvlighthouse.com.au) solar insolation calculator. Automated weather data on the Harding Icefield were downloaded from RAWS Harding Icefield (http://www.raws.dri.edu/cgi-bin/rawMAIN.pl?akAHAR). The data that support the findings of this study are available from the corresponding author upon request.