Impact of Former Peat Extraction Field Afforestation on Soil Greenhouse Gas Emissions in Hemiboreal Region

Abstract

The reduction of greenhouse gas (GHG) emissions and climate change mitigation are global issues. Peatlands in Europe are widely distributed in the Nordic–Baltic region, and Baltic countries are some of the largest peat suppliers for horticulture in Europe. However, there is no sustainable substitute for peat in the horticulture industry. Therefore, it is necessary to identify suitable re-cultivation types for former peat extraction fields, because knowledge about the effect of re-cultivation on annual carbon and GHG budgets is limited. Ecosystem GHG (CO₂, CH₄, N₂O) exchange measurements, environmental parameter assessment and sampling in the study were conducted in a hemiboreal vegetation zone for 24 consecutive months in former peat extraction fields with different re-cultivation management strategies (land use types). The aim of the study was to assess the influence of diverse re-cultivation management strategies on the GHG emissions of former peat extraction fields. The most suitable re-cultivation management is afforestation with Scots pine (Pinus sylvestris) in order to obtain the lowest annual CO₂eq values and ensure additional carbon sequestration in living tree biomass. The developed linear mixed-effect models showed a good model fit (R²CO₂ = 0.80, R²CH₄ = 0.74) for the analyzed land use types, and thus can be used for CO₂ and CH₄ emissions estimation.

1. Introduction

Climate change and the global increase in greenhouse gas (GHG) concentrations in the atmosphere are a worldwide environmental issue [1]. Legal frameworks to mitigate climate change have been developed internationally (UNFCC, Kyoto Protocol, Paris Agreement) [2,3], regionally (EU Green Deal) [4] and nationally [5]. The main target of the regulations is to reduce the rise in air temperature well below 2 °C and to reduce GHG emissions in the future. Furthermore, based on these regulations, climate neutrality must be achieved by 2050. Every year, the involved countries report GHG emissions and removals from the Land Use, Land Use Change and Forestry (LULUCF) sector. Countries without regionally specific emission factors determine the impact on climate by using default emission factors from the IPCC guidelines [6,7]. Moreover, from 2026, the European Parliament requires countries to include wetland management (including peat extraction) in the LULUCF accounting. Therefore, practical information about the sustainable and climate-smart management of former peat extraction fields is a necessity.

Peatlands across European Union member states (EU27) cover roughly 284,000 km² [8], of which 33,000 km² are protected territories under Natura 2000. The peat soil distribution of Europe indicates a strong northern bias, where most of the peat soils are distributed in the Nordic–Baltic region, and Baltic countries are some of the largest peat suppliers for horticulture in Europe [9]. The target of the Green Deal is to reduce GHG emissions; therefore, aims for peat conservation and soil protection are discussed at the EU level, because peat is classified as a fossil resource and GHG emissions arise from the extracted peat itself (off-site emissions during its decomposition or combustion). However, peatlands are of high economic importance, because the horticulture industry and forestry are still dependent on peat, as no sustainable alternative material to peat substrates for seedlings is yet available for complete substitution. The progress of the peat industry towards climate neutrality, promoting research and innovation, as well as the re-cultivation of territories, is one of the goals of climate change mitigation targets. In order to fulfil the obligations of the LULUCF sector, it is necessary to identify suitable types for the re-cultivation of peat extraction fields and to use them in order to reduce the negative carbon footprint and achieve climate neutrality goals.

An often mentioned peatland restoration measure is the rewetting of drained peatlands by closing drainage systems to restore the water table to previous conditions [10]. Furthermore, peatlands in Europe are often converted into agricultural or forestry land [11]. The afforestation of former peat deposits is another type of re-cultivation that reduces GHG emissions and ensures carbon storage in tree biomass, compared to maintaining former peat extraction fields in their current state. Forests and soils (especially organic soils) are significant carbon pools [12,13,14]; therefore, they can play an important role in climate change mitigation targets. Organic soils can both remove and emit GHG, thus contributing to the global atmospheric GHG concentration, and have a huge impact on the reported GHG emission levels in the LULUCF sector, including emissions and CO₂ removals from wetlands, forest lands, croplands, grasslands and agricultural lands [15,16]. Whether organic soil with a regulated groundwater level is a net source or a sink of CO₂ emissions is determined by the site fertility, dominant tree species and depth of drained soil layer [17,18]. In the context of achieving climate goals, it is essential to carry out research within which it is possible to determine the most suitable re-cultivation practice of former peat extraction fields to reduce GHG emissions and favor climate mitigation goals. Therefore, the aim of the study is to assess the influence of diverse re-cultivation management strategies on the GHG emissions of former peat extraction fields.

2. Materials and Methods

The study was carried out in former peat extraction fields located in a hemiboreal vegetation zone in Latvia [19]. The study was conducted for 24 consecutive months, from December 2016 to November 2018, in peat extraction fields located across the territory of Latvia. In total, twelve fields (objects) with four different management regimes (land use types) after peat extraction (

Table 1. Study objects by land use type and management approach.

Land Use	Object Code	Objects	WGS84 Coordinates
			Lat.	Long.
Peat extraction field	Control	Cena mire	56.83119	23.99491
		Kaigu mire	56.71854	23.57486
		Silgulda mire	57.32854	27.39501
Abandoned peat extraction field with vegetation of herbs and dwarf shrubs	Vegetation	Silgulda mire	57.31271	27.39117
		Cena mire	56.82605	23.97721
		Cepļa mire	57.22004	26.47574
Afforested area (pine stand)	Pinus	Cepļa mire	57.21654	26.47937
		Kaigu mire	56.74538	23.60054
		Lambārte mire	56.50447	24.31264
Afforested area (birch stand)	Betula	Lielsala mire	57.35495	22.32452
		Silgulda mire	57.31000	27.40589
		Plece mire	56.72750	21.51994

) were studied. In each land use type, three objects were selected and one sample plot (R = 12.62 m) was established in each.

The first group comprised peat extraction fields (Control) that had been left intact after peat excavation, and fields had effective drainage systems. The second group contained abandoned peat extraction fields (Vegetation) covered with poor ground vegetation (dwarf shrubs and herbs), and groundwater was close to the peat surface and without a drainage system. The remaining groups were afforested former peat extraction fields with Scots pine (Pinus sylvestris L.) and silver birch (Betula pendula Roth). Stands corresponded to the Myrtillosa mel. forest type [20] and the groundwater level was lowered (with drainage systems). In all objects, the dominant peat layer was raised bog (Sphagnum), but now the top layer is fen or transitional mire peat.

2.1. Greenhouse Gas Sampling

The ecosystem GHG exchange—CO₂, CH₄, N₂O—measurements were performed using the darkened cylindrical chamber method [21]. The chamber consists of two parts—a chamber made from PVC (H—30 cm, D—50 cm, Volume—65 L) and a collar, which is inserted into the soil (ca. 10 cm deep). In total, five collars (measurement points) were randomly dispersed in a 500 m² circular sample plot (R = 12.62 m) with a 10–15 m distance from groundwater measurement wells. Gas samples were collected using a syringe connected to a tube inserted into the chamber. Gas samples from the chamber were collected in evacuated 100 mL bottles. From each chamber, four GHG samples were taken with 20 min intervals within one hour (at 0; 20; 40; 60 min). Gas sampling was performed once every month in each object. Samples were marked and placed in previously prepared sample boxes. The obtained gas samples were transported to the Climate Change Laboratory of the Department of Geography of the University of Tartu in order to determine the concentrations of greenhouse gases (CO₂, CH₄, N₂O) in ecosystems. Analyses were performed with a Shimadzu GC-2014 gas chromatograph, equipped with an electron capture detector, flame ionization detector and Loftfield autosampler [22].

Environmental parameters were collected during gas exchange measurements in order to obtain factors affecting GHG emissions (Figure 1a,b). Measured environmental parameters were air temperature, soil temperature at 4 depths (5; 10; 20; 30 cm), soil moisture content, groundwater level and samples, water electrical conductivity and water pH. Soil temperature and soil moisture were measured in the same place every month, at one measurement point per object (sample plot). Soil samples were collected in all objects at the initial stage (one time in each object) with a probe. Soil samples were collected from seven depths (0–10; 10–20; 20–30; 30–40; 40–50; 50–100; 100–150 cm) in three replications per object. The volume of each soil sample was 100 cm³. Soil sample preparation and analysis were performed according to the ICP Forests guidelines (

Table 2. Soil characteristics of study objects at all sampling depths (0–150 cm).

Study Object	Ctotal, g kg−1	Ntotal, g kg−1	Ptotal, g kg−1	C/N	pH
Control	570.80 ± 9.70 a	10.45 ± 1.14 a	0.19 ± 0.04 a	64.51 ± 6.05 a	3.28 ± 0.17 a
Vegetation	538.41 ± 7.92 b	9.15 ± 0.76 a	0.19 ± 0.02 a	59.34 ± 3.94 a	2.88 ± 0.07 b
Pinus	537.49 ± 3.77 b	10.07 ± 0.89 a	0.23 ± 0.06 a	60.40 ± 5.10 a	2.68 ± 0.08 c
Betula	528.69 ± 14.25 b	17.31 ± 1.96 b	0.38 ± 0.07 b	40.69 ± 6.67 b	3.50 ± 0.22 a

Note: average ± 95% confidence interval; a,b,c—significant differences.

). The groundwater level was measured each month with a manual ruler as the distance from the ground surface level to water level, in two groundwater wells per object, inserted 2 m deep into the soil. Water samples were collected from two groundwater wells once every two months in each object (in total, 4 times in each object in a season). Sampling, storage, processing and analysis were performed according to the ICP Forests methodology.

To describe the stand characteristics, sample plot measurements were performed in afforested study objects according to the National Forest Inventory methodology [23] for the assessment of tree species, dominant tree height (H), diameter at breast height (DBH) and basal area (BA) measurements (

Table 3. Taxation indices of afforested study objects. DBH—diameter at breast height, H—tree height, BA—basal area.

Study Object	Object	Dominant Tree Species	Age, Years	Trees per ha−1	DBH, cm	H, m	BA, m2 ha−1	Growing Stock, m3 ha−1	C Sequestrated, t C ha−1 year
Pinus	Cepla_7	Pine	42	3920	7.8	9.2	12.8	70.0	0.5
Pinus	Kaigu_7	Pine	72	1500	16.6	24.1	28.3	366.2	1.4
Pinus	Lamb_7	Pine	127	960	16.6	14.2	20.4	135.3	0.3
Betula	Liels_8	Birch	40	1220	13.6	14.7	18.9	155.0	1.2
Betula	Plece_8	Birch	42	1940	14.6	16.2	28.1	257.3	1.8
Betula	Silg_8	Birch	35	1880	10.6	14.7	14.9	113.3	1.0
Vegetation	Cena_3	Birch	NA	700	3.6	4.9	0.7	2.3	NA
Vegetation	Cepla_3	Pine	NA	180	10.6	7.9	1.8	8.8	NA
Vegetation	Silg_3	Birch	NA	660	7.0	6.8	2.0	8.0	NA

Note: NA—information not available.

). No recent human intervention or management in study sites had been carried out for the past 35–40 years. To assess the aboveground carbon (C) stock of Pinus and Betula land use types, the biomass equations developed by Liepiņš et al. [24] and dominant tree species C content coefficients developed by Bārdule et al. [25] were used to calculate the annually sequestrated C content by dividing C stock by stand age.

2.2. Data Analysis

GHG flux was calculated by using the linear regression slope of gas concentration changes in the chamber. If the R² value of the acquired linear regression was less than 0.7, the result was classified as an outlier due to errors during measurements (poorly ventilated chamber, fan decline, etc.). The assessed slope was further expressed as the GHG flux from the area of soil:

$$flux=\frac{M}{R}\frac{P}{T}\frac{V}{t}\frac{slope}{A}$$

(1)

where flux is the soil GHG flux (µg GHG m² h⁻¹); M is the molar mass of GHG (g mol⁻¹); R is a universal gas constant (m³ Pa K⁻¹ mol⁻¹); P is the assumption of air pressure inside the chamber (101,300 Pa); T—air temperature (K); V—chamber volume (0.063 m³); t—time period between first and last GHG flux sampling (0.5 h); slope—slope of the hourly GHG concentration changes inside the chamber; A—collar area (0.1995 m²).

CO₂ emissions were expressed as heterotrophic soil respiration (to exclude ground vegetation impact on emissions) by multiplying the total soil CO₂ flux of vegetation, Pinus and Betula fields with the heterotrophic respiration coefficient (43% of the total respiration) presented by Berglund et al. [26]. Later in the analysis, only the heterotrophic soil CO₂ flux was used.

Data’s normal distribution was checked with the Shapiro test. Statistical significance between average flux values was assessed with the Wilcoxon signed-rank test for non-parametric variables. T-test was used to assess the statistical significance of average values for parametric variables. Correlation was assessed using the Pearson correlation test. All analyses were performed with a 95% confidence interval. Study sites were divided into two categories based on the groundwater level (GWL) due to different emission responses with the change in groundwater level (GWL > 40 cm classified as deep; GWL < 40 cm classified as shallow).

Annual CO₂, CH₄ and N₂O emissions were calculated by summarizing the average monthly emissions in each year, including emission hotspots in the calculation. CO₂eq was calculated by multiplying each GHG with the global warming potential (GWP; CO₂ = 1; CH₄ = 27; N₂O = 273) on a 100-year scale, according to the IPCC [27].

Linear mixed-effect models were developed to predict soil CO₂, CH₄ and N₂O emissions by using the soil temperature (T_soil), groundwater level, groundwater level category (GWL_kat) and land use type (Land_use) as fixed factors, and object (repeated measurement component) and cycle (time component) as random factors. The interaction of continuous and categorical fixed effects was determined using regression models. All data of CH₄ and N₂O emissions were normalized by preserving the original sign (positive or negative) and applying log transformation. Linear models were developed by calculating each land use type’s average CO₂, CH₄ and N₂O emissions from 5 sample chambers in each measurement month to minimize the flux hotspot’s influence on model outliers, and we checked for linear normal distribution and residual vs fitted values. For the CO₂ flux model, only measurements from April to November were included, due to very low emissions in winter months (close to 0 or negative measurements obtained). The best model fit was assessed by comparison of the R² and Aikake information criterion (AIC) values. The significance of fixed effects was estimated by Wald χ² (model ANOVA). All statistical analyses were performed with RStudio [28]. Study data are available on demand from the authors.

3. Results

For all land use types, CO₂-C flux (µg C m² h⁻¹) followed a seasonal trend, where higher emissions could be observed during the vegetation season from April till November, with low (close to 0) emissions during the winter season (Figure 2a). Lower emissions than in all other land use types were observed in control and vegetation fields during the summer season, but, overall, the trend and magnitude was similar between the analyzed land use types. The average heterotrophic soil CO₂-C flux in the control land use type was 18,386.74 ± 2140.93 µg C m² h⁻¹ (average ± 95% confidence interval), in the vegetation land use type 19,598.24 ± 2006.89 µg C m² h⁻¹ and in the Pinus land use type 23,515.34 ± 2062.55 µg C m² h⁻¹, and the highest flux was observed in the Betula land use type, 25,544.41 ± 2062.56 µg C m² h⁻¹. Significant differences between the average CO₂-C flux of two season’s measurements were observed between the control and Pinus, Betula land use types (p < 0.001), and between the vegetation and Pinus, Betula land use types (p < 0.001); however, the average CO₂-C flux between afforested (Pinus and Betula) peat extraction fields was insignificant, and the flux between control and vegetation land use types was also insignificant.

CH₄-C flux (µg C m² h⁻¹) remained low throughout all seasons; however, vegetation fields with CH₄-C flux hotspots were observed during July and August in 2016, and Betula fields with CH₄-C flux hotspots in September 2016 and May 2017 (Figure 2b). The average soil CH₄-C flux in the control land use type was 132.53 ± 148.51 µg C m² h⁻¹, in the vegetation land use type 324.10 ± 139.21µg C m² h⁻¹ and in the Pinus land use type 15.64 ± 143.07 µg C m² h⁻¹, but in the Betula land use type, it was 320.87 ± 143.07 µg C m² h⁻¹. Significant differences (p < 0.001) in the average CH₄-C flux have been observed between control and Pinus, between Pinus and vegetation and between Pinus and Betula. Due to additional emission hotspots, the median could be a better tool to describe the CH₄ emissions in each land use type. The median of CH₄-C in the control plot was −0.312 µg C m² h⁻¹, in vegetation 24.652 µg C m² h⁻¹, in Pinus −12.59 µg C m² h⁻¹ and in Betula −16.01 µg C m² h⁻¹.

N₂O-N flux (µg N m² h⁻¹) followed a similar trend as CH4-C flux—mainly insignificant differences, with exceptional months of emission hotspots (Figure 2c). Higher N₂O-N flux was observed for Betula fields compared to other land use types; however, additional months of high N₂O-N flux were observed also for the control land use type. The average soil N₂O-N flux in the control land use type was 5.29 ± 2.03 µg N m² h⁻¹, in the vegetation land use type 0.51 ± 1.91µg N m² h⁻¹, in the Pinus land use type −0.56 ± 1.96µg N m² h⁻¹ and in the Betula land use type 7.58 ± 1.96µg N m² h⁻¹. Statistically significant differences (p < 0.001) in the average N₂O-N flux have been observed between all land use types, but not between the vegetation and Pinus land use types, and between the Betula and control land use types.

The annual soil heterotrophic CO₂-C flux (t ha⁻¹ yr⁻¹) was the lowest in the control, followed by the vegetation, Betula and Pinus land use types (

Table 4. The annual soil CO₂-C, CH₄-C, N₂O-N flux and CO₂eq in study objects.

Study Object	CO2-C, t ha−1 yr−1	CH4-C, kg ha−1 yr−1	N2O-N, kg ha−1 yr−1	CO2eq, t ha−1 yr−1
Control	1.46 ± 0.05 a	10.39 ± 0.59 a	0.44 ± 0.04 a	1.86 ± 0.08
Vegetation	1.62 ± 0.06 b	27.07 ± 1.92 b	0.04 ± 0.01 b	2.37 ± 0.11
Pinus	2.30 ± 0.07 c	−1.41 ± 0.06 c	−0.04 ± 0.01 c	2.26 ± 0.07
Betula	2.23 ± 0.06 c	25.11 ± 2.25 b	0.61 ± 0.03 d	3.08 ± 0.13

Note: average ± 95% confidence interval; a,b,c,d—significant difference indicator.

). Differences were statistically significant between the control and all analyzed land uses, and vegetation and all land use types. The annual CH₄-C flux (kg ha⁻¹ yr⁻¹) was the highest in the vegetation and Betula land use types, but the lowest (negative) in the Pinus land use type. Differences between the annual CH₄-C flux were significant between all analyzed land use types, except for vegetation and Betula. The annual N₂O-N flux (kg ha⁻¹ yr⁻¹) was the highest in the Betula land use type, followed by the control and vegetation and the lowest in the Pinus land use type. Furthermore, N₂O-N flux was close to zero in Pinus and vegetation fields. Significant differences were observed between all analyzed land uses. Overall, the annual CO₂eq (t ha⁻¹ yr⁻¹) was the lowest in the control, followed by the Pinus and vegetation and the highest in the Betula land use type.

3.1. CO₂ Emission Model

A linear mixed-effect model for heterotrophic soil CO₂-C respiration was developed using the soil temperature (T_soil), groundwater level (GWL) and land use type (Land_use) as fixed factors. The model’s conditional R² was high—0.76 (Figure 3,

Table 5. The linear mixed-effect model of CO₂-C flux variable significance (model coefficients presented in Table S1).

Variable	χ2	p-Value
(Intercept)	5.6188	<0.05
T_soil	132.3021	<0.001
Land_use	21.6026	<0.001
GWL	7.8972	<0.01

). The estimated CO₂-C flux and the actual CO₂-C flux fit was high (R² = 0.80). All variables included in the model were statistically significant (p < 0.001).

3.2. CH₄ Emission Model

A linear mixed-effect model for soil CH₄-C respiration was developed using the groundwater level (GWL), groundwater level category (GWL_kat), land use type (Land_use) and interaction between all mentioned variables as fixed factors. The model’s conditional R² was high—0.70 (Figure 4). The model’s estimated CH₄-C flux and the actual CH₄-C flux fit was high (R² = 0.74). Significant variables in the model were land use type and its interaction with GWL_kat, and the interaction between GWL, GWL_kat and land use type (

Table 6. The linear mixed-effect model CH₄-C flux variable significance (model coefficients presented in Table S2).

Variable	χ2	p-Value
(Intercept)	1.33	0.25
Land_use	13.65	<0.01
GWL	0.53	0.47
GWL_kat	0.44	0.51
Land_use:GWL	4.80	0.19
Land_use:GWL_kat	10.21	<0.05
GWL:GWL_kat	0.007	0.93
Land_use:GWL:GWL_kat	10.52	<0.05

).

3.3. N₂O Emissions

A linear mixed-effect model for N₂O-N flux could be developed; however, none of the environmental parameters measured and included in the model showed a significant response; thus, an equation was not developed. Therefore, an assumption that N₂O-N emissions are affected by different factors besides CH₄-C and CO₂-C was made.

4. Discussion

The after-use management practices of former peat extraction fields may be very diverse in terms of climate change mitigation targets: firstly with carbon sequestration capacity, secondly with impact on climate and nature (soil emission magnitude) and thirdly with establishment costs and economic viability. The average soil temperature in the analyzed study sites differed significantly between afforested sites (Pinus and Betula) and sites without tree cover (control and vegetation) (Figure 1b), due to tree canopy cover that reduced the solar radiation’s influence on the soil temperature [29]. Furthermore, groundwater level differences between sites were notable (Figure 1a). The highest groundwater level was observed for the vegetation field, where no ditch (drainage) system had been established, but the groundwater level on all other fields fluctuated throughout the measurement campaign. The peat chemical content also significantly differed between the land use types (Table 2).

Tree planting in former peat extraction fields appears as one of the best options for re-cultivation; however, much is dependent on the selected tree species due to differences in biomass growth and carbon sequestration [30,31]. Carbon sequestration in tree aboveground biomass ensures additional C annual sequestration that ranges from 0.3 to 1.4 t C ha⁻¹ year (average 0.7 t C ha⁻¹ year) for Pinus fields, and from 1.0 to 1.8 t C ha⁻¹ year (average 1.3 t C ha⁻¹ year) in Betula fields (Table 3). Higher average C sequestration can be observed for Betula fields; however, the stand age in Pinus fields was notably higher. Soil heterotrophic CO₂-C flux followed a similar trend throughout the measurement campaign, and the analyzed field average and annual CO₂-C emissions can be ranked in increasing order: Control < Vegetation < Pinus < Betula (Figure 2a, Table 4). The lowest CO₂-C emissions in control sites can be explained by the lack of ground vegetation and tree cover; thus, there was no carbon input in peat from litter, fine roots and microbial processes [32]. Other land use types have several other biological processes influencing soil heterotrophic respiration, such as carbon loss from the decomposition of litter detritus and soil organic matter by microorganisms [33]. Comparing the average CO₂-C emissions from our analyzed fields with other studies of peat extraction fields, the emission range is similar, ranging from 25,000 to 29,000 µg C m² h⁻¹ [34], compared to the average values of 18,000–25,500 in our study (Figure 2a). However, a comparison of heterotrophic CO₂-C emissions from forests on drained peatlands in Finland to our obtained results for afforested sites showed lower CO₂-C flux values than in Finland [35]. Soil CH₄-C flux throughout the measurement seasons remained low (close to 0), with exceptional hotspots for different fields and months. Overall, the lowest (negative) average and annual CH₄-C flux was observed for the Pinus field, followed by low emissions for the control field and Betula field and the highest for the vegetation field (Figure 2b, Table 4). The highest CH₄-C flux for the vegetation field could be explained by the lack of drainage systems and high groundwater level for these fields, which promote soil CH₄-C emissions [36,37]. The relatively high annual CH₄-C flux for Betula sites can be explained by the high emissions in different hotspots and measurement months, but the median of CH₄-C flux for the Betula land use type was negative (−16.01 µg C m² h⁻¹), indicating a similar amount of methane accumulation in the soil compared to the Pinus land use type. The highest average, median and annual N₂O-N flux were observed for Betula fields, followed by the control; however, in vegetation fields, the emissions were low (close to 0), and in Pinus fields, the average N₂O-N emissions were even negative (Figure 2c, Table 4). In the comparison of GHG CO₂eq values, the lowest emissions were observed for control fields, followed by Pinus, vegetation and Betula (Table 4). The results of annual CO₂eq indicate Pinus fields as one of the most suitable re-cultivation practices, due to low CO₂eq emissions, negative annual CH₄ and N₂O flux and with additional carbon sequestration in living tree biomass during tree growth (Table 3). Moreover, afforestation with Betula could be a feasible option due to carbon sequestration in the biomass; however, the CH₄-C emission hotspots and high N₂O-N flux increase the annual CO₂eq values.

The linear mixed-effect models indicate differences in the analyzed factors influencing GHG emissions. The main CO₂-C influencing factors based on the linear mixed-effect model are the soil temperature, groundwater level and land use type (Table 5). The significance of the soil temperature on CO₂-C emissions has been reported in various studies [38,39]; however, the groundwater level and land use types reflect the differences in the analyzed re-cultivation practices. Overall, the developed CO₂-C flux model had a good linear fit (R² = 0.80) to the actual CO₂-C measurement data (Figure 3, Table S1). The model of CH₄-C flux determined that the main influencing factors are the groundwater level, land use type and the interaction between groundwater level, land use type and groundwater level category (Table 6). Several studies have reported the importance of the groundwater level for GHG emissions. A high water table increases methanogenesis and favors methanotrophy; thus, a high groundwater level increases the soil CH₄-C emissions [35,40,41], as observed in our study, where the highest CH₄-C fluxes were observed in the land use type without drainage systems (vegetation land use type). The estimated model of CH₄-C flux to the actual CH₄-C emissions showed a good fit (R² = 0.73) to the developed model (Figure 4, Table S2). The relatively high model fit allows us to use the developed models for the GHG flux prediction of peat extraction fields and the emissions under the field’s further re-cultivation with ground vegetation or forests. However, N₂O-N emissions could not be estimated with the currently measured environmental factors; therefore, a model was not developed. N₂O-N emission modelling was unsuccessful in this study, possibly due to differences in the principal processes causing N₂O-N emissions from the soil (alternative electron acceptors inhibiting denitrification) between the analyzed land use types [42].

5. Conclusions

The most suitable secondary use of former peat extraction fields is afforestation, especially with Scots pine, in order to favor climate change mitigation goals, as a positive impact can be reached in terms of low annual CO₂eq and additional carbon sequestration in tree aboveground biomass. The main factors influencing CO₂ emissions are the soil temperature, groundwater level and land use type. However, the main CH₄-emission-influencing factors are related to the groundwater level, land use type and the interaction between these variables. Based on the obtained influencing factors, the developed linear mixed-effect models showed a good model fit for the analyzed land use types, and thus can be used for CO₂ and CH₄ flux estimation.