Introduction

The highly reflective snow surfaces of the world contribute substantially to maintaining the energy balance of the earth. Clean snow surfaces reflect as much as 99% of incoming solar radiation, protecting our atmosphere from warming1. Even a small-scale reduction in snow cover extent or snow albedo will drastically reduce the protective properties of the cryosphere and contribute to a warming climate2,3.

A variety of microorganisms grow in snowpacks across the cryosphere, including bacteria4,5, annelids6,7, chytrids8, fungi9, and algae10,11,12,13,14. Snow algae bloom on the surface of the snowpack during the summer months when there is sufficient interstitial water to provide the necessary habitat for the algae. They develop colored photoprotective pigments to protect the cells from intense solar radiation present at the surface of the snowpack15,16,17. These photoprotective pigments, primarily astaxanthin, color the snow red where they bloom and cause the albedo of the snow to decrease by around 20%18,19,20. Quantifying the impacts of this reduced albedo across glaciers and snowfields would be prohibitive with in situ studies. With the increasing capabilities of uncrewed aerial vehicles (UAV), we can evaluate the prevalence and albedo influence of snow algae across large, inaccessible areas.

Previous studies have utilized UAVs to evaluate a variety of research questions in the cryosphere. The advances in Structure-from-Motion (SfM) technology allow researchers to construct high-resolution digital elevation models (DEMs) from UAV imagery. This approach has been applied as a means to evaluate snow and glacier ablation21,22,23,24, assess glacier dynamics25,26,27,28,29, and retrieve snow depth30,31,32,33,34,35,36,37. While there have been studies that use the UAV to study glacier algae on the Greenland Ice Sheet38,39, the remote detection of snow algal blooms in mid-latitudes using UAVs have received little attention.

One of the first studies to remotely estimate snow algae was conducted by the National Air and Space Administration/Jet Propulsion Laboratory Airborne Visible Infrared Imaging Spectrometer (AVIRIS) flown on a twin otter in the Sierras40. The current research on the remote detection of snow algae has largely been focused on using satellite images40,41,42,43,44,45. One of the most prevalent indices developed thus far is the red (610–680 nm) to green (500–590 nm) SPOT satellite reflectance ratio developed by Takeuchi et al.45. The Takeuchi index has been applied to other satellites, including Landsat46 and Sentinel-2A44. However, satellite images are less likely to have closely corresponding ground validation data and usually have limited spectral or spatial resolutions sufficient to capture the snow algae. Current and upcoming spaceborne imaging spectroscopy missions, such as the Italian Space Agency’s PRecursore IperSpettrale della Missione Applicativa, Hyperspectral PRecursor of the Application Mission (PRISMA), the German hyperspectral Environmental Mapping and Analysis Program (EnMAP) and NASA’s Surface Biology and Geology (SBG) mission, feature instruments with high spectral resolution, such as 30 m ground sampling distance, which is sufficient to detect algae in snow47. Along with the high spectral resolution of the recent and upcoming spaceborne imaging spectrometers described above, hyperspectral sensors mounted on UAVs may be able to detect subtle spectral features such as algae absorption in snow. For example, there is potential to use hyperspectral data to retrieve physical quantities in order to develop snow models that incorporate algal cell absorption, which could then facilitate the analysis of temporal evolution of snow algae blooms. This study seeks to downscale the red/green index utilized in previous satellite remote sensing of snow algae studies to the UAV level and validate the approaches with coupled in situ snow algae data. The new methods and approaches described here could be applied to the plethora of hyperspectral sensors that are coming on the market at more approachable prices to researchers.

Principal component analysis (PCA) is a proven remote sensing approach to differentiate spectral groups from multispectral imagery. The technique has been utilized to retrieve the spatial distribution of snow cover from MODIS and Sentinel-1 data48, to map avalanche debris from Sentinel-1 images49, and to conduct wetland change detection using a combination of Landsat MSS and TM images50. PCA allows researchers to identify land objects by utilizing the information in all bands of data instead of using a select few in an index approach. This band reduction approach can be a great help when analyzing hyperspectral and multispectral data51. While PCA has been thoroughly documented in other applications, to the author’s knowledge, it has yet to be applied to detect snow algae in the cryosphere. By utilizing the information in all bands of multispectral imagery through the PCA approach we may be able to extract snow algae coverage with greater precision and with the potential to guide future spectral classifications through the examination of the loading patterns in the principal components (PCs). By evaluating the band loading patterns of the PCs that are most useful in distinguishing snow algae we can design a more effective spectral index.

The snowpack of the North Cascades has been declining as a result of climate change with the average spring snowpack projected to decrease by 38–46% by 205052. This reduction in snow cover will cause glacier ice to be exposed to solar radiation for longer periods of time, exacerbating glacier melt. Snow algae are likely contributing to this increased snowmelt by decreasing the albedo of the snow19,53,54. Yet, they are not currently included in regional watershed melt models. In addition, the habitat extent of snow algae is expected to expand with a warming climate and a lengthened melt season55. Understanding and evaluating the extent and impacts of the snow algae at a large scale will aid our ability to predict the stability of the North Cascades snowpack.

Radiative transfer models have been previously established for clean snow, such as the Two-streAm Radiative TransfEr in Snow (TARTES)56. However, modeling the influence of biologically active snow algae on the albedo of the snow has only recently been established. The SNow, ICe, and Aeosol Radiative (SNICAR) model3 has been well validated as a tool to model the spectral albedo of snow with light-absorbing impurities and has been used in the International Panel on Climate Change assessments57,58,59,60. This model has been updated as the first of its kind to include the contribution of snow algae to albedo in a unified code base61.

In this study, we utilize the high-spectral and -spatial resolution MicaSense Dual camera system mounted on a UAV to map snow algae in a basin in the North Cascades under differing bloom conditions and employing two separate approaches with the first approach being exploratory in nature and the second being based on established spectral properties. In the first approach, we reduced the 10-band multispectral imagery into three components by conducting a principal components analysis and classifying the snow algae using the resulting components. In this way, we were able to utilize the information in all ten bands of the original imagery to classify the snow algae and determine which bands were important in the classification. Second, we tested a spectral indexing approach using only the red and green bands to map the snow algae and correlate the cell concentration to the proposed index. This second approach utilized previously established spectral properties of snow algae to guide the spectral classification. Our in situ snow algae data were then used to model the spectral albedo of the snow algae using the SNICAR model and compare the modeled spectra to the UAV-derived spectra. Finally, we quantified the radiative forcing of the snow algae across the entire study area using the UAV-derived snow algae map and calculated the resulting impact on snowmelt.

Results and discussion

Bloom conditions: cell counts and pigment concentrations

UAV flights and coincident ground sampling campaigns were conducted with the DJI M210 and MicaSense Dual Camera on July 2, 2021 and July 30, 2021. Four additional ground sampling campaigns occurred on 26 June, 9 July, 18 July, and July 22, 2021. At the time of the additional four ground campaigns, a MAPIR near-infrared camera was flown on a DJI Phantom 3. However, the MAPIR camera does not have the resolution needed to spectrally discriminate snow from snow algae (Supplemental Fig. 1). Therefore, we focus here only on the MicaSense Dual Camera results. The snow algae cell concentrations, determined by Guava Flowcytometry from the field samples, were almost eight times higher during the first flight on 2 July vs. 30 July (average cell concentrations 431,000 ± 398,000 cells/mL and 55,000 ± 60,000 cells/mL, respectively.) There were two outlying data points in the 30 July cell concentration results that increased the error in the second survey (Supplementary Table 1). Moreover, chlorophyll A concentrations were 29 times higher at the time of the 2 July survey than the 30th July survey and based on ground-based sampling and analysis of pigments and cell concentrations between the two flights, the 2 July flight was during a peak in bloom intensity (Fig. 1).

Fig. 1: Average cell and pigment concentrations.
figure 1

Average cell (a) and primary and secondary pigment concentrations of (b) Astaxanthin, (c) Lutein, (d) Chlorophyll A, (e) Chlorophyll B, and (f) Beta Carotene (n = 10 per date) collected where snow algae were present in Bagley Basin between June 26, 2021 and July 30, 2021. Boxes and lines represent the quartiles of the data with the dots being outlying data as defined as points beyond 1.5 times the interquartile range.

Full size image

All sample pigment concentrations decrease between the two MicaSense overflight dates, 2 July and 30 July, with the average concentration of lutein decreasing the most in the samples from 51.23 ± 38.16 µg/L on 2 July to 1.73 ± 3.23 µg/L on 30 July (Supplementary Table 2). There were numerous non-detect samples for lutein in the second survey that reduced the average concentration and increased the standard deviation (Supplementary Table 1). Normalizing the pigment concentrations relative to chlorophyll a, the cells have a much higher astaxanthin to chlorophyll a ratio on the second survey than on the first survey (Fig. 2). This shift is also associated with a decrease in the ratio of lutein to chlorophyll a while the ratios of chlorophyll b and beta carotene to chlorophyll a remain similar between the two surveys (Fig. 2).

Fig. 2: Distribution of snow algae pigment concentrations relative to chlorophyll a for each survey date.
figure 2

Boxes and lines represent the quartiles of the data where the dots are outlying data as defined as points beyond 1.5 times the interquartile range.

Full size image

The increase in the ratio of astaxanthin to chlorophyll a, the decrease in the ratio of lutein to chlorophyll a, and the associated decrease in snow algae cell concentration suggest that the bloom is reducing its photosynthetic pigment production and increasing its photoprotective pigment production in the later survey. This could indicate that the algae bloom during the early-season survey of 2 July was considerably more intense than during the second survey when the peak bloom has passed. This shift in pigments could also be due to a change in species composition, however, this type of analysis was outside the scope of this study. There also could potentially be a difference in the selected sampling locations between the two dates that is being reflected in the pigment shift, though we attempted to control for this as much as possible.

Principal component thresholding development

The ordination process reduced the variation in the 10-band imagery to explain 99.67% of the variation with the first PC, 0.30% with the second PC, and 0.01% with the third for the 2 July survey (Supplementary Table 5). For the 30 July survey, the ordination reduced the 10 bands with the first PC explaining 99.78% of the variation, the second PC explaining 0.19%, and the third PC explaining 0.01% (Supplementary Table 5). This PCA approach was an exploratory analysis to identify which bands can be used to explain the most variation in the data and to sufficiently distinguish snow algae without limiting ourselves to the typical snow algae bands.

The first PC of the 2 July survey equally weights all ten original MicaSense bands and acts essentially as an index of the overall image brightness. The second PC has negative correlations with the red edge and near-infrared bands and positive correlations with the rest of the bands (Table 1). This PC contrasts the brightness in the visible bands with the red edge and near-infrared wavelengths. This contrast is likely accounting for the variation between the snow pixels and the rest of the image where most of the variation in the snowpack is seen in the visible wavelengths. Plotting these first two PCs against each other using the training points provided clear separation between the snow algae and the snow classes (Fig. 3).

Table 1 Loadings of each MicaSense band for principal components 1, 2, and 3 from the July 2, 2021 and July 30, 2021 surveys.
Full size table
Fig. 3: Principal component classification points.
figure 3

Classification points for the (a) July 2, 2021 and (b) July 30, 2021 surveys with their associated principal component values. Dashed lines represent the threshold values used for snow algae classification, for 2 July these values were PC1 = −90,000 and PC2 = −2500 and for July 30, 2021 these values were PC1 = 50,000 and PC3 = −500.

Full size image

The ordination of the 30 July survey created a first PC that is composed of all 10 bands equally weighted, representing the overall image brightness (Table 1). The second PC contrasts a strong correlation with the near-infrared band and the red edge bands that is balanced against the visible bands. This PC compares the visible wavelengths with the infrared wavelengths while weighting the near-infrared and the blue bands the strongest and in opposite directions. Note that while the sign of the loadings for the first and second PCs for 2 July and 30 July are reversed, the loading patterns are the same and that is ultimately most important in the ordination process. The third PC weights the blue band, green-2 band, and red edge-3 band against the red-1, red-2, and red edge-1 bands (Table 1).

For the 2 July flight, the first two PCs provide the best distinction between the snow algae and the clean snow (Fig. 3 and Supplementary Fig. 4). Based on the training data, we set PC thresholds to classify snow algae as pixels that have a PC1 value greater than −90,000 and a PC2 value less than −2500 (Fig. 3). This index indicated that 1% of the surface snow was covered by snow algae (Supplementary Table 4).

For the 30 July survey, the first and the third PC demonstrate the best separation of the snow algae from the clean snow. Plotting the first and the second PC against each other revealed substantial overlap of the clean snow and the snow algae training points while the first and the third PC provided sufficient separation to develop classification thresholds (Supplementary Fig. 4). Based on the training points, pixels with a value less than 50,000 for the first PC and less than −500 for the third PC were classified as snow algae (Fig. 3). This classification resulted in a larger separation between the snow algae and snow classification points for the 2 July survey (Fig. 3). This index estimated that 2.06% of the remaining surface snow was covered by snow algae (Supplementary Table 4).

PCA offers a valuable way to extract spectrally unique features from multispectral imagery, yet it has not been widely used in snow algae remote sensing applications, due to limited coincident ground observations. This approach is most useful when there is associated classification training data. Without the training data it would be difficult to apply this approach at large scales unless the snow algae bloom is so distinct from the rest of the snowpack that it can be visually confirmed in the imagery. Here, snow algae cell counts, pigments, and associated GPS data enabled the development of principal component thresholds used to map the snow algae with the multispectral imagery.

The PCA method provided plenty of separation between the Snow Algae and the Snow and Other classes on the 2 July survey. Here, “Other” is used as a bulk term to incorporate other light-absorbing particles such as mineral dust and black carbon, which were not quantified in this study. During this survey, most of the study area was snow-covered at an average depth of about 2 m. Where the snow algae were not blooming, the snow was considerably clean with very few visible impurities like dirt, dust, or other debris. In comparison, on the 30 July survey, the snowpack thickness was reduced with an average depth of less than 1 m. As the snow melts, impurities in the snowpack aggregate at the surface, lowering the snow albedo as the season progresses. Though the ground-measured impurity content from the snow algae samples remained the same between the two surveys (Supplementary Table 3), the surrounding snow could have increased in impurity content and contributed to the less distinct separation between the Snow Algae and the Snow classes in the PCA ordination process for the 30 July survey.

The PCAs for the two surveys were independent ordinations and the difference between the results of each could point to biological shifts in the algal blooms. An evolution of the algal species composition could result in a difference in absorption features for the snow algae due to the associated change in relative pigment concentrations. In addition, the maturation and stabilization of the snow algae bloom towards the later end of the season could present different spectrally than young, early-season blooms and, therefore, impact the ordination results.

The PCA approach was independent from the optimized red/green band index development and intended to be exploratory. The relatively less effective classification developed based on the ordination of the 30 July flight could also be due to the reduced snowpack coverage of the basin and increased vegetation coverage. The ordination process encompasses as much of the variation in the 10-band orthomosaic into as few principal components as possible. The orthomosaic of the 2 July survey was composed almost completely of snow and rocks, with some areas of water and minimal patches of vegetation. Whereas the 30 July orthomosaic consisted primarily of vegetation, rock, and water with less than half the image being snow-covered. For the 30 July survey, the ordination process was driven by the variation in the vegetated portion of the survey area instead of the variation in the snow-covered portion of the image. Though we achieved sufficient separation of the training point classes by analyzing the entire survey area, future studies could attempt to mask the portions of the image that are not relevant to the analysis prior to ordinating the data for further distinction.

Optimized red/green band index development

The second approach, to develop an optimized red/green band index, was a separate approach from the exploratory PCA approach described above. On July 2, 2021 survey, there is a distinct separation between snow algae and clean snow with multiple band combinations providing sufficient separation for snow algae classification (Supplementary Fig. 2). While there could potentially be more than one applicable band combination, we proceeded with the red (642–658 nm) to green (524–538 nm) ratio since that one cleanly distinguished the snow algae classification points from the others (Fig. 4) and falls within the established red (610–680 nm) to green (500–590 nm) band ratio previously developed for snow algae mapping with the SPOT satellite45.

Fig. 4: MicaSense reflectance in Band 5 (642–658 nm) versus Band 3 (524–538 nm).
figure 4

MicaSense reflectance in Band 5 (642–658 nm) versus Band 3 (524–538 nm) for the classification training points on (a) July 2, 2021 and (b) July 30, 2021. The dashed line represents the equation: Band 5 = 1.02(Band 3) + 0.015.

Full size image

A snow mask was developed based on the classification points (Fig. 4) and applied to the image so that only the snow-covered portions of the study area were analyzed. Snow-covered area was defined as pixels that had a red band reflectance greater than 0.3:

$${R}_{{\lambda }_{B5}} \, > \, 0.3$$
(1)

Similarly, an optimized red/green band index was developed based on the training points. The original red/green band index classifies snow algae as any pixel with a red/green ratio greater than 1.0245. When we apply the original red/green band index to our training points, this captures the “clean snow” or snow algae-free classified training points, as well as the snow algae training points. We then use our classification points, derived from coincident ground sampling, to produce an optimized red/green band index. This enables us to then distinguish algae-free snow from snow with algae. These pixels are identified based on their reflectance (R) at a specific wavelength (λ) in a particular band (B). For example, here in the red (\({R}_{{\lambda }_{B5}}\)) and the green band (\({R}_{{\lambda }_{B3}}\)) as follows:

$$\frac{{R}_{{\lambda }_{B5}}-0.015}{1.02\times {R}_{{\lambda }_{B3}}}\, > 1$$
(2)

We then apply the optimized red/green band index (Eq. (2)) in order to quantify the area of snow that has algae (total snow algae-covered area) in the UAV survey image orthomosaics (Supplementary Table 4). Equation (2) detects the presence of snow algae based on the snow algae optical properties. Snow algae pigments absorb more in the green range than in the red resulting in an index value greater than 1. Clean snow absorbs similarly in both the red and the green wavelengths causing the index to be equal or less than 1. The optimized red/green band index (Eq. (2)) seemed to be less effective in separating the snow algae from the rest of the snow on the 30 July survey than on the 2 July survey (Fig. 4).

This process revealed that on July 2, 2021, 1.16% of the snow in the study area was covered by snow algae for a total of 1352.18 m2 (Supplementary Table 4). For the 30 July survey, there is a visible reduction in snow coverage (Fig. 5) from 116,226.95 m2 of snow on 2 July to 40,593.57 m2 on 30 July (Supplementary Table 4). However, snow algae covered a slightly higher percentage of the snow, 1.37%, on the 30 July survey for a total of 556.07 m2 of snow algae (Supplementary Table 4).

Fig. 5: MicaSense true color orthomosaics derived from UAV imagery.
figure 5

MicaSense true color orthomosaics derived from UAV imagery on (a) July 2, 2021 and (b) July 30, 2021. Areas in pink represent snow algae detected using the Eq. (2) optimized red/green band index. Results for the PCA mapping are similar to the optimized red/green index results displayed here.

Full size image

The wavelengths of the MicaSense Dual camera system align closely with those of the Sentinel 2A and Landsat 8 satellites, providing the future ability to directly upscale the UAV-derived and ground-validated indices to satellite platforms. The similarity in spectral bands and the ability to closely couple ground validation samples with survey flights could greatly advance the capabilities of satellite remote sensing for detecting spectrally subtle features like snow algae. The MicaSense camera systems have been employed by previous studies that used the spectral capabilities to map a variety of research focuses, including mapping submerged vegetation62, wetland functions63, arctic vegetation64, and glacier algae38,39. However, to the authors’ knowledge, this is the first dedicated application of the UAV and the MicaSense Dual camera system to map snow algae.

The 10 bands of the MicaSense imagery all gather useful information but parsing out which band combinations are superior depends on a robust ground validation dataset. For this study, we sought to explore the band combinations that were already established in the literature to assess their efficacy from a UAV platform with coincident ground truth data. Our study found the most success by using an optimized red/green index, originally developed for the SPOT satellite45. Based on our ground data, we added in a red band reflectance threshold and fine-tuned the index by incorporating multipliers into the ratio. Previous studies that have used the original red/green ratio have classified snow algae pixels as those that have a ratio greater than one46. Our optimized red/green index classifies snow algae pixels as those that have a red band reflectance greater than 1.02 times the green band reflectance plus 0.015 (Eq. (2)). The improved spectral classification of the optimized red/green index is enabled by the increased spectral resolution of the MicaSense camera bands vs. the SPOT satellite, for which the original red/green index was developed. This increased spectral resolution further helps to distinguish algae-free snow pixels from snow algae pixels, particularly late in the melt season when the surrounding snow algae blooms contain other impurities that are spectrally similar to the snow algae pigments, especially within the spectral resolution of the SPOT satellite. Again, this methodology points towards the potential for increased spectral classification of snow impurity types with hyperspectral sensors.

The optimized red/green index provided a clear distinction in the Snow and Snow Algae training points on the 2 July survey, likely because of the comparatively clean, early melt season snow that is spectrally distinct from any dark snow algae patches. The average inorganic impurity concentration was the same for both survey dates indicating that the impurities were not playing a role in the efficacy of the classification approach and instead the water content of the snow could be responsible (Supplementary Table 3). The optimized index did not clearly separate debris-covered snow from snow algae. This could provide a challenge for future studies that seek to classify snow algae in heavily debris-covered areas. However, the applications of imagers with more bands, such as hyperspectral imagers, may have the ability to further discern snow algae in debris-covered snow.

Based on the 100 randomly located points, both the principal component thresholding and the optimized red/green index approaches indicated that there was a slightly higher percent of snow algae covering the snow in the 30 July survey compared to the 2 July survey. The percentage of snow algae coverage was quite similar for the principal component thresholding and the optimized red/green index approach for the 2 July survey, 1.01% and 1.16%, respectively. However, they slightly differed in their results on the 30 July survey where the principal component thresholding approach determined that 2.06% of the snow surface contained snow algae compared to 1.37% determined by the optimized red/green index approach.

The difference between the results of the two approaches supports the dependency of index development on the snow algae bloom conditions, where early-season peak intensity blooms are more consistently detected than late-season post-peak blooms. Future studies should acknowledge the influence that the algal bloom state and the age of the snowpack have on these remote sensing indices when attempting to map snow algae without ground data and further coupled studies are needed to refine the temporal evolution of these parameters.

Optimized red/green band index correlations

When regressing the optimized red/green band index with the snow algae cell concentrations there is a nonsignificant relationship between the index value and the natural log of the snow algae cell concentration, with the index value only explaining 16% of the variation in the cell concentration (F1,17 = 3.232, P = 0.09, r2 = 0.16) (Fig. 6). Separating out the two survey dates, the 2 July index explains more variation in the natural log of the cell concentration in comparison to examining both surveys together (30% and 16%, respectively), but the relationship is still nonsignificant (F1,7 = 3.038, P = 0.13, r2 = 0.30) (Fig. 6). The relationships are much stronger with the 30 July survey index. Examining only the 30 July data, there is a significant relationship between the optimized red/green band index and the natural log of the snow algae cell concentration, where the index explains 86% of the variation in the cell concentration (F1,8 = 50.23, P < 0.001, r2 = 0.86) (Fig. 6).

Fig. 6: Correlation between the optimized red/green band index with the natural log of the cell concentration by date.
figure 6

The black dashed line represents the linear regression of all the points where ln(y) = 4.488x + 6.188 (F1,17 = 3.232, P = 0.09, r2 = 0.16), the red solid line represents the linear regression of just July 2, 2021 points where ln(y) = 8.185x + 3.038 (F1,7 = 3.038, P = 0.13, r2 = 0.30), and the blue solid line represents the linear regression of just July 30, 2021 points where ln(y) = 5.0575x + 4.4519 (F1,8 = 50.23, P < 0.001, r2 = 0.86). Gray shading represents the result of local polynomial regression fitting for each date.

Full size image

There is also a significant positive correlation between the optimized red/green band index and the ratio of astaxanthin to chlorophyll a, however, the index only explains 24% of the variation in the ratio of astaxanthin to chlorophyll a (F1,17 = 5.30, P = 0.034, r2 = 0.24) (Supplemental Fig. 5). This relationship is not improved by analyzing the two survey dates separately.

We developed the regression model to predict algal concentration from the UAV imagery with the goal of estimating the algal abundance over the entire survey area. However, this study showed that the ability of the optimized red/green index to predict snow algal concentration varies greatly depending on the time of the year. During the early summer season, the algal blooms have much greater range in cell and pigment concentrations that does not correlate well with the proposed index. Yet, later in the summer season, the algal blooms seem to stabilize with lower variation in the cell and pigment concentrations allowing for the proposed index to have greater predictive power. The difference in the predictive power of the index could also be due to a shift in algal species composition in the snowpack and an associated change in relative pigment concentrations. When both surveys are combined, the optimized red/green index only explains 16% of the variation in the natural log of cell concentration—not enough to confidently apply the index at large scales without associated ground data. Though, future studies could collaborate with numerical modeling to simulate snow algal blooming and further develop this methodology65.

The results of this study demonstrate that imagery and sampling must be closely coupled when developing algorithms to map snow algae. Over the three weeks between our two surveys, the predictive ability of the index shifts greatly. While the optimized red/green index can be used to capture the presence/absence of snow algae, the variation in cell and pigment concentration reveals the need for the development of multiple indices to predict snow algae bloom intensity and/or the need for high-resolution spectral data from hyperspectral sensors to enable quantitative retrieval of algal concentration based on other methods, such as, model inversion. In addition, ingestion of these observational data into numerical models that simulate snow algal blooming (i.e., Onuma et al.65) could be powerful in order to further extrapolate radiative transfer impacts on snowmelt across larger domains, such as at the watershed level.

Albedo modeling

The SNICAR-modeled spectral albedo varies greatly with cell concentration in the visible part of the spectrum (Fig. 7). As to be expected, increasing snow algae cell concentration is associated with a decrease in albedo within the visible portion of the spectrum, 350–750 nm. There is a much greater variation in the spectral albedo and cell concentration within the first survey on 2 July, as compared to the 30 July survey. Associated with the decrease in cell concentration between the survey dates, the average spectral albedo increased within the 205–750 nm range of the spectrum shifting from 0.56 ± 0.17 on 2 July to 0.82 ± 0.09 on 30 July (Fig. 7).

Fig. 7: Spectral profile of the snow algae samples as simulated by the SNICAR model.
figure 7

Solid lines represent the samples from July 2, 2021 survey and dashed lines represent the July 30, 2021 survey. The darker red the line, the higher concentration of snow algae cells were in the sample. The gray bars align with the ten MicaSense camera bands with the darker bars noted with arrows corresponding to the green (524–538 nm) and red (642–658 nm) bands that were used in the optimized red/green band analysis.

Full size image

Although not directly comparable, the spectrally integrated snow algae modeled albedo simulated with the SNICAR model resembles the shape of our directional surface reflectance measured by the MicaSense camera on both UAV survey dates (Fig. 8). On the 2 July survey, the shape of the snow algae curve is much more exaggerated in the SNICAR output than in the MicaSense data. The snow algae SNICAR curve shape on July 30, 2021 survey more closely matches that of the MicaSense, but with about a 0.25 vertical shift in reflectance (Fig. 8).

Fig. 8: MicaSense measured directional reflectance values.
figure 8

MicaSense measured directional reflectance values for the Snow Algae, Snow, and Other classes on (a) July 2, 2021 and b July 30, 2021. Lines represent the local polynomial regression with shaded areas representing the standard error.

Full size image

The 10-band spectral resolution of the MicaSense imagery provides insight into the observed directional reflectance shifts of the snowpack from July 2, 2021 to July 30, 2021 under similar observation conditions. The Snow Algae classification points tend to absorb more strongly in the 400–600 nm range and reflect more in the 600–900 nm range (Fig. 8). In comparison, the Snow class reflects more in the lower wavelengths and absorbs more in the higher wavelengths (Fig. 8). To capture the snow algae spectral signature, the mean directional reflectance of the 444–560 nm MicaSense bands was calculated and compared to the mean reflectance in the 650 to 842 nm bands. In both survey dates, the mean directional reflectance of the Snow Algae class increases from 0.56 ± 0.10 on 2 July and 0.55 ± 0.12 on 30 July in the lower wavelengths to 0.69 ± 0.10 on 2 July and 0.62 ± 0.07 on 30 July in the higher. Comparatively, the reflectance of Snow class decreases from 0.91 ± 0.11 on 2 July and 0.69 ± 0.09 on 30 July in the low wavelengths to 0.83 ± 0.10 on 2 July and 0.62 ± 0.08 on 30 July in the high wavelengths (Fig. 8). There is more separation between the Snow and the Snow Algae classes in the mean directional reflectance of the 650 to 842 nm bands in the 2 July survey than compared to the 30 July survey where they are both 0.62 (Fig. 8).

Radiative forcing

Modeling the radiative forcing of the snow algae revealed an average IRF of 236.56 ± 79.55 W m−2 on 2 July and 88.86 ± 40.96 W m−2 on July 30, 2021 (Table 2). Taking the MicaSense index estimated snow algae-covered area for each survey date and the IRF of the snow algae, we calculated an average daily RF by snow algae of 18,290 ± 6150 MJ on 2 July and 2671 ± 1231 MJ on 30 July. If 0.334 MJ melts 1 kg of snow66, and we assume a wet snow density of 600 kg m3 that is typical this time of year, then we calculate that the presence of snow algae caused 91 ± 31 m3 of snow to melt on 2 July and 13 ± 6 m3 of snow to melt on 30 July.

Table 2 Instantaneous radiative forcing of snow algae, dust, and other light-absorbing impurities across the crysophere from this study and previous studies.
Full size table

In comparison to previous studies of snow algae radiative forcing, the calculated IRF values of the first survey in this study were greater than twice those calculated in previous studies of red snow algae (Table 2). The maximum IRF reported for snow algae in previous studies was from green snow algae in Antarctica and was still less than the average IRF of the first survey in this study19. The average IRF of snow algae calculated in this study is more comparable to the IRF of dust in the European Alps and the San Juan Mountains, U.S.67,68.

Our RF calculations demonstrate the substantial impact that snow algae have on snowmelt in mid-latitude snow through their snow-darkening effects. The snow algae patches can produce an IRF of up to 360 W m−2 during peak bloom, far greater than previous observations of red snow algae at higher latitudes that measured up to 186 W m−2 in the Antarctic Peninsula19 and 88 W m−2 in Alaska69. This suggests that snow algae might play a larger role in snowmelt at mid-latitudes than at higher latitudes. The mid-latitude range may provide more optimal snow algae growing conditions due to the greater amounts of incoming solar radiation and longer melt season, though further studies should examine this trend with standardized sampling and analysis.

We calculated the RF of the snow algae over the coinciding survey day, indicating that over the month of July snow algae contributed between 91 and 13 m3 of snowmelt each day. Assuming the snow algae RF scales linearly over the 29 days between the surveys, we can integrate the snowmelt over each day and estimate that there was a total of 1508 ± 536 m3 of snow algae-caused snowmelt in the 0.1-km2 basin. DEM differencing reveals a total 8,000,000 m3 of snowmelt in the study area between the two surveys. While only 1% of the snow in the study area was covered by snow algae, the snow algae contributed to 0.02% of the total snowmelt in the basin. However, the timing of snowmelt and feedbacks between temperature (i.e., heatwaves which are becoming more common in the PNW), snow algal bloom conditions, and snowmelt should be furthered explored. More frequent flights capturing the bloom conditions between the two aerial surveys, such as what were conducted with the insufficient MAPIR camera, could help us fill in these gaps. In addition, partitioning of the feedbacks and relative contributions of snow algae vs. other impurities should also be explored. Given the importance of seasonal snowmelt to downstream ecosystems and water resources across the North Cascades, RF attributed to snow algae should be considered in future watershed melt models as they may also play a role in the timing of peak snowmelt, as demonstrated by large dust loadings in SW Colorado70.

Conclusion

Here, we show the ability to remotely detect and map snow algae in mid-latitudes using UAVs affixed with a multispectral camera and that remote detection optimization differs based on bloom conditions. The two approaches presented in this study show great potential to remotely detect snow algae with a UAV. The PCA approach demonstrates the potential of multispectral and hyperspectral imagery to detect snow algae by using the information in all bands of data instead of a select few. In addition, the PCA approach allowed us to explore which bands may be useful to incorporate in future indices that may not have been previously considered. These methods could be applied to current and upcoming high-resolution spaceborne imaging spectroscopy missions, along with hyperspectral sensors mounted on UAVs in order to advance our ability to detect subtle spectral features such as algae absorption in snow and feed them into physically based radiative transfer models.

Both the principal component thresholding and the optimized red/green band indexing approaches provide clear separation between snow algae and the rest of the snowpack for the first survey date. By the second survey date, both approaches had less distinct separation between the classes that complicated the classification. In addition, the variation in snow algae cell and pigment concentration enabled a comparison of the approaches for differing bloom intensities. However, given our finding that the success of the spectral classification depends heavily on algae bloom intensity, careful consideration should be given when applying algorithms to other UAV and satellite imagery, demonstrating the need for further coupled ground validation.

Combining our mapping results with the ground data and spectral modeling, we calculated the RF for snow algae and scaled this across the mapped snow algae extent to estimate the cumulative impact on RF. The RF implications of snow algae have been documented previously in Antarctica19 and Alaska69, yet this is the first study to evaluate snow algae RF in the mid-latitude range of the North Cascades. Our results demonstrate the potential to map snow algae and assess the RF over expansive areas of the cryosphere using UAV technology.

The influence of snow algae on albedo in the cryosphere is not currently included in water resource and global climate models. Considering the substantial RF of snow algae demonstrated in this study, there is need for the inclusion of snow algae in both. Being able to predict the temporal and spatial distribution of snow algae could improve melt models and better prepare us for a warming future. These findings are also applicable to other mountain ranges with seasonal snow that feed large basins and have snow algae present in spring and summer, such as the Colorado Rockies and Sierra Nevada in California, among others. Not much is known about the extent of snow algae habitat, however, snow algae blooms are expected to expand in their range with warming temperatures55. Being able to reliably map snow algae will greatly enhance our understanding of the RF implications caused by snow algae across the cryosphere.

Methods

Study site

The Bagley Lakes are a set of snowmelt-fed lakes located in the Mount Baker-Snoqualmie National Forest (48.854°N, −121.692°W, 1277 m above sea level). The basin area remains snow-covered until late in the summer and has consistent blooms of snow algae every year. Only red snow algae were present during our surveys. The basin covers an area of 0.1 km2 consisting of a flat marshy expanse and a large lake nestled in a cirque with steep rocky sides (Supplementary Fig. 3). Beneath the larger, upper lake is a connected, smaller lake that flows into a tributary of the North Fork Nooksack River. The snowmelt-fed Nooksack River supports agricultural irrigation, drinking water, and salmon spawning in the region. The local Tribes have relied on the snowmelt from the North Cascades to replenish the rivers and support the salmon populations since time immemorial. After the winter of 2020 to 2021, the area had an early May snowpack thickness of 4 m71. The basin receives its snowpack both through direct precipitation and through avalanching events from the steep walls of the basin that deliver snow and debris to the bottom of the basin. The area is reasonably accessible from the road during the late summer months, making it an ideal site for repeat surveys. A portion of the study area lies within an avalanche runout zone and, due to numerous winter avalanching events, there debris is common in this area. For this analysis, we chose to remove the avalanche runout zone from the study due to the steepness of the terrain, as well as the complex spectral properties that would make it unlikely to detect snow algae in the patches of snow within the debris-covered area. The lake area was also removed from the analysis since the photogrammetric process does not perform well over water surfaces and due to our inability to safely retrieve validation samples from the lake when it was snow-covered. This also focused our algorithms specifically on snow-covered with snow algae.

UAV surveys

Over the summer of 2021, three survey flights were completed of the Bagley Lakes basin. The first two surveys were flown to collect data during the snow algae bloom, one on July 2, 2021 and a second survey on July 30, 2021. The third survey was flown on September 24, 2021 to acquire snow-free ground elevation data. The UAV surveys were conducted with a DJI Matrice 210 equipped with a MicaSense Dual camera system comprised of both the RedEdge MX and the RedEdge MX Blue sensors. Combined, the two sensors acquire data in the coastal blue (430–458 nm), blue (459–491 nm), green-1 (524–538 nm), green-2 (546– 574 nm), red-1 (642–658 nm), red-2 (661–675 nm), red edge-1 (700–710 nm), red edge-2 (711–723 nm), red edge-3 (731–749 nm), and near-infrared (813–871 nm) (www.micasense.com). The cameras captured RAW (12-bit) images every 1 s. Prior to UAV takeoff, the MicaSense cameras captured an image of the MicaSense camera calibration panel to later convert the imagery to reflectance. The DJI Matrice 210 was flown as close to solar noon as possible where the entire survey area was under direct sunlight and the UAV followed a flight plan designed and executed in the DJI Pilot application. Flight altitude was set at 90 m above ground level with a camera view angle of 90° and the flight lines were set for a 75% sidelap and frontlap, generating a ground resolution of 7 cm/pixel. The survey area flight took ~20 min to complete with the Matrice and required 2–3 batteries.

Prior to takeoff for each survey flight, a base station was established with an EMLID Reach RS2 RTK GNSS receiver in RTK mode that logged its location for the entirety of the field survey and communicated GPS corrections to the rover GPS. The base was stationed where it would have a clear view of the sky and within range of the rover to transmit corrections. Targets to be used as Ground Control Points (GCPs) and Ground Validation Points (GVPs) were spread out over the entire survey area and the GPS coordinates were collected with the rover EMLID Reach RS2 RTK GNSS receiver in RTK mode.

Positional accuracy

The orthomosaic image accuracy is derived from the ground control points (GCP) and ground validation points (GVPs) for each survey date. On the 2 July survey, there were 12 total targets, 9 that were used as control points and 3 as validation points. This resulted in a total GVP accuracy of 0.37 m. The 30 July survey only used 10 targets instead of 12 that were divided into 7 GCPs and 3 GVPs. Even though there were only 7 GCPs to guide the SfM image processing, the total error of the GVPs was less than that of the 2 July survey at 20.69 cm. The orthomosaics produced in the SfM process generated pixel sizes of 7 cm. Since the XY error is less than the pixel size, 1.64 cm for 2 July and 1.36 cm for 30 July, the lateral error is encompassed within the pixel.

Sample collection

Immediately following each UAV survey, snow algae blooms were identified, and the GPS coordinates were acquired with the rover GPS receiver. For the first survey, 9 snow algae samples were collected (n = 9) and 10 snow algae samples were collected for the second survey (n = 10), both with associated GPS data. The samples were collected opportunistically across a range of bloom intensities dispersed throughout the study area to characterize the heterogeneity of the snow algae bloom intensities and to test the proposed indices. In addition to the coincident samples collected after the flights on 2 July and 30, additional ground samples were collected on June 26 (n = 5), July 9 (n = 10), July 18 (n = 10), and July 22 (n = 10) 2021. The samples were collected from a 0.1 m by 0.1 m square to a depth of 2 cm using a metal spatula and stored in WhirlPak bags. The samples were immediately placed in a black trash bag to prevent continued photosynthesis. The coordinates of the center of the sample locations were acquired and the samples were melted in the dark at ambient temperature over 24 h. The next day, the melted samples were processed for pigment composition (three replicate filters), cell concentration, and ash-free dry mass (AFDM). The pigment and ash-free dry mass samples were vacuum filtered through pre-ashed and -weighed 0.47-µm glass fiber filters. The filters were folded in half, wrapped in tin foil, and stored at −20 °C until they were ready to be processed for further analysis.

Pigment analysis

The frozen filters (triplicates per homogenized sample) were processed for pigment composition by high-pressure liquid chromatography (HPLC) and quadrupole time of flight (QTOF) following a modified method based on Remias and Lutz72. The filters were shock frozen in liquid nitrogen for 10 min before being ground in a mortar and pestle until disintegrated. The ground filter was then resuspended in 1 mL dimethylformamide (DMF) and shaken on a vortex mixer for 10 min with 1.0-mm glass beads. The samples were centrifuged for 10 min at 4700 rpm and the resulting supernatant was extracted and filtered through a 0.2-µm Nylon filter. The internal standard, trans-beta-Apo-8’-carotenal, was added to the sample to a constant concentration.

The samples were stored frozen and run within 24 h of processing on an Agilent HPLC 1290 Infinity 2 with an Agilent Advanced Bio 6545XT LC-QTOF and Agilent Eclipse + C18 RRHD. Two running buffers were used: buffer A was isopropyl alcohol:acetonitrile at 90:10 with 0.1% formic acid, and buffer B was nanopure water:acetonitrile at 60:40 with 0.1% formic acid. The buffers were run at a flow rate of 0.35 mL/min with the following gradient: from 0 to 12 min 80% buffer A and 20% buffer B, from 12 to 15 min 30% buffer A and 70% buffer B, from 15 to 15.1 min 5% buffer A and 95% buffer B, from 15.1 to 20 min 80% buffer A and 20% buffer B.

QTOF force settings were set to 2000 nozzle voltage, 3500 capillary voltage, sheath gas flow rate of 11 L/min, sheath gas temperature of 350 °C, drying gas flow rate of 10 L/min, drying gas temperature of 325 °C. The charge source was a Dual AJS ESI. We acquired MS/MS with a collision energy of 6.5 (m/z)/100 + 2. Samples were quantified for the following pigments and their associated masses: beta carotene 536.87 g/mol, lutein 568.43 g/mol, astaxanthin 597.39 g/mol, chlorophyll a 893.54 g/mol, and chlorophyll b 929.50 g/mol. The presence of the pigments were verified by the presence of qualifying mass fragments averaged across replicates. Measurements were taken from distinct samples with replicates used for error quantification.

Cell concentration

Melted snow algae samples were homogenized and a 12 mL aliquot of whole sample was collected. The sample was fixed with 50% glutaraldehyde to 2% final concentration and stored at −20 °C. The samples were run on a Guava easyCyte flow cytometer equipped with a blue laser. Bead checks were conducted at the beginning of each sample run to ensure consistency between run dates. Four replicates of 200 µL were analyzed for total of 800 µL per sample and cells with high red fluorescence and high forward scatter were counted as red snow algae as validated by optical microscopy. The sample was run through the Guava for 120 s. Measurements were taken from distinct samples with each separate sample measured four times in 200 µL increments.

Ash-free dry mass

The AFDM of the snow algae samples was determined following the method described in ref. 73. The frozen filters were dried for 1 h at 104 °C and left to cool in a desiccator before being weighed. The filters were then combusted in the muffle furnace for 1 h at 550 °C and left to cool in the desiccator before being weighed again. The AFDM (Supplementary Table 2) was calculated as the difference between the weight of the filter before and after combustion divided by the sample volume to reveal the AFDM as mg organic material per liter. Measurements were taken from distinct samples with replicates used for error quantification.

Image processing

The flight images were refined to only include in-survey images. Images acquired during takeoff, landing, and transit were removed from the analysis. The MicaSense images were calibrated to surface reflectance based on the calibration panel images within the Agisoft Metashape Professional Software version 1.6.4 following MicaSense processing procedures74. The altitude of the images was adjusted with an open-source Python script to represent meters above sea level based on the launch elevation and the flight altitude (Agisoft LLC, 2017)75. The Image Quality, or sharpness, was estimated and blurry images, defined as those with a quality below 0.5, were removed from the rest of the analysis to improve photogrammetric processing following the guidelines in the Agisoft Metashape Professional Edition User Manual (Agisoft LLC, 2021)76. Each set of images was aligned separately with high accuracy, 500,000 key point limit, and 0 tie point limit. All points with a reprojection error greater than 0.5 pixels and any obvious outlier points were removed from the sparse point cloud.

The GCPs and GVPs were manually marked in at least 6 images where the target was most visible. The target coordinates were loaded, and the sparse point cloud was updated. The camera positions were optimized before constructing the dense point cloud with high-quality and aggressive depth filtering. The DEM was constructed based on the dense point cloud and the orthomosaic was constructed based on the DEM.

Principal components analysis

The first approach reduced the 10-band MicaSense orthomosaic to three principal component bands using the Forward PCA Rotation New Statistics & Rotate tool in the ENVI software version 5.6. The efficacy of the ordination was tested with 100 randomly located points that were visually classified as Snow or Other (i.e., rock, water, vegetation) in addition to the in situ snow algae locations that were collected in the field. These points were plotted with respect to the first three principal component bands to assess which combination of bands best distinguishes the snow algae from clean snow. Threshold limits were developed based on these points and applied to the entire orthomosaic.

Red/green band index development

In the second approach, snow algae detection indices were developed based off the MicaSense bands for both the 2 July and the July 30, 2021 surveys. The 100 training points for each day were used to assess the efficacy of the indices and determine which MicaSense band pairs separated the snow algae from the snow the best for both survey dates. The developed indices were chosen based on literature-proposed indices and based on the bands that most closely aligned to those of satellite platforms.

Albedo and radiative forcing

The albedo of the snow algae-covered snow and clean snow was simulated using the SNICAR model. The snow algae cell and pigment concentration data collected in the field was used as input parameters for the model. Our pigment dry mass fractions, calculated following Flanner et al.61, fell below the allowable range set by the SNICAR model. To fit SNICAR model parameters, the percent dry mass of each pigment was calculated and that percent was taken of 0.05 µg to fit in the allowable dry mass range. To fit SNICAR model parameters, pigment dry mass fractions were scaled by a factor of 0.05 to reach the low end of the allowable range. Thus far, the SNICAR parameterization has been established based on the primary glacier algae pigments, as opposed to snow algae pigments. After conversations with the developer of SNICAR, we decided to scale up our values to be usable inputs to the current version of the model (M. Flanner, Personal communication, February 25, 2022). Average cell diameter was set at 15 µm based on optical microscope analyses conducted on snow algae samples collected from the same study area the previous year. The solar zenith angle was estimated as 40° on 2 July and 43–50° on 30 July using the NOAA Solar Zenith Calculator based on the sample date, time, and location (gml.noaa.gov/grad/solcalc/azel.html) (Supplementary Table 3). Snowpack thickness was derived by differencing the UAV generated DEMs from the snow algae UAV surveys with the snow-free UAV survey and snowpack density was set to 600 kg/m3 based on the literature values for the North Cascades region77. The snow grain shape was set to spheroids, snow grain effective radius was set to 100 µm and the rest of the parameters, including the dust parameters, were set as the default or zero. As the output of this model, the broadband snow albedo is calculated as well as the albedo and fraction of incident irradiance within each spectral band at 10 nm spaced intervals. Clean snow albedo was modeled using the same parameters described above, but without any of the snow algae inputs.

The instantaneous radiative forcing (IRF) of the snow algae was calculated as:

$${IRF}\approx \mathop{\sum }\limits_{350}^{850}{E}_{d}\left(\lambda \right)\left({R}_{{clean}}\left(\lambda \right)-{R}_{{algae}}\left(\lambda \right)\right)\triangle \lambda$$
(3)

where IRF is calculated as the sum of the wavelength-specific downward flux, Ed, multiplied by the wavelength interval and the difference between clean snow reflectance, Rclean, and snow algae-covered snow reflectance, Ralgae, over the wavelengths 350–850 nm. The downward flux was retrieved from https://www.pvlighthouse.com.au (last access: February 2, 2022)19,69.