Future methane fluxes of peatlands are controlled by management practices and fluctuations in hydrological conditions due to climatic variability

Peatland management practices, such as drainage and restoration, have a strong effect on boreal peatland methane (CH₄) fluxes. Furthermore, CH₄ fluxes are strongly controlled by local environmental conditions, such as soil hydrology, temperature and vegetation, which are all experiencing considerable changes due to climate change. Both management practices and climate change are expected to influence peatland CH₄ fluxes during this century, but the magnitude and net impact of these changes is still insufficiently understood. In this study, we simulated the impacts of two forest management practices, rotational forestry and continuous cover forestry, as well as peatland restoration, on hypothetical forestry-drained peatlands across Finland using the land surface model JSBACH (Jena Scheme for Biosphere–Atmosphere Coupling in Hamburg) coupled with the soil carbon model YASSO and a peatland methane model HIMMELI (Helsinki Model of Methane Buildup and Emission for Peatlands). We further simulated the impacts of climatic warming using two RCP (Representative Concentration Pathway) emission scenarios, RCP2.6 and RCP4.5. We investigated the responses of CH₄ fluxes, soil water-table level (WTL), soil temperatures and soil carbon dynamics to changes in management practices and climate. Our results show that management practices have a strong impact on peatland WTLs and CH₄ emissions that continues for several decades, with emissions increasing after restoration and clearcutting. Towards the end of the century, WTLs increase slightly, likely due to increasing precipitation. CH₄ fluxes have opposing trends in restored and drained peatlands. In restored peatlands, CH₄ emissions decrease towards the end of the century following decomposition of harvest residue in the top peat layers despite increasing WTLs, while in drained peatland forests sinks get weaker and occasional emissions become more common, likely due to rising WTLs and soil temperatures. The strength of these trends varies across the country, with CH₄ emissions from restored peatlands decreasing more strongly in southern Finland, and forest soil CH₄ sinks weakening most in northern Finland.

1 Introduction

Boreal peatlands are considerable sinks of carbon, storing approximately 270–370 Pg of carbon (Turunen et al., 2002). However, they are also a large source of methane (CH₄) (Turetsky et al., 2014; Abdalla et al., 2016), the second-most-important greenhouse gas after carbon dioxide (IPCC, 2013). In total, CH₄ emissions from northern peatlands form approximately 20 % of global wetland emissions and are a notable source of uncertainty in the global methane budget (Saunois et al., 2016). The magnitude of CH₄ fluxes varies strongly depending on local factors, particularly soil water-table level (WTL) and soil temperature, controlling CH₄ emissions (Christensen et al., 2003; Turetsky et al., 2014). While pristine peatlands are strong sources of CH₄, drained and managed peatlands can turn to small sinks of CH₄ when the lowered WTL facilitates oxic conditions and CH₄ oxidation in a thicker surface peat layer. However, drained peatlands also simultaneously often turn into sources of CO₂ due to increased oxic soil respiration (Ojanen et al., 2010; Korkiakoski et al., 2019). On the other hand, restoration of drained peatlands to wetlands can reverse these changes, turning peatlands to sources of CH₄ and sinks of CO₂ (Wilson et al., 2016). In order to understand these trade-offs, the role of CH₄ fluxes in the peatland carbon balance and the climate impacts of peatland management, it is important to study peatland CH₄ fluxes in managed peatlands for several decades after the harvest and restoration management processes.

CH₄ is produced in peatlands by microbes when soil carbon compounds decompose in anoxic conditions, typically in the water-logged peat layer (e.g. Lai, 2009). From the anoxic layer, CH₄ is transported to the atmosphere directly via plants and through the soil layers above through diffusion and ebullition. Microbes also oxidize CH₄, largely in the oxic soil layer above the WTL (Xu et al., 2016). WTL is thus an important factor controlling the CH₄ flux from peatlands, as it controls the thicknesses of both anoxic and oxic layers where CH₄ is produced and oxidized, respectively. However, previous studies have found that in wetlands where the WTL stays constantly high, other factors become more important in controlling the variation in CH₄ flux (Olefeldt et al., 2013; Turetsky et al., 2014). Soil temperature has been found to be a significant controller of peatland CH₄ flux, impacting soil microbial activity (Bubier et al., 1995; Zhao et al., 2016). Additionally, vegetation properties such as plant productivity and species composition control both methane production and transport (Dorodnikov et al., 2011; Turetsky et al., 2014).

While pristine peatlands are a source of CH₄, they are typically a net sink of carbon, as soil organic matter gradually accumulates in the deep anoxic peat layers where decomposition is very slow (Turunen et al., 2002; Nilsson et al., 2008). However, globally, approximately 15 Mha of boreal peatlands have been drained for forestry (Paavilainen and Päivänen, 1995; Päivänen and Hånell, 2012), which has a considerable impact on the ecosystem carbon balance (Ojanen et al., 2013; Korkiakoski et al., 2019; Mäkipää et al., 2023). Nearly one-third of these peatlands (4.7 Mha) is located in Finland. The most common forest management option in peatland forests is rotational forestry, which involves clearcutting the forest and requires ditches to control the WTL in order to maintain a sufficiently low WTL for forest production after forest harvesting (Paavilainen and Päivänen, 1995; Nieminen et al., 2018). While lowering the WTL decreases CH₄ emissions and may turn the soil into a small CH₄ sink (Ojanen et al., 2010), it simultaneously increases the soil CO₂ emissions by enhancing soil organic matter decomposition (Ojanen et al., 2013; Korkiakoski et al., 2019). This can turn the forest into a source of atmospheric CO₂ or decrease its sink (Ojanen et al., 2010, 2013; Hommeltenberg et al., 2014). This effect is particularly strong in nutrient-rich peatlands that receive their water from both precipitation and groundwater and that support faster forest growth but also faster decomposition of soil carbon (Meyer et al., 2013; Ojanen and Minkkinen, 2019). In nutrient-poor peatlands, which receive water solely from precipitation and thus have limited nutrient availability, this effect tends to be smaller, although in forestry-drained peatlands, fertilization also plays a role in nutrient availability and its impacts (Ojanen et al., 2013; Minkkinen et al., 2018). Furthermore, clearcutting and ditch maintenance impose a heavy nutrient and carbon load on local water bodies, decreasing water quality considerably (Nieminen, 2004). Therefore, management options that do not involve regular clearcuts or ditch maintenance have been suggested in order to mitigate the climatic and environmental impacts of rotational forestry in boreal peatlands. In this study, in addition to rotational forestry, we consider continuous cover forestry and restoration of wetlands as alternative management solutions, as these have been put forward as methods to mitigate the climatic and environmental impacts of peatland forestry (Nieminen et al., 2018; Günther et al., 2020).

Continuous cover forestry is an alternative forest management option to rotational forestry (Nieminen et al., 2018). The precise harvesting methods vary, but in effect, only part of the forest stand is removed at one time, leading to a heterogeneous forest structure. Continuous cover forestry decreases the need to maintain ditches, as the continuous forest cover upholds a reasonably low WTL through evapotranspiration (Pothier et al., 2003; Leppä et al., 2020). It can also improve the peatland carbon balance, as increased CH₄ and CO₂ emissions following clearcut harvests can be avoided (Korkiakoski et al., 2020, 2023).

To fully restore the ecohydrological conditions and functions of wetlands, restoration of drained peatlands by rewetting has been performed (Menberu et al., 2016; Günther et al., 2020). In drained peatland forests, this typically means at least reducing drainage by blocking ditches, as well as removing or reducing tree cover, which should lead to a significant rise in the WTL and a gradual return of wetland vegetation (Tarvainen et al., 2013; Maanavilja et al., 2015; Menberu et al., 2016).

Besides changes in the management of drained peatlands, the rapidly changing climate will also influence future CH₄ fluxes. In Finland, mean annual temperatures are rising at twice the rate of global averages (Mikkonen et al., 2015). Precipitation and also evapotranspiration are expected to increase, causing changes in hydrological conditions as well (Ruosteenoja and Jylhä, 2021). Previous studies have shown varying responses of peatland CH₄ fluxes to warming, largely depending on the simultaneous changes in local WTLs (Turetsky et al., 2008; Laine et al., 2019; Peltoniemi et al., 2016). However, the combined effects of peatland management practices and climate change on CH₄ fluxes are insufficiently understood.

Understanding the combined current and future effects of climate and land-use change in any ecosystem is challenging due to the various processes and feedback loops between different variables. To answer this challenge, several complex process-based models have been developed over the years. These types of models aim to describe the relevant biogeochemical and physical processes of a system, thus simulating how the system functions and how it responds to changes in the surrounding conditions. While these models can be useful in, for example, testing hypotheses and understanding future conditions, they can also contain considerable sources of uncertainty due to their complexity. Furthermore, considering the level of detail in how the simulation is created is important. Larger-scale processes and variations in climatical regions can be studied with regionally averaged data, while more local processes require specific site-level information. Careful consideration of the reliability of the results and how closely the simulated results describe real-life conditions is necessary to interpret model studies correctly. Despite their drawbacks, these models are necessary tools to understand the impacts of human actions on ecosystems, such as drained peatlands (Lehtonen et al., 2023).

In this study, we have simulated what would happen to CH₄ fluxes in a forestry-drained, nutrient-rich peatland during the 21st century if in 2020 a management decision was made to (1) continue rotational forestry, (2) shift to continuous cover forestry or (3) restore the area to a wetland. The simulations were run throughout Finland at a regional level under two climate scenarios, RCP2.6 and RCP4.5. We used the land surface model Jena Scheme for Biosphere–Atmosphere Coupling in Hamburg (JSBACH), which is driven by daily climate model data and simulates hydrology, vegetation dynamics and carbon balance, as well as the response of the ecosystem to climatic changes, extreme events and management options. We coupled JSBACH with the YASSO soil carbon model and Helsinki Model of Methane Buildup and Emission for Peatlands (HIMMELI) to simulate relevant processes in peatland carbon cycling. In this study, we investigated the simulated CH₄ flux as well as soil WTL, temperature and soil carbon pools to understand (1) the combined effects of peatland management practices and climate change on future CH₄ fluxes and (2) how CH₄ fluxes vary across Finland.

2 Model description and methodology

2.1 Land surface model JSBACH

JSBACH (Reick et al., 2013) is the land surface model of the Max Planck Institute for Meteorology Earth System Model (MPI-ESM) (Giorgetta et al., 2013) that simulates terrestrial energy, hydrology and carbon fluxes. Sub-grid-scale heterogeneity is described through different vegetation types, in JSBACH called plant functional types (PFTs), which are represented in the model through separate tiles within each grid cell. These are linked with a set of properties, such as phenology type or albedo, that relate the PFTs to the processes accounted for in JSBACH (Reick et al., 2013). A detailed description of the whole model can be found in Reick et al. (2021). In this work we made site simulations with only one PFT per site, assigning the extratropical evergreen PFT, which corresponds reasonably well with Finnish peatland forests, for forests and peatland vegetation for restored and pristine peatland runs. We used a version of JSBACH3 that has been connected with the YASSO soil carbon model (JSBACH-PEAT; Goll et al., 2015) and that we applied to peatlands with WTL dynamics following Kleinen et al. (2020). Additionally, we coupled JSBACH-PEAT with HIMMELI to simulate peatland methane dynamics (Raivonen et al., 2017). Then, to simulate forest management on peatlands, we coupled JSBACH-PEAT with another model set-up of JSBACH that can account for forest growth within the forest PFTs in a similar way to Nabel et al. (2020) and is hereafter referred to as JSBACH-FOM. The peat layer is transferred between model setups. Below, we discuss the model parts that have been modified and are most relevant for this study (Fig. 1a).

Figure 1 Schematics of (a) how the models JSBACH, YASSO and HIMMELI are linked to each other and (b) the experiment timeline in each model run. The acronyms in the figure are T_a – air temperature; P_a – air pressure; Pr – precipitation; RH – relative humidity; SW and LW – shortwave and longwave radiation; WS – wind speed; NPP – net primary production; LAI – leaf area index; WTL – water-table level; and C pools – carbon pools.

2.1.1 Soil water-table level control

The vertical soil moisture dynamics (vertical diffusion, gravitational drainage, and water inputs and losses at different layers) in JSBACH are described through a one-dimensional Richards equation (Reick et al., 2021),

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(D{\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}\right)+{\displaystyle \frac{\partial K}{\partial z}}+S,\end{array}$$

where D(z,t) is the soil water diffusivity; K(z,t) the soil hydraulic conductivity; and S(z) the water inputs and losses from precipitation, evapotranspiration, surface runoff and drainage. Li et al. (2024) implemented soil-moisture-dependent WTL in the peatland-YASSO in a drained peatland to partition the peat soil into anaerobic and aerobic fractions. We further coupled the peatland forest evapotranspiration with WTL and allowed WTL movement deeper down in the total peat column, which is important in drained peatlands. WTL was estimated from the peat column water volume that is controlled by liquid precipitation, snowmelt, evapotranspiration (ET) and run-off. The formulation of WTL used for the pristine wetland was made following the approach in Wania et al. (2009), and the potential evapotranspiration was used to estimate ET. In our implementation for drained peatlands, the actual simulated forest ET was used to drive the water balance. This change was made to account for the impact of forest growth on the ET. The range of WTL was increased to 0.95 m to allow WTL to fluctuate in a deeper layer than the 0.3 m of the pristine wetland set-up. Moreover, the minimum fractional water content was increased to 0.65 from its default value of 0.25, implying that in our formulation it does not represent the original physical definition of the water holding capacity of sphagnum peat but rather has to be considered a tuning parameter for adjusting the WTL variability to its observed level and range of variability.

2.1.2 YASSO soil carbon model

To simulate soil organic carbon and its decomposition, we used the YASSO soil carbon model that has been coupled with JSBACH (Goll et al., 2015). YASSO divides soil carbon into slowly and rapidly decomposing pools consisting of carbon originating from either woody or non-woody plant structural parts. The four rapidly decomposing pools (i.e. AWEN; acid-soluble, water-soluble, ethanol-soluble and non-soluble) take into account the chemical composition of the litter input that in turn depends on the PFTs. Each pool consists of an aboveground and a belowground part and in addition, there is the slowly decaying pool, called the catotelm pool. The soil carbon in each pool is determined by Eq. (2) (Viskari et al., 2022):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial x\left(t\right)}{\partial t}}=M(\mathit{\theta },c)x(t{)}^T+b(t),\end{array}$$

where the operator M is the product of decomposition and mass fluxes between compartments, and b(t) is the litter input to the soil. θ represents the parameters driving decomposition and c the factors controlling decomposition. The model parameters have been determined based on litter decomposition measurements worldwide in order to obtain realistic heterotrophic respiration rates for a range of conditions. In JSBACH-FOM, the model parameters follow the calibration done in forest soils by Tuomi et al. (2009) (Table A1 in the Appendix). In JSBACH-PEAT, parameters follow the parametrization used by Li et al. (2024) (Table A2 in the Appendix). A more detailed description of the model and parameter optimization can be found in Tuomi et al. (2009).

In the peatland implementation, the WTL further divides soil into oxic and anoxic layers. The decomposition of the anoxic fraction is slowed down from the oxic reference decomposition (Tables A1 and A2). The slowly decaying pool is the bottom layer and the fast-decaying layer on top of that consists of belowground and aboveground AWEN pools, in that order. The thickness of the anoxic and oxic layers is determined by the carbon contents of the pools and their bulk densities. The oxic and anoxic fractions of each carbon pool are determined by WTL depth in the soil.

2.1.3 JSBACH-FOM

To simulate forest growth and harvest, we used JSBACH-FOM. It accounts for the age of the forest PFTs and has a different control of maximum forest leaf area index (LAI) compared to JSBACH-PEAT, where maximum LAI is only dependent on the plant functional type. In JSBACH-FOM, it is dependent on available leaf biomass, which is used to simulate forest regrowth (see Nabel et al., 2020; Wey et al., 2022). The change in total vegetation carbon is used as a basis for the calculation of the growing forest. The number of trees per area, or stem number, is calculated from the total vegetation carbon, assuming that the forests are in a self-thinning state. The biomass per individual tree can then be used together with allometric relationships to derive the maximum LAI of the forest. The implementation of forest harvesting and forest growth based on a maximum LAI is described in detail by Nabel et al. (2020). The setup of the model parameters was adjusted to be suitable for our simulations and is explained in detail in Appendix B. The harvest is done for each tile when the forest age reaches the preset rotation time. Then, the forest is clearcut at the beginning of a year, and the age and size restart from zero. The forest stand carbon pools are redistributed due to harvesting. The harvested carbon, making up 77 % of the aboveground woody pool, is removed from the calculations. Half of the green carbon, both the above- and belowground vegetation parts, goes into the YASSO AWEN litter pools according to predefined fractions. The other half is sent directly to the catotelm pool. Similarly, half of the belowground woody carbon, accounting for 30 % of the growing stand total woody carbon, is distributed into the belowground litter pools, while the remaining half goes to the catotelm pool. In addition, the remaining 23 % of the aboveground woody carbon, accounting for the aboveground growing stand, is distributed into the belowground litter pools (50 %) and the catotelm pool (50 %). The redistribution of the cut forest stand carbon to the soil carbon pools in the context of selection harvests is done in a similar way.

2.1.4 HIMMELI methane model

To simulate CH₄ fluxes in peat soils, we used HIMMELI, which has been developed to simulate the buildup, transport and oxidation of CH₄ in peat soils (Raivonen et al., 2017). The concentration of CH₄ at soil depth z is described through Eq. (3):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}{\mathrm{C}}_{{\mathrm{CH}}_{\mathrm{4}}}(t,z)& =-{\displaystyle \frac{\partial }{\partial z}}{F}_{{\mathrm{CH}}_{\mathrm{4}}}-Q_{\mathrm{plt},{\mathrm{CH}}_{\mathrm{4}}}\\ & -Q_{\mathrm{ebu},{\mathrm{CH}}_{\mathrm{4}}}+R_{{\mathrm{CH}}_{\mathrm{4}}}-R_{\mathrm{O}},\end{array}\end{array}$$

where $F_{{\mathrm{CH}}_{\mathrm{4}}}$ is the diffusive flux of CH₄ in peat, $Q_{\mathrm{plt},{\mathrm{CH}}_{\mathrm{4}}}$ is the transport rate via plant roots and $W_{\mathrm{ebu},{\mathrm{CH}}_{\mathrm{4}}}$ via ebullition, $R_{{\mathrm{CH}}_{\mathrm{4}}}$ is the production rate of CH₄, and R_O is the oxidation rate. HIMMELI is driven by soil temperature, WTL, LAI of aerenchymatous plants and the rate of anaerobic soil respiration. HIMMELI simulates several microbial processes and transport pathways of CH₄, CO₂ and O₂, mainly CH₄ production and oxidation, aerobic respiration, ebullition, gas diffusion within the peat layer, and transport in the plant aerenchyma. Unlike in the HIMMELI version used in JSBACH-PEAT, in JSBACH-FOM, gas transport within plant aerenchyma was not included, as there are very few vascular plants with aerenchyma in drained peatlands (Laiho et al., 2003; Päivänen and Hånell, 2012). The methane model parameters used with JSBACH-FOM and pristine land simulations are given in Table A3.

We used a HIMMELI version that was modified by Li et al. (2024) to better suit the simulation of CH₄ uptake as well, which is relevant in a drained peatland. The modified version differs from the original in how the concentrations of compounds in soil layers are treated in the case of lowering the water table.

2.2 Experiment design

The two setups of JSBACH used in this study differ in terms of biomass dynamics and soil carbon model parametrization. In JSBACH-FOM, the PFT was extratropical evergreen, and in JSBACH-PEAT, it was wetland vegetation that described vegetation typical to natural peatlands. Peatland parameters from Hagemann and Stacke (2015) were used to describe the soil properties, e.g. soil porosity, saturated hydraulic conductivity, field capacity and wilting points, and saturated moisture potential. Maximum root depth was set to 1.6 m because the model top soil layers tend to occasionally become relatively dry and thus unrealistically limit photosynthesis.

The simulations were run forcing JSBACH with the regional EURO-CORDEX daily resolution climate data from the centre point of each mainland administrative district in Finland (Fig. 2) (Jacob et al., 2014). The regional CORDEX models are forced by coarse-resolution global CMIP5 climate models. The EURO-CORDEX models are validated and bias-corrected by Finnish observations (Fig. 2) (Räisänen and Räty, 2013; Räty et al., 2014). Three climate models (CanESM2, MIROC5 and CNRM-CM5) based on two emission scenarios, RCP2.6 and RCP4.5, were used in order to better understand climate-related uncertainty in the results. For RCP2.6, we only used the models MIROC5 and CNRM-CM5, as not all required drivers were available in CanESM2.

Figure 2 Map of the Finnish mainland regions and the locations of the chamber measurement sites.

The initial state of JSBACH can be adopted from observations or produced in a spinup run, where selected state variables are usually taken to an equilibrium state under a given climate. For certain very slowly evolving state variables, such as carbon storage in pristine peatlands, a spinup to equilibrium would require an unrealistically long time, and spinup runs have to be interrupted prior to equilibrium. Thus, to account for soil carbon accumulation, we used a spinup of 10 000 years, the time period since the last ice age, running only the soil carbon processes of the model system using climate data from the years 1900–1930 together with the net primary production (NPP) produced with the respective climate and CO₂ concentration of the year 1900. The resulting carbon stock from the spinup run was approximately 110–120 kg m⁻² (Table A4) and the peat depth was 3.4–3.9 m, which is within the range expected from previous studies that have estimated the total carbon stock of peatlands in Finland (Turunen et al., 2002; Juutinen et al., 2013).

To account for the impact of the transient changes in the climate and the CO₂ concentration on the system state, we continued a full pristine wetland run from 1900 to 1950 with increasing CO₂ concentration (Fig. 1b). In 1950, the model version was changed to JSBACH-FOM, conifer seedlings were planted and the model was run until 2020 when the forest was 70 years old, with 640 trees per hectare. JSBACH-FOM does not explicitly simulate ditches, but model drainage and runoff do remove water from the soil. In 2020, each simulation was split into three management options that were continued until the end of the century. The first option was rotational forestry, in which the forest was harvested every 70 years, i.e. in 2020 and 2090. The harvest in this option was performed as a clearcut, where all the living stand was removed, and the harvest residue was relocated to the YASSO pools.

In the second option, continuous cover forestry, we performed selection harvesting by removing 50 % of the woody biomass every 15 years (Juutinen et al., 2021). The harvest was simulated using JSBACH-FOM without applying the clearcut–growth cycles controlled through FOM. Instead, the state of the stand carbon storage and the consequent changes in the soil carbon storage were modified in simulation restarts timed to take place at constant time intervals. In each restart, a pre-determined fraction of the stand biomass was removed, with identical relocation fractions of the biomass sent to the wood products and soil carbon pools as in the baseline clearcut case. Effectively, the manipulation returns the stands to an earlier growth phase of the stand growth curve.

The third option was restoration to an undrained peatland, in which there was a clearcut in 2020, and the simulation was continued using JSBACH-PEAT with wetland vegetation cover until 2100.

2.3 Flux evaluation data

To validate the simulated methane fluxes, we used manual soil chamber measurements of CH₄ fluxes from 5 forestry-drained, 27 restored and 6 pristine peatland sites in Finland that were established for various research projects (Fig. 2). Measurement points had intact ground vegetation and tree roots; thus, the measurements include all components of CH₄ flux between the forest floor and the atmosphere.

Fluxes were measured using portable greenhouse gas measurement devices (LI-7810, LI-COR; M-GGA-918, ABB and Gasmet DX4015, Gasmet), except for the oldest measurements for which gas samples were taken using syringes and analysed in the laboratory of the Natural Resources Institute Finland with a gas chromatograph equipped with a flame ionization detector (FID) for methane. For all measurements, a similar opaque round metal chamber (height 30 cm, diameter 31.5 cm) equipped with a fan for air mixing was utilized. Measurements were carried out on a biweekly-to-monthly interval during the snow-free season (May–October) in 2021. A single measurement lasted either 2–5 min (portable devices) or 25 min (manual sampling for laboratory analysis).

In addition to the CH₄ flux, WTL and soil temperature at 5 cm in depth during the flux measurements were measured. WTL was measured from a well that was a perforated plastic tube installed into the soil when establishing the study sites. It was measured either manually or utilizing the Odyssey capacitance water level loggers (Dataflow Systems Ltd). Soil temperature was measured utilizing either manual temperature probes or iButton DS1921G loggers (Maxim Integrated).

The five forestry-drained sites were located in southern to central Finland (four sites) and northern Finland (one site). Each site had a control treatment without logging and a partial harvest treatment (thinning, overstorey harvesting or strip cutting). Three of the sites had an additional clearcut treatment. The flux measurements ranged from the first to the ninth year after cuttings during the years 2016–2021, depending on the site. The experiment setup and thus also the number of measurement points varied from site to site, yet there were always several points per treatment covering the typical soil moisture variations from the vicinity of the nearest ditch to the mid-strip. All the forestry-drained sites had been drained for several decades before the measurements. They were originally drained for practical forestry purposes, with a typical ditch spacing of 40 m and a ditch depth of ca. 1 m. The sites represented a wide range of site types from oligotrophic to eutrophic peatlands.

The restored and pristine sites were located in southern and central Finland. The restored sites were rewetted, by either damming or filling up the ditches between 1993 and 2020. The sites had previously been drained for forestry; thus, they had similar ditch spacing and ditch depth as the drained sites in this study. During restoration, ditch banks were cleared of trees, but no other tree stand management was done. CH₄ flux measurements were done during 2021 from several points per site. At the restored sites, measurement points were located on strips and on formerly filled or dammed ditches. The sites again represented a wide range of site types from oligotrophic to eutrophic.

As the aim of this paper was not to calibrate the model to these sites specifically but to understand the average drained peatland fluxes across Finland, the measurements were mainly used to check that the simulated values fell within the expected ranges and that the responses to major environmental variables controlling CH₄ fluxes were similar.

2.4 Data analysis

In addition to the simulated CH₄ fluxes, we investigated management impacts on peatland WTL, carbon pools in the upper (fast-decaying) peat layers, LAI and soil temperature in order to understand the controls of CH₄.

To estimate trends in CH₄ fluxes and the environmental variables controlling them between 2020 and 2100, we used the Mann–Kendall trend test, which is used to determine whether there is a monotonic upward or downward trend in a time series (Mann, 1945; Kendall, 1948). The test was calculated for the continuous cover forestry and restoration options and for each climate model and region separately. The trends were not calculated for rotational forestry, as the clearcutting in 2090 interrupts any linear trends in the time series. The magnitude of the trend was estimated using Sen's slope estimator, which is used to quantify significant linear trends in a time series (Sen, 1968). Both were calculated using the “trend” package, version 1.1.5. in R (Pohlert, 2023).

3 Results

3.1 Model validation

In both forestry-drained and restored peatlands, the range of simulated CH₄ fluxes was mostly within the measured values, and the responses to soil temperature and water-table level changes were similar in measured and simulated fluxes (Fig. 3). In restored peatlands, both models and measurements showed emissions increasing with increasing soil temperatures (Fig. 3a and Table A5). In forested peatlands, measurements and model results showed a weak sink of CH₄ when WTL was low (below −15 cm) and that the sink weakened with higher WTL (Fig. 3b and Table A6). The measurements showed occasional emissions regardless of WTL, while the model results showed emissions only in the clearcut option when WTL was high. The CH₄ sinks were strongest in the mature forest (control sites) in both simulations and measurements. Average simulated WTL after restoration matched well with measured WTL, although there was more variation in the measured WTL, particularly right after the restoration.

Figure 3 Comparisons between modelled and measured daily CH₄ fluxes and water-table level. Panel (a) shows CH₄ fluxes according to soil temperature classes in restored peatlands, and panel (b) shows CH₄ fluxes according to water-table level classes in forestry-drained peatlands. Panel (c) shows the WTL before and after restoration.

3.2 Environmental controls

Clearcutting and restoration had the largest impacts on the WTL and fast-decaying carbon pools but had no impact on summer soil temperatures (Fig. 4). Summer WTL rose on average by 30 cm right after the clearcut and by 20 cm after restoration, while there was only a slight rise in WTL after continuous cover forestry (Fig. 4). After the clearcut, WTL remained higher than in the continuous cover forestry option for 20–30 years and then decreased before the next clearcut in 2090. The higher WTL after a clearcut compared to restoration was due to differences in the LAI – in JSBACH-PEAT, LAI recovered rapidly after the restoration, while after the clearcut, forest LAI took considerably longer to grow, keeping transpiration rates lower during the first years. Fast-decaying pools increased following all management practices due to harvest residue but started to decrease due to the decomposition of the harvest residue in the rotational forestry and restoration options.

Figure 4 The impact of forest management on summer WTL, soil temperature (−22 cm below ground) and fast-decaying carbon pools under RCP scenario 4.5. Each boxplot shows the variation created by yearly variation and differences in the three climate models. Two regions, southwest Finland and Lapland, are presented to show regional variability across Finland. The first column (2010–2019) shows the situation in a mature forest before the beginning of the different management scenarios. The next two columns (2020–2024 and 2025–2029) show the situation during the first decade of the different management scenarios, and the last four columns (2040–2049, 2060–2069, 2080–2089 and 2090–2099) show the development of the variables during the later decades.

Towards the end of the century under RCP4.5, there was a slight positive trend in WTL in the continuous cover forestry option, particularly in northern and central Finland, although the trend was not significant in all models, and the estimated Sen's slopes varied considerably between models (Table A8). In the restored option, Uusimaa and other southern regions had a slight positive trend in WTL as well, but western coastal regions had a slight negative trend that was also not significant in all models (Table A7). Fast-decaying pools had a clear negative (positive) trend throughout Finland in the restoration (continuous cover forestry) option but the trend was stronger (weaker) in southern and central Finland. Summer soil temperatures increased throughout the country with little regional variation in both management options.

Figure 5 Average annual CH₄ flux in (a) restored peatlands and (b) under forest management scenarios (rotational and continuous cover forestry) during the 21st century under RCP scenario 4.5 in southwestern Finland and Lapland. Each boxplot shows the variation in the CH₄ flux created by yearly variation and differences in the three climate models. The first column (2010–2019) shows the situation in a mature forest before the beginning of the different management scenarios. The next two columns (2020–2024 and 2025–2029) show the situation during the first decade of the different management scenarios, and the last four columns (2040–2049, 2060–2069, 2080–2089 and 2090–2099) show the development of CH₄ flux during the later decades.

3.3 The impact of management practices and climate on CH₄ fluxes

3.3.1 Restoration

Before the first harvest in 2020, the mature forest was a small sink of CH₄ with an average annual sink of 2 kg (C) ha⁻¹ in southwestern Finland and 1.8 kg (C) ha⁻¹ in Lapland (Fig. 5). Following the first harvest in 2020, the restoration option turned into a strong source of CH₄ (20–45 kg (C) ha⁻¹ a⁻¹). The emissions in southwestern Finland were on average twice as large as in Lapland during the first decade after restoration. After the first decade, emissions in southern Finland decreased to the same level as emissions in Lapland.

During the 21st century, there was a strong negative trend in CH₄ emissions, leading to smaller emissions by 2100 (Fig. 6 and Table A7). The slope of the trend was largest in eastern Finland and smallest in Lapland, causing Lapland to have higher emissions in later decades compared to other regions (Fig. 5).

Figure 6 Regional variation in the CH₄ flux over Finland under RCP4.5. Panel (a) shows the mean annual flux averaged over three climate models in the three management options. Panel (b) shows the average Sen's slope, describing a linear trend in a time series in restored peatlands and in the continuous cover forestry management option. The dashed lines indicate how many of the three climate models estimated a statistically significant (p<0.05) trend. No trend was calculated for rotational forestry, as the clearcutting in 2090 would disrupt any linearity in the trend.

3.3.2 Continuous cover forestry

In continuous cover forestry, the start of harvesting in 2020 had very little impact on the CH₄ sink, which remained approximately −1 to −2 kg (C) ha⁻¹ a⁻¹ (Fig. 5). The sink was stronger in southern Finland. In later decades, there was a slight positive trend in the CH₄ in the northern and eastern parts of Finland under RCP4.5, indicating a weakening CH₄ sink (Fig. 6). However, there was also considerable variability between the climate models (Table A8).

3.3.3 Rotational forestry

Rotational forestry had a stronger impact on the CH₄ fluxes compared to continuous cover forestry, with the average sink weakening to 1 kg ha⁻¹ a⁻¹ in both southwestern Finland and Lapland under RCP scenario 4.5 (Fig. 5). In the first decades after clearcutting, there was little regional variation in the CH₄ sink, but after 2030, the sink decreased in southern Finland, while there was little change in Lapland throughout the century (Fig. 6a).

After the clearcut in 2020, the soil was on average a CH₄ sink despite occasional emissions during summer (Fig. 7). However, after the clearcut in 2090, the sporadic emissions increased, causing the soil to turn into a source of methane for several years under RCP4.5 in most parts of the country (Fig. 5). The number of days with emissions during summer increased the most in the southern and central parts of the country, where over one-third of the days during summer months had emissions compared to fewer than 10 d after the 2020 harvest. In Lapland, the increase in emission days was smallest, and on average, the soil stayed a CH₄ sink. Soil temperatures were on average 1–2° warmer in 2090 than in 2020, while WTL was only slightly higher, by 1–2 cm. Anoxic conditions in the fast-decaying upper carbon pools increased considerably after the 2090 clearcut, likely due to the slightly higher WTL and increases in the fast-decaying pool.

Figure 7 Differences in CH₄ emissions and their environmental controls following the clearcuts in 2020 and 2090 under RCP4.5. Panel (a) shows the number of days during a year with any CH₄ emissions; panel (b) shows mean temperatures during June, July and August; panel (c) shows the number of days when any part of the fast-decaying soil carbon pools had anoxic conditions; and panel (d) shows the mean summer WTL.

3.4 The impact of climate scenarios

In the forestry-drained peatlands, there were stronger regional differences under RCP2.6 compared to RCP4.5, with a stronger CH₄ sink in both management options in southern Finland compared to northern Finland (Fig. A1). Similarly, regional differences in WTL were stronger under RCP2.6 (Fig. A2). In southern Finland, the average WTL stayed below −30 cm throughout the century, while the average WTL in northern Finland was approximately −20 cm. Thus, in both management options, the cumulative sinks calculated over 2020–2100 were stronger under RCP2.6 in all counties except Lapland, where the sink was weaker (Table 1). Under both management options, the cumulative sink was approximately 30–40 kg (C) ha⁻¹ stronger in the continuous cover forestry option compared to rotational forestry.

Table 1 Cumulative CH₄ fluxes between 2020 and 2100 averaged over the climate models.

In the restoration option, CH₄ emissions during the first decades after restoration were lower under RCP2.6, particularly in Lapland, where emissions were approximately half those compared to RCP4.5 (Fig. A1). Linear trends were stronger in southern and western Finland under RCP2.6 but weaker in eastern and northern Finland (Table A7). Thus, the cumulative emissions were slightly higher in southern and western Finland under RCP2.6 and lower in eastern and northern Finland (Table 1).

4 Discussion

4.1 Management impacts on CH₄ fluxes

According to our simulations, CH₄ fluxes on drained boreal peatland forests vary depending on the management practices. Considering harvesting practices, clearcutting resulted in a considerable rise in WTL and a consequent weakening of the CH₄ sink with occasional emissions as well. In previous empirical studies, the forest CH₄ sink decreased following a clearcutting, and some sites turned into sources of CH₄ during the first years following a clearcutting (Wu et al., 2011; Korkiakoski et al., 2019). However, some studies have also found no considerable changes in CH₄ fluxes following a clearcutting (Saari et al., 2009). Our results indicate that the CH₄ sink may stay weakened following the clearcut for over a decade, which is in line with observations made by Wu et al. (2011), and highlight the considerable long-term effects of harvesting practices on peatlands. In comparison, selection harvesting had only minor effects on the WTL and CH₄ sink. While this is in line with previous studies (Sundqvist et al., 2014; Korkiakoski et al., 2020), it should be noted that ditches and logging trails may still be sources of CH₄, weakening the CH₄ uptake as well after selection harvests (Korkiakoski et al., 2020). In both harvesting options, variation in the CH₄ flux between southern and northern Finland was small and seemed to follow variations in WTL, indicating that fluctuations in hydrological conditions due to climatic variability may control forest CH₄ fluxes (Ojanen et al., 2010).

Following restoration, the simulated peatlands turned to considerable sources of CH₄. This increase in CH₄ emissions has been associated with higher WTL, as well as with the recovery of peatland vegetation capable of CH₄ transport directly to the atmosphere from the anoxic soil layers (Putkinen et al., 2018; Urbanová and Bárta, 2020). Contrary to the forested peatlands, there were considerable differences in emissions between southern and northern Finland, particularly during the first decade after restoration when CH₄ emissions in southern Finland were nearly twice as large as in northern Finland. This difference was likely due to soil temperatures that were nearly 5 °C higher in southern Finland, which control CH₄ production in particular (van Hulzen et al., 1999; Turetsky et al., 2014).

4.2 Long-term trends and the impact of climate

In restored peatlands, the simulated CH₄ emissions decreased under both climate scenarios towards the end of the century. This trend was likely due to the decreasing fast-decaying carbon pool in which a large majority of the CH₄ emissions are produced in the model. This pool was high immediately after restoration due to harvest residue and started to decrease with the decomposition of the residue. The decreases in both the fast-decaying carbon pools and CH₄ emissions were particularly strong in southern Finland, where decomposition of organic material was higher due to higher temperatures. Furthermore, under RCP2.6, decomposition of the fast-decaying pools was slower due to lower temperatures throughout the century, causing the net emissions to be slightly higher in the southern and western parts of the country despite lower temperatures typically decreasing CH₄ emissions (Lai, 2009). Previous research suggests that particularly nutrient-rich peatlands may have strong CH₄ emissions immediately after restoration, followed by a decrease over time (Wilson et al., 2016). However, the decay of fast-decaying pools is possibly overestimated in the simulations and may lead to an underestimation of CH₄ emissions in the latter part of the century. Still, long-term monitoring studies following restoration of drained peatland forests are still scarce, the majority of the studies having been done in cut-over peat extraction areas (e.g. Wilson et al., 2016). Peatland plant species also differ in their CH₄ transport rates and efficiency, as well as their contribution to CH₄ formation through substrate production (Dorodnikov et al., 2011; Ge et al., 2023). Consequently, gradual changes in vegetation composition following restoration may have significant impacts on the ecosystem CH₄ flux, although this was not considered in our simulations. Nonetheless, the impact of warming on CH₄ fluxes has been shown to be conflicting in previous studies and is strongly dependent on peatland hydrology and its development (Turetsky et al., 2008; Peltoniemi et al., 2016). Our results thus highlight the need for long-term research on restored peatland forests in order to better understand the impacts of restoration and warming.

In the continuous cover forestry option, variations in the CH₄ sinks seemed to follow trends in WTL, with a weakening sink when WTL increased and a strengthening sink in areas where WTL decreased. This was particularly noticeable under the RCP2.6 scenario, where regional differences in WTL and consequently in CH₄ fluxes were stronger. However, it is important to notice that the trends in WTL were quite uncertain and depended largely on the climate model, which can be due to uncertainties in future precipitation patterns (Ruosteenoja and Jylhä, 2021). A previous study by Gong et al. (2012) found that the WTL decreased towards the end of the century in drained peatlands but also found that changes in drained peatlands were mostly smaller than in pristine peatlands. Understanding the changes in precipitation and their further impacts on peatland hydrology is essential to accurately project future greenhouse gas (GHG) fluxes and their climatic impacts.

In the rotational forestry option, the CH₄ sink started to strengthen a decade after the clearcut in most parts of Finland aside from Lapland, which is likely due to a lowering of the WTL following forest regrowth and increased evapotranspiration. In 2090, the increased emissions following the second clearcutting were likely due to temperature increases together with increased anoxic conditions in the upper peat soil. The increased anoxic conditions resulted from a slight increase in WTL and fast-decaying carbon pools, resulting in a thicker anoxic layer compared to the oxic layer. However, it should also be noted that the variation in CH₄ fluxes following the clearcutting was considerably larger than in the previous decades. Under RCP2.6, rotational forestry peatlands stayed as sinks despite similar increases in WTL, possibly due to lower soil temperatures.

4.3 Methodological limitations

The recovery of vegetation after forest management practices has been kept simplified here and may affect some of the results. In the continuous cover forestry option, the harvest removes 50 % of the total stand biomass of the forest rather than removing specific trees. In practice this means that the harvesting decreases the total green and woody biomass of the forest, but the forest growth resembles that of a mature forest rather than a mixture of young and old forest, which then leads to a rapid recovery of the forest after harvesting. This means that the impacts of selection harvesting on WTL and consequently CH₄ fluxes may have been underestimated. In the rotational forestry option, the recovery of vegetation is slow after the harvest, as the simulation of pioneer species and natural tree and shrub regeneration are not included in JSBACH. This leads to a prolonged period of very high WTL but might also underestimate the CH₄ transport through peatland species such as Eriophorum vaginatum, which might be considerable following a clearcutting (Hamberg et al., 2019). Furthermore, the simulation did not include the thinnings that are also commonly done in Finland in rotational forestry. These would likely have mediated the drop in WTL between the clearcuts and increased harvest residue in the soil, possibly slowing down the strengthening of the CH₄ sink, particularly in southern and central Finland in between clearcuts. In restored peatlands, the recovery of pristine vegetation may take a longer time after restoration compared to the very rapid recovery simulated by JSBACH.

We ran the simulations for each Finnish county to study the effect of climatic variation while keeping the required computational power at a reasonable level. While most geographic and climatic variation was adequately represented by this approach, northernmost Finland was largely represented by one single county, which means that some of the variation in northern boreal regions is likely missing. However, there are very few forestry-drained peatlands in the northernmost regions in Finland, particularly in northernmost Lapland, and official reporting such as the national forest inventory is often published at the county level (Natural Resources Institute Finland, 2023).

The functioning of peatland ecosystems and their methane fluxes can vary considerably depending on local environmental and climatic conditions such as microtopography, nutrients, hydrological conditions and past use (e.g. Lai, 2009). Thus, country-wide simulations of hypothetical (i.e. not connected to real-life sites) nutrient-rich peatlands with deep peat layers are unlikely to represent realistic conditions at any single site. Rather, this approach allows us to estimate and compare possible impacts of peatland management and climate and their feedbacks throughout a wide geographic range.

Finally, it is important to consider the limitations of complex process-based models and the reliability of their results. Due to the complex processes and limited calibration data available, variations in both input data and model structure can cause considerable deviations in the results. One important factor in limiting the uncertainty is calibration of the relevant model parameters. In this paper, we have calibrated the forestry-related parameters against forest measurements in Nordic countries (Appendix B), while CH₄ parameters have been previously estimated by Raivonen et al. (2017) and Li et al. (2024). The three climate models and two RCPs show uncertainty related to climate, which, according to Mäkelä et al. (2019), is the most considerable source of uncertainty in carbon balance simulations in land surface models. Another important factor in limiting model uncertainty is model comparison studies. While this was outside the scope of the paper, further simulation studies of CH₄ fluxes on drained peatlands are necessary to further improve our understanding of their current and future changes.

5 Conclusions

We simulated the impact of management practices and climate on CH₄ fluxes from forestry-drained boreal peatlands using the land surface model JSBACH. Our simulations showed that restoration turned peatlands into sources of CH₄, but the magnitude of emissions varied regionally, with larger emissions in southern Finland than in northern Finland. Furthermore, emissions decreased towards the end of the century as harvest residue diminished from the upper peat layers. In forested peatlands, clearcutting had a stronger weakening effect on the forest CH₄ sink compared to selection harvesting, and the effect was stronger towards the end of the century under RCP4.5. Water-table level was found to have a strong control on the CH₄ fluxes, particularly on forested peatlands.

Appendix A

Table A1 YASSO parameters for JSBACH-FOM. Decomposition parameters are given for aboveground (AG) and belowground (BG) AWEN (acid-hydrolysable, water-soluble, ethanol-soluble, and neither hydrolysable nor soluble) carbon pools. The parameters are the reference decomposition rate, anoxic decomposition modifier, peat decomposition modifier, anoxic peat decomposition factor, and anoxic and oxic decomposition rates. The proportion of soil carbon from the N to A pool and the relationship between AWEN to H (humus) and N to H are also given.

Table A2 YASSO parameters for the peatland version of JSBACH. Decomposition parameters are given for aboveground (AG) and belowground (BG) AWEN (acid-hydrolysable, water-soluble, ethanol-soluble, and neither hydrolysable nor soluble) carbon pools. The proportion of soil carbon from the N to A pool and the relationship between AWEN to H (humus) and N to H are also given.

Table A3 HIMMELI model parameters in JSBACH-FOM and JSBACH-PEAT.

Table A4 Simulated total soil carbon pool (m⁻²) on peatlands in each region at the end of the spinup in 1950.

Table A5 Simulated and measured CH₄ fluxes in restored peatlands depending on soil temperature.

Table A6 Simulated and measured CH₄ in forestry-drained peatlands depending on water-table depth (WTD).

Table A7 Linear trends and their uncertainty related to climate models of CH₄ fluxes, soil respiration, summer WTL and summer soil temperature (ST) in the Finnish counties in restored peatlands. The values show the Sen's slope parameter averaged over the climate models and their variability. Bolded values indicate significant trends (p<0.05) estimated by the Mann–Kendall trend test. For the numbered list of Finnish counties, see Fig. 2.

¹ Indicates that only one model showed a significant trend; ² indicates that two models showed a significant trend.

Table A8 Linear trends and their uncertainty related to climate models of CH₄ fluxes, soil respiration, summer WTL and summer soil temperature (ST) in the Finnish counties in continuous cover forestry peatlands. The values show the Sen's slope parameter averaged over the climate models and their variability. Bolded values indicate significant trends (p<0.05) estimated by the Mann–Kendall trend test. For the numbered list of Finnish counties, see Fig. 2. ¹ Indicates that only one model showed a significant trend; ² indicates that two models showed a significant trend.

Figure A1 Average annual CH₄ flux in (a) restored peatlands and (b) under forest management scenarios (rotational and continuous cover forestry) during the 21st century under RCP scenario 2.6 in southwestern Finland and Lapland. Each boxplot shows the variation in the CH₄ flux created by yearly variation and differences in the three climate models. The first column (2010–2019) shows the situation in a mature forest before the beginning of the different management scenarios. The next two columns (2020–2024 and 2025–2029) show the situation during the first decade of the different management scenarios, and the last four columns (2040–2049, 2060–2069, 2080–2089 and 2090–2099) show the development of the CH₄ flux during the later decades.

Figure A2 The impact of forest management on summer WTL, soil temperature (−22 cm below ground) and fast-decaying carbon pools under RCP scenario 4.5. Each boxplot shows the variation created by yearly variation and differences in the three climate models. The first column (2010–2019) shows the situation in a mature forest before the beginning of the different management scenarios. The next two columns (2020–2024 and 2025–2029) show the situation during the first decade of the different management scenarios, and the last four columns (2040–2049, 2060–2069, 2080–2089 and 2090–2099) show the development of the variables during the later decades.

Appendix B: JSBACH-FOM description

B1 Model description

JSBACH version 4 has recently been developed, which can account for a distribution of forest age classes within the forest PFTs (Nabel et al., 2020). The model we used, JSBACH-FOM, also includes the age for forest PFTs, but it was implemented in JSBACH version 3. In addition to the ability to account for the age of forest PFTs, the model differs from earlier versions in the way the maximum LAI is treated. The maximum LAI that can be reached during the growing season has previously been PFT dependent but constant. In JSBACH-FOM, the maximum LAI is dependent on the available leaf biomass needed to simulate forest regrowth and age.

In our work, the phenology was calculated using the Logistic Growth Phenology (LoGro-P) model (Böttcher et al., 2016). The LoGro-P model provides the seasonal development of LAI for each PFT based on the temperature and soil moisture. The LAI is further limited by a maximum LAI, which is either prescribed for each PFT or, as in our case, dependent on the available leaf carbon for the forest PFTs. Other PFT-dependent parameters are the phenology type (e.g. evergreen and grass) and specific leaf area (SLA).

The calculation of the growing forest starts from the total vegetation carbon. The number of trees per area, or stem number, is calculated from the total vegetation carbon, assuming that the forests are in a self-thinning state. The biomass per individual can then be used together with allometric relationships to derive the maximum LAI of the forest. The implementation of the forest growth is based on maximum LAI. The setup of the model parameters for our simulations is explained in detail in Sect. B2. The forest ageing is done at the beginning of each year. The harvest is done for each tile when the rotation time is reached. Once the criterion is reached, the forest is clearcut at the beginning of the year, and the age and size of the trees restart from zero. The carbon pools are redistributed due to harvesting. The harvested carbon, making up 77 % of the aboveground woody pool, is removed from the calculations. Furthermore, all green and reserve carbon, both above- and belowground, go into the YASSO AWEN litter pools according to predefined fractions. The belowground woody carbon is distributed into the belowground litter pools. In addition, a fraction of the aboveground woody carbon, i.e. the slash fraction, is distributed into the aboveground litter pools. The slash fraction equals 23 % of the aboveground woody carbon pool and accounts for the damage and small woody fragments left in the forest during harvest.

Briggs and Cunia (1982)Briggs and Cunia (1982)Mäkelä and Vanninen (1998)Mäkelä and Vanninen (1998)Marklund (1988)Mälkönen (1974)Drexhage and Gruber (1999)Marklund (1988)Briggs and Cunia (1982)Briggs and Cunia (1982)Hakkila (1991)Mäkelä and Vanninen (1998)Mäkelä and Vanninen (1998)Marklund (1988)

Table B1 References for biomass equations that were considered in this publication (biomass types are AG – aboveground; ST – stem; CR – crown; BG – belowground; and FL – foliage). Aboveground biomass is composed of stem and crown biomass.

Figure B1 Biomass per area as a function of stem number. The lines have been cutoff based on the valid range of dbh; in addition the model has a maximum stem number defined as ln (N)_max = 9.5. Marklund (1988) is shown as a cyan line and modified Marklund as a dashed blue line, with parameters α_{nr_ind} = 15 and β_{nr_ind} = −0.39. The grey lines are calculated using four allometric relationships of aboveground biomass and two of belowground biomass, giving in total eight estimates for the relationship between total biomass per area and stem number.

B2 Model setup

The standard setup of land surface parameters for extratropical coniferous forest was used, with some modifications, to represent Scots pine. The specific leaf area (SLA) was set to 61.62 cm² g⁻¹, according to Goude et al. (2019). Competition for resources between trees in a stand results in self-thinning. The pine forest in the model is assumed to be in a self-thinning state. The parameters needed for describing the forest growth were derived from published expressions for self-thinning and allometric relationships for Scots pine. Self-thinning of evenly aged stands can be modelled using a relationship between the quadratic mean diameter and the number of trees per area (Eq. B1) according to Reineke (1933).

$$\begin{array}{}\text{(B1)}& \ln \left(N\right)=p-q\cdot \ln \left(D_g\right),\end{array}$$

where N is the stem number (ha⁻¹) and D_g is the mean diameter at breast height (dbh) in centimetres weighted with the basal area, i.e. the quadratic mean diameter. Each tile represents only one age class of trees, and the model has only one diameter per age class, but the quadratic mean equals the arithmetic mean when the variance is zero. The intercept (p) and slope (q) of the log–log relationship were obtained from Hynynen (1993): 12.669 and −1.844, respectively. Hynynen (1993) used data for Scots pine from 19 unthinned, even-aged and monospecific plots in Finland to derive the parameters p and q. Only plots where no extensive natural disturbance had occurred during the study period were included. The data were collected between 1924 and 1989, on average six measurements over 38 years. In JSBACH-FOM, the self-thinning expression (Eq. B2) relates the stem number and the biomass per unit area.

$$\begin{array}{}\text{(B2)}& \ln \left({\mathrm{BM}}_{\mathrm{veg}}\right)={\mathit{\alpha }}_{\mathrm{nr}\mathrm{\_}\mathrm{ind}}+{\mathit{\beta }}_{\mathrm{nr}\mathrm{\_}\mathrm{ind}}\cdot \ln \left(N\right),\end{array}$$

where BM_veg is the vegetation carbon in the maximum green pool (kg m⁻²). In order to derive the slope (β_{nr_ind}) and intercept (α_{nr_ind}), allometric relationships between the dbh and the biomass of various above- (AB) and belowground (BG) components are needed.

$$\begin{array}{}\text{(B3)}& {\mathrm{BM}}_{\mathrm{ind}}=\mathrm{AG}+\mathrm{BG}\end{array}$$

The biomass per area is obtained from the biomass of a single tree, BM_ind (kg m⁻²), and the stem number.

$$\begin{array}{}\text{(B4)}& {\mathrm{BM}}_{\mathrm{veg}}={\mathrm{BM}}_{\mathrm{ind}}\cdot N/\mathrm{10}\mathrm{000}\end{array}$$

We considered allometric relationships for pine based on data from Finnish sites compiled by Zianis et al. (2005). A summary is given in Table B1. First the total biomass is obtained as a function of dbh by summing the AG and BG biomass according to Eq. (B3). When both the stem number and the biomass are expressed as a function of dbh, the coefficients α_{nr_ind} and β_{nr_ind} in Eq. (B2) can be obtained by plotting ln (BM_veg) against ln(N). We used four allometric relationships for the AB and two for the BG biomass, which gave eight different relationships between biomass per area (BM_veg) and stem number (N). These are plotted in Fig. B1, grey lines. In addition we used allometric relationships by Marklund (1988) based on data from Sweden to derive one more relationship, the cyan line in Fig. B1.

The relationships are plotted only for the valid range of the original allometric relationships (given in Table B1). There is an upper limit for the stem number in the plot, i.e. ln (N)=9.5, based on a cut-off value in JSBACH-FOM, which is used to prevent an excess number of very small trees. From Fig. B1 it can be seen that the biomass based on Marklund (1988) fits within the range of the ones derived from data from Finnish plots. We selected the parameter values based on Marklund (1988) for the simulations due to the fact that the equations are based a large sample size, hundreds of felled trees from a large geographical area. In addition, the valid range of dbh (0–45 cm) is larger than for the other allometric relationships. Marklund (1988) also provides relationships for both above- and belowground biomass. In the final simulations we used a modified relationship, the dashed blue line in Fig. B1.

Figure B2 Relationships describing the leaf carbon as a function of total biomass. Marklund (1988) is shown as a cyan line and modified Marklund as a dashed blue line with parameters α_leaf = −2.0 and β_leaf = 0.71.

Kolari et al. (2009)Palmroth and Hari (2001)Ilvesniemi et al. (2009)Kolari et al. (2009)Ilvesniemi et al. (2009)Kolari et al. (2009)Kolari et al. (2009)Kolari et al. (2009)

Table B2 Summary of the comparison of model results with observations for Hyytiälä and Sodankylä. Hyytiälä CO₂ fluxes are the average of 2002–2007. All-sided LAI was divided by 2.57 to obtain one-sided LAI. GPP refers to gross primary production, NEE to net ecosystem exchange and R_e to ecosystem respiration.

JSBACH-FOM also requires the relationship of leaf carbon (C_leaf) and the biomass of a tree (BM_ind).

$$\begin{array}{}\text{(B5)}& \ln \left({\mathrm{C}}_{\mathrm{leaf}}\right)={\mathit{\alpha }}_{\mathrm{leaf}}+{\mathit{\beta }}_{\mathrm{leaf}}\cdot \ln \left({\mathrm{BM}}_{\mathrm{ind}}\right)\end{array}$$

Again we used the allometric relationships for AG and BG biomass in Table B1 to derive the total biomass of a tree, which gave eight descriptions for the tree biomass (BM_ind). In addition we had five expressions describing the foliage biomass as a function of dbh based on data from Finnish sites (Table B1). The foliage biomass was divided by 2 to obtain C_leaf. These expressions were used to derive 40 different relationships between leaf carbon (C_leaf) and tree biomass (BM_ind) according to Eq. (B5). These are plotted in Fig. B2 as grey lines. The relationships are again plotted only for the valid range of the original allometric relationships (given in Table B1). We also used allometric relationships by Marklund (1988) based on data from Sweden to derive one set of α_leaf and β_leaf, the cyan line in Fig. B2. We decided to use these in JSBACH-FOM due to the wide range of dbh where the relationships are valid. According to Fig. B2, the relationship from the Marklund (1988) data fits within the range of the ones derived from data from Finnish plots. Based on test simulations, we modified the Marklund (1988) coefficients for the final simulations, the dashed blue line in Fig. B2. In the model, the maximum LAI is calculated from leaf carbon per tree using the SLA. The increase in the maximum LAI is stopped when then the increase in maximum LAI per year is less than 1 %.

The parameters derived from the Marklund (1988) biomass equations for pine described in Sect. B2 were adjusted manually in order to get a better agreement between simulated and observed GPP. The relationships are plotted in Figs. B1 and B2, and relevant model results compared to observations are shown in Table B2.