Quantifying the Impact of Climate Change and Forest Management on Swedish Forest Ecosystems Using the Dynamic Vegetation Model LPJ‐GUESS

Abstract

Boreal and temperate forests are undergoing structural, compositional and functional changes in response to increasing temperatures, changes in precipitation, and rising CO₂, but the extent of the changes in forests will also depend on current and future forest management. This study utilized the dynamic vegetation model LPJ‐GUESS enabled with forest management (version 4.1.2, rev11016) to simulate changes in forest ecosystem functioning and supply of ecosystem services in Sweden. We compared three alternative forest policy scenarios: Business As Usual, with no change in the proportion of forest types within landscapes; Adaptation and Resistance, with an increased area of mixed stands; and EU‐Policy, with a focus on conservation and reduced management intensity. LPJ‐GUESS was forced with climate data derived from an ensemble of three earth system models to study long‐term implications of a low (SSP1‐2.6), a high (SSP3‐7.0), and a very high (SSP5‐8.5) emissions scenario. Increases in net primary production varied between 4% and 8% in SSP1‐2.6, 21%–25% in SSP3‐7.0 and 25%–29% in SSP5‐8.5 across all three forest policy scenarios, when comparing 2081–2100 to 2001–2020. Increased net primary production was mediated by a higher soil nitrogen availability and increased water use efficiency in the higher emission scenarios SSP3‐7.0 and SSP5‐8.5. Soil carbon storage showed small but significant decreases in SSP3‐7.0 and in SSP5‐8.5. Our results highlight differences in the predisposition to storm damage among forest policy scenarios, which were most pronounced in southern Sweden, with increases of 61%–76% in Business‐As‐Usual, 4%–11% in Adaptation and Resistance, and decreases of 7%–12% in EU‐Policy when comparing 2081–2100 to 2001–2020.

Introduction

Climate change is influencing forest ecosystems in multiple ways, causing changes in productivity, shifts in the natural distributions of species, and alterations in the frequency and intensity of climate-related disturbances (Anderegg et al., 2022; Arneth et al., 2020; Jönsson & Lagergren, 2018; Seidl et al., 2017). This is leading to changes in the supply of ecosystem services (ES), in other words, the benefits that society can derive from ecosystems (Millennium Ecosystem Assessment, 2005). An increased understanding of how these material and immaterial natural values are affected by human-induced land-use changes and changes in climate is therefore urgently needed (Costanza et al., 1997; IPBES, 2019; Potschin & Haines-Young, 2011).

Boreal forests hold about one-third of the total global forest carbon (C) stock, and about two-thirds of the ecosystem C is stored below ground in soils (Pan et al., 2011). These forests have experienced a decline in the net C uptake from the 1990s to the 2010s due to changes in abiotic and biotic disturbance frequencies, gradual soil warming and altered harvesting regimes (Pan et al., 2024). In northern Europe, climate warming is expected to cause an increased occurrence of forest fires, pathogen and insect-related disturbances, and droughts (Ruosteenoja et al., 2018; Seidl et al., 2017). Although changes in the future wind speeds remain uncertain, storm damage may also be amplified at higher latitudes as an outcome of a reduction in the number of days with frozen soil during winters (Jönsson et al., 2015; Peltola et al., 1999). However, gradual climate warming in cold-limited regions may cause a prolonging of the vegetation season, and increase C assimilation and net primary production, potentially offsetting some of the decline in C uptake (Wang et al., 2023).

In Sweden, boreal and temperate forests cover 68% of the land surface (SNFI, 2024). The Swedish annual production of timber is among the highest within the EU, and in 2015, 94% of the annual net increment was harvested (Forest Europe, 2020). The demands on forests to provide renewable resources for a growing bio-economy are increasing in Sweden (SNFP, 2018). The Swedish National Forest Program, launched in 2018, promotes practices to increase productivity to further stimulate the sustainable production of timber, pulp and biofuels. However, the current approach to forest management in Sweden has been criticized for retaining a disproportionately strong orientation toward achieving production outcomes, despite legislation demanding equal prioritization of environmental concern and production (Lindahl et al., 2017). Intensive forest management strategies emphasizing timber production may reduce the forests' potential to provide other benefits, and provisioning ES are often traded off against regulating/supporting ES (Pohjanmies et al., 2017; Sing et al., 2018). The 2022 assessment of the national Environmental Quality Objective “Sustainable forests” revealed negative trends for biodiversity in Swedish forests, primarily due to a loss of old-growth forests with high biodiversity values, and fragmentation of habitat (SFA, 2022a). Competing demands for timber production have caused the area of previously unfelled pristine forest in Sweden to decrease between 2003 and 2019 (Ahlström et al., 2022). The number of different ES provided in forests in boreonemoral Scandinavia generally increases with stand age, which emphasizes the importance of protecting old-growth stands to maintain multifunctional forests (Jonsson et al., 2020). Solutions such as prolonging rotations and reducing harvesting intensities are also associated with multiple additional benefits, such as improved C storage and climate change mitigation potential (Law et al., 2018; Skytt et al., 2021). Despite the negative trends for old-growth forests, habitat quality in retention patches and set-asides in managed production forests in Sweden is improving, exemplified by the nationwide gradual increase in deadwood from 6.2 m³ ha⁻¹ in 2000 to 9.4 m³ ha⁻¹ in 2020 (SNFI, 2024). Further improvements to habitat size and quality may also result from the EU Law on Nature Restoration which was passed in 2024, due to the binding requirements to restore 20% of degraded land and sea areas in the EU by 2030 (European Environment Council, 2024).

The far-reaching consequences of climate change and of management decisions on forest structure and functioning are challenging to assess, especially when considering outcomes in the intermediate to far future (Belyazid & Giuliana, 2019; Pilli et al., 2022). Dynamic vegetation models (DVMs) are tools which simulate terrestrial vegetation dynamics as a function of climate input and environmental data (Littell et al., 2011). DVMs are suitable for addressing questions that concern intermediate to long-term changes in biogeochemical cycles, ecosystem structure and functioning, and their implications for ES and societal benefits at local to regional and global scales (Gregor et al., 2022; Islam et al., 2024; Olin et al., 2015; Zhu et al., 2016). DVMs also capture transient changes in vegetation over time, including interannual or seasonal variation in growth or in the exchange of C and water between the vegetation and the atmosphere (Bergkvist et al., 2023; Islam et al., 2024).

Here we use the DVM LPJ-GUESS (v4.1.2, rev11016) to assess the possible impact of climate change on Swedish forest ecosystems over the course of the 21st century and quantify the potential effects on a range of indicators of forest ES and ecosystem functioning (EF). The chosen model version incorporates updated functionality to simulate a range of forest management systems, including the clear-felling system, suitable for management of even-aged coniferous monocultures. The simulation of uneven-aged silvicultural systems includes continuous cover forestry with single-tree selection, characterized by a constant maintained canopy cover and regular removals of the largest trees (Lindeskog et al., 2021). LPJ-GUESS has previously been evaluated for its capacity to reproduce the structure of European production forests (Lindeskog et al., 2021) and the standing volume in southern, central, and northern Sweden (Bergkvist et al., 2023), and model results have shown agreement with observations.

This study aims to investigate the plausible effects of changing climate conditions and management practices on production forests in Sweden. We consider the impact of a low emissions scenario (SSP1-2.6), a high emissions scenario (SSP3-7.0), and a very high emissions scenario (SSP5-8.5). The consequences of forest management practices over the 21st century are assessed within three forest policy scenarios, representing (a) a continuation of contemporary management practices, (b) an increased focus on adaptation to reduce the harmful effects of climate change, and (c) shifting practices toward conservation and reduced management intensity. Specifically, our objectives are to: (a) analyze the trends in C uptake and emissions over the course of the 21st century, and to determine the impacts on (b) forest ecosystem function, (c) forest ecosystem services and on (c) the predisposition to storm damage, comparing the outcomes at the end of the century (2081–2100) to a contemporary historical period (2001–2020). We hypothesize that the high (SSP3-7.0) and very high (SSP5-8.5) emission scenarios will enhance both net primary production and heterotrophic respiration rates, with major connotations for multiple ecosystem functions and ecosystem service supply at the end of the century (Holmberg et al., 2019; Wang et al., 2023). Furthermore, we hypothesize that the predisposition to storm damage by the end of the century will increase the most in the scenario representing a continuation of contemporary management, due to the higher proportions of storm-sensitive conifers in this scenario.

We also perform a species-specific calibration of two parameters governing modeled tree height to diameter ratio, to address previous suggestions for improving modeled forest structure (Bergkvist et al., 2023). The calibration is presented in Appendix A with an evaluation of model results against observational data from the Swedish National Forest Inventory (SNFI) to validate modeled height, diameter, stand density, and standing volume for the main tree species in the study domain (SNFI, 2022b).

Material and Methods

Ecosystem Model Description

LPJ-GUESS is a dynamic vegetation model which simulates terrestrial vegetation based on external input of climate and environmental data (Smith et al., 2001, 2014). Vegetation structure and composition is modeled explicitly for plant functional types (PFTs) that either represent a functional group of species or a specific tree, shrub or herbaceous species (Hickler et al., 2012). The model has been applied in a wide variety of studies on ecosystem functioning, vegetation dynamics and on C and nitrogen (N) exchange on both global, regional, and local scales (Ahlström et al., 2012; Gustafson et al., 2021; Jönsson & Lagergren, 2018; Wårlind et al., 2014). In this study we utilize the forestry-enabled version of LPJ-GUESS (Lindeskog et al., 2021) in which forest management is explicitly expressed. Moreover, we use the European version of the model in which PFTs have been parameterized as tree species (Hickler et al., 2012; Lindeskog et al., 2021). We simulate the coniferous species Scots pine (Pinus sylvestris L.), Norway spruce (Picea abies L. Karst.) as well as the broadleaved birch (Betula pendula Roth.) which together represent almost 92% of all standing stock on productive forest land (defined as forest land with a mean annual increment of 1 m³ ha⁻¹ or higher) (SNFI, 2024). Additionally, we also simulate Pedunculate oak (Quercus robur L.) and European beech (Fagus sylvatica L.), two broadleaved species associated with high recreational, esthetic and biodiversity values in southern Sweden.

Vegetation in LPJ-GUESS is simulated dynamically in replicate patches that typically represent an area of 1000 m², where individual cohorts of the selected PFTs compete for resources. External, spatially explicit climate data required for running the model include precipitation, surface downwelling short-wave radiation, air temperature, atmospheric C dioxide concentration [CO₂], N deposition, and fixed soil texture input data. Plant physiological processes such as stomatal regulation, photosynthesis, respiration and phenological development are represented mechanistically (Sitch et al., 2003; Smith et al., 2014). Modeled photosynthesis is dependent on the air temperature, surface downwelling shortwave radiation, soil water and N content, and on CO₂ concentration. Stomatal conductance may be reduced during periods of water stress, which inhibits leaf CO₂ uptake, exerting a limiting influence on photosynthesis. This mechanistic linkage also allows for an increased uptake of CO₂ per unit of water transpired in a higher CO₂ environment. For summergreen PFTs, leaf flushing occurs when a growing degree-day threshold has been passed, which is defined at the level of individual PFTs. Leaf shedding occurs 210 days after full canopy closure (Smith et al., 2014). Height and diameter growth of individuals result from annual allocation of net primary production (NPP) to roots, sapwood, and leaves and are constrained by PFT-specific allometric relationships. The growth of vegetation is influenced by competition among individuals for resources of sunlight, nutrients and water (Smith et al., 2014). Mortality results from self-thinning, low growth efficiency, environmental stress or from age. Mortality may also occur from patch-destroying disturbances with a return time commonly set to 100–400 years. In this study, these were enabled in all grid cells during model spin-up to gradually build up stable soil C and N pools. They were also enabled in protected areas and set-asides during the historical and transient simulations, but not in managed forests (see 2.3). When natural vegetation is simulated, establishment is determined by PFT-specific bioclimatic limits, the light level at the forest floor and PFT-specific maximum establishment densities determined by shade tolerance traits (Sitch et al., 2003; Smith et al., 2001).

Decomposition of soil organic matter (SOM) and the resulting mineralization of N is dynamically simulated within 11 SOM pools with varying C to N ratios and rates of decay, based on the approach utilized in the CENTURY model (Parton et al., 1993; Smith et al., 2014). Decomposition of SOM results in heterotrophic respiration (Rh) to the atmosphere. The Rh rate and the N mineralization rate depends on the decay rate of SOM and the C:N ratio of each pool, on simulated near-surface (25 cm depth) soil moisture and on near-surface (25 cm depth) soil temperature. The amount of available N in the soil depends on the rate of mineralization, the demand of the receiver SOM pool, and the plant uptake of N. If the plant demand for N outweighs the soil supply, then growth is reduced and more C is allocated to fine roots the following year in order to facilitate increased N uptake (Smith et al., 2014). Soil biological N fixation is modeled mechanistically whereas atmospheric N deposition is determined from external input data. Leaching of N varies with the soil sand fraction and with the percolation of water from the soil, but N losses may also occur as a result of gaseous emissions and during wildfire disturbances (Smith et al., 2014). LPJ-GUESS requires a spin-up period to generate stable vegetation compositions (commonly 500 years) and C and N pools (maximum 50,000 years), where detrended climate data for the first 30 years of climate data (1850–1879 in the present application, see below) is repeated with constant CO₂ concentration and N deposition rates that correspond to pre-industrial values.

The added functionality of forest management (FM) incorporates even-aged FM systems, where trees are established after clear-felling either by planting or through natural regeneration at the patch level. The FM module also includes uneven-aged management systems such as single-tree selection based on a set target diameter as well as shelterwood systems (Lagergren & Jönsson, 2017; Lindeskog et al., 2021). A user-defined management scheme specifies the tree species (PFTs) composition at the initiation of a simulation at the scale of individual stands. Clear-felling and thinning events cause the extraction of biomass from a simulated patch. These may be set to occur either at specified simulation years, or may alternatively be triggered by stand density and diameter size thresholds. Thinning includes the removal of biomass from below or from above, based on a user-defined setting. A module determining forest predisposition to storm damage, presented as an annual index value for simulated stands, has also been developed (Lagergren et al., 2012).

Climate, Deposition and Soil Data

We used climate data simulated in the Coupled Model Intercomparison Project Phase 6 (CMIP6) at a spatial resolution of 0.5ᵒ ×  0.5ᵒ during all transient simulations as model input (Figure 1 and Table B1) (Eyring et al., 2016). The climate forcing data were provided as output from the three earth system models MRI-ESM2.0 (Yukimoto et al., 2019), EC-Earth3-Veg (Döscher et al., 2021) and GFDL-ESM4 (Dunne et al., 2020). The data were supplied by the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) phase 3b (Lange & Büchner, 2021) and have been bias-corrected according to standardized procedures for ISIMIP climate data (Lange, 2019). Daily totals of precipitation, surface downwelling shortwave radiation and mean air temperature spanning the years 1850–2100 were used as input to LPJ-GUESS. The input data included forcing for three alternative future climate change trajectories SSP1-2.6, SSP3-7.0 and SSP5-8.5 that continue from the historical period ending in 2014 and run until year 2100. Annual atmospheric CO₂ concentrations for 1850–2100 were taken from CMIP6 (Meinshausen et al., 2020) and modeled annual N deposition rates for SSP1-2.6, SSP3-7.0 and SSP5-8.5 during 1850–2099 were provided by the National Center for Atmospheric Research (NCAR) (Hegglin et al., 2018). Soil input data for texture were based on the WISE 3.0 (Wide-field Infrared Survey Explorer, version 3) data set which represents mineral forest soils, in other words, histosols were excluded from the study (Batjes, 2005).

[IMAGE OMITTED. SEE PDF]

Setup of Forest Management During Transient Simulations

LPJ-GUESS enabled with forest management was used to simulate natural and managed forests in all grid cells covering productive forest land in Sweden (mountainous regions were excluded) (Figure 2). Norway spruce, Scots pine and Silver birch were simulated within all regions. The two nemoral tree species beech and Pedunculate oak were simulated in southern Sweden. Due to the relatively limited silvicultural interest in these latter two species in Sweden, our simulations assume that they are not established in managed forests in central Sweden in the future even if changing climate conditions may allow it.

[IMAGE OMITTED. SEE PDF]

Protected and set aside areas were simulated as natural unmanaged forests with a mix of different naturally occurring tree species. Random stand-replacing disturbances with a return time of 100 years were enabled in simulations of unmanaged forests to represent natural forest dynamics. As the occurrences or magnitudes of these disturbances are not influenced by the resistance or resilience of simulated forest types or altered by any specific climate conditions, their generic representation was considered unsuitable to visualize future potential damage in the range of different forest types included. For this reason, they were disabled in managed production forests. Instead, we here represent an increasing risk of disturbances by quantifying the predisposition to storm damage (see 2.5).

Monocultural stands were initiated with management regimes representing current recommended practices, with one pre-commercial thinning followed by one to three thinnings, depending on forest type and region during the nineteenth and twentieth centuries (Table 1). In forest types managed with clear-felling the thinning regimes prioritized the removal of small diameter trees (thinning from below). The first thinning in Scots pine monocultures and in Norway spruce monocultures in each region was set to occur at the approximate age when the stands reached heights between 10–13 m and 12–15 m, respectively, following recommendations (SFA, 2015).

Table 1 Applied Forest Management Within Each Region of Sweden for the Simulated Forest Types

Forest type	Region	Planting density (stems ha−1)	Pre-comm. Thin.	1st thin.	2nd thin.	3rd thin.	4th thin.	Final felling (years)
Scots pine monoculture	Southern	2700	12 (15%)	45 (30%)	60 (20%)	-	-	75
Central	2500	19 (15%)	57 (30%)	76 (20%)	-	-	95
Northern	2300	21 (15%)	63 (30%)	84 (20%)	-	-	105
Norway spruce monoculture	Southern	2700	7 (20%)	42 (25%)	56 (15%)	-	-	70
Central	2500	9 (20%)	45 (25%)	72 (15%)	-	-	90
Northern	2300	15 (20%)	50 (25%)	75 (15%)	-	-	100
Silver birch monoculture	Southern	5000	9 (35%)	27 (25%)	36 (20%)	-	-	45
Central	4500	11 (35%)	33 (25%)	44 (20%)	-	-	55
Northern	4000	13 (35%)	39 (25%)	52 (20%)	-	-	65
Mixed coniferous forest (spruce-pine)	Southern	1400/1400	6 (20%)	20 (variable)	40 (variable)	60 (variable)	-	75
Central	1300/1300	9 (20%)	20 (variable)	40 (variable)	60 (variable)	80 (variable)	90
Northern	1200/1200	15 (20%)	20 (variable)	40 (variable)	60 (variable)	80 (variable)	100
Mixed forest (spruce-birch)	Southern	2300/4000	9 (20%)	20 (variable)	40 (variable)	-	-	60
Central	2000/3500	12 (20%)	20 (variable)	40 (variable)	60 (variable)	-	75
Northern	1700/3000	-	20 (variable)	40 (variable)	60 (variable)	80 (variable)	90
Oak monoculture	Southern	5000	12 (15%)	Thinnings every 10th year throughout the rotation (variable)	120
Beech monoculture	Southern	5000	10 (45%)	20 (45%)	40 (25%)	60 (25%)	80 (20%)	100
Shelterwood (Scots pine)	Southern	5000	10 (20%)	40 (25%)	90 (92%)	-	-	100
Central	4500	11 (20%)	44 (25%)	99 (92%)	-	-	110
Northern	4000	15 (20%)	48 (25%)	108 (92%)	-	-	120
Single tree selection (Norway spruce)	Southern	2700		Continuous removals every 10th year (10%)	-
Central	2500		Continuous removals every 15th year (10%)	-
Northern	2300		Continuous removals every 20th year (10%)	-
Unmanaged (set-aside)	All regions	-	-	-	-	-	-

Mixed stands of Norway spruce-Silver birch (hereafter spruce-birch) and Norway spruce-Scots pine (hereafter spruce-pine) were simulated with a selective cutting function which maintained an even proportion of woody biomass between the species in the mixture through reoccurring thinning extractions every twentieth year over the course of a rotation (Lindeskog et al., 2021). In southern Sweden, both Pedunculate oak and beech were simulated as even-aged stands (Table 1). Fully layered Norway spruce stands were simulated with single tree selection management, with harvest of the largest trees in the stand occurring every tenth to twentieth year depending on region, following established recommendations (SFA, 2014) (Table 1). Additionally, a shelterwood system was simulated as an alternative for Scots pine monocultures.

Each simulated stand was replicated within 6 patches in each grid cell. The model was run in cohort mode, which implies that individuals within each patch retained the same size for a given age. The stands were initiated through planting or natural regeneration at different years during the nineteenth and twentieth centuries, so that at each simulated time step, the stands of even-aged forest types differed in age by 5 years. For example, Norway spruce monocultures in southern Sweden with a rotation period length of 70 years were represented by 14 stands of different ages, ranging from 1 to 70 years (Table 1). The final felling age gradually increased from southern to northern Sweden for all forest types, hence the area of younger forests was higher in the more productive southern region compared to in central or northern Sweden. The simulations assumed that the rotation lengths were not altered over time.

Forest Policy Scenarios

All forest types were simulated one-at-a-time in the model for the entire 21st century (2001–2100). The forest types were used to create three different forest policy scenarios, representing alternative compositions of the forest landscape. Output was generated for each forest type in each grid cell. The outputs from each simulated forest type were weighted based on the assigned proportions of forest area within each forest policy scenario (Figure 3). This resulted in differing relative contributions between forest types to the overall mean for each simulated output variable in each policy scenario. The weighting approach assumed that all forest types exist in every grid cell in the proportions that are given for the region that the grid cell is situated in. Hence, we did not take local variations in the cover of forest types into account within specific areas in a region.

[IMAGE OMITTED. SEE PDF]

The first forest policy represents the contemporary forest composition in Sweden, denoted as Business As Usual (BAU). In order to form the BAU landscape, the simulated forest types (Table 1) were weighted based on the areal coverage of each forest type in each of the three regions in Sweden in 2020 (SNFI, 2022a). BAU assumes that the area of each forest type within each region stays constant over the entire simulation period (2001–2100) (Figure 3). In BAU, the proportion of the forest landscape subjected to clear-felling varies from 84% (northern Sweden) to 86% (southern Sweden). Data on the area of protected forest land were taken from the Swedish Forest Agency (SFA), and aggregated from the county to the regional level for southern, central and northern Sweden, which determined the proportion of set-asides and protected areas in each region (SFA, 2022b).

Additionally, two alternative forest landscapes were created for the far future (2081–2100) by assigning different weights to each of the simulated forest types (Figure 3). The transition from the contemporary forest landscape in 2020 to a future where the landscape composition has changed requires assumptions on the speed of transition and the successful transformation of stands. Some stands require a long time to achieve a successful transformation, for example, when converting Norway spruce monocultures into a stand managed with selection felling (Lundqvist, 2017). In the analysis of these two alternative policy scenarios we decided to focus specifically on the far future (2081–2100), as it allows for the assumption that the transformation of those forest types would be completely successful during that time, which might not be the case by mid-century. The scenarios for the far future intend to highlight potential implications of changes in forest policy and to form a basis for comparison to BAU.

The Adaptation and Resistance (AR) scenario was created to depict a development where risk-spreading strategies become more prevalent in response to climate change-related hazards. Establishing mixed forests is one of the most frequently undertaken measures of adapting to climate change among private forest owners in Sweden (Blennow, 2012). In AR, the proportion of even-aged monocultural stands decreases in the landscape at the end of the 21st century and mixtures of spruce-birch and pine-spruce increases (Figure 3). In AR, the proportion of production forests managed with clear-felling ranges from 80% (northern Sweden) to 87% (southern Sweden).

The EU-Policy (EUPOL) forest scenario represents a significant shift in the composition of the forest landscape compared to BAU, where the emphasis on timber production is reduced to favor alternative values such as habitat protection and cultural ecosystem services (Figure 3). It reflects an alternative future where the Swedish forest policy aligns to a larger extent with the EU Biodiversity Strategy for 2030 (European Commission, 2020) and with the EU Nature Restoration Law (European Environment Council, 2024). The Biodiversity Strategy aims to increase habitat availability through restoration and through improvement of habitat connectivity within forests, with specific aims to improve the quality and quantity of dead wood. In EUPOL, we assumed an increase in the area of protected and set-aside forest land from 9% to 18% in southern Sweden, from 11% to 22% in central Sweden, and from 12% to 24% in northern Sweden (Figure 3). The forest area managed with continuous cover forestry (CCF) is also approximately doubled within each region compared to in BAU. In EUPOL, between 60% (northern Sweden) and 74% (southern Sweden) of the landscape is subject to clear-felling.

Indicators for Ecosystem Services and Ecosystem Functioning

The model output of LPJ-GUESS was used to calculate indicators of seven ES and six EFs, as well as an indicator of the forest predisposition to storm damage for each forest policy scenario. In each forest policy the indicators were weighted based on the proportion of different forest types (Figure 3). The ES assessed in this study via the indicators include both regulating, maintenance and provisioning services in forests (Table 2). These were determined based on the CICES classification of ecosystem services (Haines-Young & Potschin, 2018).

Table 2 Indicators for Ecosystem Services (ES), Ecosystem Functioning (EF) and for the Predisposition to Storm Damage

	Ecosystem service	Ecosystem service description	Section	Indicator name	Indicator description	Unit	Abbreviation
ES	Decomposition	Nutrient release and maintenance of soil quality	Regulation and Maintenance	N mineralization	Annual rate of net N mineralization in soil	kg N ha−1 year−1	N min
ES	Maintenance of habitat	Maintenance of ecological conditions required for sustaining forest-dwelling species	Regulation and Maintenance	Coarse woody debris	Annual addition of coarse woody debris to forest floor	kg C m−2 year−1	CWD
ES	Regulation of freshwater quality	Maintenance of the chemical quality of freshwaters	Regulation and Maintenance	Leached N	The net change in the annual amount of leached N from soil	kg N ha−1 year−1	-
ES	C sequestration in biomass and soil	The reduction of atmospheric CO2 by uptake in ecosystems	Regulation and Maintenance	Net Ecosystem Production	Annual total ecosystem C uptake	kg C m−2 year−1	NEP
ES	C storage in soil	The reduction of atmospheric CO2 by storage in ecosystems	Regulation and Maintenance	Soil C	C stored in soil	kg C m−2	-
ES	C storage in biomass	The reduction of atmospheric CO2 by storage in ecosystems	Regulation and Maintenance	Biomass C	C stored in roots, leaves and wood of trees	kg C m−2	-
ES	Supply of timber, pulp and biofuels	Harvestable biomass from mature forest stands	Provisioning	Potential harvest C	Amount of C available for harvest for timber, pulp and biofuel production	kg C m−2	-
EF	-	-	-	Annual net primary production	Annual net primary production	kg C m−2 year−1	NPP
EF	-	-	-	Annual heterotrophic respiration	Annual heterotrophic respiration	kg C m−2 year−1	Rh
EF	-	-	-	Water Use Efficiency	Monthly gross primary production/monthly transpiration (Seasonal: May to Sept)	g C kg−1 H2O	WUE
EF	-	-	-	C Use Efficiency	Annual net primary production/Annual gross primary production	-	CUE
EF	-	-	-	Biomass C to soil C ratio	Biomass C/soil C	-	-
EF	-	-	-	Rh to litter C input ratio	Annual heterotrophic respiration/Annual total C litter input	-	-
-	-	-	-	Predisposition to storm damage	Storm sensitivity index	m3 ha−1	Storm SI

The average change in an indicator was used to highlight potential changes in ES, EF or in predisposition to storm damage, here defined as the difference between its historical value (average value for 2001–2020) and its projected value at the end of the century (average value for 2081–2100).

Six indicators quantified changes in ES related to the regulation and maintenance of forest ecosystems (Table 2). Annual net N mineralization was used as a proxy to describe the ecosystem service of decomposition, and the average annual rate of added coarse woody debris to forests was used as an indicator for the maintenance of habitat. Freshwater quality status was indicated by modeled annual rate of N leaching. Model output on stocks of soil C and tree biomass C were analyzed to determine changes in ecosystem C storage.

The C sequestration capacity of forest ecosystems was modeled as the annual average net ecosystem productivity (NEP), estimated as: 1 $$\text{NEP}=\text{NPP}-\text{Rh}$$

Hence, a positive NEP indicates a net ecosystem C sink. Additionally, NPP and Rh were used to highlight changes in EF with regard to C cycling (Table 2). Seasonal mean water use efficiency was calculated for May to September, based on average values for each month for simulated gross primary production (GPP) and transpiration rate of vegetation (T) as: 2 $$\text{WUE}=\frac{\text{GPP}}{T}$$

One indicator provided information on the annual harvest of stem wood C, which represents the supply of ES relating to timber, pulp and biofuel production. For monocultures, the change in harvested stem wood C was calculated as the difference between stands felled in 2016 and 2096 (years during which all forest types were subjected to harvest). For stands of Norway spruce managed as CCF, changes in harvest were calculated as differences between the sum of all removals during the course of a standard rotation period for Norway spruce in the respective region (see Table 1). An indicator of storm sensitivity developed for Swedish forest conditions was used in order to estimate the predisposition to storm damage. The storm sensitivity index shows the potential storm damaged fraction of a forest stand that would result in case of a severe storm (Lagergren et al., 2012). The index is estimated as a product of several predisposing factors such as species, height to diameter relationship, differences in height within and between patches and time since thinning (Jönsson et al., 2015; Lagergren et al., 2012).

Differences between forest policy scenarios and emission scenarios were tested for statistical significance using a one-way analysis of variance (ANOVA) followed by Tukey's post-hoc test. If the assumptions of the model were not fulfilled (normality of residuals and homogeneity of variances), the non-parametric Kruskal-Wallis test and Dunn's Multiple Comparison test was used instead (Dunn, 1964; Kruskal & Wallis, 1952). Additionally, principal component analysis (Jolliffe & Cadima, 2016) was used to visualize the influence of using different ESMs on the model output.

Results

Trends in C Uptake and Emissions From 2011 to 2100

The long-term trends in LPJ-GUESS simulations representing BAU in combination with SSP3-7.0 and SSP5-8.5 indicated deviations from historical rates of C uptake and emissions, with gradually increasing rates of NPP and Rh in all regions (Figure 4), suggesting that changes of the largest magnitude will manifest at the end of the century in these high emission scenarios. The trends were less pronounced for NPP in the SSP1-2.6 simulation, with a peak in all three regions during mid-century. NPP was 9%, 12% and 17% higher in southern, central and northern Sweden, respectively, by 2060 compared to 2001 (Figure 4). The levels of NPP and Rh in AR and in EUPOL were notably different from BAU during 2081–2100. The deviation in NPP, comparing BAU to AR or EUPOL during 2081–2100 was generally larger in central and in northern Sweden than in southern Sweden (Figure 4). The mean Rh was consistently higher for EUPOL than for BAU and AR in all regions regardless of emission scenario (Figure 4). The average regional NEP remained positive over the studied period for both SSP1-2.6, SSP3-7.0 and SSP5-8.5. In general, the outcome of changing the forest scenario on NEP were minor, but with a trend toward a lower C uptake in EUPOL and a higher C uptake in AR (Figure 4).

[IMAGE OMITTED. SEE PDF]

End-Of-Century Changes in Ecosystem Functioning

In SSP1-2.6, annual average NPP increased significantly in AR (6%) and in EUPOL (8%) (p < 0.001) but not in BAU (4%) when 2081–2100 was compared to 2001–2020 (Figure 5). The increase in annual NPP ranged from 21% in BAU, to 25% in EUPOL (p < 0.001) in SSP3-7.0, and from 25% in BAU to 29% in EUPOL in SSP5-8.5 (Table 3). The trends for NPP varied with latitude, and gains were most pronounced in northern Sweden in all three emission scenarios (Figure 5). The end-of-century rates of Rh in SSP1-2.6 increased significantly in AR (5%) and in EUPOL (8%) (p < 0.01) but not in BAU (3%) (Table 3). Rh gains ranged from 20% in BAU to 26% in EUPOL in SSP3-7.0 (p < 0.001), and from 25% in BAU to 31% in EUPOL in SSP5-8.5 (p < 0.001). The substantially higher atmospheric CO₂ content in SSP3-7.0 and in SSP5-8.5 produced a fertilizing effect on photosynthesis, which reduced the trade-off between C uptake and transpiration, and enhanced WUE by between 30% and 34% in SSP3-7.0 (p < 0.001), and by 33%–39% in SSP5-8.5 (p < 0.001). The WUE increase was more modest in SSP1-2.6, ranging from 7.5% in EUPOL to 10% in AR (p < 0.001). In all three emission scenarios the biomass C to soil C ratio showed the largest gains in EUPOL due to the expansion of set-asides, protected areas and of stands managed with CCF, but the ratio also consistently increased in AR and BAU.

[IMAGE OMITTED. SEE PDF]

Table 3 Average Values for Six Indicators of Ecosystem Function, Seven Indicators of Ecosystem Services, and One Indicator of Predisposition to Storm Damage in Sweden

Indicator	Historical	SSP1-2.6	SSP3-7.0	SSP5-8.5
BAU	BAU	AR	EUPOL	BAU	AR	EUPOL	BAU	AR	EUPOL
NPP (kg C m−2 year−1)	0.51	0.53	0.55	0.55	0.62	0.64	0.64	0.64	0.65	0.66
Rh (kg C m−2 year−1)	0.43	0.44	0.45	0.46	0.51	0.52	0.54	0.53	0.54	0.56
NEP (kg C m−2 year−1)	0.09	0.09	0.10	0.09	0.11	0.11	0.11	0.10	0.11	0.10
C use efficiency	0.45	0.44	0.44	0.45	0.40	0.40	0.41	0.39	0.39	0.40
Biomass C (kg C m−2)	5.56	6.73	6.80	7.32	7.47	7.54	8.07	7.62	7.68	8.21
Potential harvest C (kg C m−2)	4.43	5.67	5.67	5.00	6.04	6.05	5.32	6.13	6.14	5.39
Soil C (kg C m−2)	12.2	11.9	11.9	11.8	11.8	11.8	11.7	11.8	11.7	11.6
Biomass C to soil C ratio	0.51	0.64	0.65	0.71	0.72	0.73	0.79	0.74	0.75	0.81
Rh to litter C input ratio	1.48	1.46	1.45	1.43	1.47	1.47	1.45	1.49	1.49	1.47
Coarse woody debris (kg C m−2 year−1)	0.07	0.08	0.07	0.08	0.09	0.09	0.10	0.10	0.09	0.10
Leached N (kg N ha−1 year−1)	6.26	5.21	5.09	5.15	6.79	6.63	6.66	7.41	7.17	7.21
N mineralization (kg N ha−1 year−1)	17.9	19.6	22.2	25.8	24.4	27.9	32.4	26.3	29.7	34.5
Water Use Efficiency (g C kg−1 H2O)	8.0	8.7	8.8	8.6	10.7	10.7	10.3	11.1	11.0	10.6
Predisposition to storm damage (m−3 ha−1)	27.3	40.2	34.3	32.9	40.8	34.3	33.5	42.0	34.9	33.8

Changes in Ecosystem Services

The mean change in NEP, comparing end-of-century rates to historical rates, was positive and significant for AR and EUPOL in all three emission scenarios (p < 0.05), but not for BAU in SSP1-2.6 when considering Sweden as a whole (Table 3). Greater rates of N mineralization (p < 0.001) positively influenced NPP (Figure 5; Figure B2, Appendix B) which drove the observed increases in NEP. The average gains were around 0.01 kg m⁻² year⁻¹ in SSP1-2.6, 0.02–0.03 kg C m⁻² year⁻¹ in SSP3-7.0 and approximately 0.02 kg C m⁻² year⁻¹ in SSP5-8.5 (Table 3). The rates of uptake were higher in southern Sweden compared to in northern Sweden (Tables B2 and B4, Appendix B).

The C stored in tree biomass consistently increased across the country in all emission and policy scenarios when comparing 2081–2100 to 2001–2020 (p < 0.001) (Figure 6 and Table 3). The largest gains were observed in the EUPOL policy, with increases ranging from 1.8 kg C m⁻² in SSP1-2.6 to 2.7 kg C m⁻² in SSP5-8.5. Under BAU the average increases in tree biomass C were 1.2 kg C m⁻² in SSP1-2.6, 1.9 kg C m⁻² in SSP3-7.0 and 2.1 kg C m⁻² in SSP5-8.5. The AR policy produced similar outcomes to BAU, with gains varying from 1.2 kg C m⁻² in SSP1-2.6 to 2.1 kg C m⁻² in SSP5-8.5. Greater spatial variation in gained woody biomass in EUPOL in all three emission scenarios resulted from a larger proportion of set-asides, where stochastic stand-replacing disturbances with a return time of 100 years were enabled, influencing ecosystem structure (Figure 6).

[IMAGE OMITTED. SEE PDF]

As an outcome of increased NPP, and given that the lengths of rotation periods were not altered, the model suggested potential gains in harvestable biomass in all three emission scenarios (p < 0.001). In SSP1-2.6 they ranged from 13% in EUPOL to 28% in both BAU and AR (Figure 7 and Table 3). Under SSP3-7.0, increases varied from 20% in EUPOL to 36% in BAU and AR, and under SSP5-8.5 from 22% in EUPOL to 38% in BAU and AR (Figure 7; Figure B4, Appendix B).

[IMAGE OMITTED. SEE PDF]

Soil C decreased significantly at the end of the century in SSP3-7.0 and in SSP5-8.5. Soil C losses ranged from −0.4 to −0.6 kg C m⁻² in EUPOL (p < 0.01), −0.4 to −0.5 kg C m⁻² in AR (p < 0.001), and −0.3 to −0.4 kg C m⁻² in BAU (p < 0.01) (Figure 7 and Table 3). These reductions in soil C occurred despite the simulated increases in CWD input, and were driven by higher Rh, stimulated by increased soil temperatures and similar or higher precipitation levels (Figure 7). Although soil C losses were observed also in SSP1-2.6, the trends were not significant in BAU, EUPOL or in AR (Table 3).

The SSP1-2.6 trajectory also suggested a possible general improvement of freshwater quality, as the rates of N leaching were significantly lower in all regions and forest policy scenarios (p < 0.001). The reduced N leaching resulted from lower N availability in the soil pool due to reduced deposition of nitrate (NO₃⁻) in all regions (Figure B1, Appendix B) in conjunction with a high continuous plant N demand (Figure 6). Similar or slightly higher N leaching rates as in the historical climate (2001–2020) in SSP3-7.0 and in SSP5-8.5 indicated that the ecosystems were not N saturated, which suggested that all of the N produced as a result of increased N mineralization or input via fixation and deposition was utilized by the vegetation (Figure 6).

Annual input of coarse woody debris increased significantly in EUPOL and in BAU (p < 0.001), but not in AR in SSP1-2.6. In SSP3-7.0 average increases were 42% in BAU, 35% in AR, and 50% in EUPOL (p < 0.001). For SSP5-8.5, input increases were 46% in BAU, 38% in AR, and 52% in EUPOL (Figure 7; Figure B4, Appendix B).

Changes in Predisposition to Storm Damage

The trends for predisposition to storm damage showed similar but slightly more severe outcomes in SSP3-7.0 and in SSP5-8.5 compared to in SSP1-2.6. Under BAU, the predisposition to storm damage increased by 47% (SSP1-2.6), 49% (SSP3-7.0) and 54% (SSP5-8.5) comparing 2081–2100 to 2001–2020. In AR, the corresponding increases ranged from 26% (SSP1-2.6 and SSP3-7.0) to 28% (SSP5-8.5), while in EUPOL, they ranged from 20% (SSP1-2.6) to 22% (SSP3-7.0) and 24% (SSP5-8.5).

While the increases in storm predisposition were significant in central and northern Sweden across all policy and emission scenarios (p < 0.05, Tables B2 and B3, Appendix B), the trends varied for forest policy scenarios in southern Sweden. In this region the storm predisposition index under BAU showed a 61%–76% rise in potential volume loss in the event of a major storm at the end of the century, with severity increasing with higher climate impact. These results were mainly related to the observed biomass increase in storm-sensitive Norway spruce monocultures, which were more common in BAU compared to in AR or EUPOL. In southern Sweden, AR resulted in a 4%–11% higher predisposition to storm damage, whereas EUPOL showed a 7%–12% lower predisposition when 2081–2100 was compared to 2001–2020 (p < 0.05) (Table B4 & Figure 6).

Discussion

Overview

Changes in various ecological functions and in the provisioning of ecosystem services can be expected in Sweden due to alterations in climate over the 21st century (Jönsson et al., 2015; Jönsson & Lagergren, 2018; Pilli et al., 2022). Here we have quantified outcomes for both SSP1-2.6, SSP3-7.0 and SSP5-8.5 using an ecosystem model which combines advanced representations of vegetation dynamics, demography and biogeochemistry with sophisticated and flexible representations of the influence of management on forest composition and structure. The evaluation performed indicated agreement between the modeled and observed forest structure for a number of observed variables along a north to south gradient in Sweden (Figures A2-A6, Appendix A). The modeled average biomass C for the historical period (2001–2020) of 5.56 kg C m⁻² agreed well with previous modeled estimates for tree C storage in production forests in Sweden (Mazziotta et al., 2022). Estimates of NEP of 0.02–0.08 kg C m⁻² (average 2012–2021) for the Scandinavian region have previously been generated with an ensemble of DVMs, where 11 out of 16 models accounted for wood harvest (Friedlingstein et al., 2022). These results indicate a slightly lower C uptake than found within our study, as our findings for NEP showed an average annual C sink in Sweden of 0.09 kg C m⁻² within the historical climate (average for 2001–2020, Table 3).

We attribute the resulting increased biomass accumulation in all regions in BAU to climate change and increased CO₂, as our set management regimes maintained the same management and rotation period length in 2081–2100 as in 2001–2020. Although SSP1-2.6 represents a trajectory where climate change is mitigated to some extent at the end of the century, the positive effects of mid-century climate on forest growth was still evident within older forest age classes in the landscape in 2081–2100. In other words, even within this low-emissions scenario, the mid-century peak in temperature and CO₂ contributed to the noted significant difference in end-of-century biomass C when compared to historical biomass C. Our results regarding end-of-century changes in C stock of between 2.1 and 2.7 kg C m⁻² for SSP5-8.5 are consistent with Anderegg et al. (2022) where an average gain in C biomass of about 3 kg C m⁻² was predicted for boreal Scandinavia, comparing the period 2081–2100 to 1995–2014. Our simulations have provided results for an “optimal outcome,” with quantifications of changes in indicators for ES and EF for undamaged stands within the landscape. Both the magnitude and frequency of disturbances are expected to change in the future, with impacts likely to profoundly affect both the C storage and sequestration capacity of forests (Aldea et al., 2024; Jönsson et al., 2009; Ruosteenoja et al., 2018; Seidl et al., 2014, 2016).

The gradual warming in SSP3-7.0 and in SSP5-8.5 resulted in increased microbial activity and decomposition of soil organic matter, which caused a higher N mineralization on an annual basis. These model results are consistent with findings from incubation experiments performed on samples from boreal forest soils (Xiao et al., 2023). The simulated increase of annual inflow of nitrogen-rich CWD in our study seemed to have increased decomposition rather than to stabilize soil C. In cold-limited boreal regions, soil metabolic activity has shown to be highly sensitive to increased soil temperature (Barcenas-Moreno et al., 2009). The modeled increase in Rh was proportionally greater in northern Sweden, as a result of a more pronounced climate warming in higher latitudes (Loranty et al., 2018). In addition to higher Rh in the high and very high emission scenarios, the model projected end-of-century increases in mean annual NPP of 21%–25% in SSP3-7.0 and of 25%–29% in SSP5-8.5, supporting our first hypothesis. Although NPP has been shown to increase in higher CO₂ environments, the response is also mediated by N availability (McCarthy et al., 2010), which increased considerably in our simulations toward the end of the 21st century (Figure 5). The growth responses to higher CO₂ will therefore likely vary at the site scale and depend on site fertility (Oren et al., 2001). Moreover, increased air temperatures will influence the growth rates and competitive strength of native tree species in Sweden. The change in air temperature under SSP3-7.0 and SSP5-8.5 scenarios in southern Sweden indicated an exceedance of the temperature optimum for growth of Norway spruce saplings by mid-century (Figure 1) (Marchand et al., 2023). Reductions in soil moisture can be expected in southern Sweden with increased climate impact due to higher temperatures and evapotranspiration (Ruosteenoja et al., 2018), which may be further disadvantageous to established stands of Norway spruce. However, although the gradual warming in southern Sweden will favor the more drought-adapted Scots pine over Norway spruce, the silvicultural interest in the former species has been limited due to high browsing pressures from ungulates and lower potential profitability (Lodin et al., 2017).

Implications of Business As Usual (BAU)

Our findings for the BAU scenario suggest significant gains in potential harvests in Sweden. However, the predisposition to storm damage also showed proportional increases of 47% in SSP1-2.6, 49% in SSP3-7.0 and of 54% in SSP5-8.5, highlighting larger risks of losses of C and economic value (Table 3). The predisposition to storm damage was significantly higher in BAU in comparison to in the two other forest policy scenarios, most notably in southern and in central Sweden, which confirmed our second hypothesis. Norway spruce monocultures account for 37% of the total forest area in southern Sweden, and this species has shown higher sensitivity to disturbances such as storms and drought compared to other commercially important tree species in that area (Aldea et al., 2024; Felton et al., 2016). Higher stand heights and total biomass toward the end of the rotations in Norway spruce monocultures compared to in spruce-birch mixtures also contributed to increased predisposition to storm damage. As most storms occur during autumn and winter in Sweden, reductions in the duration of frozen soil is likely to exacerbate storm damage in coniferous forests, which can be expected with gradual climate warming (Jönsson et al., 2015; Peltola et al., 1999). During conditions where N mineralization is stimulated and more N is available for uptake, C investment in root systems may also be reduced compared to during N limited conditions (Noormets et al., 2015). This may increase the ratio of aboveground woody biomass to belowground root biomass (Vicca et al., 2012) which also could have potential implications for the risk of wind throw in storm-sensitive stands. Moreover, greater proportions of storm-damaged forests amplify risks of Spruce bark beetle (Ips typographus) attacks, which may also reduce or offset ecosystem C sinks (Jönsson et al., 2015).

Implications of the EU-Policy Alternative (EUPOL)

The larger proportions of set-asides and CCF in EUPOL were intended to represent landscapes with a higher degree of naturalness, providing older trees and more within-stand structural heterogeneity, along with high levels of dead wood, which has been shown to support greater species richness and species abundance (Brockerhoff et al., 2008; Carnus et al., 2006; Häkkilä et al., 2021). This scenario was created in order to visualize a future where the area of protected forest land increases, and management shifts toward practices of reduced intensity, endorsed by the European Commission in the EU Biodiversity Strategy 2030 and EU Forest Strategy 2030 as a means to increase biodiversity, C storage and C uptake (European Commission, 2020, 2021). We identified the greatest end-of-century gains in C biomass in EUPOL, with average increases in tree biomass C of 32% in SSP1-2.6, of 45% in SSP3-7.0 and of 48% in SSP5-8.5. The proportionally larger gains in tree biomass C in EUPOL compared to in BAU or AR mainly depended on the expansion of unmanaged areas and set-asides which favored the development of older forests and larger sized trees. Our findings corroborate previous modeled scenarios in both Sweden and Finland which have indicated that these areas provide substantial gains in forest biomass C (Mazziotta et al., 2022; SFA, 2022c; Trivino et al., 2023). In SSP1-2.6, the potential harvestable stemwood C increased by 13% despite the greater proportions of set-asides and protected areas. These results suggest that increased productivity from climate change may potentially alleviate the trade-off severity between high timber yields and C storage in some areas. However, the future growth-reducing impacts of natural disturbances are uncertain and therefore also their adverse effects on observed growth increases.

EUPOL showed a significant decrease in the predisposition to storm damage in southern Sweden in SSP1-2.6 and in SSP3-7.0 (Table B4, Appendix B), and the lowest increase in predisposition in central Sweden (Table B3, Appendix B). This suggests that increasing the proportion of set-asides and protected areas in these regions is a viable alternative to achieving synergistic benefits between C storage and multiple additional ES, including the maintenance of habitat (CWD) and nutrient release from decomposition, given that the area of annually felled forest does not increase over time. Our findings corroborate Mazziotta et al. (2022) who similarly showed increases in local scale synergies between ES in set-asides both during medium and high projected climate change in Sweden.

Implications of the Adaptation and Resistance Alternative (AR)

The AR scenario represented a future where mixed forests of spruce-birch, mixed coniferous forests of spruce-pine and broadleaved monocultures form between 39% and 54% of the forest landscape. Stakeholders in Sweden associate these forest types with lower risk of disturbance-related damage and economic losses (Hallberg-Sramek et al., 2023) and non-industrial private forest owners in Sweden have indicated that they would like to have larger proportions of them in the future landscape (Bergkvist et al., 2024). A higher tree species diversity in mixed stands can enhance functional trait diversity, influencing within-stand growing conditions and resource availability (Depauw et al., 2024). In Sweden, mixtures of spruce-pine and spruce-birch have been shown to provide additional benefits to biodiversity, esthetic beauty, recreation and water quality regulation when compared to monocultures of Norway spruce (Felton et al., 2023). Existing studies indicate similar outcomes for mixtures of spruce-birch and spruce-pine as for spruce monocultures in terms of wood production capacity, but the relatively limited amount of field experiments in Scandinavia prevent firm conclusions regarding productivity differences (Felton et al., 2016, 2023). We noted a generally lower predisposition to storm damage in the AR scenario when compared to BAU, primarily in southern and in central Sweden (Tables B3 and B4, Appendix B). Simultaneously, the tree biomass C increase was similar to observed for the more storm-sensitive scenario BAU (Table 3). Damage from storms, biotic pests and diseases may also be reduced in tree species mixtures as the more vulnerable tree species occur less frequently in the mix compared to in a monoculture (Carnus et al., 2006; Felton et al., 2023). Although mixed stands provide a range of benefits, birch or pine in mixtures with spruce have also shown to be susceptible to browsing damage from ungulates in Sweden (Felton et al., 2023), a factor that was not considered in our simulations.

Limitations of the Study

Although our approach accounted for future variation in landscape composition, we assumed no change in the age distribution of stands in each grid cell over time. Hence, forest stands retained the same rotation period length in the future changed climate as during the historical period. Further model development should aim to provide a setting where rotation periods are altered for a forest type in response to changed productivity, resulting from changed climate conditions over the course of the simulation period. With such a setting enabled, harvested C would exhibit greater increases over time with increased forest NPP. In that case, the production forest tree biomass C stock would not show similar increases as in this study, which would have implications for the inflow of CWD to the soil, soil C levels, and the predisposition to storm damage.

Our findings regarding changes in C stocks, harvest potentials and NEP should be interpreted as the best possible outcome given these assumptions, as the simulations did not account for ecosystem C losses due to disturbances such as fires, storms or biotic damage (fungal, insect) in managed forests. However, climate change will likely amplify ecosystem C losses through increasing both the frequency and magnitude of disturbances, which has been suggested to completely offset or even turn forests from C sinks into sources (Gregor et al., 2022; Ma et al., 2012). Indications of the major challenges associated with extreme weather conditions and future increased CO₂ emissions have already been observed in Canada, where forest fires caused losses of about 15 Mha of forest in 2023 (Jain et al., 2024). The risk of fire is generally higher in areas with low forest fragmentation and a high proportion of conifers, as these provide fuels with a greater ignition potential compared to less fragmented forest landscapes with higher proportions of deciduous stands (Felton et al., 2023; Krawchuk et al., 2006). The application of fire suppression techniques will be essential to limit the future size of wildfires, which in Sweden have been relatively effective in areas with denser populations and well-developed road networks (Pinto et al., 2020). In this study, climate related risk was indicated by the predisposition to storm damage, with generally increased predisposition to damage with tree age. Complex interactions between storm damage, drought stress and spruce bark beetle attacks and other biotic agents were not accounted for. Harvest was simulated for managed stands, and a generic representation of stand-replacing disturbances was possible to include for set-asides and unmanaged areas in order to represent natural forest dynamics. Prolonging or shortening the return time of this setting, to mimic climate change impacts on the occurrence of natural disturbances, would further influence the buildup of C storage within these areas.

The model results presented in this study are averages of output from model runs with three differing sets of climate data input, originating from the ESMs MRI-ESM2.0, EC-Earth3-Veg and GFDL-ESM4. Differences in the effective climate sensitivity among the ESMs influenced the model output of LPJ-GUESS, where MRI-ESM2.0 and GFDL-ESM4 showed greater similarities in the representation of NPP, Rh, Biomass C, and NEP compared to EC-Earth3-Veg, which tended to produce higher values of NPP and Rh in all regions (Figure B5, Appendix B). We utilized the three available emission scenarios provided by the ISIMIP repository, which correspond to both low, high and very high continued emissions (Lange & Büchner, 2021). In 2022, the observed annual atmospheric CO₂ growth rate was higher than in any year of the past 800,000 years, lending credibility to the inclusion of these higher emission scenarios in our model study (Friedlingstein et al., 2022).

Although the area of mixtures of spruce-birch, spruce-pine and of spruce managed as CCF changed in AR and in EUPOL compared to in BAU, the transitions were not explicitly simulated within LPJ-GUESS. Instead the transitions were implicitly assumed to take place during the intermediate period 2021–2080, at a rate where all stands were assumed to have been successfully transformed by 2081. During the stand conversion period from Norway spruce monocultures to CCF the growth rate of stands can be reduced by over 30%, which has implications for timber provisioning during this transition period (Drössler et al., 2014). The C balance of a recently transformed stand in 2081–2100 will also be different from in stands where the transformation was completed at a much earlier time. However, these effects on the C balance were not possible to include in this study. Eddy-covariance measurements have indicated that stands located on ditched or drained peat soils may have considerably higher soil respiration rates than observed for mineral soils, which may offset the NEP and turn stands from sinks into sources (Bergkvist et al., 2023; Lagergren et al., 2008, 2019). Histosols were however not represented within our study, as all simulations were performed for forests on mineral soils.

Conclusions

Climate change places new demands on forest owners and managers to improve adaptive capacity and to consider alternatives to traditional forestry practices in Sweden. However, uncertainty regarding the extent of climatic change, and a lack of detailed projections of future potential outcomes and their implications for forests, can make such decisions challenging. Here we have explored end-of-century impacts of both a low (SSP1-2.6), a high (SSP3-7.0) and a very high (SSP5-8.5) emissions scenario. We have also contrasted three alternative forest policy directions in Sweden and assessed outcomes for ecosystem functioning, ecosystem services and for forest storm vulnerability. As hypothesized, the model results indicated a considerably higher NPP in the range of 21%–29% in the high and very high emissions scenarios at the end of the century in Sweden. In the SSP1-2.6 scenario of effective climate mitigation, NPP peaked during mid-century, but statistically significant increases in NPP still occurred in 2081–2100 when compared to 2001–2020 in the AR and EUPOL policy scenarios. The findings indicated consistent increases in forest biomass C in all emission scenarios in Sweden, under the assumption that the rotation period lengths in stands were not altered over time. C gains were the greatest in the more conservation-oriented alternative EUPOL, which also showed considerably higher net N mineralization rates and provisioning of coarse woody debris, but also higher soil respiration rates. In line with our hypothesis, reduced proportions of conifers in the AR and EUPOL forest policy scenarios reduced the predisposition to storm damage when compared to BAU. The AR forest scenario also provided similar gains in biomass C and potential harvest C as in BAU, as well as higher net N mineralization rates. Although our model results are associated with uncertainty, the direction and persistence of the trends may provide relevant information to stakeholders and forest managers regarding both potential benefits and risks. A more in-depth analysis of trade-offs and synergies of simulated outcomes for the different forest policy alternatives in future studies could provide further information and aid decision making. We can conclude that increasing the proportions of set-asides and protected areas provides synergistic benefits between C storage and a range of additional benefits, including improved maintenance of habitat and higher nutrient release from decomposition. It would be valuable to further explore the scenario-specific extent and magnitude of natural disturbances and impacts on ecosystem C gains, as well as the influence of shortened rotations for reducing storm damage predisposition.

Appendix

Appendix A - Supplementary Material for Model Calibration and Evaluation

Here we present details on the set-up of simulations during the evaluation of model performance, and outcomes of the evaluation for stand diameter, height, standing volume and stand density.

A1

Calibration of Allometry

Our aim was to improve model estimates of stand height and diameter by calibrating the two allometric parameters k_allom2 and k_allom3. These two parameters modulate the height-diameter ratio of simulated trees in LPJ-GUESS; A1 $$H={k}_{\text{allom}2}\times {D}^{{k}_{\text{allom}3}}$$ Where D is the diameter (m) and H is the height (m). Our calibration aimed to provide new estimates for k_allom2 and k_allom3 for each of the included tree species in this study in order to better represent observed growth patterns and species-specific traits in the model. Measurements of diameter and height originating from Swedish forests were fitted to a log-linear equation to derive new estimates for k_allom2 and k_allom3 (Huang & Titus, 1992). Scots pine (Pinus sylvestris L.) and European beech (Fagus sylvatica L.) data were available for individual trees (Jucker et al., 2022). The beech data originated from southern Sweden, and the Scots pine data from southern and central Sweden. For Norway spruce, birch (Betula spp.) and Pedunculate oak (Quercus robur L.) only plot-level means were available for calibration of the parameters. The model was run in cohort mode for even-aged stands, which implies that individuals of the same age were identical in size within a patch. For this reason, plot-level means of observed height and diameter were considered reliable for the calibration. The Norway spruce and birch data originated from temporary plots collected during 2007–2021 (SNFI, 2022b) and for oak, from long-term experiments on Visingsö, southern Sweden (Carbonnier, 1975). The birch data did not distinguish between the two birch species present in Sweden, Silver birch (Betula pendula Roth) and Downy birch (Betula pubescens Ehrh.).

A2

Evaluation of Model Performance

Model results with updated allometry estimates were validated for the years 2006–2015 in five different areas of Sweden representing both temperate and boreal zones, in order to capture a gradient in climate conditions (Figure A1). The results were compared against plot-level chronosequences of mean stand height, quadratic mean diameter, stand density, standing volume and height-diameter ratio, originating from the Swedish National Forest Inventory. The comparisons of simulations were made against observations based on stand age. Updated PFT parameters developed for Sweden for Norway spruce and Scots pine, which have shown improved agreement of simulated productivity against observations (Bergkvist et al., 2023), were utilized in the simulations of managed forests both for the evaluation and during transient simulations (see Section 2.3). Monocultures of beech, birch and oak were represented using the original vegetation parameters. Stands were initiated during the nineteenth and twentieth centuries by clear-felling of natural forest and re-planting, allowing one full rotation to occur before the re-establishment of the simulated stands used in the evaluation (Table A2).

[IMAGE OMITTED. SEE PDF]

Table A1 Calibrated Parameter Values for the Two Allometric Parameters in the Height:Diameter Relationship in LPJ-GUESS ($$H={\mathit{k}}_{\mathit{allo}\mathit{m}\mathit{2}}\times {D}^{{\mathit{k}}_{\mathit{allo}\mathit{m}\mathit{3}}}$$, Equation 1)

Species	kallom2	kallom3
Norway spruce	83.6 (76.9–90.9)	0.989 (0.939–1.039)
Scots pine	51.1 (45.4–57.6)	0.821 (0.763–0.879)
Birch spp.	43.4 (36.6–51.2)	0.638 (0.562–0.713)
Pedunculate oak	35 (34.1–36.0)	0.499 (0.481–0.517)
Beech	33.7 (31.4–36.2)	0.446 (0.397–0.495)

Table A2 Forest Management in Stands Utilized in the Evaluation of Height, Diameter, Standing Volume and Stand Density

Forest type	Area	Planting density (stems ha−1)	Pre-comm. Thin.	1st thin.	2nd thin.	3rd thin.	4th thin.
Scots pine monoculture	Dalarna	2500	19 (15%)	57 (30%)	76 (20%)	-	-
Skellefteå	2300	21 (15%)	63 (30%)	84 (20%)	-	-
Norway spruce monoculture	Kronoberg	2500	9 (20%)	45 (25%)	72 (15%)	-	-
Dalarna	2300	15 (20%)	50 (25%)	75 (15%)	-	-
Birch monoculture	Kronoberg	5000	9 (35%)	27 (25%)	36 (20%)	-	-
Västergötland	4500	11 (35%)	33 (25%)	44 (20%)	-	-
Dalarna	4000	13 (35%)	39 (25%)	52 (20%)	-	-
Oak monoculture	Skåne	5000	12 (15%)	Thinnings every 10th year throughout the rotation (variable)
Beech monoculture	Skåne	5000	10 (45%)	20 (45%)	40 (25%)	60 (25%)	80 (20%)

A3

Setup of Forest Management for County-Level Evaluations

Managed monocultures of Scots pine and Norway spruce were simulated with vegetation parameters developed for Sweden (Bergkvist et al., 2023). Random patch-destroying disturbances were disabled during the model runs in order to model the structure of undamaged forests. Coniferous monocultures, birch and beech stands were managed with a pre-commercial thinning followed by 1–4 thinnings (Table A2). Oak stands were established 1896 and onwards, with extraction of wood occurring every tenth year through selective felling. Stands were simulated in 24 patches within all areas.

To generate standing volume, the stem biomass output from LPJ-GUESS, given in kg C m⁻², were first converted to dry matter using constant conversion of C to dry matter of 0.51 kg C per kg dry matter for Norway spruce and Scots pine, and of 0.48 kg C per kg dry matter for birch (Aalde et al., 2006). Then biomass conversion factors determined for age classes were used to assess the stem wood volume (Lehtonen et al., 2004). The wood volume for pedunculate oak and beech was estimated using volumetric functions based on simulated height and diameter.

A4

Results of the Model Evaluation

Model results with updated allometry indicated slight improvements to the modeled height to diameter ratio, most noticeable in Norway spruce monocultures (Figure A2; Table A1). A correlation analysis between simulations with updated allometry and observations did not result in higher correlations when compared to a similar analysis between simulations with original allometry and observations, either for simulated diameter, standing volume, height or density (Figures A3–A6).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Estimated stand height for each simulated forest type showed generally high correlations with observations across a latitudinal gradient from southern to northern Sweden, ranging from 0.77 (oak monocultures in Skåne) to 0.95 (birch monocultures in Västergötland) (Figure A5). Simulated stand diameter displayed a higher correlation with observations in coniferous stands of Scots pine and Norway spruce (R = 0.80 to 0.92) compared to in birch stands (R = 0.54 to 0.74) (Figure A3). The model tended to underestimate the observed standing volume at sites with high standing volumes, and produce higher values than observed for sites with low standing volume (Figure A4). Observed abundant natural regeneration in older broadleaved monocultures reduced the correlation between observed and simulated stand density, as natural regeneration was disabled in our simulations (Figure A6).

Appendix

Appendix B - Supplementary Material for Transient Simulations (2001–2100)

Here we present additional information regarding trends for nitrogen deposition (2011–2100), detailed information on changes in climate and in indicator values for different regions of Sweden. This appendix also includes heatmap visualizations for SSP3-7.0, and a principal component analysis intended to highlight differences between ESMs in their influence on four indicators related to the C cycle.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Table B1 20-Year Mean Annual Air Temperature, Annual Precipitation, and Annual Shortwave Radiation for the Historical Period (2001–2020) and the Far Future (2081–2100) in Southern, Central and Northern Sweden

Region	Air temperature (ᵒC)	Precipitation (mm)	Shortwave radiation (W m−2)
Southern Sweden
Historical period (2001–2020)	7.3	731	115
Far future (2081–2100)
SSP1-2.6	7.7	796	119
SSP3-7.0	10.1	819	116
SSP5-8.5	10.9	781	117
Central Sweden
Historical period (2001–2020)	4.6	787	106
Far future (2081–2100)
SSP1-2.6	5.0	871	109
SSP3-7.0	7.7	980	106
SSP5-8.5	8.5	959	106
Northern Sweden
Historical period (2001–2020)	1.0	722	93
Far future (2081–2100)
SSP1-2.6	1.6	768	95
SSP3-7.0	4.6	859	91
SSP5-8.5	5.4	843	90

Table B2 Average Values for Six Indicators of EF, Seven Indicators of ES, and One Indicator of the Predisposition to Storm Damage in Northern Sweden

Indicator	Historical	SSP1-2.6	SSP3-7.0	SSP5-8.5
BAU	BAU	AR	EUPOL	BAU	AR	EUPOL	BAU	AR	EUPOL
NPP (kg C m−2 year−1)	0.44	0.47	0.48	0.49	0.57	0.58	0.59	0.58	0.60	0.61
Rh (kg C m−2 year−1)	0.37	0.40	0.40	0.42	0.47	0.48	0.49	0.49	0.50	0.51
NEP (kg C m−2 year−1)	0.07	0.07	0.08	0.07	0.09	0.10	0.10	0.09	0.10	0.10
NPP to GPP ratio	0.46	0.45	0.46	0.47	0.42	0.42	0.44	0.41	0.41	0.43
Biomass C (kg C m−2)	4.81	5.93	6.13	6.71	6.78	7.02	7.65	6.96	7.20	7.81
Potential harvest C (kg C m−2)	3.81	4.83	4.97	4.47	5.26	5.41	4.83	5.35	5.50	4.91
Soil C (kg C m−2)	9.6	9.4	9.3	9.2	9.3	9.2	9.1	9.3	9.2	9.1
Biomass C to soil C ratio	0.50	0.63	0.66	0.73	0.73	0.77	0.84	0.75	0.79	0.86
Rh to litter C input ratio	1.42	1.39	1.40	1.39	1.38	1.39	1.37	1.39	1.41	1.39
CWD C input (kg C m−2 year−1)	0.05	0.06	0.06	0.07	0.08	0.08	0.09	0.08	0.08	0.09
Leached N (kg N ha−1 year−1)	5.39	4.67	4.57	4.68	5.73	5.59	5.64	6.18	6.02	6.06
Net N mineralization (kg N ha−1 year−1)	13.9	15.1	16.1	18.3	18.5	19.9	22.3	20.0	21.6	24.2
Water Use Efficiency (g C kg−1 H2O)	8.4	9.1	9.1	8.9	11.7	11.6	11.1	12.3	12.2	11.6
Predisposition to storm damage (m−3 ha−1)	25.9	35.8	35.9	36.6	35.7	35.6	37.2	36.6	36.1	37.3

Table B3 Average Values for Six Indicators of EF, Seven Indicators of ES, and One Indicator of the Predisposition to Storm Damage in Central Sweden

Indicator	Historical	SSP1-2.6	SSP3-7.0	SSP5-8.5
BAU	BAU	AR	EUPOL	BAU	AR	EUPOL	BAU	AR	EUPOL
NPP (kg C m−2 year−1)	0.55	0.57	0.58	0.59	0.65	0.66	0.68	0.67	0.68	0.70
Rh (kg C m−2 year−1)	0.46	0.47	0.48	0.49	0.54	0.55	0.57	0.56	0.57	0.60
NEP (kg C m−2 year−1)	0.09	0.10	0.10	0.10	0.11	0.11	0.11	0.10	0.11	0.10
NPP to GPP ratio	0.44	0.42	0.43	0.44	0.38	0.39	0.40	0.37	0.37	0.38
Biomass C (kg C m−2)	6.27	7.50	7.56	8.13	8.25	8.29	8.92	8.38	8.40	8.99
Potential harvest C (kg C m−2)	4.96	6.35	6.36	5.56	6.69	6.73	5.86	6.78	6.79	5.91
Soil C (kg C m−2)	15.4	15.1	15.0	14.9	14.9	14.9	14.8	14.9	14.9	14.8
Biomass C to soil C ratio	0.50	0.61	0.61	0.66	0.68	0.69	0.74	0.69	0.70	0.75
Rh to litter C input ratio	1.49	1.47	1.47	1.45	1.49	1.50	1.48	1.51	1.52	1.50
CWD C input (kg C m−2 year−1)	0.08	0.08	0.08	0.09	0.10	0.10	0.12	0.11	0.10	0.12
Leached N (kg N ha−1 year−1)	6.62	5.28	5.17	5.15	6.85	6.66	6.61	7.78	7.54	7.56
Net N mineralization (kg N ha−1 year−1)	19.6	21.3	23.9	28.6	26.5	29.9	36.0	29.3	32.8	39.4
Water Use Efficiency (g C kg−1 H2O)	7.7	8.4	8.4	8.1	9.9	9.9	9.6	10.0	10.0	9.6
Predisposition to storm damage (m−3 ha−1)	26.5	39.7	32.9	30.7	40.0	32.4	30.8	40.4	32.8	31.3

Table B4 Average Values for Six Indicators of EF, Seven Indicators of ES, and One Indicator of the Predisposition to Storm Damage in Southern Sweden

Indicator	Historical	SSP1-2.6	SSP3-7.0	SSP5-8.5
BAU	BAU	AR	EUPOL	BAU	AR	EUPOL	BAU	AR	EUPOL
NPP (kg C m−2 year−1)	0.63	0.64	0.65	0.65	0.71	0.73	0.72	0.73	0.74	0.74
Rh (kg C m−2 year−1)	0.51	0.50	0.52	0.53	0.58	0.59	0.60	0.59	0.60	0.62
NEP (kg C m−2 year−1)	0.12	0.14	0.14	0.13	0.14	0.14	0.12	0.13	0.14	0.12
NPP to GPP ratio	0.44	0.43	0.42	0.42	0.39	0.38	0.38	0.38	0.37	0.37
Biomass C (kg C m−2)	6.40	7.65	7.46	7.81	8.13	7.87	8.14	8.26	7.97	8.26
Potential harvest C (kg C m−2)	5.21	6.73	6.44	5.55	7.01	6.73	5.81	7.11	6.84	5.90
Soil C (kg C m−2)	14.4	14.0	14.0	14.0	14.0	14.2	14.1	13.8	13.8	13.7
Biomass C to soil C ratio	0.55	0.68	0.66	0.69	0.73	0.70	0.72	0.75	0.73	0.76
Rh to litter C input ratio	1.58	1.59	1.53	1.50	1.66	1.62	1.59	1.68	1.64	1.61
CWD C input (kg C m−2 year−1)	0.10	0.10	0.09	0.09	0.12	0.11	0.11	0.12	0.11	0.11
Leached N (kg N ha−1 year−1)	7.70	6.24	6.06	6.13	8.91	8.74	8.80	9.57	9.19	9.25
Net N mineralization (kg N ha−1 year−1)	24.6	27.3	32.9	38.6	34.6	42.3	49.6	36.7	43.6	51.0
Water Use Efficiency (g C kg−1 H2O)	7.4	8.4	8.5	8.4	9.4	9.6	9.5	9.5	9.6	9.4
Predisposition to storm damage (m−3 ha−1)	31.0	49.9	32.4	27.3	52.1	33.5	28.5	54.6	34.5	28.9

Acknowledgments

The authors would like to thank Sam Rabin and Adrian Gustafson for aggregating climate data for the study. This work was made possible by FORMAS, grant number 2019-01968. D.W. and P.A.M. acknowledge financial support from the Strategic Research Area MERGE (Modeling the Regional and Global Earth System—), and F.L. acknowledges the Horizon project AVENGERS.

Conflict of Interest

The authors declare that they have no competing interests.

Data Availability Statement

LPJ-GUESS version 4.1 was used to generate the results of this study. LPJ-GUESS is an open-access software available for use under Mozilla Public License v. 2.0 (no registration required) (Nord et al., 2021). Input data for nitrogen deposition for 1850–2099 were provided by the National Center for Atmospheric Research (NCAR), available for use under Creative Commons CC BY-SA 4.0 (no registration required) (Hegglin et al., 2018). Tree-level observational data from the Tallo database (v1.0.0) were used to calibrate two parameters in the model. The data are available under Creative Commons 4.0 (no registration required) (Jucker et al., 2022). Climate input data were provided by the Inter-Sectoral Impact Model Intercomparison Project version 3b (ISIMIP3b, v1.1), available under CC0 1.0 Universal Public Domain Dedication (CC0 1.0) (no registration required) (Lange & Büchner, 2021). Observational data on stand structure from the Swedish National Forest Inventory, collected between 2007 and 2021, were used to evaluate the model. This data is freely available for download (no registration required) (SNFI, 2022b). A data set containing the code to optimize the allometric parameters, to plot the figures, and the details of the underlying data analysis is provided under Creative Commons Attribution 4.0 International (no registration required) (Bergkvist, 2024).

References

Aalde, H., Gonzalez, P., Gytarsky, M., Krug, T., Kurz, W. A., Lasco, R. D., & Verchot, L. (2006). Generic methodologies applicable to multiple land‐use categories, 4. In IPCC Guidelines for National Greenhouse Gas Inventories (pp. 1–59). Institute for Global Environmental Strategies IGES.

Ahlström, A., Canadell, J. G., & Metcalfe, D. B. (2022). Widespread unquantified conversion of old boreal forests to plantations. Earth's Future, 10(11). [DOI: https://dx.doi.org/10.1029/2022EF003221]