Carbon dioxide and methane exchange of a patterned subarctic fen during two contrasting growing seasons

The patterned microtopography of subarctic mires generates a variety of environmental conditions, and carbon dioxide (CO₂) and methane (CH₄) dynamics vary spatially among different plant community types (PCTs). We studied the CO₂ and CH₄ exchange between a subarctic fen and the atmosphere at Kaamanen in northern Finland based on flux chamber and eddy covariance measurements in 2017–2018. We observed strong spatial variation in carbon dynamics between the four main PCTs studied, which were largely controlled by water table level and differences in vegetation composition. The ecosystem respiration (ER) and gross primary productivity (GPP) increased gradually from the wettest PCT to the drier ones, and both ER and GPP were larger for all PCTs during the warmer and drier growing season 2018. We estimated that in 2017 the growing season CO₂ balances of the PCTs ranged from −20 g C m⁻² (Trichophorum tussock PCT) to 64 g C m⁻² (string margin PCT), while in 2018 all PCTs were small CO₂ sources (10–22 g C m⁻²). We observed small growing season CH₄ emissions (< 1 g C m⁻²) from the driest PCT, while the other three PCTs had significantly larger emissions (mean 7.9, range 5.6–10.1 g C m⁻²) during the two growing seasons. Compared to the annual CO₂ balance (−8.5 ± 4.0 g C m⁻²) of the fen in 2017, in 2018 the annual balance (−5.6 ± 3.7 g C m⁻²) was affected by an earlier onset of photosynthesis in spring, which increased the CO₂ sink, and a drought event during summer, which decreased the sink. The CH₄ emissions were also affected by the drought. The annual CH₄ balance of the fen was 7.3 ± 0.2 g C m⁻² in 2017 and 6.2 ± 0.1 g C m⁻² in 2018. Thus, the carbon balance of the fen was close to zero in both years. The PCTs that were adapted to drier conditions provided ecosystem-level resilience to carbon loss due to water level drawdown.

1 Introduction

Northern mires have sequestered substantial amounts of atmospheric carbon (C) since the last glacial period. The C storage of these peat soils has been estimated to be 415 ± 150 Pg of C (Hugelius et al., 2020), which adds up to about 30 % of the global soil C. This C storage has accumulated through the photosynthetic fixation of carbon dioxide (CO₂) by mire vegetation, which in the long term has been larger than the release of C through plant respiration and peat decomposition. In the short term, however, the C balance of a mire can switch from a sink to a source, as the annual C accumulation rate is sensitive to variations in moisture conditions and temperature (Alm et al., 1999; Bubier et al., 2003; Lindroth et al., 2007; Olefeldt et al., 2017) and to the length of the snow-free period (Aurela et al., 2004; Lund et al., 2012). Understanding the annual variability in peatland C dynamics is essential, as the subarctic and arctic regions warm rapidly, 2–3 times as fast as the rest of the world (Masson-Delmotte et al., 2018). This is projected to result in increased evapotranspiration and altered precipitation patterns, affecting in turn the C balance of mires (Tarnocai, 2006).

Subarctic mires endure long winters and relatively short growing seasons and have near-zero mean annual air temperatures. The net response of ecosystem C exchange to annual weather conditions depends on multiple processes and is thus hard to predict. In particular, the timing of soil thaw and snowmelt has been shown to impact the growing season length and consequently the annual net CO₂ uptake of mires, because earlier springs advance the timing of bud burst (Aurela et al., 2004; Grøndahl et al., 2008). In addition, warm springs increase microbial activity, ecosystem respiration (Lafleur et al., 2005) and, where anoxia prevails, methanogenesis and methane (CH₄) emissions to the atmosphere (Kim et al., 1999).

Ecosystem respiration, gross primary production and net CH₄ production in mires depend strongly on water table level (WTL), which determines the depth of oxic and anoxic layers in the peat. Microforms with varying WTL, e.g. hummocks and hollows, create different habitats for plant and microbe communities (Wieder et al., 2006). Peat decomposition is greater in the oxic layer above the WTL (Silvola et al., 1996), and deeper oxic layers also enhance plant production as roots require aeration (Bubier et al., 2003). Consequently, the net CO₂ exchange in northern mires not only varies among different sites (Bubier et al., 1998; Frolking et al., 1998; Lindroth et al., 2007) and from year to year (Aurela et al., 2004; Lund et al., 2012) but also among the microforms and plant communities within a site (Bubier et al. 1998; Alm et al. 1999; Heikkinen et al., 2004). However, there does not seem to be a clear universal spatial pattern in relation to changing moisture conditions (Alm et al., 1999; Waddington and Roulet, 2000; Laine et al., 2007a; Strack and Waddington, 2007; Maanavilja et al., 2011; Korrensalo et al., 2019). In contrast, CH₄ emissions are generally larger from wet than dry plant communities, as they depend on the balance between microbial production in anoxic conditions and oxidation above the WTL (Saarnio et al., 1997; Segers, 1998; Alm et al., 1999; Heikkinen et al., 2004; Wieder et al., 2006; Laine et al., 2007b), but there is additionally marked small-scale spatial variation related to nutrient and substrate availability and plant species composition (Svensson and Rosswall, 1984; Kettunen, 2003; Christensen et al., 2004; Dorodnikov et al., 2011).

High temperatures with water level drawdown have been reported to decrease CO₂ uptake (Chivers et al., 2009; Munir et al., 2014) and CH₄ emissions (Peltoniemi et al., 2016; Olefeldt et al., 2017) of mires. During the summer of 2018, a large-scale heatwave and drought took place in north-western Europe, including northern Finland (Lehtonen and Pirinen, 2019a, b). Rinne et al. (2020) found that this drought reduced CO₂ uptake and CH₄ emissions, as compared to a reference year, on most of the Fennoscandian mires studied. However, the magnitude of this effect varied among the mires, as did the duration and severity of the drought.

In this study, we examine the CO₂ and CH₄ exchange between the subarctic Kaamanen fen and the atmosphere during two contrasting years (2017 and 2018). The site was included in the synthesis of Rinne et al. (2020) that applied spatially averaged eddy covariance (EC) data. However, the annual variation in the plant-community-scale C exchange has not been investigated previously, and it is unknown how different communities react to water level drawdown during drought events. Thus, our specific objective is to study the small-scale CO₂ and CH₄ flux variation in response to moisture conditions and compare this to the ecosystem-scale response. We (1) utilise the EC flux measurement technique to detect the ecosystem-scale variation in C exchange, (2) study how C exchange varies spatially and temporally among different plant communities by using manual flux chamber measurements and (3) determine the main environmental factors controlling the C fluxes by means of a linear mixed-effects model. We will also put our results in the context of the earlier gas exchange data from this measurement site.

2 Materials and methods

2.1 Study site

The study took place at a patterned mesotrophic fen at Kaamanen in northern Finland (69^∘8.435^′ N, 27^∘16.189^′ E; 155 m a.s.l.). The average annual mean temperature at the Inari Ivalo weather station, 59 km south of Kaamanen, during the 30-year reference period of 1981–2010 was −0.4 ^∘C, and the corresponding mean annual precipitation sum was 472 mm (Pirinen et al., 2012). A major part of the fen is patterned with strings and flarks (or hummocks and hollows, respectively). Hummocks of a few metres' width surround wet flarks and form strings, up to a few tens of metres in length that sprawl through the fen. The strings are elevated from the water table by 0.3–0.8 m and contain ice lenses that can remain frozen until late summer (Aurela et al., 2001). The site has no permafrost, even though it is located 300 km north of the Arctic Circle, and the nearest isolated permafrost palsas are currently about 50 km north of the site. The peat depth within the study area is 1–2 m (Piilo et al., 2020).

The vegetation within the patterned fen, which is the focus of our measurements, can be divided into four main plant community types (PCTs): (1) Ericales–Pleurozium on the dry string tops (ST), (2) Betula nana–Sphagnum on the string margins (SM), and (3) Trichophorum tussock (TT) and (4) flark (F) communities in the wet hollows (Fig. 1, Table 1) (Maanavilja et al., 2011). The areal coverage of different land cover types at the site has been estimated with high spatial resolution remote sensing by Räsänen and Virtanen (2019) (Table 1). In addition to the low vegetation, a few sporadic Betula pubescens and Pinus sylvestris trees are growing on the driest parts of the fen.

A stream flows through the fen from north to south, but there is also nearly continuous surface water flow across the peatland along the elevation gradient. The fen is flooded during the spring thaw, but the magnitude of flooding varies annually. The string tops, however, are separated from the water table of rest of the fen (Fig. 1) and therefore mainly receive nutrients from precipitation instead of the lateral inflow.

The C exchange of the fen has been studied extensively during the past few decades (Aurela et al., 1998, 2001, 2002, 2004; Hargreaves et al., 2001; Laurila et al., 2001; Heikkinen et al., 2002; Maanavilja et al., 2011; Kross et al., 2016; Rinne et al., 2020; Piilo et al., 2020). The ecosystem–atmosphere CO₂ exchange has been measured with the EC technique continuously since 1997 and CH₄ exchange since 2010. On average, the fen has been estimated to be a small (approximately −15 g C m⁻² yr⁻¹) atmospheric C sink during the measurement periods 1997–2002 and 2011–2016 (Aurela et al., 2004, and unpublished data), but the C exchange varies among the microforms: the driest plant communities act as the weakest net CO₂ sink (Heikkinen et al., 2002; Maanavilja et al., 2011) and CH₄ source (Heikkinen et al., 2002).

Table 1 Vegetation composition of the four main fen plant communities and their areal coverage inside a 200 m radius around the eddy covariance tower.

Figure 1 The four dominant plant communities (Table 1) within the main eddy covariance source area and the other land cover types at the Kaamanen site. The land cover map is centred at the EC tower. The schematic cross section of the fen microtopography shows the average position of the studied plant community types (based on Maanavilja et al., 2011).

2.2 Flux measurements

2.2.1 Chamber measurements

A total of 17 chamber flux measurement plots were chosen to represent the four PCTs: five plots for the SM communities and four plots each for the ST, TT and F communities. Eight of the 17 aluminium collars (60 cm × 60 cm) were installed during the first days of June 2017, during the soil thawing, to accompany the collars that were already installed previously in 2006. The overall vegetation condition and species composition inside the old collars were checked to match the new study plots. The collars were positioned within 50 m from the instrument booth (marked with a star in Fig. 1).

The CO₂ and CH₄ fluxes between the ecosystem and atmosphere were measured with manual flux chambers biweekly during the growing seasons of 2017 and 2018: six times between 12 June and 11 October 2017, and seven times between 31 May and 4 September 2018. All chamber measurements were conducted between 09:00 and 16:00 local standard time.

The measurements were conducted with a transparent polycarbonate chamber (width × depth × height of 60 cm × 60 cm × 30 cm in 2017 and 60 cm × 60 cm × 40 cm in 2018). The chamber was connected with a 50 m long inlet tube (Teflon, inside diameter 3.1 mm) to a closed-path infrared gas analyser (Picarro G2401, Picarro Inc, USA) to detect the changes in CO₂, CH₄ and H₂O mixing ratios in its airspace. The chamber air was mixed with a battery-driven fan. The chamber closure time for each measurement was 2 min. Air temperature inside the chamber and soil temperature at the 10 cm depth was measured at each plot on the chamber flux measurement days with IKES Pt100 sensors. Soil moisture was measured at the plots with a ML3 ThetaProbe sensor (Delta-T Devices Ltd., England) simultaneously with the chamber flux measurements. Photosynthetic photon flux density (PPFD) was measured during the chamber closures with PQS1 PAR Quantum sensor (Kipp&Zonen, Delft, the Netherlands) on top of the chamber.

The flux was determined by measuring the CO₂ and CH₄ mixing ratio change within the chamber, first in ambient light, then in one or two reduced light conditions and lastly in complete darkness; the amount of incoming solar radiation was reduced by 40 %–50 %, 75 %–90 % and 100 %, respectively. The chamber was lifted off the collar between the measurements to restore the ambient gas concentration inside the chamber. Due to a small lag, until the data can be accepted for CO₂ and CH₄ flux calculation, only the data from the final 1.5 min of the 2 min closure time were used.

The CO₂ and CH₄ fluxes were calculated as

$$\begin{array}{}\text{(1)}& F_x={\displaystyle \frac{p\times M_x\times V}{R\times T\times A}}\times {\displaystyle \frac{\mathrm{d}{c}_x}{\mathrm{d}t}},\end{array}$$

where p is atmospheric pressure, M_x is the molecular mass of CO₂ (44.01 g mol⁻¹) or CH₄ (16.04 g mol⁻¹), R is the universal gas constant (8.314 J mol⁻¹ K⁻¹), T is the mean air temperature during chamber closure, V is the chamber volume, A is the chamber base area, and $\frac{\mathrm{d}{c}_x}{\mathrm{d}t}$ is the mean CO₂ or CH₄ mixing ratio change in time calculated with linear regression based on ordinary least squares. The mixing ratio is expressed with respect to dry air, so no correction for water vapour dilution was necessary. A micrometeorological sign convention was used: a positive flux indicates a flux from the ecosystem to the atmosphere (emission), and a negative flux indicates a flux from the atmosphere into the ecosystem (uptake).

For estimating the PCT-specific CH₄ flux time series from the chamber measurements, a mean flux was calculated for each of the 17 chamber plots from the two to four chamber closures that were conducted during each measurement day. The measurement data were screened to ensure that they were unaffected by CH₄ ebullition events and disturbances induced by a closing chamber. The criteria for discarding measurements were as follows: CH₄ mixing ratio was > 5 ppm at the start of the closure, or the normalised root mean square error of the linear regression fit was > 0.02, or there was an obvious non-linearity in the time series. The number of rejected/total data for each PCT was $\mathrm{27}/\mathrm{148}$ (F), $\mathrm{4}/\mathrm{150}$ (TT), $\mathrm{41}/\mathrm{188}$ (SM) and $\mathrm{4}/\mathrm{147}$ (ST).

2.2.2 Eddy covariance measurements

The EC measurements were conducted on a tower 5 m above the mean fen surface. The EC system consisted of a three-dimensional anemometer (USA-1, METEK Meteorologische Messtechnik GmbH, Germany), a closed-path infrared gas analyser for CO₂ and H₂O mixing ratios (LI-7000, LI-COR Biosciences, USA) and a laser-based gas analyser for CH₄ mixing ratio (RMT-200, Los Gatos Research, USA). The heated 6 m inlet tubes for the gas analysers were mounted 0.3 m below the centre of the anemometer sound paths. The inner tube diameter and the flow rate were 3.1 mm and 6 L min⁻¹, and 8 mm and 15 L min⁻¹ for the LI-7000 and RMT-200, respectively.

The EC data were sampled at 10 Hz, and standard methods were used to calculate half-hourly turbulent fluxes (Aubinet et al., 2012). Block averaging and a double rotation of the coordinate system were applied first (McMillen, 1988). Water vapour fluctuations affecting CO₂ mixing ratios were compensated for (Webb et al., 1980), but a similar procedure was not necessary for temperature (Rannik et al., 1997). Flux losses due to high-frequency signal attenuation were corrected for using methods detailed by Moore (1986) and Tuovinen et al. (1998).

The half-hourly averaged data were screened, and the data were accepted on the basis of the following criteria: relative stationarity < 100 % (Foken and Wichura, 1996), number of recorded data per 30 min > 17 400, number of signal spikes per 30 min < 360, mean CO₂ mixing ratio of 340–550 ppm. A western wind direction sector (260^∘–315^∘), within which the ecosystem changed from fen to pine forest at a distance of 100 m, was excluded. In addition, periods of insufficient turbulence were discarded with the friction velocity limit of 0.1 m s⁻¹.

2.2.3 Abiotic and biotic environmental measurements

Additional meteorological variables measured close to the EC tower included air temperature and humidity at 3 m height (Vaisala HMP 230), global and reflected radiation (Kipp&Zonen CM7), and downward and upward photosynthetic photon flux density (Kipp&Zonen PQS 1). Water vapour pressure deficit (VPD) was calculated from air temperature and relative humidity according to Jones (2013). Soil temperature profiles were measured in both a string (at −10, −30, −50, −75 and −105 cm) and a flark (at −10, −30 and −50 cm) (IKES Pt100 sensors). The data were collected continuously by data loggers as 30 min averages. A soil temperature time series was generated for each chamber plot by adopting the −10 cm flark temperatures for F, TT and SM and the −10 and −30 cm string temperatures for ST, and adjusting them to match the plot-specific soil temperatures on the chamber flux measurement days.

The water table level relative to the peat surface was measured from perforated tubes placed next to each chamber measurement plot. During early summer, when there were still ice lenses inside strings, the WTL of string tops was measured as the depth of an ice lens or the meltwater overlaying the lens. These measurements were conducted simultaneously with the chamber flux measurements.

Plant species coverage and mean height were measured biweekly in each collar for estimating the leaf area index (LAI). Harvested samples of different species groups (deciduous shrubs, evergreen shrubs, forbs, graminoids and mosses) were collected during peak summer 2017 in 57 sampling plots of 50 cm × 50 cm, and LAI was measured from these samples with an A4 scanner. With the help of the samples, empirical relationships between LAI and species group coverage and height were established with ordinary least squares linear regressions (Juutinen et al., 2017). Biweekly LAI was then estimated for the collars with these empirical relationships.

The phenology of the fen vegetation was also tracked utilising daily phenocamera images taken with StarDot Netcam SC 5 digital camera. The camera was placed at 3 m in a weather-proof housing on a pole facing the north, and the viewing angle of the camera was adjusted to 45^∘. Images were automatically taken every 30 min, and daytime images during the growing season were used for the image analysis. The processing of the images was executed with FMIPROT, a software designed for image processing for phenological and meteorological purposes (Linkosalmi et al., 2016; Tanis et al., 2018). The green chromatic coordinate (GCC) was used as a greenness index and calculated as

$$\begin{array}{}\text{(2)}& \mathrm{GCC}={\displaystyle \frac{\mathrm{\Sigma }G}{\mathrm{\Sigma }G+\mathrm{\Sigma }R+\mathrm{\Sigma }B}},\end{array}$$

where ΣG, ΣR and ΣB are the sums of green, red and blue channel indices, respectively, of all pixels comprising a region of interest (ROI). In addition to a more general view, ROIs were defined separately for flark and string PCTs by grouping F and TT areas and SM and ST areas together, respectively. The growing season start dates were determined based on soil temperatures measured at 10 cm depth in strings and flarks, with temperatures rising over +3 ^∘C indicating the start of the growing season. The end dates were defined in the strings from the timing of soil freezing at 10 cm depth and in flarks from the appearance of continuous snow cover seen in the daily phenocamera images.

To estimate peat bulk density and C and N concentrations for the different PCTs, peat samples were collected from 40 typified plots placed at distances of 25 to 150 m from the EC tower in cardinal, intercardinal and secondary intercardinal directions (placement of plots is described in Räsänen and Virtanen, 2019). At each plot, a sample of approximately 5 cm × 5 cm × 5 cm was cut out of the peat at 0–5 cm depth (i.e. straight under the litter layer, where plant structures are still discernible) and at 15–20 cm depth. The samples were dried for bulk density estimates and ground using a ball mill, and their C and N concentrations were analysed using a CNS-2000 analyser (LECO Corporation, Saint Joseph, MI, USA). Peat C and N content (mg cm⁻³) was then calculated using the bulk density and C and N concentrations. The pH was measured in the field using a sample of water collected from the bottom of a 30 cm deep hole at each plot.

2.3 Partitioning and gap filling of CO₂ fluxes

2.3.1 Environmental response functions

To fill the gaps in the collected CO₂ flux data, both EC and chamber based, and to analyse the processes controlling net ecosystem exchange (NEE), the CO₂ fluxes (F_NEE) were partitioned to two opposite flux components:

$$\begin{array}{}\text{(3)}& F_{\mathrm{NEE}}=F_{\mathrm{GPP}}+F_{\mathrm{R}},\end{array}$$

where F_GPP is the negative flux due to gross primary productivity (GPP), which represents the CO₂ uptake by the vegetation through photosynthesis, and F_R is the positive flux due to ecosystem respiration (ER), which describes the release of CO₂ to the atmosphere through autotrophic and heterotrophic processes.

Gaps in the F_NEE time series were filled with parameterised values that were estimated separately for F_GPP and F_R. The dependence of F_GPP on solar radiation was parameterised by a rectangular hyperbola (e.g. Whiting, 1994):

$$\begin{array}{}\text{(4)}& F_{\mathrm{GPP}}={\displaystyle \frac{\mathrm{PPFD}\times \mathit{\alpha }\times {\mathrm{GP}}_{\max }}{\mathrm{PPFD}\times \mathit{\alpha }+{\mathrm{GP}}_{\max }}},\end{array}$$

where PPFD is the measured photosynthetic photon flux density, α is the initial slope between F_GPP and PPFD, and GP_max is the theoretical maximum gross photosynthesis rate.

The respiration flux was parameterised by an exponential dependence on temperature (Lloyd and Taylor, 1994):

$$\begin{array}{}\text{(5)}& F_{\mathrm{R}}=R_{\mathrm{10}}\times e^{E_{\mathrm{0}}\left(\frac{\mathrm{1}}{T_{\mathrm{0}}}-\frac{\mathrm{1}}{T_{\mathrm{s}}-T_{\mathrm{1}}}\right)},\end{array}$$

where R₁₀ is the base respiration rate at 10 ^∘C, E₀ is the activation energy, T₀=56.02, T₁=227.13, and T_s is the soil temperature at the 10 cm depth.

The calculations and analyses were made with the Python programming language (Python Software Foundation, version 2.7, https://www.python.org, last access: 1 October 2020) with the NumPy (http://www.numpy.org/, last access: 1 October 2020) and SciPy (http://www.scipy.org/, last access: 1 October 2020) libraries.

2.3.2 Gross primary productivity parameterisation for chamber-based fluxes

For each of the chamber measurement plots, the GPP flux time series was calculated over the two growing seasons by utilising the dependence of photosynthesis rate on solar radiation. To estimate the F_GPP time series from the data from the six and seven measurement days in 2017 and 2018, respectively, a time-invariant, PCT-specific $k=\mathit{\alpha }/{\mathrm{GP}}_{\max }$ parameter was estimated. This was done by pooling the daily data from all plots of the same PCT, including only those data sets of different shading levels in which the highest PPFD exceeded 800 µmol m⁻² s⁻¹, and fitting the light response curve (Eq. 4) to these data. The PCT-specific k was calculated from the fitted α and GP_max values as the variance-weighted mean k, discarding the fits that had a relative error greater than 100 %. The k values obtained were around 0.002 (Table 2), which is a typical value for mesotrophic boreal fens (Bubier et al., 1999).

Table 2 PCT-specific k=α $/$ GP_max parameter (± standard error).

After determining the PCT-specific k values, GP_max was estimated for each plot measurement by fitting Eq. (4) (modified to incorporate k) to the corresponding CO₂ flux and PPFD data. GP_max was then linearly interpolated between the measurement days assuming GP_max of 0 at the start and end of the growing season. Finally, the half-hourly F_GPP time series were calculated for each plot from the PCT-specific k parameter and the time series of the plot-specific GP_max parameter and PPFD measurements.

2.3.3 Respiration parameterisation for chamber-based fluxes

A respiration flux time series was calculated for both growing seasons, separately for each PCT, by adopting the dark chamber measurements of CO₂ flux as respiration data and relating them to the T_s measured simultaneously at the plot. The PCT-specific E₀ values were obtained by fitting Eq. (5) to all respiration and temperature data of each PCT (Table 3).

Table 3 PCT-specific activation energy (E₀) (± standard error).

The PCT-specific R₁₀ values were then calculated for each measurement day by using Eq. (5) with these E₀ values, respiration data and the corresponding half-hour mean soil temperatures. The chamber-based CO₂ flux measurements did not cover the growing season start and end days. Therefore, the R₁₀ values for those days were estimated with Eq. (5) by utilising the CO₂ fluxes measured with the EC technique during a 2-week period around the growing season start and end dates. Continuous time series of half-hourly R₁₀ for each PCT were obtained by linear interpolation. Finally, the F_R time series over the growing seasons were calculated for each plot from the PCT-specific E₀ parameter and the time series of the PCT-specific R₁₀ and the plot-specific soil temperature time series.

2.3.4 Gap filling of CO₂ eddy covariance fluxes

The EC flux measurement time series had gaps due to equipment failures and quality control filtering applied during the post processing of data. The gaps in the CO₂ flux data were filled by modelled F_GPP and F_R values, calculated with parameterised Eqs. (4) and (5), respectively, which were fitted to measurements. The fitting was performed in a moving window of at least 5 and 15 d, for F_GPP and F_R, respectively, long enough to result in at least 30 half-hourly observations. The modelling was performed in two steps. First, the F_R was parameterised with the nighttime data (PPFD < 30 µmol m⁻² s⁻¹), and second, by utilising obtained respiration parameters (E₀ and R₁₀) the GPP parameters (α and GP_max) were obtained by fitting F_NEE to all available data. In total, 62 % and 63 % of the CO₂ flux data were gap filled in the time series of 2017 and 2018, respectively.

2.4 Gap filling of CH₄ fluxes

The CH₄ flux ($F_{{\mathrm{CH}}_{\mathrm{4}}}$) time series collected with chambers were gap filled separately for each plant community type by assuming an exponential temperature dependence (Marushchak et al., 2016):

$$\begin{array}{}\text{(6)}& F_{{\mathrm{CH}}_{\mathrm{4}}}=a_{\mathrm{10}}\times Q_{\mathrm{10}}^{\left(T_{\mathrm{s}}-T_{\mathrm{0}}\right)/T_{\mathrm{0}}},\end{array}$$

where a₁₀ is the CH₄ flux at 10 ^∘C, Q₁₀ is the temperature coefficient, T_s is soil temperature at 10 cm depth, and T₀=10 ^∘C. Including water table level as an additional explanatory variable to the model did not improve the model.

The Q₁₀ coefficient was determined for each PCT from the data set of measured daily CH₄ fluxes and soil temperatures from each plot that included both growing seasons (Table 4).

For each of the measurement days, the PCT-specific a₁₀ values were calculated by using Eq. (6) with the previously determined Q₁₀ coefficients and measured soil temperatures. For days in the growing season beginning and end, where there were no chamber measurements, the a₁₀ values were estimated by using the CH₄ fluxes measured with the EC technique during a 2-week period in the start and end of the growing season. The a₁₀ data were linearly interpolated between the measurement days to obtain a continuous half-hourly time series. The CH₄ flux time series for each PCT were calculated for both growing seasons by using the time-invariant Q₁₀ coefficients, a₁₀ time series and the continuously measured soil temperature at each plot.

Table 4 PCT-specific Q₁₀ coefficient (± standard error).

For filling gaps in the filtered EC time series of CH₄ fluxes, a simple moving average interpolation of the half-hour fluxes was used. A moving average window of ± 1, 2, 4, 8, 16 or 32 d was used depending on the length of the gaps in data. In total, 64 % and 70 % of the CH₄ flux data were gap filled in the time series of 2017 and 2018, respectively.

2.5 Estimating flux uncertainty

The uncertainty of the annual CO₂ and CH₄ balances obtained from the EC-based fluxes were estimated by taking account of the most significant error sources. The random error estimate included the statistical measurement error (E_meas) and the error caused by gap filling of missing data (E_gap) (Räsänen et al., 2017):

$$\begin{array}{}\text{(7)}& E_{\mathrm{meas}/\mathrm{gap}}=\sqrt{{\sum }_i{\displaystyle \frac{{\left(F_{i,\mathrm{obs}}-F_{i,\mathrm{mod}}\right)}^{\mathrm{2}}}{n_{\mathrm{obs}}}}}\sqrt{n_{\mathrm{obs}/\mathrm{gap}}},\end{array}$$

where F_obs is the half-hourly CH₄ or CO₂ flux that remained after all the filtering procedures, F_mod is the corresponding fitted value and $n_{\mathrm{obs}/\mathrm{gap}}$ is the number of observed or gap-filled data. This provides a conservative error estimate for E_meas, and for E_gap it includes the effect of random variability on the model fits (Aurela et al., 2002).

Additionally, the annual error due to friction velocity filtering (E_ustar) was estimated by recalculating the annual EC-based CO₂ and CH₄ balances with modified data sets that were screened with two additional friction velocity limits (0.05 and 0.15 m s⁻¹). E_ustar was calculated as the average deviation from the annual balance calculated with the optimal friction velocity limit (0.1 m s⁻¹) (Aurela et al., 2002).

The total uncertainty of the annual EC-based CO₂ and CH₄ balances was calculated as

$$\begin{array}{}\text{(8)}& E_{\mathrm{tot}}=\sqrt{E_{\mathrm{meas}}^{\mathrm{2}}+E_{\mathrm{gap}}^{\mathrm{2}}+E_{\mathrm{ustar}}^{\mathrm{2}}}.\end{array}$$

The uncertainties of the PCT-specific chamber-based fluxes (F_GPP, F_R and $F_{{\mathrm{CH}}_{\mathrm{4}}}$) were estimated by combining the uncertainty due to the estimated parameters (GP_max, R₁₀ and a₁₀) and the flux variation among the plots of each PCT. The uncertainty of F_NEE was calculated assuming that the uncertainties of F_GPP and F_R are independent. The uncertainty estimate was calculated for each half hour in the time series and further for the monthly and growing season sums.

We assumed that random errors in the response function fitting parameters are independent, so the standard error of the function f (either F_GPP, F_R or $F_{{\mathrm{CH}}_{\mathrm{4}}}$) was calculated as

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}_f=\sqrt{{\left({\displaystyle \frac{\partial f}{\partial a}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_a^{\mathrm{2}}+{\left({\displaystyle \frac{\partial f}{\partial b}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_b^{\mathrm{2}}},\end{array}$$

where a and b are R₁₀ and E₀ for F_R, GP_max and k for F_GPP and a₁₀ and Q₁₀ for $F_{{\mathrm{CH}}_{\mathrm{4}}}$, and σ_a and σ_b denote their standard errors, respectively.

2.6 Linear mixed-effects models

The effect of environmental variables on ER, GPP, NEE and CH₄ fluxes was evaluated with the linear mixed-effects (LME) model that was fitted by maximum likelihood. The chamber flux measurement data of the 13 measurement days from both years were used in this analysis, and all four PCTs were pooled together. Both logarithmically transformed and non-transformed response variables were tested, and the final model was chosen based on model residual plots. The F_R data were transformed logarithmically. Normalisation of F_R to 10 ^∘C was tested, but it did not improve the model performance. The F_NEE and F_GPP data were normalised to a common radiation level of PPFD of 1200 µmol m⁻² s⁻¹ (denoted as F_NEE1200 and F_GPP1200), which represent a near-optimal radiation level for photosynthesis in northern ecosystems (Laurila et al., 2001). For the CH₄ fluxes, daily mean values were used with a logarithmic transformation. Normalisation of $F_{{\mathrm{CH}}_{\mathrm{4}}}$ to 10 ^∘C was tested, but it did not improve the model.

The following fixed explanatory variables were tested in the models for F_R and F_NEE1200: T_s, WTL, GCC and total vascular LAI; for F_GPP1200: WTL, GCC, daily maximum VPD and vascular LAI; and for CH₄ flux: soil temperature (at −10 cm), WTL, GCC and LAI divided into four plant groups (deciduous shrubs, evergreen shrubs, forbs and graminoids). Highly cross-correlated explanatory variables (Pearson correlation coefficient r> $\left|\mathrm{0.7}\right|$) were excluded from the models (Table A5), and the number of variables for the final model was reduced with a backward stepwise procedure by minimising Akaike's information criterion value. In all regressions, the measurement plot was included as a random effect. To evaluate the relative impact of each explanatory variable, the standardised regression coefficients were calculated, while the marginal coefficient of determination ($R_{\mathrm{m}}^{\mathrm{2}}$) was used to quantify how much the fixed effects explain of the variance of the response variable. Data analyses were conducted in R (R Core Team, 2019) with the packages nlme (Pinheiro et al., 2019), MASS (Venables and Ripley, 2002) and reghelper (Hughes, 2020).

3 Results

3.1 Environmental conditions

Compared to long-term statistics measured at the Inari Ivalo weather station, the growing season of 2017 at Kaamanen fen had an average temperature but high precipitation sum, while the growing season 2018 was warm and dry (Fig. 2). The annual average temperature was close to the long-term value (−0.4 ^∘C) in both years (−0.6 ^∘C in 2017 and 0.4 ^∘C in 2018), but in 2018 the monthly means of May, July and November were clearly higher than the reference values by 4.4, 5.4 and 6.5 ^∘C, respectively (Fig. 2a). On average, the April–October period was considerably warmer in 2018 (7.7 ^∘C) than in 2017 (5.5 ^∘C) (Figs. 2 and 3a). In 2018, the daily mean temperatures rose to 10 ^∘C already in early May, while in 2017 such temperatures were not recorded until early June. In 2018, the daily mean exceeded 20 ^∘C on 12 d, and the maximum daily mean temperature of 25.9 ^∘C was recorded on 18 July. In comparison, the daily mean temperatures never rose over 20 ^∘C during 2017, the maximum being 19.3 ^∘C (28 July). In August, September and October, the differences between the years were not as large as in July.

The annual precipitation sum was higher in 2017 (499 mm) than the long-term average (472 mm) and the precipitation in 2018 (473 mm). However, there were large differences in the monthly precipitation sums between the years (Fig. 2b). In 2017, there was hardly any rain in May, while the precipitation sum of the summer months of June–August was 32 % higher than the 30-year average. On the other hand, the precipitation sum of July 2018 was only 34 mm, which is less than half of the 30-year average. This dry spell was followed by a rainy August and September.

A marked difference between the study years was observed in VPD (Fig. 3b), which serves as an indicator for drought events (Lindroth et al., 2007; Aurela et al., 2007). We defined drought as the period during which the daily maximum VPD exceeded 2 kPa. In 2017, this limit was not exceeded, while in 2018 it was exceeded in total on 13 days between 2 July and 1 August, with a maximum of 3.1 kPa observed on 18 July (Fig. 3b). The average daily maximum VPD during this period was 1.69 kPa, while during the same period in 2017 it was 0.94 kPa. The corresponding mean air temperatures were 18.5 and 14.2 ^∘C, respectively (Fig. 3a).

Figure 2 Monthly (a) mean air temperature and (b) precipitation sum in 2017 and 2018 at the Kaamanen fen, and 30-year averages (Pirinen et al., 2012) measured at the Ivalo weather station (68^∘36^′ N, 27^∘25^′ E; 59 km south of Kaamanen). The monthly values measured at Ivalo in 2017 and 2018 are marked with diamonds.

Soil temperatures (Fig. 3c) varied from day to day less than the air temperature. The temperatures measured within the string top (ST) communities were lower than those of flark (F), Trichophorum tussock (TT) and string margin SM until late summer due to the presence of ice lenses in strings. In F, TT and SM, the average soil temperature during May–September was 9 ^∘C in 2017 and 11 ^∘C in 2018, while in ST it was 6 ^∘C in 2017 and 8 ^∘C in 2018. The difference in the maximum daily soil temperatures between the years, about 2 ^∘C for F, TT and SM and 3 ^∘C for ST, was not as prominent as in the air temperatures.

The growing season, as defined by the 3 ^∘C soil temperature limit, began 2 weeks earlier in 2018 than in 2017 (Sect. 2.2.3) (Table 5). In spring, when the snow is melting, the string plant communities (SM, ST) are first exposed to direct sunlight, and therefore the growing season began a few weeks earlier in those plant communities. The growing season ended as the peat soil froze and continuous snow cover was established in October.

The greenness index GCC showed a clear seasonal pattern (Fig. 3d). The GCC variation during May–June 2017 and in May 2018 coincided with the soil temperature rise and snowmelt. The flark plant communities (F and TT) with mostly sedge and moss vegetation had a lower GCC than the string plant communities (SM and ST), which had a more diverse vegetation composition (Table 1). The drought impact on vegetation was clearly visible in the field in July 2018, and accordingly a higher maximum GCC was recorded in 2017 than in 2018.

Figure 3 (a, b) Daily mean air temperature, (c, d) daily maximum vapour pressure deficit, (e, f) PCT-specific daily mean soil temperatures at 10 cm (F indicates flark, TT indicates Trichophorum tussock, and SM indicates string margin) or 20 cm depth (ST indicates string top), and (g, h) greenness index of the flark (F, TT) and string (SM, ST) areas. The drought period (2 July–1 August 2018) is denoted with shading. In the green chromatic coordinate (GCC) plots (g, h), the “s” boxes indicate the start of growing season dates and the “e” boxes the growing season end.

Table 5 Growing season start and end dates.

There were significant microtopography-related differences among the PCT-specific WTL data (Fig. 4a, b). The peat surface of flarks was usually barely submerged, while the Trichophorum tussocks stuck out a few centimetres from the water. While WTL rose somewhat inside the strings, it remained at a depth of about 10 and 50 cm in the SM and ST communities, respectively. During the early growing season, the water table depth in ST was bounded by the ice lens depth. During July 2018, WTL dropped in all communities, which matched the drought period observed as an increased VPD. The LAI was systematically lowest in the wettest plant community type (i.e. F) and highest in the ST community type (Fig. 4c, d).

Figure 4 (a, b) Mean water table level and (c, d) average vascular leaf area index in the chamber plots in 2017 and 2018. F indicates flark, TT indicates Trichophorum tussock, SM indicates string margin, and ST indicates string top. Error bars represent the standard deviation within a plant community type. The drought period in July 2018 is denoted with shading.

Table 6 Mean (± standard deviation) peat pH measured at 30 cm depth and bulk density, C and N content and C:N ratio within the top 20 cm peat layer (measured from 0–5 and 15–20 cm depth) for each plant community type.

3.2 Ecosystem-level CO₂ and CH₄ fluxes

During the winter (October–April), the ecosystem respiration determined the magnitude of CO₂ exchange and the fen was a source of CO₂. As the spring advanced in April and May GPP gradually got started, but respiration still dominated the NEE, and the switch to a CO₂ sink took place on 16 June 2017 and 30 May 2018, the difference in timing reflecting the onset of the growing season (Fig. 5a, c, d, Table 5).

The largest difference between the years in the cumulative CO₂ balance was recorded in early July, when the fen had accumulated 30 g C m⁻² more in 2018 than in 2017 (Fig. 5d). This difference was generated in June, mostly due to the larger GPP in 2018. However, the increased cumulative uptake was offset later between 20 July and 9 August 2018, when GPP decreased and the fen momentarily became a CO₂ source at the end of the drought period (Fig. 5a).

In both years, the fen turned from CO₂ sink to a source in early September (Fig. 5a) even though the plants kept photosynthesising until October (Fig. 5d). After this switchover, the trajectories of the cumulative CO₂ balances were similar in 2017 and 2018 (Fig. 5d), leading to a small annual sink (< 10 g C m⁻²) in both years (Table 7).

The fen acted as a CH₄ source throughout the year (Fig. 5e, f). In winter, emissions were low (0.006 g C m⁻² d⁻¹), after which a distinctive CH₄ emission pulse was observed in May of both years. This pulse lasted about 1 week and took place 2 weeks earlier in 2018 than in 2017 due to the different time of the snowmelt (Fig. 5e).

The drought period in July 2018 reduced the CH₄ emissions compared to the previous year (Fig. 5e, f). However, the drought did not affect the CH₄ fluxes until 22 July, i.e. 21 d since the drought started. The reduction in CH₄ emissions continued until the end of August, when the drought was already over.

The annual carbon balance of the fen, estimated with the EC measurements of CO₂ and CH₄ fluxes, was close to zero in both years (Table 7).

Figure 5 (a) Daily CO₂ flux, (b) cumulative CO₂ flux, (c) respiration flux, (d) gross primary productivity, (e) CH₄ flux and (f) cumulative CH₄ flux measured with the eddy covariance technique during 2017 and 2018. Negative values denote an ecosystem sink and positive values denote a flux to the atmosphere.

3.3 Factors affecting CO₂ and CH₄ exchange

According to the LME models, the main environmental factors controlling the C exchange at the fen were WTL, GCC and vascular LAI for F_NEE1200 and F_GPP1200, i.e. the radiation-normalised F_NEE and F_GPP (Table 8). Increases in GCC and vascular LAI were associated with larger CO₂ sink, and the positive regression coefficient of WTL indicated that drier growing locations were larger CO₂ sinks. F_R increased with increasing soil temperature, GCC and vascular LAI, and reduced with increasing WTL. Higher CH₄ emissions were associated with wetter and warmer conditions, and higher GCC and graminoid LAI. Only the total vascular LAI was needed to explain the variation in CO₂ flux components, while for the CH₄ flux LAI had to be partitioned, because the graminoid and forb LAIs showed an opposite effect on the flux.

Soil temperature was highly correlated with air temperature ($R_{\mathrm{m}}^{\mathrm{2}}=\mathrm{0.88}$) and VPD ($R_{\mathrm{m}}^{\mathrm{2}}=\mathrm{0.78}$), GCC was highly correlated with air temperature ($R_{\mathrm{m}}^{\mathrm{2}}=\mathrm{0.76}$) and VPD ($R_{\mathrm{m}}^{\mathrm{2}}=\mathrm{0.74}$), while WTL was highly correlated with evergreen LAI ($R_{\mathrm{m}}^{\mathrm{2}}=\mathrm{0.75}$) (Table A5).

3.4 Plant-community-level CO₂ and CH₄ exchange

The growing season ecosystem respiration sums varied among the PCTs (Fig. 6a, Table A2): it was smallest in flarks and increased with the microtopographical altitude of the PCT. The growing season respiration sums differed significantly between the years, with the exception of SM (Fig. 6a), being higher in 2018 than in 2017 in the F (by 52 g C m⁻²), TT (by 79 g C m⁻²) and ST (by 109 g C m⁻²) communities.

Similarly to ER, the GPP sums increased gradually from the wettest F to the dry SM and ST communities (Fig. 6b, Table A2). The growing season GPP sum increased significantly from 2017 to 2018 in the F (by 33 g C m⁻²), SM (by 108 g C m⁻²) and ST (by 137 g C m⁻²) communities but not in TT.

Neither the PCTs nor the years differed significantly in their net ecosystem exchange of CO₂ (Fig. 6c). Additionally, the growing season NEE balances were similar in 2017 to 2018, because ER and GPP sums had a similar increase. This was also observed in the ecosystem-scale fluxes (Table 7). It appears that the flark communities (F and TT) shifted towards being a CO₂ source to the atmosphere, while the string communities (SM and ST) shifted towards being a sink of CO₂.

Table 7 Annual CO₂, CH₄ and carbon (CO₂+ CH₄) balances measured with the eddy covariance technique in 2017 and 2018.

Table 8 Standardised regression coefficients (± standard error) for the explanatory variables produced by linear mixed-effects models for chamber-based F_NEE1200, F_GPP1200, F_R and $F_{{\mathrm{CH}}_{\mathrm{4}}}$.

The ST communities had distinctly the lowest CH₄ emission. The growing season CH₄ balance of the F, TT and SM communities were similar in 2018, but in 2017 SM had a slightly larger CH₄ emission than F and TT (Fig. 6d, Table A2). CH₄ emissions differed significantly between the years only in the F communities, where the growing season balance was 2.2 g C m⁻² higher in 2018 than in 2017.

Figure 6 Mean growing season flux sums of (a) ecosystem respiration, (b) gross primary productivity, (c) net ecosystem exchange and (d) CH₄ flux of the flark (F), Trichophorum tussock (TT), string margin (SM) and string top (ST) plant communities. Negative NEE denotes an ecosystem CO₂ sink. Error bars represent the 95 % confidence intervals.

Seasonal variation in the ER, GPP and CH₄ fluxes was observed for all PCTs, with the maximum monthly exchange consistently taking place in July (Fig. 7). For NEE, however, the timing of the largest monthly exchange varied among the PCTs.

With the exception of October, the monthly mean respiration rates were higher in 2018 than in 2017, even though the difference for the TT, SM and ST communities was not significant (at the 5 % significance level) during the middle of the growing season (Fig. 7a, b, Table A3). The monthly mean GPP rates in the early growing season (May–June, and for the F community also in July) were also higher in 2018 than 2017 for all PCTs (Fig. 7c, d, Table A3).

The difference between the years was less clear in the monthly NEE. No significant differences were observed for F, whereas the NEE of TT and SM was higher in 2018 in May and April, respectively. In ST, a higher net emission could be detected in April and September 2018, and a higher net uptake in June and October 2018 (Fig. 7e, f, Table A3).

The differences in the monthly CH₄ emissions between 2017 and 2018 were not significant (Fig. 7g, h, Table A3). However, it seems that, while the flark plant communities F and TT had higher CH₄ emissions in 2018 than in 2017 in the first half of the growing season (May–July), the SM community had a lower CH₄ emission in 2018 than in 2017 in the latter part of the season (August–October).

Figure 7 Mean monthly flux sums of (a, b) ecosystem respiration, (c, d) gross primary productivity, (e, f) net ecosystem exchange and (g, h) CH₄ flux of the flark (F), Trichophorum tussock (TT), string margin (SM) and string top (ST) plant communities. Negative NEE denotes an ecosystem CO₂ sink. Error bars represent the 95 % confidence intervals.

4 Discussion

Using the C flux measurements with the EC technique in 2017 and 2018, we conclude that there were two phenomena that had a substantial effect on the annual C balance of the Kaamanen fen: the 2-week difference in the growing season start and the 1-month drought period in 2018. Additionally, with the chamber measurements, we could capture the variation in C exchange among different plant communities.

4.1 Spatial variation of C exchange among plant communities

Vegetation composition of plant communities is adapted to the prevailing moisture and nutrient conditions, which can vary greatly across the fen microtopography (Tables 1 and 6, Fig. 4). The separation of string tops from the water table is reflected in their shrub-dominated plant community composition, which in turn creates a soil of low pH and high C:N ratio. The ST PCT had the lowest pH of 4.6 (Table 6), i.e. a value approaching those found in bogs (pH < 4.2), while the three other PCTs had a pH (5.7–5.9) typical of mesotrophic fens (Wieder et al., 2006). Similarly, the C:N ratio was higher in the string than flark communities. Thus, the strings can be described as nutrient-poor bog-like islands within a mesotrophic fen (Aurela et al., 1998; Maanavilja et al., 2011). The string plant communities (SM, ST) with woody ericaceous shrubs had a relatively large LAI, while the flark communities dominated by graminoids (F, TT) had a smaller LAI (Fig. 4c, d). As expected, we found that GPP closely followed the quantity of photosynthesising plant material (Figs. 6 and 7, Table A4) (e.g. Alm et al., 1999; Bubier et al., 2003; Munir et al., 2014; Korrensalo et al., 2019). This was the case both for respiration and GPP: flarks with the highest WTL had the lowest exchange rates, which gradually increased along with increasing height above the water table (Figs. 4 and 6a, b, Table 7). This dependence was also confirmed by statistical modelling (Table 8), and the observation is in accordance with previous findings for northern fens with microtopography (Alm et al., 1997; Strack et al., 2006; Maanavilja et al., 2011). Maanavilja et al. (2011) estimated the growing season respiration sums of the Kaamanen fen to be approximately 50, 100, 250 and 225 g C m⁻² for F, TT, SM and ST, respectively, and the growing season GPP sums to be approximately 100, 140, 290 and 250 g C m⁻² for F, TT, SM and ST, respectively, which were in general lower than our estimates (Fig. 6a, b, Table A2).

In general, the differences in NEE balances between the PCTs were less clear than those in ER and GPP (Fig. 6c, Tables A2 and A4). NEE balances have earlier been found to vary substantially in mires with a hummock-hollow microtopography, with either hollow (Maanavilja et al., 2011; Schneider et al., 2012) or hummock (Bubier et al., 2003; Riutta et al., 2007) communities being larger net CO₂ sinks, or all communities being small CO₂ sinks (Strack et al., 2006). Maanavilja et al. (2011) estimated that in the Kaamanen fen the ST communities acted as the smallest growing season CO₂ sink (−10 g C m⁻²), while F, TT and SM were fairly similar sinks (−50, −40 and −40 g C m⁻², respectively). This differs from our finding that in 2017 the CO₂ balances ranged from −20 g C m⁻² (TT) to 64 g C m⁻² (SM), while in 2018 all PCTs were small CO₂ sources (10–22 g C m⁻²), which could be due to differing meteorological conditions. During the measurements of Maanavilja et al. (2011) in 2007, the growing season was not as rainy as in 2017 nor had drought events similar to 2018, and the monthly air temperatures and precipitation sums were close to the 30-year averages (Fig. 2).

In our data, CH₄ emissions were significantly higher from F, TT and SM than ST. In F and TT, the high emissions can be explained by the anoxic conditions that derive from the high WTL and are favourable for CH₄ production by archaea as well as by the graminoids of these communities that allow effective CH₄ transfer to the atmosphere (Ward et al., 2013). In SM, WTL is not as high, and a possible explanation for its high CH₄ emission is the higher plant biomass, which provides higher substrate availability for CH₄ production (Korrensalo et al., 2018). The low CH₄ emissions from the ST communities reflect the low WTL, inefficient CH₄ production and the lack of CH₄ transport routes through graminoid plants, which leads into higher CH₄ oxidation (Saarnio et al., 1997; Marushchak et al., 2016). Heikkinen et al. (2002) estimated that the monthly average CH₄ fluxes at this site from mid-June to the end of September 1995 were 2.3 g C m⁻² on flarks, 2.0 g C m⁻² on lawns (corresponding to TT) and 0.18 g C m⁻² on strings. These values are at the high end of the range of our estimates for F (mean 1.5, range 0.7–3.1 g C m⁻² month⁻¹) and TT (mean 1.7, range 0.8–3.3 g C m⁻² month⁻¹) for the same months. Heikkinen et al. (2002) did not report high CH₄ emissions for SM (mean 2.1, range 1.0–3.5 g C m⁻² month⁻¹ in the present study), most likely due to a different classification of string communities and lacking measurements from the margins. Our data for string tops (mean 0.11, range 0.08–0.14 g C m⁻² month⁻¹) are closer to the string emissions measured by Heikkinen et al. (2002). The vastly differing emissions from SM and ST suggest that the division between string margin and top is important for correctly estimating the ecosystem-scale CH₄ emissions as both PCTs cover a notable area of the fen (Table 1).

4.2 Temporal variation of C exchange

4.2.1 Annual variation

We estimated that the annual carbon balance of our fen was similar in 2017 and 2018 (Table 7), even though the meteorological conditions differed considerably between the years. The effect of an earlier onset of growing season in 2018 on the CO₂ balance was counterbalanced by the drought event later in that year. Aurela et al. (2004) have reported a mean annual CO₂ balance of −22 g C m⁻² (1997–2002) for the same fen, with variation between −4 and −53 g C m⁻². Our annual CO₂ balances were −8.5 and −5.6 g C m⁻² in 2017 and 2018, respectively, which are on the low side compared to Aurela et al. (2004). We assume that these differences mostly arise from the timing of snowmelt and spring temperatures. Hargreaves et al. (2001) estimated the annual CH₄ emissions of the fen to be 5.5 g C m⁻² (measurements conducted in 1995, 1997 and 1998), while our estimates of 7.3 and 6.2 g C m⁻² in 2017 and 2018, respectively, are slightly higher. These balances do not include the lateral aquatic transfer of dissolved organic C and particulate C through the fen ecosystem. Aurela et al. (2002) estimated, based on Sallantaus (1994) and Kortelainen et al. (1997), that the leaching of total organic carbon was 7.5 g C m⁻² yr⁻¹.

Considering the different PCTs, we found higher growing season sums of both ER and GPP in 2018 than in 2017 for all PCTs (Fig. 6a, b; Table A2). However, we did not observe significant changes in the growing season NEE (Fig. 6c, Table A2). CH₄ emissions differed significantly between the years in F only, with the growing season balance being 2.2 g C m⁻² higher in 2018 than in 2017 (Fig. 6d, Table A2).

Differences in ER, GPP and CH₄ emission sums between the years are in the same direction in all plant communities, i.e. higher values in 2018 especially during the early growing season months. From July onward, however, the changes in CH₄ exchange (Fig. 7, Table A3) were not as uniform. For most of these differences between the years, the spring weather emerges as a potential explanation.

4.2.2 Effects of spring timing

The higher ER and GPP rates during the early growing season of 2018, observed at both the ecosystem (Fig. 5c, d) and plant community (Fig. 7a–d, Table A3) level, were most likely due to the higher temperatures and earlier growing season start in April–June 2018 (Fig. 2a, Table 5), which advanced the activity of both plants and soil microbes (Jones, 2013). We observed the effect of earlier warm conditions in soil temperatures, GCC and LAI (Figs. 3e–h and 4c, d), which were in turn found to explain the variation in ER and GPP fluxes (Table 8). Warmer temperatures have been reported to increase soil respiration in boreal mires, while the impacts on GPP and NEE have varied more (Chivers et al., 2009; Ward et al., 2013). Early snowmelt and warm spring temperatures have been suggested to increase the annual net CO₂ uptake of northern mires (Aurela et al., 2004; Sagerfors et al., 2008), and our results support this observation: net CO₂ uptake increased during the early growing season at both the ecosystem (Fig. 5a) and plant community (Fig. 7e, f) scale. Of the PCTs, the string top communities showed a particularly large increase (Table A3).

The springtime increase of CH₄ emissions started during the thaw in May with a pulse of CH₄ that was stored below the snowpack (Fig. 5e), as has been commonly observed in subarctic mires (Friborg et al., 1997; Panikov and Dedysh, 2000; Hargreaves et al., 2001; Gažovič et al., 2010). The spring pulse in May, with a peak magnitude of 0.05 g C m⁻² d⁻¹, accounted for approximately 0.5 g C m⁻² of CH₄ emissions in both years, but it occurred a few weeks earlier in 2018 (Fig. 5e). This pulse occurred simultaneously with the soil temperature rise in the flarks (Fig. 3e, f). The pulses correspond to 6 %–7 % of the annual emissions, which is within the range of 3.5 %–11 % found in previous studies on subarctic fens (Friborg et al., 1997; Hargreaves et al., 2001) but larger than the proportion of < 3 % that Rinne et al. (2007) observed at a boreal fen. Even though the springs differed and the flark plant communities F and TT seemed to have higher CH₄ emissions in May–July in 2018 than in 2017 (Table A3), the differences in the monthly CH₄ balances between 2017 and 2018 remained small (Fig. 7g, h, Table A3). This suggests that, unlike the CO₂ exchange responses, the spring weather variation mostly affects the timing but not the magnitude of CH₄ exchange in northern mires.

4.2.3 Effects of summer drought

Our results show that the C exchange of the fen was significantly affected by the drought that took place in July 2018 (Fig. 3b, d). The drought was observed as a water level drawdown by 5–20 cm, higher-than-average temperatures and an elevated VPD. These anomalies likely caused drought stress in plants (Alm et al., 1999). The EC measurements suggest that this event decreased net CO₂ uptake by decreasing GPP rapidly (Fig. 5a, d), most likely because the plants regulated their stomatal openings and gas exchange as water availability decreased.

The drought impact was less obvious in chamber measurements, i.e. at the PCT level, but a bit surprisingly it appears that CO₂ uptake decreased in the wet flark communities F and TT and increased in the dry string communities SM and ST during the drought. The CH₄ emissions from SM also seemed to decrease following the reduced methanogenesis in the newly formed oxic layer (Deppe et al., 2010), but this took place after the drought in August–October (Fig. 7g, h, Table A3). In spring, a lack of precipitation would not affect water availability as much as in midsummer, as the meltwater from snow maintains a high water table level. While it is known that respiration depends on WTL (Alm et al., 1997; Christensen et al., 1998; Bubier et al., 2003), it appears that the effect of WTL on temporal ER variation is specific to a mire and a microform type. For instance, Strack et al. (2006) reported a significant increase in ER during drought events on a poor fen with microtopography, while Deppe et al. (2010) found no significant effect of WTL fluctuations on CO₂ exchange on an ombrotrophic bog and alpine wetland, and Aurela et al. (2007) found a nonmonotonic link between ER fluxes and WTL on a sedge fen.

The drought decreased the ecosystem-scale CH₄ emissions temporarily (Fig. 5e). The drought-induced decrease is likely due to the reduced volume of anoxic peat, while the oxic zone increased correspondingly, thus reducing CH₄ production and increasing CH₄ oxidation (Strack and Waddington, 2008; White et al., 2008; Deppe et al., 2010). However, we also observed differences among the PCT-specific responses; the SM communities reacted most to the drought, with smaller August emissions in 2018 than in 2017 (Fig. 7g, h, Table A3). In F, the growing season CH₄ emissions were higher in 2018 than in 2017, most likely due to the higher soil temperatures in 2018, as CH₄ emissions increase with increasing peat temperature until an optimal methanogenesis temperature within 20–30 ^∘C (e.g. Dunfield et al., 1993; Laine et al., 2007b). In the other PCTs, no systematically higher emissions were observed in 2018, probably since the temperature-induced enhancement was offset by the lower WTL during the drought. At the ecosystem level, the difference in the annual CH₄ balance remained minor (7.3 vs. 6.2 g C m⁻²). This indicates that the microtopography typical to northern mires may be a factor that increases the stability of their CH₄ emissions with respect to abrupt environmental changes such as drought and heatwaves.

In an earlier study, the timing of snowmelt was shown to be a key parameter controlling the annual CO₂ balance of the Kaamanen fen (Aurela et al., 2004). In our study, however, the higher C sequestration during the longer growing season in 2018 was offset by the drought in July. This indicates that extreme weather events can substantially affect the annual C balance of northern fens through affecting their CO₂ and CH₄ exchange. However, the drought, which covered all of north-western Europe, did not affect the CO₂ and CH₄ exchange at Kaamanen as much as it did at more southern mires during the same period (Rinne et al., 2020). The water availability of fens that have more or less continuous flow through them have greater resilience to water level drawdown during droughts than ombrotrophic bogs. Also, here the fen microtopography seems to increase the resistance to C loss at the ecosystem scale, as each of the plant communities are adapted to different environmental conditions. This could be seen with the string plant communities that acted as CO₂ sinks during the drought, thus decreasing the overall drought impact on the fen, which most likely is a direct consequence of string top species being already adapted to drier environments.

Heatwaves are predicted to become more frequent in the subarctic region as the climate warms (Masson-Delmotte et al., 2018). However, the impact of heatwaves on the C exchange of northern mires strongly depends on local soil moisture conditions. While drought leads to diminished C sequestration, warming accompanied by sufficient precipitation is likely to support the long-term peat accumulation in the subarctic, non-permafrost mires (Loisel et al., 2020). On the other hand, the vegetation composition and biomass production on these fens are susceptible to lowering water table level (Mäkiranta et al., 2018). Therefore, the functioning of subarctic fens may undergo substantial changes if the water balance changes concurrently with the warming climate.

5 Conclusions

We studied CO₂ and CH₄ exchange of a subarctic fen and found both sensitivity and inherent resilience in their response to meteorological variation. Even though meteorological and environmental conditions differed in many ways between the two measurement years, our EC data showed that the annual C balance of the fen did not differ markedly between the years. While the relatively early onset of the growing season in 2018 strengthened the CO₂ sink, this gain was counterbalanced by a later drought period. Variations in water table level, soil temperature and vegetation characteristics (leaf area and greenness) explained the majority of the variation in the ecosystem-level C exchange. These environmental factors also varied among the PCTs, which was reflected in their widely differing CO₂ and CH₄ fluxes. The flark and Trichophorum tussock communities had lower ER and GPP than string margins and tops. Even though the string margin and top PCTs were similar in terms of CO₂ exchange, CH₄ emissions from string margins were notably larger than those from string tops. In 2017, they even clearly exceeded the emissions from flarks and Trichophorum tussocks. The mean ER and GPP fluxes of all PCTs were higher during the warmer 2018 growing season than in 2017, while the changes in NEE and CH₄ fluxes were lesser.

The characteristic microtopography present at the Kaamanen fen generates a wide range of environmental conditions. This sustains a diversity of adapted plant and microbe communities, making the fen resilient to C loss during extreme meteorological events. However, if drought events become more common, the long-term impacts on ecosystem functioning may be more drastic than what was observed in our study.

Appendix A: Abbreviations, statistical tests and analysis

Table A1 List of abbreviations.

Table A2 PCT-specific growing season sums of ER, GPP, NEE (the difference between ER and GPP) and CH₄ flux. Statistically significant differences between 2017 and 2018 are indicated with an asterisk (Z test, p< 0.05).

Table A3 PCT-specific monthly ER, GPP, NEE and CH₄ flux sums. Statistically significant differences between 2017 and 2018 are indicated with an asterisk (Z test, p< 0.05).

Table A4 Comparison of the growing season ER, GPP, NEE and CH₄ flux sums between the plant community types. Statistically significant differences (Z test, p< 0.05) are indicated with an asterisk.

Table A5 Squared correlation coefficients ($R_{\mathrm{m}}^{\mathrm{2}}$) between the variables included in the LME model analysis.