Prior to the acquisition of the UAV imagery, hemispherical conical reflectance factor (HCRF) measurements were collected on snowfields in Hardangervidda, Southern Norway, between 8 July and 6 August 2023. The data are presented in Chevrollier et al. (2024). The measurements were collected using an ASD FieldSpec 4 spectroradiometer (spectral range 0.35–2.5 μ m) and a black tripod, following the methodology of Cook et al. (2017). Each measurement was acquired with the bare fiber (field of view of 25 ° ) in 10 replicates, immediately (<10 s) after a reference spectrum was measured using a calibrated Spectralon panel. All measurements were performed at nadir view (viewing zenith angle θ = 0 ° ) with the tripod oriented toward the Sun to avoid shadow effects. Most (80%) of the measurements were taken for a solar zenith angle between 38 ° and 50 ° , and the remaining (20%) were taken for a solar zenith angle below 59 ° . The spectra were corrected for the panel spectral response and the step at 1 μ m caused by the misalignment of the SWIR and NIR sensors in the ASD FieldSpec instrument, following Painter (2011). The step at 1 μ m was, however, not visible on any spectrum, probably because the instrument was warmed up for at least an hour prior to the measurements, and the selected surfaces were homogeneous enough. Five spectra were removed from the analysis because clear calibration errors were detected.

A DJI Mini 3 camera was then deployed to capture RGB imagery of the snowfields. The Mini series are the cheapest and lightest ( < 250 g) of the DJI UAVs and are, therefore, most readily affordable for integration into monitoring programs. The integrated RGB camera has a 12-megapixel sensor, with a 82.1 ° field of view. The imagery used in this study was directly saved in JPEG format, and no further processing was applied before analysis. The Spectralon calibration panel used for the HCRF measurements was placed on the snowfield for each image so that each image was normalized to the incoming irradiance to retrieve reflectance. All images were acquired with the camera at nadir view, and the UAV in the stationary mode directly above the middle of the snowfield. The images of the entire snowfields had a ground resolution of 1.3–3.0 cm, while the images acquired to analyze smaller scale surface patterns had a ground resolution of 0.4–0.7 cm.

A tree-based gradient boosting binary classifier was built for binary snow detection using the scikit-learn machine learning python package (Pedregosa et al., 2011). The model takes the UAV R, G, and B reflectance as input and returns whether the pixel is snow or not. Surface areas were manually labeled on 27 different UAV images, representing different surface types and resolutions, with half of the data points labeled as snow and the other half labeled as non-snow (surrounding vegetation and rocks), to ensure an equal representation of both classes during training and testing, for a total of 1.6 10 7 data points. In addition, 80% of these data points were fed into the model for training, whilst the remaining 20% was used for testing. The model was trained with a hold-out validation set of 10% of the training dataset and was set to stop the training after the error during validation did not change for five consecutive epochs. The hyperparameters of the classifier were tuned to minimize the error (percentage wrongly classified), and the validation set served as a control to avoid overfitting. The learning rate was set to 0.5, the booster set to “exact”, and the other parameters were kept as the scikit-learn default values (scikit-learn version 1.5.1). The error on the training and testing set was similar and very low (0.02%), with as many false positives as false negatives, indicating that the model performed very well in delineating snow surfaces from the UAV imagery despite the complex surface features linked to the presence of light-absorbing particles. A watershed algorithm using scikit-image (van der Walt et al., 2014) was finally applied to remove the pixels wrongly identified as snow on land. For the remaining features such as branches, manual masks were outlined and applied to the images using the OpenCV-Python package (Bradski, 2000).

The circular surfaces associated with the HCRF ground measurements were visually identified on the higher-resolution UAV imagery using a series of pictures acquired during the acquisition of ground measurements and manually delineated on ImageJ (Schneider et al., 2012) using the footprint size of each measurement, calculated from the distance between the fore-optic and the snow surface. The R, G, and B reflectance of the UAV imagery corresponding to each HCRF ground measurement were then extracted, and the ground measurements were rescaled to R, G, and B bands using the spectral response of DJI Phantom Pro 4 (Burggraaff et al., 2019) as no spectral response could be provided by the manufacturer for DJI mini 3. Ground red-to-green (RG) band ratio indices were then computed from ground HCRFs, and the associated airborne RG band ratio indices were computed from the average RG index of the circles on the UAV imagery. This ratio was chosen because it has been used in the past with airborne (e.g., Healy and Khan, 2023) and spaceborne (Gray et al., 2021; Ganey et al., 2017; Engstrom et al., 2022) data to detect the presence of red algae, motivated by a strong dip in the surface reflectance between the chlorophyll and carotenoid pigment absorption features of the algae. This analysis relied on the Python packages scikit-image (van der Walt et al., 2014), pandas (Pandas development team, 2024; Wes, 2010), Matplotlib (Hunter, 2007), SciPy (Virtanen et al., 2020), PIL (Umesh, 2012), GeoPandas (den Bossche et al., 2024), Rasterio (Gillies et al., 2013), and NumPy (Harris et al., 2020).

An ordinary least square regression was fitted between the ground and airborne indices, showing a strong correlation between the 129 UAV and ground RG indices (r = 0.96, Figure 1). This confirms that the RG indices captured by the UAV can closely reproduce ground RG indices. The bias and the deviation from the 1:1 line could be explained by several factors: 1) the RGB spectral response used to calculate the ground RG indices may not be exactly the same as that of the sensors onboard DJI mini 3; 2) the distribution of data in the UAV pixels for each circle is unlikely normally distributed; hence, the mean may not be a good predictor; 3) the spatial response of the FieldSpec sensor is not perfectly homogeneous, i.e., the pixels in the middle of the footprint may be more important; 4) the JPEG compression and lack of further processing may have introduced small biases in the airborne reflectance data; and 5) small changes in atmospheric conditions including cloud cover may have occurred between the acquisition of the ground and airborne measurements. These factors are not investigated here as the correlation from Figure 1 serves only as a validation of the UAV signal and is not used as a statistical model in the analyses.

For each of the 129 surfaces pairing ground and UAV RG indices, the light-absorbing particle properties were available from Chevrollier et al. (2024). Specifically, Chevrollier et al. (2024) retrieved the red algae, dust, and black carbon abundance, as well as their associated broadband albedo-reducing effect by inversion of hyperspectral reflectance using deep learning and radiative transfer modeling. These properties are, therefore, not direct measurements at the surface but retrievals from model inversions. This physics-based approach accounts for the effect of the changing snow grain size on the albedo-reducing effect caused by the absorbing particles. A positive BBA reduction corresponds to a darkening effect and is calculated by differentiating the BBA of the surface with and without algae. Causal statistical models were then built between each of the surface variables and the RG indices of the UAV imagery using Bayesian inference with Markov Chain Monte Carlo (MCMC) sampling, implemented with PyMC (Abril-Pla et al., 2023). This method enables the retrieval of probability distributions rather than single values as statistical model parameters, providing estimates of the uncertainties in the model predictions. Probabilistic inference requires the definition of prior distributions for the model parameters and the distribution of the likelihood function. Using Bayes’ theorem, posterior distributions are then computed from the prior and likelihood distributions, providing a refined estimate of the model parameters. Posterior distributions are finally sampled using MCMC to generate probability distributions for the model parameters that preserve the full statistical information contained in the posterior distributions. In this case, the algal concentration was fitted to the RG ratio with a second-order polynomial of the form y = a x 2 + b x + c , and the algal albedo reduction effect was fitted to the RG ratio with a linear model of the form y = a x + b , similar to other studies (e.g., Engstrom et al., 2022). The prior distributions on the slope and intercept for the algal albedo reduction model were normal, with a mean of 0 and a variance of 1. The prior distributions on the parameters a , b , and c in the algal abundance model were represented using Gaussian distributions with means of 100, −150, and 20, respectively, and a variance of 2. In both cases, the likelihood to maximize was represented by a Student-T distribution with ν = 3 degree of freedom, following standard practice for robust inference (Gelman et al., 1995). The variance of likelihood is necessarily positive; hence, we chose a Half-Cauchy distribution prior. Each posterior distribution was sampled 1,000 times with 2,000 iterations over four Markov chains. The median of the sampled posteriors was used to model algal abundance and albedo reduction from the UAV RG ratios, and the 95% credible interval (high-density interval with α = 0.95) was computed on the samples of the posteriors in order to estimate the uncertainty of the model predictions, plotted as shaded envelopes in Figure 2. The model validity domains were constrained to observed RG ratios producing physical (non-negative) algal properties, i.e., values between 0.099 and 0.1404 for the albedo reduction model, and to RG ratios between 0.0995 and 0.1404 for the algal abundance model. The pixels on the imagery with a RG ratio value outside these intervals could have been masked out, but the main effect would have been the removal of areas with clean snow and, therefore, an overestimation of algal coverage and impact. Before applying a given model, the pixels on the UAV imagery with RG ratios below or above the model validity domain were therefore set to the domain boundaries. The 2.9% of the pixels with a reflectance higher than 1 on at least one of the three bands was set to 1. The difference in the results between correcting the reflectance to 1 and removing the pixels was negligible (<100 cells mL − 1 absolute difference on the mean algal abundance and < 1 × 10 ⁻⁴ absolute difference on the mean algal albedo-reducing effect).

The UAV RG indices correlated very well with the algal abundance (r = 0.93) and the algal-driven albedo reducing effect (r = 0.83) but did not correlate with the dust and black carbon properties (r < 0.14, Figure 2). The RG ratio was used to map dust concentrations in the Alps deposited from a Sahara dust event (Di Mauro et al., 2015) as the Sahara dust abundance correlated with the ratio. The lack of correlation in our study is probably due to relatively lower concentrations of dust on the snowfields in the studied region and most likely due to different optical properties between the dust from the Sahara and the dust from the Hardangervidda region.

In order to calculate the algal properties on the visible bloom areas solely, surface pictures taken before each ground measurement were used to select the RG ratios associated with a visible bloom and determine a conservative threshold of 1.029. Specifically, the threshold was selected so that no bloom was visible on any surface picture associated with a lower RG ratio than this threshold.

The daily radiative forcing (RF; W m − 2 ) caused by the algae was then calculated by multiplying the BBA reducing with the daily (24 h averaged) shortwave incoming radiation measured with a four-component radiometer (CNR4, Kipp & Zonen, Netherlands) at the local weather station (Pirk et al., 2023). The algal-driven melt ( M algae , kg m − 2 ) was calculated by multiplying the RF with the energy of fusion of ice ( Δ H f = 334 kJ kg − 1 ). M algae was subsequently multiplied by the area to derive the melt in kg of snowmelt or divided by the snow density (here assumed to be 600 kg m − 3 ) to obtain the melt in meter equivalent.

The algal abundance maps were all rescaled at 3 cm resolution and numerically sampled pseudo-randomly to evaluate two ground sampling strategies in terms of how representative of the entire snowfield area ground samples of algal counts can be. The first strategy consisted in taking several small samples with a fixed footprint but distributed spatially across the snowfield, and the second consisted in increasing the sample footprint from a fixed location. In the first case, mean algal concentrations were calculated for 1,000 random points with an increasing sampling footprint from 0.001 to 4 m 2 . Hence, for each footprint size, 1,000 concentrations were calculated. In the second case, mean algal concentrations were calculated by repeating 1,000 times a random sampling of up to 30 samples with a fixed footprint of 12 × 12 cm. Hence, for each total number of samples, 1,000 concentrations were calculated. Then, the probability of obtaining a representative average concentration for a given sampling strategy (for example, 1 m 2 samples) was computed by comparing the 1,000 concentrations to the mean concentration of the entire patch with a given error. Hence, when half of the 1,000 random points were within the threshold error, the probability of a representative sampling was 0.5 for this given strategy. Here, different threshold errors were tested around 25%, which is in the range of error reported for the Neubauer hemocytometers (https://www.emsdiasum.com/docs/technical/datasheet/68052-14) typically used for field-measured cell counts.

The snow algal blooms were intense enough to reduce the overall albedo of each snowfield surface ( Δ BBA total ; Figures 3, 4). The albedo-reducing effect of the blooms was on average 0.012 ± 0.005 over the ∼ 9,740 m 2 of the snow-covered area analyzed and increased up to 0.028 ± 0.004 at the scale of a 710- m 2 snowfield (Figure 4). The darkening effect of red snow algae is often reported on the part of the snow surface that is covered with a visible bloom only, which represented 13% of the total area (see Methods). For each snowfield, the algal albedo-reducing effect of the algae on visible bloom areas was significantly higher than the estimates computed on the entire snow-covered area ( Δ BBA bloom ; Figures 3, 4), in particular for the snowfields with a main localized bloom. On average, Δ BBA bloom amounted to 0.045 ± 0.003 over the entire snow-covered area analyzed. This estimate is in good agreement with the albedo-reducing effects of 0.032 ± 0.006 and 0.045 ± 0.009 calculated for two glaciers in British Columbia in summer 2020 (Engstrom et al., 2022). The small difference can be explained by a combination of factors. For example, the algal concentration, distribution, and pigmentation differ between sites and time of imaging, and the presence of other LAPs varies between sites, resulting in different algal albedo-reducing effects for a given LAP concentration (Skiles and Painter, 2018; Kaspari et al., 2014; Chevrollier et al., 2024).

Given the large temporal variability in the incoming shortwave irradiance in the studied region (Supplementary Figure S1), the impact of red algae on snowmelt varied greatly during our observational period. On a cloudy day, the radiative forcing of the algae was negligible for all snowfields, whilst on a sunny day, the daily radiative forcing was 4.3 ± 1.8 W m − 2 over the total imaged area and reached up to 10 ± 1 W m − 2 on the snowfield where the algal albedo-reducing effect was the highest. Overall, the illumination conditions were mostly cloudy over the studied period, and the average algal radiative forcing was 2.2 ± 0.9 W m − 2 , causing 5,500 ± 2,300 kg of daily algal-induced snowmelt. This corresponds to ∼ 3 mm-equivalent of daily algal induced-melt within the total bloom area, amounting to ∼ 9 cm-equivalent during the entire observational period if the effect of the algae remained constant.

The coefficient of variation of the calculated algal concentration (standard deviation divided by the mean) ranged from 100% to 400%. The total coverage and intensity of the blooms therefore varied between and within the snowfields, confirming the high heterogeneity of snow algal blooms documented in other studies (Thomas, 1972; Thomas and Duval, 1995; Takeuchi et al., 2006; Lutz et al., 2017; Van Hees et al., 2023). When deep suncups were present on the snowfields, which are bowl-like structures forming spontaneously on the snow during the melting season in temperate alpine areas (Post and LaChapelle, 2000), the algae appeared to accumulate within them (Figures 5A, B). Differences in algal accumulation within suncups were also noticeable, and the algae could either be present in the entire depression (Figure 5A), mostly in the bottom of the depression (Figure 5B), or in the channel-like features where suncups are connected to each other (Figure 5C). In cases where the surface appeared smoother, yet still featuring small-scale ridges and depressions, algal distribution did not seem to be linked to surface roughness (Figure 5D).

In order to further characterize bloom heterogeneity and investigate the potential implications for the representativeness of surface samples collected on the ground, numerical sampling experiments were designed. Specifically, the pixels from the UAV imagery were extracted in two different ways to represent different ground sampling strategies. First, the pixels were extracted from an increasingly large footprint from a single fixed point, and second, the pixels were extracted from an increasing number of 12 × 12-cm footprint squares distributed across the snowfield (see Methods). For all snowfields, the second strategy was more efficient at capturing a representative average algal concentration of the entire snowfield (Figure 6). For a total area of extracted pixels 10 times lower, the probability of obtaining a representative concentration was at least twice higher with the second strategy, and the representativeness of the samples barely improved with increasing the footprint size from a fixed point. This was verified even for the snowfields with a main localized bloom (4 and 6; Figure 6). These results indicate that the scale of the heterogeneity of the snowfields always superseded the maximum sampled area tested of 4 m 2 . Representativeness was defined as being within a 25% error from the average concentration of the entire snowfield (see Methods), but the second strategy remained more efficient at capturing a representative concentration even for a different representativeness threshold (Supplementary Figure S2). The probability of representativeness, however, decreased with a stricter threshold, highlighting the limitations of using a discrete number of samples in the field to evaluate the concentration of an entire snowfield with precision. In addition, sampling various random points at the surface of a snowfield manually without destroying its natural structure is, at best, challenging, and at worst, impossible. The results hence emphasize the utility of UAV imagery to study algal blooms non-destructively, in particular in the context of monitoring, when sampling through time is required. However, thorough validation of the abundance retrievals from the UAV imagery against field-measured cell counts is still missing, in order to validate the absolute values of the retrieved concentrations.

Elucidating the processes behind algal spatial distribution was beyond the scope of this study, but we demonstrated the potential of a light-weight, low-cost, and easy-to-deploy UAV to map and describe algal surface distribution that can be used to monitor the blooms at high temporal and spatial resolution. Such datasets would enable the monitoring of large areas that cannot be accessed, circumventing the representativeness bias associated with ground sampling, and helping to elucidate major questions in the field of bio-albedo, such as the role of meteorological and snow conditions on the spatial distribution of algal blooms (Roussel et al., 2024; Chen et al., 2023).