# Quantified positive radiative forcing at a greening Canadian boreal-Arctic transition over the last four decades

Florent Dominé, Arthur Bayle, Maria Belke-Brea, Esther Lévesque, Ghislain Picard

<figure><img src="image_1.png" /></figure>

## To cite this version:

Florent Dominé, Arthur Bayle, Maria Belke-Brea, Esther Lévesque, Ghislain Picard. Quantified positive radiative forcing at a greening Canadian boreal-Arctic transition over the last four decades. Remote Sensing of Environment, 2025, 322, pp.114715. 10.1016/j.rse.2025.114715. hal-05344788

HAL Id: hal-05344788

https://hal.science/hal-05344788v1

Submitted on 3 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Quantified positive radiative forcing at a greening Canadian boreal-Arctic transition over the last four decades

Florent Dominé<sup>a,b,c,1,*</sup>, Arthur Bayle<sup>d,e,*†</sup>, Maria Belke-Brea<sup>a,b</sup>,
Esther Lévesque<sup>b,e</sup>, Ghislain Picard<sup>f</sup>

## ARTICLE INFO

Edited by Jing M. Chen

Keywords:

Forest-tundra ecotone

Arctic greening

Lichen

Dwarf birch

Black spruce

Spectral albedo

Shortwave forcing

## ABSTRACT

Climate warming in northern and Arctic regions drives vegetation growth and shifts species distribution. In northern Quebec's Boreal-Arctic transition (forest-tundra ecotone), this is seen in the replacement of lichen by shrubs, primarily dwarf birch. These changes impact surface albedo, contributing to climate forcings with broad consequences. This study measures vegetation changes in Tasiapik valley near Umijuaq, Quebec, using a combination of (1) hyperspectral data (347–2400 nm) collected from 62 vegetation assemblages, including lichen, dwarf birch, willow, and spruce, to calculate broadband albedo, and (2) remote sensing data from Landsat satellites over 1984–2023. By combining these data, the proportion of vegetation type for each pixel was determined at the beginning and end of the 40-year period. The areal coverage of six main vegetation types was quantified over the 9.25 km² valley. The most significant change was lichen replacement by dwarf birch with lichen understory, leading to an albedo reduction from 0.233 to 0.168 and a summer shortwave forcing of 11.17 W m⁻². At the valley scale, the spatially-averaged summer forcing was 2.16 W m⁻² when considering all observed vegetation changes. These values, lower than those in previous Norwegian studies, highlight the spatial variability of shortwave forcing due to lichen replacement. We observed that the vegetation change producing the greatest positive radiative forcing also caused the strongest greening. This suggests that Landsat-based greening may be used as a proxy for surface albedo change on an Arctic scale. This unique combination of ground and satellite data allows quantification of a direct, first-order effect of Arctic shrubification.

## 1. Introduction

The forest-tundra ecotone (FTE) (Payette et al., 2001) is a critical circumpolar transition region between the Arctic tundra and the boreal forest. It extends over 13,400 km at the southern edge of the Arctic (Walther et al., 2019). Its latitudinal extent can reach several hundred km, which is the case in northern Quebec (Gamache and Payette, 2005). Climate change is driving the FTE northward, as manifested by some expansion of trees (Beck et al., 2015; Gamache and Payette, 2005), and also by shrub expansion, as observed in northern Quebec (Ropars and Boudreau, 2012; Tremblay et al., 2012), western Canada (Lantz et al.,

2013), Alaska (Tape et al., 2006) and Siberia (Frost and Epstein, 2014). As vegetation transitions from tundra to shrubland to forest, the broadband albedo (i.e., the albedo integrated over the solar spectrum) decreases (Beringer et al., 2005) and the resulting shortwave (SW) forcing positively feeds back on climate. Consequences include not only atmospheric warming on land surfaces but also sea ice decline because of the greenhouse effect due to extra atmospheric moisture generated by enhanced evapotranspiration (Bonfils et al., 2012; Chen et al., 2022; Lafleur and Humphreys, 2018).

While monitoring albedo changes at large scale is a challenge, changes in vegetation are routinely monitored by remote sensing, often

<footnote><sup>a</sup> Takuvik Joint International Laboratory, Université Laval (Canada) and CNRS-INSU (France), Québec City, Canada</footnote>

<footnote><sup>b</sup> Centre d'Études Nordiques, Université Laval, Québec City, Canada</footnote>

<footnote><sup>c</sup> Department of Chemistry, Université Laval, Québec City, Canada</footnote>

<footnote><sup>d</sup> Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, LECA, Grenoble, France</footnote>

<footnote><sup>e</sup> Département des sciences de l'environnement, Université du Québec à Trois-Rivières, Trois-Rivières, Québec, Canada</footnote>

<footnote><sup>f</sup> Université Grenoble Alpes, CNRS, IRD, Grenoble INP, IGE, Grenoble, France</footnote>

<footnote>* Corresponding authors.
E-mail addresses: florent.domine@gmail.com (F. Dominé), arthur.bayle.env@gmail.com (A. Bayle).</footnote>

<footnote n="1">¹ These authors should be considered joint first authors</footnote>
using the normalized difference vegetation index, NDVI, (Boelman et al., 2011; Juszak et al., 2014; Rees et al., 2020) which is considered a proxy for the amount of photosynthetically active vegetation. For example, data from Landsat satellites, starting in 1984, yield the NDVI evolution of Arctic regions, from which trends in greening (i.e., the increase in photosynthetically active vegetation) are deduced (Bayle et al., 2022; Berner and Goetz, 2022; Berner et al., 2020; Crichton et al., 2022; Delbart et al., 2005; Edwards and Treitz, 2017; Ju and Masek, 2016; Pattison et al., 2015). More specifically, the evolution of the FTE near Umiujaq in northern Quebec has been investigated using Landsat-derived NDVI (Beck et al., 2015; Lemay et al., 2018). In this region, prostrate vegetation is often comprised of lichen, a vegetation type with a high albedo (Aartsma et al., 2020; Finne et al., 2023; Payette and Delwaide, 2018; Reinhardt et al., 2022). Its replacement with erect vegetation such as dwarf birch (Tremblay et al., 2012) can be expected to lead to large changes in broadband albedo and therefore to a significant modification of the surface energy budget.

Quantifying the impacts of these vegetation changes on the SW energy budget over large areas and extended timescales presents several challenges. First, it requires detecting spatially comprehensive vegetation change at scales suited to the vegetation assemblages of the FTE, typically ranging from decameter to meter scales (Assmann et al., 2020; Dobbert et al., 2021; Putkiranta et al., 2024; Yang et al., 2023). Second, it necessitates the knowledge of the broadband albedo of each species or more practically of each vegetation assemblage affected by these changes. Broadband albedo is not strictly speaking an intrinsic property of a surface because it depends on variables such as the solar zenith angle (SZA) (Yang et al., 2020), cloudiness (Wolf et al., 2024), terrain slope (Aartsma et al., 2020) and surface roughness (Cierniewski et al., 2017). A rigorous calculation of broadband albedo therefore requires among other factors SW spectral data, i.e., albedo as a function of wavelength, so that it can be calculated for different spectral distributions, caused e.g., by changing cloudiness or SZA. Third, it needs monitoring of the incident shortwave radiation at the site of interest and ideally its spectral distribution. To our knowledge, no study has simultaneously addressed these three challenges. Investigations of the spatial distribution of Arctic vegetation assemblages using medium- to very high-spectral-resolution remote sensing have been conducted (Putkiranta et al., 2024; Yang et al., 2023) but without addressing the SW energy budget issue. The SW energy budget change due to vegetation changes has been addressed (Aartsma et al., 2020; Beringer et al., 2005; Finne et al., 2023) but only for changes in specific assemblages, without determining changes in their spatial distribution or the SW energy budget over a given geographic area. While SW radiation monitoring datasets are available for the high Arctic (Boike et al., 2019; Domine et al., 2024) and at the FTE (Domine et al., 2024), they are not commonly integrated into comprehensive analyses of vegetation dynamics and their SW climate forcing.

To address these gaps in local studies, we combined field work, Landsat remote sensing data and existing radiation monitoring data near the Inuit community of Umiujaq on the Hudson Bay coast in north-eastern Quebec. This site, situated in the FTE, features a diverse mosaic of lichen tundra, shrub tundra, and open forest within less than 1 km (Lackner et al., 2021). Our actions consisted of (i) recording in the field the spectral albedos (347–2400 nm) of the dominant vegetation assemblages in the area. We collected 62 spectra, from which broadband albedo and NDVI values were calculated, and which were used to obtain spectral signatures that were used to determine the spatial distribution of vegetation assemblages using remote sensing: (ii) analyzing Landsat spectral measurements, providing information on areal changes in major vegetation types at high spatial resolution over a 40-year period. To achieve this, we developed machine learning classification models trained with ground-measured spectral albedos and applied probabilistic spectral unmixing methods to evaluate subpixel vegetation changes; (iii) measuring the spectral characteristics of downwelling SW radiation under various cloudiness and SZA conditions and combining

those with broadband SW radiation monitoring (Domine et al., 2024), the hyperspectral field spectroscopy obtained in (i) and with four decades of high-resolution unmixed vegetation mapping, to estimate SW climate forcing at a spatially comprehensive scale. Furthermore, to connect the mapping of SW forcings with the more commonly studied mapping of NDVI, we computed Landsat-based greening trends over the last four decades and compared their magnitude with vegetation transitions and radiative forcings. By integrating these approaches, this study provides novel insights into the spatial and temporal dynamics of vegetation changes and associated greening in one location of FTE and their implications for Arctic climate feedback.

## 2. Methods

Our analysis aims to estimate the radiative forcing at the scale of a valley located at the boreal-arctic transition over the last 40 years by mobilizing ground-based spectral data to inform satellite remote sensing. The procedure is structured in three parts, with (1) a detailed description of the ground-based spectral data, (2) our mapping of vegetation type trajectories and associated greening, and (3) the calculation of radiative forcing by combining satellite remote sensing and ground-based spectral data, with the complete workflow shown in Fig. 1.

### 2.1. Site description

The Inuit community of Umiujaq in Nunavik, Northern Quebec (56.55°N. 76.55°W), is on the east coast of Hudson Bay (Fig. 2). The greater area consists of basaltic plateaus eroded to form rocky cuestas with little vegetation and valleys where the soil consists of marine silty sediments and sand beaches raised by isostatic uplift after the melting of the Laurentide icesheet (Bhiry et al., 2011; Fortier et al., 2020). The main area of interest is Tasiapik Valley, east of Umiujaq, towards Tasiujaq Lake (Fig. 2). A secondary study site is near the coastal plain. The upper soil layer near the coast and in the upper part of Tasiapik valley is mostly sand, while silt is encountered in the lower part of the valley leading to Tasiujaq Lake (Lemieux et al., 2020). The vegetation covers consist mostly of lichen tundra, shrub tundra and stunted krummholz spruce (Gagnon et al., 2019). Lichens found most frequently in the present study were *Cladonia stellaris* Opiz, *C. rangiferina* L., *Alectoria ochroleuca* Hoffm. Massal, and *Stereocaulon paschale* L. The dominant shrub species is dwarf birch, *Betula glandulosa* Michx., which is colonizing the lichen tundra (Lemay et al., 2018) (Fig. 3). Other shrub species of interest to this study are willows (predominantly *Salix planifolia* Pursh.), blueberries *Vaccinium* sp. and *Rhododendron greenlandicum* Oeder. Black crowberries, *Empetrum nigrum* L., a non-erect species, was also observed. Mountain alder *Alnus viridis* (Chaix) D.C. ssp. *crispa* (Aiton) Turrill grows only along stream beds and in some humid hollows. Finally, black spruce, *Picea mariana* Mill. is present in stunted krummholz form that rarely exceeds 3 m in height, often much less (Fig. 3a, b). In the forest patch visible in Fig. 3b, however, spruce may rise to 4 m high.

### 2.2. Ground albedo measurements

Spectral albedo was measured with a field portable spectroradiometer (HR-1024 model, Spectra Vista Corporation) between 12 and 17 September 2015, between about 10:00 and 16:30 local time each day. Most of the time lighting conditions were 80 % to fully overcast, with short periods up to 60 % clear sky. Fully clear sky conditions were not encountered, so that most of our spectra were recorded under white sky conditions or close to that. Solar zenith angles (SZA) were between 52.4° and 66.5°, but given the mostly overcast conditions, these values had little impact. The light collector was an integrating sphere located at the end of a 2-m metal rod resting on a tripod (Fig. 3a) at a height of about 1.2 m. Since the sphere is a cosine receptor, 71 % of the radiance
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Flow chart illustrating the full procedure developed to investigate radiative forcings over the last four decades derived from the combination of satellite remote sensing and ground-based spectral data. Datasets are represented by rectangles and processing is represented by ovals. Green rectangles represent initial datasets while red rectangles represent main results of our analysis. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

entering the sphere comes from a circle 1.2 m in radius below the sphere and the remaining from the outer area. For measurements above 40 cm shrubs, this radius is reduced to 0.8 m. Albedo was determined from successive measurements of the downwelling and upwelling radiation, by rotating the metal rod 180° between both measurements. The horizontality of the sphere was checked with an electronic inclinometer attached next to the sphere. A photodiode monitored variations in the solar irradiance. Data where irradiance varied by >1 % between both measurements were rejected, as in (Belke-Brea et al., 2020). Light detectors were a Si photodiode and 2 InGaAs photodiodes. The spectral resolution varies between 1 and 1.5 nm for the Si detector (346.5–982 nm), 3.5 to 4 nm for the first InGaAs detector (982–1883 nm) and 2.5 nm for the second InGaAs detector (1883–2514 nm). For each spectrum, a running mean average was used to improve the signal/noise ratio. Averaging was over a 12 nm range from 346.5 to 1800 nm, and 80 nm beyond that. For wavelengths greater than 2400 nm, data quality was often poor and unreliable for albedo calculation, so that data over the 346.5–2400 nm range are presented. The impact of neglecting the 2400–2514 nm range is negligible because only 0.35 % of downwelling solar energy falls within that range and broadband albedos calculated by omitting this range is likely more reliable.

Broadband albedo, i.e., the ratio of total reflected SW energy over total incident SW energy, was calculated over the 300–2400 nm range. The contribution of the 300–346.5 nm range to albedo was estimated as follows. Based on the SMARTS atmospheric radiation model (http://www.nrel.gov/grid/solar-resource/smarts.html, last accessed on 6 June 2023), we evaluated that this range contributes 1.6 % to the total shortwave downwelling energy. The albedo of numerous vegetation types in the 300–340 nm range is about 90 % of that at 350 nm according

to (Feister and Grewe, 1995). The contribution of the 300–346.5 nm range was assessed based on the hypotheses that (i) this range contributes 1.6 % of incident irradiance; (ii) the average albedo over that range is 0.9 x albedo at 346.5 nm.

Albedo was measured over different surfaces dominated by (a) lichen; (b) spruce; (c) dwarf birch with lichen understory; (d) dwarf birch with moss understory; (e) willows; (f) varied prostrate vegetation including *Stereocaulon paschale*, *Empetrum nigrum*, mosses, *Carex* sp., and various shrubs <15 cm height, mostly *Betula glandulosa* and *Vaccinium* sp. (Fig. 4). No spectra of alders were recorded because of their low spatial coverage and because the location of small enough individuals to allow measurements were found in complex topographic environments not suitable for our setup. The locations of the measurements are indicated in Fig. 2. After each spectral measurement, photographs of the vegetation were taken to allow the assignment to a vegetation type.

To account for the variation in the spectral distribution of radiation with cloudiness, we calculated broadband albedo under two different lighting conditions (i) partially clear sky conditions, using as downwelling spectral irradiance distribution the average of five spectra measured on 17 September over lichen and spruce; (ii) fully overcast conditions, using as downwelling irradiance the average of six spectra measured on 12–13 September over lichen, birch with lichen, birch with moss and spruce. The spectral downwelling irradiances used are shown in Figure S1, showing that under overcast conditions the incident radiation is depleted in IR wavelengths relative to the partially blue-sky conditions The broadband albedo was calculated as:
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Map of the area near Umiujaq, Quebec, on the east coast of Hudson Bay with the 15 cm aerial image from 2010 (Provencher-Nolet et al., 2014) as a background. The approximate location of the various vegetation types measured are indicated. Two points of interest are shown, the Umiujaq village and the valley pass. The Tasiapik valley is depicted in yellow polygon. The centers of the 500 m × 500 m MODIS pixels used to compare white-sky and black-sky albedo products (section 2.4) are also shown. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Various vegetation types near Umiujaq. a) stunted spruce with lichen patches on 17 September 2015; b) overview of the lower part of the valley with Lake Tasiujaq on the right. Lichen tundra, dwarf birch and the spruce forest patch are visible, 10 September 2019; c) view of the pass, showing the mix between dwarf birch and lichen, and Hudson Bay in the background, 10 September 2019; d) dwarf birch patches expanding on lichen tundra, 19 August 2013. The spruce thicket is probably a relic from a warmer period a few hundred years ago (Payette and Morneau, 1993).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4. The six vegetation types investigated by spectroscopy and photographed while making measurements. a) lichen tundra. The main lichens visible here are Cladonia stellaris and Alectoria ochroleuca. Carex sp. grasses are also visible. b) Stunted black spruce. c) Dwarf birch tundra with lichen understory, mostly Cladonia stellaris. d) Dwarf birch tundra with moss understory. e) Willows. f) Varied prostrate vegetation, here Stereocaulon paschale, Empetrum nigrum, mosses, Carex sp., Vaccinium sp., Rhododendron groenlandicum and Betula glandulosa.</figcaption></figure>

$$Broadband\ albedo = \frac{\sum_{i} A(\lambda) I_{i}(\lambda_{i}) d\lambda}{\sum_{i} I_{i}(\lambda_{i}) d\lambda} \quad (1)$$

where A is the albedo at wavelength λ, as shown in Fig. 5a and I(λ) dλ is the incident energy over the wavelength range dλ as determined by the spectrometer resolution, as shown in Figure S1.

During our field measurements, fall colors of deciduous plants were starting to appear. Overall, we estimate that about 5 % of deciduous leaves had fall colors. Spruce and lichen were not affected.

## 2.3. NDVI calculations from ground spectral data

NDVI uses red (R) and near-infrared (NIR) spectral ranges. The literature reports the use of different spectral ranges for R and NIR, depending on the satellite used when NDVI is determined from remote sensing data (Berner et al., 2020; Blok et al., 2011; Juszak et al., 2014), or on the instrument used when it is determined from field data (Boelman et al., 2011; Jespersen et al., 2023; May et al., 2022). Here we

base our NDVI calculations on bands centered on the Landsat R and NIR bans, at 655 and 865 nm:

$$NDVI = (NIR_{865} - R_{655}) / (NIR_{865} + R_{655}) \quad (2)$$

Our spectra are smoothed over 12 nm, so that the actual ranges used are R(649–661) and NIR(859–871). In comparison, Landsat ranges are 630–680 nm for R and 845–885 nm for NIR, but wavelengths at the center of the band contribute more to the energy received, so that our ranges are very close to those of Landsat.

## 2.4. Effect of lighting conditions on broadband albedo

The spectra recorded in the field were mostly under overcast (i.e., white-sky) conditions. We use those data to estimate broadband albedo under both white-sky and blue-sky conditions. As mentioned above, the difference in spectral distribution of radiation between both conditions has been considered. However, albedo also depends on the direction of incoming radiation (Peltoniemi et al., 2005; Yang et al., 2020). It is diffuse (isotropic) under white-sky conditions while it is a mixture of
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Optical data for the main six vegetation types investigated. (a) Average spectral albedo for each vegetation type, with Landsat bandwidths along the X-axis; (b) broadband albedo; (c) NDVI. For (b) and (c), error bars are standard deviations. Lichen12, Lichen13 and Lichen17 refer to spectra measured on 12, 13 and 17 September with respective dominant species *Cladonia stellaris* and *C. rangiferina* (12), *C. stellaris* and *Alectoria ochroleuca* (13), and *C. stellaris* only (17). (b) and (c) also include the average values for all lichens measured. The spectrum was not included in (a) to avoid cluttering.</figcaption></figure>

diffuse and direct (directional) under blue-sky conditions. Contrary to Landsat, MODIS has a wide swath allowing multiple viewings at numerous angles, from which white-sky albedo and black-sky (fully direct) albedo were calculated. We obtained both these MODIS albedo products (MCD43D3) for four 500 m pixels in Tasiapik valley (Table S1 and Fig. 2). Two pixels had lichen and birch with lichen understory as dominant vegetations; one pixel was in the spruce forest patch; the last pixel was on rock to test the impact of rock outcrops possibly present in the other three pixels. Since black-sky albedo refers to 100 % direct illumination, the MODIS comparison informs us on an upper limit of the error considered.

## 2.5. Processing of Landsat surface reflectance data

We utilized all Tier 1 data available for our study site from 1984 to 2023 for the months June, July, August and September in Landsat Collection 2 provided by the U.S. Geological Survey (USGS) and hosted on Google Earth Engine (GEE). The Tier 1 data products analyzed include surface reflectance from Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper + (ETM+) and Landsat 8 Operational Land Imager (OLI). We exclusively chose images with an average cloud cover of less than 80 %, as scenes with high cloud cover can compromise the accuracy of geometric calibration. The C Version of Function of Mask (CFmask) was applied to categorize each pixel as clear (land/water), snow, cloud, adjacent to cloud, or cloud shadow (Zhu et al., 2015; Zhu and Woodcock, 2012). Pixels affected by snow, cloud, adjacent to cloud or cloud shadow were excluded from the analysis. Angular effects arising from variation in viewing and solar geometry throughout the Landsat time series have been identified as a significant source of variation in retrieved directional reflectance, independent of

ground-related changes (Nagol et al., 2015). Correcting for these effects is crucial due to the variability in the sun-surface sensor geometry, described as the Bidirectional Reflectance Distribution Function (BRDF), at the time of acquisition, which varies both spatially and temporally. (Roy et al., 2016) introduced a comprehensive method using a fixed set of parameters based on the RossThick-LiSparse BRDF model (Schaaf et al., 2002). This method enables the normalization of the entire Landsat time series to nadir (0° viewing zenith angle) and to a constant SZA, reducing BRDF-related variations in reflectance. We applied the correction directly in GEE using the script developed by Erik Lindquist and Daniel Wiell and available in Earth Engine Developers Groups. Our correction implementation involved the use of global coefficients and an optimal, normalized SZA set constant per location (Zhang et al., 2016). (Berner et al., 2020) identified systematic radiometric discrepancies among Landsat 5 TM, 7 ETM+, and 8 OLI satellites. They proposed a method based on phenological proximity to cross-calibrate Landsat 5 to Landsat 7, and Landsat 8 to Landsat 7, utilizing years when both satellites were available. Landsat 7 served as a benchmark due to its temporal overlap with the other two sensors. We performed cross-sensor calibration to correct for systematic radiometric discrepancies in Landsat 5 TM, 7 ETM+, and 8 OLI satellites spectral bands using the methods and coefficients for the Arctic of (Berner et al., 2023; Berner et al., 2020). We used this surface reflectance data for two applications. (1) We used all the data from 1984 to 2023 and from June, July, August and September to estimate greening trends in our study area as it is an efficient approach widely used to assess Arctic vegetation change. (2) We used data from 1984 to 1988 (Period 1) and from 2019 to 2023 (Period 2) to produce two vegetation classification maps that will allow us to obtain a more precise trajectory of vegetation types changes.
## 2.6. Diachronic vegetation classification and validation

We used Landsat surface reflectance and ground spectral data to produce two vegetation classification spanning four decades of vegetation change over our study site. Ground spectral data was used to produce “endmember” spectral signals for the six vegetation types. We computed the albedo over the Landsat 5 TM and 8 OLI bands in order to obtain the spectral signature of each vegetation type as seen by Landsat 5 TM and Landsat 8 OLI. We differentiated between the two sensors because only Landsat 5 TM images are available in the 80s (for P1) and only Landsat 8 OLI images are available in the 2020s (for P2). Based on vegetation type spectral signatures, we calibrated two random forest (RF) classifiers, one per sensor. These RF models were used to predict the distribution of vegetation types at the valley scale using Landsat bands reflectance values as predictive variables. We relied on the out-of-bag classification accuracy to select the best RF model. Landsat 5 TM BRDF and cross-sensor calibrated reflectance of each band (blue, green, red, NIR, SWIR1 and SWIR2) from summer of 1984 to 1988 were averaged to obtain one raster per band that is considered to be representative of the first period. The same procedure was applied to Landsat 8 BRDF and cross-sensor calibrated reflectance for the year 2019 to 2023. The classifications were not produced at an annual time step due to insufficient data at the beginning of the Landsat time series. We predicted the probability of each vegetation types at Landsat pixel scale using RF models and averaged reflectance values as predictive variables over the valley and attributed the vegetation types with the highest probability for the two periods. Vegetation trajectories were assessed using a Sankey diagram that allows to depict the type-by-type vegetation change.

We validated the P2 classification using photo-interpreted samples. We randomly sampled 500 points within the study area and photo-interpreted vegetation types based on a mosaic of 15 cm resolution color aerial photographs acquired during an overflight in 2010. The mosaic was produced from digital photographs, the protocol for which is detailed in (Provencher-Nolet et al., 2014). One of the co-authors (E.L.) not involved in the remote sensing-based vegetation classification procedure performed the photo-interpretation independently. Points were positioned on the very high-resolution image and each sample was assigned to one of six vegetation types. If the sample was located on a road, quarry, cast shadows or watercourse, the point was removed from the dataset. This corresponded to 30 samples over the 500 points. The remaining 470 samples were distributed as 13, 148, 59, 85, 117 and 48 for “Lichen”, “Birch-Lichen”, “Birch-Moss”, “Willow”, “Spruce” and “Varied low” vegetation types. We evaluated the quality of the P2 model by computing overall accuracy and Kappa, and balanced accuracy and F-score by vegetation types.

## 2.7. Robust greening trends estimates from Landsat

We computed NDVI using BRDF and cross-sensor calibrated Red and NIR bands and computed the annual maximum during the growing season. (Berner et al., 2020) demonstrated that $NDVI_{max}$ estimates are dependent on the number of observations during the growing season. (Bayle et al., 2024) showed that, because Landsat observations increase over the time series, the dependence of $NDVI_{max}$ on sampling frequency can lead to an overestimation of greening trends. (Berner et al., 2020) and (Berner et al., 2023) showed that this sampling bias can be partially corrected by modeling phenology on a pixel-by-pixel basis and adjusting the estimates of $NDVI_{max}$ based on mean phenology. Hence, we applied phenological modeling to adjust NDVI values before computing $NDVI_{max}$. We relied on the Harmonic Analysis of Time Series (HANTS) reconstruction method as used by (Choler et al., 2025). This procedure was only applied on NDVI observations with values higher than 0.15 to eliminate unvegetated pixels as the correction relies on phenological modeling. The complete procedure, relying on state-of-the-art methods allowed us to derive a robust annual $NDVI_{max}$ series for our study site. To

obtain slope and p-values estimates pixel-wise, we first determined the slope using a non-parametric Theil-Sen slope estimator (Sen, 1968) as implemented in the mblm package in R (Komsta, 2019) and performed a rank-based Mann-Kendall trend test (Kendall, 1975). Greening slopes and proportion of pixels exhibiting significant greening were compared among vegetation type of P1. This comparison allows to understand which vegetation changes is leading to the highest greening signal.

## 2.8. Estimating radiative forcings from probabilistic unmixing

The outcome of random forest models is not a class label but a vector of class membership probabilities which are afterward collapsed into a single label to be displayed as a map (as done in section 2.6). Despite the high resolution of Landsat images that minimize mixed classes within one pixel compared to, e.g., MODIS, mixed pixels are the norm in tundra vegetation (Assmann et al., 2020). Collapsing probabilities into single classes is convenient to understand and visualize changes in vegetation over time. However, if some classes are much more abundant than others, it is common that area estimates (and thus area change estimates) based on counting pixels overestimate abundant classes and underestimate scarce classes (Canters, 1997; de Bruin, 2000) which might have implications for radiative forcings estimates. It has been proven that random forest class probabilities are consistent estimates of the occurrence probabilities (Malley et al., 2012) and thus that the total area of each vegetation types can be estimated by summing of the pixels corresponding class membership probabilities multiplied by the pixel area (Sales et al., 2022). The average albedo of each pixel was calculated by weighting the albedo of the “pure” vegetation types by their respective proportions as determined by the random forest class probabilities (section 2.6). We use the downwelling shortwave radiation data of (Domine et al., 2024), obtained in the upper Tasiapik valley from 2013 to 2021 and found that a change in 0.01 of albedo resulted in a radiative forcings of $1.718 \text{ W m}^{-2}$. An example of these data is shown in Figure S2. Radiative forcings were computed for each pixel of the valley and the average summer radiative forcings due to vegetation type changes over the last four decades was estimated.

# 3. Results

## 3.1. Spectral albedo

Fig. 5a plots the averages of spectral albedo measured in the field for each vegetation type under mostly overcast (white-sky) conditions, since blue-sky conditions were not encountered. Lichens were found to have the highest albedos in the visible and in the short-wave infrared (SWIR, 1000 to 2500 nm) with values always exceeding 0.12 at 650 nm and sometimes exceeding 0.4 at 1670 nm, much greater than all other spectra. However, lichen albedos were very dependent on the dominant lichen species. On 17 September, in the middle part of Tasiapik valley, the species studied were dominated by *Cladonia stellaris* which has the highest albedo of lichens for which species data are available (Kuusinen et al., 2020). On 12 September, in the upper part of Tasiapik valley, the dominant species were *C. stellaris* and *C. rangiferina*, the latter species having a noticeably lower albedo (Kuusinen et al., 2020). On 13 September, near the coast, *Alectoria ochroleuca*, which in the visible appears almost as bright as *C. stellaris* (albedo of 0.375 over the 400–700 nm range vs 0.415 for *C. stellaris*, (Gauslaa, 1984; Rautiainen et al., 2024)), was also present. We therefore present lichen data separately for these three days in Fig. 5a. The other four vegetation types have fairly similar spectra and albedo values. The most notable difference is the greater albedo of willows in the 800–1100 nm range (NIR).

These spectral data were used to calculate broadband albedo and NDVI values of each vegetation types illustrated in Fig. 4. Values are listed in Table 1 with their respective standard deviations. Fig. 5b and c are graphical representations of Table 1, for easier visualization. The lichen vegetation type clearly stands apart from the other vegetation
<caption>**Table 1**</caption>

Broadband (BB) albedo and NDVI of the dominant six vegetation types. The number of spectra n and the standard deviations σ are indicated.

<table><thead><tr><th>Vegetation type</th><th>n</th><th>BB albedo<br>Clear sky</th><th>σ</th><th>BB albedo<br>Overcast</th><th>σ</th><th>NDVI</th><th>σ</th></tr></thead><tbody><tr><td>Lichen 12</td><td>9</td><td>0.200</td><td>0.009</td><td>0.191</td><td>0.009</td><td>0.396</td><td>0.034</td></tr><tr><td>Lichen 13</td><td>5</td><td>0.237</td><td>0.010</td><td>0.223</td><td>0.010</td><td>0.351</td><td>0.021</td></tr><tr><td>Lichen 17</td><td>4</td><td>0.303</td><td>0.032</td><td>0.292</td><td>0.031</td><td>0.386</td><td>0.009</td></tr><tr><td>Lichen all</td><td>18</td><td>0.233</td><td>0.048</td><td>0.222</td><td>0.047</td><td>0.381</td><td>0.035</td></tr><tr><td>Varied low</td><td>8</td><td>0.175</td><td>0.022</td><td>0.166</td><td>0.021</td><td>0.723</td><td>0.027</td></tr><tr><td>Birch-Lichen</td><td>7</td><td>0.168</td><td>0.019</td><td>0.162</td><td>0.019</td><td>0.672</td><td>0.037</td></tr><tr><td>Birch-Moss</td><td>9</td><td>0.165</td><td>0.024</td><td>0.159</td><td>0.023</td><td>0.757</td><td>0.020</td></tr><tr><td>Willow</td><td>8</td><td>0.183</td><td>0.014</td><td>0.179</td><td>0.014</td><td>0.801</td><td>0.009</td></tr><tr><td>Spruce</td><td>12</td><td>0.166</td><td>0.041</td><td>0.165</td><td>0.041</td><td>0.810</td><td>0.028</td></tr></tbody></table>

types, with higher broadband albedos and much lower NDVI values.

Lichen NDVI values are about half the values of all the other vegetation types (Table 1 and Fig. 5c). This is explained by its high red albedo, clearly visible in Fig. 5a, while its NIR albedo is similar to or lower than those of other vegetation types, expect for Lichen 17, but the red albedo is so high that it more than compensates for this to produce a low NDVI. This shows clearly in the Red-NIR plot of Figure S3, where lichen is well separated from all other vegetation, implying that its identification with remote sensing is easy and reliable.

Within each vegetation type, there are significant variations in albedo. All individual spectra are shown in supplementary material Figs. S4 to S9 and the variations within each vegetation type are also discussed in detail there, based on photographs taken along with each spectrum. Briefly, variations in spectral albedo are simply explained by variations in vegetation assemblages and structure. For example, variations in albedo for erect vegetation is explained by variations in canopy structure, with the thickest canopies having the highest albedos. This latter observation is consistent with those in the Canadian low Arctic of (Blanken and Rouse, 1994) and of (Williamson et al., 2016) who reported that for willow and dwarf birch, the highest albedo measured were for the thickest canopies with the highest leaf area index (LAI). Overall, there is therefore a good correspondence between visual vegetation assemblages and structure and spectral data.

## 3.2. Vegetation dynamics

Changes in vegetation types over the last four decades were investigated by producing two vegetation maps using satellite imagery and ground spectral data. We validated the Random Forest-based vegetation map of the second period using 470 randomly sampled photo-interpreted points based on very high-resolution images. We found an overall accuracy of 0.77 [95 % CI: 0.73–0.81] and a Kappa of 0.71 suggesting good quality in predicting vegetation types at the valley scale. The F-score was greater than 0.6 for all vegetation types with up to 0.861 and 0.828 for Spruce and Birch-Lichen respectively (Table S3). The lowest values were found for Birch-Moss and Varied low types, 0.603 and 0.673 respectively. In general, specificity was found to be high for all types (>0.9) while sensitivity was intermediate for Lichen (0.600), Birch-Moss (0.532) and varied low (0.641), suggesting, for example, that some pixels classified as Lichen should have been classified as Birch-Lichen, but pixels that should be classified as Lichen are systematically classified as Lichen. Hence, the overall surface of Lichen, Birch-Moss and Varied low are probably slightly overestimated in our classification. Overall, our random forests provide high quality estimates of vegetation types distribution at the valley scale.

Fig. 6 shows that during P1 (1984–1988), Lichen occupied a total of 2.33 km², corresponding to 25.2 % of the Tasiapik valley, while it only covered 0.36 km², 3.9 %, during P2 (2019–2023). Most of the Lichen was converted to the Birch-Lichen class and in very limited proportion to the varied low class. Overall, the Birch-Moss, Spruce, Willow and Varied low classes exhibit increases in areal cover over the last four decades while Lichen and Birch-lichen exhibit a decrease (Fig. 6). The trajectories of vegetation types are spatially structured (Fig. 6). For example, the Spruce forested areas in P1 show clear lateral expansion for both patches located in the middle and lower parts of the valley. Most transitions from Lichen to Birch-Lichen classes occurred near the pass (NW), with large patches of lichen present during P1 now showing spectral characteristics of Birch-Lichen, indicating encroachment of lichen by birch. Most of the Birch-Moss increase comes from transformation of Birch-Lichen areas. This is explained by the growth of the birch shrubs, which increase soil humidity and decrease irradiance at ground level, favoring moss at the expense of lichen (Gagnon et al., 2019). Classification probabilities for the two periods are shown in Figure S11.

In addition, we investigated changes in satellite-derived greenness over the last decades and its correspondence with changes in vegetation types (Fig. 7). Examples of greening trends plots over the 1984–2023

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** (left) Distribution of the main vegetation assemblages in period 1 (1984–1988, top) and period 2 (2019–2023, bottom). (right) Sankey diagram of the trajectories of the main vegetation assemblages between both periods. The surfaces for each vegetation assemblages are indicated. The valley is 9.25 km², excluding anthropized areas (roads and quarry). A map of vegetation changes is shown in Figure S10.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Greening in Tasiapik valley. (a) Landsat-based greening trends from 1984 to 2023 based on NDVI<sub>max</sub> with histogram depicting the distribution of all pixels independently of their significance. (b) Bar graphs showing the distribution of greening trends as a function of initial vegetation type, during P1. (c) Pie charts showing the proportion of significant greening trends (p-value <0.05) for each initial vegetation class, during P1.</figcaption></figure>

period are shown in Figure S12. Over the Tasiapik valley, 73 % of the pixels exhibit a significant greening trend (p-value <0.05). Greening trends are heterogeneous in their spatial distribution (Fig. 7a). Greening trends show large variations, depending mostly on the initial vegetation type during P1 (Fig. 7b). The strongest greening is found for Lichen (0.004 NDVI/year), and then for Birch-Lichen (0.003 NDVI/year). Almost all lichen pixels show significant greening trends (95 %), while it only concerns 40 % of the spruce class (Fig. 7c). Fig. 6 and Fig. 7 both lead to the conclusion that the main vegetation change is the

replacement of Lichen by Birch-Lichen. On-site observations over the 2012–2023 period show that this proceeds by the lateral expansion of existing bushes by clonal growth, clonal growth of individuals up to 2 m away from the edge of the main bush, or new plants by seed germination, identified by the absence of roots linking them to other individuals (Fig. 8).

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Birch expansion over lichen tundra. (a) laterally expanding bushes; (b) infilling by clonal growth with new individuals visible in foreground; (c) small plant from seed germination; (d) larger plant by seed germination. Photographs taken near the pass at the top of Tasiapik valley on 2 September 2018.</figcaption></figure>
### 3.3. Radiative forcing

The numerous types of vegetation changes shown in Fig. 6 lead to changes in spectral and broadband albedo. Fig. 5 shows that lichen replacement by birch has the strongest impact on albedo and we quantify here the resulting shortwave radiative forcing of this process as an example. This forcing is relevant only for the snow-free period, typically early June to late October. Our spectral measurements may apply neither before leaf-out nor to the fall when leaves change color and start falling, so we compute the forcing for the 15 June to 14 September period, hereafter called summer, for the years 2013 to 2021. We use the downwelling shortwave radiation data of (Domine et al., 2024), obtained in the upper Tasiapik valley from 2013 to 2021. An example of these data is shown in Figure S2. We use a broadband albedo of 0.233 for lichen (average of all values) and 0.168 for birch with lichen understory (Table 1). The calculated average summer forcing is $11.17 \pm 0.70 \text{ W m}^{-2}$ for the Lichen to Birch-Lichen transition when considering pure pixels. The error bar is the standard deviation in yearly values. The value increases to $11.69 \text{ W m}^{-2}$ if birches keep growing, as lichen is then replaced by moss (Gagnon et al., 2019) and the broadband albedo decreases to 0.165 (Table 1).

Based on spectral unmixing, we estimated per pixel albedo for the two periods and computed the difference. By converting the difference in albedo to summer radiative forcing, we obtained the spatial distribution of the forcing (Fig. 9). At the valley scale the mean forcing is $2.16 \text{ W m}^{-2}$ over the summer. The strongest positive forcing is at the top of the valley, where Birch-Lichen replaces Lichen. The strongest negative forcing is in humid areas, where Willows expand in height and laterally at the expense of Birch-Lichen and Birch-Moss, with a forcing around $-3 \text{ W m}^{-2}$ (Table S2). In terms of energy, the average daily summer forcing for the Lichen to Birch-Lichen transition is $0.47 \text{ MJ m}^{-2} \text{ d}^{-1}$. Averaged over the valley, this value drops to $0.10 \text{ MJ m}^{-2} \text{ d}^{-1}$. At valley scale, 76% of the valley exhibit positive radiative forcings.

It is interesting to explore the correlation between greening trends and radiative forcing values, regrouped by intervals of $2 \text{ W m}^{-2}$ (Figure S13). On average, greening trends were positive across all radiative forcing classes, with some negative values, also visible in Fig. 7, occurring for radiative forcing values between $-3$ and $+3 \text{ W m}^{-2}$. The lowest greening trends correspond to radiative forcings near 0, representing a lack of vegetation transition between the two periods. Negative radiative forcings (i.e., increases in albedo) mostly concerned pixels with enhanced willow cover, with concomitant greening and increased albedo (Figs. 7 and 9), in line with the observations of

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Shortwave radiative forcing for each pixel in the valley between both periods. The histogram of values and the mean valley-scale value are also shown. The pie chart shows the proportion of pixels exhibiting positive and negative radiative forcings.</figcaption></figure>

(Blanken and Rouse, 1994) and (Williamson et al., 2016), as mentioned in section 3.1. Conversely, transitions resulting in positive radiative forcings showed a positive correlation between greening and radiative forcings, with variations depending on the specific vegetation transition. The Lichen-to-Birch-Lichen shift stands out, exhibiting the highest positive radiative forcing ($11.17 \text{ W m}^{-2}$) and the strongest greening trend ($>0.005 \text{ NDVI y}^{-1}$).

### 3.4. White-sky vs black-sky albedo

The differences in MODIS white-sky and black-sky albedo products for the four pixels shown in Fig. 2 are reported in Table S1 for summer 2015 (15 June to 14 September), when our field measurements were performed. Table S1 shows that the average albedo difference is $<0.01$ for all four pixels investigated, and even as low as 0.0035 for the pixel Lichen 1. The average for the three vegetated pixels is 0.0069. This again is an upper limit because black-sky conditions are fully directional. Given the fairly high latitude of Umiujaq (56.55°N) and the resulting large SZA, blue-sky conditions have a large fraction of diffuse radiation. The (white-sky) – (blue-sky) difference can reasonably be expected to be significantly less than the values of Table S1, and we estimate that they are probably close to half this value, i.e., 0.0035.

## 4. Discussion

### 4.1. Implications

The consequences of shrubification extend far beyond local effects, influencing broader systems such as sea ice coverage through enhanced evapotranspiration (Chen et al., 2022). Accurately quantifying these effects necessitates detailed investigations. In this study, we combined spectral data from satellite remote sensing with ground measurements to examine not only the impact of shrubification but the broader vegetation dynamics occurring within the forest-tundra ecotone. Our analysis revealed an albedo change of approximately 0.065 for the Lichen to Birch-Lichen transition, corresponding to a summer shortwave forcing of $11.17 \text{ W m}^{-2}$. When all vegetation changes are considered comprehensively at the valley scale, the shortwave forcing is reduced to $2.16 \text{ W m}^{-2}$. This finding was possible by adopting a holistic approach for understanding vegetation-induced radiative forcings, which is only achievable through the integration of high-resolution satellite remote sensing and advanced spectral unmixing techniques.

Other studies have investigated the shortwave forcing due to shrubification. (Aartsma et al., 2020) measured the broadband albedo of ten lichen plots and of another nearby ten shrub plots dominated by *Betula nana* in southern Norway. The average lichen albedo was 0.255, close to our 0.233 value, but their average shrub value, 0.132, was much lower than our shrub values, 0.165 to 0.183. Their albedo difference was thus 0.124, almost twice our 0.065 value, so that the resulting shortwave forcing was bound to be much greater. The work of (Finne et al., 2023) confirms that vegetation changes in Norway have a greater impact on the surface energy budget than in Northern Quebec. These authors compared shrub and lichen albedos in southern and northern Norway. The classification they used is different to ours as they mixed data with values from Svalbard, making detailed comparison difficult. However, they also had solar irradiance monitoring at their sites and calculate that the mean summer radiative forcing was $1.61 \text{ MJ m}^{-2} \text{ d}^{-1}$ in southern Norway and $1.28 \text{ MJ m}^{-2} \text{ d}^{-1}$ in northern Norway, when pale ecorticate lichen such as *Cladonia stellaris* was replaced by shrubs dominated by *Betula nana*. Our value for the Lichen to Birch-Lichen transition is 0.47 MJ m$^{-2}$ d$^{-1}$. In their case, the difference with our values seems to be more due to a higher lichen albedo than to a low shrub albedo, since all their lichen broadband albedo values were above 0.3 and they even report a high value of 0.389. This is in fact not totally surprising as their dominant lichen species was *C. stellaris*. (Reinhardt et al., 2022) measured experimental albedos on non-natural assemblages from pure
lichen and noted the high broadband albedo of *C. stellaris*, 0.364 ± 0.019. In both these studies the shortwave forcing was much greater than our value, showing that forcings due to similar vegetation changes can show large spatial variations, stressing the need for numerous studies in different locations.

Albedo effects of vegetation changes are relevant to climate modeling. (Chen et al., 2022) used an Earth system model to compare Arctic climate between the mid-Holocene and the pre-industrial period. Vegetation expansion in the mid-Holocene led to a surface albedo reduction which produced a forcing of 9.13 W m⁻² in spring and summer. This value is close to ours for the complete lichen replacement by shrubs but much higher than our mean valley value, considering the observed replacement over 40 years. A direct comparison is however difficult in the absence of details on species.

The forcing values we measured therefore appear to be lower than other available values. This suggests large spatial variability on the impact of vegetation change, at least between Norway and Northern Quebec. The need for similar studies in a large variety of environments must therefore be promoted for accurately quantifying the effects of vegetation change in Arctic and boreal regions. We however stress that our study allowed results at the valley scale while previous ones rather focused on a specific replacement only and their values may not be extrapolated to a larger scale. Our results were possible thanks to the use of detailed hyperspectral data on vegetation assemblages measured in the field, combined with the use of visible, NIR and SWIR band data provided by satellite remote sensing. This allowed the identification of vegetation types proportions for each pixel and the calculation of the shortwave forcing at the valley scale. Such a method seems to be rather unique. A similar classification in Tasiapik valley was performed using color photographs obtained with an airplane and using just three spectral bands (red, blue, green) (Provencher-Nolet et al., 2013). This allowed the determination of 4 classes where the dominant vegetation types lichen, herbs, shrubs and spruce were confirmed with field validation. Maps were obtained in 1994 and 2010 covering a 16-year interval. Their maps are overall quite consistent with ours, although their shorter time intervals inevitably produce differences with our Fig. 6 and the changes they show are less intense. As in this study, the main change they observed was the replacement of lichen by shrubs. However, they were not able to differentiate between shrub types as was possible with the addition of NIR and SWIR bands. Their study was focused on vegetation changes and did not attempt to quantify the resulting shortwave forcing. Nevertheless, it is reassuring to note the proximity of the trajectories identified in both studies.

Interestingly, the vegetation transition that results in the highest shortwave forcings is the one producing the highest greening magnitude (Fig. 7, Fig. 9 and Figure S13 Greening is first and foremost a relative measure of the evolution of plant cover and is highly dependent on the initial state. In the case of dwarf birch encroachment on lichen, vegetation with a particularly low NDVI (0.379 ± 0.036) is gradually replaced by high NDVI vegetation (0.686 ± 0.040). Conversely, the transition from Birch-Lichen (0.686 ± 0.040) to Spruce (0.829 ± 0.028) necessarily results in lower greening. The sensitivity of greening to shrub encroachment on lichen that we demonstrate here could explain that the forest-tundra ecotone is showing high-magnitude greening (Berner and Goetz, 2022; Ju and Masek, 2016) compared to other regions. Our study hence suggests that greening computed from satellite remote sensing can inform surface albedo changes and radiative forcings at the forest-tundra ecotone. Nonetheless, we show that the forcing that results from shrub encroachment on lichen is very dependent on the lichen species involved (Fig. 5) while NDVI appears similar among lichen species. Shrub encroachment onto various lichen species would therefore result in similar greening magnitude but different radiative forcings, depending on the lichen species. Mapping of lichen at species scale would improve modeling of radiative forcings in relation to vegetation dynamics. Also, our vegetation types probabilities show significant spatial heterogeneity, with uniform vegetation patches extending over

small areas (Figure S11). We show that vegetation dynamics (greening) and their consequences (radiative forcing) are also played out at these fine scales, as suggested by (Assmann et al., 2020). Our analysis therefore demonstrates the importance of using decametric-scale satellite remote sensing associated with spectral unmixing approaches to capture ongoing environmental changes in Arctic ecosystems in response to climate change.

Nonetheless, opposite patterns were also found with both high greening trends and negative radiative forcings (Figure S13). This corresponds to either the growth of existing willows or to the expansion of willows at the expense of other shrubs (Fig. 6), which have a lower albedo than willow (Fig. 5b). This illustrates that, while increased vegetation greenness is typically linked to positive radiative forcings (Figure S13), certain vegetation transitions can lead to both increased greenness and negative radiative forcings. This underscores the importance of considering species-specific broadband albedo in addition to greenness for investigating radiative feedback due to Arctic shrubification.

In summary, the impact of shrubification on the surface SW radiation budget is important and in summer appears highly spatially variable. The large-scale impacts of shrubification are potentially enormous as it generates long-distance feedbacks. Being able to accurately quantify the shortwave forcing is therefore one essential aspect. Since shortwave forcing is highly dependant on the change in species, the identification of these species must be performed. We have shown here that high-resolution hyperspectral data over the whole SW range has the potential to allow detailed species identification over large areas, using all the optical bands of satellites. This furthermore allows the determination of the shortwave radiative forcing for each pixel. The ground determination of the optical characteristics of vegetation, combined with satellite remote sensing, therefore offers an excellent potential for the determination of shortwave radiative forcings due to vegetation change. An improved accuracy may be expected with a refined identification of lichen species.

## 4.2. Limitations and sources of error

Our work relies mainly on our ability to correctly measure the spectral albedo of vegetation types representative of the valley and the changes in the areas covered by these same vegetation types over the last 40 years. This combination of spectral data from outer space and the ground is necessarily imperfect, and we propose to discuss some of its limitations and sources of error.

This study relies on spectral albedo measurements from 62 vegetation assemblages, selected to represent the diverse vegetation types across the valley. While every effort was made to capture the variability in key vegetation categories, such as lichen assemblages—which play a critical role in determining valley-scale radiative forcing—this representation inevitably has limitations. For lichens, we specifically sampled both the brightest and darkest types (Figs. 5 and S3) to encompass the full range of variability. With 18 lichen surfaces measured, statistical calculations suggest that this sample size (Figure S4) is sufficient to capture the observed albedo range of a factor of 1.8. Although this provides a strong foundation, uncertainties inherent in such sampling remain. In calculating radiative forcing, we aggregated the albedo values of various lichen species and used the average as a reference value. This aggregation introduces potential biases, as the dominant species present across the valley could lead to either over- or underestimation of radiative forcing. The broadband albedo of lichens varies significantly, from 0.186 in mixed mats of *C. stellaris* and *C. rangiferina* mats with sedge, cranberry shrub, and moss cover to 0.342 in almost pure *C. stellaris* patches (Figure S4). This variability results in substantial differences in radiative forcing estimates: assuming the valley would be entirely covered by *C. stellaris* mats, the valley-scale radiative forcing would reach 5.93 W m⁻². Conversely, assuming complete coverage by mixed mats yields a much lower radiative forcing of 0.53 W m⁻² Based on
photographs taken during ground spectral data sampling, our aggregated spectra approach aligns well with the actual composition of vegetation in Tasiapik valley. However, this method may limit the extrapolation of our findings to larger spatial scales, where the regional distribution of lichen species may differ significantly. A more detailed accounting of species composition at broader scales would be necessary to test the wide-scale applicability of these results.

Another source of uncertainty lies in the measurement of spruce albedo. Due to logistical constraints, we measured short spruce trees, which display notable albedo variations (Figure S9). Most spruce in the valley, however, are taller than those sampled. Despite this, short spruces are typically younger and are significant contributors to spruce expansion and related albedo changes over time. Therefore, while the omission of tall spruce may limit the scope of our measurements, it is unlikely to significantly affect our conclusions regarding radiative forcing trends. The seasonal timing of measurements also presents a potential limitation. Data collection occurred in late summer when approximately 5 % of deciduous leaves had already started to turn yellow. This could have slightly affected reflectance measurements. However, evergreen species such as spruce and lichen remained unaffected. Furthermore, previous studies (e.g., (Cipar et al., 2008)) have demonstrated the relative stability of optical properties and in particular indices such as NDVI throughout the summer-to-fall transition. Given the small proportion of deciduous leaves impacted and the early timing of our work, the effect on our results is likely minimal and can reasonably be neglected.

Lighting conditions during measurements introduce an additional source of uncertainty. The SW forcing values we calculated were based on “clear-sky” (in fact mostly but not purely blue-sky) albedo values, which differ from white-sky conditions, frequent near Umiujaq. Difference in albedo values arise from variations in the spectral distribution of radiation (Figure S1) and directional illumination (Peltoniemi et al., 2005; Yang et al., 2020). The spectral distribution used in our forcing calculation (2.16 W m⁻²) corresponds to a blue-sky scenario with some white-sky influence, which likely reflects average summer conditions at Umiujaq. Errors from this assumption are expected to be small. However, angular distribution differences also warrant consideration. For the three MODIS vegetation pixels analyzed (Table S1), the average albedo difference between white-sky and black-sky conditions is 0.0069. Given that blue-sky is a mix between black-sky and white-sky, and that our actual summer conditions, are a mix between blue-sky and white-sky, we estimate that the error we make in our sky conditions are at most half this 0.0069 value. Applying this correction would reduce the albedo difference between lichen and birch-lichen (0.065) by approximately 9 %, directly impacting the radiative forcing calculation. Since this process is the primary driver of the forcing estimate, we assign a 9 % uncertainty to this value due to unaccounted illumination variations.

To estimate the change in cover of the different vegetation types we calibrated Random Forest models and employed a probabilistic spectral unmixing method, thus limiting the biases inherent in a pixel-counting approach (Sales et al., 2022). Nevertheless, our approach is based on the Landsat time series which, while offering a better spatial resolution more suited to forest-tundra ecotones (Assmann et al., 2020), is made up of different sensors that can hamper a good estimate of vegetation dynamics. Although we have deployed methods to limit uncertainties (Roy et al., 2016), we have not directly validated the 1:1 correspondence between ground and satellite spectral measurements. To do so, we would have had to identify pixels with homogeneous vegetation over 30 × 30 m or carry out several spectral ground measurements within one single Landsat pixel (see (Wong et al., 2024) for similar conceptualization). Such an approach represents an interesting prospect for improving the combined use of ground and satellite measurements and provides assessment of radiative forcings over larger areas.

Finally, the quality of a land cover classification is highly dependent on its premises and promises. This means that the choice of classes we wish to map (*promises*) must be in line with those that our training data

capture (*promises*). In our case, we considered classes as close as possible to the reality of the spectral measurements made on the ground, but this necessarily comes with limitations. As our validation results show, there is an inequality of quality, with, for example, the Varied low vegetation type with the lowest sensitivities but high specificities (Table S3). This can be explained by the diverse floristic assemblage that makes up this class, accentuating spectral diversity and class variance. The vegetation types Lichen, Birch-Lichen, Spruce and Willow are more resolved, giving better results. Increasing the diversity of vegetation measured spectrally on the ground would improve vegetation classification and perhaps capture finer-grained vegetation dynamics.

Overall, all these limitations demonstrate the inherent difficulty of measuring radiative forcing linked to vegetation dynamics in these heterogeneous environments, where changes are played out on fine scales, and the importance of controlled field measurements to inform satellite images. Ultimately, our study underscores the value of integrating high-resolution field hyperspectral data with satellite remote sensing to bridge critical gaps in understanding the spatial variability and climate implications of vegetation dynamics, paving the way for more accurate quantifications of radiative forcings and their incorporation into broader climate models.

## Author contribution

FD and AB designed research. FD obtained funding. MBB performed field measurements with assistance from FD. FD, AB and MBB analyzed data. FD and AB wrote the paper. EL and GP provided comments and expert advice.

## CRediT authorship contribution statement

**Florent Dominé:** Writing – original draft, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Arthur Bayle:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Maria Belke-Brea:** Writing – review & editing, Methodology, Investigation, Data curation.
**Esther Lévesque:** Writing – review & editing, Methodology, Investigation, Data curation.
**Ghislain Picard:** Writing – review & editing, Methodology, Investigation.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This work was funded by the BNP-Paribas foundation (APT project to FD), the Natural Sciences and Engineering Research Council of Canada (Discovery Grant and Northern Research supplement programs to FD) and the French Polar Institute (IPEV, program 1042 to FD). AB acknowledges a grant from MELCCFP, Québec, Canada and its hosting at Université Grenoble Alpes (UGA), Laboratoire d'Écologie Alpine (LECA), Grenoble, France. We thank Mathieu Barrère for help with field measurements. We thank the community of Umiujaq for their warm welcome. We are grateful to the Centre d'Études Nordiques for the use of their research base at Umiujaq.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.rse.2025.114715.
## Data availability

The data obtained in this study will be made available on the Zenodo repository upon paper acceptance. For review purposes the data of Fig. 5a and of Figure S1 are available at: https://www.dropbox.com/scl/fi/b17ljk7hcn1tw3uvptaj/ZenodoFile.txt?rlkey=gflg470ry78x5m5y175hy7bq3&dl=0

The data of Fig. 6 (Rasters of land cover), Fig. 7 (Rasters of greening slope and p-values) and Fig. 9 (Raster of radiative forcings) are available at: https://zenodo.org/records/14998695.

## References

Aartsma, P., Asplund, J., Odland, A., Reinhardt, S., Renssen, H., 2020. Surface albedo of alpine lichen heaths and shrub vegetation. Arct. Antarct. Alp. Res. 52 (1), 312–322. https://doi.org/10.1080/15230430.2020.1778890.

Assmann, J.J., Myers-Smith, I.H., Kerby, J.T., Cunliffe, A.M., Daskalova, G.N., 2020. Drone data reveal heterogeneity in tundra greenness and phenology not captured by satellites. Environ. Res. Lett. 15 (12), 125002. https://doi.org/10.1088/1748-9326/abbf7d.

Bayle, A., Roy, A., Dedieu, J.-P., Boudreau, S., Choler, P., Lévesque, E., 2022. Two distinct waves of greening in northeastern Canada: summer warming does not tell the whole story. Environ. Res. Lett. 17 (6), 064051. https://doi.org/10.1088/1748-9326/ac74d6.

Bayle, A., Gascoin, S., Berner, L.T., Choler, P., 2024. Landsat-based greening trends in alpine ecosystems are inflated by multidecadal increases in summer observations. Ecography n/a(n/a), e07394. https://doi.org/10.1111/ecog.07394.

Beck, I., Ludwig, R., Bernier, M., Levesque, E., Boike, J., 2015. Assessing permafrost degradation and land cover changes (1986-2009) using remote sensing data over Umiujaq, Sub-Arctic Quebec. Permaf. Periglac. Process. 26 (2), 129–141. https://doi.org/10.1002/ppp.1839.

Belke-Brea, M., Domine, F., Barrere, M., Picard, G., Arnaud, L., 2020. Impact of Shrubs on Winter Surface Albedo and Snow Specific Surface Area at a Low Arctic Site: In Situ Measurements and Simulations. J. Clim. 33 (2), 597–609. https://doi.org/10.1175/jcli-19-0138.1.

Beringer, J., Chapin, F.S., Thompson, C.C., McGuire, A.D., 2005. Surface energy exchanges along a tundra-forest transition and feedbacks to climate. Agric. For. Meteorol. 131 (3), 143–161. https://doi.org/10.1016/j.agrformet.2005.05.006.

Berner, L.T., Goetz, S.J., 2022. Satellite observations document trends consistent with a boreal forest biome shift. Glob. Chang. Biol. 28 (10), 3275–3292. https://doi.org/10.1111/gcb.16121.

Berner, L.T., Massey, R., Jantz, P., Forbes, B.C., Macias-Fauria, M., Myers-Smith, I., Kumpula, T., Gauthier, G., Andreu-Hayles, L., Gaglioti, B.V., Burns, P., Zetterberg, P., D'Arrigo, R., Goetz, S.J., 2020. Summer warming explains widespread but not uniform greening in the Arctic tundra biome. Nat. Commun. 11 (1), 4621. https://doi.org/10.1038/s41467-020-18479-5.

Berner, L.T., Assmann, J.J., Normand, S., Goetz, S.J., 2023. LandsatTS': an R package to facilitate retrieval, cleaning, cross-calibration, and phenological modeling of Landsat time series data. Ecography 2023 (9), e06768. https://doi.org/10.1111/ecog.06768.

Bhiry, N., Delwaide, A., Allard, M., Begin, Y., Filion, L., Lavoie, M., Nozais, C., Payette, S., Pienitz, R., Saulnier-Talbot, E., Vincent, W.F., 2011. Environmental change in the great White River region, Hudson Bay: five decades of multidisciplinary research by Centre d'études nordiques (CEN). Ecoscience 18 (3), 182–203. https://doi.org/10.2980/18-3-3469.

Blanken, P.D., Rouse, W.R., 1994. The role of willow birch forest in the surface-energy balance at arctic treeline. Arct. Alp. Res. 26 (4), 403–411. https://doi.org/10.2307/1551802.

Blok, D., Schaepman-Strub, G., Bartholomeus, H., Heijmans, M.M.P.D., Maximov, T.C., Berendse, F., 2011. The response of Arctic vegetation to the summer climate: relation between shrub cover, NDVI, surface albedo and temperature. Environ. Res. Lett. 6 (3), 035502. https://doi.org/10.1088/1748-9326/6/3/035502.

Boelman, N.T., Gough, L., McLaren, J.R., Greaves, H., 2011. Does NDVI reflect variation in the structural attributes associated with increasing shrub dominance in arctic tundra? Environ. Res. Lett. 6 (3), 035501. https://doi.org/10.1088/1748-9326/6/3/035501.

Boike, J., Nitzbon, J., Anders, K., Grigoriev, M., Bolshyanov, D., Langer, M., Lange, S., Bornemann, N., Morgenstern, A., Schreiber, P., Wille, C., Chadburn, S., Gouttevin, L., Burke, E., Kutzbach, L., 2019. A 16-year record (2002–2017) of permafrost, active-layer, and meteorological conditions in the Samoylov Island Arctic permafrost research site, Lena River delta, northern Siberia: an opportunity to validate remote-sensing data and land surface, snow, and permafrost models. Earth Syst. Sci. Data 11, 261–299. https://doi.org/10.5194/essd-11-261-2019.

Bonfils, C.J.W., Phillips, T.J., Lawrence, D.M., Cameron-Smith, P., Riley, W.J., Subin, Z. M., 2012. On the influence of shrub height and expansion on northern high latitude climate. Environ. Res. Lett. 7 (1), 015503. https://doi.org/10.1088/1748-9326/7/1/015503.

Canters, F. (1997). Evaluating the uncertainty of area estimates derived from fuzzy land-cover classification. Photogramm. Eng. Remote. Sens., 63(4), 403–414. <Go to ISI://WOS:A1997WR92400011.

Chen, J., Zhang, Q., Kjellström, E., Lu, Z.Y., Chen, F.H., 2022. The contribution of vegetation-climate feedback and Resultant Sea ice loss to amplified Arctic warming

during the mid-Holocene. Geophys. Res. Lett. 49 (18). https://doi.org/10.1029/2022GL098816. Article e2022GL098816.

Choler, P., Bayle, A., Fort, N., Gascoin, S., 2025. Waning snowfields have transformed into hotspots of greening within the alpine zone. Nat. Clim. Chang. 15, 80–85. https://doi.org/10.1038/s41558-024-02177-x.

Ciernowski, J., Ceglarek, J., Karnieli, A., Królewicz, S., Każmierowski, C., Zagajewski, B., 2017. Predicting the diurnal blue-sky albedo of soils using their laboratory reflectance spectra and roughness indices. J. Quant. Spectrosc. Radiat. Transf. 200, 25–31. https://doi.org/10.1016/j.jqsrt.2017.05.033.

Cipar, J., Cooley, T., Lockwood, R., 2008, Aug 13. Summer to autumn changes in vegetation spectral indices of deciduous trees. Proceedings of SPIE [remote sensing and modeling of ecosystems for sustainability v]. In: Conference on Remote Sensing and Modeling of Ecosystems for Sustainability, San Diego, CA.

Crichton, K.A., Anderson, K., Charman, D.J., Gallego-Sala, A., 2022. Seasonal climate drivers of peak NDVI in a series of Arctic peatlands. Sci. Total Environ. 838, 156419. https://doi.org/10.1016/j.scitotenv.2022.156419.

de Bruin, S., 2000. Predicting the areal extent of land-cover types using classified imagery and Geostatistics. Remote Sens. Environ. 74 (3), 387–396. https://doi.org/10.1016/S0034-4257(00)00132-2.

Delbart, N., Kergoat, L., Le Toan, T., Lhermitte, J., Picard, G., 2005. Determination of phenological dates in boreal regions using normalized difference water index. Remote Sens. Environ. 97 (1), 26–38. https://doi.org/10.1016/j.rse.2005.03.011.

Dobbert, S., Pape, R., Löffler, J., 2021. How does spatial heterogeneity affect inter- and intraspecific growth patterns in tundra shrubs? J. Ecol. 109 (12), 4115–4131. https://doi.org/10.1111/1365-2745.13784.

Domine, F., Sarrazin, D., Nadeau, D.F., Lackner, G., Belke-Brea, M., 2024. Meteorological, snow and soil data, CO2, water and energy fluxes from a low-Arctic valley of northern Quebec. Earth Syst. Sci. Data 16 (3), 1523–1541. https://doi.org/10.5194/essd-16-1523-2024.

Edwards, R., Treitz, P., 2017. Vegetation greening trends at two sites in the Canadian Arctic: 1984–2015. Arct. Antarct. Alp. Res. 49 (4), 601–619. https://doi.org/10.1657/AAR0016-075.

Feister, U., Grewe, R., 1995. Spectral albedo measurements in the UV and visible region over different types of surfaces. Photochem. Photobiol. 62 (4), 736–744. https://doi.org/10.1111/j.1751-1097.1995.tb08723.x.

Finne, E.A., Bjerke, J.W., Erlandsson, R., Tommervik, H., Stordal, F., Tallaksen, L.M., 2023. Variation in albedo and other vegetation characteristics in non-forested northern ecosystems: the role of lichens and mosses. Environ. Res. Lett. 18 (7). https://doi.org/10.1088/1748-9326/ace06d. Article 074038.

Fortier, R., Banville, D.-R., Lévesque, R., Lemieux, J.-M., Molson, J., Therrien, R., Ouellet, M., 2020. Development of a three-dimensional geological model, based on quaternary chronology, geological mapping, and geophysical investigation, of a watershed in the discontinuous permafrost zone near Umiujaq (Nunavik, Canada). Hydrogeol. J. 28 (3), 813–832. https://doi.org/10.1007/s1004q-020-02113-1.

Frost, G.V., Epstein, H.E., 2014. Tall shrub and tree expansion in Siberian tundra ecotones since the 1960s. Glob. Chang. Biol. 20 (4), 1264–1277. https://doi.org/10.1111/gcb.12406.

Gagnon, M., Domine, F., Boudreau, S., 2019. The carbon sink due to shrub growth on Arctic tundra: a case study in a carbon-poor soil in eastern Canada. Environ. Sci. Communicat. 1 (9), 091001. https://doi.org/10.1088/2515-7620/ab3cdd.

Gamache, I., Payette, S., 2005. Latitudinal response of subarctic tree lines to recent climate change in eastern Canada. J. Biogeogr. 32 (5), 849–862. https://doi.org/10.1111/j.1365-2699.2004.01182.x.

Gauslaa, Y., 1984. 2. Infrared and visible reflectance in different lichen species and its ecological significance. Ecography 7 (1), 13–22. https://doi.org/10.1111/j.1600-0587.1984rb01100.x.

Jespersen, R.G., Anderson-Smith, M., Sullivan, P.F., Dial, R.J., Welker, J.M., 2023. NDVI changes in the Arctic: functional significance in the moist acidic tundra of northern Alaska. PLoS One 18 (4). https://doi.org/10.1371/journal.pone.0285030.

Ju, J.C., Masek, J.G., 2016. The vegetation greenness trend in Canada and US Alaska from 1984-2012 Landsat data. Remote Sens. Environ. 176, 1–16. https://doi.org/10.1016/j.rse.2016.01.001.

Juszak, I., Erb, A. M., Maximov, T. C., & Schaepman-Strub, G. (2014). Arctic shrub effects on NDVI, summer albedo and soil shading. Remote Sens. Environ., 153(0), 79–89. doi:https://doi.org/10.1016/j.rse.2014.07.021.

Kendall, M.G., 1975. Rank Correlation Methods, 4th edition. Charles Griffin.

Komsta, L., 2019. Median-based linear models. R package version 0.1.21. https://cran.r-project.org/web/packages/mblm/index.html.

Kuusinen, N., Juola, J., Karki, B., Stenroos, S., Rautiainen, M., 2020. A spectral analysis of common boreal ground lichen species. Remote Sens. Environ. 247, 111955. https://doi.org/10.1016/j.rse.2020.111955.

Lackner, G., Nadeau, D.F., Domine, F., Parent, A.-C., Leonardini, G., Boone, A., Anctil, F., Fortin, V., 2021. The effect of soil on the summertime surface energy budget of a humid subarctic tundra in northern Quebec, Canada. J. Hydrometeorol. 22 (10), 2547–2564. https://doi.org/10.1175/jhm-d-20-0243.1.

Lafleur, P.M., Humphreys, E.R., 2018. Tundra shrub effects on growing season energy and carbon dioxide exchange. Environ. Res. Lett. 13 (5). https://doi.org/10.1088/1748-9326/aab863. Article 055001.

Lantz, T.C., Marsh, P., Kokelj, S.V., 2013. Recent shrub proliferation in the Mackenzie Delta uplands and microclimatic implications. Ecosystems 16 (1), 47–59. https://doi.org/10.1007/s10021-012-9595-2.

Lemay, M.-A., Provencher-Nolet, L., Bernier, M., Lévesque, E., & Boudreau, S. (2018). Spatially explicit modeling and prediction of shrub cover increase near Umiujaq, Nunavik. Ecol. Monogr., 0 (0). doi:https://doi.org/10.1002/ecm.1296.

Lemieux, J.-M., Fortier, R., Murray, R., Dagenais, S., Cochand, M., Delottier, H., Therrien, R., Molson, J., Pryet, A., Parhizkar, M., 2020. Groundwater dynamics
within a watershed in the discontinuous permafrost zone near Umijujaq (Nunavik, Canada). Hydrogeol. J. 28 (3), 833–851. https://doi.org/10.1007/s10040-020-02110-4.

Malley, J.D., Kruppa, J., Dasgupta, A., Malley, K.G., Ziegler, A., 2012. Probability machines consistent probability estimation using nonparametric learning machines. Methods Inf. Med. 51 (1), 74–81. https://doi.org/10.3414 ме00-01-0052.

May, J.L., Oberbauer, S.F., Unger, S.L., Simon, M.J., Betway, K.R., Hollister, R.D., 2022. Shading decreases and delays NDVI and flowering of prostrate Arctic shrubs. Arctic Sci. 8 (3), 967–978. https://doi.org/10.1139/as-2020-0043.

Nagol, J.R., Sexton, J.O., Kim, D.-H., Anand, A., Morton, D., Vermote, E., Townshend, J. R., 2015. Bidirectional effects in Landsat reflectance estimates: is there a problem to solve? ISPRS J. Photogramm. Remote Sens. 103, 129–135. https://doi.org/10.1016/j.isprsjprs.2014.09.006.

Pattison, R.R., Jorgenson, J.C., Raynolds, M.K., Welker, J.M., 2015. Trends in NDVI and tundra community composition in the Arctic of NE Alaska between 1984 and 2009. Ecosystems 18 (4), 707–719. https://doi.org/10.1007/s10021-015-9858-9.

Payette, S., Delwaide, A., 2018. Tamm review: the north-American lichen woodland. For. Ecol. Manag. 417, 167–183. https://doi.org/10.1016/j.foreco.2018.02.043.

Payette, S., Morneau, C., 1993. Holocene relict woodlands at the eastern Canadian treeline. Quat. Res. 39 (1), 84–89. https://doi.org/10.1006/qres.1993.1010.

Payette, S., Fortin, M.-J., Gamache, I., 2001. The subarctic forest–tundra: the structure of a biome in a changing climate. BioScience 51 (9), 709–718. https://doi.org/10.1641/0006-3568(2001)051[0709:TSFTTS]2.0.CO;2.

Peltoniemi, J.I., Kaasalainen, S., Näränen, J., Rautiainen, M., Stenberg, P., Smolander, H., Smolander, S., Voipio, P., 2005. BRDF measurement of understory vegetation in pine forests: dwarf shrubs, lichen, and moss. Remote Sens. Environ. 94 (3), 343–354. https://doi.org/10.1016/j.rse.2004.10.009.

Provencher-Nolet, L., Bernier, M., Levesque, E., 2014. Quantification of recent changes to the forest–tundra ecotone through numerical analysis of aerial photographs. Ecoscience 21 (3–4), 419–433. https://doi.org/10.2980/21-(3-4)-3715.

Putkiranta, P., Räsänen, A., Korpelainen, P., Erlandsson, R., Kolari, T.H.M., Pang, Y., Villoslada, M., Wolff, F., Kumpula, T., Virtanen, T., 2024. The value of hyperspectral UAV imagery in characterizing tundra vegetation. Remote Sens. Environ. 308, 114175. https://doi.org/10.1016/j.rse.2024.114175.

Rautiainen, M., Kuusinen, N., Majasalmi, T., 2024. Remote sensing and spectroscopy of lichens. Ecol. Evol. 14 (3), e11110. https://doi.org/10.1002/ece3.11110.

Rees, W.G., Golubeva, E.I., Tutubalina, O.V., Zimin, M.V., Derkacheva, A.A., 2020. Relation between leaf area index and NDVI for subarctic deciduous vegetation. Int. J. Remote Sens. 41 (22), 8573–8589. https://doi.org/10.1080/01431161.2020.1782505.

Reinhardt, S., Aartsma, P., Skoyen, K., Renssen, H., 2022. Shrub encroachment interacts with environmental variation to reduce the albedo of alpine lichen heaths: an experimental study. Nord. J. Bot. 2022 (3). https://doi.org/10.1111/njb.03314.

Ropars, P., Boudreau, S., 2012. Shrub expansion at the forest-tundra ecotone: spatial heterogeneity linked to local topography. Environ. Res. Lett. 7 (1), 015501. Article 015501. https://doi.org/10.1088/1748-9326/7/1/015501.

Roy, D.P., Zhang, H.K., Ju, J., Gomez-Dans, J.L., Lewis, P.E., Schaaf, C.B., Sun, Q., Li, J., Huang, H., Kovalskyy, V., 2016. A general method to normalize Landsat reflectance data to nadir BRDF adjusted reflectance. Remote Sens. Environ. 176, 255–271. https://doi.org/10.1016/j.rse.2016.01.023.

Sales, M.H.R., Bruin, S., d., Souza, C., & Herold, M., 2022. Land use and land cover area estimates from class membership probability of a random Forest classification. IEEE Trans. Geosci. Remote Sens. 60, 1–11. https://doi.org/10.1109/TGRS.2021.3080083.

Schaaf, C.B., Gao, F., Strahler, A.H., Lucht, W., Li, X., Tsang, T., Strugnell, N.C., Zhang, X., Jin, Y., Muller, J.-P., Lewis, P., Barnsley, M., Hobson, P., Disney, M., Roberts, G., Dunderdale, M., Doll, C., d'Entremont, R.P., Hu, B., Roy, D., 2002. First operational BRDF, albedo nadir reflectance products from MODIS. Remote Sens. Environ. 83 (1), 135–148. https://doi.org/10.1016/S0034-4257(02)00091-3.

Sen, P.K., 1968. Estimates of the regression coefficient based on Kendall's tau. J. Am. Stat. Assoc. 63 (324), 1379–1389. https://doi.org/10.1080/01621459.1968.10480934.

Tape, K., Sturm, M., Racine, C., 2006. The evidence for shrub expansion in northern Alaska and the Pan-Arctic. Glob. Chang. Biol. 12 (4), 686–702. https://doi.org/10.1111/j.1365-2486.2006.01128.x.

Tremblay, B., Levesque, E., Boudreau, S., 2012. Recent expansion of erect shrubs in the low Arctic: evidence from eastern Nunavik. Environ. Res. Lett. 7 (3), 035501. Article 035501. https://doi.org/10.1088/1748-9326/7/3/035501.

Walther, C., Huttich, C., Urban, M., Schmullius, C., 2019. Modelling the Arctic tundra ecotone using ALOS PALSAR and optical earth observation data. Int. J. Appl. Earth Obs. Geinf. 81, 195–206. https://doi.org/10.1016/j.jag.2019.05.008.

Williamson, S.N., Barrio, I.C., Hik, D.S., Gamon, J.A., 2016. Phenology and species determine growing-season albedo increase at the altitudinal limit of shrub growth in the sub-Arctic. Glob. Chang. Biol. 22 (11), 3621–3631. https://doi.org/10.1111/gcb.13297.

Wolf, K., Jäkel, E., Ehrlich, A., Schäfer, M., Feilhauer, H., Huth, A., Weigelt, A., Wendisch, M., 2024. Impact of clouds on vegetation albedo quantified by coupling an atmosphere and a vegetation radiative transfer model. EGUSphere 2024, 1–30. https://doi.org/10.5194/egusphere-2024-3614.

Wong, R.E., Berner, L.T., Sullivan, P.F., Potter, C.S., Dial, R.J., 2024. Pixel walking along the boreal forest–Arctic tundra ecotone: large scale ground-truthing of satellite-derived greenness (NDVI). Glob. Chang. Biol. 30 (6), e17374. https://doi.org/10.1111/gcb.17374.

Yang, J.X., Li, Z.H., Zhai, P.M., Zhao, Y.Y., Gao, X.Q., 2020. The influence of soil moisture and solar altitude on surface spectral albedo in arid area. Environ. Res. Lett. 15 (3). https://doi.org/10.1088/1748-9326/ab6ae2. Article 035010.

Yang, D., Morrison, B.D., Hanston, W., McMahon, A., Baskaran, L., Hayes, D.J., Miller, C. E., Serbin, S.P., 2023. Integrating very-high-resolution UAS data and airborne imaging spectroscopy to map the fractional composition of Arctic plant functional types in Western Alaska. Remote Sens. Environ. 286, 113430. https://doi.org/10.1016/j.rse.2022.113430.

Zhang, H.K., Roy, D.P., Kovalskyy, V., 2016. Optimal solar geometry definition for global long-term Landsat time-series bidirectional reflectance normalization. IEEE Trans. Geosci. Remote Sens. 54 (3), 1410–1418. https://doi.org/10.1109/TGRS.2015.2480684.

Zhu, Z., Woodcock, C.E., 2012. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sens. Environ. 118, 83–94. https://doi.org/10.1016/j.rse.2011.10.028.

Zhu, Z., Wang, S., Woodcock, C.E., 2015. Improvement and expansion of the Fmask algorithm: cloud, cloud shadow, and snow detection for Landsats 4–7, 8, and sentinel 2 images. Remote Sens. Environ. 159, 269–277. https://doi.org/10.1016/j.rse.2014.12.014.