Ecosystem respiration during snowmelt and soil thaw leads to a rare annual CO₂ net loss in a boreal fen

Although boreal peatlands play a critical role in the global carbon cycle, their year-round carbon dioxide (CO₂) dynamics, and particularly the contribution of the non-growing season, remain poorly constrained in annual balance estimates. Using 17 years (2005–2021) of eddy covariance measurements from a fen in southern Finland, we first quantified the magnitude, timing, and interannual variability of CO₂ fluxes. We then examined in greater detail the non-growing season, specifically the non-productive season defined by the net productivity of the system. We assessed the flux drivers during different periods of the year, with particular emphasis on soil temperature dynamics and the role of thermal legacy effects. On average, the non-productive season accounted for 60 % of the year (226 ± 27 d), ranging from mid-September to late April, and offset 57 % (±33 %) of the subsequent productive season's CO₂ uptake. Emissions declined from autumn to spring, with the highest carbon emissions occurring across September–December and the lowest in January–February. Soil temperature, both concurrent and lagged up to four months, was the main control of CO₂ fluxes during November–December and spring thaw, while photosynthetically active radiation (PAR) dominated during the onset of the non-productive season. Variability in annual CO₂ balances was large, and in two years (2016 and 2018) the fen switched from a net CO₂ sink to a source. Finally, we focused on 2016 in detail: an exceptional six-week CO₂ release during April–May released 84 g C m⁻², offsetting 38 % of the following productive season's CO₂ uptake. This event was linked to unusually warm late-autumn soils, minimal snow insulation, and subsequent rapid surface freezing, which likely enhanced CO₂ accumulation and stimulated CO₂ release during thaw. Our results demonstrate that short-lived but intense events during the non-productive season can determine the annual peatland CO₂ balance and therefore significantly affect the annual carbon budget of boreal peatlands.

1 Introduction

Boreal and sub-Arctic peatlands store between 400 Gt (Gorham, 1991; Yu et al., 2011) and 1000 Gt of carbon (Nichols and Peteet, 2019), contributing up to 30 % of the global soil carbon pool (Friedlingstein et al., 2022). These ecosystems have had a net cooling effect on the atmosphere over the past millennia, as carbon uptake by plants has outweighed losses from heterotrophic respiration (Hugelius et al., 2020). However, global warming – especially at high latitudes, where temperatures are rising three to four times faster than the global average (AMAP, 2021; Rantanen et al., 2022) – may alter the balance between carbon sequestration and loss in boreal and sub-Arctic regions. Elevated temperatures have already intensified carbon cycling by enhancing both productivity and respiration in these ecosystems (See et al., 2024).

Carbon dynamics outside the growing season, i.e., late autumn, winter, and early spring, play a critical yet understudied role in shaping the annual carbon balance of boreal peatlands. These non-growing season periods have become increasingly important as warming trends intensify during the coldest months (Rantanen et al., 2023). While fluxes are typically low in mid-winter, cumulative emissions can offset a significant proportion of growing season uptake, with up to three-quarters of summer sequestration negated in some years based on in situ measurements in the boreal–Arctic zone (Arndt et al., 2023; Virkkala et al., 2022). These emissions are expected to increase even further in the future (Natali et al., 2019), linked with decreasing snow depth and snow cover duration (Pongracz et al., 2024).

Carry-over effects from previous seasons can strongly influence CO₂ fluxes during the non-growing season. For example, high microbial activity during warm summer and autumn can lead to gas build-up in the soil, which is then subsequently released during soil freeze-up or thaw (Raz-Yaseef et al., 2017; Sullivan et al., 2012). Pulses of high carbon release during freezing and thawing have been observed for both methane (Kübert et al., 2026; Mastepanov et al., 2013) and CO₂ (Raz-Yaseef et al., 2017; Arndt et al., 2020; Wang et al., 2023). These events, though infrequent, have been shown to have the potential to reduce annual carbon uptake by up to 46 % (Raz-Yaseef et al., 2017), highlighting the importance of understanding carbon exchange across the full annual cycle.

Despite their non-trivial implications for net carbon sink strength, non-growing season processes remain underrepresented in carbon flux studies, largely due to logistical and environmental challenges. Data from the ABCFlux database show that the majority of flux measurements still occur during the growing season, with June to August accounting for 32 % of observations, compared to just 18 % for winter (Virkkala et al., 2022; Virkkala et al., 2025; See et al., 2024). Eddy covariance (EC) sites provide high-temporal-resolution measurements of atmosphere–biosphere exchanges at the ecosystem scale and are the only method capable of directly quantifying net carbon fluxes across these areas. However, relatively few EC sites operate continuously year-round (Pallandt et al., 2022), leaving major gaps in our understanding of winter emissions and seasonal transitions (Hugelius et al., 2024).

This study utilizes a unique 17-year (2005–2021) time series dataset (Alekseychik et al., 2024) of the net ecosystem exchange of CO₂ (NEE) from the Siikaneva fen to understand the overall interannual variability and examine the role of non-growing season periods in the annual carbon balance of a boreal fen. Given the scarcity of long-term, year-round measurements, this study offers a rare dataset for assessing long-term carbon exchange in these ecosystems.

The aim of this study is to assess the non-growing season carbon exchange and identify periods within the non-growing season that contributed significantly to the annual carbon balance. We specifically ask:

1. How much does the non-productive season contribute to the annual CO₂ balance in the Siikaneva fen, and what are the seasonal patterns of CO₂ fluxes during the non-growing season?

2. What are the primary environmental drivers of CO₂ fluxes during the non-growing season? Specifically, do prior soil thermal conditions predict non-growing season fluxes, particularly during high-emission periods?

We hypothesize that the non-productive season contributes substantially to the annual CO₂ balance, with distinct seasonal patterns characterized by higher emissions in autumn and spring thaw periods. We expect soil temperature to be the primary driver of non-growing season CO₂ fluxes, with soil thermal conditions from prior seasons potentially affecting CO₂ exchange, especially during late autumn and spring soil thaw.

2 Methods

2.1 Site description

Siikaneva wetland (Fig. 1) study site is an oligotrophic fen within a large aapamire complex located in southern Finland (61°50^′ N, 24°12^′ E, 162 m above sea level), with ongoing measurements since 2005. The site has been part of the ICOS research infrastructure since 2017 and is classified as a Class 2 ecosystem station. Peat depth at the site ranges from 2 to 6 meters. The terrain is generally flat, with vegetation dominated by peat mosses (Sphagnum balticum [Russow] C.E.O. Jensen, S. majus [Russow] C.E.O. Jensen, S. papillosum Lindb.), sedges (Carex rostrata Stokes, C. limosa L., Eriophorum vaginatum L.), and Rannoch rush (Scheuchzeria palustris L.) (Korrensalo et al., 2022). On one site the fen is neighboured by Scots pine forest growth on mineral soil. The mean annual temperature for the period 1971–2000 was 3.3 °C, and the total annual precipitation was 713 mm, based on data from the Hyytiälä weather station, located 5 km from Siikaneva (Drebs et al., 2002).

Figure 1 Photo composite featuring Siikaneva in September (left), December (middle) and July (right). Photos by Maximilian King and Karoliina Särkelä.

2.2 Non-growing season definitions: non-productive season and thermal non-growing season

Many different definitions exist for the growing season and non-growing season, which can be based on environmental variables such as air temperature (referred to as the thermal growing season when persistently above 5 °C (Kollo et al., 2023; Ruosteenoja et al., 2016)), snow cover period (Arndt et al., 2020; Rafat et al., 2022), or thresholds in gross primary production (Böttcher et al., 2014). In this study, we define growing season and non-growing season based on the ecosystem carbon exchange (NEE), following one of the definitions suggested by Körner et al. (2023), and specifically referred to as productive or non-productive season. The productive season corresponds to the period when the ecosystem acts as a net carbon sink (daily NEE < 0), and the non-productive season when it is a net carbon source (daily NEE > 0). We used a threshold of three days (following the sink period definition by Aurela et al., 2004): the productive season starts on the first of three consecutive days with negative NEE and ends on the first of three consecutive days with positive NEE.

Additionally, we evaluated thermal growing season and non-growing season We defined the thermal growing season for each year as the period between the day of year when the cumulative sum of daily mean air temperature above a 5 °C threshold (air temperature °C −5 °C) reached its minimum, corresponding to the point before sustained daily temperatures consistently exceeded 5 °C, and the day of year when this cumulative sum reached its maximum, corresponding to the point after which daily temperatures no longer contributed to net accumulation above the 5 °C threshold (Ruosteenoja et al., 2016).

2.3 Eddy covariance (EC) measurements

The net ecosystem exchange (NEE) of CO₂ at the Siikaneva fen site was measured using the eddy covariance (EC) technique (Baldocchi et al., 2003), providing continuous, high-frequency estimates of ecosystem–atmosphere CO₂ exchange. Measurements have been conducted since 2005, with the flux tower installed at a height of 2.7 m above the peat surface. 80 %–90 % footprint was used for flux calculations, which varies between tens of meters to a maximum of 200 m depending on atmospheric conditions. Between 2005 and 2015, the EC system consisted of a Metek USA-1 three-dimensional sonic anemometer (Metek GmbH, Elmshorn, Germany) and a LI-COR LI-7000 CO₂ $/$ H₂O gas analyzer (LI-COR, Lincoln, NE, USA). In late 2015, the system was upgraded to a Gill HS-50 anemometer (Gill Instruments Ltd., UK) and a LI-COR LI-7200 enclosed-path gas analyzer, resulting in a data gap between September 2015 and February 2016. The processed and gapfilled data were available at the Finnish metadata catalog Etsin. Raw 10 Hz data were processed in EddyUH (Mammarella et al., 2016), including despiking, double coordinate rotation, sonic temperature correction, frequency response correction, and were calculated using mixing ratios. Fluxes were filtered for turbulence intermittency and atmospheric stability and corrected for CO₂ storage below the measurement height. Further details on EC data processing can be found in Mammarella et al. (2016). Missing NEE data was gapfilled with the sum of modelled gross primary productivity (GPP) and ecosystem respiration (Reco). GPP was modelled as a function of photosynthetically active radiation (PAR) and air temperature, with parameters optimized for the preceding two-week period Photosynthesis was constrained to zero when air temperature dropped below 0 °C. Reco was modelled as a function of air temperature using site- and period-specific coefficients. The equations for the flux partitioning are provided in Kulmala et al. (2019) with site-specific modifications – for the flux timeseries in Siikaneva fen, air temperature was used to model the Reco as soil temperature records were inconsistent across the study period. Processed and gap-filled hourly NEE data were aggregated to daily sums and converted from µmol CO m⁻² s⁻¹ to g C m⁻² d⁻¹. Negative NEE values indicate net ecosystem carbon uptake.

2.4 Ancillary measurement

We reviewed a range of ancillary environmental variables alongside CO₂ fluxes, including soil temperature, photosynthetically active radiation (PAR), vapor pressure deficit (VPD), precipitation, water table depth, and snow depth.

Soil temperature was measured using Campbell 107 thermistors (Campbell Scientific, Logan, UT, USA) from 2005 to 2016, and UMS TH3-s temperature profile probes (UMS GmbH & Co. KG, Willmars, Germany) from 2017 to 2021. The latter represents the average of five microsites within the EC footprint. Data from 5 and 50 cm depths (2005–2016) and from 5 and 45 cm depths (2017–2021) were used, specifically from a lawn microsite representative of the EC footprint. For consistency, temperatures at 45 and 50 cm are both referred to as “soil temperature at 50 cm depth”. Photosynthetically active radiation (PAR, µmol m⁻²) was measured with a Li-Cor Li-190R quantum sensor (LI-COR, Lincoln, NE, USA), but data from 2009–2015 were excluded due to sensor malfunction. Data from the Hyytiälä forestry site (5 km to the northeast of Siikaneva) were used for this period, measured with a Li-Cor Li-190SZ quantum sensor (LI-COR, Lincoln, NE, USA). Vapor pressure deficit (VPD, kPa) was calculated using the fCalcVPDfromRHandTair function in the R package ReddyProc (Wutzler et al., 2018), based on air temperature and relative humidity from a Rotronic HC2 sensor (Rotronic AG, Bassersdorf, Switzerland). Water table depth (WTD) was measured with Druck PDCR1830 and Campbell CS451 pressure transducers (Campbell Scientific). WTD is expressed in centimeters (cm) relative to the surface, with negative values indicating depths below ground. Precipitation (mm) was recorded using an ARG-100 tipping bucket rain gauge (2005–2016), which underreports sleet and snow. In 2017, it was replaced by an OTT Pluvio2S weighing rain gauge (OTT HydroMet, Kempten, Germany). Instrument upgrades in late 2015 caused data gaps across several variables (dates vary by parameter). To assess autumn and winter 2015 conditions, air temperature and snow depth data were retrieved from the Hyytiälä Forestry Station, 5 km from Siikaneva. Air temperature there was measured with a Pt100 sensor inside a custom shield. A linear regression comparing air temperatures at the two sites yielded a slope of 0.98 and R² of 0.98, validating the use of Hyytiälä's temperature data.

2.5 Freeze–thaw status classification

The freeze–thaw state of the study area was classified using the Sentinel-1 IW GRD SAR (Synthetic Aperture Radar) data available from autumn 2014 onwards. These data have been processed by Google Earth Engine into calibrated and ortho-rectified grids of Sigma Nought backscatter values with 10 m pixel size. Median backscatter time series (VV, VH, and VH $/$ VV ratio) were extracted for a 200 m zone around the Eddy Covariance mast and the broader Siikaneva Fen. The 200 m zone represents ∼ 75 % flux footprint area (Riutta et al., 2007). Data from both ascending and descending orbits were used. The temporal resolution was, on average, two acquisitions per week. Due to changes in liquid water content, soil freezing typically reduces the measured backscatter intensity by a few decibels (Hallikainen et al., 1985; Ulaby et al., 1982; Cohen et al., 2021), and the presence of wet snow further reduces the backscatter by an additional few decibels (Nagler and Rott, 2000; Luojus et al., 2007)

Freeze–thaw periods were identified through visual inspection of backscatter trends, supported by in situ measurements (air and soil temperature, soil water content, and snow depth). Freezing onset was marked by a ≥2 dB drop in backscatter in VV polarization; the frozen period continued as long as all orbits consistently showed low values. Thaw onset was defined either by a sharp drop indicating wet snow or by a rise in backscatter consistent with thawed soil; thaw continued as long as all orbits showed higher backscatter values typical of thawed conditions.

2.6 Data analysis

We quantified the relative importance of environmental drivers controlling the carbon sink or source strength, expressed as daily net ecosystem exchange (NEE), using a machine-learning random forest regression algorithm. Random forest is a machine-learning technique widely used for complex multi-regression analyses (e.g., López-Blanco et al., 2017; Wei et al., 2022) that builds an ensemble of decision trees by repeatedly sampling random subsets of the training data (Breiman, 2001; Pedregosa et al., 2011). Each decision tree partitions the explanatory variables (covariates) to create clusters of data, within which regression models predict the response variable (NEE). Variable importance is then quantified as the percentage of times a variable is used in the decision splits across all trees, reflecting its contribution to explaining variability in NEE.

Our set of explanatory variables included photosynthetically active radiation (PAR), vapor pressure deficit (VPD), water table depth (WTD), and soil temperature. All variables and NEE were aggregated to daily means. Daily data were then binned into two-month periods: January–February, March–April, May–June, July–August, September–October, and November–December. These periods loosely follow the approximate periods when the system turns into a net source, i.e., the onset of the NGS (September–October); when the surface freezes (November–December); when the surface is typically frozen (January–February); when the surface thaws (March–April); when the system turns into a net sink of CO₂, i.e., the onset of the GS (May–June); and when it is consistently a net sink of CO₂ (July–August). Variable importance was then assessed for each model, with standard deviations calculated to evaluate uncertainty.

Days with fewer than 12 hourly NEE observations were excluded from the analysis to ensure data quality. Due to incomplete PAR records at the Siikaneva site for 2009–2015, PAR data from the nearby Hyytiälä Forest Station (5 km to the northeast of the Siikaneva fen) were used as a proxy. Due to strong collinearity between air temperature and soil temperature, we did not include air temperature in the analysis. Precipitation was excluded as an explanatory variable, since the instrument measuring precipitation could not reliably detect snow and sleet. Therefore, such precipitation data could introduce biases in the analysis by underestimating the amount of precipitation during winter compared to summer.

To investigate the potential influence of soil thermal conditions over different temporal scales and legacy effects, three separate model frameworks were constructed, each differing in how soil temperature was aggregated:

• a.

  Daily mean soil temperature at 5 cm depth, representing contemporaneous conditions.

• b.

  Daily running mean of 60 d of soil temperature at 5 cm, representing the conditions in the last two months.

• c.

  Daily running mean of 60 d of soil temperature at 5 cm, lagged by 60 d, representing the conditions 2–4 months prior. For example, the soil temperature data point used as a driver for the daily NEE on 1 April is the mean soil temperature between 2 December and 30 January.

3 Results

3.1 Duration of the non-productive and thermal non-growing season

On average, the non-productive season accounted for 60 % (226 ± 27 d) of the year, beginning on 17 September (±19 d) and ending on 30 April (±15 d). The non-productive season started latest in 2011 on 26 October and earliest in 2018 on 30 July (Fig. 2a) and ended the earliest in 2006 on 18 March and latest in 2016 on 2 June (Fig. 2a). The non-productive season that started in 2018 and ended in 2019 was the longest NGS (278 d) in the whole study period, due to an exceptionally early start of the season.

On average, the thermal non-growing season accounted for approximately 53 % of the year (195 ± 17 d), beginning on 13 October (±10) and ending on 26 April (±9), therefore being 31 d longer than the non-productive season. The mean absolute difference between the non-growing season periods based on the two definitions was 27 d for the onset and 12 d for the end. The largest delay between the onset of the thermal growing season and the start of the productive season occurred in spring 2016 (thermal growing season onset 27 April vs. productive season onset 3 June), whereas the earliest onset of the non-productive season relative to the thermal non-growing season occurred in 2018 (30 July non-productive season onset vs. 22 October thermal non-growing season onset).

Figure 2 (A) Daily net ecosystem exchange of CO₂ (NEE; g C m⁻² d⁻¹) in 2005–2021 and the non-productive season (black rectangular outline) periods centered in the winter. Each row represents one year, beginning in July and ending in June of the following year. Colors indicate daily mean NEE values, where the yellow–green gradient reflects negative values (net carbon sink) and the yellow–pink gradient reflects positive values (net carbon release). Black rectangular outlines mark the NGS periods, defined as the period beginning after the first three consecutive days of positive NEE and ending after the first three consecutive days of negative NEE. On average, the NGS started on 17 September and ended on 30 April. (B) Daily mean air temperature for 2005–2021, with the thermal non-growing season centered on winter. Each row represents one year, beginning in July and ending in June of the following year. Black rectangular outlines mark the thermal non-growing season following the definition of Ruosteenoja et al. (2016). The season boundaries are determined from the annual cumulative sum of daily temperature deviations from a 5 °C threshold: the growing season begins when the cumulative sum reaches its minimum and ends when it reaches its maximum, indicating that temperatures remain persistently above or below 5 °C. On average, the thermal non-growing started on 13 October and ending on 26 April.

3.2 Interannual and interseasonal CO₂ fluxes

The fen acted as a net CO₂ sink during the study period, with an average annual uptake of −51 ± 39 g C m⁻². In two years, 2016 and 2018, the fen was a net source of CO₂, releasing 28 and 21 g C m⁻², respectively (Fig. A1a). However, when looking at the pairs of full non-productive and productive seasons (Fig. 3a), only productive season in 2016 and the preceding non-productive season from 2015–2016 were a net source of CO₂. The mean cumulative carbon release during the non-productive seasons was 49 ± 15 g C m⁻² (Fig. 3a), offsetting on average 57 % of the following productive season carbon uptake (−100 ± 31 g C m⁻²) (Fig. 3a). Interannual variability was high (±33 %). The lowest cumulative carbon emissions over the non-productive season were observed in 2012 (26 g C m⁻²), while the highest were recorded in 2016 (84 g C m⁻²) (Fig. 3a). The mean daily CO₂ flux during non-productive season ranged from 0.14 g C m⁻² in 2012 to 0.33 g C m⁻² in 2016. Mean during the whole study period was 0.21 ± 0.048 g C m⁻².

During the non-productive season, CO₂ emissions generally declined from autumn to spring (Fig. 3a and b). The transitional months of October and April, marking the shift into and out of the growing season, showed the highest interannual variability (Fig. 3a and b). Emissions peaked in November (9.03 ± 2.07 g C m⁻²), followed by October (8.09 ± 3.60 g C m⁻²) and December (7.9 ± 2 g C m⁻²; Fig. 3b). Mean CO₂ emissions during October–December were 75 % higher than those in January–March (8.4 vs. 4.8 g C m⁻²; Fig. 3b). April and May of 2016, both part of the non-productive season that year, exhibited the highest monthly net CO₂ losses and the largest deviations from the long-term mean, followed by August and September in 2018 (Figs. 3a and 2b). Fluxes during that 2016 period exceeded the 95th percentile for 32 nearly consecutive days (see Fig. S5 in the Supplement). All other periods with three or more consecutive days above the 95th percentile occurred between September and December, with event durations of up to five days (see Fig. S5).

The strongest CO₂ sinks occurred in June and July (−33 ± 9.6 and −33 ± 11 g C m⁻², respectively; Fig. 3B). Monthly variability was higher during periods of net CO₂ sink (May–September; Fig. 3b) than during the period of CO₂ source (October–April; Fig. 3b).

There was a significant data gap in the EC measurements from 22 September 2015, to 25 February 2016. Notably, 2016 had the highest annual carbon emissions (1 January– 31 December, Fig. A1a). To ensure that this finding was not an artifact of the extended gap-filled period, annual balances were recalculated for all years, excluding data from 1 January to 25 February (Fig. A1b). Since 2016 still showed the highest annual balance (meaning a net source of CO₂), we conclude that the gap-filling procedure did not significantly bias the results. Additionally, we evaluated the gapfilling performance after the data gap and during the anomalously high flux values. The adjusted R² between the gap-filled NEE and the highest quality measured data from spring 2016 (26 February–30 May) was 0.78, indicating a generally good agreement (see Fig. S7).

Figure 3 (a) Non-productive season and productive season net ecosystem exchange (NEE) sums of CO₂ for each full non-productive season and the following productive season, with the net balance of the combined period indicated by black points. Negative values indicate a carbon uptake, and positive values indicate carbon loss to the atmosphere. We use the terms productive season and non-productive season to refer to periods when the ecosystem acts as a net carbon sink (productive) or a net carbon source (non-productive), based on daily NEE values, following the definition for a sink period by Aurela et al. (2002). The start of the non-productive season is defined as the first of three consecutive days with positive NEE, and the productive season begins on the first of three consecutive days with negative NEE. The non-productive season sum of CO₂ is calculated for the full season, starting in autumn and ending in spring. (b) Monthly NEE sums across the study years, with the mean monthly sums (black line) and standard deviation (shaded area) shown for the period 2005–2021.

3.3 Respiration release in spring 2016

Since 2016 exhibited the highest cumulative CO₂ emissions during the non-growing season (84 g C m⁻²) and was one of the few years with an annual net carbon loss (28 g C m⁻²), coupled with unusually high net emissions of CO₂ in April and May (Fig. 3b), we examined this year in greater detail to better understand the underlying drivers and the conditions leading up to and during the peak emission period (Fig. 4).

The late autumn of 2015 was notably warm, with November-December reaching the highest mean temperature in the 2005–2021 record at 1.69 °C, 3.22 °C above the long-term average (Fig. 4c). A sharp temperature drop in January 2016 drove surface soil temperatures to their lowest levels recorded from 2005 to 2021, with air temperatures plunging from above zero to −25 °C over a two-week period (Fig. 4c). The cold spell lasted nearly three weeks, and the second-lowest air temperature in the study period (−25.2 °C) was recorded. Snow depth was not measured for most of the winter in early 2016 at Siikaneva. Measurements at Hyytiälä (5 km to the northeast of Siikaneva) recorded no snowpack at the beginning of the cold spell. During the cold spell, snow depth only increased from 0 to 4 cm – well below the 2016–2021 average of 16 cm. Surface soils froze, but the deeper soil layers remained warmer throughout the winter (Fig. 4b) due to the warm late autumn. From January to March 2016, the mean soil temperature at 50 cm depth was about 0.85 °C higher than the long-term at 2.95 °C. The soil temperature showed unique fluctuations at 50 cm depth but not in shallower layers (Fig. S2).

Complementary satellite-derived freeze–thaw status classification indicated the beginning of thawing for the 200 m radius from the EC mast on 27 March, marking the start of snowmelt (see Table S2 in the Supplement). At this point, CO₂ flux increased sharply while soil temperatures at 5 cm remained at 0 °C (Fig. 4a–b), together with CH₄ flux (see Fig. S8). The CH₄ flux decreased after around two weeks, while elevated CO₂ emissions persisted throughout the entire thawing period (see Fig. S8). Once the soil temperature at 5 cm began to rise more rapidly in the beginning of May, a second, more pronounced peak in CO₂ emissions followed, fluctuating in response to the warming topsoil (Fig. 4a–b). In April 2016, daily emissions (positive NEE) were nearly four times higher than the monthly average: 0.53 g C m⁻² d⁻¹ compared to the typical 0.13 g C m⁻² d⁻¹. Cumulative emissions in April and May were substantial, offsetting 38 % of the total carbon uptake observed during the subsequent productive season.

To assess whether the high positive values of NEE were related to a delayed onset of photosynthetic activity, we evaluated GPP (see Fig. S9). Overall, GPP in spring 2016 closely followed the long-term mean and did not exhibit anomalous behaviour. The onset of photosynthetic activity, defined as the first occurrence of five consecutive days with GPP exceeding 1 g C m⁻² d⁻¹, occurred on 6 May, which is five days earlier than the 17-year average and within the range of interannual variability. This timing coincided with the decline in NEE and the increase in CH₄ emissions (see Fig. S8), suggesting the onset of biological activity. However, the NEE stayed positive indicating a net release of CO₂ until June 3rd, when it turned permanently negative (see Fig. 2b).

Figure 4 (a) Daily means of Net Ecosystem Exchange of CO₂, (b) soil temperature at 5 and 50 cm depth, and (c) air temperature, snow depth, and water table depth (expressed as centimeters from the surface, where positive values indicate water level above the surface) from 1 November 2015, to 1 June 2016. In each panel, the solid line shows observed values, the dotted line represents the mean for each day of the year (2005–2021, excluding 2015–2016), and shading indicates the standard deviation. Air temperature was measured 9 km away at the Hyytiälä Forest Station due to a data gap at the Siikaneva site in autumn 2015. The blue shading highlights the period of high flux values observed in spring 2016.

3.4 Variable importance controlling NEE

We used machine learning random forest models to evaluate the relative importance of environmental drivers controlling CO₂ exchange, with a particular focus on the role of soil temperature over different time scales. The pronounced respiration pulse observed in spring 2016 (Fig. 4a) suggested that soil thermal conditions from the preceding autumn and winter (Fig. 4b) might influence carbon emissions during soil thaw. To investigate this, we compared three temporal representations of soil temperature as predictors: concurrent daily soil temperature (Fig. 5a), a 60 d running mean (Fig. 5b), and a 60 d running mean lagged by two months (Fig. 5c). The 60 d window was chosen to capture potential cross-seasonal effects of soil temperature, such as the influence of summer conditions on autumn respiration and late-autumn temperatures on winter CO₂ efflux.

Random forest models that were used to evaluate the relative importance of environmental drivers controlling CO₂ exchange revealed that during soil thawing and snowmelt in March–April (see Fig. S4), fluxes were influenced by both concurrent soil temperature (Fig. 5a) and the thermal history (Figs. 5b and 4c) of the preceding four months. Concurrent soil temperature dominated fluxes during the highest non-growing season emissions in November–December (38 % ± 5.2; Fig. 5a). At the onset of the non-growing season (September–October), CO₂ exchange was primarily driven by PAR (58 % ± 3.1 %; Fig. 5a). During January–February, when soils were typically frozen or snow-covered (see Fig. S4), no single driver clearly dominated (Fig. 5a). During the productive season, fluxes were mainly driven by PAR (43 % ± 3.4 %; Fig. 5a) and WTD (36 % ± 3.2 %; Fig. 5a). The soil thermal conditions in the last 60 d were also an important driver during the onset of the GS (Fig. 5b).

Overall, model performance was moderate to high (R² = 0.38–0.71) and varied depending on the soil temperature representation. The best-performing model for November–December and January–February periods used concurrent soil temperature (R²= 0.38 and R²= 0.53, respectively; model used in Fig. 5a), for March–April and May–June the 2–4 month lagged soil temperature performed best (R²= 0.52 and R²= 0.67, respectively; model used in Fig. 5c), and for July–August the 60 d running mean showed the highest explanatory power (R²= 0.71; model used in Fig. 5b), while differences between representations were small in September–October (R²= 0.62–0.65).

Figure 5 Relative importance of photosynthetically active radiation (PAR), vapor pressure deficit (VPD), water table depth (WTD), and soil temperature at 5 cm in predicting daily net ecosystem exchange (NEE), based on random forest models. All variables and NEE were aggregated to daily means. Shading represents standard deviation. Due to incomplete PAR records at the Siikaneva site, data from the Hyytiälä Forest Station were used for 2009–2015. Subplots differ in the soil temperature aggregation used: (a) Daily mean soil temperature at 5 cm; (b) Mean soil temperature at 5 cm, averaged over the two months preceding each flux measurement; (c) Mean soil temperature at 5 cm, averaged over a two-month window and lagged by 60 d (i.e., a lagged two-month mean).

4 Discussion

Our aim is to examine how much the non-productive season contributed to annual CO₂ exchange at the Siikaneva fen and to describe the seasonal patterns of CO₂ fluxes during this period. Additionally, we investigated the main environmental factors driving CO₂ emissions specifically focusing on concurrent and lagged soil temperature as a driver. The year 2016 emerged as a key focus because it featured an unusually large CO₂ release during the spring thaw.

4.1 Annual ecosystem carbon exchange and its seasonal variation through the year

On average, the studied fen ecosystem functioned as a CO₂ sink, with an annual net carbon uptake of −51 ± 39 g C m⁻². This is more than twice the mean of boreal wetlands based on a recent synthesis (−17 g C m⁻² yr⁻¹; Virkkala et al., 2021), indicating it is a highly productive wetland site. Emissions during the non-productive season offset approximately 60 % of the carbon sequestered during the following productive season (Fig. 3A), a proportion comparable to that observed in other wetland sites (Virkkala et al., 2021; Aurela et al., 2002; Wang et al., 2023; López-Blanco et al., 2018; Yao et al., 2022), although the studies differ in the definition used for the non-growing season.

The non-productive season lasted on average from mid-September through the end of April (Fig. 2a), and was a month shorter than the thermal non-growing season (Fig. 2). The 5 °C threshold used to define the thermal growing season was a better predictor of the onset of the growing season than of its end, indicating a closer link between temperature conditions and spring onset than autumn termination. The two years, 2016 and 2018, with a positive annual CO₂ balance (i.e., when the ecosystem acted as a net source of CO₂ rather than a sink) showed strong deviations between the thermal and productivity-based season definitions. In 2016, the non-productive season persisted unusually late into early June (ending on 3 June), whereas in 2018 it began exceptionally early, already on 30 July. The unusually early onset of the non-productive season in 2018 has previously been linked to extreme summer drought conditions (Rinne et al., 2020). However, it is important to note that this early start of the non-productive season does not necessarily indicate a lack of photosynthesis or actual plant growth, but rather a net carbon loss at the ecosystem level.

Emissions generally declined from autumn to spring (Fig. 3b): mean fluxes from October to December (8.4 g C m⁻² month⁻¹) were nearly twice those observed from January to March. Autumn is typically the period with the highest emissions (Byrne et al., 2022; Commane et al., 2017), which can potentially shift the system to an annual net carbon source. However, at boreal sites, growing season uptake generally offsets autumn emissions, even under current warming trends (See et al., 2024). PAR emerged as a strong driver during September–October, a period when the system typically shifted from productive to non-productive. This transition is highly sensitive to declining light availability, which influences the onset and progression of senescence. As a result, day-to-day variability in carbon uptake is largely determined by PAR. During the peak emission period from November to December (Fig. 3b), carbon exchange was mainly controlled by soil temperature (Fig. 5a), which is typically the dominant driver of non-growing season emissions (Natali et al., 2019).

Emissions were lowest in January and February (Fig. 3b), with low interannual variability and high uncertainty in concurrent environmental drivers (Fig. 5a). This uncertainty may partly arise from a temporal lag between environmental conditions, such as soil temperature, and the observed fluxes, as snow and ice can restrict gas diffusion, allowing CO₂ to build up in the soil and snowpack before being released later (Martz et al., 2016; Morgner et al., 2010). The 60 d averaged soil temperature was also not identified as an important driver (Fig. 5b), indicating that midwinter fluxes are not consistently higher or lower following warmer or colder prior soil conditions.

From March to April, random forest analysis highlighted both concurrent soil temperature and the preceding months' conditions (e.g., the 60 d running mean and December–January soil temperatures) as important predictors of fluxes (Fig. 5). This indicates a legacy effect of late-autumn and early-winter temperatures on spring emissions, consistent with the hypothesis that CO₂ accumulates during frozen conditions and is released upon thaw (Raz-Yaseef et al., 2017; Sullivan et al., 2012). In addition, thermal conditions in autumn and winter may influence microbial activity and substrate availability during thaw.

4.2 Spring 2016 respiration enhanced by winter conditions and snowmelt

In 2016, an exceptional respiration pulse during April–May dominated the annual CO₂ budget (Fig. 4a). Lasting about six weeks, this event offset 37 % of the subsequent productive season uptake and produced the highest annual CO₂ loss of the record. Comparable pulses have been reported from other cold-climate wetlands (Wang et al., 2023) and permafrost-affected tundra sites (Raz-Yaseef et al., 2017; Arndt et al., 2020). Raz-Yaseef et al. (2017) and Wang et al. (2023) associated the pulse with the rapid release of gases that accumulated beneath snow and ice during winter. A similar event was reported by Arndt et al. (2020), who further suggested that snowmelt infiltration and oxygen supply can stimulate microbial activity, enhancing respiration during thaw. Our results add to this evidence but are, to our knowledge, the first to document such an event in a non-permafrost site. Overall, our observations indicate that the 2016 events combined the release of CO₂ accumulated during winter with an additional increase in respiration during spring thaw, together enhancing total emissions.

During snowmelt, CO₂ and CH₄ emissions rose sharply while surface soil temperatures remained near 0 °C, suggesting release of stored gases from the snowpack and soil, possibly aided by ice cracking (Raz-Yaseef et al., 2017). The conditions during the previous autumn and winter were favourable for gas buildup in soil. Phenocamera observations show the presence of an ice layer during winter (see Fig. S10), which appeared more extensive than in other years and likely formed after a sudden cold spell in January following a warm autumn. This ice layer would have restricted gas diffusion, allowing CO₂ and CH₄ produced at depth to accumulate. Warm conditions in late autumn 2015 and early winter kept deep soil layers warmer than average throughout the winter months (Fig. 4b), enhancing CO₂ and CH₄ production during the non-productive season and increasing the pool of gases available for release once transport pathways reopened. The random forest analysis confirmed a soil temperature legacy effect, with late-autumn and early-winter conditions influencing fluxes during spring thaw (Fig. 5c).

After snow cover disappeared, CH₄ fluxes declined while CO₂ emissions remained elevated, indicating a shift from transport-dominated fluxes to enhanced respiration. This increase in respiration was likely driven by the decomposition of labile organic matter made available through freeze–thaw damage to the organic matter in the soil. Soil freeze-up can create mechanical disturbance and break down biomass in the soil, which could promote faster decomposition in the following spring (Brooks et al., 2005; Byun et al., 2021). Freezing can also damage the winter microbial populations, which releases easily degradable carbon and nutrients to the surviving community, boosting respiration post thaw (Brooks et al., 2005). In addition to responses related to soil freezing, infiltration of oxygen-rich snowmelt water has been shown to stimulate aerobic microbial activity at depth (Arndt et al., 2020), further amplifying respiration during thaw.

Given the severity of the cold spell in January and the formation of a surface ice layer, it is possible that vegetation was affected by these extreme conditions, potentially delaying the onset of photosynthetic activity and contributing to the positive high NEE observed in spring. To evaluate whether elevated NEE was driven by increased respiration or reduced photosynthetic uptake, we examined GPP dynamics (Fig. S9). As GPP followed the long-term average and exceeded 1 g C m⁻² d⁻¹ already on 6 May, four weeks before NEE turned permanently negative, this indicates that the onset of photosynthetic activity was not delayed. Therefore, the persistently high positive NEE values observed in spring were not primarily driven by reduced photosynthetic uptake but instead reflect elevated ecosystem respiration. Additionally, the thermal growing season started on 27 April (Fig. 2b), one day after the long-term mean, indicating that spring thermal conditions were within the typical range and therefore unlikely to have constrained the onset of ecosystem activity or reduced early-season carbon uptake.

Beyond the biological explanations for the 2016 event, a physical process may also be relevant. Campeau et al. (2021) found that weakened thermal stratification in autumn promotes turbulent diffusion that has the potential to release large amount of CO₂ stored in deep porewater to the surface. In our data, temperature fluctuations at 50 cm depth (see Fig. S2) during winter while the upper layers remained frozen and stable suggest hydrological mixing within saturated peat. The subsequent CO₂ pulse coincided with a period of weak thermal stratification, when the surface and deeper soil layers reached similar temperatures, conditions that could have supported the rapid release of stored CO₂. The pulse may therefore have been supported not only by microbial processes but also by transient physical mixing and mobilization of deeper CO₂. Although direct evidence of porewater CO₂ dynamics is lacking, this temporal correspondence suggests that biological and physical pathways may have interacted to produce the unusually high CO₂ efflux in 2016. More broadly, these findings highlight the need to consider CO₂ accumulation in the peat profile and both vertical and lateral transport as a fully integrated net ecosystem carbon balance (López-Blanco et al., 2025).

Such events appear to be infrequent but can have outsized impacts on the annual carbon balance. In 17 years of monitoring, 2016 was the only year when such a pronounced spring pulse was observed. The event suggests that no single environmental driver can account for its magnitude; rather, several conditions (e.g., warm soils at depth, surface freeze-up) likely needed to coincide for it to occur. To better understand these dynamics, direct monitoring of peat gas buildup during the non-growing season would be valuable for determining whether spring pulses primarily reflect the release of gases accumulated under frozen conditions or enhanced microbial respiration initiated during thaw.

4.3 Vulnerability of northern wetlands CO₂ balance to NGS changes

The dominance of non-productive season processes in affecting the annual CO₂ balance means that warming late autumns, earlier spring thaws, or increased frequency of freeze–thaw events could substantially alter the net sink strength of boreal fens. Projected late-season warming may prolong microbial activity and amplify NGS losses (Natali et al., 2019) unless offset by concurrent growing season gains (See et al., 2024).

In this time series, two years (2016 and 2018) shifted the fen from a net CO₂ sink to a source (Figs. 3 and A1). In 2018, a summer drought reduced uptake during the productive season and caused an earlyonset of the non-productive season (Rinne et al., 2020). While this extended the non-productive season to its longest duration in the record, total productive season emissions remained smaller than in 2016, demonstrating that short-lived but intense events during the non-productive part of the year can outweigh more gradual seasonal anomalies in determining the annual carbon balance.

Representation of episodic emissions during spring thaw and autumn freeze is currently lacking in process-based models and future projections (e.g., Watts et al., 2023). Given that high, episodic carbon release during the non-productive season resulted in the highest annual net loss of CO₂ in this study, models aiming to predict peatland carbon–climate feedbacks must explicitly represent the thermal legacy effects from autumn to spring and soil gas accumulation and release mechanisms during freezing and thawing. Without such process-level representation, predictions may underestimate the potential for short-term non-growing season events to transform high-productivity peatlands from strong sinks into net sources of CO₂.

5 Conclusions

Our study highlights the importance of non-growing season carbon exchange in shaping the annual carbon balance of boreal peatlands. While the fen generally acted as a CO₂ sink, episodic events – such as the exceptional spring pulse in 2016 – demonstrated that infrequent disturbances can strongly influence the annual carbon balance.

Our results point to the significance of soil thermal history, potential peat gas accumulation and increased respiration following soil freezing as key mechanisms driving these pulses. These processes remain poorly constrained and are not yet well represented in ecosystem models, but they are critical for predicting peatland carbon exchange. Long-term, continuous measurements proved particularly valuable: our 17-year record not only captured infrequent but influential events such as the respiration pulse in 2016 but also provided the broader context to recognize their rarity and evaluate their contribution to interannual variability.

Autumn accounted for the largest share of emissions, emphasizing the role of shoulder seasons in offsetting growing-season uptake. Consequently, warmer autumns and predicted reduced snow cover are likely to enhance CO₂ efflux by prolonging the period of unfrozen soils and stimulating late-season microbial activity.

Looking forward, warming autumns and changes in snow cover will alter soil thermal conditions, which can affect both concurrent CO₂ fluxes during the non-growing season and fluxes later in the season. Understanding these processes – including soil thermal legacies and potential subsurface gas buildup – is therefore essential for assessing the vulnerability of peatland carbon storage under future climate change.

Appendix A

Figure A1 (a) Cumulative sums of gap-filled daily sums of net ecosystem exchange of CO₂ (NEE), starting from the day of year (DoY) 1 and (b) starting from DoY 57. As 2016 had the highest annual balance, but there was a significant measurement gap in DoY 1–57 in 2016, we calculated the annual sums for all years excluding this period. 2016 remained to have the highest balance, and the relative ranking in the annual balances remained similar. We therefore concluded that the positive balance of 2016 was not a result of the gapfilling period.