Introducing shrubs enhances the representation of high-latitude vegetation and carbon cycling in the ORCHIDEE land surface model

Arctic–Boreal terrestrial ecosystems are rapidly changing under amplified high-latitude warming, including widespread expansion of shrubs, with consequences for regional carbon and energy balances. Yet, high-latitude vegetation diversity and vegetation–climate interactions remain under-represented in many global land surface models. In ORCHIDEE, the land surface component of the IPSL Earth system model, high-latitude vegetation is represented primarily as boreal trees or grasslands, omitting explicit shrubs. Here, we implement three high-latitude shrub plant functional types (PFTs) (tall deciduous, low deciduous, and evergreen dwarf shrubs) in ORCHIDEE (revision 9269). Following literature recommendations, this classification combines phenology and stature to capture key functional contrasts while keeping the number of new PFTs limited. The implementation builds on ORCHIDEE's existing woody vegetation scheme by recalibrating a targeted set of parameters controlling allometry, carbon allocation, recruitment, mortality and phenology. Parameter values are constrained using synthesised pan-Arctic observations to obtain regionally representative shrub traits. Shrub spatial distributions are prescribed with updated PFT maps that combine ESA CCI products with Arctic and regional shrub mapping information. The resulting shrub PFTs reproduce observed ranges of shrub size and biomass allocation across the Arctic–Boreal domain. Introducing shrubs reduces simulated total aboveground biomass in the Arctic–Boreal region from 54 to 46.7 Pg C (−13.5 %) and mean annual gross primary productivity from 498 to 481 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (−3.4 %) over the simulated period 1992–2020, with a stronger reduction in the tundra region (4.6 to 3 Pg C (−34.8 %); and 334 to 289 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (−13.5 %)), increasing agreement with benchmarking datasets. A key strength of our implementation is its simplicity, as it builds on ORCHIDEE’s existing woody vegetation framework. In addition, the use of synthesised pan-Arctic observations provides regionally representative observational constraints, making the methodological choices transferable beyond ORCHIDEE. Overall, this work provides a data-constrained shrub representation in ORCHIDEE with minimal added process complexity and establishes a foundation for future development of shrub-climate interactions and dynamic shrubification processes.

1 Introduction

High-latitude regions are warming two to four times faster than the global average under anthropogenic climate change (AMAP, 2021; Rantanen et al., 2022), driving changes in ecosystem structure and function across boreal and Arctic tundra ecosystems (Parmesan et al., 2022). One of the most prominent ecological responses in the Arctic is the expansion of shrubs across tundra ecosystems, often referred to as shrubification (Berner et al., 2020; Myers-Smith et al., 2011). Through a complex set of plant–soil–atmosphere interactions, shrubification processes alter tundra carbon and energy balances, with potential feedback effects on local to global climate (Mekonnen et al., 2021; Heijmans et al., 2022). Shrub expansion can enhance carbon uptake through increases in productivity and biomass, but also increase carbon losses through ecosystem respiration and litter inputs, modify albedo via snow-shrub interactions, alter soil moisture and temperature, and contribute to permafrost degradation (Mekonnen et al., 2021). Due to the many and complex interacting processes involved, the net climate impact of shrubification and its future trajectory remain uncertain (Myers-Smith et al., 2011; Loranty and Goetz, 2012). Reducing such uncertainties requires both long-term field observations and improved process representation in ecosystem models. Consequently, coupled Earth System Models (ESMs) require credible representations of high-latitude vegetation characteristics and dynamics – including shrubs – to simulate observed change and project future scenarios.

However, shrub functional diversity is not always sufficiently represented in global land surface models (LSMs), the terrestrial land surface components of ESMs. Shrub tundra vegetation exhibits strong fine-scale heterogeneity with many structural variations and differing ecosystem functions. Since LSMs operate at much coarser spatial resolution and under computational constraints, plant species are typically aggregated into plant functional types (PFTs) based on similar characteristics and ecosystem function (Wullschleger et al., 2014). For high-latitude shrubs, the most common distinction is between evergreen and deciduous shrubs (Chapin et al., 1996; Wullschleger et al., 2014). However, this categorisation has been criticized as insufficient to capture divergent shrub responses to environmental change, motivating calls for more detailed, trait-informed representations of tundra vegetation in LSMs (Alton, 2011; Wullschleger et al., 2014; Mekonnen et al., 2021). This highlights a central modelling trade-off: PFT classifications must remain general enough for global models while still capturing the functional contrasts that control vegetation dynamics and shrub-climate interactions. To reduce existing uncertainties and biases, the choice of shrub categorisation and key characteristics should be grounded in current ecological knowledge and synthesised field observations.

In response, there have been efforts to improve the representation of high-latitude vegetation in ESM land components, including introducing explicit high-latitude shrub PFTs in LPJ-GUESS (Wolf et al., 2008) and CLASSIC (Meyer et al., 2021), and updating shrub parameterisations in JULES (Harper et al., 2016, 2018). Nevertheless, shrub representations remain uneven across models, and dedicated shrub PFTs are still absent from some widely used LSMs, including ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems), the global land surface component of IPSL-CM.

In its present version (revision 9269), ORCHIDEE only represents high-latitude ecosystems as boreal forest, grassland or bare soil. As a result, regions that are shrub-dominated in reality are often prescribed as forest in the model. The lack of tundra plant diversity in the current model version leads to an inaccurate representation of tundra ecosystems and inherent feedback processes, together with biases in simulated high-latitude carbon stocks and CO₂ fluxes (Zhu et al., 2015).

Previous attempts to address this limitation in ORCHIDEE were made by Druel et al. (2017), who implemented Arctic shrubs, grasses, and non-vascular plants in an earlier version (r1322). Since then, ORCHIDEE has undergone major structural changes, including a completely new carbon allocation scheme with tree demography (Naudts et al., 2015), addition of the nitrogen cycle (Vuichard et al., 2019), and a new multi-layer snow scheme (Wang et al., 2013). As a result, the implementation by Druel et al. (2017) is no longer compatible with the current model structure and could not be maintained in the present ORCHIDEE tag 4.3.

In this study, we introduce shrub PFTs into ORCHIDEE (r9269) to improve the representation of high-latitude vegetation and carbon cycling. To this aim, we: (1) define a classification of shrub types and associated traits constrained by synthesised pan-Arctic observations, (2) implement three new shrub PFTs in ORCHIDEE using existing woody-vegetation functionality, (3) develop updated annual PFT maps that include the new shrub PFTs across the Arctic–Boreal domain, and (4) evaluate the impact of shrub inclusion on high-latitude aboveground biomass and productivity fluxes using independent benchmarking datasets.

2 Methods

Introducing shrubs into a global land surface model requires simplifying shrub diversity into a limited set of characteristics that can be represented with few PFTs. In this study, we specifically focus on the role of high-latitude shrubs in carbon cycle processes in ORCHIDEE LSM.

First, we identified the central characteristics of shrubs and what distinguishes them from trees. Then, we established a meaningful categorisation into shrub PFTs (see Sect. 2.1.1), and determined calibration targets for each type based on observational data (Sect. 2.1.2). Following these data-based shrub characteristics, shrub PFTs were implemented in ORCHIDEE by selecting and optimising a set of PFT-specific parameters (Sect. 2.2). The spatial distributions of the shrub PFTs was prescribed with updated PFT maps (see Sect. 2.3). Finally, pan-Arctic simulations were evaluated against independent observation-based data products (Sect. 2.4).

2.1 Shrub PFT design and observational constraints

We adopted a data-driven strategy in which shrub PFT definition and calibration targets were guided by established ecological understanding and constrained by synthesised observations across the high-latitude region. This makes the shrub implementation approach transferable to modelling frameworks beyond ORCHIDEE, and increases region-scale representativeness of shrub characteristics.

2.1.1 Shrub characterisation and classification

Shrubs are low-stature woody plants – both deciduous and evergreen – adapted to limited growing conditions, such as in arid or cold regions (CAVM Team, 2024). In high-latitude regions, shrubs typically have lower height and stem diameter and therefore lower biomass relative to boreal trees. Furthermore, they are characterised by a multi-stemmed growth form, and a comparatively large belowground allocation to shallow, laterally extensive root systems (Myers-Smith et al., 2015b). Together with rapid regrowth and tolerance of stem loss, these traits may contribute to the ability of shrubs to persist further north than trees (Myers-Smith et al., 2011; Götmark et al., 2016; Treml et al., 2019). These common characteristics and survival mechanisms distinguish shrubs from trees and guide our shrub implementation by indicating which parameters and outputs should be targeted (cf. Sect. 2.2).

Beyond these shared traits, shrubs span a wide range of sizes and growth forms across the high-latitude region. Observed shrub stature ranges from prostrate dwarf shrubs (< 20 cm), through erect dwarf and low shrubs (≈ 20–50 cm), to tall shrubs that can exceed 2 m (Myers-Smith et al., 2015b). Shrub stature and composition vary with climate and local conditions (Yang et al., 2026), but this diversity cannot be resolved in a large-scale model like ORCHIDEE. Shrubs must therefore be grouped into PFTs that capture the functional contrasts most relevant for ecosystem functioning, responses to environmental change and climate interactions (Wullschleger et al., 2014).

Based on evidence from field studies and model experiments we have chosen to distinguish shrub types by phenology and size. We implemented three high-latitude shrub PFTs: tall deciduous shrubs, low deciduous shrubs and evergreen dwarf shrubs. The classification into deciduous and evergreen shrub PFTs follows tundra PFT recommendations and implementations in other models (Chapin et al., 1996; Wullschleger et al., 2014; Dobbert et al., 2021; Meyer et al., 2021). The approach was refined by adding a distinction by height, because shrub stature influences canopy structure, productivity, and snow interactions, and phenology alone is insufficient to capture observed divergent growth responses to climate change (Bonfils et al., 2012; Wullschleger et al., 2014; Saccone et al., 2017; Mekonnen et al., 2021).

The new shrub PFTs in ORCHIDEE represent (i) tall deciduous shrubs, such as willow (Salix spp.) and alder (Alnus spp.), which are often dominant under shrubification (Myers-Smith et al., 2011; Mekonnen et al., 2021); (ii) low deciduous shrubs, mainly representing dwarf birch (Betula nana), which are likewise expanding, but differ in their ecological function compared to tall deciduous shrubs through their lower stature (Bonfils et al., 2012; Sweet et al., 2014); and (iii) evergreen dwarf shrubs, which are expanding under different conditions, with distinct effects on ecosystems and climate feedbacks (Vowles et al., 2017; Vowles and Björk, 2019; Myers-Smith et al., 2011).

This categorisation captures functionally distinct trait combinations while keeping the number of PFTs manageable for data-driven calibration of a global model. The proposed PFTs represent dominant tundra shrub species, align with available calibration data, and support regionally representative parameter values.

Raynolds et al. (2019)Berner et al. (2015)Berner et al. (2024)Wang et al. (2016); Iversen et al. (2015)

Table 1 Data-based characterisations and calibration constraints for the three shrub PFTs, given as representative values and typical ranges. (^∗ Data sources for height: Myers-Smith et al. (2015b); Elmendorf et al. (2012); Raynolds et al. (2019); Macander et al. (2017); Walker and Maier (2008); Wolf et al. (2008).)

2.1.2 Observational constraints

Typical characteristics to constrain the calibration of the shrub PFTs (in Sect. 2.2) were compiled from various synthesized pan-Arctic field observations and literature. Since Sect. 2.1.1 identified low height and biomass, growth form with multiple stems, and a large fraction of belowground biomass as the most distinct shrub characteristics, we focused on observations of height, diameter, above- and belowground biomass and carbon flux measurements (see Table 1). Contrary to the approach taken in other shrub model implementation studies (e.g. Meyer et al., 2021), we chose not to calibrate the model against observations collected at a single field site. Instead, we relied on field observations and literature from across the high-latitude region, to capture the most representative characteristics for the simulated shrubs. This approach was used to enhance the chances that the shrub implementation reflected region-scale patterns rather than local tuning.

The three shrub PFTs are distinguished by height into a tall shrub PFT (aim for ≈ 1.5 m), a low shrub PFT (≈ 0.5 m), and a dwarf shrub PFT (≈ 0.2 m, Table 1). These values represent a compromise between different threshold values from literature (Myers-Smith et al., 2015b; Elmendorf et al., 2012), different shrub mapping products (Raynolds et al., 2019; Macander et al., 2017; Walker and Maier, 2008) and the LPJ-GUESS LSM (Wolf et al., 2008) (see Sect. A3 for further details).

The desired diameters were derived from shrub height based on allometric relationships established by Berner et al. (2015), and are in line with a small survey we conducted in Kobbefjord, Greenland (Table A2).

Figure 1 Locations of observational data sites. Orange triangles: Six EC CO₂ flux sites from FLUXNET and AmeriFlux. Purple dots: 331 shrub-dominated aboveground biomass sites in 24 different areas from Berner et al. (2024). Red rectangle: Model grid cell used for parameter optimisation. Green shading indicates the Arctic–Boreal region as defined in Olefeldt et al. (2021a, b), with the tundra (light green) and boreal (darker green) subregions (Dinerstein et al., 2017).

Ranges for aboveground biomass (AGB) were derived from the Arctic Plant Aboveground Biomass Synthesis Dataset (Berner et al., 2024), containing measurements of peak summer biomass collected at 636 field sites across the Arctic region between 1998 and 2022. Only sample plots with at least 80 % shrub AGB were selected, resulting in 331 sites located in 24 different areas (see Fig. 1) with a median (interquartile range) biomass of 142.5 (83.25–265.0) g C m⁻². The dataset does not provide a categorisation of shrub types, but we used the vegetation description for each plot to extract specific aboveground biomass estimates for the three shrub PFTs by filtering for species or references to their height (see Sect. A3 for further details). We summarised biomass by field site (median across plots) before calculating statistics across sites and converted from dry biomass to carbon assuming 50 % C content. This resulted in the following AGB calibration targets: tall deciduous shrubs (14 sites): 583.5 (325.25–671.3125) g C m⁻², low deciduous shrubs (65 sites): 220 (130–479) g C m⁻², and evergreen dwarf shrubs (238 sites): 131 (74–218) g C m⁻².

Fractions of belowground biomass were determined from Wang et al. (2016) and Iversen et al. (2015). Average fractions of belowground biomass were calculated to be 68.3 ± 9.2 % (mean ± standard deviation) for shrub-dominated sites (> 80 % of AGB) from the Wang et al. (2016) dataset, and Iversen et al. (2015) report 73 % for deciduous and 76 % for evergreen shrubs, which were combined into a target of approximately 70 % (60 %–80 %) of biomass allocated belowground for all three simulated shrub PFTs.

Table 2 Eddy-covariance CO₂ flux measurement sites located in shrub tundra locations from FLUXNET and Ameriflux. Columns: IGBP = IGBP vegetation classification (CSH = closed shrubland, OSH = open shrubland). Veg = Vegetation classification by Oehri et al. (2022) derived from site descriptions based on CAVM (Walker et al., 2005) (S2 = Low-shrub moss tundra, P2 = Prostrate/Hemiprostrate dwarf-shrub lichen tundra, G4 = Tussock sedge, dwarf-shrub, moss tundra). E = elevation [m]. T = Mean annual temperature [°C]. P = mean annual precipitation [mm]. GPP = median (interquartile range) annual GPP [$\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$]. GPP estimated using nighttime flux partitioning method.

Data documenting the carbon fluxes between the atmosphere and shrub vegetation were extracted from the FLUXNET2015 (Pastorello et al., 2020) and AmeriFlux (Chu et al., 2023) datasets. We selected sites classified as open or closed shrubland (OSH/CSH), located above 50° N and with the full FLUXNET-standard post-processing available. Based on the available vegetation descriptions, further sites were excluded, including burn sites with fire succession vegetation, and one tussock-tundra site (US-ICt), resulting in six selected shrub sites (cf. Table 2). It is important to note that due to the nature of EC flux measurements, the recorded fluxes represent the aggregate of the entire ecosystems around the tower within a variable footprint, and the data does not specify the contributions of individual PFTs. Therefore, the sites were not assigned to specific shrub PFTs, but used instead as a benchmark for shrub tundra ecosystems to calibrate and evaluate all three ORCHIDEE shrub PFTs against. Notably, the available vegetation descriptions indicate that dwarf and low shrubs seem to be more prominent than tall shrubs at most sites (e.g. van der Molen et al. (2007); NEON (National Ecological Observatory Network) (2026); cf. Table 2). Across all six sites, a range of average daily gross primary productivity (GPP) of 0.73 (0.18–2.45) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (median (interquartile range)) and average annual GPP of 336 (209–413) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ was extracted. The GPP statistics were calculated using nighttime flux partitioning, after excluding negative GPP estimates and years without positive values.

2.2 Implementation of shrub PFTs in ORCHIDEE

2.2.1 ORCHIDEE model description

We used ORCHIDEE (revision 9269), the land surface component of IPSL-CM (Peylin et al. (2026), see Krinner et al. (2005), Naudts et al. (2015), Vuichard et al. (2019) for descriptions of previous model versions, and supplementary Sect. A1 for an overview of the version used here). ORCHIDEE simulates exchanges of energy, water, and greenhouse gases and explicitly represents terrestrial C and N cycles. ORCHIDEE r9269 builds on ORCHIDEE tag 4.3, which forms the basis for upcoming model applications, including CMIP-7 FAST-TRACK and TRENDY (Peylin et al., 2026).

In ORCHIDEE, vegetation is represented by PFTs, and many parameters are PFT-specific to capture functional differences among ecosystems. The baseline ORCHIDEE r9269 configuration represents global vegetation using 15 PFTs, mainly focused on different types of forest ecosystems, and includes no dedicated shrub PFTs (see Table A1). High-latitude vegetation is represented by five PFTs representing boreal forests, grassland and bare soil. The spatial occurrence of PFTs is prescribed with PFT maps, based on satellite-derived maps of land cover classes (Harper et al., 2023), further combined with information from the land use harmonisation database (LUH; Hurtt et al., 2011) to reconstruct the historical evolution of PFTs back to 1850. Due to the lack of shrub PFTs, ORCHIDEE simulates boreal forest (80 %) and boreal grass (20 %) at high-latitude locations where shrub cover has been observed.

To improve ORCHIDEE's representation of high-latitude ecosystems, we introduced three high-latitude shrub PFTs – tall deciduous shrubs, low deciduous shrubs, and evergreeen dwarf shrubs – into ORCHIDEE (r9269). Shrubs were implemented using existing woody-PFT functionality by recalibrating a targeted set of parameters controlling allometry, recruitment and mortality, carbon allocation, and phenology. Deciduous shrubs were derived from the boreal summergreen tree PFT (PFT 8), and evergreen dwarf shrubs from the closest phenological analogue temperate broad-leaved evergreen tree PFT (PFT 5).

Table 3 ORCHIDEE PFT-parameters altered to implement shrub PFTs based on tree PFTs, their function in the model, and the direction of change for shrubs (− = decrease, + = increase).

2.2.2 Model experiments

Calibration of the shrub PFT parameters was performed in two steps. First, manual sensitivity tests were conducted to identify influential parameters and value ranges (Sect. 2.2.3). Then, the selected parameters were optimised against observational calibration targets using the ORCHIDEE data assimilation system ORCHIDAS (Sect. 2.2.4). Sensitivity tests and optimisation were performed separately for each shrub PFT at grid cell level, primarily using the 2° × 2° grid cell around Toolik Lake, Alaska (67.6–69.6° N, 148.6–150.6° W). This area was found suitable both because it includes a well-studied diverse shrub tundra environment, and because initial ORCHIDEE simulations across the high-latitude region showed that it represents average growing conditions for the region. Different grid cells were tested during the development process, and parameterisation proved to be only weakly sensitive to grid cell-choice. While parameterisation and optimisation were limited to a single grid cell for computational reasons, they were based on regionally representative targets derived from synthesised pan-Arctic observations.

Site-level simulations consisted of a 340 year spinup (using cyclical 1901–1920 forcing), followed by a historical simulation (1900–2020). By using a semi-analytical spinup method, soil C/N pools could be brought to equilibrium with only 340 years spinup (Lardy et al., 2011). Atmospheric forcing conditions were supplied from CRU JRA v2.4 reanalysis, in line with the TRENDY model intercomparison experiment protocol (University of East Anglia Climatic Research Unit and Harris, 2023; Sitch et al., 2024).

2.2.3 Parameter selection and changes

Manual sensitivity tests at site-level determined a suitable set of parameters to target and the direction of their change to capture shrub characteristics starting from tree PFTs (see Table 3). To match the data-based shrub characteristics, the main aims were to decrease height, diameter and biomass, increase the belowground biomass allocation, and keep plant stands stable throughout those changes. The sensitivity analyses were performed by varying parameter values and evaluating their impact on relevant output variables. In some cases, parameters were varied individually using incremental adjustments. In other cases, multiple parameters were modified simultaneously, for example to account for changes in the balance between recruitment and mortality. The sensitivity analyses were partly iterative, requiring certain parameters to be re-adjusted after others had been constrained in order to maintain internal process consistency. Once model behaviour converged towards the predefined calibration objectives, formal optimisation with ORCHIDAS was initiated. In the subsequent paragraphs, the text in italics refers to the parameter names in ORCHIDEE r9269.

Plant allometry

Plant height and diameter were reduced by decreasing the initial diameter of saplings (qmd_init) and the reference plant height at a diameter of 1 m (pipe_tune2). In the current ORCHIDEE version, height can either be controlled statically through (pipe_tune2), or adjusted dynamically based on precipitation. Shrub PFTs use the static height, since precipitation is not the dominant driver for high-latitude shrub growth (Myers-Smith et al., 2015a).

Recruitment and mortality

Besides the influence of allometric parameters, the average stand diameter, height and biomass are the result of a balance between recruitment and mortality. We increased recruitment to stimulate self-thinning mortality and maintain dense, low-stature shrub stands.

In ORCHIDEE, natural mortality of forests is simulated through a self-thinning mechanism based on site carrying capacity (Reineke, 1933; Yoda et al., 1963):

$$\begin{array}{}\text{(1)}& d^{\text{ind},max}={\left({\displaystyle \frac{d^{\text{qmdia}}}{\text{alpha\_self\_thinning}}}\right)}^{\frac{\mathrm{1}}{\text{beta\_self\_thinning}}},\end{array}$$

where $d^{\text{ind},max}$ is the maximum stand density (maximum number of individuals of a given quadratic mean diameter (d^qmdia)), that a site can support, and “alpha_self_thinning” and “beta_self_thinning” are PFT-specific parameters. Self-thinning occurs when actual density exceeds the maximum ($d^{\text{ind},max}$). To allow higher shrub densities, we increased “alpha_self_thinning” for shrub PFTs and prescribed it as constant.

We increased natural mortality through increased recruitment. Recruitment was increased by tuning the PFT-specific parameters (“recruitment_alpha”, “recruitment_beta”) controlling yearly recruits (d^ind,rec; m⁻²) as a function of light availability ($f^{\text{Pgap},\text{trees},\text{season}}$, calculated as in Naudts et al., 2015) and stand density (d^ind):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& d^{\text{ind},\text{rec}}\\ & ={\text{recruitment\_alpha}}^{\left(\begin{array}{c}\log (\text{recruitment\_beta}\cdot d^{\text{ind}})\\ \times \sqrt{\log (max(\mathrm{100}\cdot f^{\text{Pgap},\text{trees},\text{season}},\mathrm{1}\left)\right)}\end{array}\right)}\end{array}\end{array}$$

We also reduced the prescribed height of recruits (“recruitment_height”) to introduce smaller saplings. Increased recruitment and self-thinning shifted stands towards smaller plants. Increased mortality emulates difficult Arctic growing conditions, and higher density emulates multi-stemmed shrub growth and maintains realistic biomass.

Carbon allocation

Belowground biomass fraction was increased by decreasing “alloc_min” and “alloc_max”, reducing allocation to aboveground sapwood and investing more in belowground organs and coarse roots.

We increased the leaf-to-sapwood area ratio through the parameters “k_latosa_min”, “k_latosa_max”, to stabilise growth after biomass reductions, and reduced specific leaf area (SLA, the ratio of leaf surface area to mass) through “sla_init”, to control leaf area index (LAI).

Phenology

For low deciduous shrubs, we advanced growing-season onset by recalibrating the offset parameter in the GDD–NCD phenology model. The GDD-NCD model calculates budburst as a function of growing degree days (GDD) and chilling days (NCD) (Chuine, 2000; Hänninen and Kramer, 2007; Orlandi et al., 2004):

$$\begin{array}{}\text{(3)}& {\text{GDD}}^{\text{thres}}={\displaystyle \frac{\text{gddncd\_ref}}{e^{\text{gddncd\_curve}\cdot \text{NCD}}}}-\text{gddncd\_offset},\end{array}$$

where GDD^thres is the temperature-sum threshold for budburst, NCD is the number of chilling days accumulated over the chilling period, and “gddncd_ref”, “gddncd_curve” and “gddncd_offset” are model parameters. These parameters are not PFT-specific and were originally calibrated against satellite phenology data at global scale (Botta et al., 2000), with the default values “gddncd_ref”: 964, “gddncd_curve”: 0.0058, “gddncd_offset”: 12.8. We increased “gddncd_offset” only for low deciduous shrubs, to lower the onset threshold and improve agreement of simulated growing season timing with shrub tundra observations.

Others

For evergreen dwarf shrubs, originally based on a temperate tree PFT, we adapted key physiological parameters to reflect boreal conditions. Specifically, the reference temperature used to calculate critical leaf age (“leaf_age_crit_tref”) was reduced and nitrogen use efficiency was increased via the parameter “nue_opt”, to stabilise growth in nitrogen-limited high-latitude locations, in line with estimated values reported by Kattge et al. (2009). “nue_opt” represents the photosynthetic nitrogen use efficiency, defined as the carboxylation capacity at 25 °C ($V_{max}^{\mathrm{25}}$) per leaf nitrogen content, based on Kattge et al. (2009), and is used for calculation of $V_{max}^{\mathrm{25}}$ after accounting for sugar loading and leaf efficiency.

2.2.4 Parameter optimisation (ORCHIDAS)

After identifying influential parameters and plausible ranges (Sect. 2.2.3), shrub PFT parameter values were formally optimised against the observational targets determined in Sect. 2.1.2 using the ORCHIDEE data assimilation system ORCHIDAS. ORCHIDAS is an optimisation framework developed at Laboratoire des Sciences du Climat et de l'Environment (LSCE) for the ORCHIDEE model (ORCHIDAS team, 2024; see for example Bastrikov et al., 2018). We used a random-based Monte Carlo approach (Genetic Algorithm (GA); Goldberg, 1989) to identify parameter sets that minimise a cost function J(x) between model outputs (M(x)) and observational target ranges (Y). The observational targets (Y) included height, peak summer AGB, minimum fraction of belowground biomass, mean GPP, peak summer leaf area index, and age (ref. Table A4). The cost function can be written as:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}J\left(x\right)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\right(Y-M\left(x\right){)}^{\mathsf{T}}{R}^{-\mathrm{1}}(Y-M(x\left)\right)+(x-x_{\mathrm{b}}{)}^{\mathsf{T}}\\ & \times B^{-\mathrm{1}}(x-x_{\mathrm{b}})]\end{array}\end{array}$$

and it also measures the misfit between the parameter values x and prior information on them, x_b. Both mismatch terms are weighted by the prior error covariance matrices on observations R and parameters B, respectively (ORCHIDAS team, 2024).

Figure 2 Methodoloy used to include three shrub PFTs into the annual PFT maps for ORCHIDEE (1992–2020). Based on annual data by Harper et al. (2023), supplemented with shrub cover from the Circumpolar Arctic Vegetation Map (CAVM; Raynolds et al., 2019), and re-partitioned following functions derived from Macander et al. (2017) and Macander et al. (2022). Different approaches taken for the Boreal and Tundra subregions.

Similar to the sensitivity analysis, the optimisations were performed at grid cell level, mainly using the Toolik Lake grid cell. For each shrub PFT, a separate optimisation was performed, with prior best-guess parameter values and bounds derived from the initial sensitivity analyses. The GA was configured to use a population pool with at least 30 candidate parameter sets evolved over 30 generations following the principles of genetics and natural selection, and stopped after 30 iterations or when the cost function stabilised, which gives around 30 × 30 = 900 independent model evaluations across each optimisation run. Further details on the optimisation setup can be found in the appendix (Sect. A4).

2.3 Prescribed vegetation distributions (PFT maps)

In ORCHIDEE r9269, the spatial occurrence of PFTs is prescribed annually through PFT maps, derived from ESA CCI satellite-based land cover products (1992–2020) (Harper et al., 2023). Prior to the implementation of shrub PFTs in ORCHIDEE, the high-latitude shrub cover in Harper et al. (2023) was redistributed at grid cell level to boreal tree PFTs (80 %) and boreal grass (20 %), resulting in biased model outputs over the high-latitude region. We updated ORCHIDEE PFT maps to include three high-latitude shrub PFTs by combining multiple mapping products and applying simple partitioning rules (see Fig. 2).

In a first processing step, the total shrub cover was set, based on Harper et al. (2023) shrub cover. Since the high-latitude shrub mapping in Harper et al. (2023) has known caveats, including a general underestimation of shrub cover in the Arctic due to methodological constraints (Harper et al., 2023), we supplemented it with shrub cover from the Circumpolar Arctic Vegetation Map (CAVM) (Raynolds et al., 2019) in the Arctic tundra region. Where shrub cover on the Circumpolar Arctic Vegetation Map exceeded shrub cover on Harper et al. (2023), the shrub fraction was increased at the expense of boreal grass.

Subsequent processing steps partitioned the total shrub fraction into deciduous and evergreen shrubs. Within the CAVM domain (Arctic tundra), we mapped CAVM shrub classes to ORCHIDEE shrub PFTs as follows: P1 and P2 (prostrate shrub tundra) and S1 (erect dwarf shrub tundra) were assigned to evergreen dwarf shrubs (PFT 18), while S2 (low shrub tundra) was assigned to deciduous shrubs (PFTs 16 and 17). To avoid discontinuities introduced by the dominant-class mapping approach in CAVM, we gap-filled the deciduous fraction using a latitude-dependent relationship. The fraction of deciduous shrubs of total shrub cover, $f_{{\text{ds}}_{\mathrm{t}}}$, was derived from a linear regression of the fraction of S2 to total CAVM shrub cover between 59 and 74° N:

$$\begin{array}{}\text{(5)}& f_{{\text{ds}}_{\mathrm{t}}}=\left\{\begin{array}{ll}-\mathrm{0.0106}\times \text{lat}+\mathrm{0.904},& \text{if S2}=\mathrm{0}\\ f_{\text{S2}},& \text{if S2}>\mathrm{0}\end{array}\right.\end{array}$$

where lat is latitude in degrees north and f_S2 is the deciduous shrub fraction derived directly from mapped S2 cover.

South of the CAVM domain, shrub cover was prescribed by combining the four shrub types that are mapped in Harper et al. (2023) and redistributing the total shrub fraction between deciduous and evergreen shrubs, to correct the unequal distribution of shrub types in Harper et al. (2023). The deciduous shrub fraction south of the CAVM domain, $f_{{\text{ds}}_{\mathrm{b}}}$, was set to approximately 70 %, increasing over time to reflect observed shrubification trends, using a relationship derived from Macander et al. (2022):

$$\begin{array}{}\text{(6)}& f_{{\text{ds}}_{\mathrm{b}}}=\mathrm{0.6879}+\mathrm{0.000405}\times \left(\text{year-1992}\right)\end{array}$$

In the final processing step, deciduous shrubs were further partitioned by their height. Across the entire Arctic–Boreal domain, deciduous shrubs were divided between tall deciduous shrubs (PFT 16) and low deciduous shrubs (PFT 17) as a function of latitude to represent decreasing shrub height with increasing latitude (Pajunen, 2009). Latitude itself is not a driver of shrub height but was used as a proxy for growing season length and summer temperatures which correlate well with vegetation height (Gould et al., 2003; Myers-Smith et al., 2015a). We used Macander et al. (2017) data, which map tall (> 1.5 m) and low deciduous shrubs (0.2–1.5 m) in Arctic Alaska, aggregated the mapped region in 0.2°-wide latitudinal bands, and derived a linear relationship for the tall-shrub fraction of deciduous shrub cover, f_tds, from 12.48 % at 68.6° N to 7.47 % at 70.2° N:

$$\begin{array}{}\text{(7)}& f_{\text{tds}}=-\mathrm{0.03125}\times \text{lat}+\mathrm{2.26845}\end{array}$$

Values were clipped to [0,1], implying a northern limit for tall shrubs where the function reaches zero.

Finally, all PFT fractions were normalized to ensure their sum is equal to 1.

2.4 Model evaluation

2.4.1 High-latitude model simulations

To evaluate the new shrub PFTs and their impact on high-latitude carbon cycling, we ran regional simulations north of 50° N at 2°×2° spatial resolution. Like the site-level simulations (Sect. 2.2.2), the model was spun up for 340 years using cyclical 1901–1920 CRU JRA v2.4 forcing (University of East Anglia Climatic Research Unit and Harris, 2023), but followed by a transient simulation for 1861–1900 with increasing atmospheric CO₂ concentrations. Subsequently, the TRENDY protocol (Sitch et al., 2024) was used to set-up a historical simulation for 1900–2020 forced by CRU JRA v2.4. From 1992, vegetation fractions were prescribed annually using PFT maps. For the simulation years before 1992, vegetation was set constant to the 1992 distribution. We ran a baseline simulation with ORCHIDEE r9269 with 15 PFTs and the default ORCHIDEE PFT maps (ORCHIDEE-15), as well as a simulation including the new shrub PFTs (ORCHIDEE-18), with the updated PFT maps (Sect. 2.3), but otherwise identical settings. The simulated years 1992–2020 were used for model evaluation.

2.4.2 High-latitude benchmarking data

The outputs of the regional simulations were evaluated against high-latitude data products to benchmark aboveground biomass and gross primary productivity. The evaluation was limited to the Arctic and Boreal (Arctic–Boreal) region, defined like the BAWLD (Boreal-Arctic Wetland and Lake Dataset) domain in Olefeldt et al. (2021a). This domain was subdivided into boreal and tundra zones following Dinerstein et al. (2017) for refined evaluation (see Fig. 1).

Simulated aboveground biomass across the Arctic–Boreal region was evaluated against two field- and satellite-based data products, by Spawn et al. (2020) for 2010 and Liu et al. (2015) for 1993–2012, remapped and converted to the 2° × 2° resolution used for the ORCHIDEE simulations using conservative remapping. Due to data availability (Spawn et al., 2020), 2010 was used as an example year for evaluation.

To ensure a robust evaluation of regional GPP, we benchmarked our simulations against an ensemble of three independent data-based products, all remapped to the ORCHIDEE grid. These were deliberately chosen to represent distinct methodological approaches: (1) the FLUXCOM product (Jung et al., 2019), representing a machine-learning approach that upscales FLUXNET eddy-covariance observations using remote-sensing and meteorological predictors; (2) the Virkkala et al. (2021) dataset, representing a statistical modelling approach to upscale synthesised CO₂ flux data; and (3) CARDAMOM (Hugelius et al., 2024; López-Blanco et al., 2019), representing a model-data-fusion approach that assimilates remote sensing and gridded datasets into a structurally simple ecosystem carbon model, allowing for observation-constrained parameter retrieval. By utilising this structurally diverse ensemble, we mitigate the risk of evaluating our model against the biases or scaling assumptions inherent to any single upscaling method.

Table 4 Selected PFT-specific model parameters, their function in ORCHIDEE, original values in the underlying forest PFTs and recalibrated values for the new shrub PFTs. The given values are the final values after parameter optimisation, except the parameter values marked with an asterisk, which were set manually. “–” indicates that a parameter value was not changed from the baseline. PFT 8, 16 and 17 are all based on metaclass (MTC) 8, while PFT 5 and 18 are based on MTC 5. All PFT-specific parameters not listed here were kept at the same values for shrub PFTs as in the underlying forest PFTs (8 and 5, for deciduous and evergreen, respectively).

Table 5 Simulated shrub structural attributes and biomass allocation for ORCHIDEE shrub plant functional types across the Arctic–Boreal region given as median (interquartile range) (1992–2020).

Figure 3 Simulated (green, ORCHIDEE-18, 1992–2020) and observed (light brown, Berner et al., 2024) aboveground biomass for three shrub PFTs across the Arctic–Boreal region. Observed biomass from Berner et al. (2024) was filtered to match the ORCHIDEE shrub PFT classification. Numbers in parentheses indicate the number of field sites included after filtering. Boxplots show median (center line), interquartile range (boxes) and whiskers extend to the farthest data point lying within 1.5× the interquartile range. A complete visualisation including outliers is shown in Fig. B1.

3 Results

3.1 Evaluation of shrub PFTs compared to field observations

Introducing three shrub PFTs and optimising key parameters against observational targets produced shrub stands with realistic structure and allocation dynamics across the Arctic–Boreal region. Up to 14 parameters were re-calibrated per shrub PFT to capture their growth form and carbon allocation. Table 4 summarises the optimised parameters and their values for the baseline tree PFTs and the new shrub PFTs after optimisation. The parameters with the highest change from tree PFTs include the controls on plant size, both directly (PIPE_TUNE2, up to 89 % reduction) and indirectly via recruitment (RECRUITMENT_BETA, up to 100 % reduction), which were significantly decreased to capture the difference in stature between tree and shrub PFTs.

The simulated shrub characteristics (Table 5) fall within the range determined from reported in situ field observations (see Table 1). For instance, median heights (1992–2020) are: 1.37 m (tall deciduous), 0.47 m (low deciduous) and 0.23 m (evergreen dwarf), and basal diameters 4.84 cm (tall deciduous), 2.23 cm (low deciduous) and 0.47 cm (evergreen dwarf), closely matching the literature-based aims established for the shrub PFTs. For all three shrub PFTs, a fraction of ∼ 64 %–66 % is allocated to belowground biomass, in line with ecological knowledge (Table 5). Overall, our estimates provide a physically plausible shrub baseline for evaluating high-latitude carbon stocks and flux dynamics in regional simulations.

3.1.1 Aboveground biomass

Simulated aboveground biomass differs between shrub PFTs, and closely aligns with Berner et al. (2024) observations (Fig. 3). Tall deciduous shrubs in ORCHIDEE (median (IQR): 420.48 (335.8–535.6) g C m⁻²) fall within the range of measured AGB at sites dominated by tall willow or alder shrubs (583.50 (325.2–671.3) g C m⁻²). Simulated low deciduous shrubs (183.74 (155.0–224.1) g C m⁻²) closely match observations at dwarf birch-dominated sites (219.50 (130.4–478.6) g C m⁻²). Modelled evergreen dwarf shrubs show similarly low AGB (126.77 (105.0–147.3) g C m⁻²) as observed at field sites dominated by evergreen dwarf shrubs species (131.25 (73.8–217.5) g C m⁻²). For each shrub type, observed variability from Berner et al. (2024) is larger than simulated, which is expected both because of the limited information available to accurately filter the observational data, and the limitations of model PFTs to fully represent observed heterogeneity.

Figure 4 Observed (black) CO₂ fluxes at six shrub-tundra EC sites and simulated (color) shrub CO₂ fluxes at four ORCHIDEE grid cells (columns, where the coordinates above the figures represent the centre of the 2° × 2° grid cells). Rows represent CO₂ fluxes: gross primary productivity (GPP), ecosystem respiration (R_eco), and net ecosystem exchange (NEE). Negative values indicate uptake of CO₂ from the atmosphere by the vegetation, while positive values denote release. Shown is the average seasonal cycle and standard deviation (shading) in daily resolution, smoothed with a 10 d running mean. Goodness of fit is measured with root mean square error normalised by amplitude (nRMSE). Last row: annual sum of GPP simulated for all high-latitude PFTs in ORCHIDEE (green) and observed at the EC sites (brown). Note that where two EC sites are located in one ORCHIDEE grid cell, the observed fluxes shown represent different time frames and the ORCHIDEE fluxes are shown for the combined time frame.

3.1.2 Carbon fluxes

Comparing simulated seasonal cycles of CO₂ fluxes to eddy-covariance observations at six shrub-tundra sites shows that the evergreen dwarf shrub PFT (PFT 18) best matches the observed fluxes, the low deciduous shrub PFT (PFT 17) scores reasonably well, and the tall deciduous shrub PFT (PFT 16) overestimates fluxes and mismatches timing (Fig. 4). This is in line with the vegetation descriptions of the EC sites pointing to dominance of low and dwarf shrubs. Modelled biases from gross primary productivity (GPP) and ecosystem respiration (R_eco) propagate into simulated net ecosystem exchange (NEE), causing the largest mismatch in net carbon balance at some sites, while cancelling each other out at others. Despite their limitations to fully reproduce the observed seasonal cycle of CO₂ fluxes, the shrub PFTs capture total annual GPP significantly better than the baseline high-latitude tree and grass PFTs (Fig. 4, see annual GPP). The diurnal and seasonal EC data were not used in the sensitivity test or the optimisation, making this a true data-model evaluation.

In terms of plant productivity (GPP), both the low deciduous (PFT 17) and evergreen dwarf shrub (PFT 18) PFTs capture the observed timing and magnitude of the average annual cycle reasonably well at all three sites (nRMSE: 0.15–0.36 (PFT 18); 0.21–0.42 (PFT 17)), while the tall deciduous shrub PFT (PFT 16) delays growing season onset and overestimates peak growing-season uptake substantially (nRMSE 0.34–0.81) (Fig. 4). A consistent mismatch across the three shrub PFTs is the overly abrupt spring onset of GPP, not matching the slow onset suggested by the field observations. In particular, the deciduous shrub PFTs experience a rapid increase in photosynthetic activity at the onset of the growing season – a behavior that can be traced to model functionality where leaf onset for deciduous PFTs occurs within one single day. As a consequence of the rapid onset, simulated GPP of PFT 17 and 18 peaks slightly earlier in the growing season than observed. The magnitude of peak C uptake is matched realistically by PFT 18 at five out of six sites, while the low and especially the tall deciduous shrubs overestimate peak GPP.

For our specific modelling purposes in this paper, more urgent than matching timing and magnitude, is to correctly reproduce total annual carbon uptake. Annual sums are best reproduced by the low deciduous and dwarf evergreen shrub PFTs, whereas boreal trees, grasses and tall shrub PFTs tend to overestimate GPP, and underestimate it outside of their optimal growth range (Fig. 4, last row). Over the entire tundra region, median (IQR) annual GPP (1992–2020) is 587 (414–791) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for PFT 16, 453 (357–561) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for PFT 17 and 321 (233–403) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for PFT 18; the latter matching the 336 (209–413) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ estimated at the six EC flux sites best.

With respect to respiratory loses (R_eco), the dwarf evergreen shrub PFT shows the best fits across sites (nRMSE: 0.12–0.30), whereas both deciduous shrub PFTs overestimate growing-season respiration, reflected in higher nRMSE values (0.17–0.74 (PFT 17), 0.29–0.94 (PFT 16); Fig. 4). This can likely be attributed to overestimated maintenance respiration, while the simulated heterotrophic respiration outside of the growing season matches the observations well (cf. Fig. B2).

In ORCHIDEE, NEE reflects the net balance between GPP and R_eco, so biases in both gross fluxes propagate into NEE (Fig. 4). As a result of the underestimation of spring-time GPP, the shrub PFTs simulate weaker net C uptake than observed early in the year. The observed growing season NEE is matched well by the simulated evergreen dwarf shrubs in terms of magnitude, but the timing is shifted, with too early onset and end (nRMSE values of 0.19–0.38). For both deciduous shrub PFTs, the fit between simulated and observed NEE is less optimal (PFT 17: 0.22–0.45; PFT 16: 0.30–0.76). Both overestimate peak C uptake compared to four out of six site observations and misrepresent the seasonal pattern by simulating two peaks during the growing season. The first peak results from the strong onset of GPP, while the second arises from an abrupt decline in (maintenance) respiration at the end of the growing season. Due to the strong fluctuations and the underestimation of carbon uptake (GPP) outside of the growing season, the simulated total annual NEE is close to neutral for the shrub PFTs at all four grid cells (between -29 and 40 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), while the observed annual NEE shows a larger variation, from 197 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ at RU-Vrk (denoting net loss of CO₂ to the atmosphere) to -142 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ at RU-Cok (cf. Fig. B3).

3.2 Updated vegetation distribution

The updated PFT maps used in this study allocate on average 9.9 % (2.14 Mkm²) of the Arctic–Boreal study area > 50° N (∼ 21.75 Mkm²) to the three shrub PFTs during the 1992–2020 period. 1.26 Mkm² (∼ 5.8 % of the land area analysed) were allocated on average to evergreen dwarf shrub cover, 0.64 Mkm² (∼ 3 %) to low deciduous shrub cover, and 0.25 Mkm² (∼ 1.1 %) to tall deciduous shrubs. Introducing explicit shrub PFTs prevented the reassignment of shrub fractional cover to boreal tree and grass PFTs in ORCHIDEE PFT maps, thereby removing erroneously prescribed forest cover in tundra regions, and established the spatial baseline for subsequent changes in simulated carbon stocks and fluxes.

Figure 5 Change in fractional cover of high-latitude PFTs in ORCHIDEE after the introduction of three shrub PFTs. Bottom row: Fractional cover of the three new shrub PFTs. Maps are based on the vegetation distribution in 2010. Numbers indicate respective total PFT cover over the study area in 2010.

Tall deciduous shrubs are located sparsely across the Arctic–Boreal region but are most abundant across the southern latitudes, with hotspots in southwestern Alaska and western Siberia (Fig. 5). While low deciduous shrubs display a similar spatial pattern, evergreen dwarf shrubs dominate further north and are widespread across Arctic tundra regions, reaching up to 67 % fractional vegetation cover in individual grid cells. Boreal tree and grass cover classes have been reduced compared to the baseline ORCHIDEE-15 PFT simulation, and are no longer prescribed in place of shrub tundra (Fig. 5).

Figure 6 Simulated and observation-based pan-Arctic aboveground biomass. (a) Simulated and observation-based (from Spawn et al. (2020) and Liu et al., 2015) total pan-Arctic AGB in 2010. Difference maps between ORCHIDEE with 18 PFTs and the previous (ORCHIDEE-15PFT) model version, as well as the observation-based products. Hatching indicates agreement with Spawn et al. (2020) standard error interval. Numbers indicate total AGB across the region (top row), and total difference (bottom row). (b) Simulated (blue/green) and observed (brown) total pan-Arctic AGB 1992–2020.

3.3 High-latitude aboveground biomass

Across the Arctic–Boreal region, adding shrub PFTs with realistic biomass reduced the overestimation of AGB previously caused by prescribing boreal forest in shrub-dominated landscapes. Our simulations bring simulated carbon stocks into the range of observation-based estimates and improve the spatial agreement of biomass patterns compared to two independent AGB datasets (Fig. 6).

The combination of implemented shrub PFTs with a new PFT map has reduced the simulated aboveground biomass by 13.5 % from 54 Pg C in ORCHIDEE-15 to 46.7 Pg C in ORCHIDEE-18 over the entire high-latitude region (1992–2020), now falling into the range of 34.3 ± 16.1 Pg C estimated by Spawn et al. (2020) for 2010, and closely matching the 45.2 Pg C estimate by Liu et al. (2015) for 1993–2012 (cf. to 53.1 and 45.8 Pg C for ORCHIDEE-15 and -18, respectively, for 1993–2012). Simulated AGB consistently increases over the historical time period, but the trend is slightly higher than estimated by Liu et al. (2015) (Fig. 6b). Spatially, agreement with Spawn et al. (2020) improved in 2010 to 633 out of 989 land grid cells, corresponding to 56.9 % of the domain, falling within one standard error, compared to 479 out of 989 land grid cells (42.9 % of the area) for ORCHIDEE-15 (see Figs. 6a and B4). The biomass reductions are largest in western Siberia and southwestern Alaska (see Fig. 6a), identical to the two areas with largest replacement of boreal forest with shrub cover (cf. Fig. 5). Thus, the reduction in estimated AGB is stronger if considering the tundra region separately, with 34.8 %, from 4.6 Pg C in ORCHIDEE-15 to 3 Pg C in ORCHIDEE-18 (1992–2020), now falling between the estimates by Spawn et al. (2020) (2.8 ± 1.1 Pg C in 2010) and Liu et al. (2015) (4 Pg C 1993–2012). Overall, adding shrub PFTs to ORCHIDEE reduces the bias in tundra carbon stocks and significantly increases model agreement with benchmarking data.

Figure 7 Average annual GPP across the Arctic–Boreal region. (a) Distribution in ORCHIDEE-18 and different upscaled (Virkkala et al. (2021); FLUXCOM (Jung et al., 2019)) and data-constrained (CARDAMOM (Hugelius et al., 2024)) products, as well as difference maps with ORCHIDEE baseline and data-based products, in 2010. Numbers indicate area-weighted mean annual GPP across the region (top row), and total difference (bottom row). (b) Simulations (blue/green) and upscaled/data-constrained estimates (brown) of average annual GPP over the boreal and arctic region.(Virkkala et al. (2021): area-weighted mean of ensemble median of upscaled estimates (1992–2015); FLUXCOM: area-weighted mean and mean absolute derivation (2001–2015), CARDAMOM: area-weighted mean of model median and 95 % confidence interval (2001–2019)).

3.4 High-latitude carbon fluxes

Mean annual GPP across the Arctic–Boreal region (1992–2020) is reduced by 3.4 %, from 498 to 481 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ through the introduction of shrub PFTs, moving simulations toward the FLUXCOM estimate from upscaled flux observations (Jung et al. (2019): 350 ± 50.7 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (mean ± mean absolute derivation, 2001–2015); cf. to ORCHIDEE-18 (2001–2015): 490 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) and an estimate using the data-constrained data-assimilation framework CARDAMOM (Hugelius et al. (2024): 407 (388–426) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (weighted mean (95 % CI), 2001–2019); cf. to ORCHIDEE-18 (2001–2019): 491 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; Fig. 7b). Conversely, the upscaled Virkkala et al. (2021) estimate is closer to the higher GPP estimation from ORCHIDEE-15 (512 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Virkkala, mean 1992–2015); cf. to 495 and 478 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (ORCHIDEE-15 and ORCHIDEE-18 1992–2015, respectively)). However, this high region average masks a spatially heterogeneous estimate by Virkkala et al. (2021), supporting the higher ORCHIDEE-15 estimate over some parts of the domain, and the lower ORCHIDEE-18 estimate over others. Spatially, all three data-based products agree that ORCHIDEE overestimates GPP over most of the Eurasian side of the domain (cf. difference maps in Fig. 7a). On the North American side, the difference is more heterogeneous and varies between the data products. FLUXCOM and CARDAMOM point to an overestimation in the higher latitudes and an underestimation at the southern edge of the domain, while the Virkkala et al. (2021) estimate indicates that ORCHIDEE underestimates GPP across Canada and Alaska. Compared to the baseline ORCHIDEE version, GPP was reduced primarily in the tundra region, concentrated largely in the areas with high evergreen dwarf shrub cover (ORCHIDEE-18–ORCHIDEE-15 difference map in Fig. 7a; cf. Fig 5).

Figure 8 Simulated (blue/green) and upscaled/data-constrained (brown; Virkkala et al., 2021; Jung et al., 2019; Hugelius et al., 2024) average annual GPP, like in Fig. 7b, but only for the Arctic tundra region.

Consequently, a separate evaluation of the tundra region shows a stronger improvement, with a 13.5 % reduction in mean annual GPP from 334 (ORCHIDEE-15) to 289 (ORCHIDEE-18) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (1992–2020) and a closer yet imperfect agreement with FLUXCOM and CARDAMOM benchmarking data (FLUXCOM: 171 ± 45 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; CARDAMOM: 215 (201–231) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), while the Virkkala et al. (2021) estimate points to an underestimation (364 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$; Fig. 8).

Overall, including shrubs reduced total Arctic–Boreal GPP, even more so in the tundra region, and improved agreement with 2out of 3 data-based evaluation products.

Haverd et al. (2018)Reick et al. (2021)Harper et al. (2018)Meyer et al. (2021)Lawrence et al. (2019)Rezsöhazy et al. (2025)Wolf et al. (2008)Sulman et al. (2021)

Table 6 Overview of high-latitude shrub PFTs in selected land models of ESMs. Parentheses indicate shrub PFTs explicitly designed for outside of the high-latitude region. Information on shrub PFT height thresholds is indicated in parentheses where known.

4 Discussion

4.1 Shrub PFT classification and parameterisation

The new ORCHIDEE shrub PFT classification distinguishes shrub types by phenology (deciduous or evergreen) and plant height, following literature recommendations (Wullschleger et al., 2014; Mekonnen et al., 2021). Our classification captures the differing carbon cycle and ecosystem functions of the most prominent shrub types (Mekonnen et al., 2021) and lays the foundation for implementing distinct shrub-climate interactions and feedbacks. This categorisation balances real-world functional diversity with the computational constraints of a global LSM.

Our classification and parameterisation are more detailed than the shrub implementations in most other global LSMs, which often only differentiate between deciduous and evergreen shrubs and therefore are not able to capture growth-form-dependent differences and climate interactions (Bonfils et al., 2012; Wullschleger et al., 2014; Saccone et al., 2017; Mekonnen et al., 2021). A comparison of shrub implementations in different LSMs (Table 6; see Wullschleger et al. (2014) for further ecosystem models) shows that some models only include one generic shrub PFT in high-latitudes (e.g. CABLE2.4, CLM5, JSBACH3.2), while others at best differentiate between deciduous and evergreen (e.g. CLASSIC, JULES4.6). ELM and LPJ-GUESS, known for their high ecological realism (Martín Belda et al., 2022), include similar categorisations to the one presented in our study, differentiating between evergreen and deciduous shrubs and different heights (Wolf et al., 2008; Sulman et al., 2021). Our choice made for ORCHIDEE captures the same distinctions (evergreen-deciduous and height differences) and thus supports explicit representation of growth-form–dependent shrub–climate feedbacks (e.g., via phenology and canopy structure), while keeping the PFT space manageable given the current model complexity.

Contrary to the site-level calibrations in many other shrub implementations (e.g. Wolf et al., 2008; Meyer et al., 2021; Sulman et al., 2021), the ORCHIDEE shrub PFT classification and parameterisation are constrained by synthesised observations spanning the entire high-latitude region, increasing its transferability and limiting site-specific biases (Alton, 2011; Kuppel et al., 2014). This is in line with recent work by Murphy et al. (2025), highlighting the limitations of model parameterisations based on plot-level observations from a single site. Thus, the methodological choices, including the shrub PFT classification scheme, calibration targets and vegetation distribution, are applicable beyond the ORCHIDEE model.

While parameter selection and optimisation were mainly carried out in a single grid cell due to computational limitations, the results indicate that parameterisations can be considered regionally representative, since simulated shrub characteristics match synthesised observations not only in the optimisation grid cell, but across the entire region. Parameterisation in a single grid cell was further deemed sufficient for the targeted variables like height and biomass, because in ORCHIDEE they are less sensitive to local environmental conditions than the observed variability indicates (cf. the relatively narrow ranges of height, diameter and biomass over the entire region given in Table 5, and the smaller variability in AGB in ORCHIDEE compared to observations in Fig. 3). Notably, CO₂ flux-variables, with a stronger variability based on site conditions, were carefully evaluated at several locations to avoid site-specific biases. While this limited reproduction of spatial variability justifies the methodological choices in the present work, it calls for an enhanced representation of environmental dependencies in ORCHIDEE shrub parameter calculations in future work.

Our shrub implementation follows a classical fixed PFT approach, instead of a more flexible trait-based PFT implementation, where PFT traits vary in response to environmental conditions. Dynamic trait-based approaches have been proposed as a useful mechanism to capture functional diversity within generalised PFTs and dynamic responses to climate change (Wullschleger et al., 2014; Verheijen et al., 2013; Bodegom et al., 2011) and several studies suggest that understanding and implementing trait variations of high-latitude shrubs in models can improve projections of changes in ecosystem function and climate change feedbacks (Myers-Smith et al., 2019; Iturrate-Garcia et al., 2020; Eldridge et al., 2011; Yang et al., 2022). However, a trait-based implementation is challenging given the limited availability of shrub-plant trait observations (Myers-Smith et al., 2019; Thomas et al., 2020; Iturrate-Garcia et al., 2020). Our chosen shrub PFT implementation is sufficient to capture the role of different shrub types in present-day high-latitude carbon uptake and storage, as demonstrated by our results. However, future developments aiming to improve plant phenology timings and magnitudes, and capture shrub-climate interactions and future expansion may benefit from adding trait-based variations to the existing shrub PFTs. Recent developments of ORCHIDEE already introduced dynamic parameter calculations, of e.g. dynamic tree height and mortality (Peylin et al., 2026), that may be extended to shrubs. One such example could be for the tundra–treeline ecotone, where the structural distinction between tall deciduous shrubs and small trees can be fluid, particularly for taxa such as mountain birch (Betula pubescens ssp. tortuosa) and green alder (Alnus viridis) (Mauri and Caudullo, 2016; Óskarsdóttir et al., 2022). These may occur in both shrub-like and tree-like forms depending on climatic conditions. In such cases, dynamic parameterisation could, in principle, allow gradual transitions between shrub and tree growth forms.

Overall, the shrub PFT categorisation by phenology and height captures key functional contrasts and is sufficiently constrained to support regional carbon-cycle evaluation (Figs. 6–8), while leaving trait-continuous responses as a clear next development step.

4.2 Shrub implementation in ORCHIDEE

A key strength of this new shrub functionality lies in its relative simplicity: it introduces shrub PFTs by relying largely on existing processes already used for woody vegetation, rather than adding new functional modules to an already complex model. Our experiment demonstrates the adaptability of the ORCHIDEE framework.

In ORCHIDEE, all woody PFTs share core process formulations for photosynthesis, respiration, and carbon allocation to leaves, stems, and roots, while functional differences are expressed through PFT-specific parameters controlling plant height, leaf area index, allometry, and phenology (Krinner et al., 2005; Naudts et al., 2015). Shrub PFTs were implemented relying on this shared functionality, maintaining process consistency, and minimizing added model complexity, but adjusting parameter values that control key differences in ecosystem functioning between shrubs and trees. This highlights that the carbon allocation scheme of the recent ORCHIDEE version, based on the pipe model (Shinozaki et al., 1964), is generic enough to represent trees and shrubs with only a few parameter changes.

Key shrub–tree contrasts that are directly targeted in our calibration and evaluation include: morphological differences, which have been addressed through an adaptation of the growth form to low heights and stem diameter and higher density of stems; and adaptability of shrubs to limited growth conditions, which have been addressed through reduced sizes, as well as earlier leaf onset in PFT 17 and higher nitrogen use efficiency in PFT 18 (cf. Sect. 2.1.1).

Furthermore, shrubs differ from trees in their interactions with high-latitude environments. While the newly implemented shrub PFTs improve high-latitude vegetation representation and carbon cycling in ORCHIDEE, further developments are needed to explicitly capture shrub-climate interactions and feedbacks. Central feedback processes to consider when modelling shrubs include their interactions with snow, impacts on soil thermal regimes and permafrost dynamics, as well as albedo changes (Wullschleger et al., 2014; Mekonnen et al., 2021). For instance, snow cover influences shrub growth, by insulating and protecting shrub buds and tissues from frost, and shaping shrub height and growth form (Myers-Smith et al., 2011; Wheeler et al., 2014). Conversely, shrub presence increases snow accumulation and reduces its compaction, with impacts on soil temperatures and permafrost through the thermal insulation of the snowpack (Myers-Smith et al., 2011; Domine et al., 2016; Heijmans et al., 2022). Depending on shrub height and stature, shrub branches protruding above the snowpack can reduce winter albedo and accelerate snow melt (Belke-Brea et al., 2020; Shuman et al., 2025). These dynamics may potentially represent feedbacks on shrub growth and climate warming, and will thus be important to capture in modelling efforts of future high-latitude ecosystem-climate interactions in ORCHIDEE. After implementation of shrub-climate interactions, the next step would be the representation of dynamic shrub expansion processes through competition with other high-latitude PFTs.

Overall, our implementation demonstrates that shrub-appropriate PFT traits can be represented within ORCHIDEE’s existing woody framework, reducing biases in high-latitude carbon cycling (Figs. 6–8), while laying the groundwork for dynamic shrub-climate interactions.

4.3 Vegetation distribution

As vegetation fractions are prescribed in ORCHIDEE, the PFT maps condition simulated regional carbon stocks and fluxes. Updating these maps to include three shrub PFTs reduces the extent of forest-prescribed tundra and introduces spatially explicit shrub fractions that contribute to the biomass and GPP changes reported in Sect. 3.3 and 3.4.

Our updated PFT maps combine the ESA CCI–based PFT product (Harper et al., 2023) with additional shrub cover information (Raynolds et al., 2019) and apply empirical partitioning rules derived from local to sub-regional shrub mapping (Macander et al., 2017, 2022) to distribute shrub cover among the three ORCHIDEE shrub PFTs (Sect. 2.3). These additional steps address the known limitations of high-latitude shrub cover representation in Harper et al. (2023), and provide the information needed to distinguish between the three ORCHIDEE shrub PFTs. While this approach improves shrub representation in high latitudes compared to Harper et al. (2023) alone, it introduces its own set of uncertainties that should be considered when interpreting regional results.

First, partitioning shrub cover among the three shrub PFTs relies partly on relationships derived from Alaska/Yukon mapping products that are applied at pan-Arctic scale (Sect. 2.3). This extrapolation is necessary due to the lack of pan-Arctic mapping products with comparable detail in their shrub categorisation, but it may misrepresent regional differences in shrub composition. Second, shrub type definitions differ among observational and remote-sensing products and do not map uniquely onto the three ORCHIDEE shrub PFTs. For example, Macander et al. (2017) define low deciduous shrubs as 0.2–1.5 m tall, and tall shrubs as > 1.5 m, whereas our tall/low cut-off lies between ∼ 0.5 and 1 m; this mismatch may lead to an underestimation of tall-shrub fractions.

Independent evaluation of shrub cover fractions remains challenging because no pan-Arctic benchmark product provides comparable detail on shrub type and height. Using a reconstruction of tall shrub cover across 11 regions of the Siberian tundra (Frost and Epstein, 2014), we obtain a mean tall shrub cover of 13.9 %, supporting that tall shrubs may be underestimated in our PFT maps (mean 1.3 %, max 10.2 %). In the absence of unified benchmarks, we evaluate shrub distributions indirectly through their impact on regional AGB and GPP patterns, and interpret model improvements as the combined effect of shrub traits and shrub spatial distributions.

Overall, the updated PFT maps provide a necessary spatial foundation for representing high-latitude shrubs in ORCHIDEE and contribute to improved AGB and GPP relative to benchmarking data. However, uncertainty in shrub distributions remains a primary source of uncertainty in regional pan-Arctic carbon simulations.

4.4 Improved Model performance

Including shrub PFTs improves high-latitude carbon cycling in ORCHIDEE, reducing both simulated AGB and GPP across the Arctic–Boreal region and improving model performance compared to benchmark data.

Comparing simulated GPP with different data products highlighted a general agreement between simulated and observed GPP, but also underscored uncertainties in the upscaled observational datasets (Sect. 3.4). Across the Arctic–Boreal region, FLUXCOM (Jung et al., 2019) and CARDAMOM (Hugelius et al., 2024) provide lower GPP estimates than ORCHIDEE (350 (FLUXCOM, 2001–2015) and 407 (CARDAMOM, 2001–2019) $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, compared to 498 and 481 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, for ORCHIDEE-15 and ORCHIDEE-18, respectively, 1992–2020). In contrast, the Virkkala et al. (2021) upscaling yields estimates higher than ORCHIDEE-18 (512 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, 1992–2015; Fig. 7b). This spread among upscaled products implies that absolute GPP bias cannot be uniquely diagnosed from any single regional benchmark, and our simulated estimates fall within the range of recent regional data-constrained estimates derived from independent approaches.

Virkkala et al. (2021)Virkkala et al. (2021)Virkkala et al. (2021)Virkkala et al. (2021)Virkkala et al. (2021)Virkkala et al. (2021)

Table 7 Mean (standard deviation) annual ($\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) and seasonal ($\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{month}}^{-\mathrm{1}}$) GPP simulated in ORCHIDEE and observed/upscaled by different sources. (FLUXCOM (Jung et al., 2019), CARDAMOM (Hugelius et al., 2024), Virkkala et al. (2021) upscaled and observed, ABCflux database (Virkkala et al., 2022)). Seasons: spring = MAM, summer = JJA, autumn = SON, winter = DJF.

^a Virkkala et al. (2021) report growing season values (not JJA summer) with unknown length, in $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{period}}^{-\mathrm{1}}$, here divided by 3 assuming the growing season corresponds approximately to the three summer months (JJA). Values taken directly from the paper, reported over a larger area than examined in this study, resulting in an expected small overestimation. ^b Total annual values aggregated from reported seasonal means.

Differentiating between the tundra and boreal ecoregions (Dinerstein et al., 2017) and by seasonal fluxes, as well as adding synthesized site level observations to the comparison, provides a clearer picture of GPP biases (Table 7). In the boreal region, where the reduction in ORCHIDEE-18 relative to ORCHIDEE-15 is small (0.47 %), alignment between data-based products is high, indicating a modest positive bias in ORCHIDEE. A comparison with synthesized FLUXNET observations from Virkkala et al. (2021) further supports the overall high agreement between simulated GPP from both ORCHIDEE-versions and boreal zone site observations, with a slight improvement in ORCHIDEE-18. The FLUXNET synthesis in the ABCflux database suggests higher GPP values in the boreal region, mainly based on higher estimates outside of the growing season (Virkkala et al., 2022). The main disagreement among data-based products originates in the tundra region. There, the Virkkala et al. (2021) upscaled estimate and the synthesized flux observations of the ABCflux database point to higher GPP estimates than in ORCHIDEE-15, while the other products (FLUXCOM, CARDAMOM, Virkkala et al. (2021) synthesized observations) indicate an overestimation bias in both ORCHIDEE-15 and ORCHIDEE-18. Focusing on the growing season reveals a stronger and more consistent signal across all datasets: all data products agree that ORCHIDEE-15 overestimates GPP (100 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{month}}^{-\mathrm{1}}$) in tundra summer (JJA), and that the introduction of shrubs reduces this bias (85 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{mo}}^{-\mathrm{1}}$), increasing agreement with upscaled and synthesized data products (see Table 7). In contrast, outside of the growing season, ORCHIDEE underestimates GPP relative to the ABCflux database and FLUXCOM. The stronger reduction in simulated GPP in the tundra region (13.5 %), particularly during the summer, can be explained with a higher shrub fraction in the ORCHIDEE PFT maps in the tundra region, combined with the fact that differences in simulated GPP among ORCHIDEE PFTs are most pronounced during the growing season. Overall, multiple benchmarks indicate that shrub PFTs reduce GPP toward observation-based ranges, especially in the Arctic tundra subregion.

Hugelius et al. (2024)Hugelius et al. (2024)Ramage et al. (2024)Virkkala et al. (2021)Virkkala et al. (2021)Virkkala et al. (2021)Virkkala et al. (2021)

Table 8 Mean high-latitude NEE, total (Tg C yr⁻¹) and per area ($\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) simulated in ORCHIDEE and estimated by different sources (Hugelius et al., 2024; Ramage et al., 2024; Virkkala et al., 2022, 2021). Numbers in parentheses denote uncertainties, given as standard deviation (ORCHIDEE, ABCflux) or 95 % confidence interval (Hugelius et al., 2024; Ramage et al., 2024).

Reported Arctic–Boreal NEE estimates in the literature span a wide range (−587 (−862,-312) Tg C yr⁻¹ in the top-down estimate of (Hugelius et al., 2024) to 12 (−606,661) Tg C yr⁻¹ (Ramage et al., 2024)). Both ORCHIDEE setups simulate NEE at the higher end of this range (Table 8). In the boreal zone, both ORCHIDEE simulations are very similar and may overestimate NEE relative to observation-based estimates, whereas in the tundra, ORCHIDEE-18 simulates a reduced carbon uptake more consistent with observations. However, improvements in process representations governing R_eco and the onset of the growing season are needed to increase confidence in simulated NEE.

This study represents the first implementation and evaluation of shrub PFTs in the trunk version of ORCHIDEE (based on v4), which forms the basis for upcoming model applications, including CMIP-7 FAST-TRACK and TRENDY (Peylin et al., 2026). Accordingly, diagnostics directly constrained by the new plant representation were prioritised, with a particular focus on shrub dynamics during the growing season. We therefore emphasise GPP and C stocks (e.g. AGB), which are mechanistically linked to the processes and parameters targeted in this study (phenology, canopy development, allometry and allocation), rather than the full suite of processes governing ecosystem respiration and net carbon balance (e.g. soil carbon turnover, litter and heterotrophic respiration), which were not explicitly addressed during this implementation stage.

Ecosystem respiration is generally reported to be more uncertain than GPP, because it is more difficult to constrain observationally and depends on flux partitioning and multiple temperature- and substrate-sensitive processes (López-Blanco et al., 2018, 2019). Consequently, NEE, defined as the balance between GPP and R_eco, may appear unbiased due to compensating errors in photosynthesis and respiration. In our simulations, limitations in the seasonal dynamics of GPP – particularly an overly abrupt leaf onset for deciduous shrubs – propagate into R_eco and thus into NEE, complicating process-level attribution. For this reason, we report R_eco and NEE for completeness, but use GPP and biomass as the primary, more interpretable benchmarks for assessing whether the shrub PFTs improve high-latitude vegetation functioning at this stage of model development.

5 Conclusion

High-latitude vegetation has been under-represented in the global land surface model ORCHIDEE, contributing to biases in simulated carbon stocks and CO₂ fluxes across the Arctic–Boreal region. Relying largely on existing model functionality and synthesised pan-Arctic observations, we introduced three high-latitude shrub PFTs (tall deciduous, low deciduous, and evergreen dwarf shrubs) into ORCHIDEE LSM (r9269) and updated prescribed vegetation distributions to include shrubs.

Our shrub PFT classification combines two key dimensions (phenology and stature) to capture divergent shrub ecosystem functions while keeping the framework tractable for global and coupled simulations. Shrub traits are constrained using synthesised pan-Arctic observations to favour regionally representative parameter values rather than site-specific tuning. We prescribed shrub cover distributions using updated PFT maps that combine ESA CCI–based products with Arctic shrub information and sub-regional partitioning approaches.

Our implementation produces shrub PFTs with realistic size, biomass, and C allocation characteristics and improves simulated carbon stocks and gross primary productivity across the Arctic–Boreal region compared to regional benchmark datasets, with the strongest gains in tundra areas, where shrubs are more dominant and replace previously prescribed forest. This represents an important improvement of ORCHIDEE's representation of high-latitude vegetation and carbon cycling.

Important priorities for future work include improving seasonal dynamics of photosynthesis and respiration (and hence net ecosystem exchange), reducing uncertainties in shrub-cover mapping and shrub-type partitioning, and representing shrub-climate interactions involving snow, soil thermal regimes, and permafrost. Building on the baseline established here, implementing dynamic shrub expansion and trait variability will be essential for projecting future Arctic–Boreal ecosystem–climate feedbacks. More broadly, our results highlight the value of combining regionally synthesised observations with land surface model development to improve high-latitude vegetation representation in Earth system models.

Appendix A: Supplementary information for methods

A1 ORCHIDEE model description

ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) is the land surface component of the French ESM IPSL-CM (Boucher et al., 2020) developed by the Institute Pierre Simon Laplace (IPSL). ORCHIDEE can be run in a coupled set-up interacting with the IPSL atmospheric, ocean and ice sheet models, or as a stand-alone model, offline, as done in this work. Inputs required for the simulations include atmospheric forcing (air temperature, humidity, wind speed, surface pressure, and radiation), land surface characteristics (PFT distribution and dominant soil type), and initial conditions for energy, water, carbon, and nitrogen pools. These inputs are sourced from global datasets such as atmospheric reanalysis products, land cover maps, and soil texture databases. ORCHIDEE simulates carbon, nitrogen, energy and water cycling within the terrestrial biosphere, and the exchanges along the soil–vegetation–atmosphere interface (Krinner et al., 2005; Vuichard et al., 2019). Different biomes are represented by fifteen generic PFTs, including forests, crops and grasslands. They are differentiated based on climate zone (tropical/temperate/boreal), leaf habit (evergreen/deciduous), leaf type (broadleaf/needleleaf) and photosynthetic pathway for grasses and crops (C₃/C₄) (Box, 1996) (see Table A1). Woody vegetation is represented using a dynamic, vertically discretised canopy scheme that enables explicit simulation of tree-scale demographic and growth processes, including light penetration allowing photosynthesis to be calculated at each canopy level, sapling establishment and density-dependent self-thinning (Naudts et al., 2015). Three diameter classes within each woody PFT allow simulation of different plant sizes within a stand. Mortality due to fire, wind, bark beetles and forest management are implemented in ORCHIDEE, but deactivated in this work (Chen et al., 2018; Marie et al., 2024; Yue et al., 2014, 2018). The terrestrial carbon cycle, including processes of photosynthesis, respiration, soil carbon dynamics and plant carbon allocation and phenology, is simulated in the STOMATE module. Carbon assimilation and stomatal conductance are described following the concepts of Farquhar et al. (1980) and Ball et al. (1987), as implemented through the analytical solution proposed by Yin and Struik (2009). Photosynthetic processes are resolved at a half-hourly temporal resolution, whereas the allocation of carbon and nitrogen to vegetation compartments (namely leaves, fruits, roots, above-and belowground sapwood and heartwood, labile, and carbohydrate reserve pool), litter, and soil pools is simulated on a daily time step. The allocation of carbon to leaves, roots and wood follows allometric relationships, based on the pipe model theory (Shinozaki et al., 1964). Soil mineral properties are prescribed at the grid cell scale based on the 12 class USDA soil taxonomy (Forbes et al., 1987). The hydrological and biophysical processes of water and energy exchange are described by the surface–vegetation–atmosphere transfer scheme SECHIBA (de Rosnay and Polcher, 1998; Ducoudre et al., 1993). Energy fluxes in the model are calculated uniformly at the grid cell scale, without explicit representation of sub-grid heterogeneity, with the exception of snow processes, for which a distinct energy balance is evaluated over the snow-covered fraction of each cell. In contrast, hydrological processes are simulated separately for different vegetation components within a grid cell, employing three distinct soil columns representing bare soil, tall vegetation (e.g. trees), and short vegetation (e.g. grasses). Vertical discretisation differs among model components: soil hydrology is represented using an 11-layer profile extending to a depth of 2 m, with layer thickness increasing with depth according to a geometric spacing. Thermal processes and soil carbon and nitrogen dynamics share the same grid, but extended downward by seven additional layers to a total depth of approximately 90 m.

A2 ORCHIDEE PFTs

Table A1 ORCHIDEE plant functional types. PFTs 2 to 9 represent forest ecosystems. Bold: high-latitude PFTs; italics: shrub PFTs added in this study.

A3 Calibration data and targets

Height threshold values

Section 2.1.2 reports the final target heights for three ORCHIDEE shrub PFTs: a tall shrub PFT, including plants of 0.5 to 3 m height (typical shrub height being 1.5 m), a low shrub PFT, with a typical height of 0.5 m, and a dwarf shrub PFT covering heights below 30 cm (typical height: 20 cm). These are derived as a compromise among literature definitions, mapping products, and LPJ-GUESS. For transparency, the height thresholds in the sources considered are summarised below:

• Myers-Smith et al. (2015b): dwarf shrubs < 20 cm; erect dwarf/low shrubs up to 50 cm; tall shrubs up to 3 m.

• Elmendorf et al. (2012) (maximum potential height): dwarf shrubs < 15 cm; low shrubs 15–50 cm; tall shrubs > 50 cm.

• Circumpolar Arctic Vegetation Map (CAVM; Raynolds et al., 2019): prostrate/hemiprostrate dwarf shrubs < 15 cm; erect dwarf shrubs < 40 cm; low shrubs 40–200 cm.

• Alaska shrub mapping (Macander et al., 2017): dwarf shrubs up to 20 cm; low deciduous shrubs up to 1.5 m; tall deciduous shrubs > 1.5 m.

• Toolik Lake vegetation map (Walker and Maier, 2008): low shrubs 40–100 cm; tall shrubs > 2 m.

• LPJ-GUESS LSM (Wolf et al., 2008): dwarf shrubs up to 20 cm; low shrubs 20–50 cm; tall shrubs 50 cm–2 m.

Field survey used as an independent diameter check (Kobbefjord, Greenland)

Section 2.1.2 mentions a small survey of tall willow shrubs (Salix glauca) in Kobbefjord, Greenland, used to confirm that the diameter targets derived from Berner et al. (2015) are realistic. The measurements are provided in Table A2).

Table A2 Results of a small survey of tall willow shrubs in Kobbefjord, Greenland. Height, basal diameter (BD), main stem diameter at different heights (D_1.0, D_1.0), crown area (CA) and number of branches were measured/counted.

Table A3 Filter procedure applied to the vegetation-description fields in the Arctic Plant Aboveground Biomass Synthesis Dataset (Berner et al., 2024) to obtain PFT-specific shrub subsets. Counts refer to the shrub-dominated subset (≥ 80 % shrub AGB) before site-level aggregation. Biomass values shown are the site-level median across plots, summarised across sites (median and interquartile range) after conversion to g C m⁻² (50 % C content).

Aboveground biomass (Berner et al., 2024): exact filtering logic for shrub-PFT subsets

Section 2.1.2 describes the overall selection of shrub-dominated plots from the Arctic Plant Aboveground Biomass Synthesis Dataset (Berner et al., 2024), that compiles 32 individual datasets spanning 636 field sites across the Arctic region. The observations contain measurements of aboveground biomass of five different PFTs, including shrubs, mostly collected in peak summer season in 1998 to 2022, given in grams of oven-dried aboveground live biomass per square meter of ground surface. We extracted shrub biomass from the dataset by selecting the sample plots dominated by shrubs, using a selection criterion of at least 80 % of aboveground biomass being shrub biomass, well aware that this causes a slight underestimation of biomass compared to the model assumption that 100 % of a simulation site is covered by shrubs. We only considered shrub biomass from the resulting 677 plots at 331 field sites located within 24 larger areas. Following Berner et al. (2024), to avoid biases from sites with many plots, we calculated the median to group all sample plots at each field site, before calculating statistics across all field sites. To make the biomass data comparable to ORCHIDEE output, we converted the data to g C m⁻² assuming that half of the dry biomass is carbon. We used the vegetation descriptions contained in the dataset for each plot to extract specific aboveground biomass estimates for the three shrub PFTs. Table A3 provides additional implementation details on the string-based filtering of plot vegetation descriptions used to approximate shrub types. For tall deciduous shrubs we selected plots that contain the word tall or references to willow (salix) or alder (alnus) in their vegetation description, and excluded plots characterised as low or dwarf, resulting in 32 plots at 14 sites at 7 locations, with a median biomass of 583.5 (325.25–671.3125) g C m⁻², which is considerably higher than the overall shrub biomass. For low deciduous shrubs we filtered the vegetation description for a mention of low, not large, and Birch (betula or betula nana), finding 204 plots at 65 sites at 15 locations and a median biomass of 220 (130–479) g C m⁻². Finally, filtering for evergreen dwarf shrubs returned data from 410 plots at 238 sites at 15 locations, by filtering for the word dwarf and excluding plots with mentions of birch (betula), willow (salix), alder (alnus), tall, and other selected criteria, resulting in a median biomass of 131 (74–218) g C m⁻².

A4 Parameter optimisation

Optimisation was performed on site-level simulations at the Toolik Lake grid cell. The simulations were restarted from the end of a 340-years long spin-up simulation with prior best-guess parameter values, and then run with CRU JRA meteorological forcing from 2001–2010 for 110 years. The first 100 years served as additional spinup to allow the model to adjust to changed parameter values, and the final 10 years to calculate the cost function to minimise. If the cost function did not converge after 30 × 30 iterations, further iterations were added.

Table A4 Optimisation target ranges for the three shrub PFTs and the output variables: Height [m], peak summer aboveground biomass [g C m⁻²], minimum belowground biomass fraction, mean gross primary productivity [$\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$], peak growing season leaf area index, and age [years].

Optimised parameters included: pipe_tune2, alpha_self_thinning, recruitment_beta, recruitment_alpha, alloc_max, recruitment_height, k_latosa_min, k_latosa_frac, slainit. Different parameter combinations were optimised for the different shrub PFTs.

The following model output variables were optimised against observation-based targets: height, peak summer AGB, minimum fraction of belowground biomass, mean GPP, peak summer leaf area index, age. The optimisation targets are defined as ranges of variation for a set of model variables and are listed in Table A4. During the optimization the cost-function value contribution (see Eq. 4) from a given variable was assumed to be zero (i.e., considered to be acceptable) if the simulated variable time-series falls within the prescribed range. Whenever the simulated time series fell outside the prescribed range of acceptable variability, the cost function was defined as the deviation between the simulated values and the mean of the specified minimum and maximum bounds. Consequently, the contribution to the cost function increases with the magnitude of the deviation from the prescribed range.

Appendix B: Supplementary information for results

Figure B1 Simulated (green, ORCHIDEE-18, 1992–2020) and observed (light brown, Berner et al., 2024) aboveground biomass for three shrub PFTs across the Arctic and Boreal region. Observed biomass from Berner et al. (2024) was filtered to match the ORCHIDEE shrub PFT classification. Numbers in parentheses indicate the number of field sites included after filtering. Boxplots show median (center line), interquartile range (boxes) and whiskers extend to the farthest data point lying within 1.5× the interquartile range. Complete visualisation including outliers of Fig. 3.

Figure B2 Individual respiration fluxes of the ORCHIDEE shrub PFTs at the same grid cells as in Fig. 4 (columns). Rows represent respiration fluxes: heterotrophic respiration, maintenance respiration and growth respiration. Positive values indicate release of CO₂ from the vegetation to the atmosphere. Shown is the average seasonal cycle and standard deviation (shading) in daily resolution, smoothed with a 10 d running mean.

Figure B3 Simulated (green) and observed (brown) annual sum of CO₂ fluxes at four ORCHIDEE grid cells and six EC sites (as in Fig. 4). Rows represent CO₂ fluxes: gross primary productivity (GPP), ecosystem repiration (R_eco), and net ecosystem exchange (NEE). Negative values indicate uptake of CO₂ from the atmosphere by the vegetation, while positive values denote release.

Figure B4 Simulated (ORCHIDEE with 15 PFTs) and observation-based (from Spawn et al. (2020) and Liu et al., 2015) pan-Arctic AGB in 2010. Difference maps between ORCHIDEE with 15 PFTs and the observation-based products. Hatching indicates agreement with Spawn et al. (2020) standard error interval. Numbers indicate total AGB across the region (top row), and total difference (bottom row).