Impacts of shrub coverage for arctic ecosystem carbon uptake and storage

Although shrubs employ different water and carbon-use strategies compared to trees and are expanding rapidly across warming tundra and grassland, they are not adequately represented in global land surface models. In this study, we incorporated two types of shrubs, deciduous and evergreen, into the QUINCY nutrient-enabled terrestrial biosphere model, which features a state-of-the-art treatment of soil nutrient dynamics and carbon exchange. We investigated changes in carbon (C) fluxes and storage due to shrub cover, it's response to climate and CO₂ fertilization effect, as well as the role of nitrogen availability.

With this new implementation, shrubs exhibited reasonable seasonal cycle of gross primary production (GPP) at 50 % of the Arctic study sites. The model achieved mean R² values of 0.5 on average, when compared with field measurements and remote sensing products, which represents an improvement for evergreen shrubs compared with the old model parametrization. However, the model-measurement agreement is lower at sites with heterogeneous local vegetation, which is not captured by the model hydrology or represented by the shrub plant functional type. The modelled gross primary production of shrubs is similarly distributed compared to needleleaf evergreen forest, with no significant difference in the median. Compared to graminoids, shrubs' carbon fluxes are 40 % higher. Shrubs produce substantial above-ground biomass, albeit lower than needleleaf trees, and exhibit distinct phenological patterns. While CO₂ fertilization generally benefits all plant types, shrubs appear to have particularly robust growth under elevated CO₂ concentrations.

Providing the plants with unlimited nitrogen doubled gross primary production at most sites and improved the model-measurement agreement by 15 % whereas the modelled deciduous shrubs substantially reduce their nitrogen sources more than evergreen shrubs. A similar effect occurred when nitrogen and carbon contents were initialized based on permafrost profiles, which partly alleviated nitrogen limitation in the model. Overall, this study demonstrates the impact of evergreen and deciduous shrub vegetation on the Arctic carbon cycle and ongoing changes. It recommends including shrub plant functional types in global Earth system models. However, the strong nitrogen limitation of Arctic shrub productivity when using standard model parameterizations suggests that global models underestimate the Arctic's contribution to global land carbon.

1 Introduction

The high latitude ecosystems north of 60° N are vulnerable to climate change. This region is warming at around three times the global rate (Giesse et al., 2024). Arctic soils contain vast amounts of organic carbon (C) in the permafrost (frozen soil). Under climate warming and the thawing of permafrost, this carbon could be released to the atmosphere as carbon dioxide or methane inducing a positive climate feedback loop (Schuur et al., 2015, 2022). An increase of atmospheric temperature prolongs the growing season and deepens the active layer, causing major shifts in vegetation growth and cover in the Arctic (Mekonnen et al., 2021, and references therein), (Liu et al., 2025).

Arctic shrub vegetation is short-statured (mostly under 2 m tall) and often prostrate, it can grow beyond the cold limits of trees due to its adaption to cold environments. Shrubs exhibit faster wood production than trees, enabling them to complete their growth within a growing season that is 20 %–50 % shorter (Treml et al., 2019; Körner, 2012), which is a key adaptation to the Arctic conditions, and intense competition for limited light and nutrients beyond the treeline. Shrubs are classified as either evergreen or deciduous species. Deciduous shrubs grow faster than evergreen species, since their small leaf mass per unit area enables them to intercept more light for a given investment of leaf C than evergreen shrubs (Parker et al., 2021; Mekonnen et al., 2018).

Heijmans et al. (2022) found that shrubs are becoming increasingly dominant in the Arctic. Current estimates suggest that they comprise approximately 40 % of the plant biomass in the Polar Arctic (Orndahl et al., 2025). The expansion rate has been estimated at around 38 % from 1982 to 2014 and 27 % from 2000 to 2020 in the high-northern latitudes, due to increased biomass in existing stands and new establishment (Chen et al., 2020; Andreu-Hayles et al., 2020; Gui et al., 2025). Studies (e.g. Meyer et al., 2021; Elmendorf et al., 2012), mainly report an enhanced coverage of tall deciduous shrubs, such as birch, willow, and alder. This process is often termed as “shrubification”. Recently, the expansion of evergreen dwarf shrubs, which are ecologically distinct from tall deciduous shrub species, has also been observed in the Arctic (Vowles and Björk, 2019). However, these vegetation changes are, however, spatially heterogeneous and are controlled by a mixture of environmental factors other than temperature, such as changes in soil wetness and plant-plant interactions (García Criado et al., 2025; Martin et al., 2022). The ability of shrubs to respond more to warming than e.g. graminoids is linked to their capacity to develop woody tissue, thereby increasing the height and width of their canopy (Bret-Harte et al., 2001, 2002).

Changes in Arctic vegetation cover trigger several feedback loops in the atmosphere. These include decreasing the albedo, thereby increasing the surface temperature, and impacting soil-plant interactions. Examples of this are deepening the active layer, higher mineralisation rates, and improving litter quality. These changes affect the ecosystem's C balance (Bonfils et al., 2012; Meyer et al., 2021; Mekonnen et al., 2021). Furthermore, shifts in Arctic vegetation potentially affect permafrost thaw by altering surface temperatures (Heijmans et al., 2022). Buckeridge et al. (2010) conclude that increased nutrient inputs via enhanced litter inputs from smaller shrub species subsequently cause an increase in the growth of tall shrub tundra species.

Field studies show that shrubs impact soil nutrient cycling in various ways, as they are associated with high leaf litter quality and quantity. This means they contain more nitrogen, which enhances decomposition rates (Martin et al., 2022; Moore et al., 2024, and references therein). This is in contrast to boreal trees. The vertical distribution of soil carbon (C) is determined by nutrition allocation, litter quality, and plant root depth; for example, shrubs store more C in deep soils than graminoids (Du et al., 2024). Besides climate, several other factors influence shrub establishment and growth. The acceleration of nutrient cycling, driven by higher soil temperatures, greater litter inputs, and greater litter turnover, has been found to play a key role in vegetation change. However, despite the importance of these processes the impact of nutrient dynamics on shrub growth and C cycling is neglected in most commonly used vegetation models (Miller and Smith, 2012).

The magnitude of the impact of the rapid shrub expansion in Arctic landscapes on the current and future Arctic C balance is uncertain. Given the importance of shrubs in these landscapes shrub plant functional types (PFTs) must be included in land-biosphere models to more accurately represent C dynamics in the heterogeneous Arctic landscapes. Accordingly, we here employ the QUINCY terrestrial biosphere model (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system v1.0) (Thum et al., 2020) which includes the main features of Arctic ecosystems, such as soil C storage, nutrient limitation and vertical discretisation, to realistically simulate C dynamics in the Arctic and physical processes relevant in high latitudes (Lacroix et al., 2022). Unlike other major community land models, QUINCY represents vegetation with PFTs based on physiology rather than geography, overcoming the need for regional PFT definitions, such as boreal trees and limiting the number of parameters required to a minimum (Thum et al., 2020). To assess the effects of shrub expansion on ecosystem C dynamics across the Arctic, we implemented two shrub PFTs into QUINCY (Sect. 2.1 and 2.3). Using this new model configuration, we performed simulations that took into account coverage of three different vegetation types: boreal forest, graminoids, and shrubs. Additionally, we conducted sensitivity simulations with constant CO₂, with nitrogen supply set to levels that do not limit carbon cycling, and with nutrient initialisation in the permafrost (Sect. 2.2), respectively. These simulations were conducted at 16 different Arctic sites. Through these simulations and two observational datasets (Sect. 2.4), we aimed to answer the following research-specific questions:

1. How much carbon do shrubs absorb and store, and does this align with observed data? – Sect. 3.1

2. How does the response of carbon uptake and storage in shrub-covered ecosystems differ from that in boreal forests and graminoid-dominated tundra when subjected to combined climate warming and CO₂? – Sects. 3.2 and 3.2.1

3. What role does nitrogen cycling play for shrub growth and C storage, particularly in response to increased nitrogen availability from permafrost thaw? – Sect. 3.3

2 Methods

2.1 QUINCY model

The QUINCY terrestrial biosphere model (Thum et al., 2020) simulates the coupled carbon, nitrogen (N) and phosphorus (P), as well as water and energy cycles. It represents physical processes relevant to high latitudes, such as snow dynamics, soil freeze-thaw, and adaptive rooting depths to layers in non-frozen layers (Lacroix et al., 2022). The standard version of the model includes 12 different PFTs with flexible parameterizations for leaf characteristics, phenology, and other growth and turnover dynamics. The concept of PFTs groups all plant species into a limited number of classes, thus allowing for a much more traceable interpretation of the plant interactions (Poulter et al., 2015). The approximately 40 parameterizations representing the different plant characteristics and functionalities in QUINCY are based on plant trait databases (Kattge et al., 2020). Here, we have added additional PFTs to simulate Arctic shrubs.

In QUINCY, vegetation consists of structural pools (leaves, sapwood, heartwood, coarse roots, fine roots, and fruit) and a fast-overturning, respiring non-structural pool (labile), as well as a seasonal, non-respiring, and non-structural storage pool (reserve). The canopy comprises several layers (a multilayer canopy scheme), including separate representations of photosynthesis and canopy conductance for sunlit and shaded leaves separately within each layer. Vegetation processes as well as nutrient uptake and availability respond to local environmental conditions on a process-specific timescale, dynamically enabling a certain memory effect (Thum et al., 2020). Plants partition the required C between growth, respiration and storage, acclimatising to meteorological conditions and nutrient availability.

The soil is explicitly discretized vertically into 15 layers of increasing thickness down to a depth of 9.5 m, and this applies to soil temperature, texture, moisture, and the biogeochemical soil pools. Soil biogeochemical pools are represented by five organic pools: structural; polymeric; and woody litter; and fast and slow-degrading soil organic matter; as well as solute concentrations of ammonium (NH₄), nitrate (NO₃) and phosphorus (PO₄). Each of these organic pools contains C, N, and P (Thum et al., 2020). The model specially accounts for the nitrogen cycle in soil by simulating nitrogen inputs through symbiotic (microbial) fixation in higher plants, as well as free-living (non-symbiotic) N fixation in all parts of the ecosystem. It also accounts for N deposition, and the transformation processes of mineralization, nitrification, and denitrification, as well as the transport between soil layers. N is always accounted for as ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$.

2.2 Site-level simulations setup and forcings

Site-level QUINCY model simulations were conducted for sixteen Arctic sites (Fig. 1 and Table S1 in the Supplement) from 1901 to 2010. Data on shrubs and GPP estimates were provided by the ABCflux dataset version 2 for these sites. In the reference simulation, each site is represented by one of the two shrub PFTs according to the vegetation reported at each site (see Table S1). The simulations include nitrogen dynamics, as mineral nitrogen is often known to limit plant growth in the Arctic (Mekonnen et al., 2018; Martin et al., 2017; Prager et al., 2017). Different model configurations were chosen to quantify the differences between ecosystems covered by graminoids and NE trees at the same sites. A second set of simulations disentangles the effects caused by climate drivers, CO₂ fertilization and nutrient limitation (Table 1).

Figure 1 Map of the 16 chosen stations: BON: Bonanza, US; TVC: Trail Valley Creek, CA; SAM: Sammaltunturi fell, FI; ISK: Iskoras, NO; YAK: Yakutsk, RU; STR: Stordalen, SE; ABI: Abisko, SE; RAI: Raisduoddar, NO; XHE: Neon Healy, US; WLF: Wolf Creek, CA; LET: Lettuoso, FI; SIL: Siikaneva, FI; IGA: Igarka, RU; SEI: Seida, RU; FAI: Fairbanks, US; IMN: Imnavait, US. The black edge indicates the needleleaf and the white edge the broad-leaved shrub.

Friedlingstein et al. (2022)Friedlingstein et al. (2022)Friedlingstein et al. (2022)Friedlingstein et al. (2022)Friedlingstein et al. (2022)

Table 1 Description of the conducted simulations (1901–2010).

The standard initialization of QUINCY does not account for the elevated levels of soil organic C and N characteristic of permafrost soils. These contents cannot be built up in a standard QUINCY spin-up either, as this would require a simulation of the type used for paleo-timescales that reproduces at least the last interglacial and glacial periods. This is neither practical for our purposes, nor feasible given the lack of the required model forcing data. In this study, we therefore adapted the initialization procedure by adjusting the C and N contents, as well as the distribution of the initialization, when a site is identified as permafrost based on the definition and methods in Palmtag et al. (2022). The permafrost distribution curve is initialized using a Weibull function with a shape parameter of 0.6 instead of 2, and the minimum relative deep permafrost concentration is set at 0.3 of the surface soil concentration according to the estimated distribution. Soil organic matter is initialized with a single permafrost content of 121 000 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for the first 3 m, which is the vertically integrated average C. Permafrost nitrogen is initialized using a C:N ratio of 14. It should be noted that the applied approach does not account for spatial heterogeneity in permafrost soil content and distribution during initialization.

The meteorological forcing that drives the model is derived from the University of East Anglia Climatic Research Unit Japanese Reanalysis (CRU-JRA; Harris, 2019), which is a reanalysis of incoming solar radiation, surface temperature, humidity, precipitation, surface pressure, and surface winds. This dataset is available with a spatial resolution of 0.5°×0.5° and at a temporal resolution of 6 h. Data for our selected sites were extracted from the CRU-JRA dataset according to their geographic coordinates (Zaehle and Friend, 2010). Furthermore, the soil is represented by 15 layers, with increasing thickness down to a depth of 9.5 m. The model starts the transient simulations after a spin-up period of 1000 years using cyclic forcing data from 1901–1930. The atmospheric CO₂ concentrations are prescribed according to Friedlingstein et al. (2022).

2.3 New shrub implementation

We implemented two new shrub PFTs: broadleaved deciduous (BD) and needleleaf evergreen (NE). The two shrubs are based on the respective default parametrization of the woody PFT, but due to the physiological and phenological differences important aspects concerning the initial leaf N, plant allometry, root distribution and wood turnover times are specified based on measurement studies and trait databases (e.g. TRY, Kattge et al., 2020). The new parameters are listed in Table 2.

(Kattge et al., 2020)(Kattge et al., 2020)Zhang et al. (2013)Thum et al. (2020)Zhang et al. (2013)Zhang et al. (2013)He et al. (2025)Ma et al. (2021)Fu et al. (2014)Fu et al. (2014)Fu et al. (2014)

Table 2 New parameters for NE and BD shrubs replacing the default needleleaf evergreen PFT. For simplicity, the latter one is referred to as “Tree”.

2.4 Observational data

The ABCfluxnet v2 data, cited in the references Virkkala et al. (2026), Leffler et al. (2026), was recently compiled to provide a dataset of in situ Arctic–boreal terrestrial ecosystem CO₂ flux measurements. The synthesized GPP is given as cumulative flux densities aggregated at the monthly timescale. These methods were employed during the growing season. We selected sites where shrub growth has been observed at multiple points in time. The 16 study sites are located across the northern USA, Canada, Finland, Norway, Sweden and Siberia, as shown in Fig. 1. These sites are in wet, moist or dry tundra, fens, bogs or boreal forests and are vegetated by either a mixture of shrubs, trees and grasses or only shrubs (see Table S1). The measurements were averaged and used across different years to minimize the impact of climate variability. The data set includes eddy-covariance (EC) measurements, which generally represent a 100 m radius footprint of the measurement tower. For some sites, only chamber data exists (Table S1) which can be considered the most highly resolved, plot-scale GPP estimate of a specific vegetation/soil unit.
To further evaluate the simulated GPP, we also use data extracted from the global MODIS/Terra Level-4 products of MOD17A2H (Collection 6) (Running et al., 2015) and MOD17A2HGF (Collection 6.1) (Running and Zhao, 2021) at the site coordinates. Both products provide a cumulative eight-day composite of GPP values at a spatial resolution of 500 m and have been interpolated. GPP is estimated by applying a radiation-use efficiency framework that integrates absorbed photosynthetically active radiation (APAR) with the biome-specific radiation conversion efficiency (ε), which is constrained by meteorological stress. While the MOD17A2H product contains gaps due to the quality of the input data, the MOD17A2HGF product reduces data gaps by employing additional temporal gap-filling based on quality control flags. Both products are subject to uncertainties relating to the quality of the input data, land cover classification, and model assumptions. With a GPP time series that is more continuous in space and time, the MOD17A2HGF dataset introduces an additional source of uncertainty in the form of interpolation, which dampens the signals caused by inter-annual variability in climate or ecosystem disturbances. Extensive validation against eddy covariance measurements shows that MOD17 series products can moderately capture the interannual and seasonal variability in GPP, with R² over 0.6 across different biomes and climate regimes (Zhu et al., 2018; Tang et al., 2015; Wang et al., 2017). We extracted GPP time series from the datasets for the different sites using the nearest neighbour approach.

3 Results and Discussion

3.1 Modelled vegetation productivity vs. observations

The modelled shrubs capture a reasonable monthly mean seasonal cycle of GPP at half of the validated sites (six out of twelve, with a 50 % measured seasonal cycle), when compared to the field measurements in Fig. 2. Modelling with the new NE shrub PFT improves the mean R² value to 0.46 across nine sites, compared to the default (tree) PFT (see the dashed line in Fig. 2 and Table S2 in the Supplement). In contrast, simulations with the BD shrub (7 sites) mostly underestimate the measured seasonal cycle of GPP; the default (tree) PFT therefore agrees better with the measurements (mean R²=0.6). For further evaluation, we compared our model's GPP outputs to GPP MODIS satellite products. MODIS data is widely used to evaluate the carbon cycle in global terrestrial biosphere models and in the Coupled Model Intercomparison Project (CMIP) models (50 % use the JSBACH surface scheme) (Hu et al., 2022), showing a high agreement in the high latitudes. In this study, the MODIS GPP products show a higher correlation with the model, with mean R² values of 0.71 for NE shrubs and 0.50 for BD shrubs (Table S2) than compared to the site measurements. NE shrubs generally improved the agreement with the MODIS data compared to the NE default (tree) PFT (R²=0.59). The MODIS products generally align with field measurements, though they diverge significantly at two sites. At the ISK site (Norway), the GPP of MODIS and the (CN) model agree quite well, whereas the GPP measured in the field is much larger. This may be due to the specifics of the local climate and landscape at the site coordinates used. Specially, half of the field measurements at the ISK site are taken in wet tundra areas (Pirk et al., 2024) whereas however the model soil hydrology tends to be rather dry. At the LET (Finland) site, which has a significant proportion of forest (e.g. Heydel and Tackenberg, 2017) both observational datasets show a higher GPP which the model does not reflect because it represents shrubs only (in this simulation). Another possible source of the model bias is the underestimation of soil moisture at sites with a higher wetness than 25 %, such as the Stordalen (STR, Sweden) site (Holmes et al., 2022). Overall, two out of sixteen sites (Lettuoso (LET) and Stordalen (STR)) are considered not to be representative of the simulated shrub PFTs due to their heterogeneous ecosystems.

Figure 2 Multi-year mean seasonal cycle of GPP [$\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{month}}^{-\mathrm{1}}$] by measurements (field data and MODIS satellite data: MODGF (MOD17A2HGF) and MOD2H (MODIS17A2H) and model with new shrub PFT, tree PFT and graminoid PFT at 16 sites (BON: Bonanza, US; TVC: Trail Valley Creek, CA; SAM: Sammaltunturi fell, FI; ISK: Iskoras, NO; YAK: Yakutsk, RU; STR: Stordalen, SE; ABI: Abisko, SE; RAI: Raisduoddar, NO; XHE: Neon Healy, US; WLF: Wolf Creek, CA; LET: Lettuoso, FI; SIL: Siikaneva, FI; IGA: Igarka, RU; SEI: Seida, RU; FAI: Fairbanks, US; IMN: Imnavait, US) at individual time periods where field data is available.

All in all, the model's performance depends on the chosen configuration, the evaluated shrub PFT, local heterogeneity of vegetation coverage. While this comparison is an important benchmark for the model and its new parametrization, it is limited to monthly GPP values. The structural and phenological plant characteristics can be only discussed quantitatively in the later analysis of this study due to limited data availability and time constraints.

3.2 Present day C uptake and storage

The simulated (multi-site) average annual GPP flux at shrub sites is 371 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, ranging from 132 to 735 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for the time period 1990–2010 (Fig. 3a). These estimates (and their partitioning) are consistent with the data synthesized by Virkkala et al. (2025) who report an average (area-upscaled) Arctic GPP of 300 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for tundra and 416 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for the Arctic permafrost region, although this accounts for all tundra vegetation species (graminoids, shrubs, mosses, and lichen). The annual GPP for shrubs has a much wider interquantile range than for tress and graminoids, i.e. it is more variable (Fig. 3b). The distribution is more symmetrical, with shorter tails and a flatter peak and a lower skewness and kurtosis. This shows the better adaptation of shrubs to the cold Arctic climate, as it has been observed elsewhere (e.g. Mekonnen et al., 2021). Thus, the median (not considering the high values) of the new shrub and tree PFTs are not significantly different, while the average annual GPP increases by 16 % due to the implementation of shrubs. Additionally, Meyer et al. (2021) report a substantially higher GPP at shrubs than at trees when simulating a Canadian site with the CLASSIC model, which aligns well with the site measurements. Shrub dominance instead of graminoid cover would increase the annual mean and median GPP by 40 %, respectively.

Figure 3 Simulated multi-year (1990–2010) annual mean GPP [$\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$] per PFT for 16 sites in the Shrubs, Trees and graminoid (grass) simulations, with the values of two most recent years marked with “+” and (b) the distribution across all sites (BON: Bonanza, US; TVC: Trail Valley Creek, CA; SAM: Sammaltunturi fell, FI; ISK: Iskoras, NO; YAK: Yakutsk, RU; STR: Stordalen, SE; ABI: Abisko, SE; RAI: Raisduoddar, NO; XHE: Neon Healy, US; WLF: Wolf Creek, CA; LET: Lettuoso, FI; SIL: Siikaneva, FI; IGA: Igarka, RU; SEI: Seida, RU; FAI: Fairbanks, US; IMN: Imnavait, US). The first eight bars (BON-RAI) show the results at NE shrubs, the later ones at BD shrubs.

The GPP change of the three PFTs is determined by different N pathways, and the strength of the correlation varies not only across the PFTs, but also between sites. For example, the GPP at NE trees and shrubs increase much more strongly in response to ${\mathrm{NH}}_{\mathrm{4}}^{+}$ uptake than at graminoids, whereas comparing graminoids and trees to BD shrubs reveals a similar slope (Fig. S1 in the Supplement). Thus, among the shrub species, the model simulates faster growth in response to ${\mathrm{NH}}_{\mathrm{4}}^{+}$ uptake for deciduous shrubs, as also reported by Chen et al. (2020). For all shrubs, the change in ${\mathrm{NH}}_{\mathrm{4}}^{+}$ uptake explained 50 %–70 % of the GPP change compared to trees (see R² in Fig. S2 in the Supplement). For NE shrubs, N fixation and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ net mineralization also showed a significant correlation with the GPP changes (Fig. S2). Most of the increase in GPP at evergreen shrubs could be additionally explained by higher biological nitrogen fixation (BNF) compared to the default PFT. Here, the absolute numbers approximately half the magnitude of the plant ${\mathrm{NH}}_{\mathrm{4}}^{+}$ uptake (Fig. S3 in the Supplement). Indeed, Davies-Barnard and Friedlingstein (2020) and Schore et al. (2023) report substantially higher asymbiotic BNF in shrubs than in trees, driven by bacterial activity. The GPP changes in deciduous shrubs are usually lower than in evergreen shrubs, which may explain the low correlation. The GPP differences from graminoids are generally less correlated with changes in the N uptake pathways.

The site-averaged foliage (projection) cover of shrubs (i.e. the percentage coverage of the ground area covered by leaves when viewed from above), is very similar to that of trees (4 % higher for shrubs), indicating slightly higher self-thinning in shrubs (Fig. S8 in the Supplement).

Shrub cover reduced the carbon-use efficiency (NPP$/$GPP), at all sites, from 0.64 (graminoids) and 0.68 (trees) to 0.63 at NE and 0.57 at BD shrubs, indicating lower efficiency in converting GPP to biomass overall (Fig. S5a in the Supplement). This is generally consistent with satellite estimates of the NPP$/$GPP ratio for the three different vegetation types reported by Zhang et al. (2013). Unlike needleleaf evergreen shrubs, most deciduous shrubs had lower net primary production (NPP) than trees (Fig. S5a) following the growth respiration (Fig. S5b in the Supplement). Higher growth respiration indicates a faster growth rate and higher seasonal leaf turnover (Ge et al., 2017).

The new shrub PFTs produced reasonable total aboveground C values, ranging from 0.7 to 2.5 kg C m⁻² (Fig. 4a, site-mean=1.3 kg C m⁻²) between 1990–2010. Given that approximately 50 % of C in plants is found in the aboveground part (Ma et al., 2018) these values translate into biomass estimates consistent with measurements of 0.4–1.5 kg C m⁻² (Poley et al., 2020; Zhao et al., 2021; Sulman et al., 2021; Orndahl et al., 2025).

Figure 4 Multi-year (1990–2010) annual mean total vegetation C, values of two most recent years marked with “+” and (b) distribution of multi-year annual mean total soil C per PFT in [kg C m⁻²] at all sites (BON: Bonanza, US; TVC: Trail Valley Creek, CA; SAM: Sammaltunturi fell, FI; ISK: Iskoras, NO; YAK: Yakutsk, RU; STR: Stordalen, SE; ABI: Abisko, SE; RAI: Raisduoddar, NO; XHE: Neon Healy, US; WLF: Wolf Creek, CA; LET: Lettuoso, FI; SIL: Siikaneva, FI; IGA: Igarka, RU; SEI: Seida, RU; FAI: Fairbanks, US; IMN: Imnavait, US). The first eight bars (BON-RAI) show the results of NE shrubs, the later ones of BD shrubs.

Most of the simulated C here was stored in the living wood of the shrubs followed by the storage in coarse roots and heartwood (Fig. S4 in the Supplement). ). Shrubs produce less sapwood than trees, but have a faster production rate (Fig. S6a in the Supplement). According to e.g. Treml et al. (2019) this is due to the narrower and fewer xylem cells in shrubs. As explained by Parker et al. (2021), Mekonnen et al. (2018) deciduous shrubs grew faster than evergreen species (Fig. S6a), since their small leaf mass per unit area enabled them to intercept more light for a given investment of leaf C than evergreen shrubs. Furthermore, the modelled shrubs stored more C in coarse roots than trees, increasing the coarse root:sapwood ratio by a factor eight. This is consistent with findings by Wang et al. (2016) and a global dataset (including 530 shrub land sites) by Ma et al. (2021). Thus, the total aboveground C storage of shrub vegetation is lower overall than that of trees (at 13 out of 16 sites, Cliff's $\mathit{\delta }=-\mathrm{0.39}$). Indeed, Jobbágy and Jackson (2000) report that trees allocate more nutrients aboveground than shrubs and graminoids.

While the distribution of annual mean vegetation carbon at shrub sites overlaps with that of trees (Fig. 4a), the mean biomass at sites with shrub coverage is 17 % lower than the mean tree biomass, as indicated by the lower carbon use efficiency of shrubs (Fig. S5b). The shrub distribution has a much larger interquartile range and is more skewed towards higher values than tree distribution, which points to the higher variability of shrubs. The total vegetation C of graminoids is significantly lower (effect size: Cliff's δ=1.0) than in both woody PFTs when considering distribution, median and range, as also found by Lett et al. (2004). Graminoids store relatively more C in fine roots and the labile pool (almost 40 %), but woody plants have significantly larger wood and coarse root storage.

The higher fraction of coarse root C stored by shrubs determines the increase in total soil C (Fig. 7a), reaching 14.3 kg C m⁻² compared to other PFTs (vs. 13.2 for trees and 11.1 kg C m⁻² for graminoids) at 81 % of the stations. This small (6 %) increase in soil C due to shrub coverage compared to trees is consistent with observational studies reporting greater belowground than aboveground C with shrub vegetation coverage (Ma et al., 2021). It should be noted that the soil type is a primary factor that can explain both increases and decreases in soil organic content, as a meta-analysis of different shrub species found globally (Li et al., 2016). Parker et al. (2021) even report a substantial C loss by shrub expansion, which is consistent with other studies. At a landscape scale, Zhao et al. (2023) found an increase in soil organic carbon with shrub expansion in the Arctic. Our modelling study shows that, evergreen shrubs tend to conserve more C than deciduous shrubs (first eight bars in Fig. 4), as reported by Mekonnen et al. (2018). In comparison with graminoids, soil C increases by 18 % as substantially more C is stored in coarse roots of shrubs. However, since the model does not yet consider the long-term C storage in permafrost soils, the impact on soil C is only discussed briefly here. Given that shrubs exhibit higher root nutrient absorption rates than trees, as supported by Chen et al. (2020), who analyzed 170 root samples from Arctic shrubs, and as corroborated by our simulations (Fig. S6b in the Supplement) significant impacts on nitrogen cycling can also be expected (Martin et al., 2022). Overall, we demonstrate that the newly modelled shrubs produce reasonable C fluxes and storage while reproducing important physiological processes (e.g. wood formation and root allocation) that differ from those of trees.

3.2.1 Transient changes in uptake and storage

Our transient simulations from 1901 to 2010 show an average annual GPP change of 60 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for shrubs, with a significant increase since 1960 (Fig. 5a), the increase in GPP over time for trees and graminoids is 30 %–50 % weaker.

Figure 5 Multi sites average 15-years rolling mean total changes (a, b) and CO₂ fertilization and climate effect (c, d) in GPP relative to 1901–1911 [$\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$] carried out with the N dynamic and C-only model configuration, respectively.

Figure 6 Multi-year (1990–2010) average seasonal cycle of (a) growing season length, (b) GPP change (Shrubs-ShrubsConly), (c) change of N uptake and net mineralization, including ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Shrubs-ShrubsConly) at particular stations.

Initializing N and C in the permafrost (simulation: ShrubsPM) increased the soil C pool by 10–60 kg C m⁻² (Fig. 7a, 4 %–18 %), towards better agreement with the permafrost C estimates of 20–150 kg C m⁻² (0–3 m), where the majority of the soil C is stored (Hugelius et al., 2014). It also improved the total soil N content from 1.5 to 7.4 kg N m⁻² (Fig. S12 in the Supplement). Palmtag et al. (2022) report estimates of around 5 kg N m⁻², when extrapolated to our modelled soil depths of 9 m using the same vertical distribution.

Figure 7 Multi-year (1990–2010) average total soil C (a) and total soil nitrogen (b) at the 16 different sites estimated by the three simulations: Shrubs, ShrubsConly and ShrubsbPM.

We deduct from the Shrubs and ShrubsNof simulations the CO₂ fertilization effect on shrub GPP across the study sites, with a steep rise in the CO₂ effect on GPP since the 1960s (Fig. 5c) following increased atmospheric CO₂ concentrations (Intergovernmental Panel on Climate Change, 2023). This effect is more pronounced at sites where nitrogen dynamics have a low impact. For example, at YAK (Yakutsk, Russia), where there is no difference between C-only and CN simulations, we found an GPP increase of up to 40 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, caused by CO₂ fertilization alone. Most shrub-dominated sites are strongly nitrogen-limited, which results in a small average CO₂-fertilization effect (Fig. 5a) across all sites between 1960–2010. This estimate aligns with the results of the QUINCY model study for the high-Arctic tundra by Lacroix et al. (2022), as well as with a number of field N fertilization experiments involving tundra and mixed shrub coverage (DeMarco et al., 2014; Shaver et al., 2001). For NE trees, positive changes in GPP outweigh negative changes at only a few sites over time, leading to a small overall change (Fig. 5a). This is consistent with their lower N uptake compared to shrubs. Thus, the newly modelled shrubs demonstrate higher growing potential (N uptake and GPP) under increased CO₂ in comparison to NE trees (default parameters). The CO₂-fertilization effect (on GPP) is higher the less nitrogen-limited the sites are (e.g. the YAK site in Russia). Indeed, it has been found that shrubs deepen their roots in response to increasing CO₂ (Mekonnen et al., 2021).

After 1960, the climate effect clearly becomes dominant due to the significant temperature increase in the Arctic, as reported, for example, by Wei et al. (2018). The multi-site average GPP increases until 45 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ at shrub-dominated sites relative to 1901–1910, which is approximately 12 % of the absolute GPP. The foliage projection cover of trees rises over time, similarly to that of shrubs (Fig. 3). In contrast, the GPP of graminoids shows only half of the increase. Until 1960, the climate impact of all three PFTs exhibited a 20-year cycle, driven by variations in dominant climate factors. Among all three PFTs, shrubs respond most to climate effects in the model, e.g. to increasing temperature since 1960. This feature has been cited as a reason for the shrub expansion in the Arctic, as described in e.g. Mekonnen et al. (2021).

As expected, alleviating N limitation in the (ShrubsConly, ShrubsNOf) simulations resulted in a substantially higher CO₂ fertilization effect (on average 58 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, Fig. 5b). This can be explained by the change in (plant) biological nitrogen fixation (BNF) in plants during an unlimited N supply. BNF increases with temperature (Yu and Zhuang, 2020) and the model formulation used here favors elevated CO₂ (Meyerholt et al., 2016). Thus, BNF increases by 5 % yr⁻¹ with rising CO₂ (Fig. S10 in the Supplement), meaning that by 2010 the CO₂ fertilized C-only simulation had 30 % more BNF than the non-fertilized simulation. BNF is an important N source for plants, whereas the contribution of symbiotic (bacteria) and free-living organism is uncertain (Davies-Barnard et al., 2020). In our C-only model used here plant BNF is unlimited. Thus, a simulation that is twice the size of the observations in the Arctic (see Meyerholt et al., 2016) compensates the absence of other N sources, such as symbiotic BNF (gray line in Fig. S10).

The effect of climate on the total vegetation C (Fig. S9 in the Supplement) is similar to that on GPP, with a sharp increase around 1970 to a change of 0.2 kg C m⁻² (8 % of the total amount). Indeed, numerous studies, such as Elmendorf et al. (2012), have reported an increase in Arctic ecosystem productivity due to climate warming. The CO₂ fertilization has a minor effect on shrubs and graminoids, and even leads to a reduction in tree vegetation C of about 5 %. Comparing CN and C-only simulations for shrubs reveals the same climate effect, as the climate forcing is the same, but a much higher CO₂ fertilization effect occurred in the C-only simulation (0.5 kg m⁻², 25 %). Norby et al. (2010) and McCarthy et al. (2010) confirm that local N availability influences the strength and the sign of the CO₂ fertilization effect on Arctic trees. Changes in vegetation C also affect the foliage projection cover of trees and shrubs. CO₂ fertilization has no significant impact on modelled shrub foliage projection cover (Fig. S8). For trees, however, we see a decrease in foliage cover when accounting for CO₂ fertilization over the last 20 simulation years (Fig. S8, orange compared to coral line) which results in higher self-thinning and thus fewer trees per leaf area. Similar to shrubs, graminoids show only a modest response to CO₂ fertilization (Fig. 5a).

3.3 Role of nitrogen availability

Here, we compare the C-only (ShrubsConly) and CN (Shrubs) simulations, i.e. simulations in which the nitrogen supply is maintained at levels that do not constrain plant growth vs. an explicit simulation of nitrogen input and output processes, as well as the limitation of plant growth due to nitrogen availability. Accounting for nitrogen limitation significantly reduced the GPP at all sites in shrub cover simulations (Figs. 6b and S7 in the Supplement), most notably for deciduous shrubs. These reductions ranged from 3 $\mathrm{g}\mathrm{C}{\mathrm{month}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$ (2 %, not significant) for the evergreen shrub at the YAK site in Russia to 150 $\mathrm{g}\mathrm{C}{\mathrm{month}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$ (75 %) for the deciduous shrub at the IGA site in Russia (Fig. 6a). We demonstrate that the magnitude of GPP change (Fig. 6b) depends on the length of the growing season (Fig. 6a), which impacts the nitrogen supply to the ecosystem through net N mineralization, fixation (not shown) and uptake (Fig. 6c). At sites with deciduous shrubs and a three-month growing season only (at the SEI, IMN, IGA sites; whereas SEI and IGA shrubs behave very similar), these main N sources declined in the CN simulations, leading to a significant decrease in GPP (Fig. 6b). In contrast, shrubs at other sites, e.g. at the WLF site in Canada, showed at least an increase in N uptake at the start of the growing season under the CN configuration. When both N uptake terms and the use of N storage increase (at the sites YAK, STR and ABI, where the shrubs at STR and ABI behave similarly), the reduction in GPP is only slight (Fig. 6b and c). We demonstrated that nitrogen supply strongly limits most simulated shrubs, which is consistent with other studies (Street and Caldararu, 2022). Other than modelled evergreen shrubs, the deciduous shrubs have a higher N demand and respond more rapidly to N additions, as has been reported for alpine shrubs by Zhou et al. (2023). The C-only model, which allows for unlimited biological N fixation, performs better than the CN model when compared to measurements (Table S2). This suggests that the N dynamics in the model may be overly restrictive. Supporting this hypothesis, we extracted 89 measurements of the C$/$N ratio at Betula nana leaves and 36 measurements of the C$/$N ratio at Betula glandulosa leaves were extracted from the TRY database (Kattge et al., 2020). The C$/$N ratio estimated by the C-only configuration at seven stations (with deciduous species) is 24.0 (median), which is consistent with the measured median values (26.06±2.0, 22.56±2.1). However, the median value of the CN model (29.06) instead overestimates the measured C$/$N ratio (i.e. it underestimates the leaf N concentration), thus supporting our hypothesis of overly too strong N limitation in the CN model. We recommend to improve the model by refining key processes such as BNF, to better reflect Arctic conditions.

Higher soil N levels, originating from the improved soil initialization, stimulated primary productivity. This increases the C storage in the coarse root, sapwood and heartwood pools of the plants. Consequently, total vegetation C increased by 1 kg C m⁻² (93 %) on average (Fig. S11 in the Supplement).

For most of the sites with NE shrubs, the effects of the enhanced soil N supply (Fig. S7) are clearly reflected in the annual GPP, which increased by 73–290 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ compared to the standard model simulations (Shrubs). The additional nutrients led to a significant increase in the leaf-to-root mass ratio (Fig. S13 in the Supplement), which indicates significant alleviation of nitrogen limitation. This pattern has been widely reported for other vegetation types and causes the shrubs to allocate more growth to aboveground biomass than to roots. The latter is needed to take up more nutrients. In contrast, GPP for BD shrubs is less strongly affected (Fig. S7), and the leaf-to-root mass ratio at BD shrub sites do not change significantly (Fig. S13b in the Supplement). This feature highlights the contrasting nitrogen sensitivities and acquisition strategies of NE and BD shrubs. Indeed, Martin et al. (2022) report a rapid increase in GPP and leaf-to-root mass ratio in NE shrubs with higher N availability whereas BD shrubs only show a slight response to increased N due to the conservative economics and higher respiration efficiency (Street and Caldararu, 2022).

The change in GPP due to the PM initialization varies significantly over time. From the initialization the GPP difference rapidly increases by 10 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ to a level which is more or less constant until 1960. With increasing CO₂ the GPP difference increases by 30 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ until 2010.

4 Conclusions

We incorporated broadleaf deciduous and needleleaf evergreen shrubs into the QUINCY nutrient-enabling terrestrial biosphere model using empirically based parameters to enhance the representation of Arctic carbon cycle processes. We simulated the growth of shrub, tree and graminoids separately at 16 Pan-Arctic sites separately to assess the impact of Arctic shrub cover on carbon uptake and storage. We found that:

• Modelling shrubs using empirically based parameters enabled a reasonable simulation of gross primary production (GPP), achieving R² values of 0.5 on average when compared with field data and MODIS satellite observations. This approach also accurately captured well-known physiological processes that are distinct from those of trees.

• At the same simulated sites, shrubs exhibited greater growth than NE trees and graminoids, leading to a mean increase in GPP of 16 % or 40 %, respectively.

• Climate-induced effects dominated over other driving factors when simulating enhanced shrub growth from 1901 to 2010. The results imply that this is due to the nitrogen limitation, which was strengthened by CO₂ fertilization, but weakened by a warmer climate. This enhanced nitrogen fixation and soil nitrogen turnover.

• The initialization of soil nitrogen concentrations, parametrized from observation products, enhanced shrub growth, especially at NE sites, improving the model-measurement comparison by 15 %. It also affected the increased C uptake over 1901–2010.

Overall, our results demonstrate that explicitly presenting shrubs significantly changes simulated carbon uptake and storage. This underlines the importance of including evergreen and deciduous shrub plant functional types (PFTs) in global land surface models to accurately predict ongoing changes in the Arctic C cycle. In the light of the significant impact of nitrogen limitation on ecosystem C fluxes as revealed here, we recommend improved constraints on the soil N cycle in low-temperature conditions and its changes in a warming climate. We also call for additional measurements needed to parametrize the PFT-specific traits, such as the biomass allocation to roots, which influence shrubs' efficiency in accessing soil nitrogen (Martin et al., 2022) and are likely to become increasingly important with permafrost thaw.