Introduction

Tundra vegetation is experiencing distributional and compositional changes as a response to climate warming, with the expansion of woody plants to treeless tundra (including growth and infill of existing woody plants, i.e., shrubification) being among the most fundamental phenomena1,2,3,4. However, this vegetation redistribution and increase in biomass can be constrained by cryogenic land surface processes (LSP), which are presently active over ca. 25% of the Earth’s terrestrial areas that are exposed to seasonal and/or permanent ground freezing5,6,7,8,9,10. In arctic and alpine tundra, cold climatic conditions prevail due to high latitude and/or altitude, and plants and other lifeforms are often near their thermal tolerance limits11,12. Here, LSP are expected to mainly impact vegetation by causing mechanical disturbances. Repeated ground freezing and thawing can destabilize the topsoil (here the maximum rooting depth of vegetation is found to be <30 cm for most tundra plants13) thus preventing the formation of continuous vegetation cover, further constraining forest-tundra ecotone dynamics12. In addition, LSP mediate other local conditions; for example, soil moisture and growing season length are strongly constrained by snowpack dynamics in snowbed environments (i.e., nivations), thus affecting species distributions, community compositions and vegetation productivity in tundra ecosystems14. Moreover, cryogenic mixing of soil horizons may significantly contribute to carbon and nutrient cycling15,16.

Recent studies suggest a widespread loss of active LSP (i.e., stabilization of soil surface and/or disappearance of snow beds) due to climate change17,18. If realized, this could result in feedbacks where the warmer climate generates an increase in vegetation biomass and uptake of atmospheric CO2, which are then further accelerated by the stabilization of ground surfaces. In addition, changes in vegetation may also impact LSP activity19. Not only can changes in vegetation profoundly impact tundra ecosystems’ functioning by altering productivity and biotic interactions, they can also feed back into the global climate system through changes in albedo, evapotranspiration, permafrost and carbon cycling7,20,21,22,23,24. Yet, the effect of multiple LSP in constraining future tundra biomass over broad environmental gradients has remained unexplored and consequently unaccounted for in simulation models concerning tundra biomass.

Here, we use comprehensive field-quantified and remotely sensed data from northern Europe to incorporate the impacts of multiple LSP on spatial predictions of tundra above-ground biomass (AGB; Fig. 1). We first establish a statistical relationship with field-measured AGB (433 field sites) and satellite-derived normalized difference vegetation index (NDVI; Landsat 8) to produce spatially extensive estimates of AGB over 2917 observations sites. Then using the larger data we focus on the effects of three dominant LSP on the estimated AGB that create typical surface-geomorphological features of arctic and alpine tundra landscapes5; patterned ground, frost boils and hummocky terrain created by cryoturbation (i.e., frost-churning), lobes and terraces occurring at the mountain slopes produced by solifluction (i.e., slow mass wasting mostly due to cyclic freezing and thawing of soil), and nivation processes comprising weathering and fluvial processes indicated by nivation hollows and snowbed environments. Our study area represents broad regions of the forest-tundra ecotone25. It ranges from low-lying sub-Arctic coniferous and birch forests to barren arctic-alpine tundra, characterized by intense cryogenic soil activity. The study area has strong climatic gradients (mean annual air temperature range ca. −6.5 to 1.5 °C, mean annual precipitation sum range ca. 400–900 mm) and is mainly representative of sporadic and discontinuous permafrost5,26 (Fig. 1). Within this climatic domain representing arctic-alpine tundra, the presence of the three LSP and their role in ecosystem functioning is widely documented5,6,11. Moreover, the discontinuity of permafrost, relatively high ground temperatures and low ground-ice content (compared to colder tundra subzones) can indicate that the soil thermal response t. o climatic changes is likely to be rapid27. Accordingly, a recent study by Aalto et al., (2017; ref. 17) suggested that the study area is likely to experience substantial near-term reductions in the distribution of multiple cryogenic LSP. This, in turn, has the potential to alter landscape patterns and dynamics, and affect features such as slope movements, seasonal thaw depths and surface hydrology18,28.

Fig. 1: The study area and observation sites, and their climatic coverage.
figure 1

a Location of the study area (black box) in respect to the circum-Arctic region. b Normalized difference vegetation index (NDVI) overlain by observation sites (n = 2917) and type of LSP (Nivat = nivation, Solif = solifluction, Cryot = cryoturbation). The black rectangle indicates the model prediction domain. c–e show field photos of typical surface features of LSP: cryoturbation (frost-churning), solifluction (solifluction lobes) and nivation (snow accumulation sites), respectively. f the relationship between AGB measured in-situ (sampling regions are labeled as “BM region 1” and “BM region 2”, total number of measurements is 433) and NDVI; the red line depicts a fitted generalized least squares function used to estimate AGB for the LSP sites (Methods). The kernel smoothed NDVI distributions in the presence of any of the three LSP (n = 1440) and absence of all LSP (n = 1477) are shown in the background. g Permafrost (PF) conditions at each LSP observation site according to a recent ground thermal modeling26. The PF categories continuous, discontinuous, sporadic, and PF absent are associated with the presence of at least one LSP, whereas the last category (PF absent, LSP absent) does not include any LSP.

Full size image

Here, we first quantify the LSP effect on spatial patterns of above-ground biomass in arctic alpine tundra, and then investigate the potential implications of these relationships under warming climate. We deploy a machine learning framework (boosted regression trees, BRT) to statistically relate the estimated AGB to the current activity of LSP at a fine spatial scale29,30 while controlling for the effects of key climatic and local environmental gradients, and to account for spatial correlation structures in the data (Methods). The use of machine learning allows us to model complex, non-linear relationships, and interactions among predictors. BRT also provides an estimate of each predictors’ relative influence in a multivariate modeling setting. We quantify the effects of LSP on AGB by predicting AGB over the study area at high spatial resolution using model configurations that both neglect and account for the effects of the three LSP on AGB. These analyses allow us to compare the locations where LSP are either present or absent, after controlling for other environmental conditions. We expect that in the studied arctic-alpine system (mainly sporadic-discontinuous permafrost, low ground ice content) LSP will demonstrate secondary importance to AGB after the climatic factors.

We show that three surveyed LSP collectively reduced AGB in the study area on average by nearly one-third (30%, 123.0 g m−2). This effect was significant over various landscape positions and was especially pronounced in snowbed environments (mean AGB reduction 57.3%). Our results suggest that LSP are important in shaping future patterns of tundra biomass, as long as cryogenic ground activity is retained by climate warming.

Results

LSP effects on above-ground biomass

Sites with active LSP present were associated with significantly lower NDVI and consequently estimated AGB compared to sites with no LSP activity (Supplementary Table 1). We found that the current AGB patterns were most strongly related to growing degree days (GDD, mean effect size = 653.8 g m−2 [509.8–797.7 g m−2, 95% range of variation]), potential incoming solar radiation (RAD, 430.7 g m−2 [115.8–842.0 g m−2]) and soil chemistry (pH, 306.8 g m−2 [290.8–323.3 g m−2]; Supplementary Fig. 1). Within the LSP realm, the total effect of the three LSP on AGB was −123.0 g m−2 (−30.0%; Fig. 2a) on average (−77.4 to −168.7 g m−2). There were significant negative effects of LSP on AGB (P ≤ 0.001, paired one-sided t-test, n = 500) across the AGB gradient in our study area; the presence of active cryoturbation decreased AGB by an average of −86.6 g m−2 (−26.3, −64.1 to −109.2 g m−2) compared to sites without cryoturbation, but otherwise similar climatic, topographic and soil conditions. On mountain slopes, which are important transition zones between forest and tundra biomes12, the presence of active solifluction reduced AGB by an average of −21.5 g m−2 (-21.7%, -7.1 – -35.9 g m-2). In high-elevation snowbed environments the presence of nivation led to a minor absolute (−14.9 g m−2, −6.2 to −23.6 g m−2), but a large relative average decrease (−57.3%) in AGB. These results support our argument that cryogenic soil disturbance is likely to impact the accumulation of vegetation biomass in cold regions over a wide range of landscape positions.

Fig. 2: The effect of cryogenic land surface processes (LSP) in the biomass models.
figure 2

a The bar plots show the effect sizes of the observed LSP on estimated AGB after controlling for the effects of other environmental predictors. The lines indicate 95% confidence intervals based on boosted regression trees (BRT) modeling over 500 bootstrap samples. Asterisks indicate statistically significant differences between sites with LSP absent and LSP present (P ≤ 0.001, paired one-sided t-test, n = 500). b Prediction error of the AGB models in terms of root mean squared error (RMSE) between observations and each of the model predictions (calibrated first without LSP [Base model] and then with LSP [LSP model]). The boxplots depict median, 1st and 3rd quartiles, and 95% range of variation over 500 cross-validation rounds. Asterisks indicate statistically significant model improvement (P ≤ 0.01, paired one-sided t-test, n = 500). c Predictors’ relative influence in the above-ground biomass (AGB) models. Here, the LSP includes the summed influences of cryoturbation, solifluction and nivation (numerical results are presented in Supplementary Table 2). The black lines indicate the standard deviation over 500 BRT model runs with bootstrap sampling. GDD growing degree days, WAB water balance, Tfeb February mean air temperature, RAD potential incoming solar radiation, RAC residual autocovariate for the LSP model (i.e., RACLSP), TWI topographic wetness index, Peat peat cover.

Full size image

We tested the effect of LSP on AGB by using two separate BRT model configurations. The first model did not include the LSP (hereafter “Base model”); the second model included the LSP (“LSP model”). The LSP model showed a significant improvement in the predictive performance (P ≤ 0.01, paired one-sided t-test, n = 500) in terms of root mean squared error (RMSE; mean RMSEBase model = 332.1 g m−2 (186.8–572.4 g m−2, 95% range of variation), mean RMSELSP model = 316.8 g m−2 (175.7–547.8 g m−2); Fig. 2b). This significant decrease in RMSE indicates that LSP should be considered in biomass models in order to better predict local AGB patterns in arctic and alpine tundra. The relative influence of individual LSP remained low on average (Supplementary Table 2), but their combined influence of ca. 6.4% was ranked as the fourth highest after GDD, pH and WAB (Fig. 2c). This finding complies with our expectation that LSP can be considered as secondary factors controlling the AGB patterns in these systems.

Our data and analyses are associated with a considerable amount of uncertainty as indicated by large variability around mean estimates and prediction error. Such uncertainties are likely to reflect landscape-level heterogeneity in considered environmental variables not captured by our data, and uncertainty related to the NDVI-model used to derive AGB estimates for the LSP sites. In addition, the deployed visual approach for mapping active LSP is likely to contain uncertainty. Whereas active cryoturbation can be readily identified based on indications of frost-heaving and change in vegetation cover, active solifluction may be more challenging to identify without direct velocity measurements. This is because the maximum depth of soil movement can be deeper than the rooting depth enabling vegetation to grow on active solifluction features19. This may partly explain our finding of a relatively smaller reduction in AGB over solifluction sites compared to cryoturbation or nivation sites (Fig. 2a). In addition, nivation processes are commonly active in late-lying snowbeds which have relatively consistent year-to-year patterns31, making such snowbeds good indicators of the presence of active nivation processes28. However, it is likely that some of the found effects of nivation on AGB derives from other mechanisms such as snow cover limiting energy input. Moreover, our statistical modeling framework assumes that above-ground biomass in a given area is at equilibrium with its environment and that the relevant environmental variables have been adequately sampled. Our input data showed clear spatial autocorrelation structures for all predictors. For climate predictors the spatial autocorrelation was found to extend longer horizontal distances than for local topography and soil-related predictors, and LSP (Supplementary Figs. 2–3). Such data structures were explicitly accounted for in the modeling, and the consequent analyses were not confounded by spatial autocorrelation (Supplementary Figs. 4–5; see Methods).

Future AGB trajectories under climate warming

In tandem with the projected increase of AGB20, the studied active LSP are expected to widely decrease as a response to climate change17 (Supplementary Figs. 6–8), with largely unknown consequences for tundra vegetation. To examine this, we simulated potential future AGB trajectories over the study area using the two model configurations (i.e., Base model and LSP model). We predicted AGB at high spatial resolution and under several climate warming scenarios constrained by the corresponding predicted changes of LSP (Methods). The results suggest that LSP can significantly hinder the increase in AGB, leading to future trajectories that are suppressed compared to the predictions that do not account for LSP (Fig. 3). Here, we stress that our data and methods do not enable us to account for two-way soil-vegetation interactions and feedbacks or factors such as herbivory, and thus, these future projections should be interpreted with caution because of the inadequacies of the data on which they are based.

Fig. 3: Cryogenic land surface processes suppress tundra biomass trajectories under climate warming.
figure 3

The above-ground biomass (AGB) predictions over the study area are based on two boosted regression trees (BRT) model configurations: Base model (climate + topography and soil) and LSP model (Base model + cryoturbation, solifluction, and nivation). In a–c, cryoturbation, solifluction, and nivation are independently considered, respectively; in d, all LSP are considered together. The AGB were extracted from areas suitable for LSP occurrences at each temperature scenario (\(\triangle\)T), where monthly air temperatures (averages of 1981–2010) were modified from +0.5 °C to +8.0 °C at 0.5 °C interval. The lines depict the means and the shaded areas represent 50% of the predictions. The crosses indicate statistically non-significant (P > 0.05) differences between the predictions (a paired one-sided Wilcoxon’s test).

Full size image

Even after a considerable warming and loss of LSP, the models accounting for the LSP effects project lower AGB values than the Base models, meaning that LSP potentially still restricts the accumulation of phytomass under future climates. These LSP effects are particularly evident over slope and snowpack environments. For example, under the +3 °C warming scenario the LSP model predicts on average 71.7 g m−2 (−21.1%) lower AGB (P ≤ 0.001, paired one-sided Wilcoxon’s test), compared to the Base model at sites maintaining solifluction activity in the future. Similarly, the corresponding mean reduction for nivation sites due to disappearance of active nivation would be 75.0 g m−2 (−45.8%, P ≤ 0.001) compared to prediction by the Base model.

However, it is reasonably to assume that the interactions between cryogenic soil processes and vegetation may differ across tundra vegetation types6,12,32. In snowbed environments the large effect of nivation is partly related to the fact that in these environments the AGB is generally low. Observational and experimental studies across the Arctic suggest that the warming may have contrasting effect on community compositions depending on the starting vegetation type and the responsiveness to climate warming can vary greatly9,33. This is likely to be true for communities characterized by active LSP. Here, we cannot separate AGB responses to LSP across vegetation types, but our results suggest that LSP are one additional factor that could cause the contrasting vegetation responses to warming within and between landscapes and regions. Notably, cryoturbated soils can act as favorable new microsites for seed germination and seedling establishment e.g., for tall shrubs as observed in other tundra regions32,34 as repeated disturbances limit biotic competition and heave nutrients from deeper soil layers. After LSP activity has stopped or at least markedly reduced, the biomass accumulation can thus be rapid. Then also biogeomorphic feedbacks can come to play; once LSP activity has decreased to a certain level, an increasing vegetation cover can accelerate the stabilization of the cryogenic processes19. Thus, LSP may be one important mechanism of how climate is controlling local ecotone (e.g., treeline) dynamics by facilitating vegetation establishment and altering competition35,36,37. This further underlines the multifaceted role of LSP in shaping tundra AGB, and the potentially complex and nonlinear AGB-LSP responses in the future8,9,14,36,37,38.

We further estimate that cryogenic soil disturbance currently prevents 5.5% (0.12 Tg) of aboveground carbon storage from being sequestered in the study area, assuming 50% of the AGB is carbon39. In tundra landscapes, differences in AGB translate to pronounced spatial heterogeneity in short-term carbon stocks, especially at poleward aspects where the microclimate can sustain active LSP long after macroclimatic conditions become unsuitable14 (Fig. 4). Here we suppose that the main impact of LSP on carbon stocks derive from limiting vegetation establishment and growth. However, as the majority of carbon in the tundra ecosystems is stored in seasonally or permanently frozen soils40, further research is needed to quantify the role of LSP in controlling below-ground carbon stocks. For example, cryoturbation is known to be an important mechanism to transport organic material from the soil surface into deeper layers, thus increasing long-term carbon storage40. Our results indicate that an increase in AGB—possibly regarding tall and long-lived woody vegetation22—could be delayed due to cryogenic soil disturbance (Figs. 3 and 4). Therefore, LSP can indirectly mediate biophysical feedbacks such as changes in ground reflectance and latent heat fluxes by impacting tundra ecosystem structure (e.g., vegetation height and volume, and canopy closure) and functioning (nutrient cycling and soil moisture variations)20,23. Plant species and functional types have different capabilities to cope with (or even take advantage of) physical disturbance41. Thus, loss of cryogenic soil activity should lead to a plant succession that may take time to reach a new equilibrium in an overall slow-reacting tundra vegetation.

Fig. 4: Above-ground biomass (AGB) patterns in tundra are constrained by cryogenic land surface processes (LSP).
figure 4

The close-ups (a–c Base model, d–f LSP model) show the predicted AGB in a topographically heterogeneous landscape (ca. 69°03’N, 20°51’E) at present climate (1981–2010), and under two mean annual air temperature warming scenarios (i.e., +1.5 °C and +3.0 °C), respectively. In a–c, the contours depict elevation from sea level; in d–f, the gray hatched polygons encompass the areas suitable for the occurrence of LSP (i.e., the LSP realm) at each warming scenario. The boxplots indicate the AGB variation (black line depicting median, the whiskers extending to 1st and 99th percentile) within the LSP realm. Asterisks indicate statistically significant differences between the two model steps (P ≤ 0.001, paired one-sided Wilcoxon’s test, n = 1000 randomly chosen pixels).

Full size image

Our result suggesting nearly one-third reduction in AGB due to LSP at the study area underline the fact that cryogenic land surface dynamics must be carefully considered when assessing current and future biomass patterns over cold-climate regions. We further postulate that recent cryospheric alterations, such as decrease of LSP and reduced seasonal snow cover7,42 could partly explain the observed greening trends of the tundra regions2,9. Our simulations suggest that the decrease in soil disturbance could trigger non-linear feedbacks where AGB increase (e.g., shrub expansion to tundra) is first delayed (i.e., diverging AGB trajectories between the models in Fig. 3), but later accelerated following ground surface stabilization after disappearance of active LSP. Therefore, cryogenic disturbances can cause a short-term reduction in AGB, but over a longer time period this effect is likely to diminish. Nevertheless, the results suggest that the LSP still affect AGB patterns at the end of 21st century. Notably, as opposed to our study area, in poorly drained tundra regions associated with continuous and ice-rich permafrost (e.g., Northern Siberia, Alaska, and Canada) climate warming can increase LSP activity e.g., due to promoting the development of permafrost thaw slumps43 and ice-wedge degradation44. Therefore, large uncertainties are related to the timing and spatial characteristics of AGB-LSP responses due to associated geo-ecological lags related to ground thermal response, altered land surface dynamics, and rates of vegetation establishment.

These poorly explored responses of tundra vegetation to climate warming can have far-reaching consequences for global atmospheric conditions via albedo, evapotranspiration feedbacks, and carbon balance20,45,46,47,48. Thus, our results challenge the projections of cold-region vegetation dynamics in the near and intermediate future (50–100 years)20,45,49, and provide a solid foundation for investigating vegetation-LSP relationships over other parts of tundra with different environmental conditions. We show that further development of cryogenic ground components is urgently required for the next generation ecosystem and land surface models over large spatial domains.

Methods

Study area

The study area (78 000 km2) is located between 68–71°N and 20–26°E, with strong climatic gradients, ranging from wet maritime to relatively dry continental, over tens of kilometers. The landscape of this climatically sensitive high-latitude region has been affected by multiple glaciations in the past. It includes the Scandes Mountains near the Arctic Ocean and low-relief areas to the south and east. The majority of the region (52%) is underlain by sporadic permafrost. Continuous and discontinuous permafrost are limited to the highest mountains of the study area (2% and 7%, respectively)17,26. This large proportion of sporadic, typically warm and shallow permafrost in the study area indicates that ground thermal response to climate warming can be rapid27. Our data do not cover low-relief plateaus of continuous permafrost (similar to northern Siberia and Alaska), where the generally high ice content of soil may lead to different and enhanced LSP responses under climate warming (e.g., ice wedge degradation and surface ponding) with altered AGB feedbacks43,44.

LSP observations

The data consist of 2917 study sites (each 25 m × 25 m) and includes previously combined observations (both in-situ [n = 581] and remote-sensing [n = 2336]) of the active surface features of three cryogenic LSP common in the area: cryoturbation, solifluction, and nivation. These LSP are mainly associated with seasonal freeze–thaw processes. Cryoturbation (i.e., frost churning) is a general term for soil movement caused by differential heave, and it creates typical surface features such as patterned ground, frost boils and hummocky terrain5. Solifluction is the slow mass wasting of surficial deposits through frost creep and permafrost flow, where gravitation causes frost-heaved soil to settle downwards during the summer thaw, creating features of lobes and terraces50. In addition, solifluction also includes gelifluction which is a mass wasting process caused by high porewater pressure in unconsolidated surface debris creating similar lobes and terraces5,50. We use the term nivation to collectively designate various weathering and fluvial processes which are intensified and depicted by the presence of snowbeds (which in general are melting in mid-July – late-August) and nivation hollows28,51. We expect the presence of such a snowbed to be an indication of active nivation processes, since in these environments the year-to-year spatial snow patters are fairly consistent31.

The rationale behind LSP sampling is described in previous geomorphic studies which served as a basis for the used protocol52,53. Due to the large study domain, study objectives (focus on distribution of active surface features, not on activity itself) and modeling resolution (50 m × 50 m), we used a visual method to estimate the presence/absence of the mapped LSP. We used high-resolution aerial photography (spatial resolution of 0.25 m−2) and targeted field surveys (GPS accuracy ~5 m; Garmin eTrex personal navigator) to construct the LSP dataset. A binary variable (1 = presence, 0 = absence) was assigned to each LSP to indicate their evident activity (or absence). The activity/absence of the LSP was visually estimated based on the evidence in ground surface, indicated by e.g., frost-heaving, cracking, microtopography (e.g., erosional and depositional forms), soil displacement indicative to a process form (e.g., solifluction lobes, patterned ground), changes in vegetation cover and late-lying snow. Such indicators average the LSP activity over several years. Even small areas with slight indication of activity were considered active processes. However, such a protocol based on a visual assessment is susceptible for incorrect activity classification; solifluction may be active despite having a complete vegetation cover19 and the presence of late-lying snowbed, although being a good indication28, does not necessarily mean that active nivation processes are present.

Remotely sensed vegetation index

For obtaining remotely sensed vegetation index for the study area, we employed a maximum-value compositing approach. We downloaded all available clear sky (less than 80% land cloud cover) Landsat OLI 8 images overlapping the study area from June to September between 2013 and 2017 (total of 1086 scenes) from the United States Geological Survey (USGS) database (http:\\earthexplorer.usgs.gov\). Images were USGS surface reflectance products, which were preprocessed (georeferencing, projection, and atmospheric corrections) by USGS54. Landsat-8 satellite is the latest addition to the Landsat mission that has provided repeated land surface information globally since the 1970’s and is the most commonly used fine-scale satellite system for vegetation mapping. The native resolution of the Landsat OLI sensor is 30 m for the spectral bands used in the image processing steps of this study.

Normalized difference vegetation index (NDVI), a widely used spectral index to estimate the amount of green vegetation, was calculated as55:

$$({{{{{\rm{\rho }}}}}}{{{{{\rm{NIR}}}}}}-{{{{{\rm{\rho }}}}}}{{{{{\rm{red}}}}}})/({{{{{\rm{\rho }}}}}}{{{{{\rm{NIR}}}}}}+{{{{{\rm{\rho }}}}}}{{{{{\rm{red}}}}}})$$
(1)

where ρNIR and ρred are the surface reflectance for their respective Landsat bands, 0.851–0.879 \(\mu\)m and 0.636–0.673 \(\mu\)m.

USGS provides pixel-based quality assessment bands for all surface reflectance products. These bands were used to mask clouds, snow, water, and other low-quality pixels from the individual NDVI scenes. Additionally, if the NDVI images still had unphysical values over 1 or under -1, these pixels and their surroundings of 100 m radius were excluded. We determined maximum values for each 30 m resolution pixel of the study area individually. After masking cloud, snow, and water from the scenes, obvious scattered erroneous NDVI values remained in some scenes. Therefore, we excluded the values outside the pixel-based 95% percentile prior to maximum composite.

The CFmask cloud detection algorithm that is used to generate the quality assessment band has clear difficulties in distinguishing small snow patches from clouds. As such, a large portion of late-lying snow beds were repeatedly and incorrectly classified as clouds. Moreover, the CFmask algorithm creates buffers around the cloud pixels54, hence much information was lost around the snow patches that were incorrectly identified as clouds. After these processing steps, some pixels around the extreme late-lying snow beds had still too low number of NDVI records to provide reliable NDVI values for the maximum composite. To fill these small and scattered gaps in the initial maximum NDVI composite, we selected 74 mostly cloud-free scenes between August and September. For these 74 scenes, we manually digitized cloud masks to exclude cloud-contaminated pixels with high certainty. Moreover, every pixel must have passed the following quality checks to be included in the gap-filling composite: not classified as water in the USGS quality assessment band; normalized difference snow index (NDSI) value less than 0.4, and blue band reflectance less than 0.1 (to exclude snow); reflectance of red band between 0.03 and 0.4 (second check for water and snow, and deepest shadows); NDVI between 0 and 0.4 (lower threshold to exclude snow and water contamination; higher threshold to exclude erroneous values, as very late snowbed habitats always have very limited vegetation cover). Additionally, if the NDVI images had unphysical values over 1 or under -1, these pixels and their surroundings (200 m radius) were excluded. Pixels in the 74 selected images which passed these checks, were then used to create a secondary maximum NDVI composite that was used to fill the gaps in the initial maximum NDVI composite. The secondary composite comprised 0.4% of the pixels in the final composite. Among all 2917 LSP observation sites, 2.9% were located within the gaps in the initial maximum NDVI composite, and thus received their maximum NDVI values from the secondary NDVI composite.

In the used Landsat data, the nivation sites were not covered by snow, but instead were associated with generally lower AGB values as nival processes affect the vegetation’s structure and composition (Supplementary Table 1).

Above-ground biomass data

Above-ground biomass (AGB) reference data were collected from two regions, with a total of 433 sites that represent an area of > 4000 km2 (Supplementary Fig. 9). The first dataset (hereafter BM region 1; centering to ca. 69°N, 21°E) was collected between 2008 and 2011, and the second dataset (BM region 2; centering to ca. 70°N, 26.2°E) between 2015 and 2017. Both study regions are representative of an arctic and alpine treeline ecotone and include data from mountain birch forest to barren oroarctic tundra56,57.

The BM region 1 dataset consists of 309 field sites (each 10 m × 10 m), which are located around eight different massifs covering a wide range of environmental conditions (Supplementary Figs. 9–10). Sampling was performed in transects to cover various aspects of the slope (i.e., topoclimatic conditions), starting from the foothill of the mountain, and ending at the summit. A plot was systematically established at every 20 m increase in elevation and recorded with a GPS device. Four clip-harvest biomass samples (20 cm × 20 cm) were taken 5 m from the plot center in every cardinal direction. Two samples were used in bare mountaintops (north, south). The clip-harvest samples were dried for 48 h at +65 °C, and dry weight was recorded. The sample biomass values were converted to g m-2 and the average sample value was calculated for each site (Supplementary Fig. 9). The original BM region 1 dataset contains forest and treeline plots, but these were excluded from the final analyses due to an incomparable tree sampling strategy with BM region 2, which could introduce uncertainty into biomass estimates.

The BM region 2 data were collected from three different massifs having an elevation range from 120 m to 1064 m (Supplementary Fig. 9). The biomasses were sampled from 102 sites (each 24 m × 24 m in size) that were chosen using a stratified sampling to cover gradients of thermal radiation (potential incoming solar radiation), soil moisture (topographic wetness index, TWI) and vegetation zone (forest, treeline, and alpine zones). Radiation and TWI were calculated from a 10 m digital elevation model (DEM, provided by the National Land Surveys of Finland and Kartverket, the Norwegian mapping authority), and assigned to one of three classes based on observation percentiles (breaks at 20% and 80%) leading to total of 27 strata. Vegetation zones were digitized based on aerial imagery. After the first field survey, 22 sites were added to account for vegetation types that were not sufficiently represented by the GIS-based stratification. Thus, the total sample size of the BM region 2 dataset is 124 AGB sites.

The same clip-harvest sample protocol was used as in BM region 1; additional samples were also taken from 12 m in every cardinal direction, thus each site had eight AGB samples (Supplementary Fig. 9). Trees with diameter at breast height (DBH) greater than 20 mm were measured from a 900 m2 circular plot, which corresponds to the size of the NDVI product resolution. Large stems (DBH > 80 mm in the forest and 40 mm at the treeline) were measured from the whole plot, whereas smaller stems were measured from five subplots. Specifically, the center subplot was 100 m2, and the four subplots located at 8 m to every cardinal direction were each 12.5 m2. For the subplot observations, we used a plot expansion factor (900/150 = 6) to generalize the observations for the whole plot assuming a homogeneous forest structure i.e., each subplot stem represents six trees within the 900 m2 plot. A total of 98% of the measured stems were mountain birch (Betula pubescens ssp. czerepanovii), making it the most abundant species in the area. For predicting stem biomass, we used the average of three allometric equations58,59,60, in order to reduce the uncertainty related to the transferability of an individual allometric model. In addition, Populus tremula (1% of the observations) were found on low-altitude south-facing slopes, and Salix caprea (1%) in moist, nutrient-rich sites. Species-specific models61,62 were used to estimate their respective stem biomasses. Individual pines (Pinus sylvestris) were scattered in the area but were not present in any of the sampled plots.

The plots of above-ground tree biomass were converted to g m−2 and added to the mean clip-harvest AGB to obtain the total vascular plant AGB for each site. The BM region 1 and BM region 2 datasets were combined, and the NDVI value was extracted from the site center coordinates.

Spatial autocorrelation (SAC) is a common property of any spatial dataset and means that observations are related to one another by the geographical distance63. SAC in the model residuals violates the independence assumption commonly required by statistical models and can lead to inflated hypothesis testing and biased model estimates64. To investigate whether the plot-scale AGB data are spatially autocorrelated, we calculated semivariogram which describes the spatial dependency between the observations as a function of distance between the point pairs65. Semivariogram were calculated as:

$${{{{{\rm{\gamma }}}}}}(h)=\frac{1}{2N(h)}\mathop{\sum }\limits_{i=1}^{{N}_{h}}{\left(Z\left({s}_{i}\right)-Z\left({s}_{i}+h\right)\right)}^{2}$$
(2)

where N(h) denotes the number of data pairs within distance h, and \(Z\left({s}_{i}\right)\) is an observation (or model residual) in location i. For the calculation, we used R package gstat66 (version 2.0-0). A visual inspection of the semivariogram indicated spatial autocorrelation at short distances (“AGB” in Supplementary Fig. 11). Therefore, for the NDVI-AGB conversion, we used a generalized least squares modeling (GLS, as implemented in R package nlme67 [version 3.1-137]) that can explicitly account for SAC in the data. For the modeling, the AGB values were log(x+0.1) transformed. The GLS, where AGB was modeled as a function of NDVI, were fitted assuming an exponential spatial correlation structure:

$$\gamma \left(h\right)={c}_{0}+c\left(1-{e}^{-h/a}\right)$$
(3)

where \({c}_{0}\) is the difference between the intercept and origin (i.e., the “nugget” parameter in geostatistics), \(c\) is the amount of variance (i.e., the “sill”) and a represents the distance of spatial dependency (i.e., the “range”). The fitted GLS was as follows:

$${{\log }}\left({AGB}\right)=-1.038629+9.725572\times {NDVI}$$
(4)

The estimated spatial correlation parameters were c0 = 0.516, c = 0.484 and a = 260.605, indicating that the distance of spatial autocorrelation extends to ca. 261 m. The semivariogram for the model residuals indicated a notable reduction in the amount of spatial autocorrelation compared to the AGB data (Supplementary Fig. 11). The fitted model explained 70.6% of the deviance in the data. When the predicted values were converted back to the response scale, the model explained 60.5% of the deviance. Therefore, for the subsequent analyses we use the above-ground biomass estimated by the model.

Environmental predictors

In addition to LSP, we used climate, topography, and soil predictors to model AGB. Gridded monthly average temperatures and precipitation data (1981–2010; spatial resolution 50 m × 50 m) based on a large collection (> 950) of Fennoscandian meteorological stations were used in a spatial interpolation scheme17. Three climate predictors—growing degree days (GDD, °C, base temperature 5 °C), mean February air temperature (Tfeb, °C) and water balance (WAB, mm)—were calculated from the gridded climate data. WAB is the difference between total annual precipitation and potential evapotranspiration (PET), which was estimated from the monthly air temperature and precipitation data68:

$${PET}=58.93\times {T}_{{above}0^\circ C}/12$$
(5)

These climatic predictors were selected to represent different aspects of climate that are critical for tundra vegetation: heat requirements, cold tolerance and moisture availability. In addition, two local scale topographic predictors were calculated from a DEM (spatial resolution of 50 m × 50 m, provided by the National Land Survey Institutes of Finland, Norway, and Sweden): topographic wetness index69 (TWI, a proxy for soil moisture) and potential annual direct solar radiation70 (MJ cm-2 a-1). Slope angle was initially considered as a potential predictor for AGB but was later omitted due to the strong correlation with TWI (-0.93, P ≤ 0.001). We also calculated peat cover (%) from a digital land cover classification71. Here, the native resolution of 100 m was resampled at 50 m to match the resolution of the climatic and topographic predictors, using nearest-neighbor interpolation. The binary peat cover variable was transformed to a continuous scale using a spatial mean filter of 3 × 3 pixels52. Finally, the topmost soil layer of a global gridded soil database72 was used to obtain pH data. Again, the original resolution of 250 m was also resampled to 50 m resolution using bilinear interpolation.

Our fine-scale data revealed strong environmental gradients over the 78,000 km2 study area (Supplementary Table 1), most of which were only moderately inter-correlated (Spearman’s correlation coefficient < [0.7]; Supplementary Table 3). All predictors showed clear spatial autocorrelation structures (Supplementary Figs. 2–3).

Climate data for future warming scenarios

To investigate tundra AGB change along climate warming, we prepared climate data by modifying monthly air temperatures from +0.5 °C to +8 °C at 0.5 °C intervals. To account for the unequal warming among different months, we calculated monthly adjustment factors (Supplementary Table 4). The relative changes in current monthly air temperatures were compared to those projected by an ensemble of 23 global climate models derived from the archives of Coupled Model Intercomparison Project Phase 5 (CMIP5)73 for the period of 2070–2099, assuming a representative concentration pathway (RCP) 8.574. The monthly precipitation data needed to be adjusted to account for the air temperature increase. The amount of adjustment was determined by regressing the projected change in monthly precipitation (\(\triangle\)Prec, %) by the CMIP5 models against the projected change in the corresponding monthly air temperature (\(\triangle\)T, °C; Supplementary Fig. 12). The resulting relationship \(\triangle\)Prec = −0.445 + 3.85 \(\times\) \(\triangle\)T was used to adjust the monthly precipitation data at each temperature scenario. Subsequently, all climate predictors (i.e., GDD, Tfeb and WAB) were recalculated for each corresponding warming scenario.

Future LSP data

We used the recently developed statistical modeling approach17 to predict future LSP distributions at the study area. In brief, this was done by relating each investigated LSP to climatic, topographic, and soil predictors using an ensemble of four statistical techniques ranging from regression (generalized linear models, generalized additive models) to machine learning (boosted regression trees, random forest). The LSP predictions of cryoturbation, solifluction and nivation over the study area (spatial resolution 50 m × 50 m) were produced for each corresponding warming scenario (i.e., from +0.5 °C to +8 °C at 0.5 °C intervals).

Statistical AGB modeling

Prior to modeling, the AGB values were log-transformed to approximate a normal distribution. All statistical analyses were conducted using the R statistical computing environment75. We used boosted regression trees (BRT, as implemented in the R package “gbm”76) to link AGB to the climatic, topo-edaphic and LSP predictors using all 2917 observations and assuming Gaussian error distribution. BRT is a sequential machine-learning method that combines a large number of iteratively fitted regression trees to a single model with improved prediction accuracy29. BRTs automatically incorporate interactions between predictors, are capable of modeling highly complex non-linear systems, and provide an estimate of the predictors’ relative influence in the models. We optimized BRTs using an “out-of-bag” estimator and setting the interaction depth to six, learning rate to 0.01 and bag fraction to 0.75. In order to control for potential differences caused by choice of modeling technique, we also ran parallel analyses using random forest30 (RF) with the maximum number of trees set to 500. To assess whether accounting for LSP improves the AGB predictions, the models were fitted using two configurations:

$${AGB}={GDD}+{WAB}+{T}_{{Feb}}+{TWI}+{Peat}+{RAD}+{pH}+{RA}{C}_{{base}}$$
(6)

[Base model]

$${AGB}= \,{GDD}+{WAB}+{T}_{{Feb}}+{TWI}+{Peat}+{RAD}+{pH}+{RA}{C}_{{LSP}} \\ +{{{{{\boldsymbol{Cryoturbation}}}}}}{{{{{\boldsymbol{+}}}}}}{{{{{\boldsymbol{Solifluction}}}}}}{{{{{\boldsymbol{+}}}}}}{{{{{\boldsymbol{Nivation}}}}}}$$
(7)

[LSP model]

where the RACBase and RACLSP refer to residual autocovariates (RAC). Following the ref. 77, these were calculated from the BRT residuals and then added back to the models as additional predictors to account for SAC. This was done despite a visual inspection of the BRT residuals without the RACs did not reveal evident signs of SAC (Supplementary Fig. 4). However, after calculating semiovariograms over very short distances (≤ 5 km) SAC patterns in the residuals could be detected (Supplementary Fig. 5). The two fitted exponential functions \(\gamma \left(h\right)=0.707(1-{e}^{-h/162.385})\) for the Base model step and \(\gamma \left(h\right)=0.699(1-{e}^{-h/158.249})\) for the LSP model step indicated that the SAC extends to ca. 162 m and 158 m, respectively. Consequently, the RACs were calculated using the function autocov_dist() as implemented in R-package spdep78 (version 1.1-3). For calculating the RACBase and RACLSP we used the above-mentioned neighborhood radius of 162 m and 158 m, respectively, with an inverse distance weighting scheme and symmetric weight matrix (style=“B” in the autocov_dist -function, following ref. 79).

The predictive capability of the models was assessed with a repeated cross-validation (CV) scheme, where the models were fitted 500 times at each round using a random sample of 99% of the data (out of 2917, with no replacement) and subsequently evaluated against the remaining 1%. A distance-block (hereafter h-block, i.e., the minimum distance between calibration and evaluation data) of 158 m (panel b in Supplementary Fig. 5) was specified to omit model calibration data located at the immediate vicinity of the evaluation data, which could lead to overly optimistic CV-statistics due to spatial autocorrelation80. The CV procedure yielded an average of 2873 observations available for model calibration, and 30 for evaluation per CV round. After each CV run the predicted and observed AGB were compared in terms of root mean squared error (RMSE) for both model steps.

To quantify the effects of LSP on AGB, the LSP models were used to predict AGB (using all 2917 observations and 500 bootstrap sampling-rounds) under two conditions. The first included that other environmental predictors (i.e., climate, topography, and soil) and RACLSP were fixed to their median values (constrained to an environmental envelope allowing each LSP occurrence), and each LSP is in turn set to absent (0). The second condition was that environmental predictors and RACLSP were fixed to their median values as mentioned above, but each LSP is in turn set to present (1). The total effects of LSP on AGB over the study area were determined by adding the mean effect sizes of each individual LSP as follows:

$${Total}\,{effect}= \,\left({AG}{B}_{{cryot}0}-{AG}{B}_{{cryot}1}\right)+\left({AG}{B}_{{solif}0}-{AG}{B}_{{solif}1}\right) \\ +\left({AG}{B}_{{nivat}0}-{AG}{B}_{{nivat}1}\right)$$
(8)

Effect sizes for continuous predictors were calculated over the calibration data by comparing the predicted minimum and maximum AGB values while fixing other predictors to their median values.

To obtain spatial estimates, we predicted AGB over the study area at a spatial resolution of 50 m × 50 m, at each temperature scenario based on two model configurations: Base model and LSP model, using all observations. For the predictions, the RACs were fixed to their median values (i.e., zero). This was done in order to avoid additional interpolation over the study area that could increase uncertainty in the consequent analyses. Although the AGB predictions were produced over the entire study area, all comparisons of AGB values between the model steps were limited to areas suitable for the occurrence of LSP (i.e., the LSP realm)17. As the AGB predictions by BRT were in close agreement with those of RF (Spearman’s rank correlation coefficient = 0.96), we only present the BRT results.