Synthesis of greenhouse gas emission factors for forest organic soils in the Köppen–Geiger Dfb climate zone

Greenhouse gas (GHG) emissions from organic soils are a key component of land-use-related emissions, particularly in countries with large areas of organic soils. Temperate-zone forest soils remain less studied in GHG research than boreal soils. However, recent work has expanded coverage in the northeastern temperate region, which, under the Köppen–Geiger climate classification, shares key climatic characteristics with the southern boreal region (Dfb). This study synthesised updated GHG flux data to evaluate carbon balance and emissions from forest organic soils in the Dfb zone, stratified by drainage status, nutrient availability, and dominant tree species. Such stratification revealed CO₂ source-sink patterns, which encourage the ecological relevance of using these categories for data aggregation. The dominant tree species reflected nutrient status: drained coniferous and deciduous stands have been reported as CO₂ sources, emitting 0.03 ± 0.55 and 0.47 ± 0.29 t CO₂-C ha⁻¹ yr⁻¹, respectively, though soils tended to shift toward CO₂ sinks in stands older than 25 years. In contrast, undrained soils have generally been observed to function as CO₂ sinks, although not necessarily in all sites. However, this stratification was less informative for CH₄ and N₂O. CH₄ fluxes were primarily determined by water table level rather than by other site variables, whereas N₂O showed a tendency toward elevated emissions in deciduous stands, irrespective of drainage status.

1 Introduction

Organic soils, when predominantly peat-based, contain large stocks of organic carbon (C) and nitrogen (N), where widespread, are considered to be disproportionately large sources of land-use greenhouse gas (GHG) emissions (Lake et al., 2025). As a result, they have attracted growing attention from stakeholders in sustainable soil C management and GHG emission mitigation, thereby stimulating efforts to improve the accuracy of GHG emission estimates (e.g., Alm et al., 2023; Tiemeyer et al., 2020; Jauhiainen et al., 2023). The Intergovernmental Panel on Climate Change (IPCC) guidelines characterise emissions from organic soils based on stratification by boreal, temperate, and tropical climate zones (Eggleston et al., 2006). These zones are geographically broad, thus encompassing substantial variability in climate, vegetation, hydrology, nutrient availability, and geomorphology conditions, all of which can influence the magnitude of GHG emissions (He and Roulet, 2023). Compared with the extensive data available for the boreal zone, studies on GHG fluxes from organic soils in the temperate climate zone have been relatively scarce. For example, the most recent yet now outdated default IPCC emission factor (EF) for carbon dioxide (CO₂) from drained organic soils in temperate forest land was synthesised from only eight study sites, whereas the boreal EFs were based on data from 184 sites (Hiraishi et al., 2014). A more recent EF synthesis incorporated data published during the following decade, increasing the number of temperate sites to 35 and boreal sites to 218 (Jauhiainen et al., 2023). However, most temperate and boreal sites remained concentrated in Sweden and Finland, respectively. In the years following this synthesis, the number of study sites in the Baltic States, classified as part of the temperate zone according to the IPCC, has increased substantially, with directly measured organic-soil GHG emissions now reported from 56 forest sites (Butlers et al., 2025, 2022). This recent expansion has considerably improved data availability and geographical coverage for the northeastern part of the temperate zone.

The Baltic States, together with the southern parts of Finland and Sweden, fall within the Dfb climate zone (cold winters, no dry season, warm summers) according to the Köppen–Geiger climate classification (Beck et al., 2023). Given that the Dfb zone partially overlaps with the temperate zone and that most temperate-zone data fall within this overlap, using Dfb-specific data to elaborate EFs for the respective region is more appropriate than aggregating data of a much broader temperate zone. Combining data within narrower climatic ranges increases the spatial density of aggregated data, thereby supporting the development of EFs specific to additional geographical strata, such as forest types, drainage conditions, or nutrient status, thus enhancing the accuracy of soil emission upscaling.

In this study, we compiled existing annualised net CO₂, CH₄, and N₂O flux data from drained and undrained (naturally wet) forest land sites with organic soils across the Baltic States reported by Butlers et al. (2022, 2025) and Kamil-Sardar et al. (2026), supplemented with forest site data from the Köppen–Geiger Dfb climate zone aggregated by Jauhiainen et al. (2023). Our objective was to synthesise GHG EFs stratified by drainage status (drained, undrained), dominant tree species, and nutrient status, for the Dfb climate zone, to identify the most suitable grouping options for upscaling using Latvia as a case study.

2 Material and methods

2.1 Data compilation

Peer-reviewed annualised site-specific net GHG flux values, estimated using the chamber method, were compiled for the synthesis of organic soil CO₂, CH₄ and N₂O EFs (Hiraishi et al., 2014) for forest land. The addressed net CO₂ fluxes represent soil CO₂ balance: sum of soil CO₂-C loss by heterotrophic respiration and C input to soil by litter. The compiled dataset comprised study results from drained and undrained forest sites in the Baltic states (Butlers et al., 2022, 2023, 2025; Kamil-Sardar et al., 2026) supplemented by a subset of aggregated drained site values from Jauhiainen et al. (2023) overlapping with the Köppen–Geiger Dfb climate zone (cold winter, no dry season, warm summer, Beck et al., 2023), for which nutrient status was explicitly reported (Fig. 1). Dominant tree species in sites were coniferous, such as Norway spruce (Picea abies (L.) H.Karst., n=32) and Scots pine (Pinus sylvestris L., n=36), and deciduous species such as birch (Betula pendula Roth and Betula pubescens Ehrh., n=22) and alder (Alnus glutinosa (L.) Gaertn., n=18). We also included values assigned to other species, mixed stands, clearcuts, or cases where the tree species were not specified (n=25). Where tree species were not reported, the data were used only to derive tree species–independent EFs and were excluded from species-specific EF derivation. As a result, 297 annualised GHG flux values from 131 sites were compiled (Butlers, 2026). Fluxes can be traced to published peer-reviewed articles, 93 site-specific CO₂, CH₄ and N₂O fluxes were acquired by contacting article authors, which were validated against originally published summary results (Butlers et al., 2022, 2023). In addition to flux values, we compiled available site-specific data on quantitative soil properties, forest stand characteristics, and annual C inputs to soil by litter. Availability of these supplementary data was not an inclusion criterion. Study inclusion was based solely on the flux and site-characterisation criteria described above, including identification of drainage and soil nutrient status.

Figure 1 Locations of the sites in the Köppen–Geiger Dfb climate zone from which GHG flux data were collected and included in the study. The figure excludes three sites located in Canada.

Soil properties such as C, calcium (Ca), magnesium (Mg), phosphorus (P), and organic matter (OM) content were available for 34 % of sites, whereas N, pH, C:N ratio, and soil bulk density were available for 59 % of sites. Overall, at least some soil property data were available for 93 % of sites. Stand age was the most frequently reported stand parameter, available for 43 % of sites. Only 19 % of sites provided a complete set of stand variables (age, growing, stock, mean tree height, diameter at breast height, and basal area); these were also the sites for which annual C input to soil via litter (from above- and belowground ground vegetation, tree fine roots, and foliar litter) was available.

2.2 Data quality control and harmonisation

Quality control of the aggregated fluxes followed the approach described by Jauhiainen et al. (2023). This procedure assessed whether the originally reported peer-reviewed fluxes represented annual net soil GHG fluxes or those could be derived by applying supplementary procedures, including annualising warm-season measurements, separating heterotrophic from total soil respiration, or accounting for litter inputs. Where required, Jauhiainen et al. (2023) had already performed the necessary annualisation and harmonisation procedures and assigned weights to individual estimates according to their relative reliability. A weight of 1.0 was assigned to estimates without major identified limitations. In contrast, a weight of 0.5 was assigned where uncertainty was higher, for example because of limited spatial replication or greater reliance on supplementary data when deriving the annual flux estimate. Jauhiainen et al. (2023) provide detailed descriptions of adjustments applied to individual flux estimates.

In this synthesis, we used the estimates compiled by Jauhiainen et al. (2023) directly and retained their original weights. For additional fluxes compiled in this study, we applied the same quality-assessment criteria. We verified that the originally reported values represented annualised net soil GHG fluxes suitable for direct comparison with the Jauhiainen et al. (2023) dataset. Because no additional data manipulation was required, we assigned a weight of 1.0 in all cases.

Before further processing of the compiled dataset, we ensured that the units were harmonised. Additional harmonisation was needed for classifying soil nutrient status. The assigned soil nutrient status (rich, moderate) was harmonised by aligning the classification schemes used across the datasets. In the dataset of Jauhiainen et al. (2023), sites labelled as nutrient-poor were combined with the Baltic data classified as Myrtillosa turf. mel. forest site types and grouped under the “moderate” nutrient class. Similarly, sites labelled as nutrient-rich in Jauhiainen et al. (2023) and those classified as Oxalidosa turf. mel. and Dryopterioso–caricosa forest types in the Baltic data were grouped under the “rich” class.

2.3 Emission factor synthesis and statistical analysis

Using the assigned weights, for each EF aggregation level, we calculated the weighted mean with the associated standard error (SE) and 95 % confidence interval (CI). In addition, weighted medians with weighted 25th and 75th percentiles (IQR) are presented. In the text, uncertainty is expressed as the SE.

Differences in fluxes, grouped by drainage status, dominant tree species or type (coniferous, deciduous), and nutrient status, were assessed using Kruskal–Wallis tests. Mixed-effects models were used to evaluate the influence of species and nutrient status on GHG emissions while accounting for site-level random effects (R Core Team, 2023). These statistical methods were also applied to flux metadata, such as soil physico-chemical properties, to further strengthen the ecological rationale for data grouping aimed at elaborating EFs suitable for upscaling.

3 Results

3.1 Soil CO₂ fluxes

No significant differences in soil annual net CO₂ fluxes were detected among dominant tree species or between coniferous and deciduous species type groups when stratified by drainage and nutrient status (p=0.32–0.84). Nutrient status had somewhat greater predictive power than species type: in mixed-effects models, nutrient status explained ∼ 13 % of the variance in net CO₂ fluxes (marginal R²=0.13, χ²=11.6, p<0.001), while species type accounted for only ∼ 3 % within nutrient-rich sites (marginal R²=0.03, χ²=1.9, p=0.17). However, robust comparison was constrained by the data structure, as deciduous stands occurred exclusively on nutrient-rich sites. Nevertheless, the data demonstrated the dominant tree species type reflecting nutrient availability, as notable differences between coniferous and deciduous stands were observed for Ca, $\mathrm{C}/\mathrm{N}$ ratio, K, Mg, N, P, and pH (Fig. A1). These observations indicated that both soil nutrient status and dominant tree species can be used to stratify data for EF calculation, depending on patterns observed in flux data.

Overall, mean net CO₂ emissions were 0.22 ± 0.18 t CO₂-C ha⁻¹ yr⁻¹ in drained sites and −0.60 ± 0.37 t CO₂-C ha⁻¹ yr⁻¹ in undrained sites. The undrained sites were represented by nutrient-rich stands of Norway spruce, birch, and alder, with the highest mean removals observed in birch stands (−1.20 ± 0.76 t CO₂-C ha⁻¹ yr⁻¹). However, undrained soils at some sites across all dominant species categories have also acted as CO₂ sources (Fig. A2). The mean net CO₂ flux of drained soils in coniferous and deciduous stands was 0.03 ± 0.55 and 0.47±29 t CO₂-C ha⁻¹ yr⁻¹, respectively. Scots pine and Norway spruce were equally represented within the drained coniferous category (Fig. 2), and statistical analysis indicated that differences between fluxes of these categories are not significant. However, a tendency toward lower fluxes, with a mean of −1.13 ± 0.52 t CO₂-C ha⁻¹ yr⁻¹, was observed in nutrient-moderate sites (Fig. 2), which were dominated by Scots pine (80 % of observations). In contrast, drained nutrient-rich coniferous sites dominated by Norway spruce (67 % of observations) were a source of 0.55 ± 0.26 t CO₂-C ha⁻¹ yr⁻¹, similar to what was observed in drained deciduous stands, which were present in nutrient-rich sites only. Consequently, drained nutrient-rich soils were, on average, a CO₂ source of 0.51 ± 0.17 t CO₂-C ha⁻¹ yr⁻¹. Drained alder stands appeared as CO₂ sinks; however, their fluxes did not differ statistically from other nutrient-rich sites.

Figure 2 Mean soil annual net CO₂ fluxes by soil drainage and nutrient status, and dominant tree species. Error bars indicate a 95 % confidence interval; n indicates the number of sites represented. A common lowercase letter indicates statistically insignificant differences in average annual net CO₂ fluxes. For numerical values, see Table B1.

Among stand variables, age showed the strongest relationship with net CO₂ flux in drained soils, indicating a tendency for soils to shift from a CO₂ source to a sink beyond approximately 25 years of stand age (Fig. 3). No relationship was observed between gross CO₂ emissions and stand variables or annual litter. Gross CO₂ flux, however, showed a weak decreasing trend with decreasing soil pH and N concentration, although the correlations were too weak to quantify. Net CO₂ fluxes showed no consistent relationship with water table level, as both sinks and sources occurred across the range of water table conditions. Annual cumulative C input to soil from tree fine roots and foliar litterfall increased linearly with stand age (slope = 0.046, intercept = 0, R²=0.9). No relationship with fluxes was found for above- or belowground vegetation, likely due to limited data availability.

Figure 3 Variation in soil annual net CO₂ fluxes depending on stand age (all sites with known stand age combined).

3.2 Soil CH₄ fluxes

CH₄ fluxes showed less pronounced effects of nutrient richness and an apparent influence of species (Fig. A3). In drained sites, soils acted as CH₄ sources in Scots pine stands (6.3 ± 2.6 kg CH₄ ha⁻¹ yr⁻¹) and as CH₄ sinks in Norway spruce stands (₋4.5 ± 0.8 kg CH₄ ha⁻¹ yr⁻¹), irrespective of nutrient status. Soils in drained, nutrient-rich birch stands were approximately CH₄-neutral, while alder stands showed on average higher and more variable CH₄ emissions under both drained and undrained conditions (Fig. 4). In undrained Norway spruce and birch sites, mean CH₄ emissions were 9.3 ± 8.0 and 23.3 ± 18.0 kg CH₄ ha⁻¹ yr⁻¹, respectively. However, due to the CH₄ emission hotspot in alder stands, mean undrained-soil CH₄ emissions reached 92.88 ± 78.05 kg CH₄ ha⁻¹ yr⁻¹.

Figure 4 Mean soil annual CH₄ fluxes by soil drainage and nutrient status, and dominant tree species. Error bars indicate a 95 % confidence interval; n indicates the number of sites represented. A common lowercase letter indicates statistically insignificant differences in average annual CH₄ fluxes. For numerical values, see Table B2.

However, these patterns reflect variation in water table levels among the respective study sites rather than species effects. For instance, alder stands in drained sites appear as comparatively larger CH₄ sources because their mean water table level was 16 cm below the surface. Moreover, when the alder emission hotspot in undrained sites is excluded, mean CH₄ emissions decrease from 178.2 to 12.8 kg CH₄ ha⁻¹ yr⁻¹. Similarly, Scots pine stands appear biased toward CH₄ sources, as half of the sites have mean water table levels above ∼ 20 cm. In the compiled dataset, a linear regression between mean water table level and annualised CH₄ fluxes identified −20 cm as a threshold separating source and sink conditions, with CH₄ uptake increasing by 0.102 kg CH₄ ha⁻¹ yr⁻¹ for each 1 cm decrease in water table level. The mean CH₄ fluxes for sites with mean water table level below and above 20 cm were −3.91 ± 0.41 and 70.88 ± 44.30 kg CH₄ ha⁻¹ yr⁻¹, respectively. For comparison, drained and undrained sites were on average sources of 1.17 ± 1.58 and 92.88 ± 78.05 kg CH₄ ha⁻¹ yr⁻¹.

3.3 Soil N₂O fluxes

Soil N₂O flux patterns suggested that species type may influence flux behaviour, as excess emissions were observed in deciduous stands regardless of drainage conditions, though these excesses were more pronounced under drained conditions (Figs. A4, 5). On average, N₂O fluxes in deciduous stands were 2.82 ± 0.68 kg N₂O ha⁻¹ yr⁻¹, irrespective of nutrient or drainage status. Due to the presence of outlier (excess) fluxes, statistical differences among deciduous stands could not be identified. Mean fluxes in drained deciduous sites were 5.58 ± 1.87 kg N₂O ha⁻¹ yr⁻¹, whereas in undrained sites they were 3.51 ± 1.35 kg N₂O ha⁻¹ yr⁻¹ or on average 37 % lower. In all conifer stands, emissions were on average 1.58 ± 0.35 kg N₂O ha⁻¹ yr⁻¹.

Figure 5 Mean soil annual N₂O fluxes by soil drainage and nutrient status, and dominant tree species. Error bars indicate a 95 % confidence interval; n indicates the number of sites represented. A common lowercase letter indicates statistically insignificant differences in average annual N₂O fluxes. For numerical values, see Table B3.

3.4 Upscaling: Latvian Case

The synthesised EFs are intended for application within the Dfb climate zone; however, the level of EF aggregation can be adjusted for specific regions to reflect local conditions, such as nutrient status and the presence of dominant tree species. Here, we demonstrate the application of the synthesised EFs for upscaling tailored to Latvia's specific conditions. Latvia is situated in the hemiboreal region according to the vegetation zone classification (temperate–boreal transition; Ahti et al., 1968) and in the IPCC-defined cold temperate, moist climate region (Eggleston et al., 2006), which overlaps with the Köppen–Geiger Dfb classification (Beck et al., 2023). Therefore, the synthesised EFs are applicable, although specific Latvian conditions must be taken into account.

In Latvia, the area of drained organic soils is estimated at 529 kha, representing 8 % of the state area. Around 314 kha or 59 % of drained organic soils are in forest land (Fig. A5). At the Latvia scale, the coarsest level of granularity that can be reasonably applied for emission upscaling corresponds to the heterogeneity and spatial distribution of forest site types and dominant tree species as provided by the National Forest Inventory (NFI) data (Fig. 6). NFI distinguishes eight forest site types with organic soils based on drainage status and nutrient availability, the latter classified by ground-vegetation composition (Bušs, 1981). Nutrient availability for forest site types is further grouped into three categories: nutrient-poor, nutrient-moderate, and nutrient-rich soils (Bušs, 1981). For emission upscaling, we combined nutrient-poor and nutrient-moderate areas and applied the synthesised EFs for the “moderate” category.

Figure 6 Distribution of soil nutrient status and dominant tree species across drained forest lands in Latvia (Latvian National Forest Inventory, https://www.silava.lv/images/Petijumi/Nacionalais-meza-monitorings/MRM-rezultati/Kopsavilkumi-2023-IV-cikls.xlsx, last access: 12 January 2026).

Although differences in drained soil CO₂ EFs, stratified by nutrient status and tree species, were generally not statistically significant, we applied EFs based on distinct source-sink patterns. To characterise emissions in all drained Scots pine stands, we used EFs aggregated for Scots pine under moderate and rich nutrient status. The choice was based on ecological reasoning that Scots pine does not tend to dominate on nutrient-rich soils, which was supported by NFI data showing that Scots pine accounts for less than 0.1 % of nutrient-rich forest types. This showed that it's appropriate to apply Scots pine EF, which aggregates any assigned nutrient status. A similar issue was encountered for Norway spruce, which typically occurs on nutrient-moderate soils (Latvian National Forest Inventory). However, the aggregated EF for such sites was underrepresented (3 sites, 20 %) and exhibited high flux variability (−1.91 ± 1.89 t CO₂-C ha⁻¹ yr⁻¹). Therefore, we considered it appropriate to combine all available data points and derive a single EF for all drained Norway spruce stands, irrespective of nutrient status. For deciduous stands, only nutrient-rich EFs were available. Therefore, the choice was between applying a single EF to all deciduous stands or using birch- and alder-specific EFs for their respective areas. In addition to birch and alder (grey alder, black alder), approximately 6 % of organic soils are overgrown with other deciduous species, generally aspen, oak, or ash. For less common species for which national allometric equations are unavailable, the NFI classifies these stands as birch; we followed the same approach in this study. For alder species, we applied an alder-specific EF, as it exhibited a distinct source–sink pattern compared with birch.

For undrained soils, as well as for CH₄ and N₂O when grouped by drainage status, the use of a single EF was considered most appropriate. An exception was made for drained soil N₂O, for which species–type–specific EFs were deemed more appropriate. These decisions were based on the absence of statistically significant differences in EFs or in source–sink dynamics, or on the observation that differences were driven by factors other than species or nutrient status, as demonstrated in the section on CH₄ emissions.

Using this approach, the weighted CO₂ EF for drained organic soils in Latvia is 0.14 t CO₂-C ha⁻¹ yr⁻¹ (Table 1). Applying a single EF to all deciduous stands increases the weighted value to 0.23 t CO₂-C ha⁻¹ yr⁻¹. Conversely, if a single nutrient-based EF were applied for CO₂, the area-weighted mean would be −0.71 t CO₂-C ha⁻¹ yr⁻¹, due to nutrient-rich organic soils being underweighted in Latvia's area distribution. Similarly, if, for N₂O rather than species-type-specific EFs, a single value was used, the weighted 3.80 kg N₂O ha⁻¹ yr⁻¹ would be reduced to 2.82 kg N₂O ha⁻¹ yr⁻¹. Consequently, the choice of stratification can lead to differences in country-specific EFs.

Table 1 Emission factors were identified as appropriate for upscaling GHG emissions from drained forest organic soil in Latvia.

¹ All Scots pine sites. ² All Norway spruce sites. ³ All nutrient-rich deciduous sites. ⁴ All sites. ⁵ All conifer sites.

4 Discussion

4.1 Emission factor stratification improves quality of the estimates

The main reason for developing EFs that reflect different local conditions, such as soil nutrient status, dominant tree species, or other factors, is their use in emission upscaling with improved accuracy. Applying a single EF across a wide geographic scale inevitably introduces bias if the EF does not capture the heterogeneity of conditions in the area. However, this is common practice in national GHG inventories (Lake et al., 2025), where EFs applied often represent averages derived from datasets covering broad geographic areas, without accounting for the representativeness or composition of the underlying data. These averaged EFs are subsequently applied to areas that may differ substantially in their actual composition in terms of dominant tree species and site nutrient status. Increasing the granularity of EFs represents an effective way to enhance the precision of upscaling. However, this also means that segregated EFs must be sufficiently representative of the corresponding area strata, necessitating more observations. Our study demonstrates an approach to mitigate the restraints caused by a limited number of observations. We constructed EFs for the northeastern part of the temperate zone by including a proportion of observations from the boreal zone, thereby ensuring partial overlap of data sources across regional boundaries. This overlap was justified by the demonstrated comparability of the regions, specifically by the Dfb climate zone, and enabled the acquisition of more observations without compromising representativeness. Furthermore, regional datasets are often produced by different research teams using varying study designs and equipment, which may introduce systematic biases into the results. Incorporating a portion of data from neighbouring areas into local datasets likely reduced potential biases arising from such methodological differences in EF estimation. This, in turn, technically improves the comparability of the elaborated EFs with those from the boreal zone, an essential requirement for large-scale GHG inventories that support coordinated policy decisions, such as those at the EU level. Notably, Jauhiainen et al. (2023) demonstrated that expanding and stratifying the dataset into more specific site categories improves EF precision, even if the updated mean emission values remain similar to previous estimates. Despite the uncertainties, the magnitude and direction of fluxes estimated in our study are broadly consistent with published measurements from boreal and hemiboreal forest organic soils (Jauhiainen et al., 2023), and the synthesised EFs reflect the best available knowledge for the region. In the absence of country-specific EFs, Latvia would have to rely on the IPCC default EF of 2.60 t CO₂-C ha⁻¹ yr⁻¹, based on observations from eight sites (Hiraishi et al., 2014). Alternatively, regionally refined EFs could be adopted, such as the values of 1.61 and 1.46 t CO₂-C ha⁻¹ yr⁻¹ reported by Jauhiainen et al. (2023) for nutrient-rich and nutrient-poor temperate forest organic soils, respectively, based on observations from 19 sites. However, the temperate climate region covers a broad geographical and climatic range extending to considerably warmer and drier parts of Europe, making these aggregated EFs poorly representative of Latvian conditions. The EF synthesised in the present study, 0.14 t CO₂-C ha⁻¹ yr⁻¹, is substantially lower and is more consistent with the EFs derived by Jauhiainen et al. (2023) for drained boreal forest organic soils, which range from −0.23 to 0.73 t CO₂-C ha⁻¹ yr⁻¹ depending on site nutrient status and stand productivity. The refined EF developed in this study, based on data from 81 study sites, also represents an improvement over the EF of 0.52 t CO₂-C ha⁻¹ yr⁻¹ currently applied in Latvia's national GHG inventory (Latvia's National Inventory Document, 2025).

4.2 Representativeness and selection of emission factor strata

We observed that EF aggregation and the level of upscaling granularity can lead to diverging area-weighted EF outcomes, and identifying the most optimal approach for a given region may not be straightforward. Combining datasets from multiple sources introduces challenges related to data interpretation and harmonisation. In our case, an issue arose regarding potential differences in how soil nutrient-richness categories were assigned across the source datasets, likely reflecting variations in national forest site-type classification systems. This issue became evident in our upscaling case study, as the compiled dataset contained a substantial number of coniferous stands classified as growing on nutrient-rich soils, whereas such stands account for only about 4 % of drained forest organic soils in Latvia (Latvian National Forest Inventory). Nine of ten data points for nutrient-rich Scots pine sites were obtained from the Jauhiainen et al. (2023) synthesis, suggesting that different forest site-type classification systems may not yield directly comparable indicators of nutrient status. Accordingly, nutrient-status classes should be interpreted as operational groupings rather than strictly equivalent indicators across datasets.

Despite this issue, the results suggested that nutrient availability has a stronger predictive capacity for fluxes than tree species. This finding supports previous results indicating that nutrient status plays a key role in determining GHG emissions across different ecosystems (e.g., Säurich et al., 2019). However, in our case, aggregated EFs derived specifically for conifers on nutrient-rich soils had limited practical relevance for national upscaling. Therefore, we considered flux aggregation by species to be the most appropriate approach, as nutrient status classifications in the compiled data may vary, but stand types themselves are ecologically meaningful for data stratification. This is because dominant tree species not only reflect soil properties but also shape them (Hagen-Thorn et al., 2004; Lorenz and Thiele-Bruhn, 2019; Vesterdal et al., 2013). Fast-growing, nutrient-demanding species likely represent nutrient-rich soils because in such areas they tend to outcompete other species, whereas nutrient-tolerant species prevail in nutrient-poor sites, as they prevail in conditions of low nutrient availability that limits the growth of other species. In turn, the species established further reinforces the soil conditions, including decomposition of litter and the associated microbial community composition they inhabit (Wang et al., 2016; Bai et al., 2023): nutrient-demanding broadleaves enrich soils with labile organic matter (Vancampenhout et al., 2012), while nutrient-tolerant species maintain more acidic conditions that limit nutrient availability for other species (Šantrůčková et al., 2019). This reasoning was supported by the observation that nutrient availability differs statistically among sites dominated by different species types (Fig. A1).

The decision to stratify nutrient-rich alder and birch stands, which showed an apparently distinct sink–source pattern though not statistically significant, was more challenging. Whether upscaling should rely on such stratification or on a total mean EF remains debatable, as current knowledge lacks ecological evidence to explain these sink–source dynamics. Nutrient-rich undrained stands dominated by black alder in our dataset showed increased spatial variability in CH₄ and N₂O fluxes, complicating generalisations. High inputs of labile litter and nitrogen in deciduous stands could stimulate microbial processes that produce N₂O (Zhou et al., 2022) even under wet conditions (Elberling et al., 2023; Rautakoski et al., 2024), but more data are needed to confirm this mechanism. Such difficulties will likely persist until a sufficient number of observations is available to provide a statistically robust justification for data-use decisions, or may persist even as long as EFs are used for emission upscaling. Studies have begun linking organic soil GHG fluxes to measurable site properties. For instance, under both drained and wet conditions, the magnitudes of peat CO₂ and CH₄ emissions are linked to geochemical variables such as pH, P content, and water-holding capacity (Nielsen et al., 2023). Since soil emissions are dependent on the chemical and physical properties of the soil (Rubaiyat et al., 2023), improvements can be achieved not only by expanding the emission sample set accompanied by observations of potential influencing factors, including quantitative assessments of soil properties, but also by conducting in-depth evaluations of the relationships between soil properties and forest types. Quantifying soil–forest type–species relationships would allow us to determine which level of data stratification or aggregation is ecologically rational.

Irrespective of the chosen flux stratification, application of acquired average values for upscaling to national emissions is likely to introduce biases, if the weighted impact of site characteristics, including stand variables, water table level and management practices are not accounted for. Water table levels represented at the study sites clearly demonstrated how bias can be introduced into the synthesised CH₄ EFs, creating a false impression of species impact. A water table level threshold of approximately 20 cm is well known to be a deterministic indicator of whether a site is a CH₄ source or sink (Jovani et al., 2023; Hatano, 2019). We observed that the significant differences in fluxes among drained sites with different dominant tree species were driven by the inclusion of sites with mean water tables above −20 cm. This highlights a common issue: the classification of drainage status (i.e. the presence of drainage ditches) does not necessarily imply that the ditches are functional or that the water table is actually lowered (Butlers et al., 2025). For this reason, we do not consider dominant tree species a suitable single proxy for CH₄ flux upscaling; instead, it should rely on water table data. However, because spatial water-table data, even when derived from remote sensing, have not yet been validated against ground-truth measurements (Ibrahim et al., 2024), it is more practical to use mean fluxes by drainage status for upscaling. This approach would assume that the probability of non-functional drainage ditches within the studied sites reflects the situation across the upscaled region.

4.3 Effects of stand development on emission factors

Stand development stages included in the dataset demonstrated an impact on synthesised CO₂ EFs. The issue is well illustrated by clearcuts, where the risk of soil C loss is increased, with mean fluxes of 0.90 ± 0.28 t CO₂-C ha⁻¹ yr⁻¹, compared with the overall mean for drained soils of 0.22 ± 0.18 t CO₂-C ha⁻¹ yr⁻¹. Increased GHG emissions in a drained peatland forest after clear-cutting, at least in the short term, compared with mature forests, are consistent with earlier findings (e.g., Tikkasalo et al., 2025; Korkiakoski et al., 2023). However, such assessments are challenged by the difficulty of partitioning CO₂ emissions originating from the decomposition of historical soil organic matter and deadwood, which increase substantially after harvesting due to above- and belowground residues (Menichetti et al., 2025). The compiled dataset also showed a moderate correlation (r=0.41) between stand age and soil C balance, indicating a higher likelihood of soils acting as a C source in stands younger than ∼ 25 years and suggesting an increasing sink potential of about 0.02 t C ha⁻¹ yr⁻¹ with subsequent aging. The flux data supplemented by stand age information showed a mean age of 56 years and an interquartile range of 31–73 years. This indicates a risk of underestimating emissions if the elevated fluxes from clearcuts and young stands are not properly accounted for. Using a 25-year threshold derived from regression analysis, we found that younger stands, on average, were a source of 0.61 ± 0.44 t CO₂-C ha⁻¹ yr⁻¹.

However, the soil carbon balance relationship identified with stand age only illustrates part of the issue. The true impact of stand development should be derived from relationships between soil conditions and stand variables, such as soil gross emissions, carbon inputs from litter, and their interaction. We observed that the ground vegetation component contributes approximately 30 % of the total reported carbon input to soil. However, no quantitative relationships with site variables could be identified. This component is likely highly influenced by soil nutrient status, suggesting that potential relationships are obscured by noise arising from insufficient measurement replication and unaccounted site heterogeneity. It is crucial that future research focuses on increasing the amount of data characterising soil C input dynamics during stand development to enable dynamic soil C balance modelling. At present, however, such analyses are highly constrained by the critically limited number of available observations. Similar to GHG efflux data, regional C influx measurements are largely concentrated in Finland, where biomass–litter equations dependent on stand variables have been developed (e.g., Starr et al., 2005; Muukkonen and Mäkipää, 2006). These equations, however, are unlikely to be directly applicable at a broader regional scale, making this data gap a key limitation for advancing spatially explicit soil C balance modelling. Furthermore, C fluxes, both efflux and influx, can respond to typical management activities such as tending (e.g., Yang et al., 2022; Doukalianou et al., 2022) as well as natural disturbances (e.g., Ribeiro-Kumara et al., 2020) and therefore vary over the course of a forest rotation cycle. However, data remain insufficient to accurately incorporate these impacts. Nevertheless, the compiled dataset indicates that stratification by drainage status and dominant tree species provides a practical first-order improvement for upscaling.

4.4 Drainage impact on emission factors

EFs derived from combined drained and undrained site observations show that, on average, drained organic soils are either GHG-neutral or net sources, while undrained soils can function as both sources and sinks, with uncertainties several times higher due to variable CH₄ and N₂O fluxes. Consequently, net drainage impact on GHG emissions remains unknown (Fig. 7). Uncertainty in undrained soil GHG emissions remains largely driven by a single CH₄-emission hotspot site. It is unknown whether the alder CH₄ emission hotspot is exceptional, and increasing the amount of data will continue to reduce the mean estimate and uncertainty, or another hotspot will eventually be captured.

Figure 7 Mean GHG emissions from drained and undrained organic soils expressed in CO₂ equivalents according to AR5 100-year global warming potentials (Huang et al., 2013). Error bars indicate the combined uncertainty (root mean square sum) expressed as a 95 % confidence interval.

5 Conclusions

Undrained forest organic soils generally skew toward acting as CO₂ sinks; however, irrespective of drainage status, soils have been found to be CO₂ sources, likely depending on local conditions such as soil nutritional conditions, dominant tree species, and stand development stage. Soil nutrient status, to a large extent, was reflected by the dominant tree species, indicating that both can serve as useful categories for classifying soil GHG fluxes. Drained deciduous stands are reported to be on average CO₂ sources, while drained coniferous stands are roughly CO₂ neutral. Soils in drained stands older than 25 years are more likely to act as CO₂ sinks, highlighting the importance of the forest development stage in carbon balance assessment. CH₄ emissions from undrained soils exhibited the greatest uncertainty, largely driven by variability in groundwater levels rather than forest type or nutritional conditions. This makes the development and upscaling of CH₄ emission factors especially challenging. Similarly, N₂O emissions showed high site-level variability, with elevated fluxes driven by statistical outliers in both drained and undrained deciduous stands, resulting in comparably increased fluxes. Currently available data do not allow robust quantification of drainage impacts on the total GHG balance, owing to residual flux variability that cannot yet be statistically explained. However, expanding CO₂ datasets linked to stand variables and quantitative soil physicochemical properties offers strong potential for future site-specific assessments of soil carbon balance across stand development stages. The study's findings provide empirical support for refining greenhouse gas accounting approaches for forest organic soils.

Appendix A: Appendix figures

Figure A1 Variation in soil parameters. (Ca – calcium, Mg – magnesium, K – potassium, N – total nitrogen, P – total phosphorus, CN – $\mathrm{C}/\mathrm{N}$ ratio, pH – soil pH) in drained coniferous (red) and deciduous (green) stands. Ca, K, Mg, N and P are expressed in g kg⁻¹. The bottom and top edges of the box represent the 25th and 75th percentiles, summarising the interquartile range (IQR). The whiskers extend to the smallest and largest values within 1.5 × IQR from the 25th and 75th percentiles, respectively. A solid horizontal line in the box indicates the mean and median values.

Figure A2 Variation of soil net CO₂ fluxes by drainage status, tree species type and nutrient status. The bottom and top edges of the box represent the 25th and 75th percentiles, summarising the interquartile range (IQR). The whiskers extend to the smallest and largest values within 1.5 × IQR from the 25th and 75th percentiles, respectively. A solid horizontal line in the box indicates the mean and median values.

Figure A3 Variation in soil CH₄ fluxes by drainage status, tree species type and nutrient status. The bottom and top edges of the box represent the 25th and 75th percentiles, summarising the interquartile range (IQR). The whiskers extend to the smallest and largest values within 1.5 × IQR from the 25th and 75th percentiles, respectively. A solid horizontal line in the box indicates the mean and median values.

Figure A4 Variation in soil N₂O fluxes by drainage status, tree species type, and nutrient status. The bottom and top edges of the box represent the 25th and 75th percentiles, summarising the interquartile range (IQR). The whiskers extend to the smallest and largest values within 1.5 × IQR from the 25th and 75th percentiles, respectively. A solid horizontal line in the box indicates the mean and median values.

Figure A5 Area of drained organic soils (left panel) and reported on-site GHG emissions (right panel) from drained organic soils across different land use categories (Latvia's National Inventory Document, 2025).

Appendix B: Appendix tables

Table B1 Summary statistics of CO₂ emissions(t CO₂-C ha⁻¹ yr⁻¹) stratified by soil nutrient status and dominant tree species. SE – standard error, CI – 95 % confidence interval, Var – variance, IQR – interquartile range.

Table B2 Summary statistics of CH₄ emissions (kg CH₄ ha⁻¹ yr⁻¹) stratified by soil nutrient status and dominant tree species. SE – standard error, CI – 95 % confidence interval, Var – variance, IQR – interquartile range.

Table B3 Summary statistics of N₂O emissions (kg N₂O ha⁻¹ yr⁻¹) stratified by soil nutrient status and dominant tree species. SE – standard error, CI – 95 % confidence interval, Var – variance, IQR – interquartile range.