1. Introduction
The Greenland Ice Sheet (GrIS) has experienced a significant surface darkening over the past decades (Dumont and others, Reference Dumont2014; Tedescoand others, Reference Tedesco, Doherty, Fettweis, Alexander, Jeyaratnam and Stroeve2016) owing to a reduction in surface albedo, which affects its mass balance (Wientjes and Oerlemans, Reference Wientjes and Oerlemans2010; Saito and others, Reference Saito, Sugiyama, Tsutaki and Sawagaki2016; Cook and others, Reference Cook2020). This lowering albedo, commonly known as darkening phenomenon, has particularly affected south-west Greenland over the so called ‘dark zone’ (Ryan and others, Reference Ryan2018), and it has been ascribed to a combination of factors. For example, the albedo reduction has been linked to less frequent snowfalls (caused by atmospheric blockings) and increased solar radiation, which foster the aging of snow grains (Lewis and others, Reference Lewis2021). Moreover, the accumulation of different light-absorbing particles (black carbon, mineral dust, volcanic ashes, Dumont and others, Reference Dumont2014) and biological activities promoted the reduction in albedo (Cook and others, Reference Cook2020; Wang and others, Reference Wang, Tedesco, Alexander, Xu and Fettweis2020). The dark ice extent further correlates with solar radiation (Shimada and others, Reference Shimada, Takeuchi and Aoki2016). The extent of dark ice is not only controlled by the abundance of impurities, but also by changes in the surface structures of the bare ice surface, such as cryoconite holes, which are water-filled small pits hiding the impurities into the ice (Shimada and others, Reference Shimada, Takeuchi and Aoki2016).
Over ice, different species of glacier algae reduce surface albedo (Chevrollier and others, Reference Chevrollier, Cook, Halbach, Jakobsen, Benning, Anesio and Tranter2023; Feng and others, Reference Feng, Cook, Naegeli, Anesio, Benning and Tranter2024) thanks to the presence of dark pigments, i.e. phenolic purpurogallin (purpurogallin carboxylic acid-6-O-b-D-glucopyranoside), which strongly absorbs solar radiation (Remias and others, Reference Remias, Holzinger, Aigner and Lütz2012; Yallop and others, Reference Yallop2012; Williamson and others, Reference Williamson, Cameron, Cook, Zarsky, Stibal and Edwards2019; Halbach and others, Reference Halbach2022), allowing the algae to bear enhanced radiation level (Williamson and others, Reference Williamson2020). Their growth is favoured by water and nutrients made available by snow and ice melting, leading to a positive feedback (the melt-albedo feedback) of increasing temperature and melting promoting algal blooming (Box and others, Reference Box, Fettweis, Stroeve, Tedesco, Hall and Steffen2012; Cook and others, Reference Cook2020; Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023; Halbach and others, Reference Halbach2025). In this context, the formation of a surficial porous layer of white ice, called weathering crust (Woods and Hewitt, Reference Woods and Hewitt2023; Traversa and Di Mauro, Reference Traversa and Di Mauro2024) retains water and sediments, providing an ideal habitat for algae growth (Cooper and others, Reference Cooper2018; Takeuchi and others, Reference Takeuchi, Sakaki, Uetake, Nagatsuka, Shimada, Niwano and Aoki2018; Tedstone and others, Reference Tedstone2020; Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023). Additionally, specific cyanobacteria facilitate the creation of organic matter, leading to the formation of dark-coloured aggregates on the ice, commonly known as cryoconite (Uetake and others, Reference Uetake, Tanaka, Segawa, Takeuchi, Nagatsuka, Motoyama and Aoki2016; Takeuchi and others, Reference Takeuchi, Sakaki, Uetake, Nagatsuka, Shimada, Niwano and Aoki2018; Traversa and others, Reference Traversa, Scipinotti, Pierattini, Fasani and Di Mauro2024; Dory and others, Reference Dory2025). The porosity of this peculiar kind of ice and the dynamics of cryoconite holes also influence the spatial distribution of glacier algae (Takeuchi and others, Reference Takeuchi, Sakaki, Uetake, Nagatsuka, Shimada, Niwano and Aoki2018). For these reasons, previous studies found that glacier algae have a stronger role in glacier darkening than mineral particles, making the biological influence on surface ice melting of particular interest in glaciology (Stibal and others, Reference Stibal2017; Cook and others, Reference Cook2020; Chevrollier and others, Reference Chevrollier, Cook, Halbach, Jakobsen, Benning, Anesio and Tranter2023).
On the GrIS, algal communities are dominated by species adapted to extreme environments, such as green algae from the Zygnematales order, with evidence of the presence of the Ancylonema genre (i.e. A. nordenskioldii, a filamentous species A. alaskana, a unicellular species) or the Chlamydomonadales order (Sanguina nivaloides), and cyanobacteria such as Phormidesmis priestleyi and Chroococcaceae cyanobacterium (Onuma and others, Reference Onuma, Takeuchi, Tanaka, Nagatsuka, Niwano and Aoki2018; Williamson and others, Reference Williamson, Cameron, Cook, Zarsky, Stibal and Edwards2019; Di Mauro and others, Reference Di Mauro, Garzonio and Baccolo2020a; Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023). Among these species, the most abundant over glaciers is A. nordenskioldii (Takeuchi, Reference Takeuchi2015; Stibal and others, Reference Stibal2017; Lutz and others, Reference Lutz, McCutcheon, McQuaid and Benning2018; Takeuchi and others, Reference Takeuchi, Tanaka, Konno, Irvine-Fynn, Rassner and Edwards2019). All these algae and bacteria species present a strong growth during the glacier-melting season, dependent on the length of the melting period and available mineral dust (Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023). Their physiology includes mechanisms of protection from UV rays and optimisation of photosynthesis in low-light environments (Williamson and others, Reference Williamson, Cameron, Cook, Zarsky, Stibal and Edwards2019; Hoham and Remias, Reference Hoham and Remias2020; Di Mauro and others, Reference Di Mauro, Garzonio and Baccolo2020a).
In order to analyse the glacier algae-driven darkening over the GrIS, different approaches were carried out, especially by means of remote-sensing, e.g. by creating indices or ratios based on multi-spectral satellite data. Wang and others, (Reference Wang, Tedesco, Xu and Alexander2018, Reference Wang, Tedesco, Alexander, Xu and Fettweis2020) analysed the spatiotemporal variability of glacier algae in Greenland at the ice-sheet scale, by using data from Sentinel-3 OLCI launched by the European Space Agency (ESA) and by the MERIS spectrometer onboard ENVISAT by applying a ratio of bands 11 (709 ± 10 nm) and 9 (674 ± 8 nm) and bands 9 (709 ± 10 nm) and 7 (665 ± 10 nm), respectively, both at 300 m spatial resolution. These ratios took advantage of the radiation absorption feature located at 680 nm, usually linked to the presence of Chlorophyll-a (Takeuchi, Reference Takeuchi2002; Remias and others, Reference Remias, Holzinger, Aigner and Lütz2012). In other regions of Earth, similar approaches have been applied to retrieve algal abundance from satellite observations. For example, Takeuchi and others (Takeuchi and others, Reference Takeuchi, Dial, Kohshima, Segawa and Uetake2006) exploited SPOT-2 satellite data (20 m spatial resolution) to estimate algae abundance over snow in Alaska by applying the ratio of the red band (610–680 nm) and the green band (500–590 nm). Additionally, Di Mauro and others (Reference Di Mauro, Garzonio and Baccolo2020a) developed another spectral ratio for the European Alps based on field spectroscopy data. They proposed a reflectance ratio using Sentinel-2 band 6 (740 ± 40 nm) and band 4 (665 ± 30 nm). Other approaches included a supervised classification (random forest) on Uncrewed-Automatic Vehicle (UAV) and Sentinel-2 (Cook and others, Reference Cook2020) data or spectral unmixing for discriminating among algae and other surface impurities (Williamson and others, Reference Williamson, Cameron, Cook, Zarsky, Stibal and Edwards2019; Wang and others, Reference Wang, Tedesco, Alexander, Xu and Fettweis2020; Di Mauro and others, Reference Di Mauro, Garzonio and Baccolo2020a, Reference Di Mauro2024a; Engstrom and Quarmby, Reference Engstrom and Quarmby2023; Roussel and others, Reference Roussel2024).
An area outside of the dark zone which received particular interest in previous research is the Qaanaaq Ice Cap, located in northwest Greenland, where the darkening phenomenon due to the presence of glacier algae has been previously described (Uetake and others, Reference Uetake, Naganuma, Hebsgaard, Kanda and Kohshima2010; Sugiyama and others, Reference Sugiyama, Sakakibara, Matsuno, Yamaguchi, Matoba and Aoki2014; Anesio, Reference Anesio2024). Until now, research over the northwest of the GrIS has been carried out by means of field observations (Aoki and others, Reference Aoki, Matoba, Uetake, Takeuchi and Motoyama2014) to analyse samples of algae and cryoconite (Onuma and others, Reference Onuma, Takeuchi, Tanaka, Nagatsuka, Niwano and Aoki2018; Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023; Onuma, Fujita, and others, Reference Onuma, Fujita, Takeuchi, Niwano and Aoki2023) and their effect on ice albedo (Takeuchi and others, Reference Takeuchi, Sakaki, Uetake, Nagatsuka, Shimada, Niwano and Aoki2018). Several field campaigns have been carried out in the area since 2012 (Aoki and others, Reference Aoki, Matoba, Uetake, Takeuchi and Motoyama2014; Tsutaki and others, Reference Tsutaki, Sugiyama, Sakakibara, Aoki and Niwano2017; Nishimura and others, Reference Nishimura2023), mostly over outlet glaciers flowing from the southern side of the ice cap (e.g. Qaanaaq Glacier). Nevertheless, despite the relevance of studying such an extreme environment (the ice cap is located at around 77°N), a spatially distributed analysis of algae abundance on the ice cap in its entirety has not yet been conducted, and high-resolution remote-sensing data in this region are still underutilised. One of the main open questions in glacier algae studies is what mechanisms influence algal growth and blooms on the Greenland glaciers (Di Mauro, Reference Di Mauro2020b; Halbach and others, Reference Halbach2023). Moreover, previous research has not studied in detail the phenology of glacier algae over the GrIS. The present research aims at filling these gaps by presenting insights from field and satellite observations (Sentinel-2 at 10 m spatial resolution) over the past decade by taking advantage of the Qaanaaq area as a test site. Thus, the phenology of glacier algae was analysed and compared with variations in meteorological and glaciological parameters (such as albedo, temperature and snowfall) over the period 2016–2023.
2. Study area: the Qaanaaq Ice Cap
The study site is the Qaanaaq Ice Cap, one of the ice caps of north-western Greenland located at 77°N–69°W and extending for 273 km2 (August 2023, Fig. 1a).
Figure 1. (a) overview of the Qaanaaq Ice Cap (blue star in the overview map of Greenland, Moon and others, Reference Moon, Fisher, Stafford and Thurber2023), Sentinel-2 image acquired on 18 August 2023 in the background and ice-cap outlines in light blue and ice sheds in black. (b) study sites from the 2014 and 2023 field campaigns over Qaanaaq Glacier (black rectangle in a). (c) Zoomed-in aerial view of 2014 study sites as of UAV acquisitions in August 2023. (d) collection of spectral measurements over the Qaanaaq Glacier in 2023. (e) Light micrograph of one of the 2023 samples where Ancylonema nordenskioldii cells are clearly visible.
The ice cap, located on the Piulip Nunaa peninsula, is detached from the GrIS by a few hundreds of meters on the north-east side and includes 11 ice sheds (RGI 7.0 Consortium, 2023). The village of Qaanaaq is located close to the ice cap, at a distance of about 2 km to the south-east. Starting from Qaanaaq, several field campaigns were carried out in the past decade, especially over the Bryant and Qaanaaq Glaciers and the adjacent outlet glaciers for investigating glacier mass balance and meteorological parameters by means of an automatic weather stations (SIGMA-B AWS at 944 m a.s.l.; Aoki and others, Reference Aoki, Matoba, Uetake, Takeuchi and Motoyama2014; Tsutaki and others, Reference Tsutaki, Sugiyama, Sakakibara, Aoki and Niwano2017; Nishimura and others, Reference Nishimura2023). Over the Qaanaaq Glacier (Fig. 1b), with an area of almost 10 km2 (August 2023), previous studies have identified the presence of several phototroph blooms and cryoconite (Takeuchi and others, Reference Takeuchi, Nagatsuka, Uetake and Shimada2014, Reference Takeuchi, Sakaki, Uetake, Nagatsuka, Shimada, Niwano and Aoki2018; Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023). This site is known to host one of the highest concentrations of glacier algae over the northern GrIS (Uetake and others, Reference Uetake, Naganuma, Hebsgaard, Kanda and Kohshima2010; Anesio, Reference Anesio2024), and thus has been selected as the test site for the present research. In addition to its glaciological relevance, this glacier is also easily reachable by feet from Qaanaaq village and it was surveyed in 2014 by a campaign conducted in the context of SIGMA project (Aoki and others, Reference Aoki, Matoba, Uetake, Takeuchi and Motoyama2014) and in 2023 within an INTERACT TA project (Fig. 1d).
3. Data and methods
3.1. Data
3.1.1. Field measurements at Qaanaaq Glacier (2014 and 2023)
Surface ice samples and field spectral measurements used in the present research were collected in two different melting seasons, i.e. summers of 2014 and 2023 (Fig. 1b).
In this study, we also used glacier ice samples collected over the Qaanaaq Glacier in 2014 (32 surface-ice samples collected between 20 July and 3 August) at four sites (S1, S2, S3 and S4 at 247, 441, 672 and 772 m a.s.l., respectively), as published in Onuma and others (Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023). As for 2023, sampling activities were performed on 7–9 August, collecting 21 ice surface samples over the Qaanaaq Glacier between the glacier terminus and an elevation of about 800 m a.s.l. (passing by S2, S3 and S4 sites; S1 was in the proglacial area by that time). Surficial ice (∼1 cm) was sampled using a stainless-steel spoon and preserved in Whirl-Pak® like plastic bags (following the 2014 campaign methodology, Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023). Subsequently, the samples, once melted at room temperature at the DMI Observatory, were poured into 50 mL falcons and immediately preserved in 1% Lugol’s iodine solution (solution of potassium iodide with iodine in water). The prepared samples were then shipped to Chiba University laboratories in order to be analysed.
Simultaneously, in 2014 several spectral measurements were acquired over each sample site by means of a VIS-NIR (350–1050 nm) spectrometer (MS-720, Eiko Seiki Co., Japan; Onuma, Fujita, and others, Reference Onuma, Fujita, Takeuchi, Niwano and Aoki2023). Field spectroscopy measurements were also acquired in 2023 over each sampling site by means of a field spectrometer, the Reflectance boX—RoX (Fig. 1d), covering the wavelength range of 400–865 nm with a spectral resolution of ∼0.75 nm (Naethe and others, Reference Naethe2024). For both field observations, broadband albedo (hereafter referred as albedo) was estimated by calculating the ratio between reflected and incident spectral radiance integrated in VIS-NIR wavelengths (400–865 nm, Traversa and others,, Reference Traversa, Scipinotti, Pierattini, Fasani and Di Mauro2024). For consistency with the 2023 measurements, albedo for measurements acquired in 2014 was estimated over the same wavelength range.
3.1.2. Remote-sensing dataset
Sentinel-2 (ESA, 2025) data (tile area: ∼12 000 km2) were used as the main source of satellite imagery. Images were downloaded from the Copernicus Browser portal for 2016–2018 (https://browser.dataspace.copernicus.eu/, last access: 21 May 2025) and via Google Earth Engine (GEE) platform for 2019–2023 (since GEE provides atmospherically corrected Sentinel-2 data for the study area only since 10 June 2018). Sentinel-2 acquires data in several spectral bands, including visible, near-infrared, short-wave infrared and thermal-infrared bands, with varying spatial resolution from 10 m to 60 m. We used atmospherically corrected data (Sentinel2-L2A product in Copernicus Browser) which provides bottom of atmosphere reflectance after application of atmospheric correction (Sen2Cor processor). In GEE, the corresponding Harmonized Sentinel-2 MSI Level-2A surface reflectance product was employed. In this product, pixel values are shifted in the same range as in scenes prior to the application of PROCESSING_BASELINE ‘04.00’ on 25 January 2022 (Traversa and Di Mauro, Reference Traversa and Di Mauro2024).
We considered only images acquired during the summer months (1 June–30 September) and excluded from our dataset all images with cloud cover >30% of land surface. Additionally, clouds were evaluated on a pixel basis and cloudy pixels were excluded based on the scene classification layer (https://sentiwiki.copernicus.eu/web/s2-processing, last visit: 16 May 2025). A total of 341 images were analysed, distributed as follows: 34 in 2016, 40 in 2017, 37 in 2018, 47 in 2019, 57 in 2020, 43 in 2021, 42 in 2022 and 41 in 2023. At these high latitudes, the temporal resolution was daily from 2017 onwards, and sub daily (about 66% of acquired days in a month) in 2016 when only Sentinel-2A was operating.
In order to characterise the ice surface at higher spatial resolution, Qaanaaq Glacier was also surveyed in 2023 by means of an UAV, model DJI Mini SE (carrying a FC7203 camera, with a 4.49 mm focal length and a pixel size of 1.76 × 1.76 μm). The UAV was manually flown over S1, S2, S3 and S4 (Fig. 1c) and an additional site in between S2 and S3 at a flying altitude of 30 m above the glacier surface, which led to high-resolution RGB mosaics (∼1 cm spatial resolution) of 27 000 m2 each. We used those data to better interpret the glacier surface. In particular, we analysed Sentinel-2 time series in the areas surveyed by the UAV to also evaluate the effect of surface heterogeneity (e.g. accumulation of impurities, presence of bédières).
3.1.3. Meteorological and topographic datasets
Meteorological and topographic data used in this study were acquired from ERA5-Land reanalysis and the ArcticDEM. ERA5-Land, produced by the European Centre for Medium-Range Weather Forecasts (ECMWF), is a reanalysis dataset with 9 km global spatial resolution and provides more realistic meteorological conditions, compared to a climate simulation, by means of a data assimilation scheme using observation data (Muñoz-Sabater and others, Reference Muñoz-Sabater2021). Surface air temperature, surface pressure, dew point temperature, and total precipitation derived from the ERA5-Land hourly product were used. We derived the terrain elevation with 32 m resolution from the ArcticDEM mosaic v4.1 (Porter and others, Reference Porter2023). ArcticDEM is the high-resolution elevation dataset or Digital Elevation Models (DEM) published by The Polar Geospatial Center (PGC); it is constructed from hundreds of thousands of individual DEMs extracted from various satellite imagery (Porter and others, Reference Porter2023). The reanalysis and DEM products were downscaled spatially using the method described in the Methods section. Two ArcticDEM strips (Porter and others, Reference Porter2022) over Qaanaaq Glacier were also utilised to derive the 2023 (2–4 June 2023) slope of the glacier at a high spatial resolution of 2 m.
To validate data from ERA5-Land reanalysis, we used meteorological conditions observed from the AWS at the SIGMA-B site (Fig. 1a). This AWS was established in 2012 and provides several semi-real-time data on an hourly frequency (calculated averaging 1 min data), including air temperature (accuracy ± 0.17°C), relative humidity (acc. ± 1%), wind speed (acc. ± 0.3 m s−1), wind direction (acc. ± 3°), up- and downward shortwave radiation (acc. ± 5%), up- and downward longwave radiation (acc. ± 10%), surface air pressure (acc. ± 0.30 hPa), snow height (acc. 1 cm), snow temperature (July 2012–present; acc. ± 0.15°C), and upward and downward near-infrared radiation (July 2022–present; acc. ± 5%). The surface air temperature, relative humidity (thermo-hygrometer HMP-155) and surface pressure (barometer PTB210) were derived from the dataset quality controlled by Nishimura and others (Reference Nishimura2023). Note that the dataset does not include precipitation data due to the absence of a rainfall sensor at the SIGMA-B site. Because the period of the dataset is between July 2012 and September 2020, we used the dataset from January 2013 to September 2020 for the evaluation.
3.2. Methods
3.2.1. Laboratory analyses for algae abundance and dust weight estimation
Algal cell concentration was calculated over slides by direct cell counting using an optical microscope (BX51; Olympus, Tokyo, Japan), after having filtered through a hydrophilised PTFE membrane filter (pore size 0.45 μm; Omnipore JHWP, Merck Millipore, Japan) (Tanaka and others, Reference Tanaka2016; Onuma and others, Reference Onuma, Takeuchi, Tanaka, Nagatsuka, Niwano and Aoki2018; Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023). From each ice sample, three slides were counted for the cell quantification and the final algal concentration (cells mL−1) was calculated from the cell mean of the three counts and the filtered sample volume. During the counting, different species were detected, with particular attention to peculiar species already described in previous papers (Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023), i.e. Ancylonema nordenskioldii (Fig. 1e), Ancylonema alaskana, Sanguina nivaloides, Phormidesmis priestleyi and Chroococcaceae cyanobacterium. On the other hand, mineral-dust abundance was quantified through the combustion method (Takeuchi and Li, Reference Takeuchi and Li2008; Onuma and others, Reference Onuma, Takeuchi, Tanaka, Nagatsuka, Niwano and Aoki2018). Samples were dried (60°C, 24 h) in pre-weighed crucibles and combusted at 500°C for 3 h in an electric furnace. This step was useful to remove all organic material (Di Mauro and others, Reference Di Mauro2024a). Finally, dust abundance was calculated per liter (g L−1), as a result of the combusted sample weight (only dust remained) and sample volume.
3.2.2. Reflectance ratio identification
With the aim of defining the most accurate reflectance ratio to quantify the algal abundance from Sentinel-2 images over the Qaanaaq Ice Cap, correlations (R2) and root-mean-square error (RMSE) based on a linear regression between field spectral ratios (model predictions) and algal concentration (observations) were calculated. This variable selection approach was applied in order to identify the reflectance ratio showing the highest correlation with algae concentration through a correlation matrix (Di Mauro and others, Reference Di Mauro, Garzonio and Baccolo2020a). A similar approach was also followed by Di Mauro and others (Reference Di Mauro, Fava, Ferrero, Garzonio, Baccolo, Delmonte and Colombo2015) to quantify mineral dust abundance. Thus, we calculated all possible spectral ratios in the 400–865 nm wavelength range using both ASD and RoX data. 216 225 linear regression models were then created between these ratios and the concentration of glacier algae. The same was conducted using dust concentration as a target variable, in order to define a ratio that is correlated with algae and uncorrelated with dust (Di Mauro, Reference Di Mauro2024a).
3.2.3. Glacier phenology analysis
Once the reflectance ratio was identified, the spatial distribution of algae abundance was estimated by applying the inverse formula calculated as the linear relationship between field spectrometer measurements and algae abundance from field samples.
Moreover, different parameters were estimated from satellite data, following a similar approach as proposed by Di Mauro and Fugazza (Reference Di Mauro and Fugazza2022) for Moderate Resolution Imaging Spectroradiometer (MODIS) data over the Alpine Region. As a starting point, we calculated the broadband albedo (hereafter referred as albedo) from Sentinel-2 imagery, applying the Liang broadband-conversion algorithm (Liang, Reference Liang2001) without anisotropic correction (directional effect), as suggested in different papers in cryospheric sciences (Naegeli and others, Reference Naegeli, Damm, Huss, Wulf, Schaepman and Hoelzle2017; Traversa and Fugazza, Reference Traversa and Fugazza2021; Traversa and others, Reference Traversa, Fugazza, Senese and Frezzotti2021; Hartl and others, Reference Hartl, Covi, Stocker-Waldhuber, Baldo, Fugazza, Di Mauro and Naegeli2025). Given the interest of the present study in analysing only ice surfaces, Sentinel-2 images were first masked using Normalized Difference Snow Index (NDSI) values higher than 0.30 to extract the Qaanaaq Ice Cap (this process was applied once over the most recent bare-ice image of 2023, i.e. 18 August 2023, and all images were masked based on it in order to consider the same portion of the ice cap in the analysis). This threshold, lower than the more common 0.40 employed for ice detection (Zhang and others, Reference Zhang, Wang, Shi and Yan2019), was adjusted to consider dark ice where high concentrations of debris or impurities are present (Salomonson and Appel, Reference Salomonson and Appel2004; Stillinger and others, Reference Stillinger, Rittger, Raleigh, Michell, Davis and Bair2023), such as the Qaanaaq Ice Cap. Consequently, with the aim of excluding snow surfaces and only focusing on bare ice, an albedo-based threshold was applied on each image, only considering values lower than 0.565, which was determined as a threshold for bare-ice albedo at the ice-ablation onset for the GrIS (Wehrlé and others, Reference Wehrlé, Box, Niwano, Anesio and Fausto2021), as it lies between the typical albedo values for the weathering crust and clean ice (Cuffey and Paterson, Reference Cuffey and Paterson2010). This threshold was also considered in algae abundance estimation, excluding pixels which represented snow or surfaces other than ice. Based on the albedo threshold, different pixel-based glacier phenology parameters (Di Mauro and Fugazza, Reference Di Mauro and Fugazza2022) were estimated, i.e.: length of bare-ice season (LOiS), start of bare-ice season (SOS), end of bare-ice season (EOS), minimum of summer albedo (min(α)), mean of summer albedo (mean(α)), maximum of summer algae abundance (max(AA)), mean of summer algae abundance (mean(AA)), day of the year (DOY) of summer algae abundance maximum (max(AA)-DOY) and the length of blooming season (LObS). Unlike in Di Mauro and Fugazza (Reference Di Mauro and Fugazza2022), albedo data were not filtered as they were found to be less noisy than MODIS. To further exclude possible noisy data points or transient snow events, we calculated the SOS as the first DOY when albedo was below 0.565 for three consecutive days; similarly, we defined the EOS as the first day when albedo exceeded 0.565 for at least three consecutive days. Algae abundance metrics were calculated only during the bare ice season as defined by these metrics. For each season among the 2018–2023 period, maps of phenological metrics were thus generated. Table 1 summarises the specifications of each ice phenological parameter, and Fig. 2 graphically represents how the variables were calculated for a certain year and location.
Figure 2. Plot of albedo (α) and glacier algae abundance estimated from Sentinel-2 in 2023 over the Qaanaaq Glacier in a small patch of 5 × 5 pixels located in the range between 400 m and 500 m a.s.l. Shaded areas are based on one standard deviation. The plot schematically represents how different glacier phenology variables were calculated.
Table 1. Specifications of glacier phenology variables
On the basis of the algae abundances thus calculated, we also estimated the equivalent carbon concentration in mg L−1. The conversion was made possible by assuming the cell biovolume of 1 µL equivalent to dry weight of 0.5 mg, resulting in 1 µL corresponding to 0.25 mg of carbon (C) (Fogg, Reference Fogg1967; Takeuchi and others, Reference Takeuchi, Dial, Kohshima, Segawa and Uetake2006). Here, we calculated the average of the equivalent carbon from the mean(AA) and the total equivalent carbon at the scale of the ice cap from the sum of the pixels retrieved by max(AA).
For a better understanding of the factors spatially and temporally controlling the algae abundance variability in the Qaanaaq area, we decided to move the focus from the ice-cap scale to a local scale, drawing attention to the Qaanaaq Glacier (Fig. 1b), where the two field campaigns from 2014 and 2023 were carried out. In this context, for retrieving variations of the different metrics (algae abundance and atmospheric variables) over the Qaanaaq Glacier, 6 square polygons of 2500 m2 (5 × 5 Sentinel-2 pixels) each were used. These areas of interest (AOI) were manually identified over the glacier extent at different elevations, about every 100 m, from 200 m to 800 m a.s.l., over homogeneous surfaces where UAV surveys were carried out (e.g. Fig. 2). Therefore, these AOI allowed us to evaluate the algae abundance variations along the Qaanaaq Glacier extent, over surveyed areas, with respect to temperature and snowfall variability.
3.2.4. Downscaling of atmospheric reanalysis data
To spatio-temporally assess the relationship between meteorological conditions and glacier algal blooms, we downscaled the atmospheric variables derived from ERA5-Land using the ArcticDEM. First, the horizontal resolution of the reanalysis and topographic data were interpolated to 60 m using a bilinear method before the downscaling processing. Their coordinate system was adjusted to the Sentinel-2 coordinate system during the interpolation to compare the downscaled meteorological conditions with Sentinel-2 images directly. In the paragraphs regarding the downscaling method, the interpolated reanalysis and topographic data are referred to as E5L and ADEM, respectively.
Surface air temperature, surface pressure, humidity, and precipitation were downscaled by the physically based downscaling approach. To obtain surface air temperature with 60 m resolution for each grid cell, the air temperature was downscaled using Eq. 1 proposed by Rouf and others (Reference Rouf, Mei, Maggioni, Houser and Noonan2020):
\begin{equation}\hat T = T + \Gamma \left( {\hat Z - Z} \right),\end{equation}where a variable with a hat (ˆ) denotes the downscaled data, and without it the original data. For example, T and Z are surface air temperature (K) and terrain elevation (m) derived from E5L, respectively. The air temperature and elevation differences between a target grid cell and its eight nearest neighbours at each time step are calculated, and a line is fitted to describe the T–Z relationship. The slope of the fitted line was used as temperature lapse rate Γ for the target grid cell. Z ˆ indicates terrain elevation derived from ADEM. The downscaling methods for surface pressure and specific humidity are shown in Eqs. 2 and 3 based on Rouf and others (Reference Rouf, Mei, Maggioni, Houser and Noonan2020), respectively:
\begin{equation}\widehat {Ps} = Ps{\ }exp\left[ { - \frac{{g\left( {\hat Z - Z} \right)}}{{R{T_m}}}} \right],\end{equation}
\begin{equation}\hat q = \frac{{0.622\hat E}}{{\widehat {Ps} - 0.378\hat E}},\end{equation}where Ps, g and R are surface pressure (hPa), the gravitational acceleration (9.81 ms−1) and the ideal gas constant (287 J kg−1 K−1), respectively. Tm is the mean air temperature between the T and T^. E^ in Eq. 3 is the downscaled water vapour pressure, which was obtained from the dew point temperature of E5L downscaled by the same method as T (Eq. 1), using the formula of Bolton (Reference Bolton1980). Total precipitation rate Pr (mm s−1) was downscaled using Eqs. 4 and 5 based on Thornton and others (Thornton and others, Reference Thornton, Running and White1997) and Liston and Elder (Reference Liston and Elder2006):
\begin{equation}\widehat {Pr} = Pr\left[ {\frac{{1 + f}}{{1 - f}}} \right],\end{equation}
\begin{equation}f = {X_s}\left( {\hat z - z} \right) + {X_i},\end{equation}where f means a factor to correct total precipitation with the elevation difference. Slope Xs and intercept Xi for each grid cell every time step were obtained from the Pr–Z relationship using the same method as Eq (1). Although the maximum absolute value of f was 0.95 in Thornton and others (Thornton and others, Reference Thornton, Running and White1997), the value is set to 0.8 in this study. They reported that values too close to 1.0 will result in excessive precipitation at strong elevation gradients. Because the elevation of ADEM is finer than the elevation data with 500 m resolution they used, the elevation gradient of ADEM is stronger. For this reason, we set the lower value in the range, so that their research does not degrade accuracy, to avoid excessive precipitation after downscaling. To obtain rainfall and snowfall amounts separately, the rain-to-snow ratio was calculated from the downscaled air temperature, pressure and humidity, and was applied to the downscaled total precipitation amount. The ratio was calculated using Eqs. 6 and 7 (Yamazaki, Reference Yamazaki2001):
\begin{align}s\left( {{T_w}} \right) &= \{ 1 - 0.5exp{\left( { - 2.2\left( {1.1 - {T_w}} \right)} \right)^{1.3}}{\ }if{ }\nonumber\\ {T_w} & \lt 1.1{ }0.5exp\left( { - 2.2{{\left( {{T_w} - 1.1} \right)}^{1.3}}} \right){\ }if{\ }{T_w} \geq 1.1 , \end{align}
\begin{equation}{T_w} = 0.584\left( {T - 273.15} \right) + 0.875E - 5.32,\end{equation}where s and Tw are rain-to-snow ratio and wet-bulb temperature (°C), respectively. The snowfall and rainfall amounts are given as sPr and (1 − s)Pr, respectively.
The validation results for the downscaled air temperature, surface pressure, and relative humidity with observational data at the SIGMA-B site in the Qaanaaq Ice Cap are shown in supplemental material (Fig. S1). The downscaled rainfall and snowfall amounts are also shown in the figure, although there are no observations at the site. The validation results indicate that the temporal changes in the downscaled air temperature, surface pressure and relative humidity agree well with those in the observed conditions at SIGMA-B.
4. Results and discussion
4.1. Glacier algae and dust abundance from field observations
The two datasets from 2014 and 2023 feature strong differences in terms of glacier algae abundance along the Qaanaaq Glacier, despite a comparable elevation range. Both datasets agree in showing a strong dominance of Ancylonema nordenskioldii (82% and 93% of total algae in 2014 and 2023, respectively), followed by Ancylonema alaskana (10% and 5%, respectively) and Sanguina nivaloides (5% and 2%, respectively). In general, samples from 2014 revealed an averaged abundance of 2.2 ± 2.4 × 104 cells mL−1, with a peak around the S2 site in August, i.e. 5.4 ± 4.2 × 104 cells mL−1, doubled compared to the 2.6 ± 1.9 × 104 cells mL−1 of July. All the sites present higher values in August than in July with the following concentrations: S1 with 0.4 ± 0.5 × 104 cells mL−1 in July and 2.3 ± 2.2 × 104 cells mL−1 in August, S3 with 1.3 ± 0.6 × 104 cells mL−1 in July and 2.6 ± 1.7 × 104 cells mL−1 in August and S4 with 0.9 ± 0.9 × 104 cells mL−1 in July and 1.0 ± 0.3 × 104 cells mL−1 in August. Thus, we observed an increase in algal presence between S1 (lowest peak) and S2, followed by a decrease in S3 and S4. The same pattern was observed in 2023, but with significantly higher values of algal abundance (one order of magnitude higher), leading to an overall average of 8.1 ± 5.5 × 105 cells mL−1 in August (please note that, in 2023, observations were conducted in August only). Again, the highest abundances were found in the S2 area, showing an average of 13 ± 0.5 × 105 cells mL−1. As in 2014, the lowest values were detected in the S1 area, with an average of 2.1 ± 2.0 × 105 cells mL−1. S3 and S4 surroundings showed averages of 9.5 ± 5.5 × 105 cells mL−1 and 7.3 ± 1.7 × 105 cells mL−1, respectively. Generally, 2014 samples provided an average abundance of organic matter of 0.13 g L−1 against 2.91 g L−1 of dust abundance (4% and 96% of the total particulate matter, respectively). The ratio between organic and inorganic materials tends to remain temporally and spatially stable, with the highest level of organic abundance in August (6%) and in S3 area (6%). The same pattern was also observed in 2023, with a generally lower abundance of organic matter, which, on average, presented 0.13 g L−1 (14%) against 0.75 g L−1 (86%) of dust. Here, the greatest organic content was recorded at the S2 site (18%) and the lowest at the S4 site (11%). These results are consistent with other findings in south-west Greenland (marginal area of the dark zone), where mineral dust was estimated as 94% of the total particulate matter (McCutcheon and others, Reference McCutcheon2021).
Such differences could be ascribed to the weathering crust formation which characterised the glacier in 2014. In fact, in the 2023 campaign it was observed that surface impurities were generally spread all over the glacier surface and, conversely, in 2014, due to the weathering crust formation, surface impurities tended to accumulate in cryoconite holes inside the crust for its high porous texture (Takeuchi and others, Reference Takeuchi, Sakaki, Uetake, Nagatsuka, Shimada, Niwano and Aoki2018; Woods and Hewitt, Reference Woods and Hewitt2023). This aspect could have affected the sampled specimens, which showed a lower concentration of impurities compared to 2023 in the surficial portion of the ice due to these weathering-crust characteristics.
4.2. Field spectroscopy data and reflectance ratio
As already observed in the previous section about the algae abundance in 2014 and 2023, differences were found also with regards to the field spectra (Fig. 3).
Figure 3. Field spectra acquired in 2014 (dotted line) and 2023 (solid line) in the areas of the four sites of Qaanaaq Glacier. Reddish columns in the plot represent the spectral amplitude of bands 4 and 5 (Red and Red-Edge1) of Sentinel-2.
In 2014, an increase in the reflectance was observed between July and August measurements, in particular in locations S2 and S3 where an average increase of +0.09 and +0.08 was encountered. Focusing on August measurements, as already observed in algae abundance from the samples, S2 was the darkest site, with an average albedo of 0.44 ± 0.09 and lowest value recorded in 2014, i.e. 0.36. The other sites were as follows: S1 averaged 0.66 ± 0.07, S3 averaged 0.52 ± 0.05 and S4 averaged 0.54 ± 0.03. On the other hand, 2023 spectroscopy data revealed much lower albedo values in August, which is in accordance with a much higher algae abundance as observed in the cell count. In fact, in the S2 area, which remained the darkest analysed zone, the average of measured spectra showed an albedo of 0.19 ± 0.04, therefore 0.25 lower on average compared to 2014. All the other sites showed generally lower values in 2023, as follows: S1 average of 0.41 ± 0.09, S3 average of 0.39 ± 0.16 and S4 average of 0.20 ± 0.06. Overall, the lowest albedo was observed in the S2 area, with 0.13. Analysing the averaged spectra of the 2014 and 2023 campaigns, we observed the typical decrease of reflectance for wavelengths shorter than 750 nm, due to algae presence on bare ice (Dauchet and others, Reference Dauchet, Blanco, Cornet and Fournier2015; Cook and others, Reference Cook, Hodson, Taggart, Mernild and Tranter2017). Moreover, an absorption feature located at 680 nm was detected as well in most spectra. This feature has been usually linked to the presence of Chlorophyll-a in glacier ice (Takeuchi, Reference Takeuchi2002; Remias and others, Reference Remias, Holzinger, Aigner and Lütz2012; Di Mauro and others, Reference Di Mauro, Garzonio and Baccolo2020a). These spectral behaviours were observed in 2023, especially in S2 and S3 areas and secondarily in S4, as shown in Fig. 3 (blue and green solid lines). The higher reflectances observed in 2014 than in 2023 field campaigns could be ascribed to two main aspects: first, a clean weathering crust was observed to cover large areas of the Qaanaaq Glacier in 2014, brightening the glacier surface (Traversa and Di Mauro, Reference Traversa and Di Mauro2024). This can explain why certain specimens in the S1 area, despite a lower algae abundance in 2023 (Fig. 3), were less reflective than in 2014. Secondly, as already stated before, during the 2023 campaign, the surface impurities were observed to be generally spread all over the glacier surface and in 2014 surface impurities tended to accumulate in cryoconite holes inside the weathering crust (Takeuchi and others, Reference Takeuchi, Sakaki, Uetake, Nagatsuka, Shimada, Niwano and Aoki2018; Woods and Hewitt, Reference Woods and Hewitt2023).
With the aim of defining the highest correlations among spectral ratios and algae abundance based on field observations, the correlation matrix was calculated. The correlation matrix, for both campaigns, resulted in the highest correlations in the red and far-red portion of the spectrum (Fig. 4a).
Figure 4. Matrices of reflectance ratios coloured on the basis of (a) R2 and (b) RMSE values with algae abundances from the 2014 and 2023 campaigns. The red rectangle represents the Sentinel-2 ratio of bands 5 and 4. Scatter plots of (c) algae and (d) dust abundances estimated from field samples, correlated with corresponding reflectance ratio from field-spectroscopy measurements (averaged on the Sentinel-2 bands).
Here, the highest correlation (R2) among spectral ratios and algae abundance was detected between 695/681 or 695/687 nm wavelengths, providing an R2 of 0.90 and among the lowest RMSE values (Fig. 4b), of 163 874 and 161 587 cells mL−1, respectively. Therefore, the photosynthetic absorption of glacier algae in the red wavelengths can be exploited for their estimation from satellites. In the Sentinel-2 specific case, this is possible by taking advantage of bands 5 and 4 (centre wavelengths of 705 ± 15 and 665 ± 30 nm, respectively), which are the closest bands to the identified ratio from field measurements. The identified Sentinel-2 ratio was then calculated over the field spectroscopy measurements of the 53 samples from 2014 and 2023. The Sentinel-2 ratio presented a high correlation with algae abundances when estimated from field measurements, showing a R2 of 0.84 (Fig. 4c). The identified equation was then inverted in order to estimate glacier algae abundance (in cells mL−1) from Sentinel-2, as follows:
\begin{equation}AA\, = \,\left( {rati{o_{705/665}}\,*\,8,034,887.96} \right)\, - 7,564,737.66\end{equation}with 95% confidence intervals for the intercept (–8.53 × 106 to −6.60 × 106) and the slope (7.06 × 106 to 9.01 × 106). Values of the ratio higher than 0.941 will lead to a positive concentration of glacier algae. Lower values shall be regarded as clean ice or dirty ice with a predominance of inorganic impurities.
Moreover, this reflectance ratio was tested against dust abundances and a non-significant (p-value > 0.05) correlation (R2) was observed, i.e. 0.03 (Fig. 4d). These results support the application of the ratio of bands 5 and 4 for estimating algae abundance, avoiding dependencies from dust presence.
It is important to note that this equation can be applied to estimate algae abundance over ice surfaces when the supraglacial biological community is dominated by Ancylonema nordenskioldii species, as in the Qaanaaq area. In fact, this reflectance ratio slightly differs from previous attempts to estimate glacier algae abundance in other regions of the GrIS or other glacierised areas. In particular, despite being similar to the ratio proposed by Di Mauro and others (Reference Di Mauro, Garzonio and Baccolo2020a) for the European Alps, this new ratio takes advantage of the band 5 of Sentinel-2 instead of band 6, which has a higher wavelength of 740 ± 15 nm. Nevertheless, the present ratio is in accordance with the ratios proposed by (Wang and others, Reference Wang, Tedesco, Xu and Alexander2018; Wang and others, Reference Wang, Tedesco, Alexander, Xu and Fettweis2020) for Sentinel-3 and MERIS (709/674 nm and 709/665 nm, respectively), both based on field measurements collected in western GrIS. Especially the MERIS ratio was based on the same portions of the spectrum proposed in this paper, supporting its application at the Greenland scale.
4.3. Interannual variability over the Qaanaaq Ice Cap by means of glacier phenology metrics
The glacier phenology approach allowed us to estimate different metrics over the entire Qaanaaq Ice Cap per each summer season, from 2016 to 2023 (Fig. 5-7; Fig. S2-S7). In general, for all the analysed metrics, we observed specific patterns in the different years, due to the varying duration of snow cover over the ice. In fact, we observed that some years presented only a small portion of the ice cap with bare ice at the surface, due to long-lasting snow cover. In contrast, in some other years most of the ice cap was snow-free for a long period. In particular, summer seasons 2017, 2018 and 2021 presented more than 40% of ice-cap area with winter snow cover lasting from June to September and above all in summer 2018 75% of the area was characterised by long-lasting summer snow cover, while only marginal glaciers showed bare ice at the surface. On the other hand, summer seasons 2019, 2020 and 2023 featured less than 20% of snow-covered area, and the minimum was reached in 2023 where only portions of the ice cap over 1000 m a.s.l. were characterised by the presence of snow (6% of the whole ice cap).
Figure 5. Phenology maps retrieved from Sentinel-2 images representing the summer (June–September) minimum albedo from 2016 to 2023 over the Qaanaaq Ice Cap. All the images were masked based on the August 2023 ice-cap extent. ESRI Light Gray map in the background.
This behaviour reflects in the summer albedo of the ice cap, which shows the lowest mean(α) in 2019, 2020, 2022 and 2023 summer seasons. Particularly, albedo in 2019 shows the lowest mean(α) from 0 to 800 m a.s.l. (lowest values between 500 and 700 m a.s.l., <0.30), while 2022–2023 shows the lowest mean(α) from 800 m a.s.l. to the ice-cap top. Similar results were obtained focusing on the min(α) observed during the summer period (Fig. 5).
The lowest values were again observed in 2019 in the first 800 m of elevation, where the darkest area was observed between 400 and 700 m a.s.l., showing on average min(α) values <0.20. The overall lowest min(α) was observed in 2020 though, despite an average higher albedo than 2019 below 800 m a.s.l. On the other hand, the 2018 summer season was the brightest among the analysed years, with mean(α) and min(α) always >0.40, followed by 2017 and 2021.
The results of the analysis of mean(α) and min(α) are reflected in the LObS (Fig. 6).
Figure 6. Phenology maps retrieved from Sentinel-2 images representing the length of the blooming season (days) from 2016 to 2023 over the Qaanaaq Ice Cap, where 0 means no days of algal bloom in the summer and 75 means blooming lasting until mid-August. All the images were masked on the August 2023 ice-cap extent. ESRI Light Gray map in the background.
In fact, the 2019-summer season shows LObS values for the 0–800 m elevation range longer than a month, resulting in algal blooms lasting for most of the bare-ice conditions. In this case, the 2019 season presents the highest LObS values across the entire elevation range of the ice-cap, followed by 2020, 2022 and 2023. In fact, in these seasons the snow cover was immediately removed from the surface at lower elevations, leaving exposed ice on the surface since the beginning of June (lowest observed SOS) allowing the algae to start blooming. However, the blooming season persisted until the end of the summer in 2019 (about two months below 500 m a.s.l.), but not in 2022 and 2023, when possibly snowfalls occurred in August and September. In contrast, summer 2017 presents the lowest LObS values: areas at elevations higher than 200 m a.s.l. showed a length shorter than one month across the entire summer season, focused between the end of July (SOS) and the beginning of August (EOS). The other short bare-ice summers were 2021 and 2016, which always had about or less than one month of LObS over the ice cap. Finally, 2018 summer presented low LObS values as well (lower than one month), with the exception of the lowest elevation range (up to 200 m a.s.l.), showing more than a month of blooming.
What is observed in the spatial variability of albedo and LObS is partially reflected by the results of algae concentration. In fact, generally mean(AA) (Fig. 7) showed the highest values in the period from 2019 to 2023 (with the exception of summer 2021), when the lowest albedo and longest LObS were estimated.
Figure 7. Phenology maps retrieved from Sentinel-2 images representing the summer (June–September) mean of algae abundances (cells mL−1) from 2016 to 2023 over the Qaanaaq Ice Cap. All the images were masked on the August 2023 ice-cap extent. ESRI Light Gray map in the background.
The highest mean algae abundance was observed in summer 2020. During that summer, between 100 and 800 m a.s.l., mean(AA) showed values higher than 10.0 × 105 cells mL−1, even if high variations among spatially close pixels were observed. In contrast, summer 2019 showed a high mean(AA) at different elevations too, but lower spatial heterogeneity, with concentrations >9.0 × 105 cells mL−1 between 100 and 700 m a.s.l., exceeding the average mean(AA) of 10.0 × 105 cells mL−1 between 400 and 500 m a.s.l. When converted to equivalent carbon per L, the results are 270 and 300 mg C L−1, respectively. Similar results were also obtained for the 2022 and 2023 summer seasons. 2017 and 2018 had the lowest mean(AA), always lower than 8.0 × 105 cells mL−1 (240 mg C L−1) at all the elevations. In this context, the highest values of algal concentrations (max(AA), mostly higher than 10.0 × 105 cells mL−1) were observed between the last week of July and the first two weeks of August (max(AA)-DOY) in all the analysed years. In conclusion, the algal abundance results reflect the albedo behaviour previously described, with the exception of the 2020 summer season, which showed the highest algal abundances and thus the highest concentration of total equivalent carbon at the scale of the ice cap (2114 × 103 g C), but not the lowest min(α), which was instead recorded in 2019. Table 2 summarises these results.
Table 2. Spatial averages of mean(AA), equivalent carbon and min(α) and total equivalent carbon at the scale of the Qaanaaq Ice Cap
Another relevant output of the analyses is the geographical distribution of the algae concentration over the ice cap. In fact, in addition to the elevation differences already observed, i.e. higher abundances at lower elevations, also geographic differences were encountered. In general, the Qaanaaq Glacier can be divided into 11 ice sheds (RGI 7.0 Consortium, 2023), where a major ice divider splits the ice cap from south-east to north-west. Across the different years, higher abundances were always estimated on the ice divides of the east side and especially at the margins of the ice cap. In fact, on a spatial average, east ice divides showed at least 0.5 × 105 cells mL−1 higher than west side, with a maximum of +1.7 × 105 cells mL−1 in 2020. Outlet glaciers generally present higher algae abundances, even if not at their margins, but rather in their middle part, as already shown before for the lower elevations, e.g. at 400–500 m a.s.l. This pattern was observed in most of the Qaanaaq Ice Cap outlet glaciers, even if more pronounced over the west side. There, as shown in the next subsection, glacier terminus present brighter surfaces and general lower algae abundances in respect to middle elevation areas of outlet glaciers.
4.4. Role of meteorological and glaciological parameters on the algae abundance variations at Qaanaaq Glacier
To investigate the role of different parameters in the variability of algae abundance on Qaanaaq Glacier, we focused on specific locations along the glacier extent, from 200 m to almost 800 m a.s.l., where different UAV surveys were carried out (Fig. 8b).
Figure 8. (a) UAV view (about 60 m above the surface) of the Qaanaaq Glacier taken on 9 August 2023 as seen from site S4. Images of the Qaanaaq Glacier represented as: (b) RGB Sentinel-2 acquisition (18 August 2023), where trend analyses of algae and meteorological parameters were carried out (black squares); (c) algae abundance and (d) albedo retrieved from the Sentinel-2 18 August 2023 acquisition and (e) slope derived from strips of the ArcticDEM (2 June 2023).
Temporally, similar variations as observed at the ice-cap scale were found. In fact, 2017, 2018 and 2021 showed low variability and abundance of algae. Conversely, summer 2019 and 2023 were characterised by an exceptional algal abundance, higher than 10 × 105 cells mL−1 at the peak of the blooming season. Figure 9 displays two examples of one low algae abundance season (2018) and one of high algae abundance season (2019).
Figure 9. Temporal variation (summer 2018 on the left and 2019 on the right) of algae abundance (green lines) retrieved from Sentinel-2 acquisitions, air temperature (red lines) and snowfall (light-blue columns) from downscaled atmospheric reanalysis data. Shaded areas represent one standard deviation. Each subplot refers to a specific elevation range of the Qaanaaq Glacier, whose locations are represented in Figure 8b (black squares).
Especially during the two summer seasons 2019 and 2023, as well as in other years, the pattern showed a start of the blooming season between the end of June and the middle of July, in accordance with the observations from the Qaanaaq Ice Cap, and the peak which is mostly reached close to the middle of August. Generally, the increase in algae abundance between July and August tends to be gradual, with a sudden decrease after having reached its peak. In conjunction with this sudden decrease, meteorological observations usually present significant snowfall events (Fig. 9), which occur between the end of August and the beginning of September, possibly halting the algae blooming at the surface. From the meteorological data, we also observed a correspondence between consecutive days when surface temperature remained positive (daily averages) and no snowfall events were observed. In particular, we observed that the bloom of algae is related to a few consecutive days (e.g. five days) of daily positive temperature, and then strongly increases in abundance in relation to the number of days of positive temperature and thus length of melting season (Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023; Roussel and others, Reference Roussel2024). In this context, the two years (2019 and 2023) when highest algae abundances were observed are related to the two longest periods of consecutive positive temperature in the analysed years, ranging between 57 and 64 days. Other years presented many consecutive days of positive temperature, such as 2016 (42 days), 2020 (44 days) and 2022 (45 days). However, during these seasons, early snowfalls occurred during the blooming, without being followed by many days of positive temperature. In 2023, even if a significant amount of snow (26 mm) fell in five days around 20 August, the temperature remained positive for another 24 days, allowing the snow to melt and the algae to continue blooming. Conversely, in some other years, such as 2017, 2018 and 2021, few days of consecutive positive temperature were recorded (around one month) and many snowfall events were recorded all across the summer. Possibly, the mix of conditions respectively contributed to or opposed algae blooming over the Qaanaaq Glacier.
In this context, linear regressions (R2) were calculated between algae abundance and air temperature and algae abundance and surface albedo, using monthly (June–September) averages of these variables in the period 2016–2023. In particular, in middle range (elevations) areas, good R2 were found, especially between 400 and 600 m a.s.l. of the Qaanaaq Glacier, with a maximum R2 of 0.70 (500–600 m a.s.l.) and 0.40 (400–500 m a.s.l.) for albedo and temperature, respectively (Table 3).
Table 3. Linear regression coefficients (R2) between algae abundances (AA) and albedo (α) and air temperature (T)
The regressions are retrieved from monthly means from 2016 to 2023 for summer months (June–September) over the eight locations represented in Fig. 8b. Albedo averages were calculated by taking into account only those values <0.565 (Wehrlé and others, Reference Wehrlé, Box, Niwano, Anesio and Fausto2021).
* represents significant (p-value) regressions at the 99% confidence level;
** represents significant (p-value) regressions with 95% confidence and
*** represents not-significant regressions (p-value > 0.05)
However, the R2 of these regressions become much lower moving towards the terminus of the glacier and at elevations >700 m a.s.l., where non-significant regressions were calculated. As for snowfalls, regressions were not significant in all the analysed areas, suggesting the non-linearity of their relation with algae abundances.
Moreover, significant spatial differences were encountered along the glacier extent. In fact, the highest concentrations of algae were found in the middle portion of the glacier, between 400 m and 700 m a.s.l. (Fig. 8c-9). Conversely, the lowest concentrations were estimated especially over the terminal portion of the glacier, and also on the highest (elevation) parts of it. In these latter areas, during the brightest seasons (2017–2018), almost no algae were estimated by the satellite (Fig. 9). These findings are reflected in the observed albedo by Sentinel-2, which showed high values over these regions (Fig. 8d). Particularly, bright portions of the glaciers were observed at its margins especially at the terminus or at the sides of the outlet glacier where the ice-cap margins are located. This pattern was also observed on other neighbouring outlet glaciers and is opposite to what was observed over Alpine glaciers (e.g. Morteratsch Glacier), where the highest concentrations are displayed at the margins (Rossini and others, Reference Rossini2018; Di Mauro and others, Reference Di Mauro, Garzonio and Baccolo2020a; Millar and others, Reference Millar, Broadwell, Lewis, Bowles, Tedstone and Williamson2024). Possibly, at Qaanaaq Ice Cap, these areas are affected by a relatively higher surface slope (>10°, Fig. 8e), which favours the flow of surficial water, washing away the algae from the glacier ice. This behaviour would explain the reason behind the observed low abundance of algae, above all at the terminus of Qaanaaq Glacier.
These results are in line with previous observations in the area, where findings demonstrated that the abundance of glacier algae increases with the length of the melting season (Onuma and others, Reference Onuma, Takeuchi, Uetake, Niwano, Tanaka, Nagatsuka and Aoki2023), as well as in other areas of GrIS (Feng and others, Reference Feng, Cook, Naegeli, Anesio, Benning and Tranter2024). There, accordingly, bloom events were found to start from the end of June until August, in relation with increasing temperature and solar radiation (Shimada and others, Reference Shimada, Takeuchi and Aoki2016) and availability of liquid water at the surface. Water retention due to glacier surface roughness (depressions) was discovered to host useful nutrients (phosphorus) which set up an ideal habitat for algae development (McCutcheon and others, Reference McCutcheon2021). The presence of dark depressions was further confirmed by UAV observations, which identified rough surfaces, especially in the middle area of the Qaanaaq Glacier, where the highest concentrations of algae were observed. Moreover, the washing effect of water was previously found to reduce algae abundance (Williamson and others, Reference Williamson2020), strengthening the hypothesis behind the lower concentrations found at the margins of the outlet glaciers, where the presence of steeper surfaces can make it easier for algae to be washed away.
4.5. Limitations on glacier algae abundance estimation
Despite the relevant results obtained in the present research, thanks to the proposed remote-sensing methodology based on field measurements, cautions and considerations are further needed. First, it is important to highlight that the provided methodology to estimate algae abundance was developed only for ice surfaces and not for snow. Its application over snow surfaces could provide unreliable results. Moreover, despite being developed for glacier algae, the Qaanaaq supraglacial biological community is characterised by a strong dominance of only one species, i.e. Ancylonema nordenskiöldii. Thus, the application of this approach needs to be taken with caution when applied over glaciers with a more diversified algal community or where cryoconite granules have a stronger role in darkening the ice surface. Additionally, the methodology was based on two datasets taken by two different operators nine years apart and, despite efforts carried out in 2023 to repeat the sampling procedure, small differences could have occurred. For example, even small differences in ice sampling depth could have affected the estimations of algae and dust concentrations.
However, the results obtained in this research are in line with previous attempts in other regions of the Earth (Wang and others, Reference Wang, Tedesco, Xu and Alexander2018; Wang and others, Reference Wang, Tedesco, Alexander, Xu and Fettweis2020; Di Mauro and others, Reference Di Mauro, Garzonio and Baccolo2020a), thus suggesting the broader applicability of the methodology. In addition, we demonstrated the low dependence of the model on the presence of dust. Another point which deserves attention is the possible application of the inverse equation to estimate algae abundance from the satellite band ratio. In fact, field measurements at Qaanaaq Glacier, especially in 2023, were characterised by a high concentration of glacier algae, making the estimation less reliable when applied over lower algae concentrations and values of the band-ratio. For ratios lower than ∼0.94, the application of the inverse formula will provide negative algae abundances; also, note that the model is based on field-estimated ratios between 0.94 and 1.18, i.e. the lower and upper bounds found in this study.
Considering the application of the method to Sentinel-2 and the derived estimated algae abundance and its variability, the main limitations could be found in the temporal and spectral characteristics of this dataset. In fact, despite the high temporal resolution of Sentinel-2 at high latitudes, the high cloud cover persisting over the study area significantly reduced the temporal resolution of the analyses, providing at best 47% of investigated days between June and September in the 2020 summer season. Moreover, it is relevant to note that this high temporal resolution of Sentinel-2 is available only at high latitudes. Therefore, in e.g. ice-sheet-wide studies, where acquisitions are needed for much lower latitudes (<60° N), the application of Sentinel-2 could be limited and a trade-off with its high spatial resolution should be considered. In this context, the integration of Sentinel-3 products, which share similar spectral characteristics with Sentinel-2, but a daily temporal resolution even at lower latitudes, could benefit the research, despite the lower spatial resolution (300 m). Additionally, the accuracy of derived reflectance from the Sentinel-2 L2A product is relevant to the overall uncertainty in the estimation of algae abundance. In past studies, the accuracy of Sentinel-2 L2A derived reflectance was analysed over both dark (land and coastal areas) and bright surfaces (snow and ice), and a slight positive offset (<10%) was generally found compared to field measurements (Gorroño and others, Reference Gorroño, Guanter, Graf and Gascon2024; Naethe and others, Reference Naethe2024; Di Mauro and others, Reference Di Mauro2024b), as also observed in the present research (Fig. 3). In view of this, the application of a corresponding ratio to hyperspectral satellite data (e.g. PRISMA and EnMAP satellite missions) could improve the accuracy of the satellite-based estimations, by taking advantage of narrower spectral bands as in field spectroscopy, which showed the highest correlation with field-based algae abundance (i.e. 695/681 nm or 695/687 nm). Finally, the method suggested here could be applied to other multispectral satellites to further enhance data availability, with the aim of improving the spatial and temporal analyses (e.g. GCOM-C and MODIS products, Sentinel-3, Landsat and PlanetScope).
5. Conclusions
The study focuses on glacier algae and their glaciological role over the Qaanaaq Ice Cap, in north-western Greenland, by means of field and satellite observations. Glacier surface samples were collected during two field campaigns over the Qaanaaq Glacier, one of the outlet glaciers of the Qaanaaq Ice Cap, in 2014 and 2023 and then processed in the laboratory to estimate the algae and dust abundances. Contextually, at each sampling site, reflectance measurements were collected in the VIS-NIR portion of the spectrum by using different spectrometers. Differences were found among the two campaign samples: in fact, 2023 measurements provided much higher algae concentrations (by an order of magnitude, with overall averages of ∼2 × 104 and ∼8 × 105 cells mL−1) and much lower reflectances (−20%), when the effect of chlorophyll-a was evident with absorption at around 680 nm. In general, a strong dominance of Ancylonema nordenskiöldii species was found in both campaigns (>80%). Based on field observations, we defined the best reflectance ratio capable of estimating algae abundance, i.e. 695/681 nm or 695/687 nm, close to previous attempts in other areas (Wang and others, Reference Wang, Tedesco, Xu and Alexander2018; Wang and others, Reference Wang, Tedesco, Alexander, Xu and Fettweis2020; Di Mauro and others, Reference Di Mauro, Garzonio and Baccolo2020a). Subsequently, the ratio was applied to Sentinel-2 bands (band 5/band 4), allowing us to estimate the glacier algae abundance at the ice cap scale and to analyse their variation over time from 2016 to 2023 (summer seasons). From this analysis, thanks to a phenology approach (Di Mauro and Fugazza, Reference Di Mauro and Fugazza2022), we found strong heterogeneity in algae variation among the different seasons, where especially summers 2019, 2020 and 2023 revealed exceptional algal growth and abundance (>1 × 106 cells mL−1). Particularly, it is noteworthy that in 2020 more than two tons of equivalent carbon (2114 kg C) were estimated at the scale of the ice cap (∼270 km2). At Qaanaaq Glacier scale, we also estimated algae variability in comparison with meteorological (temperature, snowfall) and topographic (slope) parameters, observing that algae tend to grow and reach high abundances (in mid-August) when several consecutive days of positive temperature (ranging from 57 to 64 days) and corresponding snowfall hiatus take place, especially over relatively flat areas (algae presented lower concentrations on glacier margins, where steeper surfaces are found).
Further research could focus on different aspects revealed in the present research. First, the application of a similar approach based on other satellites, to test hyperspectral resolutions (PRISMA, EnMAP and CHIME) as a possible improvement in remote glacier algae estimation, or to widen the study area to the entire GrIS, thanks to the application of the identified reflectance ratio to other multispectral satellites with daily temporal resolution (e.g. GCOM-C, Sentinel-3 and MODIS products). Moreover, the application of other satellite products having a longer temporal acquisition history (e.g. MODIS) would allow estimating statistical correlations among the phenology metrics. In fact, in the present research, these correlations were not calculated in view of the few available years (eight), which would have strongly affected the statistics. Finally, given the not-linearity of relations among algae and meteorological and glaciological parameters found in this research (with few exceptions at certain elevations), future research should focus on non-linear modelling for simulating algae abundance and their growth over time.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10096.
Acknowledgements
This study was supported by the project Light-Absorbing ParticleS in the cryosphere and impact on water resourcEs (LAPSE) (No. PRIN 202283CF7) and the Ministry of education, Universities and Research (MUR) of Italy, which funded the fellowship of Dr. G. Traversa. The study was also supported by Horizon2020 (Grant Agreement No. 871120), INTERACT Transnational Access, which funded the Triple UP-scaling of Ice-Light-Absorbing particles at Qaanaaq ice cap (TUPILAQ) project, leading to the 2023 Qaanaaq field campaign, and by the ArCS II International Early Career Researchers Program funded by the National Institute of Polar Science of Japan, which supported a visiting period of Dr. G. Traversa at the Earth Observation Research Center (JAXA) and Chiba University, in the context of the Biological Darkening over the Qaanaaq ice cap (BDQ) project and by the JSPS-Kakenhi Grant-in-Aid (23221004, 23K17036, 24H00260). The authors also acknowledge Dr. T. Suzuki (JAXA) for supporting the laboratory analyses and two anonymous reviewers for their comments and suggestions which greatly improved the manuscript.
Author contributions
G.T., Y.O. and B.D.M. conceived the idea of this work. G.T. and Y.O. wrote most of the paper and collected field data, with the support of F.C.Q., and performed laboratory analyses. N.T. supported laboratory analyses for algae and dust concentrations. R.G. run the correlation matrix correlations. G.T., Y.O. and D.F. performed calculations based on remote sensing and reanalysis datasets. B.D.M. and N.T. supervised the work and supported and revised the writing of the manuscript.
Competing interests
The authors declare none.
Open access
