Winter fluxes determine the annual carbon balance of an unmanaged subarctic drained peatland

Peatlands are critical components of the global carbon (C) cycle, storing large amounts of soil organic carbon (SOC). However, drainage substantially alters their carbon exchange and hydrological functioning, often converting them into net carbon dioxide (CO₂) sources. This study presents the first year-round, ecosystem-scale Eddy Covariance (EC) assessment of CO₂ dynamics from an unmanaged drained peatland in western Iceland, originally drained in the early 1960s. Two years of continuous EC measurements were collected alongside high-resolution environmental data, including solar radiation, air and soil temperatures, soil water content, and groundwater level. Several multispectral drone flights were also conducted during the study period, which provided seasonal NDVI-based estimates of canopy greenness. The two study years differed markedly in annual weather during the growing season (GS), with 2023 GS being unusually warm and dry, while 2024 GS was cold and wet. Despite these contrasts, annual net ecosystem exchange (NEE) remained similar between the 2 years. The annual NEE was dominated by non-growing-season (NGS) respiration, which highlighted the necessity for year-round measurements. Overall, the site remained a persistent CO₂ source, emitting 4.1–4.4 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ nearly 60 years after drainage. Temperature exerted the strongest control on ecosystem respiration (R_eco), while gross primary production (GPP) responded primarily to seasonal irradiance and NDVI. A compensatory mechanism was observed during the warm year (2023) at this relatively cool site, where warming-induced increases in R_eco were offset by an enhanced GPP, resulting in a relatively stable annual NEE despite meteorological contrasts. Soil moisture and vapor pressure deficit played only minor roles under these cool and moist conditions. These findings highlight the need for continued monitoring of unmanaged drained peatlands to better quantify their contribution to regional greenhouse gas budgets.

1 Introduction

Peatlands cover only 3 % (4.23 million km²) of the Earth's land surface, yet they represent a massive global carbon (C) reservoir, storing approximately 600 Gt C (Yu et al., 2010; Xu et al., 2018; Loisel et al., 2021). This C storage is remarkably large, as peatlands store a similar amount of C as global forest soils and litter combined, despite forests covering 31 % of the Earth's land surface (FAO, 2020; Loisel et al., 2021). Functionally, these ecosystems remain a significant net C sink of approximately 0.14 Gt yr⁻¹, contributing 3 %–10 % of the total global terrestrial C sink (IPCC, 2007; Loisel et al., 2021). Peatlands function as C reservoirs because waterlogged conditions suppress decomposition, allowing organic matter (peat) to accumulate over thousands of years, often since the Last Glacial Maximum (Gorham, 1991; Yu et al., 2010). High-latitude peatlands, in particular, play an important role in the global C cycle, as their large soil C stocks may be highly sensitive to climatic warming and changes in hydrology (Loisel et al., 2021).

Drainage of natural peatlands for agriculture, forestry, or human infrastructure development is a major disturbance factor; it is estimated that the drainage has degraded 15.7 % of global peatlands, often transforming them from C sinks to persistent C sources (Joosten and Clarke, 2002; Leifeld and Menichetti, 2018; Leifeld et al., 2019). Drainage typically accelerates peat decomposition by lowering the water table and promoting an aerobic microbial oxidation of organic matter (Loisel et al., 2021). However, drainage can also alter vegetation traits (such as species composition and diversity), potentially increasing photosynthetic rates, biomass, and C inputs (Jauhiainen et al., 2019). In some cases, particularly where forests are established on drained peatlands, the increased C sequestration in biomass can temporarily offset soil carbon losses, turning the site into a net C sink (Minkkinen and Laine, 1998; Lohila et al., 2011; Bjarnadottir et al., 2021). Ultimately, the balance between drainage-induced changes in carbon inputs and outputs strongly depends on climate and management practices (Maljanen et al., 2010), highlighting the critical importance of monitoring the net carbon balance of drained wetlands in different settings.

Drained organic soils are recognized as significant sources of atmospheric greenhouse gas (GHG) emissions within national inventories under international treaties. In the 2013 IPCC Wetland Supplement (IPCC, 2014), the available studies on the GHG balances of drained peatlands at the time were used to establish default Tier 1 emission factors (EFs) to be used globally for countries lacking accurate national data. For the boreal zone, the Tier 1 EFs for CO₂ are 0.25–0.37, 0.93, 5.7 and 7.9 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for drained nutrient-poor forest/fallow lands, drained nutrient-rich forest lands, drained grasslands, and drained croplands, respectively (where drainage occurred < 50 years ago). In recent years, there have been several new studies published on the CO₂ and other GHGs balances of some drained land classes (e.g. Evans et al., 2021; Jauhiainen et al., 2023; van Giersbergen et al., 2025). Jauhiainen et al. (2023) reviewed northern European studies on drained forest peat soils and published more accurate EFs, while similar initiatives are still missing for croplands and grasslands of N-Europe; especially there are few studies existing for unmanaged drained peatlands at higher latitudes (Guðmundsson et al., 2024). In a study in Denmark, annual ecosystem respiration (R_eco) was found to increase from 4.2 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in an undrained bog to 13.3 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in a nearby drained permanent grassland (cattle grazing) (Kandel et al., 2018). The difference of 9.1 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ can, however, not be interpreted as a local EF, as neither net ecosystem exchange (NEE) nor annual CO₂ uptake were measured and therefore net carbon balance was unknown.

A major limitation of many existing peatland carbon balance studies is the incomplete representation of non-growing-season (NGS) fluxes, which are often excluded or poorly constrained due to limited winter measurements (Alm et al., 1999; Aurela et al., 2004). At high latitudes, winter periods are long and can contribute significantly to annual CO₂ emissions, even when short-term fluxes are relatively low (Aurela et al., 2002). This is particularly relevant in maritime climates, where relatively mild winter conditions can sustain continuous R_eco throughout much of the NGS (Bjornsson et al., 2007). Many studies rely on growing-season (GS) chamber measurements, requiring model-based extrapolation to estimate annual budgets (e.g. Ojanen et al., 2010; Koskinen et al., 2016). In contrast, the eddy covariance (EC) method quantifies ecosystem-scale CO₂ exchange continuously, allowing for a direct assessment of both GS and NGS contributions to the annual carbon balance (Aurela et al., 2002).

Drainage of peatlands in Iceland began in the early 20th century and was most intensive between the 1940s to 1970s, with little new drainage occurring in recent decades (Hallsdóttir et al., 2012; Guðmundsson and Óskarsson, 2014). Fully drained organic soils cover about 4196 km², consisting mostly of unmanaged grasslands (Hallsdóttir et al., 2012; Guðmundsson and Óskarsson, 2014). These areas represent the largest single source of CO₂ emissions in Iceland, estimated at ∼ 8.5 million t CO₂-C annually (∼ 65 % of reported national GHG emissions; Umhverfisstofnun, 2025). Iceland currently applies the IPCC (2014) default EFs for all drained land categories, even though more than half of these areas were drained over 50 years ago and Icelandic peatlands differ significantly from those at comparable latitudes (Arnalds et al., 2016b). Specifically, Icelandic peatlands are unusually mineral-rich due to volcanic and aeolian inputs, which increase nutrient availability and may influence carbon chemistry and turnover (Dagsson-Waldhauserova et al., 2014; Arnalds, 2015; Möckel et al., 2021, 2023).

Limited data exists on the actual annual CO₂ emissions from Icelandic drained peatlands, both unmanaged and croplands. Using peat-soil inventories of eight paired undrained and drained unmanaged 15–50-year-old site-pairs in southern Iceland, Gunnarsdóttir (2017) found an average annual loss of 1.7 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, ranging between 0.7 and 3.1 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. These estimates are much lower than the 5.7 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ IPCC (2014) default EF for boreal grasslands. In another study, Ólafsdóttir (2015) used chamber techniques to measure surface CO₂ fluxes from an uncultivated drained peatland in western Iceland, and after interpolating the annual flux found a NEE of 3.8 ± 0.65 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, which is also below the IPCC (2014) default EF.

Figure 1 Location of the LD site (unmanaged drained peatland) in W-Iceland. The red star marks the eddy covariance (EC) tower, while dark blue circles indicate buffer zones with radius of 25, 50, and 100 m around the tower. Light blue arrows illustrate water discharge pathways from the field toward the ocean, and white arrows indicate the prevailing wind direction and magnitude toward the tower.

Iceland possesses a cool temperate maritime climate at elevations below 400 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Bjornsson et al., 2007), with a growing season of approximately four months (Einarsson, 1984). Being an island in the middle of the North Atlantic Ocean, its winters are relatively mild compared to the more continental climates found at similar latitudes in Scandinavia and North America (Bjornsson et al., 2007; Ruosteenoja and Jylhä, 2022). In fact, the current relatively warm winter climate in lowland Iceland represents similar conditions as are predicted to occur at the same latitudes in Scandinavia during the latter part of this century if GHG emissions are not reduced (Hanssen-Bauer et al., 2017; Ruosteenoja and Jylhä, 2022).

The Umhverfisstofnun (2025) national inventory report emphasized that continuous flux measurements from unmanaged, drained sites in Iceland are still lacking. To address this knowledge gap, this study quantified the annual NEE of an unmanaged drained peatland in western Iceland over a two-year period using EC measurements. The study aimed to determine: (i) how close the unmanaged site was to the IPCC's (2014) default EF for drained grasslands, (ii) how strongly the NEE values and its component fluxes were influenced by different environmental conditions and (iii) how NEE and its component fluxes differed between the GS and the NGS, in the cool, maritime climate of Iceland.

2 Material and methods

2.1 Study site

The study was carried out at a 12 ha unmanaged drained grassland site at the Lækur Farm in Western Iceland (64.39° N, 21.89° W), hereafter termed the LD site. The site was artificially drained in 1961 through a network of ca. 2 m deep ditches spaced roughly 52 m apart and by putting in subterranean channels at ca. 60 cm depth between the ditches. Despite the drainage effort, the land has never been cultivated or fertilized, allowing natural vegetation succession under altered hydrological conditions. This has resulted in a dense, highly productive grass cover. There was some livestock grazing at the site prior to 1997, but since then it has remained protected. The terrain is predominantly flat, with subtle hummocks and the site drains southeast and southwest toward the ocean, which is > 120 m away from the lowest edge of the LD site (Fig. 1). The peat soil of the site is classified as Histosols, with a bulk density of 0.2 g cm⁻³, a pH of 4.6 and a $\mathrm{C}/\mathrm{N}$ ratio of 20.36 in the 0–50 cm layer.

Table 1 Summary of site characteristics and environmental parameters for the LD site (unmanaged drained peatland) in W-Iceland during 2023–2024, including climate, topography, soil, hydrology, and vegetation data.

^a Data from Hafnarfjall weather station (ID 31674; ∼ 10 km away from the study site) (Icelandic Meteorological office, 2025). ^b Data from Neðra-Skarð weather station (ID 97; ∼ 5 km away from the study site) (Icelandic Meteorological office, 2025).

The climate is cool maritime, with a mean annual temperature (MAT) of +4.8 °C and mean annual precipitation (MAP) of 1011 mm, based on 2023 and 2024 data from the nearest weather station to the site (Icelandic Meteorological office, 2025). Seasonal minimum and maximum temperatures ranged from −14 °C in winter to 20 °C in summer during the same period. Prevailing winds originate mainly from the east and northeast, averaging 4 m s⁻¹ (Table 1). The uniform microtopography of the site likely supports a relatively consistent microclimate.

The current vegetation at the site is classified as a Boreo-subalpine Agrostis grassland habitat type (Ottósson et al., 2017), with common bentgrass (Agrostis capillaris) and Arctic fescue (Festuca richardsonii) as dominating species. The peak season's plant canopy height is 20–50 cm and maximum annual aboveground biomass is ∼ 720 g m⁻² (Table 1).

2.2 Instrumentation and data collection

2.2.1 Eddy Covariance system

An EC tower was established in early January 2023 at the drained (LD) site, positioned based on predominant wind directions to optimize the flux footprint and minimize edge effects (Fig. 1). The instrumentation consisted of an open-path infrared gas analyser (LI-7500DS, LI-COR Biosciences, Lincoln, NE, USA) for real-time CO₂ and H₂O measurements, integrated with a 3D sonic anemometer (WindMaster Pro, Gill Instruments, Lymington, UK). The system recorded three-dimensional wind velocity and gas concentrations at a sampling frequency of 10 Hz. In addition to CO₂ exchange, the EC system measured latent heat flux (LE) and sensible heat flux (H), providing insight into energy exchange processes and friction velocity (u_∗). From January 2023 until May 2024, flux measurements were conducted at a height of 2.15 m above the canopy, after which the sensor height was increased to 3.5 m.

Due to the absence of grid power, the setup was powered using a hybrid power solution consisting of 200 W solar panels and a methanol fuel-cell generator (EFOY Pro 900, SFC Energy AG, Brunnthal, Germany). This configuration ensured continuous year-round operation, particularly during winter months when solar radiation is insufficient at high latitudes. The fuel cell was strategically positioned northwest of the EC system, outside the dominant wind directions, to avoid interference with flux measurements; footprint analysis and CO₂ concentration data confirmed negligible influence from this sector (further details are provided in Appendix A).

Table 2 The instrumentation and technical specifications at the LD site (unmanaged drained peatland) in W-Iceland.

2.2.2 Biomet and auxiliary sensors

A micro-meteorological tower was installed near the EC system. This tower was equipped with sensors to monitor air temperature (T_a), humidity (RH), precipitation (Rain), and solar radiation (PAR). Soil temperature (T_s) and soil moisture (SWC; with a factory calibration for organic soil) were measured at fixed depth of 10 cm (point measurement) below the surface (Table 2).

Groundwater level (WL) was recorded hourly close to the EC tower using a pressure sensor with barometric compensation (Onset, HOBO MX2001 water level logger, Bourne, MA, USA). The WL well was manually augured to a depth of 330 cm to capture the full range of seasonal water table fluctuations. Additionally, five wells were augured at different locations within the site where WL was manually measured. The WL data were post-processed to correct for sensor drift and noise from freezing conditions. WL is reported as the vertical distance from the peat surface, with negative values indicating that the water table was below ground level.

2.2.3 Site characterization and surface mapping

To characterize spatial heterogeneity and surface features at the study site, a combination of in-situ measurements and drone-based mapping was used. Peat depth was assessed using a manual peat depth probe (Eijkelkamp Soil & Water, Giesbeek, Netherlands), inserted vertically until the mineral substrate was encountered. A total of 46 measurements were collected across a systematic grid. A peat depth map was generated from the point measurements using Inverse Distance Weighting (IDW) interpolation in QGIS (v3.28.9-Firenze, QGIS Development Team). IDW was selected over geostatistical methods because the peat thickness showed significant local variability without a clear regional trend, and the sample size was better suited for a deterministic interpolation approach.

A high-resolution Digital Surface Model (DSM) was generated from drone imagery acquired during low-altitude flights. Imagery was captured using a 10-band MicaSense RedEdge-MX Dual Camera Imaging System (MicaSense Inc., Seattle, WA, USA) mounted on a DJI Matrice 300 RTK drone (SZ DJI Technology Co., Ltd., Shenzhen, China). To ensure temporal comparability, all imagery was radiometrically calibrated using a calibrated reflectance panel and a downwelling light sensor (DLS 2) for each flight. The data were processed using structure-from-motion (SfM) photogrammetry in Agisoft Metashape (v2.1.1, Agisoft LLC, St. Petersburg, Russia). This technique is well-suited for high-precision surface reconstruction in ecological studies (Westoby et al., 2012; Dandois and Ellis, 2013).

Seasonal changes in the site's Normalized Difference Vegetation Index (NDVI) were derived from multispectral imagery acquired during 14 drone flights (5 in 2023; 9 in 2024). NDVI was calculated using the standard spectral reflectance ratio:

$$\begin{array}{}\text{(1)}& \text{NDVI}={\displaystyle \frac{R_{\text{NIR}}-R_{\text{Red}}}{R_{\text{NIR}}+R_{\text{Red}}}},\end{array}$$

Where R_NIR and R_Red represent surface reflectance in the near-infrared and red bands, respectively (Rouse et al., 1974). Weekly NDVI values were estimated by fitting a smoothing spline to discrete flight data to reflect the seasonal development of canopy greenness.

2.3 EC data processing

EC data were processed using EddyPro software (v7.0.9, LI-COR Biosciences, Lincoln, NE, USA), following standard protocols for open-path sensor configurations. The processing workflow included angle-of-attack correction, double coordinate rotation, and time lag compensation. Spectral losses were compensated following Moncrieff et al. (1997). To adjust for fluctuations in air temperature and water vapor, Webb–Pearman–Leuning (WPL) density correction was applied (Webb et al., 1980). To account for instrument-related density effects, the Burba et al. (2008) correction was applied to the entire 2023–2025 dataset. Owing to the high-latitude, cool-climate setting of the Icelandic site, a temperature gradient between the ambient air and the LI-7500DS sensor head was expected to persist for much of the year. Applying the correction consistently across all seasons ensured a uniform processing workflow and avoided introducing artificial discontinuities at seasonal transitions. The correction magnitude scaled with the observed air-instrument temperature gradient and naturally became negligible during warmer periods.

To account for CO₂ accumulation below the measurement height, a storage flux correction was applied following the approach of Aubinet et al. (2001). Although the EC system height was increased from 2.15 to 3.5 m during the study period, data comparability was maintained through consistent storage corrections and the site's high degree of topographical and vegetative homogeneity.

The EC system data was integrated with the micro-meteorological sensors data (Table 2) and automatically time-aligned via the EddyPro software, enabling the computation of half-hourly averages for NEE, H, and LE. This integration also provided the metadata and diagnostic quality flags required for subsequent quality assurance and data analysis.

Following initial processing, a multi-step quality control (QC) procedure was implemented. Outliers were removed through a combination of statistical screening and visual inspection. The dataset was filtered using the turbulence and stationarity classification scheme of Mauder and Foken (2011); only data with high (0) or moderate (1) quality flags were retained. Data periods affected by instrument malfunction or signal loss were also excluded. u_∗ thresholds were determined using the moving point test implemented in the REddyProc R package (Wutzler et al., 2018). Thresholds of 0.26 and 0.29 m s⁻¹ were applied to the 2023 and 2024 datasets, respectively, to minimize underestimation of nocturnal fluxes (Papale et al., 2006). An additional filtering step was applied specifically to CO₂ fluxes. Data collected during precipitation events, dew formation, or low vapor pressure deficit (VPD) were excluded, as these conditions interfere with open-path analyser performance and can yield biologically unrealistic values (Wohlfahrt et al., 2005; Aubinet et al., 2012).

Energy balance closure (EBC) was evaluated by comparing the sum of the turbulent energy fluxes to the available energy (Foken, 2008). To account for energy temporarily stored within the ecosystem, storage terms were included in the energy balance equation:

$$\begin{array}{}\text{(2)}& H+\text{LE}=R_{\mathrm{n}}-G-S,\end{array}$$

Where H is the sensible heat flux (W m⁻¹), LE is the latent heat flux (W m⁻²), R_n is the net radiation (W m⁻²), G is the ground heat flux (W m⁻²), and S represents the sum of heat storage terms, including storage in the air column (S_a) and the soil layer above the heat flux plates (S_g).

The spatial representativeness of the EC fluxes was assessed using the Flux Footprint Prediction (FFP) model of Kljun et al. (2015) using Tovi software (v2.9.1, LI-COR Biosciences, Lincoln, NE, USA). Footprints were computed at 30 min resolution using measured u_∗, Obukhov length (L), wind direction, and lateral wind variability (σ_v), with constant surface parameters of roughness length (z₀ = 0.035 m) and displacement height (d = 0.23 m).

2.4 Data analysis and statistical approaches

2.4.1 CO₂ Flux gap-filling and partitioning

Meteorological drivers, including T_a, T_s, VPD, and PAR, were gap-filled prior to flux processing. Missing values were substituted using data from auxiliary on-site sensors and the nearest Icelandic Meteorological Office (IMO) weather station. This pre-processing follows standard FLUXNET and REddyProc procedures, which require continuous meteorological inputs for CO₂ gap-filling and flux partitioning (Wutzler et al., 2018; Pastorello et al., 2020).

The GS was defined using a daily mean T_a threshold of 5 °C, following the standard Nordic climatological approach (Carter, 1998; Førland et al., 2004). Since short warm spells in winter are common in Iceland's maritime climate and to avoid falsely identifying early or mid-winter warm events as genuine growing-season onset, the growing season onset and end was only confirmed when temperatures remained above or below 5 °C for at least 10 consecutive days.

Initial CO₂ data coverage across the study period was 78 %, mainly due to power outages. After applying quality filters (e.g., removing low-turbulence periods and sensor-related artifacts gaps lasting more than 12 h), the remaining data coverage was 65 %. To produce a continuous record of CO₂ fluxes, a two-step gap-filling approach combining Marginal Distribution Sampling (MDS) (Reichstein et al., 2005) and Extreme Gradient Boosting (XGBoost) (Chen and Guestrin, 2016) was applied.

First, MDS was performed using the REddyProc package which fills missing flux values by sampling observed data during periods with similar environmental conditions within flexible temporal windows. Environmental drivers including PAR, T_a, and VPD were used to account for seasonal and diurnal variability. Short-term gaps (≤ 24 h) were filled using MDS, while MDS estimates associated with high uncertainty (standard deviation > 4 $\mathrm{\mu }\mathrm{mol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) were excluded.

Long-term gaps and unreliable MDS outputs were subsequently re-estimated using XGBoost, a machine-learning algorithm well-suited for capturing non-linear relationships between carbon fluxes and environmental drivers in Python (v3.12.3; Python Software Foundation). To account for seasonal shifts in ecosystem response, separate models were trained for the GS and NGS using only high-quality NEE observations. Predictor variables included PAR, T_a, T_s, VPD, and time-of-year indicators to capture phenological trends.

The model was evaluated using 5-fold cross-validation on 14 610 cleaned observations. In each iteration, the model was trained on 80 % of the data and validated on the remaining 20 %, ensuring that all data points were independently used for validation. Model performance was quantified using the coefficient of determination (R²) and root mean square error (RMSE). To assess potential overfitting, training and validation performance were compared. Furthermore, SHAP (SHapley Additive exPlanations; Lundberg and Lee, 2017) values were calculated to provide a physically meaningful interpretation of variable importance and to quantify the contribution of each environmental driver to the predicted NEE. This approach effectively addressed missing or low-quality data and enabled reliable reconstruction of NEE over the entire study period.

The NEE was partitioned into gross primary production (GPP) and R_eco using a physics-guided hybrid approach. Ecosystem respiration was first parameterized from nighttime NEE with the Lloyd–Taylor temperature response function using T_s as the primary environmental driver (Lloyd and Taylor, 1994):

$$\begin{array}{}\text{(3)}& R\left(T\right)=R_{\text{ref}}\exp \left[E_{\mathrm{0}}\left({\displaystyle \frac{\mathrm{1}}{T_{\text{ref}}-T_{\mathrm{0}}}}-{\displaystyle \frac{\mathrm{1}}{T-T_{\mathrm{0}}}}\right)\right],\end{array}$$

where R_ref is R_eco at a reference temperature (T_ref = 283.15 K), E₀ is the activation energy-type parameter, and T₀ is a fitted constant (227.13 K). To account for freeze–thaw state dependence common in Icelandic peatlands, the function was fitted across two thermal regimes (cold: T_s≤ 2 °C; warm: T_s > 2 °C). Non-temperature-driven variations (calculated as the residual between observed nighttime NEE and the Lloyd–Taylor R_eco) were modelled using a constrained XGBoost algorithm. This residual model was trained on environmental drivers (T_a, T_s, VPD, SWC), diel and seasonal harmonics, and cumulative thermal indicators such as thawing and freezing degree days. To maintain biological consistency, monotonicity constraints were enforced, ensuring that predicted nighttime R_eco did not decrease with increasing temperatures. Furthermore, a temperature-based gating function was applied to prioritize the Lloyd–Taylor baseline during frozen conditions while allowing XGBoost to capture complex biological variance during warmer, more active periods. The model was evaluated using 5-fold cross-validation, and potential overfitting was assessed by comparing predictive performance metrics (R² and RMSE) between the training and validation sets. SHAP values were also calculated to assess the relative importance of model input variables.

Daytime GPP was subsequently calculated by subtracting the measured NEE from the modelled R_eco (GPP = R_eco−NEE), where negative values indicate CO₂ uptake (net-photosynthesis). During the prolonged twilight hours typical of 64° N, light sensors can become unreliable; therefore, when PAR fell below 150 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, GPP values were smoothed toward a standard light-response curve (Lasslop et al., 2010). To prevent sensor noise from biasing the carbon budget, GPP was strictly set to zero during astronomical nights (solar elevation ≤ −12°). Finally, uncertainty in GPP and R_eco estimates was quantified using bootstrap resampling of the entire processing chain (N=30; meaning that the entire workflow including Lloyd–Taylor fitting, XGBoost training, and flux partitioning was repeated 30 times), yielding time-dependent standard deviations.

Carbon flux dynamics were analysed at multiple temporal scales, including diurnal, seasonal, and annual cycles. This multi-scale approach enabled the assessment of both short-term variability and longer-term ecosystem behaviour. Annual carbon budgets were calculated for the full calendar years of 2023 and 2024, expressed in $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$. A mean annual value was derived from both years to represent the long-term site behaviour.

2.4.2 Statistical relationships

To examine how CO₂ fluxes responded to environmental conditions, analyses were restricted to periods with observed (non-gapfilled) NEE and the corresponding partitioned GPP and R_eco values. Environmental drivers (accumulated weekly PAR (PAR_acc), VPD, T_a, T_s, SWC, and WL) were likewise limited to their measured values. Weekly averages of fluxes and weekly summaries of drivers were then computed for both study years and used to evaluate relationships between environmental conditions and GPP or R_eco. In addition to individual drivers, a combined radiation–vegetation metric (PAR_acc × NDVI) was calculated, with weekly NDVI representing canopy greenness and commonly used parameter as a proxy for the fraction of absorbed PAR (fPAR) in light-use-efficiency frameworks that are typically linearly related to GPP sums (Monteith, 1972; Sellers, 1985). The PAR_acc × NDVI product is therefore an approximation to total intercepted PAR.

Linear regressions were subsequently applied in Python to quantify how strongly each environmental variable explained weekly variation in GPP and R_eco within each study year, enabling comparison of the relative influence of individual drivers and interannual differences in response strength. Model performance was evaluated using R² and associated p-values to assess statistical significance.

For the relationship between R_eco and temperature, a non-linear Q₁₀ temperature-response function following Raich and Schlesinger (1992) was applied to capture the characteristic exponential increase in metabolic activity with temperature, which linear models typically fail to represent across a broad thermal range. The function is expressed as:

$$\begin{array}{}\text{(4)}& R\left(T\right)=R_{\text{ref}}\times Q_{\mathrm{10}}^{(T-T_{\text{ref}})/\mathrm{10}},\end{array}$$

Where R(T) is ecosystem respiration at temperature T, R_ref is respiration at a reference temperature (T_ref = 10 °C), and Q₁₀ represents the factor by which respiration increases for a 10 °C rise in temperature.

Furthermore, Pearson correlation coefficients (Pearson, 1896) were calculated across the full study period to assess the overall strength and direction of linear associations between weekly flux components (NEE, GPP, and R_eco) and environmental drivers, integrating seasonal progression and interannual variability. Pearson's r was selected to characterize broad-scale relationships rather than year-specific responses.

Figure 2 Site characteristics of the LD site (unmanaged drained peatland) in W-Iceland: (a) Digital Surface Model (DSM) of the site showing the 80 % flux footprint boundaries for 2023 (instrument height of 2.15 m above mean canopy height) and 2024 (instrument height of 3.5 m); (b) peat depth map derived from spatial interpolation, overlaid with the 46 manual measurement points; (c) seasonal NDVI patterns across the 2023–2024 study period.

3 Results

3.1 Surface and subsurface characteristics

The high-resolution DSM with 7 cm ground sampling distance shows that the site is predominantly flat, with most of the area situated around 7 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Elevation ranged from 1 to 11 m, with the lowest point along the eastern boundary and the highest in the northwest, though these extremes fell largely outside the main flux footprint (Fig. 2a).

In 2023, when the measurement height was 2.15 m, the 80 % footprint averaged 50–70 m around the tower, which was smaller than expected due to the rough, hummocky surface (Fig. 2a). After the measurement height was increased to 3.5 m in 2024, the 80 % EC flux footprint extended on average 100–120 m around the tower (Fig. 2a). The 80 % footprint remained within the main study area for both years and covered relatively uniform terrain with less than 2 m elevation variation and similar vegetation composition. The site was characterized by predominantly easterly winds and high u_∗ (mostly above 0.2 m s⁻¹), indicating that mechanical turbulence dominated and atmospheric conditions were mostly neutral.

Within the footprint, subtle microtopographic features were present, including hummocks and drainage ditches with their associated spoil ridges. Surface classification indicated that ∼ 3.5 % of the footprint consisted of ditches, 9.5 % of excavated peat soil, and the remaining ∼ 87 % was the main drained unmanaged peatland grassland field (Fig. 2a).

Measured peat depths across the 46 sampling points ranged from 120 to 387 cm. Spatially, deeper peat deposits were primarily concentrated in the central and eastern parts of the site. In contrast, shallower regions were typically located near the northern and southern site boundaries. The average depth recorded from field measurements was 235 cm, while the mean depth derived from the generated peat depth map was 237 cm (Fig. 2b).

Temporal trends in NDVI followed a clear seasonal pattern, with values peaking in mid-summer (July–August), reflecting maximum canopy development and light interception. NDVI values were lowest in early spring, late autumn and winter corresponding to the dormant phases and the deciduous character of the vegetation and relative lack of moss cover in the fertile grassland (Fig. 2c).

Figure 3 Meteorological and environmental conditions at the LD site (unmanaged drained peatland) in W-Iceland. Time series of key environmental parameters measured from January 2023 to April 2025 including (a) air temperature (T_a; °C) and soil temperature (T_s; °C), (b) photosynthetically active radiation (PAR; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), (c) vapor pressure deficit (VPD; hPa), (d) precipitation (Rain; mm), (e) soil water content (SWC; m³ m⁻³), and (f) water level (WL; m). Vertical red dashed lines show the transitions between the growing seasons (GS) and non-growing seasons (NGS).

3.2 Environmental variables

The seasonality of key environmental drivers, including T_a, T_s (at 10 cm below the surface), PAR, VPD, Rain, SWC and WL, was analyzed to provide essential context for interpreting the ecosystem CO₂ fluxes (Fig. 3).

The GS began in early May and extended into early October in 2023 (3 May–5 October), reflecting a relatively long period of active plant growth. In 2024 the GS was noticeably shorter; it started later and ended earlier, covering the period from 19 May to 22 September (Fig. 3; vertical red dashed lines). This corresponds to the length of the GS being 45 % of the year in 2023 and 35 % in 2024.

Mean T_a during GS was slightly higher in 2023 (9.8 °C) than in 2024 (9.5 °C), with similar patterns during the NGS with mean T_a of 1.7 °C in 2023 and 1.6 °C in 2024. T_s followed similar trends, remaining low during NGS (mostly close to 0 °C due to freezing of surface peat layers) and reaching GS means of 10.6 °C (2023) and 10.4 °C (2024) (Fig. 3a).

PAR showed the strongest seasonal signal (Fig. 3b), with GS means of 302 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ in 2023 and 281 in 2024, while NGS values remained low. VPD remained relatively low in both years (with GS mean of 5.8 hPa in 2023 and 5.1 hPa in 2024 and even lower values during NGS), reflecting the humid maritime climate, but showed slightly stronger seasonality in 2023 (Fig. 3c).

Total annual precipitation was higher in 2024 (950 mm) than in 2023 (823 mm), an increase of roughly 15 % (Fig. 3d). During the GS, rainfall increased from 328 mm in 2023 to 391 mm in 2024, with 2023 characterized by a pronounced dry period in July–August, whereas precipitation in 2024 was more evenly distributed.

SWC reflected the precipitation pattern, with reduced values during the dry period of 2023 and consistently higher, more stable levels (∼ 0.4 m³ m⁻³) during GS 2024. Lower winter SWC in both years was associated with soil freezing, consistent with T_s values near 0 °C.

The WL fluctuated between approximately −0.5 and −0.9 m through most of the study period, close to the subterranean channel depth of −0.6 m (Fig. 3f). During GS 2023 the mean WL was −0.88 m, while GS 2024 showed a shallower mean of −0.73 m, indicating that the GS WL was on average about 17 % higher in 2024. The deepest drawdowns occurred in late GS of 2023, reaching approximately −1.1 m, while the deepest drawdowns in 2024 was only −0.86 m.

Overall, 2023 was warmer and drier compared to 2024, which was colder and had more precipitation especially during GS.

3.3 Data quality, gap filling, and partitioning

In EBC evaluation, linear regression of half-hourly H + LE against R_n − G yielded a slope of 0.79 and an intercept of 2.68 W m⁻² in 2023 and a slope of 0.74 and an intercept of 2.25 W m⁻² in 2024, indicating reasonable EBC within the expected range for peatland and wetland ecosystems (Wilson et al., 2002).

Figure 4 Time series of measured and gap-filled CO₂ flux (Net Ecosystem Exchange, NEE) at the LD site (unmanaged drained peatland) in W-Iceland from January 2023 to April 2025. The data includes the original quality-controlled observations (blue) and the continuous record reconstructed using the hybrid MDS-XGBoost gap-filling approach (black).

The final gap-filled time series (Fig. 4) combined accepted MDS estimates with XGBoost predictions, producing a smooth and internally consistent NEE dataset. Model evaluation demonstrated strong predictive performance, with R² values of 0.8 for MDS model and 0.91 for XGBoost, and RMSE below 1.36 $\mathrm{\mu }\mathrm{mol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, highlighting the high level of agreement between predicted and observed fluxes. Based on SHAP values, radiation and thermal variables emerged as the primary contributors, while seasonal and diurnal temporal indicators were additionally used by the model to capture systematic patterns in ecosystem exchange not fully resolved by meteorological drivers alone. A detailed breakdown of XGBoost feature attribution based on SHAP values is provided in Appendix B.

The gap-filled NEE was partitioned into flux components (R_eco and GPP). Evaluation metrics showed an acceptable model fit, yielding an R² of 0.59 and an RMSE of 0.913 ${\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (see Appendix B for the complete SHAP value analysis). These continuous NEE, GPP and R_eco data provided a robust foundation for analyzing seasonal dynamics and annual carbon budgets.

Figure 5 Mean diurnal cycles of carbon flux components at the LD site (unmanaged drained peatland) in W-Iceland (2023–2024 data). Diurnal patterns of net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (R_eco) are shown for (a) the growing season (GS) and (b) the non-growing season (NGS). Shaded areas represent the 95 % confidence intervals.

3.4 Diurnal NEE, GPP and R_eco

As expected, the average diurnal cycle of CO₂ fluxes showed a clear seasonal contrast, reflecting the interplay between seasonal light availability, temperature, and phenological state (Fig. 5). During GS, the ecosystem acted as a strong daytime carbon sink, with NEE and GPP exhibiting synchronized diurnal patterns. Peak GPP uptake and maximum NEE sequestration occurred at 13:00 GMT, reaching −12.1 and −6.9 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively, close to local solar noon. In contrast, R_eco peaked later at 16:00 GMT (5.4 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), reflecting thermal inertia and the delayed response to peak air and soil temperatures. The broad “shoulders” of the GPP curve indicate extended photosynthetic activity during long twilight periods at this high-latitude site.

In contrast, during the NGS, photosynthetic activity was strongly suppressed, and the ecosystem functioned as a persistent CO₂ source (Fig. 5). GPP remained low, with a maximum uptake of −1.4 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (≈ 10 % of GS peak), while R_eco exceeded GPP throughout the day, resulting in consistently positive NEE values, peaking at 1.2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ around midday. Although the diurnal amplitude of R_eco was reduced compared to the GS, it still showed a cycle, with maximum values of 2.7 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (∼ 50 % of GS peak). The sustained positive NEE during the extended NGS highlights the significant contribution of winter CO₂ emissions to the annual carbon balance.

Figure 6 Daily time series of net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (R_eco), over the study period (January 2023 to April 2025) and cumulative sum of NEE for 2023 and 2024 at the LD site (unmanaged, drained peatland) in W-Iceland. Negative cumulative NEE values show a net carbon sink, while positive values indicate a net carbon source. Vertical red dashed lines show the growing season (GS) and non-growing season (NGS) periods.

3.5 Seasonal NEE, GPP and R_eco

The daily NEE time series showed that the unmanaged drained grassland acted as a net CO₂ source during most of the year (Fig. 6a). During the NGS, NEE remained positive, fluctuating around 2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, indicating continuous CO₂ release under the long, cold and light-limited conditions. With the onset of the GS, increasing temperature and PAR drove NEE toward negative values, reflecting net carbon uptake. The strongest sink strength was observed in July (∼ −4 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), corresponding to peak photosynthetic activity.

GPP showed particularly strong seasonality (Fig. 6b). Photosynthetic uptake started to increase sharply in May and peaked around −10 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ during July. The decline in late September was also rapid, marking the end of active vegetation growth under sharply shortening day length at high latitude.

R_eco followed a similarly clear seasonal pattern (Fig. 6c). R_eco increased rapidly in late spring (May), reached peak values of 6–7 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ around August, and declined gradually during autumn (late September). Winter respiration generally remained relatively high, (∼ 2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), occasionally falling below 1 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, but never ceasing entirely.

Despite a mid-season reduction in cumulative NEE reflecting temporary carbon uptake, emissions during the NGS dominated the annual balance (Fig. 6d). The consistently increasing cumulative NEE in both years confirms the system as a persistent net CO₂ source, with slightly higher total emissions observed in 2024. Overall, these flux patterns revealed an ecosystem characterized by intense but short summer CO₂ uptake on one hand, and continuous and substantial winter respiration on the other hand.

Table 3 Annual and seasonal carbon budgets for net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (R_eco) in t CO₂-C ha⁻¹ at the LD site (unmanaged drained peatland) in W-Iceland during 2023 and 2024, partitioned into the growing season (GS), non-growing season (NGS), and the full calendar years. The GS duration was 156 d in 2023 and 127 d in 2024, while the NGS lasted 190 d in 2023 and 239 d in 2024.

3.6 Carbon budgets

The annual net C-budget (NEE) was surprisingly similar between the 2 years (Table 3), showing the drained unmanaged peatland being a C source to the atmosphere of 4.1 and 4.4 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in 2023 and 2024, respectively (approximately an 8 % difference). All the net seasonal CO₂-C loss from the unmanaged drained peatland occurred during the NGS, as the GSs were small sinks in both years.

When comparing the seasonal NEE budget between the 2 years, it should be noted that the NGS was shorter in 2023 than in 2024 (55 % versus 65 % of the year, respectively). The NGS NEE-sum was correspondingly lower in 2023 (Table 3). The GS NEE-sum was also less negative in 2023, when the seasonal temperature variation included a larger part of the late season within the defined GS.

The component fluxes (GPP and R_eco) provided further insight into these annual patterns. The annual GPP was 11 % higher in 2023, which benefited from a longer GS, but the annual R_eco was only 5 % higher in 2023 compared to 2024. The higher $\text{GPP}/R_{\text{eco}}$ ratio in 2023 than in 2024 resulted in the site being a smaller net CO₂ source in 2023 (Table 3).

While most of the GPP fluxes occurred during the defined GS, the R_eco showed a more modest intra-annual variability, with an average of 39 % and 50 % of seasonal emissions occurring during the NGS in 2023 and 2024, respectively (Table 3). Consequently, in both years nearly half of the annual respiratory CO₂ release occurred outside the GS. The average annual R_eco-sum (12.6 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) was approximately 3 times larger than the average NEE-sum (4.2 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$), showing the relative importance of the annual GPP in the site's C-balance (Table 3).

Figure 7 Environmental drivers of carbon fluxes at the LD site (unmanaged drained peatland) in W-Iceland. Relationships are between weekly aggregated fluxes and biophysical drivers for 2023 and 2024 using non gap-filled data. (a) Weekly average ecosystem respiration (R_eco; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) versus average soil temperature at 10 cm (T_s; °C) and (b) R_eco versus water level (WL; m), both restricted to the growing season (GS) due to T_s and WL winter data gaps in 2023. (c) Gross primary production (GPP; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) versus accumulated photosynthetically active radiation (PAR_acc; µmol m⁻²) and (d) GPP versus the weekly radiation–vegetation product (Accumulative PAR × NDVI), both using full annual data.

3.7 Environmental controls on CO₂ exchange

To explain the relatively small interannual differences in NEE despite contrasting weather conditions (2023 being warmer and drier while 2024 was colder and wetter), weekly aggregated R_eco and GPP values were related to measured environmental drivers (Fig. 7). For T_s and WL, the analysis was restricted to observed GS data to ensure comparability between years, as winter data were unavailable for 2023. Conversely, for PAR and NDVI, the full annual dataset was used to capture the complete seasonal development of the radiation-use efficiency of the canopy.

Weekly R_eco consistently followed a similar exponential temperature-response curve in both years relative to T_s (Fig. 7a; Q₁₀ values of 2.37 in 2023 and 2.23 in 2024). While the range of weekly temperatures was also relatively similar across both years, 2024 experienced more weeks with lower T_s and fewer peak-temperature weeks compared to 2023 (Fig. 7a). This difference in frequency of high-temperature weeks primarily explained why R_eco was significantly higher in 2023 (annual median R_eco of 4.087 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ in 2023 and 2.751 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ in 2024; Mann–Whitney U test, p < 0.0001). Additionally, a significant negative relationship between WL and R_eco was observed in 2023 (R² = 0.40, p = 0.009), whereas no such relationship existed in 2024 (R² = 0.00, p = 0.934) (Fig. 7b). Consequently, the higher cumulative R_eco in 2023 was also driven, to a lesser extent, by the lower WL, which increased the peat volume exposed to aerobic microbial decomposition.

Weekly GPP reached higher peak values in 2023 compared to 2024 (Fig. 7c). PAR_acc alone explained 72 % of the variation in GPP during 2023 but only 38 % in 2024, where several weeks exhibited a low GPP:PAR_acc ratio following a delayed start to the GS. However, when GPP was plotted against the PAR_acc × NDVI product, the relationships for both years converged and explanatory power increased to 86 % and 64 % for 2023 and 2024, respectively (Fig. 7d). 2023 included several weeks with higher PAR_acc × NDVI values than any observed in 2024. Further relationships between the R_eco and GPP with other environmental variables are provided in Appendix C.

4 Discussion

4.1 Net annual C balance

The present study site is broadly representative of unmanaged drained peatlands in Iceland, where drainage commonly leads to a shift from wetland vegetation toward graminoid-dominated grassland communities, often influenced by land use such as grazing (Arnalds et al., 2016a). These systems are classified as “Grassland on organic soil” in the Icelandic National Inventory Report (Umhverfisstofnun, 2023) and cover approximately 4200 km².

The results show that the drained unmanaged peatland remained a consistent net CO₂ source despite different environmental conditions between 2023 (warmer and drier) and 2024 (colder and wetter), with a mean annual loss of 4.2 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 3). This value was 26 % lower than the IPCC (2014) default boreal EF of 5.7 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ that is currently used in Iceland's national reporting for unmanaged drained peatlands (Umhverfisstofnun, 2025). It is, however, notably similar to the 3.8 ± 0.65 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ that was obtained by chamber measurements on unmanaged drained peatland in western Iceland (Ólafsdóttir, 2015). Further, it is a similar value to estimated mean soil organic matter losses of approximately 4 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}\text{-}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in a regional study using chamber point measurements and NDVI-driven upscaling on similar seven unmanaged drained sites in W-Iceland (Guðmundsson et al., 2026). The obtained CO₂-C loss in the present study is, however, higher than the mean annual loss estimated from the peat-stock change method in eight paired undrained and drained site-pairs in southern Iceland, which found on average annual topsoil loss of 1.7 t CO₂-C ha⁻¹ during 15–50 years of draining (Gunnarsdóttir, 2017). However, in that study a 1500 AD ash layer at 17–30 cm depth in the peat profile was used for the stock-change estimates, which may have underestimated the total CO₂ loss from the whole drained peat profile as the WLs were generally lower.

Despite the similar NEE observed in this study compared to other measurements on drained peatlands in Iceland, the annual GPP (−7.85 to −8.8 t CO₂-C ha⁻¹) and R_eco (12.2 to 12.7 t CO₂-C ha⁻¹) observed here were somewhat higher than most chamber-based estimates from other Icelandic drained sites (GPP: −4.3 to −7.4 t CO₂-C ha⁻¹ and R_eco: 7.19 to 12.3 t CO₂-C ha⁻¹) (Guðmundsson et al., 2026; Ólafsdóttir, 2015), likely reflecting the broader spatial footprint and continuous temporal coverage of the EC method.

One pertinent question is whether the 60-year duration since initial drainage contributes to the lower CO₂ emissions observed compared to the IPCC (2014) default EF. There is currently a general lack of long-term monitoring studies evaluating how long high CO₂-C emissions are sustained post-drainage (Qiu et al., 2021). Truskavetskii (2014) demonstrated that the loss of peat organic matter in croplands was linear during the first 45 years following drainage, excluding enhanced losses during the initial 7-years period. Studies extending beyond 50 years have primarily focused on how decadal peat subsidence may reduce effective drainage (WL) over time, leading to decreasing rates of organic matter loss as the peat column compacts (Rojstaczer and Deverel, 1993; Liu et al., 2026). Evans et al. (2021) analyzed CO₂ flux measurements from 16 EC towers on drained peatlands across the UK and concluded that mean annual effective WL represents the overwhelmingly dominant control on CO₂ fluxes.

Given that the subsidence at the present study site was only −0.16 m since the onset of drainage (Streeper, 2026), and the present WL generally remained deep (mainly between −0.5 and −0.7 m), the age of the drainage (60 years) is unlikely to be the primary explanation for the lower CO₂-C losses relative to the IPCC (2014) default EF for the boreal zone. Furthermore, both Gunnarsdóttir (2017) and Guðmundsson et al. (2026) investigated multiple drained peatland sites in Iceland using stock change and repeated chamber flux methods, respectively. They observed considerable variation in annual CO₂-C losses among sites that was not correlated with the age of drainage. Consequently, it is recommended that for future assessments of unmanaged drained peatlands, sites should be classified by their effective annual WL rather than the time that has elapsed since drainage.

4.2 Winter fluxes

The annual CO₂-C balance of the unmanaged drained peatland was strongly influenced by NGS emissions, which accounted for 39 % and 50 % of the annual R_eco in 2023 and 2024, respectively. Notably, positive NEE occurred only during the NGS in both years, which turned the study site into an annual C source. For comparison, Alm et al. (1999) found that winter (November-May) CO₂-C losses accounted for 23 % and 21 % of annual totals in undisturbed and drained peatlands, respectively, at latitudes of 62–65° N in Finland using a closed chamber method. In another study Tikkasalo et al. (2025) found by EC measurements over a drained peatland forest in southern Finland (latitude of 61° N), that only 28 % of the positive NEE was released during the winter season. Although the exact timing of NGS defined in the present study (October–April) differed from these reference studies (November–May), the duration was the same. This comparison suggests that the proportion of NGS CO₂ flux relative to the annual flux is higher in the present study than other studies on drained peatlands at similar latitudes in Fenno-Scandia. Climatic shifts in these continental regions are most evident during the winter, a season that has warmed markedly over the past 30 years and is projected to experience the greatest degree of warming through the end of the century (Lind et al., 2023). Recent evidence suggests that the C balance of drained peatlands in Finland has already shifted toward becoming a larger CO₂ source over the last three decades, due to warmer (winter) climates (Alm et al., 2023). In this context, Iceland may serve as a pertinent case study, illustrating how increasingly mild winters are likely to impact the annual C balance of drained peatlands at similar latitudes globally.

4.3 Environmental controls on CO₂ fluxes

The year 2023 was an extreme weather year at the field site, characterized by the coldest winter since 1995, followed by an unusually warm and dry summer (including the driest July since the 1860s) and a warm autumn (Icelandic Meteorological Office, 2025a). The 2024 winter was also among the coldest on record but was followed by an unusually wet and cold summer (the coldest since 1998) and autumn (Icelandic Meteorological Office, 2025b). Despite these contrasting meteorological conditions between 2023 (drier and warmer) and 2024 (wetter and colder), the annual CO₂-C balance remained surprisingly similar across both years (Table 3). The stabilizing mechanism that appears to have balanced the enhanced CO₂ fluxes during the warm 2023 was a strong positive response in the productivity of the unmanaged grassland vegetation. This was reflected in higher PAR_acc × NDVI (indicating enhanced canopy development and higher irradiance) and elevated GPP values during a 6-week window of the relatively short GS (Fig. 7d). Although R_eco was significantly higher in 2023 due to higher temperatures and reduced precipitation, which led to lower WL and greater peat exposure to aerobic microbial decomposition, the enhanced GPP largely offset this increase. Consequently, the annual NEE remained consistent across both years. High covariance between R_eco and GPP which led to small interannual differences in NEE has been observed in other studies (Wohlfahrt et al., 2008). Previous research has shown that the NEE of drained peatlands can be relatively stable (e.g. Wilson et al., 2016), while other studies have highlighted significant interannual variability in NEE (e.g. Bubier et al., 2003). These differences in results often depend on whether interannual driving factors consistently stimulate both C uptake and release, or primarily influence only one side of the C balance.

Unmanaged grasslands at high latitudes, particularly in Iceland, have shown a robust positive response to warming, characterized by significant increases in biomass production and NDVI (Elmendorf et al., 2012; Mortier et al., 2024). This trend has also been noted for the country as a whole as MAT has increased since the 1980s (Raynolds et al., 2015). Such a strong positive relationship between GPP (or NDVI) and temperature is not unexpected at the northern limits of the boreal/Subarctic zone, particularly where water availability or biotic factors, such as high grazing pressure, do not limit aboveground biomass.

The summer drought in 2023 resulted in a significant drawdown of WL; however, no apparent negative drought response was observed in either GPP or R_eco. In fact, R_eco showed a significantly enhanced response to the lower WL and higher T_s. This aligns with findings from 16 EC towers on drained peatlands across the UK (Evans et al., 2021), which concluded that mean annual effective water-table depth represents the overwhelmingly dominant control on CO₂ fluxes in these ecosystems. In the present study, however, the effect of lower WL was relatively small compared to the direct T_s effect on R_eco.

5 Conclusion

The present study showed that the drained Icelandic unmanaged peatland consistently released CO₂ to the atmosphere, with annual C loss being 26 % less than the IPCC (2014) default EF for drained grasslands in the northern boreal zone. The CO₂-C loss was driven by relatively high respiration fluxes outside the growing season rather than by differences in summer productivity and CO₂ emissions. Despite the clear contrast between a dry, warm year (2023) and a cool, wet year (2024), the overall carbon balance remained remarkably similar. This stability reflects the distinctive ecological setting of Icelandic peatlands: mild maritime winters that only partly inhibit microbial respiration, and cool summer climate where warm summers can greatly enhance vegetation growth and the growing season GPP. Together, these findings highlight the Icelandic peatlands function differently from many boreal or continental peatlands, where severe winters suppress respiration. As such, accurate carbon budgets for Icelandic sites require year-round flux measurements. The results also highlight the need for restoration measures that target reductions in non-growing-season respiration if long-term CO₂ losses from drained peatlands are to be mitigated.

Appendix A: Fuel cell placement and interference check

To ensure the EFOY fuel cell would not interfere with the EC system, its location was carefully selected based on the site's prevailing wind directions. As confirmed by the drone imagery (Fig. A1a), the fuel cell was positioned to the NW of the EC system, a sector from which the wind rarely blows. This strategic placement successfully mitigates potential CO₂ pollution from the device, a conclusion supported by the 2023–2024 dataset (Fig. A1b) demonstrating minimal atmospheric CO₂ concentration originating from the NW sector. Furthermore, this lack of interference is further supported by isolating data from the NGS of 2023 and 2024 (Fig. A1c), which consistently shows that CO₂ contributions from the NW remain negligible throughout the year.

Figure A1 Assessment of EFOY fuel cell placement and potential CO₂ interference. (a) Location of the EFOY fuel cell situated to the northwest (NW) of the Eddy Covariance system at the LD site (unmanaged drained peatland) in W-Iceland. (b) Atmospheric CO₂ concentrations originating from different wind sectors during the full 2023–2024 period. (c) CO₂ concentrations from different wind sectors isolated to the non-growing seasons in 2023 and 2024.

Appendix B: Machine learning interpretability

Figure B1 shows the global feature attribution of the NEE gap-filling model, quantified using mean absolute SHAP values. This metric summarizes the relative importance of each predictor in explaining the variability of the CO₂ flux residuals, where larger values indicate a stronger influence on the model's predictions. PAR emerged as the dominant predictor, reflecting its primary role in controlling photosynthetic uptake during the growing season. The high ranking of seasonal (sin_doy) and diurnal (cos_hour, sin_hour) harmonic transforms shows that the model captures recurring phenological and diurnal patterns in ecosystem exchange that are not fully resolved by instantaneous meteorological variables alone. Furthermore, the moderate importance of T_s and T_a highlights the underlying control of the thermal regime on ecosystem respiration, while the relatively low contribution of VPD suggests that atmospheric moisture demand was a secondary driver of flux variability at this high-latitude site during the study period.

Figure B2 shows the global feature attribution of the R_eco model, quantified using mean absolute SHAP values. T_s emerged as the dominant predictor, confirming the primary control of soil thermal conditions on respiration processes at this site. The importance of T_a and the seasonal harmonic transform (sin_doy) further highlighted the model's ability to capture broader thermal regimes and phenological cycles. SWC showed moderate importance, reflecting its secondary but still meaningful role in modulating respiration. In contrast, the low contribution of VPD and diurnal proxies (cos_hour, sin_hour) suggested that atmospheric moisture demand and high-frequency hourly variability were minor drivers of respiration compared to soil temperature and seasonal dynamics.

Figure B1 Mean absolute SHAP values for the XGBoost CO₂ flux gap-filling model applied to data from LD site (unmanaged, drained peatland) in W-Iceland during 2023–2025, showing the relative importance of each predictor. The results show that photosynthetically active radiation (PAR) is a strong predictor. Seasonal and diurnal harmonic proxies (sin_doy, cos_hour) and thermal variables (T_s, T_a) exhibit moderate importance, while vapor pressure deficit (VPD) contributes the least to the predicted flux variability.

Figure B2 Mean absolute SHAP values for the ecosystem respiration (R_eco) model applied to data from the LD site (unmanaged, drained peatland) in W-Iceland during 2023–2025, showing the relative importance of each predictor. The results indicate that soil temperature (T_s) is the strongest predictor of respiration. Air temperature (T_a) and seasonal proxies (sin_doy) exhibit moderate importance, followed by soil water content (SWC). Diurnal harmonic transforms and vapor pressure deficit (VPD) contribute the least to the predicted respiration variability.

Appendix C: Secondary environmental controls on ecosystem exchange

Figure C1 shows the relationships between various environmental drivers and NEE components (R_eco and GPP) calculated as weekly means for the 2023 and 2024 periods using observed (non-gap filled) data. In 2023, R_eco had significant negative relationships with SWC (R² = 0.48, p = 0.0029), whereas in 2024, this relationship was statistically non-significant (R² = 0.06, p = 0.36). A distinct negative relationship was observed between R_eco and VPD in 2023 (R² = 0.28, p = 0.016), while no relationship was found in 2024 (R² = 0.00, p = 0.85). Furthermore, the relationship between R² and NDVI was strong in both years (p < 0.0001), though it was stronger in 2023 (R² = 0.88) than in 2024 (R² = 0.70).

GPP showed a higher sensitivity to T_s in 2023 (R² = 0.79, p < 0.0001) compared to 2024 (R² = 0.61, p = 0.0006). In contrast, GPP exhibited a significant positive relationship with VPD in 2024 (R² = 0.17, p = 0.0041), while the relationship remained non-significant in 2023 (R² = 0.03, p = 0.2676). The relationship between GPP and SWC remained statistically non-significant in both years (p > 0.05). Overall, 2023 was characterized by stronger relationships with hydrological and thermal drivers, while 2024 showed more moderate relationships across most secondary environmental variables.

Figure C2 shows overall relationships between observed environmental drivers and NEE components (R_eco and GPP) calculated as weekly means for the combined 2023–2024 dataset, quantified using Pearson correlation analysis. T_s showed the strongest correlations across all flux components, with a near-perfect positive correlation with R_eco (r = 0.96) and a very strong correlation with GPP (r = 0.91). PAR and T_a were also highly correlated with carbon fluxes (p < 0.001). NDVI shows strong correlations with both GPP (r = 0.75) and R_eco (r = 0.82). WL displays a significant but more moderate correlation with R_eco (p < 0.001). VPD and SWC show non-significant correlations at the weekly timescale.

Figure C1 Environmental drivers of carbon fluxes at the LD site (unmanaged drained peatland) in W-Iceland. Relationships are between weekly aggregated fluxes and biophysical drivers for 2023 and 2024. (a) Weekly average ecosystem respiration (R_eco; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) versus soil water content at 10 cm depth (SWC; m³ m⁻³), (b) R_eco versus vapor pressure deficit (VPD; kPa), (c) R_eco versus normalized difference vegetation index (NDVI), (d) Gross primary productivity (GPP; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) versus soil temperature at 10 cm depth (T_s10; °C), (e) GPP versus VPD, (f) GPP versus SWC. The plots involving T_s or SWC are based on measured growing season (GS) data due to data gaps during the non-growing season (NGS), whereas all other plots use measured data from the full annual cycle.

Figure C2 Pearson correlation matrix (r) between observed weekly environmental drivers such as air temperature (T_a; °C), soil temperature (T_s; °C), photosynthetically active radiation (PAR; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), vapor pressure deficit (VPD; hPa), precipitation (Rain; mm), soil water content (SWC; m³ m⁻³), and water level (WL; m) and observed (non–gap-filled) CO₂ fluxes including net ecosystem exchange (NEE; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), gross primary production (GPP; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), and ecosystem respiration (R_eco; $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) for the combined 2023–2024 study period at the LD site (unmanaged drained peatland) in W-Iceland. Colores represent the strength and direction of correlations, ranging from blue (negative) to red (positive). Significance levels are shown as follows: ${}^{\ast \ast \ast }$ p < 0.001, ${}^{\ast \ast }$ p < 0.01, ^∗ p < 0.05, and ns (not significant).