A new method for estimating carbon dioxide emissions from drained peatland forest soils for the greenhouse gas inventory of Finland

In peatlands drained for forestry, the soil carbon (C) or carbon dioxide (CO₂) balance is affected by both (i) higher heterotrophic CO₂-C release from faster decomposing soil organic matter (SOM) and (ii) higher plant litter C input from more vigorously growing forests. This balance and other greenhouse gas (GHG) sinks and sources in managed lands are annually reported by national GHG inventories to the United Nations Climate Change Convention. In this paper, we present a revised, fully dynamic method for reporting the CO₂ balance of drained peatland forest soils in Finland. Our method can follow temporal changes in tree biomass growth, tree harvesting and climatic parameters, and it is built on empirical regression models of SOM decomposition and litter input in drained peatland forests. All major components of aboveground and belowground litter input from ground vegetation as well as live trees and trees that died naturally are included, supplemented by newly acquired turnover rates of woody plant fine roots. Annual litter input from harvesting residues is calculated using national statistics of logging and energy use of trees. Leaching, which also exports dissolved C from drained peatlands, is not included. The results are reported as time series from 1990–2021 following the practice in the GHG inventory. Our revised method produces an increasing trend of annual emissions from 0.2 to 2.1 t CO₂ ha⁻¹ yr⁻¹ for the period 1990–2021 in Finland (equal to a trend from 1.4 to 7.9 Mt CO₂ yr⁻¹ for the entire 4.3 Mha of drained peatland forests), with a statistically significant difference between the years 1990 and 2021. Across the period 1990–2021, annual emissions are on average 1.5 t CO₂ ha⁻¹ yr⁻¹ (3.4 Mt CO₂ yr⁻¹ for 2.2 Mha area) in warmer southern Finland and −0.14 t CO₂ ha⁻¹ yr⁻¹ (−0.3 Mt CO₂ yr⁻¹ for 2.1 Mha area) in cooler northern Finland. When combined with data on the CO₂ sink created by the growing tree stock, in 2021 the drained peatland forest ecosystems were a source of 1.0 t CO₂ ha⁻¹ yr⁻¹ (2.3 Mt CO₂ yr⁻¹) in southern Finland and a sink of 1.2 t CO₂ ha⁻¹ yr⁻¹ (2.5 Mt CO₂ yr⁻¹) in northern Finland. We compare these results to those produced by the semi-dynamic method used earlier in the Finnish GHG inventory and discuss the strengths and vulnerabilities of the new revised method in comparison to more static emission factors.

1 Introduction

Drainage of peatlands improves forest growth but also turns unmanaged land to managed land and creates anthropogenic greenhouse gas (GHG) emissions: drainage increases aerobic decomposition of soil organic matter (SOM) and leads to substantial carbon dioxide (CO₂) emissions, globally accounting for 6 % of total CO₂ emissions (Joosten, 2010). Parties and signatories of the Paris Agreement and the United Nations Framework Convention on Climate Change (UNFCCC) are committed to reporting managed land GHG annual emissions and removals by sources and sinks, respectively, from 1990 onwards. Emissions from drained peatlands are reported under the land use, land-use change and forestry (LULUCF) sector of GHG inventories.

In Finland peatlands drained for forestry comprise an important category of managed lands. Almost one-third, or 8.8 Mha (million hectares), of land area in Finland is covered by peatlands. This area also includes ca. 1 Mha thin-peated sites (< 30 cm peat layer). Of the total peatland area, 4.7 Mha is drained for forestry (Korhonen et al., 2021), 0.34 Mha for agriculture and 0.11 Mha for peat extraction (Statistics Finland, 2022). Peatland forests were mainly drained after World War II, with a peak in activity in the 1960s and 1970s (Päivänen and Hånell, 2012). This paper addresses soil CO₂ emissions from those 4.3 Mha of drained peatland forests in Finland that fulfill the FAO definition of forest land. Soil CO₂ emissions from undrained peatland forests are excluded since the water table level, and thereby the functioning of soil, has not been modified by human activity. Emissions of CO₂ from drained peatland forest soils were first globally estimated by IPCC (2006). In Finland, the emissions, or soil C balance, have been quantified by Minkkinen et al. (2002) and Ojanen et al. (2014).

Countries may use static IPCC default emission factors for reporting soil CO₂ emissions from their drained peatland forest soils. Country-specific data and improvement of GHG inventory methods beyond the IPCC default emission factors are relevant, especially for countries that have large areas of drained peatlands. Some countries have indeed opted to use country-specific emission factors; these can be based on the subsidence, i.e., fall of the surface of the organic soil, as in the Netherlands (Ruyssenaars et al., 2021), CO₂ flux measurements as in Ireland (Duffy et al., 2021) or soil C balance modeling as in the UK (Brown et al., 2021). The weakness of static emission factors is that they cannot follow temporal changes, if such changes occur, in variables that drive the soil CO₂ balance. For instance, net CO₂ emissions from peatland forest soils may gradually evolve after drainage due to a successional transition of vegetation and changes in peat decomposition and hydrology (Sarkkola et al., 2010; Straková et al., 2012). In addition, increasing global temperatures likely increase SOM decomposition rates and thus soil CO₂ emissions.

Here we present a new dynamic method for estimating the soil CO₂ balance of drained peatland forests in Finland on the timescales needed in the GHG inventory. We use both measured and modeled values of peat and litter decomposition as well as plant litter input. While water table depth (WTD) is the main factor that controls C loss from peat by heterotrophic respiration in drained wetlands (Silvola et al., 1996; Ojanen et al., 2010; Jauhiainen et al., 2019), direct data on WTD are not available nationwide for GHG inventory purposes in Finland. We therefore link decomposition of peat and litter to tree basal area per hectare (BA) using empirical regression models available for different drained peatland forest site types (Ojanen et al., 2014). BA provides a proxy for the rate of evapotranspiration, which largely controls WTD in forestry-drained peatlands (Hökkä et al., 2021; Leppä et al., 2020). In addition, BA can be used to predict litter input from trees and ground vegetation (Ojanen et al., 2014). The time series of soil CO₂ balance is then calculated using a time series of BA provided by the Finnish National Forest Inventory, NFI (Tomppo et al., 2011; Korhonen et al., 2021), combined with a time series of air temperature. Since drained peatland forests are not evenly distributed across Finland, the temperature-dependent peat and litter decomposition is modeled using spatially matched environmental variables. With this method, soil CO₂ balance can be estimated as a difference of C added to soil via litter input and CO₂ lost through SOM decomposition (Ojanen et al., 2012).

The CO₂ balance calculation method that was applied earlier in the Finnish GHG inventory for drained peatland forest soils also had dynamic components, but only in terms of belowground litter input. We call it the “semi-dynamic” or “previous” method in this paper. The estimate of SOM decomposition was static and spatially poorly representative as it relied on data collected from a limited number of field sites (Minkkinen et al., 2007). Aboveground litter input was not included in the previous calculations. Producing a revised, fully dynamic method became possible when new data and empirical models appeared that were able to link together the drivers and components of soil CO₂ balance (Ojanen et al., 2010, 2012, 2013, 2014). To illustrate the consequences of adopting the new method for the GHG inventory of Finland, we also present, for comparison, the CO₂ balance of drained peatland forest soils as estimated by the earlier semi-dynamic method (Statistics Finland, 2022).

2 Materials and methods

2.1 The concept of soil CO₂ balance

The estimation of soil CO₂ balance (CO_2 Net) is based on the annual difference between the release of CO₂ from decomposing soil organic matter (SOM), or heterotrophic soil respiration R_Het, and the C entering the soil through plant litter input, and it can be expressed as

$$\begin{array}{}\text{(1)}& {\mathrm{CO}}_{\mathrm{2}\mathrm{Net}}=R_{\mathrm{Het}}-{\displaystyle \frac{\mathrm{44}}{\mathrm{12}}}\left(I_{\mathrm{AGL}}+I_{\mathrm{BGL}}+I_{\mathrm{AGR}}+I_{\mathrm{BGR}}\right),\end{array}$$

where I_AGL and I_BGL are the annual C mass input (t or Mt of C, as appropriate in the context) of aboveground and belowground litter, respectively, from living trees and ground vegetation; I_AGR and I_BGR are the annual C input of aboveground and belowground residues, respectively, from forest harvests and naturally dying trees; and $\mathrm{44}/\mathrm{12}$ is the ratio of CO₂ and C molar masses, converting the organic C in dead organic matter to the respective mass unit of CO₂. Positive results in Eq. (1) appear when decomposition R_Het is larger than the input of litter and residues and denote net CO₂ emission from the soil to the atmosphere. Negative values denote a net increase in soil C stock. The time dimension of CO₂ balance is 1 year, and decomposition therefore includes wintertime decay as well.

R_Het can further be defined as the sum of peat decomposition and decomposition of aboveground and belowground plant litter as well as residues from harvests and trees that died naturally:

$$\begin{array}{}\text{(2)}& R_{\mathrm{Het}}=D_{\mathrm{P}}+D_{\mathrm{AGL}}+D_{\mathrm{BGL}}+D_{\mathrm{AGR}}+D_{\mathrm{BGR}},\end{array}$$

where D_P is the annual CO₂ release (t or Mt of CO₂, as appropriate in the context) from decomposing peat, D_AGL and D_BGL are the annual CO₂ release from decomposing aboveground and belowground litter, respectively, of living trees and ground vegetation, and D_AGR and D_BGR are the annual CO₂ release from decomposing aboveground and belowground harvest and natural mortality residues, respectively (Fig. 1). Autotrophic (root) respiration is not a part of the soil C stock change and is therefore not included in R_Het. To focus on CO₂ release from SOM, not confounded by autotrophic root respiration, the field plots were trenched (Ojanen et al., 2013).

Figure 1 Sources of data and workflow for the new method to produce the time series of annual soil CO₂ balance for drained organic forest land in the GHG inventory of Finland. NFI means the National Forest Inventory, FMI is the Finnish Meteorological Institute and “Luke” is the Natural Resources Institute Finland. Characteristics of the method previously applied in the GHG inventory, in comparison to the new method, are indicated in grey boxes. References to key figures, tables and equations are shown in the respective boxes. T: temperature, P: precipitation, BCEF: biomass conversion and expansion factor, FTYPE: drained peatland forest site type.

2.2 Drained peatland forest site types

One key element in our method is that litter input and decomposition are individually estimated for each drained peatland forest site type, hereafter denoted as FTYPE. In Finland, FTYPE classification was developed for guiding drainage and forest management on peatlands. It follows the general Finnish forest site type theory (Cajander, 1913) and uses species of trees and ground vegetation as indicators of moisture and fertility regimes. This typification of drained peatlands (Laine, 1989; Laine et al., 2012) is a practical tool for forestry and forest inventory in Finland. Forest growth and the rate of SOM decomposition are related to FTYPE: more fertile FTYPEs support greater tree biomass and growth, and consequently they also have higher decomposition rates (Silvola et al., 1996; Minkkinen et al., 2007; Ojanen et al., 2010).

The current typology of FTYPEs follows the principles presented in Laine et al. (2012) and includes a total of five FTYPEs: herb-rich FTYPE, Vaccinium myrtillus FTYPE, Vaccinium vitis-idaea FTYPE, dwarf shrub FTYPE and Cladonia FTYPE, listed in order of decreasing soil fertility (Table 1). Differences in fertility between FTYPEs are evident in tree and ground vegetation species composition (Table 1). We did not separate the two subtypes of the V. myrtillus and V. vitis-idaea FTYPEs (I and II with different original mire types; Laine et al., 2012) because field determination of FTYPE subtypes can be ambiguous, especially V. vitis-idaea. Drained peatland forest sites have been classified into FTYPEs in field surveys in all NFIs, but the classification shown in Table 1 has only been employed since NFI10 (carried out in 2004–2008). The FTYPEs used in NFI8 (1986–1994) and NFI9 (1996–2003) were converted to correspond to the classification in Table 1 using field-collected data on soil type, drainage status, site fertility and vegetation.

Table 1 Characteristics of drained peatland forest site types (FTYPEs) and the GHG inventory estimates of their area and proportion of all drained peatland forest area remaining forest, based on NFI12 data and following the FAO forest classification, in southern and northern Finland (n: number of sample plot centers in each FTYPE). See Fig. 1 for the southern and northern regions of Finland and the distribution of FTYPEs across the regions.

The areas and proportions of FTYPEs of all drained peatland “forests remaining forests” (i.e., forests that have not undergone another change in land use in the past 20 years) in southern and northern Finland, derived from NFI12 (2014–2018), are shown in Table 1, and the distribution of FTYPEs across Finland is illustrated in Fig. 2 (the definition of southern and northern Finland used in this paper follows the definition used in NFI). FTYPE areas have changed over the years due to changes in land use, such as new peatland drainage or reclamation of agricultural land. To provide a time series of FTYPE areas for the GHG inventory period starting from 1990 (Fig. 3), area data were collected from six NFIs, i.e., NFI8 (1986–1994), NFI9 (1996–2003), NFI10 (2004–2008), NFI11 (2009–2013), NFI12 (2014–2018) and NFI13 (2019–2023). As for all other time series that utilize NFI data, the years between NFI mid-years are interpolated. However, the years 1990–1992 in northern Finland are extrapolated because the NFI8 mid-year was after 1990, and from 2017 onwards the values are determined with the help of the not yet complete data from NFI13. The development of the Finnish NFI until NFI9 is described by Tomppo et al. (2011) and the sampling design of NFI12 by Korhonen et al. (2021).

Figure 2 Distribution of drained peatland forest site types (FTYPEs) in southern and northern Finland, as observed in NFI12. The color represents the area (kha) of FTYPE in a 20 km × 20 km grid cell, and the division of southern and northern Finland (thick black line) follows the division used in the NFI. The northernmost parts of the country (separated with a thin black line) are not included in NFI. Rhtkg: herb-rich FTYPE, Mtkg: Vaccinium myrtillus FTYPE, Ptkg: Vaccinium vitis-idaea FTYPE, Vatkg: dwarf shrub FTYPE, Jatkg: Cladonia FTYPE, listed in order of decreasing soil fertility.

Figure 3 Time series of basal area per hectare for each FTYPE and FTYPE areas in southern and northern Finland as estimated from NFI data. Rhtkg: herb-rich FTYPE, Mtkg: Vaccinium myrtillus FTYPE, Ptkg: Vaccinium vitis-idaea FTYPE, Vatkg: dwarf shrub FTYPE, Jatkg: Cladonia FTYPE, listed in order of decreasing soil fertility.

The drained peatland forest sites included in our data are those classified as forests in the Finnish FAO Global Forest Resources Assessment (FRA) (FAO, 2020): i.e., the area is ≥ 0.5 ha, trees are able to reach a height of 5 m and the canopy cover is > 10 %. Currently, NFI field surveys include FRA land classification, and the earlier datasets from the 1990s have been reclassified by NFI.

2.3 Climatic variables

The spatial distributions of FTYPEs (Fig. 2), and thereby their climatic attributes, differ due to geographic factors that originally resulted in the formation of different mires in Finland (Ruuhijärvi, 1960; Havas, 1961; Eurola, 1962). Therefore, the mean May–October air temperatures by FTYPE (approximating the growing season mean air temperature), which are needed to predict SOM decomposition in different FTYPEs (see Table 2), are calculated using real FTYPE distributions in southern and northern Finland (Fig. 2). Briefly, the May–October mean air temperature is first produced for each drained peatland forest location by combining the location data from sample plots in NFI12 with a 10 km × 10 km weather grid data provided by the Finnish Meteorological Institute (Venäläinen et al., 2005). Using the obtained data, FTYPE mean values are calculated for southern and northern Finland. For simplicity, NFI12 locations are used to calculate the FTYPE mean temperatures for all years in the 1990–2021 time series (i.e., we assume that FTYPE spatial distributions have not significantly changed in time). Lastly, the time series for both regions of the country are smoothed using 30-year rolling means (Fig. 4). Smoothed long-term means follow long-term trends, such as gradual warming, but do not bring annual fluctuation in the results and thus help in revealing the anthropogenic impact of changes in drainage area, forest development and harvests on managed land emissions.

Table 2 Empirical regression models of CO₂ release from peat and litter decomposition for FTYPEs as used in the new method. The regression models are from Ojanen et al. (2014), except that the constants for the V. myrtillus and V. vitis-idaea FTYPEs are weighted means of constants of the two subtypes; BA: tree stand basal area (m² ha⁻¹), T: mean May–October air temperature (^∘C).

Figure 4 May–October mean air temperature for each FTYPE in southern and northern Finland and mean annual air temperature and precipitation of all drained peatland forests in southern and northern Finland for 1990–2021. In all cases, the means are 30-year rolling means. Rhtkg: herb-rich FTYPE, Mtkg: Vaccinium myrtillus FTYPE, Ptkg: Vaccinium vitis-idaea FTYPE, Vatkg: dwarf shrub FTYPE, Jatkg: Cladonia FTYPE, listed in order of decreasing soil fertility.

Decomposition, i.e., CO₂ production of litter from forest harvests and trees that died naturally, is calculated using the Yasso07 model (Tuomi et al., 2009, 2011), which uses mean annual temperature, temperature amplitude (the difference between the mean temperature of the warmest and coldest month divided by 2) and annual precipitation as environmental drivers. Spatially weighed and smoothed 1990–2021 time series for these variables are produced as for May–October mean air temperature above, but since the harvest residue data are currently allocated between southern and northern Finland, not to each FTYPE, the means are calculated for southern and northern Finland, not for each FTYPE (Fig. 4, Supplement Fig. S1).

2.4 Basal area

FTYPE mean basal area per hectare, BA (Fig. 3), is used to predict peat and litter CO₂ release (see Table 2) as well as ground vegetation and arboreal fine root litter input (Tables 3, 4). The FTYPE mean BA is estimated using data on tree stem diameter (measured at breast height, DBH) of tree species groups (pine, spruce and deciduous species), collected in the sample plots of the NFI8–NFI13 inventories. As in the case of FTYPE areas, FTYPE mean BAs are interpolated for years between the NFI mid-years, extrapolated for 1990–1992 in northern Finland and derived from the not yet complete NFI13 for the years 2017–2021 to produce the time series needed in GHG inventory (Fig. 3). BAs have generally increased with time in drained peatlands, demonstrating the benefits of drainage for tree growth (Fig. 3).

Table 3 Regression models from Ojanen et al. (2014) for estimating mean arboreal fine root (≤ 2 mm diameter) biomass (g m⁻²) for FTYPEs in southern and northern Finland. BA: mean basal area (m² ha⁻¹); decid: deciduous trees; cover_shrub: mean dwarf shrub cover (% of area), taken from Table 3.

Table 4 Regression models of litter input from ground vegetation (combined dry mass of aboveground and belowground litter excluding dwarf shrub fine root litter) in different FTYPEs. Models are from Ojanen et al. (2014), except for the constants of V. myrtillus and V. vitis-idaea FTYPEs, which are area weighted means of constants of their subtypes; BA: basal area (m² ha⁻¹).

2.5 Calculating litter input

Plant litter input consists of (1) litter from living trees (excluding fine root litter), (2) arboreal fine root litter (roots of trees and dwarf shrubs with ≤ 2 mm diameter), (3) ground vegetation litter (excluding dwarf shrub fine root litter), and (4) litter originating from forest harvests and natural tree mortality. In all calculations, all litter is assumed to contain 50 % C of dry mass.

To estimate litter input from living trees, the mean biomass (Mg ha⁻¹) of each tree component (stem wood, stem bark, live branches, dead branches, foliage, base and stump, and coarse roots of > 10 mm diameter) for each FTYPE is first obtained from NFI. These biomasses are calculated by multiplying mean stem volume in each FTYPE with biomass conversion and expansion factors (BCEFs) based on biomass models of Repola (2008, 2009). BCEFs are component biomass to stem volume ratios, estimated using a subsample of trees measured in detail in each NFI for all different combinations of tree species group, FTYPE, and south and north region. Inputs of litter from living trees, consisting of senesced foliage, branches, stem and stump bark, and coarse roots (> 10 mm diameter), are then estimated by multiplying the biomass estimates and litter production rates (Table 5) of tree components. Litter production rate tells the proportion of the component mass that turns into litter in a year.

Table 5 Litter production rates of tree biomass components (i.e., the proportion of the component mass that turns into litter in a year) as derived from studies by Lehtonen et al. (2004), Muukkonen and Lehtonen (2004), Starr et al. (2005), Liski et al. (2006), and Ojanen et al. (2014).

^S Southern Finland, ^N northern Finland, ^* calculated as a mean of pine and deciduous tree estimates as no estimate for spruce was available from the literature.

To calculate arboreal fine root litter input, fine root biomass is first estimated for each FTYPE by using FTYPE mean BAs of tree groups and FTYPE mean dwarf shrub cover (Table 6) as predictors in empirical regression models (Table 3), obtained from Ojanen et al. (2014). Trees and dwarf shrubs are combined due to difficulties in distinguishing their fine roots in field samples. The dwarf shrub cover for the Vaccinium myrtillus and Vaccinium vitis-idaea FTYPEs is calculated as a weighted mean of values in I and II subtypes with the relative areas of subtypes (61.3 % and 38.7 % for the V. myrtillus FTYPE and 60.1 % and 39.9 % for the V. vitis-idaea FTYPE, respectively) as a weight factor. The fine root biomasses include an estimate for 0–20 cm rooting depth (Ojanen et al., 2013) plus 4.3 % in deeper depths (Laiho and Finér, 1996). Arboreal fine root litter input is then estimated by multiplying fine root standing biomass by the tree fine root turnover rate measured for each FTYPE (Table 6).

Table 6 Dwarf shrub cover, obtained from Ojanen et al. (2014), and turnover rate of tree fine roots (≤ 0.5 mm diameter) in FTYPEs. The SDs are calculated from average median turnover rates.

^* Value replicated from dwarf shrub FTYPE as no respective site was measured.

Turnover rates of tree fine roots in Table 6 were estimated using minirhizotrons installed in each FTYPE (except for the Cladonia FTYPE, for which the estimate of the dwarf shrub FTYPE is used). Following the minirhizotron method (Lukac, 2012), transparent acrylic tubes (9–12 tubes per site) were installed into the soil on six sites to photograph fine roots (Kari Minkkinen et al., unpublished data). After a 1-year stabilization period, all tubes were photographed 19 times during 4 consecutive years. It has been stated that minirhizotrons may underestimate root longevity in short (< 3 years) studies since stabilization after installation disturbance may take several years (Strand et al., 2008). Our justification for a 1-year stabilization is that installation to peat soil causes little disturbance (no need to dig the soil) in comparison to installation to mineral soils, as in Strand et al. (2008). Images of 1.3 cm × 2 cm (height × width) were taken with a minirhizotron camera (BTC-2; Bartz Technology, Santa Barbara, USA) and analyzed using the WinRHIZO Tron 2015b program (Régent, Canada). The diameters, lengths and depths of roots were recorded. Root longevities (one per turnover) were determined as years with the method by Kaplan and Maier (1958) using median longevity, i.e., the lifetime when 50 % of the observed roots had died. If 50 % mortality was not reached during the monitoring period, a parametric regression model with Weibull error distribution was applied to predict it using survreg (Kalbfleisch and Prentice, 2002) in the R package Survival 3.4-0. The turnover rates are based on tree fine roots with ≤ 0.5 mm diameter, and they are applied for all arboreal fine roots of ≤ 2 mm diameter.

The mean combined aboveground and belowground litter input from ground vegetation (excluding dwarf shrub fine root litter input) is calculated for each FTYPE using FTYPE mean BA as a predictor in empirical regression models (Table 4), obtained from Ojanen et al. (2014). Notably, although dwarf shrub areal cover is assumed to remain static regardless of BA in different FTYPEs (Table 3), the regression models that predict ground vegetation litter input (Tables 4, 5) have BAs of trees as predictors and thus follow changes in BA in FTYPEs.

In forest harvests, most of the collected roundwood is utilized for wood products, but some is, together with harvest residues, used as energy wood. Therefore, to calculate litter input from harvest residues, a chain of calculations is needed. (1) Of all harvested roundwood volume in Finland (Luke statistics, 2022), the proportion that is collected from drained peatland forests is first derived using the respective share of harvest area data provided by the NFI. (2) This proportion of roundwood volume is then converted to total tree biomass, and (3) using data provided by LUKE Statistics (2022) of energy wood consumption, the part of biomass that is used for energy production is deducted from the total harvested biomass. (4) The remaining biomass is finally converted to different tree biomass components using specific BCEFs. Of these components, the foliage, branches, waste wood, stumps and coarse roots are assumed to remain in the harvest sites and thus to compose the harvest residue litter input.

The yearly production of litter from trees that died naturally is derived from NFI. For these trees, litter input consists of all biomass components, including the stem wood, calculated using appropriate BCEFs. Data on litter input from harvested trees and trees that died naturally are currently available for southern and northern Finland, but not for each FTYPE separately.

Finally, litter input lacks one biomass component, i.e., tree roots in the diameter range of 2–10 mm. The reason is that neither biomass nor turnover rate estimates are available for these roots. Omission of the 2–10 mm root litter input does not affect our estimate of heterotrophic decomposition of that litter fraction as it is based on inclusive empirical emission data; see below.

2.6 Calculating soil CO₂ release, R_Het

Annual CO₂ release from decomposing peat, D_P, and decomposing aboveground and belowground litter of living trees and ground vegetation, D_AGL and D_BGL, is calculated using empirical regression models given in Ojanen et al. (2014). These models, presented in Table 2, include an FTYPE-specific constant, FTYPE May–October mean air temperature (Fig. 3) and FTYPE mean BA (Fig. 3) as dynamic predictors. The data behind the regression models include field measurements on 68 drained peatland forest sites (Ojanen et al., 2010, 2013), which cover the distribution of FTYPEs in Finland (see Fig. 1 in Ojanen et al., 2010). The measurements were originally performed on field plots with aboveground litter removed, but the CO₂ release from decomposing aboveground litter was later calculated using the Yasso07 decomposition model and added to the recorded values of CO₂ release (Ojanen et al., 2013). To focus on CO₂ release from SOM, not confounded by autotrophic root respiration, the field plots were trenched (Ojanen et al., 2013). The predicted CO₂ release is a yearly estimate, although the temperature predictor is the mean May–October air temperature (Ojanen et al., 2014).

The soil CO₂ release estimates produced by the regression models in Table 2 do not include CO₂ release from harvest residues (the empirical data for the models were collected at sites without recent harvests) or from stump and stem wood of trees that died naturally (because these parts of litter did not fit inside the CO₂ flow measurement chambers). The CO₂ release from these litter components is therefore estimated using the Yasso07 decomposition model, combined with time series of southern and northern Finland mean annual air temperature, temperature amplitude and mean annual precipitation (Fig. 4, Supplement Fig. S2). Since the Yasso07 model produces the remaining OM pool after decomposition, the CO₂ release from decomposing aboveground and belowground harvest residues as well as stump and stem wood of trees that recently died naturally, D_AGR and D_BGR, is calculated by subtracting the remaining OM pool from the inputs of these residues, I_AGR and I_BGR.

The Yasso07 model was first initialized for 50 years using constant residue input from harvests and trees that died naturally (mean of 1970–1976 input for southern Finland and mean of 1975–1977 input for northern Finland). The initialized model was then run from 1971 onwards using real annual residue input. During the 50-year initialization, the OM pool originating from harvest residues and trees that died naturally monotonically approached a plateau (Supplement Fig. S2). The rationale for using unusually short spin-up is that there is no long history for drained peatland forests in Finland (drainage activity rose after the 1940s and was most intensive in the 1970s) and the soil OM pool originating from anthropogenic activity has not been developing for long. Running a longer spin-up would create an unnaturally high harvest and natural tree mortality-based soil OM pool, leading to erroneously high CO₂ release from these slowly decomposing C stocks during the GHG inventory reporting period. However, to be aware of the effects of the short spin-up, we also initialized the Yasso07 model for 1000 years. Increasing the length of initialization from 50 to 1000 years increased CO₂ release over the years 1990–2021 annually by 2.4 %–0.6 % (mean 1.4 %) and 4.1 %–1.5 % (mean 2.8 %) in southern and northern Finland, respectively.

2.7 Uncertainty analysis

Uncertainty was assessed for each annual estimate of soil CO₂ balance, as well as for the estimate of change in balance between the years 1990 and 2021. We use the IPCC (2006) guidelines, wherein uncertainty is defined as 1.96 × SEM (standard error of the mean) and given as a percentage of the sink–source estimate. When the uncertainty of the estimate is less than 100 %, zero is not included within the 95 % confidence limits of the estimate. The accounted sources of uncertainty included NFI sampling errors in area and BA estimates, estimation error in the parameters of the models, litter production rates and dwarf shrub coverage, and uncertainty about litter production originating from living trees, harvests and natural mortality. The details of variance propagation are given in Appendix A. These and all other computations were carried out in the R environment (R Core Team, 2020).

2.8 Earlier version of the calculation method

The earlier version of the calculation method of drained organic forest soil CO₂ emissions applied in the GHG inventory of Finland (Statistics Finland, 2022) differs in many ways from the new revised method we present here. In the earlier version, the estimates of heterotrophic soil respiration in FTYPEs are based on results from two peatlands, located at 61^∘ N, 25^∘ E (covering FTYPEs from the herb-rich to dwarf shrub FTYPE) and 62^∘ N, 31^∘ E (representing the Cladonia FTYPE) (Nykänen et al., 1998; Minkkinen et al., 2007). The estimates include D_P and D_BGL only, are static in time, and do not make a difference between southern and northern Finland. Belowground litter input from living trees, ground vegetation, harvest residues and natural mortality is included, whereas aboveground litter input (and most of its decomposition) from living plants and harvests is excluded. Also, in contrast to the new revised method, the belowground litter input from ground vegetation is constant at 1.1 t CO₂-C (or 4.0 t CO₂) ha⁻¹ yr⁻¹ (Laiho et al., 2003), the arboreal fine root turnover rate is 0.85 for all FTYPEs (Liski et al., 2006) and tree leaf biomass is used as a proxy to estimate arboreal fine root biomass (Helmisaari et al., 2007).

Figure 5 Drained organic forest soil CO₂ balance for the whole country (grey area) as well as southern Finland (yellow line) and northern Finland (red line) as produced by the new fully dynamic and the previous semi-dynamic method. See Fig. 1 for the two regions of the country.

3 Results

3.1 Soil CO₂ balance

The method developed in this study produces a generally increasing trend of annual emissions from 0.2 to 2.1 t CO₂ ha⁻¹ yr⁻¹ for the period 1990–2021, equal to a trend from 1.4 to 7.9 Mt CO₂ yr⁻¹ for the entire 4.3 Mha area of drained peatland forests in Finland (Fig. 5). The difference between the years 1990 and 2021 is statistically significant (uncertainty 46 %; Table A10 in Appendix A). Across the period 1990–2021, emissions are on average 1.5 t CO₂ ha⁻¹ yr⁻¹, or 3.4 Mt CO₂ yr⁻¹, for the whole area of 2.2 Mha in southern Finland and −0.14 t CO₂ ha⁻¹ yr⁻¹, or −0.3 Mt CO₂ yr⁻¹, for the whole area of 2.1 Mha in northern Finland (Fig. 5). Considering the uncertainty, drained organic forest soils in southern Finland and in the whole country are a source of CO₂ in 2021, while the 2021 net emission in northern Finland does not significantly differ from zero (Table A9 in Appendix A).

In contrast to the results produced by the new method, the earlier method produces a decreasing trend of total emissions from 12.1 to 3.2 Mt CO₂ for the period 1990–2021 for the whole country, and the emissions are on average lower in southern than northern Finland (Fig. 5).

3.2 Total litter input and decomposition

The new method produces an increasing trend of litter input and decomposition for both southern and northern Finland (Fig. 6). In 2021, total litter input and decomposition are, per unit area, 28 % and 34 % higher in southern than northern Finland, respectively (Fig. 6). In both regions, the increasing trend in total litter input is mainly due to increasing input of aboveground tree litter and belowground arboreal litter (Fig. 6). These litter fractions also include residues from forest harvests and trees that died naturally.

Figure 6 Input of litter fractions, CO₂ release from decomposing litter and peat, and the CO₂ balance (dashed white line) for southern and northern Finland using the new fully dynamic and the previous semi-dynamic method. In the new method, ground vegetation litter consists of both aboveground and belowground litter except for dwarf shrub fine root litter, which is included in belowground arboreal litter.

Figure 7 Input and decomposition of litter from tree harvests and trees that died naturally, as well as their CO₂ balance (dashed line) for southern and northern Finland as produced by the new method. Only stem and stump wood is included for trees that died naturally.

When compared to the earlier method, the new method produces 37 % and 61 % smaller belowground arboreal litter input, per unit area, for southern and northern Finland, respectively (Fig. 6). However, total litter input is on average 62 % and 32 % higher, respectively, because aboveground tree litter and ground vegetation litter are included in the new method (Fig. 6). In 2021, the new method produces 63 % and 10 % higher CO₂ release from decomposition, per unit area, than the previous method for southern and northern Finland, respectively (Fig. 6).

3.3 Input and decomposition of litter from harvested trees and trees that died naturally

With the new method the input of litter from harvested trees and trees that died naturally, as well as the CO₂ release from their decomposition, is significantly higher, per unit area, in southern than northern Finland but increases with time in both regions (Fig. 7). The CO₂ balance of these residues is negative (thus acting as a sink) and has a decreasing overall trend in both regions (Fig. 7). In 2021, the residues form 20 % and 15 % of all litter input, and their decomposition forms 15 % and 12 % of total decomposition in southern and northern Finland, respectively (compare Figs. 6 and 7).

3.4 Comparison of FTYPEs

When litter input and decomposition (excluding litter from harvested trees and trees that died naturally) are estimated with the new method for each FTYPE individually, both litter input and decomposition have a general decreasing trend along the FTYPE fertility gradient from the most fertile herb-rich FTYPE (Rhtkg) to the least fertile Cladonia FTYPE (Jatkg) (Fig. 8). The only exception to this pattern is litter input in the Vaccinium vitis-idaea FTYPE (Ptkg), which in 2021 exceeds litter input in the more fertile Vaccinium myrtillus FTYPE (Mtkg) by 26 % and 20 % in southern and northern Finland, respectively (Fig. 8).

No clear trend of CO₂ balance emerges along the FTYPE fertility gradient (Fig. 8). Instead, net emissions are higher at both ends of the gradient than in the Vaccinium vitis-idaea (Ptkg) and the dwarf shrub (Vatkg) FTYPEs (Fig. 8). In these two FTYPEs, ground vegetation litter input and arboreal fine root litter input are relatively high in comparison to other FTYPEs (Fig. 8). In general, the share of ground vegetation litter in total litter input is higher in northern than southern Finland and increases towards the nutrient-poor end of the FTYPE fertility gradient (Fig. 8).

Figure 8 Litter input and decomposition (both exclude input and decomposition of litter from harvested trees and trees that died naturally), as well as the soil CO₂ balance (dashed line) in FTYPEs in southern and northern Finland as produced by the new method. Ground vegetation litter includes aboveground and belowground litter except for dwarf shrub fine root litter, which is included in arboreal fine root litter. Rhtkg: herb-rich FTYPE, Mtkg: Vaccinium myrtillus FTYPE, Ptkg: Vaccinium vitis-idaea FTYPE, Vatkg: dwarf shrub FTYPE, Jatkg: Cladonia FTYPE, listed in order of decreasing soil fertility.

3.5 Sensitivity of soil CO₂ balance to the main dynamic drivers

As temperature, T (and other climatic variables in Yasso07 modeling), and BA are the main dynamic drivers of soil CO₂ balance in our method, we examined their relative importance by producing three sets of scenarios, where either (1) climate, (2) BA (and the associated harvest rate), (3) both or (4) neither (the default method) were fixed to their 1990 values (Fig. 9): the first set is for the CO₂ balance of tree and ground vegetation litter input and SOM decomposition (top row of graphs, excluding harvest residues and litter from trees that died naturally); the second set is for the CO₂ balance of input and decomposition of litter from harvested trees and trees that died naturally (middle row), and the third set is for the total soil CO₂ balance (bottom row).

The CO₂ balance of tree and ground vegetation litter input and SOM decomposition (Fig. 9 top row) appears to be sensitive to the increasing T as the scenario with constant BA (i.e., only T increases) closely follows the default scenario. The effects of changes in BA are more subtle: increasing BAs until ca. 2010 (Fig. 3) decrease emissions slightly as the scenario with constant climate (only BA changes) produces less CO₂ than the scenario wherein both BA and climate remain constant. From ca. 2010 onwards BAs in FTYPEs level off or start decreasing (Fig. 3) and the BA effect turns neutral (northern Finland) or slightly positive (southern Finland) (Fig. 9 top row).

Figure 9 Soil CO₂ balances calculated with four different scenarios for southern and northern Finland and the whole country: “Default” allows all drivers to change according to data, “BA and harvest constant” fixes the values of BA and harvest rate to that in 1990, “Climate constant” fixes the growing season mean temperature used by equations of SOM decomposition and the climatic variables used by Yasso07 to values recorded in 1990, and “BA + harvest + climate constant” fixes all these variables to their 1990 levels. The top row of panels (a) shows the CO₂ balance of live tree and ground vegetation litter input as well as SOM decomposition (thus excluding harvest residues and litter from trees that died naturally); the middle row (b) shows the CO₂ balance of input and decomposition of litter from harvested trees and trees that died naturally, and the bottom row (c) shows the total soil CO₂ balance.

The CO₂ balance of harvest residues and litter from trees that died naturally (Fig. 9 middle row) is driven by changes in BA and harvest rate as the scenario with constant climate (only BA and harvest rate change) closely follows the default scenario. Changes in climate increase emissions slightly as the scenario with constant BA and harvest rate (only climate changes) produces more CO₂ than the scenario wherein BA, harvest rate and climate all remain constant.

The scenarios for the total soil CO₂ balance (Fig. 9 bottom row) combine the main characteristics of the scenarios for the two parts of the total balance: i.e., the shape and the positive slope of CO₂ emissions in the default scenario are driven by climate warming, and the main result of increasing BA and harvest rate from 1990 onwards is a lower level of net soil CO₂ emissions.

3.6 Ecosystem CO₂ balance

The ecosystem CO₂ balance in drained peatland forests is governed by the combination of soil CO₂ balance and the CO₂ sink created by living trees. When the soil CO₂ balance results of the new method are combined with tree CO₂ sink for southern Finland, the ecosystem CO₂ balance is negative; i.e., the forests are a CO₂ sink until the beginning of the 2010s when the balance turns positive (Fig. 10). This change in pattern is a product of both lower tree CO₂ sink in the 2010s and increasing soil CO₂ emissions (Fig. 10). In northern Finland, drained peatland forests are still a CO₂ sink, but the gradually decreasing tree CO₂ sink together with increasing soil CO₂ emissions also creates a decreasing trend for the ecosystem sink in the 2010s (Fig. 10).

Figure 10 Changes in tree CO₂ sink (grey bars), soil CO₂ balance (red line) and whole-ecosystem CO₂ balance (yellow line) (Mt CO₂ yr⁻¹) in drained peatland forests in southern and northern Finland across the GHG inventory reporting period (see Fig. 2 for the two regions). Note that no GHGs other than CO₂ are included in the calculated balances.

4 Discussion

4.1 The CO₂ balance of drained organic soils by the new method

The CO₂ balance of drained organic soils in Finland, as estimated by the new method we describe in this paper, has two main features (Fig. 5): first, drained peatland soils are currently a source of CO₂ in Finland, and second, the emissions are increasing with time due to the high sensitivity of SOM decomposition to increasing temperature. The mean increase of ca. 0.7 ^∘C in mean May–October temperatures across the time series in FTYPEs in both southern and northern Finland (Fig. 4) increases CO₂ emissions from soil by 8.1 Mt CO₂ in the whole country (the constant BA and harvest scenario in Fig. 9 bottom row). An increase in BA and harvest rate over the decades counteracts this pattern (the constant climate scenario in Fig. 9 bottom row), however, and the increase in net soil CO₂ emissions across the time series remains smaller (the default scenario in Fig. 9 bottom row).

Both litter input and decomposition have increasing trends in time (Fig. 6). The increasing trend in litter input is mainly driven by the general increase in BA across the FTYPEs (Fig. 4), which leads to greater litter input from living trees (evident in both aboveground tree litter input and belowground arboreal litter input in Fig. 6). Part of the increasing trend is also explained by the increasing input of litter from harvested trees and trees that died naturally (Fig. 7). Litter input from ground vegetation, instead, remains rather stable across the time series (Fig. 6). The increasing CO₂ release from decomposition is driven by both increasing BA and increasing May–October temperature (Fig. 3) as both are positively associated with decomposition in the empirical data (Table 2; Ojanen et al., 2014). Part of the increasing trend is also explained by the increasing CO₂ release from the decomposing litter of harvested trees and trees that died naturally (Fig. 7). The likely reasons why BA is positively linked to decomposition in the regression models by Ojanen et al. (2014) (Table 2) is that increasing litter input sustains larger decomposer biomass, and thus higher CO₂ production, and that higher BA – through higher evapotranspiration – maintains a lower soil water table level and a deeper oxic layer for SOM decomposition (Sarkkola et al., 2010). When the impacts on litter input and decomposition are combined, increasing BA first slightly decreases net soil CO₂ emissions (scenario of constant climate in Fig. 8 top row), but this effect disappears in the mid-2010s, when the positive trend in BA levels off or turns negative in most FTYPEs (Fig. 3).

Release of CO₂ from decomposition per area of drained peatland forest soil is significantly higher in southern than northern Finland (Fig. 6). This is due to higher temperature, but also due the generally higher BA and litter production (Figs. 3, 4, 6) in southern Finland. Input of litter from forest harvests has an increasing trend, and on the timescale of the inventory time series, harvest residues create a CO₂ sink and decrease net soil CO₂ emissions (Fig. 7). However, this does not imply that net ecosystem CO₂ emissions would decrease with harvests because the tree CO₂ sink simultaneously decreases (see trends in the 2010s in Fig. 10).

Of the FTYPEs, the dwarf shrub (Vatkg) and V. vitis-idaea (Ptkg) FTYPEs are clustered to the major water divides called Suomenselkä and Maanselkä of the southern and northern regions of the country, while other FTYPEs are more evenly distributed across the country (Fig. 2). These two common FTYPEs have relatively low nutrient status (Table 1) and both have, even without the harvest residue impact, in most of the GHG inventory time series a negative or near-zero soil CO₂ balance in the north, and the V. vitis-idaea FTYPE also in the south (Fig. 8). This is roughly in line with earlier findings (Minkkinen et al., 1999; Ojanen et al., 2014). However, the ubiquitous upwards temporal trend of soil CO₂ balance in all FTYPEs (Fig. 8) suggests that these FTYPEs are also entering the phase of net soil CO₂ emissions.

Uncertainty due to model parameters and input data were propagated into the soil CO₂ balance estimates. The parameter estimates of the peat and litter decomposition model were clearly the greatest source of uncertainty, followed by those in the fine root litter input model (including the turnover rates) and litter input from living trees. Meanwhile, uncertainty in NFI estimates of site type areas and BA had a negligible impact on annual balance estimates and a somewhat greater, but still minor, impact on the estimates of changes in emissions. This indicates that the proposed method might also be applicable in somewhat smaller regions (with less precise NFI estimates) without a drastic increase in uncertainty. Errors in annual soil CO₂ balance estimates are strongly correlated over the whole time series because the same model parameters with the same estimation errors are applied throughout the series. Since all applied models were linear in both their parameters and inputs, it was possible to handle these correlations by combining quantities including the same parameter (Appendix A). Uncertainty propagation by means of variances relies on an assumption of symmetric error distributions. This assumption is realistic since uncertainty due to all individual components in the balance estimates was 30 % or less; i.e., the coefficient of variation was at most 0.15 (Table A9; see Frey et al., 2006, Sect. 3.2.2.4). Structural uncertainty about the models could not be accounted for since large-scale validation data do not exist.

4.2 Ecosystem CO₂ balance in drained peatland forests

Peatlands in Finland were mostly drained to enhance tree growth, thus enabling forestry on peatlands. The increasing BAs in drained peatlands from 1990 onwards (Fig. 4) show that this target has largely been achieved. Therefore, when considering the net climate effects of peatland drainage, the development of the tree CO₂ sink needs to be taken into consideration (Fig. 9). In the 2010s, the trend of gradually increasing BA in drained peatland forests leveled off (Fig. 3). This is due to increasing harvests in these forests (manifested by the increasing harvest residue input in Fig. 7), which in turn is associated with the sharp turn in tree CO₂ sink and net ecosystem balance after 2010, ultimately leading to a shift from a net ecosystem sink to a source of CO₂ in southern drained peatland forests (Fig. 9). After 2010, drained peatland forests were a sink of CO₂ only in northern Finland, but when the northern tree stands mature for harvesting and industry demand for wood remains high, the net ecosystem CO₂ sink will likely be lost in northern Finland as well. Calculated for the whole country, forests growing on drained peatlands were a net sink of 0.2 Mt CO₂ in 2021.

4.3 Comparing the soil CO₂ balance predictions by the new and earlier method

The differences in the estimates of soil CO₂ balance between the new fully dynamic and the earlier semi-dynamic method are remarkable (Fig. 4). The main difference is the opposite trends in time series: at the beginning of the GHG inventory reporting period the new method produces 10.7 Mt lower soil CO₂ balance than the earlier method, but then it reaches 3.4 Mt higher balance in recent years. The mean level of emissions differs a lot in the north and less in the south (Fig. 5).

The contrasting trends of time series produced by the two methods call for explanations. In the earlier method, CO₂ release from decomposition is constant, whereas in the new method CO₂ release from decomposition increases in time as it is driven by increasing temperatures and BA (Fig. 6). Both methods predict increasing litter input in time due to increasing BA (Fig. 6). In the earlier method, however, increasing litter input combined with constant decomposition unavoidably leads to decreasing net emissions as the release of CO₂ in decomposition does not follow the increase in litter input. In the new method, both constituents of the CO₂ balance are dynamic, and the method can produce either an increasing or decreasing trend of net emissions depending on the relative magnitude of slopes in trends for litter input and decomposition. The new method predicts that decomposition increases more in time than litter input does, thus leading to an increasing trend in net soil CO₂ emissions.

Besides the difference in the estimation of decomposition, estimation of litter input has major differences between the methods. First, while the aboveground litter (including litter from living trees and ground vegetation, as well as from harvested trees and trees that died naturally) is explicitly included in both litter input and decomposition calculations in the new method, the earlier method assumes that aboveground litter input and decomposition are in equilibrium, and consequently, aboveground litter is not included in the calculation. However, findings from drained organic forest soils suggest that while increasing litter input is generally associated with increasing decomposition, litter input and decomposition are seldom in equilibrium. For instance, in Finland, litter production exceeds decomposition in the nutrient-poor FTYPEs and vice versa in the nutrient-rich FTYPEs (Fig. 1 in Ojanen et al., 2013). Other studies also stress the importance of litter input in the C balance of drained peatland forests, and despite plant-community-level litter decomposition rates being somewhat higher in drained than undrained peatland forests, increased tree litter input can still lead to a significant accumulation of new SOM, especially in conifer-dominated stands (Vávřová et al., 2009; Straková et al., 2012). A model-based assessment of forests growing on mineral soils in Finland further suggests that the proportion of C in litter input that remains in the soil C pool is ca. 4 % (Liski et al., 2006), and although these models are developed for mineral soil forests, the main principles of SOM decomposition also likely apply in drained organic soil forests. Finally, as the maturing drained peatland forests are increasingly harvested, litter production from living trees will inevitably decrease and that from harvest residues increase, making it necessary to explicitly include the input and decomposition of aboveground litter in the calculation method.

Another major change in the estimation of litter input in the new method is the FTYPE-specific turnover rates of arboreal fine roots, which constitute the major source of belowground litter. In light of recent findings, the assumption of a single turnover rate of 0.85, used for all FTYPEs in the earlier inventory method (Statistics Finland, 2022), seems incorrect. The new minirhizotron results included in the new method suggest that the tree fine root turnover rate varies markedly by FTYPE and that the turnover rate is on average lower in the nutrient-poor than in the nutrient-rich FTYPEs (Table 3). The lower turnover rate largely explains the remarkably lower arboreal fine root litter production in the new method (Fig. 6).

4.4 Evaluating the predicted soil CO₂ balance

To evaluate the soil CO₂ balance predicted by our method, it can be compared to (1) IPCC Tier 1 emission factors, (2) emission factors applied in other European countries, and (3) empirical results on soil CO₂ and C balance. As there are no alternative 30-year dynamic emission time series available, the validity of time series dynamics cannot be adequately evaluated here – this would require another comprehensive study. Instead, we will focus on the question of whether the predicted soil CO₂ balance is at a feasible level in comparison to emission factors utilized in other GHG inventories and in comparison to results presented in earlier studies.

The default IPCC emission factors for boreal drained organic soils range from 0.25 (95 % confidence interval −0.23–0.73) t CO₂-C ha⁻¹ yr⁻¹ in drained nutrient-poor forest land to 0.93 (0.54–1.3) t CO₂-C ha⁻¹ yr⁻¹ in nutrient-rich forest land (IPCC, 2014, their Table 2.1 in Chap. 2). Applying the IPCC default emission factors for the areas of nutrient-poor and nutrient-rich FTYPEs found in southern and northern Finland (Table 1) gives emission factors of 0.56 and 0.49 t CO₂-C ha⁻¹ yr⁻¹ for these regions, respectively, and an overall factor of 0.53 t CO₂-C ha⁻¹ yr⁻¹ for the whole country. The emission rate of 0.50 t CO₂-C ha⁻¹ yr⁻¹ predicted by our method for Finland for the year 2021 is slightly lower, but as the emissions with the new method have a rising trend (Fig. 5), they will likely soon exceed the default emissions.

To compare our results to emission factors used for drained organic forest soils by other European countries, we collected information from the national CRF (Common Reporting Format by the United Nations Framework Convention on Climate Change, UNFCCC) tables published in 2021. Russia, Norway, Ireland and Poland use emission factors of 0.52–0.72 t CO₂-C ha⁻¹ yr⁻¹ that are at the lower end of the 95 % CI of IPCC default method emissions for boreal nutrient-rich forest land. In temperate areas the IPCC default emission is 2.6 (2.0–3.3) t CO₂-C ha⁻¹ yr⁻¹, and temperate countries Germany, Denmark, Switzerland and the UK apply emission factors of 2–2.6 t CO₂-C ha⁻¹ yr⁻¹. As expected, the emissions predicted by our method for Finland in the northern boreal forest zone are lower than the emission factors used in countries that are partly or completely located in the temperate zone.

In a compilation of all available empirical soil CO₂ balance data in Finland, Ojanen and Minkkinen (2019) estimated emission factors for nutrient-rich (0.66 t CO₂-C ha⁻¹ yr⁻¹) and nutrient-poor (−0.19 t CO₂-C ha⁻¹ yr⁻¹) forestry-drained peat soils. Weighting these by the FTYPE areas of our study (Table 1) gives a mean emission factor of 0.37 t CO₂-C ha⁻¹ yr⁻¹. A previous country-level upscaling study by Ojanen et al. (2014) utilizing the heterotrophic soil respiration–litter production method like the one used in our study (Eq. 1) gives a maximum value (95 % CI) of 0.55 t CO₂-C ha⁻¹ yr⁻¹ for the mean emission factor. Other soil CO₂ balance (Bjarnadottir et al., 2021; Korkiakoski et al., 2023; Minkkinen et al., 2018; Uri et al., 2017; Meyer et al., 2013) and C balance (Minkkinen et al., 1999; Minkkinen and Laine, 1999; Simola et al., 2012; Lupikis and Lazdins, 2017) studies from boreal and hemiboreal northern Europe further reinforce the view that forestry-drained peat soils have relatively low CO₂ emissions; and while drainage has often led to soil C loss, nutrient-poor soils in particular may exhibit soil C sink even after drainage.

Based on the above paragraphs, we can state that the soil CO₂ balance predicted by our new method is in line with previous knowledge and within the range of previous results. Even the relatively large change from close-to-zero emissions in the 1990s to close to 10 Mt CO₂ yr⁻¹ in 2021 is within the limits of previous knowledge. According to our uncertainty assessment (Appendix A), which covers sampling error and parameter uncertainty of the original models (Ojanen et al., 2014) and the present GHG inventory time series application, this change in emissions is statistically significant. As it is virtually impossible to produce empirical country-level 30-year soil CO₂ emission time series for comparison, a comparison with mechanistic modeling results would be a viable option for future research. Mechanistic models that could couple hydrology and the C cycle and are capable of large-scale modeling of forestry-drained peat soils are currently not available, but recent developments in hydrological modeling (e.g., Hökkä et al., 2021; Laurén et al., 2021) appear promising for the development of such models.

4.5 Critical assumptions in the new method

It is well established that water table depth (WTD) is the master variable controlling decomposition in drained wetlands (e.g., Silvola et al., 1996; Ojanen et al., 2014; Jauhiainen et al., 2019), and one can argue that WTD should be included in methods that aim at predicting CO₂ release from drained peatland soils. In fact, Ojanen et al. (2010), using the data behind our new method, found that regression models between soil respiration and climate variables had a higher coefficient of determination, R², when WTD was included (see Ojanen et al., 2010, Table 4). However, because direct data on WTD, in particular for the time series starting in 1990, are not available nationwide for GHG inventory purposes in Finland, WTD cannot currently be used as a predictor of decomposition. In undrained peatlands, soil CO₂ balance fluctuates following annual variation in water saturation and temperature (e.g., Alm et al., 1997, 1999), tree growth is suppressed by anoxia, and BA remains low. Draining initiates tree growth and reduces critical high water level periods, and finally, the increasing tree growth and BA lead to lowering water table level (Sarkkola et al., 2010). The impact of WTD on decomposition is therefore implicitly included in our method through the impact of BA, which provides a proxy for the rate of evapotranspiration that largely controls WTD in forestry-drained peatlands (Hökkä et al., 2021; Leppä et al., 2020).

The in situ minirhizotron data behind the arboreal fine root turnover rates that we apply are rare, and the results we use are the first to compare arboreal fine root growth and mortality in different boreal FTYPEs. The measurements carried out with minirhizotrons focused on tree fine roots with diameter ≤ 0.5 mm, but we extrapolated the results for all arboreal roots ≤ 2 mm in diameter. Because thicker fine roots have lower turnover rates, our method may overestimate litter input of roots with a diameter of 0.5–2 mm. In terms of soil CO₂ balance, this would lead to underestimation of net CO₂ emissions. On the other hand, as the fine root turnover rates of dwarf shrubs are higher than those of trees, at least when compared to Scots pine (Kari Minkkinen et al., unpublished data), our method may in this respect underestimate the arboreal root litter input and overestimate the net CO₂ emissions.

What may impair the use of BA when estimating SOM decomposition is that relatively few data are available on CO₂ release from peatland forest soils after clear-cutting and regeneration or after other forest management options such as cuttings applied in continuous-cover silviculture. Mäkiranta et al. (2010) showed that CO₂ release from soil decreased after forest clear-cutting because of a rising water table and, on the other hand, because of the dryness of the peat soil surface. Instead, CO₂ release from decomposing harvest residues was high, and because of the removed tree stand and suffering ground vegetation, the site acted as a high source of CO₂ emissions during the whole study of 4 years. Similarly, eddy covariance observations in a drained peatland forest showed markedly elevated CO₂ emissions at least 2 years after clear-cutting (Korkeakoski et al., 2019), thus suggesting significant transitory effects of clear-cutting on soil CO₂ balance in peatland forests. However, similarly to Mäkiranta et al. (2010) study, a marked share of increased emissions likely also resulted from decomposition of fresh harvest residues in the Korkeakoski et al. (2019) study. While this harvest residue effect is included in our method in both litter input and decomposition (Fig. 7), the temporary effects of tree cut on ground vegetation are not explicitly implemented.

One fundamental assumption in our new method is that the statistical relationships between BA, ground vegetation litter input and SOM decomposition described by Ojanen et al. (2014) for different FTYPEs can also adequately predict these associations in forests under their future successional change. If future forest stands have features that were not present in the studied field sites, this assumption is not valid. This also has relevance in terms of the uncertainty analysis. A strong structural assumption in our new method is that model parameters, including the litter production and turnover rates, do not change in time. If trend-like changes do exist, the reported uncertainties of change estimates are too small (see Lehtonen and Heikkinen, 2016, for an illustration). Better understanding of temporal dynamics would be needed to make the assumptions and uncertainty assessments more comprehensive.

Our method also assumes that the spatial covariation in temperature and decomposition, found by Ojanen et al. (2014) among their field sites, is a causal relationship that can be used to describe the effects of increasing temperature on SOM decomposition under climate warming. In our time series, the May–October mean temperature rises by 0.67 ^∘C and the mean annual temperature by 1.2 ^∘C during the GHG reporting period (Fig. 3). This upward trend in temperature increases decomposition in our method (Table 2) and will reduce the C sink of forest ecosystems in drained peatlands if it is not simultaneously accompanied by increasing nutrient mineralization, plant growth and litter input. It is well established that heterotrophic soil respiration is sensitive to temperature (Silvola et al., 1996; Kätterer et al., 1998; Meyer et al., 2018), and while the response to temperature may be somewhat attenuated with time due to a gradual decrease in easily degradable C and thermal adaptation of decomposer communities and their respiration (Davidson and Janssens, 2006; Bradford et al., 2008; Bradford, 2013), it is likely that the trend of increasing SOM decomposition produced by our method represents a real outcome of climate warming. Globally the same impact of increasing temperature on soil CO₂ emissions in drained peatlands is demonstrated by the steeply increasing emission factors from high latitudes to the tropics (IPCC, 2014, their Table 2.1 in Chap. 2).

4.6 Further development needs

The impact of WTD on SOM decomposition is implicitly included in our method through BA and T, which are both proxies for the rate of evapotranspiration. However, the relationship between BA and WTD is not linear but attenuating at high BAs (Sarkkola et al., 2010). Furthermore, while ditch spacing and ditch depth have a minor role in comparison to BA, they also affect WTD (Sarkkola et al., 2010; Hökkä et al., 2021), and including models that could link ditch parameters to soil CO₂ release would likely improve our method. Models to overcome the shortcomings of large-scale WTD estimation in forest drainage are being tested (e.g., Haahti et al., 2015), but finally resolving these issues may require transforming the method into a process-based model, which would then also open avenues for addressing a wider range of environmental effects on decomposition.

Also, while our new method covers a wide range of litter input, litter input from tree roots thinner than coarse roots but thicker than fine roots, i.e., in the diameter range of 2–10 mm, is lacking. This is because there are no data on the standing biomass or litter production rate for this fraction of tree roots. In a Swedish study, Petersson and Ståhl (2006) showed that the belowground tree biomass was on average 11 % higher when roots down to 2 mm diameter were included in comparison to using Marklund's (1987, 1988) biomass equations, which produce the biomass of stump and thicker coarse roots only. Using this estimate, we could assume that we miss ca. 10 % of the tree belowground biomass (i.e., of the combined stump and coarse root biomass) and correct the tree belowground biomass estimate accordingly. However, the turnover rate of this biomass fraction is unknown.

We use 30-year rolling means of climatic variables in our models, computed to the last year of the sequence. This procedure is aimed at smoothing the effects of annual weather variation on the results and thus better reveal the mean level of soil CO₂ balance. However, when the smoothed variable has a monotonous trend in time, which is the case, for instance, for May–October or annual mean temperature that increases in time (Fig. 4), the 30-year rolling mean temperatures show past developments behind the present temperature. In our case, this means that the May–October mean air temperature, which largely determines the peat and litter decomposition (Fig. 9 top row), is constantly too low, which in turn leads to CO₂ emission estimates that are too low. The origin of applying a rolling mean is in the similar solution used in mineral soils. For a more agile smoothed value, another type of filter would be needed.

Finally, although not adopted in our method, peatlands also release organic C as dissolved and particulate organic carbon, DOC and POC into surrounding waters in their natural stage, and drainage increases the release (Nieminen et al., 2021). The IPCC 2013 Wetlands Supplement (IPCC, 2014) provides default rates of waterborne C release in the form of increased DOC leaching over the natural baseline (IPCC, 2014, their Chap. 2, Eq. 2.5). A simple application of the IPCC Tier 1 value of 0.12 (0.07–0.19) t C ha⁻¹ (IPCC, 2014, their Chap. 2, Table 2.2), multiplied by the total drained forest land area of 4.3 Mha on organic soil in Finland, gives a rough estimate of peat C loss as DOC of 0.52 Mt C. Using the IPCC default assumption that 90 % of DOC is finally released as CO₂ to the atmosphere gives an estimate of 1.70 Mt CO₂ above the soil CO₂ emissions predicted by our new method. This corresponds to a shift of about 21 % of the net soil CO₂ balance in 2021 (Fig. 5), the amount worthy of considering as part of the GHG inventory of drained peatlands. Although the fate of DOC in water courses is not entirely clear (a significant part of DOC likely ends up fueling CO₂ emissions, but the share of humus that re-sediments to rivers, lakes and seas is difficult to quantify), from the point of view of peatland forest soils, the runaway organic C reduces the peat C pool and displays one more aspect of the anthropogenic impacts of drainage worth acknowledging.

5 Conclusions

Clear merits of our method are that it includes comprehensive empirical data and models for SOM decomposition and litter production, and it links these to successional changes in tree stand characteristics, monitored by NFI, and to a regionally precise temperature regime. The method also includes realistic propagation of uncertainties in all model inputs and parameters into the annual estimates of CO₂ balance and their differences, which is a requirement for GHG inventories. Such a dynamic model that can consider the effects of climate change as well as forest management and development provides a greatly improved tool for forest policy guidance to mitigate climate change.

Appendix A: Assessment of uncertainty

Uncertainty assessments were developed for all annual estimates of soil CO₂ balance in southern Finland, northern Finland and the whole country, as well as for the estimates of change in balance between the years 1990 and 2021. The accounted sources of uncertainty included NFI sampling errors in estimates of the areas of drained peatland forest types along with basal area and biomass of trees, as well as uncertainty about the parameters of the models, litter production and turnover rates, and mean dwarf shrub coverage.

A1 Notation

Soil CO₂ balance (Mt yr⁻¹, net emission positive) in year t over forests representing site type p in region r (southern Finland or northern Finland) was estimated as Y_rpt=A_rpty_rpt, where A_rpt is the interpolated NFI estimate of the area (Mha) of site type p in region r and

$$\begin{array}{}\text{(A1)}& y_{rpt}=R_{rpt}/\mathrm{100}-\mathrm{44}\cdot \mathrm{0.5}\left(G_{rpt}/\mathrm{100}+F_{rpt}/\mathrm{100}+L_{rpt}+H_{rt}\right)/\mathrm{12}\end{array}$$

is the estimated net CO₂ exchange per area unit (Mg ha⁻¹ yr⁻¹). The components of y_rpt are the following:

$$\begin{array}{}\text{(A2)}& R_{rpt}={\mathit{\alpha }}_{\mathrm{1}}{B}_{rpt}+{\mathit{\alpha }}_{\mathrm{2}}{T}_{rpt}+{\mathit{\beta }}_{\mathrm{1},p},\end{array}$$

the prediction of peat and litter decomposition (g CO₂ m⁻² yr⁻¹), obtained using Eq. (13) of Ojanen et al. (2014, their Table A.5);

$$\begin{array}{}\text{(A3)}& G_{rpt}={\mathit{\alpha }}_{\mathrm{3}}{B}_{rpt}+{\mathit{\beta }}_{\mathrm{2},p},\end{array}$$

the biomass of litter from ground vegetation (g m⁻² yr⁻¹), as predicted using Eq. (12) of Ojanen et al. (2014, their Table A.4);

$$\begin{array}{}\text{(A4)}& F_{rpt}={\mathit{\alpha }}_{\mathrm{4}}{\mathit{\varphi }}_p\left(\sum _s{\mathit{\delta }}_s{B}_{rpst}+{\mathit{\alpha }}_{\mathrm{5}}{\mathit{\mu }}_p+{\mathit{\gamma }}_r\right),\end{array}$$

the biomass of litter from fine roots (g m⁻² yr⁻¹), as predicted using Eq. (8) of Ojanen et al. (2014, their Table A.2), corrected according to Laiho and Finér (1996) and combined with the site-type-specific turnover rates; L_rpt, the biomass of litter from living trees (excluding fine roots, Mg ha⁻¹ yr⁻¹); and H_rt, the biomass of residual organic matter from tree harvesting and natural mortality after decomposition (Mg ha⁻¹ yr⁻¹).

In the equations above, $B_{rpt}={\sum }_s{B}_{rpst}$ is the basal area of trees (m² ha⁻¹) based on NFI, s refers to tree species group, T_rpt is the time-smoothed mean May–October air temperature (^∘C) over sites of type p in region r and the remaining static components of the predictors are based on earlier studies. Model parameters

$$\begin{array}{}\text{(A5)}& {\mathit{\theta }}_R={\left[{\mathit{\alpha }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}{\mathit{\beta }}_{\mathrm{1},\mathrm{Rhtkg}}{\mathit{\beta }}_{\mathrm{1},\mathrm{Mtkg}}{\mathit{\beta }}_{\mathrm{1},\mathrm{Ptkg}}{\mathit{\beta }}_{\mathrm{1},\mathrm{Vatkg}}{\mathit{\beta }}_{\mathrm{1},\mathrm{J}\text{ä}\mathrm{tkg}}\right]}^{\prime }\end{array}$$

are based on Ojanen et al. (2014, their Table A.5),

$$\begin{array}{}\text{(A6)}& {\mathit{\theta }}_{\mathrm{G}}={\left[{\mathit{\alpha }}_{\mathrm{3}}{\mathit{\beta }}_{\mathrm{2},\mathrm{Rhtkg}}{\mathit{\beta }}_{\mathrm{2},\mathrm{Mtkg}}{\mathit{\beta }}_{\mathrm{2},\mathrm{Ptkg}}{\mathit{\beta }}_{\mathrm{2},\mathrm{Vatkg}}{\mathit{\beta }}_{\mathrm{2},\mathrm{J}\text{ä}\mathrm{tkg}}\right]}^{\prime }\end{array}$$

on Ojanen et al. (2014, their Table A.4) and

$$\begin{array}{}\text{(A7)}& {\mathit{\theta }}_{\mathrm{F}}={\left[{\mathit{\delta }}_{\mathrm{pine}}{\mathit{\delta }}_{\mathrm{spruce}}{\mathit{\delta }}_{\mathrm{deciduous}}{\mathit{\alpha }}_{\mathrm{5}}{\mathit{\gamma }}_{\mathrm{south}}{\mathit{\gamma }}_{\mathrm{north}}\right]}^{\prime }\end{array}$$

on Ojanen et al. (2014, their Table A.2). Mean dwarf shrub coverages μ_p are from Table A.3 of Ojanen et al. (2014), correction α₄=1.043 is based on Laiho and Finér  (1996) and turnover rates ϕ_p are reported in this article (Table 3). q=5 is the number of peatland forest site types, and the applied values of the site-type-specific parameters are listed in Tables 3–6. For $p\in \mathit{\{}\mathrm{2},\mathrm{3}\mathit{\}}$, the values β_1,p, β_2,p and μ_p are area-weighted averages of the values in Ojanen et al. (2014, their Tables A.5, A.4 and A.3) over the two subtypes I and II (w_MtkgI=0.615, w_PtkgI=0.6).

The region-specific annual estimates were $Y_{rt}={\sum }_p{Y}_{rpt}$, country-wide annual estimates were $Y_t={\sum }_r{Y}_{rt}$, and corresponding change estimates were $Y_{r,\mathrm{2021}}-Y_{r,\mathrm{1990}}$ and Y₂₀₂₁−Y₁₉₉₀.

A2 Uncertainty components

Uncertainty in estimates A_rpt and B_rpst due to NFI sampling was based on NFI12 (Tables A1 and A2) relying on the assumptions that (i) the relative standard errors (RSEs) ${\mathrm{Var}}^{\mathrm{1}/\mathrm{2}}\left(A_{rpt}\right)/A_{rpt}$ and ${\mathrm{Var}}^{\mathrm{1}/\mathrm{2}}\left(B_{rpt}\right)/B_{rpt}$ are unchanged across time, i.e., the same and equal to the corresponding RSEs in Tables A1 and A2 for all t, and (ii) the estimates A_rp,1990 and A_rp,2016 are uncorrelated, as are B_rps,1990 and B_rps,2021. Assumption (i) is justified by the relatively constant sampling effort of NFI across time and assumption (ii) by the fact that the estimates for the 2 years are completely based on different NFI campaigns.

Uncertainty due to estimation of parameters in Eqs. (8), (12) and (13) of Ojanen et al. (2014) was based on covariance matrices Σ_R, Σ_G and Σ_F of parameter vectors θ_R, θ_G and θ_F (Tables A3–A5) derived from Tables A.5, A.4 and A.2 of Ojanen et al. (2014).

In particular, $\mathrm{Var}\left({\mathit{\beta }}_{\mathrm{1},\mathrm{Mtkg}}\right)=w_{\mathrm{MtkgI}}^{\mathrm{2}}\mathrm{Var}\left({\mathit{\beta }}_{\mathrm{MtkgI}}\right)+w_{\mathrm{MtkgII}}^{\mathrm{2}}\mathrm{Var}\left({\mathit{\beta }}_{\mathrm{MtkgII}}\right)+\mathrm{2}{w}_{\mathrm{MtkgI}}{w}_{\mathrm{MtkgII}}\mathrm{Cov}\left({\mathit{\beta }}_{\mathrm{MtkgI}}{\mathit{\beta }}_{\mathrm{MtkgII}}\right)$, and so on.

Variances of site-type-specific estimates of dwarf shrub coverage, Var(μ_p), were similarly derived from those presented in Table A.3 of Ojanen et al. (2014), and variances of fine root turnover rates, Var(ϕ_p), were based on expert judgment (Table A6), as was that of the deep root correction α₄, Var(α₄)=0.012².

Relative standard errors of litter production estimates from living trees on drained peatland sites were estimated by region from NFI8 and NFI11 as explained in Lehtonen and Heikkinen (2016). For annual results, RSEs from NFI11 were applied to the litter series L_rpt across all site types and times (Table A7). When aggregating to the whole country and assessing the uncertainty in change between the years 1990 and 2021, correlations between regions and time points (Table A8) due to applying the same biomass models and litter production rates were also taken into account.

Relative standard error of litter originating from harvests and natural mortality was estimated similarly from NFI11 and was applied to the residual series H_rt.

Apart from correlations discussed above, all sources of uncertainty were assumed mutually uncorrelated.

A3 Impacts of uncertainty components on soil CO₂ balance estimates

Assuming negligible correlation between area estimates A_rpt, variances of region-specific annual balance estimates $Y_{rt}={\sum }_p{A}_{rpt}{y}_{rpt}$ due to NFI sampling variance in area estimates A_rpt were obtained as

$$\begin{array}{}\text{(A8)}& {\mathrm{Var}}_A\left(Y_{rt}\right)={\sum }_p{y}_{rpt}^{\mathrm{2}}\mathrm{Var}\left(A_{rpt}\right),\end{array}$$

variances of country-wide balanced estimates as

$$\begin{array}{}\text{(A9)}& {\mathrm{Var}}_A\left(Y_t\right)={\sum }_{r,p}{y}_{rpt}^{\mathrm{2}}\mathrm{Var}\left(A_{rpt}\right),\end{array}$$

and variances of the change in balance between the years 1990 and 2021 as

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{\mathrm{Var}}_A\left(Y_{r,\mathrm{2021}}-Y_{r,\mathrm{1990}}\right)& =\sum _p[y_{rp,\mathrm{2021}}^{\mathrm{2}}\mathrm{Var}\left(A_{rp,\mathrm{2021}}\right)\\ & +y_{rp,\mathrm{1990}}^{\mathrm{2}}\mathrm{Var}\left(A_{rp,\mathrm{1990}}\right)]\end{array}\end{array}$$

and

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}{\mathrm{Var}}_A\left(Y_{\mathrm{2021}}-Y_{\mathrm{1990}}\right)& ={\sum }_{r,p}[y_{rp,\mathrm{2021}}^{\mathrm{2}}\mathrm{Var}\left(A_{rp,\mathrm{2021}}\right)\\ & +y_{rp,\mathrm{1990}}^{\mathrm{2}}\mathrm{Var}\left(A_{rp,\mathrm{1990}}\right)].\end{array}\end{array}$$

To evaluate the impact of NFI sampling variance in basal area estimates B_rpst, balance estimates were rewritten as

$$\begin{array}{}\text{(A12)}& Y_{rt}=\sum _p{C}_{\mathrm{1},rpt}+\sum _{p,s}{C}_{\mathrm{2},rpst}{B}_{rpst},\end{array}$$

where intercepts

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{C}_{\mathrm{1},rpt}& =A_{rpt}\mathit{\{}{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{2}}{T}_{rpt}+{\mathit{\beta }}_{\mathrm{1},p}}{\mathrm{100}}}\\ & -{\displaystyle \frac{\mathrm{44}\cdot \mathrm{0.5}}{\mathrm{12}}}[{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{2},p}+{\mathit{\alpha }}_{\mathrm{4}}{\mathit{\varphi }}_p\left({\mathit{\alpha }}_{\mathrm{5}}{\mathit{\mu }}_p+{\mathit{\gamma }}_r\right)}{\mathrm{100}}}\\ & +L_{rpt}+H_{rt}\left]\mathit{\right\}}\end{array}\end{array}$$

do not depend on basal area estimates and

$$\begin{array}{}\text{(A14)}& C_{\mathrm{2},rpst}=A_{rpt}\left\{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{\mathrm{100}}}-{\displaystyle \frac{\mathrm{44}\cdot \mathrm{0.5}}{\mathrm{12}}}\left[{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{3}}+{\mathit{\alpha }}_{\mathrm{4}}{\mathit{\varphi }}_p{\mathit{\delta }}_s}{\mathrm{100}}}\right]\right\}.\end{array}$$

Then, for the annual region-specific balance estimates,

$$\begin{array}{}\text{(A15)}& {\mathrm{Var}}_B\left(Y_{rt}\right)=\sum _{p,s}{C}_{\mathrm{2},rpst}^{\mathrm{2}}\mathrm{Var}\left(B_{rpst}\right),\end{array}$$

and the variances for the whole country and for the change are obtained similarly (see area estimates above).

Considering uncertainty in the parameters of the peat and litter decomposition model, the balance estimates can also be expressed as

$$\begin{array}{}\text{(A16)}& Y_{rt}=C_{\mathrm{3},rt}+C_{\mathrm{4},rt}{\mathit{\alpha }}_{\mathrm{1}}+C_{\mathrm{5},rt}{\mathit{\alpha }}_{\mathrm{2}}-\sum _p{C}_{\mathrm{6},rpt}{\mathit{\beta }}_{\mathrm{1},p},\end{array}$$

where C_3,r does not depend on these parameters:

$$\begin{array}{}\text{(A17)}& {\displaystyle }{C}_{\mathrm{4},rt}=\sum _p{\displaystyle \frac{A_{rpt}{B}_{rpt}}{\mathrm{100}}},\text{(A18)}& {\displaystyle }{C}_{\mathrm{5},rt}=\sum _p{\displaystyle \frac{A_{rpt}{T}_{rpt}}{\mathrm{100}}}\end{array}$$

and

$$\begin{array}{}\text{(A19)}& C_{\mathrm{6},rpt}=A_{rpt}/\mathrm{100}.\end{array}$$

Thus, the impact of uncertainty in the parameters of the peat and litter decomposition model is

$$\begin{array}{}\text{(A20)}& {\mathrm{Var}}_R\left(Y_{rt}\right)={\mathit{c}}_{rt}^{\prime }{\mathbf{\Sigma }}_R{\mathit{c}}_{rt},\end{array}$$

where ${\mathit{c}}_{rt}={\left[C_{\mathrm{4},rt}{C}_{\mathrm{5},rt}{C}_{\mathrm{6},r,\mathrm{Rhtkg},t}{C}_{\mathrm{6},r,\mathrm{Jatkg},t}\right]}^{\prime }$. The corresponding variances for the whole country were obtained using coefficients

$$\begin{array}{}\text{(A21)}& {\displaystyle }{C}_{\mathrm{4},t}=\sum _{r,p}{\displaystyle \frac{A_{rpt}{B}_{rpt}}{\mathrm{100}}},\text{(A22)}& {\displaystyle }{C}_{\mathrm{5},t}=\sum _{r,p}{\displaystyle \frac{A_{rpt}{T}_{rpt}}{\mathrm{100}}}\end{array}$$

and

$$\begin{array}{}\text{(A23)}& C_{\mathrm{6},pt}=\sum _r{\displaystyle \frac{A_{rpt}}{\mathrm{100}}}\end{array}$$

and for the region-specific change estimates with coefficients

$$\begin{array}{}\text{(A24)}& C_{\mathrm{4},r}=\sum _p{\displaystyle \frac{A_{rp,\mathrm{2016}}{B}_{rp,\mathrm{2016}}-A_{rp,\mathrm{1990}}{B}_{rp,\mathrm{1990}}}{\mathrm{100}}}.\end{array}$$

Variances Var_G(Y) and Var_F(Y) due to parameter uncertainty in ground vegetation and fine root litter production model were derived from covariance matrices Σ_G and Σ_F and variances Var(α₄), Var(ϕ_p) and Var(μ_p) in a similar manner.

Variances in region-specific annual balance estimates due to uncertainty in litter input from living trees were estimated as

$$\begin{array}{}\text{(A25)}& {\mathrm{Var}}_L\left(Y_{rt}\right)=\sum _p{\left(A_{rpt}{\displaystyle \frac{\mathrm{44}\cdot \mathrm{0.5}}{\mathrm{12}}}\right)}^{\mathrm{2}}\mathrm{Var}\left(L_{rpt}\right),\end{array}$$

where variances Var(L_rpt) were based on RSEs in Table A7. When propagating these to the whole country, covariances $\mathrm{Cov}\left(L_{\mathrm{south},pt}{L}_{\mathrm{north},pt}\right)$, based on the NFI11 correlation 0.539 (Table A8), were accounted for, as were covariances between time points when considering the change estimates. For instance,

$$\begin{array}{}\text{(A26)}& {\mathrm{Var}}_L\left(Y_{\mathrm{2021}}-Y_{\mathrm{1990}}\right)={\mathit{c}}^{\prime }\mathbf{\Sigma }\mathit{c},\end{array}$$

where $c=(\mathrm{1},\mathrm{1},-\mathrm{1},-\mathrm{1}{)}^{\prime }$ and Σ is the covariance matrix of ${\left[Y_{\mathrm{south},\mathrm{2021}}{Y}_{\mathrm{north},\mathrm{2021}}{Y}_{\mathrm{south},\mathrm{1990}}{Y}_{\mathrm{south},\mathrm{1990}}\right]}^{\prime }$ derived using the correlation matrix of Table A8 and variances based on NFI11 for Y_r,2021 and NFI8 for Y_r,1990.

Finally, variances due to uncertainty in litter input from harvests and natural mortality,

$$\begin{array}{}\text{(A27)}& {\mathrm{Var}}_H\left(Y_{rt}\right)={\left({\displaystyle \frac{\mathrm{44}\cdot \mathrm{0.5}}{\mathrm{12}}}\sum _p{A}_{rpt}\right)}^{\mathrm{2}}\mathrm{Var}\left(H_{rt}\right),\end{array}$$

were based on RSEs in Table A7 and on the assumption that the input estimates are uncorrelated across regions and time points. To combine all sources of uncertainty, source-specific variances were simply added up:

$$\begin{array}{}\text{(A28)}& \begin{array}{rl}\mathrm{Var}\left(Y\right)& ={\mathrm{Var}}_A\left(Y\right)+{\mathrm{Var}}_B\left(Y\right)+{\mathrm{Var}}_R\left(Y\right)\\ & +{\mathrm{Var}}_{\mathrm{G}}\left(Y\right)+{\mathrm{Var}}_{\mathrm{F}}\left(Y\right)+{\mathrm{Var}}_L\left(Y\right)+{\mathrm{Var}}_H\left(Y\right).\end{array}\end{array}$$

Following IPCC guidelines (Frey et al. 2006), percentage uncertainties were defined as

$$\begin{array}{}\text{(A29)}& U_Y=\mathrm{100}\times \mathrm{1.96}\sqrt{\mathrm{Var}\left(Y\right)}.\end{array}$$

A4 Results

Uncertainties of the year 2021 (positive) balance estimates for southern Finland and the whole country are less than 100 % (Table A9), which indicates that the soils of drained organic peatlands are statistically significant sources of CO₂ according to our model predictions. However, the 2021 balance in northern Finland does not differ significantly from zero. The increase in emissions between the years 1990 and 2021 is significant in both regions and in the whole country (Table A10).

The greatest impact was due to uncertainty in the parameter estimates of the peat and litter decomposition model (41 %–67 % of variance in 2021 balance estimates, 73 %–78 % in change estimates; Tables A9 and A10). The impact of uncertainty in fine root model parameters was relatively modest (9 %–10 %) in change estimation, but somewhat greater (19 %–39 %) in annual estimates. From other sources of uncertainty, only litter input from living trees exceeded 10 % impact in any of the estimates. In particular, uncertainty in NFI estimates of site type areas and basal areas of trees had a negligible impact on annual balance estimates and contributed less than 5 % of the variance of change estimates.

The errors in annual balance estimates are positively correlated because the same model parameters with same estimation errors were applied throughout the series. For example,

$$\begin{array}{}\text{(A30)}& \begin{array}{rl}\mathrm{Corr}\left(Y_{\mathrm{2021}},Y_{\mathrm{1990}}\right)& ={\displaystyle \frac{\mathrm{Cov}\left(Y_{\mathrm{2021}},Y_{\mathrm{1990}}\right)}{\sqrt{\mathrm{Var}\left(Y_{\mathrm{2021}}\right)\mathrm{Var}\left(Y_{\mathrm{1990}}\right)}}}\\ & ={\displaystyle \frac{\begin{array}{c}\mathrm{Var}\left(Y_{\mathrm{2021}}\right)+\mathrm{Var}\left(Y_{\mathrm{1990}}\right)\\ -\mathrm{Var}\left(Y_{\mathrm{2021}}-Y_{\mathrm{1990}}\right)\end{array}}{\mathrm{2}\sqrt{\mathrm{Var}\left(Y_{\mathrm{2021}}\right)\mathrm{Var}\left(Y_{\mathrm{1990}}\right)}}}=\mathrm{0.9}.\end{array}\end{array}$$

Due to these correlations, change estimates can differ significantly from zero even if confidence intervals of annual estimates have a large overlap (Fig. A1).

Figure A1 Drained organic forest soil CO₂ balance for southern Finland, northern Finland and the whole country (black lines; equal to those in Fig. 4), as well as pointwise 95 % confidence intervals obtained by adding ± U to the balance estimate (grey lines).

Table A1 NFI12 estimates of the areas of drained peatland forest site types (FTYPEs with their Finnish abbreviations) in southern and northern Finland together with their standard errors (SEs) and relative standard errors (RSEs) due to sampling, assessed as explained in Korhonen et al. (2021, Supplement S1).

Table A2 NFI12 estimates of basal area of trees and their standard errors (SEs) and relative standard errors (RSEs) due to sampling, assessed as explained in Korhonen et al. (2021, Supplement S1).

Table A3 Covariance matrix Σ_R of the parameters of the peat and litter decomposition model derived from Ojanen et al. (2014, their Table A.5) after combining the Mtkg and Ptkg subtypes.

Table A4 Covariance matrix Σ_G of the parameters of the ground vegetation litter model derived from Ojanen et al. (2014, their Table A.4) after combining the Mtkg and Ptkg subtypes.

Table A5 Covariance matrix Σ_F of the parameters of the fine root biomass model of Ojanen et al. (2014, their Table A.2).

Table A6 The applied variances of site-type-specific fine root turnover rates, Var(ϕ_p) (expert judgment), and dwarf shrub coverages, Var(μ_p) (Ojanen et al. 2014, their Table A.3, after combining the Mtkg and Ptkg subtypes); the corresponding values ϕ_p and μ_p are given in Table 3.

Table A7 Relative standard errors (%) of litter production from living trees and from harvests and natural mortality on drained peatlands estimated from NFI11.

Table A8 Correlations of NFI estimates of litter production from living trees.

Table A9 Estimates of CO₂ release from peat and litter decomposition, net C inputs to soil converted to the units of CO₂, and soil CO₂ balance (“Net”) for the year 2021 together with the variance (“Var”) and uncertainty (“U”) of the estimates.

Table A10 Estimates of change from 1990 to 2021 in CO₂ release from peat and litter decomposition, net C inputs to soil converted to the units of CO₂, and soil CO₂ balance (“Net”) for the year 2021 together with the variance (“Var”) and uncertainty (“U”) of the change estimates.