Simulating ecosystem carbon dioxide fluxes and their associated influencing factors for a restored peatland

Restoration of drained and extracted peatlands can potentially return them to carbon dioxide (CO₂) sinks, thus acting as significant climate change mitigation. However, whether the restored sites will remain sinks or switch to sources with a changing climate is unknown. Therefore, we adapted the CoupModel to simulate ecosystem CO₂ fluxes and the associated influencing factors of a restored bog. The study site was a peatland in eastern Canada that was extracted for 8 years and left for 20 years before restoration. The model outputs were first evaluated against 3 years (representing 14–16 years post-restoration) of eddy covariance measurements of net ecosystem exchange (NEE), surface energy fluxes, soil temperature profiles, and water table depth data. A sensitivity analysis was conducted to evaluate the response of the simulated CO₂ fluxes to the thickness of the newly grown mosses. The validated model was then used to assess the sensitivity to changes in climate forcing. The CoupModel reproduced the measured surface energy fluxes and showed high agreement with the observed soil temperature, water table depth, and NEE data. The simulated NEE varied slightly when changing the thickness of newly grown mosses and acrotelm from 0.2 to 0.4 m but showed significantly less uptake for a 1 m thickness. The simulated NEE was $-\mathrm{95}\pm \mathrm{19}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ over the 3 evaluation years and $-\mathrm{101}\pm \mathrm{64}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, ranging from −219 to $+\mathrm{54}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, with extended 28-year climate data. After 14 years of restoration, the peatland has a mean CO₂ uptake rate similar to pristine sites but with a much larger interannual variability, and in dry years, the restored peatland can switch back to a temporary CO₂ source. The model predicts a moderate reduction in CO₂ uptake but still a reasonable sink under future climate change conditions if the peatland is ecologically and hydrologically restored. The ability of the CoupModel to simulate the CO₂ dynamics and its thermo-hydro-drivers for restored peatlands has important implications for emission accounting and climate-smart management of drained peatlands.

1 Introduction

Degradation of peatlands through land use change and drainage is currently estimated to emit ∼4 % of global annual anthropogenic carbon dioxide (United Nations Environment Programme, 2022). Therefore, restoring drained peatlands so that they return to carbon (C) sinks has been identified as an emerging priority for climate change mitigation (Leifeld and Menichetti, 2018). When ecologically restored successfully, peatlands can generally return to their carbon uptake function after a decade or two following the recolonization of peatland vegetation and a decrease in water table depth (Nugent et al., 2018; González and Rochefort, 2014; Richardson et al., 2023; Tuittila et al., 1999; Wilson et al., 2016; Beyer and Höper, 2015). However, the C-uptake function of restored peatlands is sensitive to climate conditions, particularly in drier years (Wilson et al., 2016). Therefore, a changing climate can potentially weaken the sink strength or even switch the restored peatlands to C sources.

In North America, about a quarter of drained peatlands that were used earlier for horticultural peat extraction have been restored by the moss layer transfer technique (MLTT) (Chimner et al., 2017; Quinty and Rochefort, 2003). Ecosystem-scale flux measurements indicate peatlands remain a CO₂ source (∼200–$\mathrm{500}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) the first few years of restoration (Petrone et al., 2001, 2003), but after a decade or two, peat vegetation recovers, and the restored bogs return to CO₂ sinks with uptake rates similar to pristine sites (Nugent et al., 2018). While the C-accumulation function can generally be fully restored within a decade or two, full restoration of the peat soil structure and ecohydrology takes a much longer time (Loisel and Gallego-Sala, 2022) with centuries to millennia required for the restored peatland to accumulate the C that was extracted. Restoration creates a novel ecosystem in transition to a rewetted steady state, and the altered ecohydrology decreases peatland ecological resilience (Kreyling et al., 2021).

The ecological function of peatlands is strongly linked to ecohydrology (Waddington et al., 2014). Clymo (1992) outlined four functional layers of pristine peatlands (i.e., green, peat litter, collapse, peat proper; Fig. 1 in his paper) and how the peat structure interacts with ecohydrology, thus regulating the growth and function of peatlands. Briefly, the bulk density of the green and peat litter layer is low, typically below 0.05 g cm⁻³. The increasing load of new growth above and the mass proportion of water, as well as the decomposition of plant material, cause the moss structure to collapse, typically increasing the bulk density gradually along the peat profile to $\sim \mathrm{0.1}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$. The result is a reduction in the space between dead leaves and stems and the soil pore sizes, increasing the capillary force for vertical water movement, thus sustaining the water supply for sphagnum mosses and the growth of the peatlands. For extracted peatlands, the MLTT gives a jump start for moss colonization at the residual catotelmic peat surface; with time, a new layer of acrotelm is formed and thickens. However, these newly regenerated mosses with low bulk density form large pores directly above the dense residual peat remaining after extraction (catotelmic peat) and do not have the negative interstitial pressures required to draw pore water, causing a capillary barrier effect (Gauthier et al., 2018, 2022). The capillary barrier decreases the ability of the new moss to draw water from the deeper compacted catotelmic peat, thus causing an overall lower surface moisture content for restored sites compared to natural peatlands (McCarter and Price, 2015). As a result, the new moss layer may become stressed quickly and even die off during dry periods. Synthesis studies have shown that vegetation colonization is much slower after restoration over warm and drier years (González and Rochefort, 2014), and data from a restored Irish extracted bog show a less resilient C-uptake function over the drier years (Wilson et al., 2016).

Under the United Nations Framework Convention on Climate Change (UNFCCC), countries with peatlands managed for extraction are required to report greenhouse gas emissions annually (IPCC, 2014). Currently, the National Inventory Report (NIR) of Canada reports emissions from restored peatland separately, and an emission factor (EF) of +2.07 t CO₂–C ${\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (positive meaning source) generated from data of three sites (all restored less than 10 years) is used (ECCC, 2021). However, the CO₂ emissions change with time as the peatland develops and gradually switches to CO₂ uptake (Nugent et al., 2018). Currently, there is a discussion that restoration can create C credits and thus could be used to offset the C emissions during the drainage phase (Tanneberger and Wichtman, 2011). Moreover, the IPCC and the Canadian NIR use default EFs (i.e., Tier 1) based on literature data (IPCC, 2019). An emission factor based on empirical observations (i.e., Tier 2) offers improvement as it is subject to the environmental conditions and the time of year the measurements were done. Yet most of the observed data are of short duration and thus cannot capture interannual variations in emissions and associated environmental variables. Process-based modeling of restored peatlands (i.e., Tier 3) can be used to determine the “representativeness” of the empirical EFs by examining the coupled hydrological–C dynamics and how they vary within and between years. He and Roulet (2023) showed that directly using literature data to generate emission factors can be biased because it does not account for seasonality and interannual climate variability.

Existing studies using models for restored peatlands are few. Lees et al. (2019) applied a satellite-based, temperature-driven gross primary productivity (GPP) model over peatland sites at various stages of restoration in the UK and Ireland and found that the model can simulate the GPP measured by eddy covariance. Premrov et al. (2021) modified the drainage function in the ECOSSE model to simulate the water table and CO₂ flux for drained and rewetted extracted bogs, but their model evaluations showed that ECOSSE still requires further development to accurately simulate the water table depth for the rewetted sites. Recently, Lippmann et al. (2023) introduced a dynamic vegetation scheme in the Peatland-VU-NUCOM (PVN) model, driven by input water table data, and evaluated the model for the measured CO₂ flux together with the vegetation competitions in two restored nutrient-rich peatlands in the Netherlands. However, none of these models consider the coupled ecohydrology and C dynamics for restored peatlands (Silva et al., 2024). Previous research showed that the CoupModel could successfully simulate peatland CO₂ dynamics associated with various land use options, e.g., drained peatlands for forestry (He et al., 2016a, b; Kasimir et al., 2021); land use change of afforested peatlands (Kasimir et al., 2018); and five European peatlands with various land uses, including restored sites (Metzger et al., 2015). Recently, the model was applied to simulate the CO₂ fluxes of a pristine continental bog (He et al., 2023b) and an active peat extraction site (He et al., 2023a). These studies provide a basis for further use of the model to simulate restored peatlands to close the land use cycle – from pristine peatlands to drainage for different land uses, culminating in final restoration.

Table 1 Data from the Bois-des-Bel peatland used for the CoupModel forcing and evaluation.

The overall aim of this study is to simulate the soil–atmosphere exchanges of heat, water, and CO₂ for a bog restored by the MLTT technique. More specifically, we aim to do the following:

1. adapt and evaluate the CoupModel to simulate net ecosystem exchange (NEE) and its hydrothermal drivers, including surface energy fluxes, soil temperature profile, and water table depth;

2. test the model sensitivity to varying thickness of newly grown mosses and the acrotelm;

3. evaluate the impact of interannual climate variability on the simulated ecosystem CO₂ flux and discuss its implications for emission factor calculation; and

4. predict the impact of future climate change on the C-uptake function of restored peatlands.

2 Site and methods

2.1 Site description

The Bois-des-Bel (BDB) peatland is located 11 km northeast of Rivière-du-Loup, Quebec ($\mathrm{47}\mathrm{°}{\mathrm{58}}^{\prime }{\mathrm{1.95}}^{\prime \prime }\mathrm{N}$, $\mathrm{69}\mathrm{°}{\mathrm{25}}^{\prime }{\mathrm{43.10}}^{\prime \prime }\mathrm{W}$). The peatland complex covers an area of 202 ha with a mean peat depth of 2.2 m. A small sector of 11 ha was extracted for horticulture peat by vacuum harvesting between 1972 and 1980. After the extraction, there was 2 m residual peat left where the top 0.8 m characterizes a Sphagnum bog peat (Lavoie et al., 2001). In the autumn of 1999, an 8.1 ha area was restored using the MLTT. The climate of the region is cool–temperate with an average long-term (according to 1981–2010 climate-normal data for Saint-Arsène) annual temperature of 3.5 °C and annual precipitation of 962 mm (ECCC, 2023). BDB is well studied, and detailed descriptions of the restoration procedure and site characteristics can be found in several publications (McCarter and Price, 2015; Strack and Zuback, 2013; Waddington and Day, 2007; Poulin et al., 2013). Nugent et al. (2018) measured the soil–atmosphere exchanges by eddy covariance between 2013–2016, i.e., 14–17 years after the restoration. In this study, we used their measured meteorological data (Table 1) for model forcing and measured water table depth, peat temperatures, and flux data for model evaluation.

2.2 Brief model description

The CoupModel (a coupled heat and mass transfer model for soil–plant–atmosphere systems) platform is a process-based model designed to simulate water and heat fluxes, along with the C–N–P cycle, in terrestrial ecosystems (Jansson, 2012; He et al., 2021). The main model structure is a one-dimensional multi-layered soil profile. Model forcing is measured weather data (Table 1). The model and technical description are freely available at http://www.coupmodel.com (last access: 20 February 2024). The CoupModel was previously applied to simulate ecohydrology and CO₂ exchanges for a pristine bog, Mer Bleue, which resembled, though with fewer trees, the BDB site before opening for extraction (He et al., 2023b). It also recently successfully simulated one ongoing peat extraction site, Rivière-du-Loup, in the same region as BDB (He et al., 2023a). The setup and model structure of the BDB simulation were thus built on the base of the upper aerobic peat layer and vegetation characteristics of Mer Bleue and the residual extracted peat layer of Rivière-du-Loup. The details of model parameter configuration for BDB are reported in Table S1 in the Supplement. Here, we report the model setup unique to the BDB site. More detailed process descriptions, model structure, and parameters are reported in He et al. (2023a) and He et al. (2023b).

2.3 Simulation design, model structure, and initial and boundary conditions

The CoupModel was used to simulate the soil vegetation processes and linked hydrology and energy flows of BDB in a 30 min time step from 14 July 2013 to 1 November 2016.

Table 2 Physical, hydraulic, and Mualem–van Genuchten coefficients for the Bois-des-Bel site.

Bulk density ρ_B, porosity θ_s, and saturated conductivity k_sat data are from McCarter and Price (2013), Gauthier et al. (2022), and Petrone (2002). Nonlinear curving fitting was run with the empirical constant $m=\mathrm{1}-\mathrm{1}/n$ with the wilting point θ_w set to 10 % for the topsoil layer and 30 % for the 40–150 cm layer (Menberu et al., 2021).

Nugent et al. (2018) conducted a detailed vegetation survey and spatial distributions at BDB in 2013, and their results show the vegetation at the site is quite homogenous across the major survey direction. These data were used to initialize the vegetation conditions in the CoupModel. The survey showed Sphagnum mosses and Polytrichum strictum cover more than 90 % of the surface with a new acrotelm thickness of ∼0.3 m, sedges (Eriophorum vaginatum and Carex spp.) cover 33 %, and ericaceous shrubs (Chamaedaphne calyculata, Rhododendron groenlandicum, Kalmia angustifolia, Vaccinium oxycoccus, and V. angustifolium) cover 39 % of the soil surface (Nugent et al., 2018). Trees (Picea mariana and Larix laricina) were few but were also beginning to expand across the site. Typha latifolia from the filled remnant ditches covers 4 % of the total site area. In our simulation, we grouped vegetation into three plant functional types or modeled vegetation layers: the first group represents the Ericaceae shrubs and the trees, which cover ∼40 % of the surface, with an assumed lowest root depth of 0.5 m. The second group represents the sedges, which cover 33 % of the surface and lowest root depth of 0.35 m. The third group represents the Sphagnum mosses and other nonvascular vegetation (Polytrichum strictum) at the soil surface, which cover 90 % of the soil surface with no roots. These three modeled vegetation layers were described in the model using the “multiple-big-leaves” concept considering dynamic competition in terms of interception of light and uptake of water. For each vegetation layer, plants were conceptually divided into leaf, stem, and coarse root and fine root. For the moss layer, the live capitulum was conceptually viewed as a leaf and the rest as a stem in the model (He et al., 2023b). C and the dynamics of the plant development (e.g., leaf area index, height) are simulated as the interactions between plant and physical driving forces (e.g., how the plant cover influences both aerodynamic conductance for both heat and momentum transfer in the atmosphere and the radiation balance at the soil surface). Since these are oligotrophic ecosystems, the influence of nutrients on C was not considered in this study. The three vegetation groups were pre-run for 14 years to spin up and reach a quasi-steady state (defined as no abrupt takeover or die-offs of one vegetation group).

For the peat soil, we simulated the first 1.8 m of peat in BDB, which includes 0.3 m of the surface newly developed acrotelm and mosses and 1.5 m of the residual extracted peat. We divided the peat soil profile into nine layers: from 0.05 m per layer at the top to 0.80 m per layer at the bottom. For each simulated layer, the peat soil water retention curve and unsaturated hydraulic conductivity were estimated by the Mualem–van Genuchten model (Mualem, 1976; van Genuchten, 1980). The physical and hydraulic properties used in this study were compiled from the measured data from BDB (Table 2). Water flow between soil layers follows Darcy's law as generalized for unsaturated flow by Richards (1931). We additionally simulated bypass flow to account for preferential water flow in the root channels and macropores by using an empirical bypass flow scheme (Jansson et al., 2004). Soil heat flow between soil layers was assumed to be mainly driven by conduction. The CoupModel solves water and heat equations simultaneously within the soil–plant–atmosphere continuum, and water and heat are coupled in a dynamic way to the plant vegetation layers: accounting for feedback interactions between the plant and the environment.

The initial conditions for water and heat were from measured data (Nugent et al., 2018). The initial condition for soil-C stocks for each soil layer was calculated from the measured bulk density and C concentration (assumed 50 %). The total C in the 1.8 m soil profile was 101.8 kg C m⁻² (Table 2). Similar to He et al. (2023a), we used two soil-C pools which differed in substrate quality and hence decomposition rate to model the impact of organic matter quality on soil respiration: labile and refractory soil C. The partitioning ratio between these two pools from Rivière-du-Loup was used for the bottom 1.5 m at BDB, while for the top 0.3 m of newly grown peat, 80 % was assumed to be in the labile pool. The decomposition rate coefficient (Table S1) and its response to temperature and water were kept the same as in He et al. (2023b).

We assumed no vertical water flow for the lower boundary condition (i.e., at 1.8 m depth) due to the very low saturated hydraulic conductivity (Table 2) and assumed a small thermal heat flow across the lower boundary condition for heat. The site was also drained laterally to the outflow at a distance of ∼200 m (Shantz and Price, 2006). The model parameter values were primarily obtained from the measured data, and where they were not available, literature values used in previous model applications were applied (Table S1).

3 Results

3.1 CoupModel evaluation for restored peatland

The CoupModel simulated the half-hourly surface energy balances well, as shown by the high agreement with the measured total radiation, sensible and latent fluxes (coefficient of determination, r²>0.7 for all; Fig. 1a–c), and surface soil heat flux (r²=0.4; Fig. 1d). However, the model tended to overestimate the sensible heat flux and underestimate the latent heat flux, particularly over the periods of spring and earlier summer, where the model simulated a smaller and delayed (∼1 month) increase in latent heat fluxes compared to the measured data (Fig. S1 in the Supplement).

Figure 1 Relationship between simulated and measured 30 min (a) total net radiation, (b) sensible heat, (c) latent heat, (d) soil surface heat flux, (e) soil temperature at 5 cm depth, (f) soil temperature at 20 cm depth, (g) soil temperature at 80 cm depth, (h) water table depth, and (i) net ecosystem exchange over the period 2013–2016 (n=56 600).

The model simulates the measured soil temperature profile over 5–20–80 cm depth generally well, with r²>0.9 for all three soil layers (Fig. 1e–g). However, the model showed difficulty in precisely simulating the soil thawing (i.e., zero curtain effect; Fig. S2 in the Supplement). The simulated temperature started to increase above zero half a month earlier than the measured data for the 20–80 cm depth in 2015 but was delayed for almost 1 month for 2016 (Fig. S1).

Model performance for water table depth was generally worse compared to the energy and temperature variables. However, the model still captured 50 % of the measured variations (r²=0.5; Fig. 1h). The CoupModel generally simulated a smaller-magnitude fluctuation compared to the measured data, and the model–data agreement was better over the summer than over the winter (Fig. 2a). For instance, large infiltration from snowmelt around May was simulated in the model every year but not represented in the measured data (Fig. 2a).

Figure 2 Measured (red) and simulated (black) 30 min (a) water table depth and (b) net ecosystem exchange.

Measured daily net ecosystem exchange data range from $\sim -\mathrm{3}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (negative indicating uptake) during July to a loss of $\sim +\mathrm{2}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ during cloudy days or shoulder seasons (Fig. 2b; note that the flux data are 30 min in Figs. 1i and 2b). The CoupModel reproduced the measured half-hourly NEE data reasonably well (r²=0.64; Figs. 1i and 2b). Nugent et al. (2018) gap-filled the BDB eddy covariance data and estimated an annual C flux of $-\mathrm{90}\pm \mathrm{10}$ (±95 % CI), $-\mathrm{105}\pm \mathrm{7}$, and $-\mathrm{70}\pm \mathrm{7}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in 2014, 2015, and 2016, respectively. The corresponding simulated annual fluxes are −89, −120, and $-\mathrm{75}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, respectively.

3.2 Sensitivity to the thickness of the newly grown mosses

We conducted a sensitivity analysis to evaluate the model responses to the thickness of the newly grown mosses (i.e., new acrotelm), which partly represents the time since the restoration. It has been argued that −100 mbar is the limiting soil moisture pressure head for sustaining moss growth (McCarter and Price, 2012). Three extra model simulations were made based on the reference run (30 cm acrotelm) with new acrotelm thicknesses of 20 cm (∼10 years after restoration), 40 cm (∼30 years after restoration), and 100 cm (hypothetical, to test the empirical threshold of −100 mbar). For the latter two model simulations, peat properties of the 20–30 cm layer in the reference run (Table 2) were assumed for the future extra 10 and 70 cm acrotelm, respectively. The vegetation was assumed to be the same as the reference run and the peat compaction due to the growth of mosses, and decomposition was not considered.

Our sensitivity analysis showed that the simulated NEE uptake increased slightly when changing the new acrotelm thickness from 20 to 40 cm but reduced (meaning less uptake) significantly for the model run with an acrotelm of 100 cm (Fig. 3). The small changes of simulated NEE can be explained by both increase in GPP and ecosystem respiration (ER) with increasing new acrotelm thickness (20–40 cm). The NPP (net primary production) of mosses shows a slight decreasing trend with increasing acrotelm thickness (Fig. 3). The reduction in CO₂ uptake in the 100 cm acrotelm thickness model run is because the model simulated that the surface mosses start to die off because they cannot take up water from the deep peat (Fig. 3).

Figure 3 The response of simulated average ecosystem CO₂ fluxes (2013–2016) to the simulated thickness of the newly grown mosses.

3.3 Impact of interannual climate variability on CO₂ uptake of restored peatlands

The BDB region shows large annual climate variability over the last 28 years from 1994 to 2021. The measured annual mean air temperature ranged from 2.6 to 5.7 °C and the annual precipitation from 633 to 1488 mm (Fig. 4a and b). This can be compared to the 30-year annual mean air temperature of 3.5±2.9 °C and the precipitation of 962 mm for the climate-normal data (1981–2010) at Saint-Arsène station (Environment and Climate Change Canada, 2023). Both annual air temperature and precipitation showed increasing trends over the measured period from 1994 to 2021, with an increasing rate of 0.03 °C yr⁻¹ for air temperature and 1.69 mm yr⁻¹ (Fig. 4), indicating possible future warmer and wetter conditions in the region. The weather over the 3 years of flux measurement (shaded cycles in Fig. 4) was generally similar to the mean climate conditions (for more discussion, see Nugent et al., 2018).

Figure 4 Variability in (a) annual mean air temperature and (b) annual precipitation between 1994 and 2021 as recorded at Rivière-du-Loup (ECCC, 2023). The shaded circles indicate the measured period of the eddy covariance tower.

We made an extra simulation (in daily time steps) with 28-year climate input based on the 2013–2016 BDB setup to represent the normal-climate variability, also including extreme years. The simulated 28-year average of CO₂ uptake was $-\mathrm{101}\pm \mathrm{64}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, ranging from a maximum uptake of $-\mathrm{219}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in 1999 to a loss of $+\mathrm{54}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in 2015 (Fig. 5). At the annual scale, CO₂ uptake seems to increase slightly with increasing air temperature, although the relationship was not statistically significant (p=0.19). Annual CO₂ flux did not show a correlation with annual precipitation, but the model simulated the BDB peatland as an atmospheric CO₂ source for 3 years (1995, 1997, and 2015), all of which had below-average precipitation.

Figure 5 Accumulated annual CO₂ flux (a) simulated with the BDB 2013–2016 setup forced by Rivière-du-Loup 1994–2021 climate data and (b) measured over 2004–2018 at Mer Bleue (He et al., 2023a), a pristine peatland in the same climate region.

We further compared the simulated flux rates with long-term measurements at Mer Bleue, a pristine shrub–Sphagnum bog within the same climate region. Over 15 years of measurements (2004–2018), Mer Bleue had an average uptake rate of $-\mathrm{108}\pm \mathrm{33}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (simulated value of $-\mathrm{90}\pm \mathrm{51}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) (He et al., 2023b), similar to the 3-year BDB uptake rate measured by the tower ($-\mathrm{90}\pm \mathrm{18}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) and the current 28-year extended simulation ($-\mathrm{101}\pm \mathrm{64}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) (Fig. 6). Therefore, after 14 years of restoration, the BDB peatland switched back to C uptake, and the uptake rate was similar to pristine sites (for more discussion, see Nugent et al., 2018). However, our model simulations additionally show that the C uptake at the restored peatland has larger interannual variability (SD of $\mathrm{64}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) compared to the pristine Mer Bleue site. Under certain dry years, the restored site can potentially switch back to C sources, while the pristine peatlands showed persistent C uptake with a smaller interannual variation (SD of $\mathrm{33}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). In other words, the restored peatlands seem to have less ecological resilience compared to the pristine peatlands.

Figure 6 Comparison of CO₂ fluxes and emission factors from the different approaches.

3.4 Impact of future climate change on CO₂ uptake of restored peatlands

We evaluate the potential impact of future climate change on the CO₂ uptake function of the restored bog using the 28-year simulation as the long-term reference run. Climate change scenarios were designed as a combination of increases in year-round air temperature for the 28-year climate data by +1 °C and +2 °C. The scenarios also changed year-round precipitation by ±10 %, the range of climate change expected for this area of Quebec (Zhang et al., 2019). Then, equilibrium model runs using the 2013–2016 BDB setup for the future climate were conducted to evaluate the potential response of C-uptake functions.

Our model simulations show that increasing air temperature will decrease the CO₂ uptake rate of restored peatlands. Increasing air temperature alone by 1 °C decreases the annual C uptake rate by 5 % compared to the reference run and ∼16 % when air temperature increased by 2 °C (Fig. S3 in the Supplement). A larger rate of CO₂ uptake decreases under the +2 °C scenario compared to the +1 °C can be explained by the simulated more pronounced water table drop (Fig. S3). Our model simulation shows a change of ±10 % in precipitation alone only influences the CO₂ flux marginally, with a reduction in uptake rate when precipitation decreases (Fig. S4 in the Supplement). The BDB region is humid (annual precipitation / potential evaporation ratio is ∼1.5–2; Hare and Thomas, 1979). Thus, a 10 % change in precipitation is predicted to influence the water table marginally (Fig. S4). We made a climate scenario with an increase in air temperature of 2 °C and reduced precipitation by 10 % (i.e., the “extreme” scenario). The restored bog still acts as a C sink overall, with a slightly reduced ($\sim -\mathrm{6}\mathit{\%}$) simulated mean uptake rate of $-\mathrm{95}\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. 7). The modified climate causes both GPP and ER to increase (Fig. 7), thus effectively canceling each other out. Our model simulations thus overall suggest the restored bogs may retain some capacity for CO₂ uptake under future climate change.

Figure 7 Simulated CO₂ fluxes under a scenario in which air temperature is increased year-round by 2 °C and precipitation is decreased by 10 %. Equilibrium model runs used the BDB 2013–2016 setup and Rivière-du-Loup 1994–2021 climate data.

4 Discussion

The current model evaluation with the dataset from the BDB site shows that the CoupModel can simulate the coupled hydrology, heat, and CO₂ fluxes of a restored peatland. The CoupModel has been applied to Mer Bleue, a pristine bog (He et al., 2023b), and Rivière-du-Loup, an active peat extraction site (He et al., 2023a). The ability of the model to simulate C dynamics associated with ecohydrology for the restored system thus closes the land use cycle and shows the model can now simulate all stages of land uses, from pristine peatlands to lands drained for extraction and ultimately restoration.

Our model evaluation highlights the model deficiencies in simulating the time of phase changes in spring. This can be partly explained by the complex processes that occur during this period. For instance, the CoupModel probably overestimated the soil frozen depth as higher heat flow was partitioned into the soil surface over May to June every year (Figs. 1d and S1). Thus, extra heat was needed for thawing in the spring and delayed the increase in latent heat fluxes (Fig. S1), temperature increase (Fig. S2), and the start of spring CO₂ uptake (Fig. 2b). Moreover, the energy balance closure calculated with measured data over the 3 years is ∼90 %, while the CoupModel has full energy conservation. Thus, uncertainties in distributing surface energy fluxes can be carried over to the soil processes (e.g., soil surface heat flux). The lower model–data agreement for soil surface heat flux is due to its residual energy flux, thus being small in flux size (i.e., 1 order of magnitude lower compared to the turbulent energy fluxes; Figs. 1d and S1).

Our model performance for CO₂ flux is similar to previous models that have been applied to restored sites, such as the ECOSSE model (Premrov et al., 2021) and the PVN model (Lippmann et al., 2023). However, the advancement of our current modeling exercise compared to the earlier studies is its capability of accurately simulating both the water table depth and C dynamics at a finer temporal resolution. The CoupModel simulates the coupled C–hydrological processes at a half-hour resolution, while a daily time step was used for the earlier models. The ability to simulate processes at a subdaily scale is particularly important for the future inclusion of CH₄ as the transport processes (e.g., ebullition) occur at a subdaily scale (Walter and Heimann, 2000). Empirical studies have shown that the water table is an important control for greenhouse gas fluxes in restored peatlands (Evans et al., 2021; Järveoja et al., 2016; Koch et al., 2023). Restoration is associated with management practices that change the hydrology of the peatlands, such as blocking the drainage ditches at the beginning of restoration. With the gradual recovery of peat vegetation and the development of peat soil structure, the water table fluctuations are further reduced, and the mean level gradually moves above the collapse layer (Shantz and Price, 2006). Therefore, following restoration, the ecohydrology and vegetation co-evolve and feed back between each other, co-regulating the overall C-uptake function of the peatland. The ability of the CoupModel to simulate the coupled processes thus has important implications for understanding the overall climate impacts of peatland restorations. Our study simulates the time frame of 14–16 years after restoration, representing a stage of fully recovered vegetation. However, the degree of vegetation recovery might vary across sites. Further, it needs to be noted that the influence of nutrients and the altered pH levels that can encourage invasive species outcompeting mosses are not considered in our study. Future modeling research should include different peatland types and cover the beginning of the restoration, thus simulating the full dynamic coupling of vegetation development, hydrology management, and peat soil development.

The extended model simulations show that restored peatlands have less resilience to climate variability than do pristine peatlands (Figs. 5 and 6). Theoretical studies have argued that bogs are complex adaptive systems based on the tight feedbacks among plant production, decomposition, and water storage represented by water table depth (Eppinga et al., 2009; Belyea and Clymo, 2001). Water table frequency distribution can be a useful measure for evaluating the success of ecohydrology restoration (Shantz and Price, 2006). The CoupModel predicts that the water table frequency distribution for BDB will gradually recover to a state of a pristine bog when the newly grown mosses at the surface reach 40 cm depth (data not shown). Our model sensitivity analysis shows that mosses cannot thrive under a 100 cm acrotelm thickness, which is in agreement with results from field studies that suggest a tension of −100 mbar as the hydrologic threshold for Sphagnum establishment (Price, 1998; Price and Whitehead, 2001). The ability of the CoupModel to reproduce such important ecohydrology regulation has implications for future model applications to evaluate the impacts of field management practices on greenhouse gas fluxes by changing boundary and lateral hydrology conditions.

The current model exercise represents a series of studies towards developing the CoupModel as an IPCC Tier 3 methodology for estimating emissions from extracted and restored peatlands (He et al., 2023a; He and Roulet, 2023). Our work to date has focused on bogs in eastern Canada but should be expanded to include bogs and poor fens in western Canada and other geographical and ecoclimate regions in the future. To date, there are few emission data from restored peatlands in Canada, and those data are snapshots, covering only sites restored within 10 years, thus explaining the current EF (+2.07 t CO₂–C ${\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) used in the Canadian NIR. We argue this EF does not reflect the temporal dynamics of greenhouse gases for restored peatlands, particularly for those sites that have full vegetation recovery (Kalhori et al., 2024). Our 28-year extended model simulations, by considering the interannual climate variation, suggest an EF of $-\mathrm{1.01}\pm \mathrm{0.64}\mathrm{t}$ CO₂–C ${\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for a bog 14–16 years post restoration by the MLTT. This should be included in the next revision of EF within the NIR of Canada. Moreover, our modeled EF is ∼ 5 times larger (meaning more uptake) than the default IPCC Tier 1 EF for temperate nutrient-poor rewetted organic soils (−0.23 with CI −0.64–$+\mathrm{0.18}\mathrm{t}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) (IPCC, 2014). The data used to generate the IPCC EF include more degraded sites in Europe and different rewetting methods. The Canadian practice of leaving a residual peat layer at the end of extraction and using MLTT for restoration seems to be beneficial for the recovery of peatland C uptake. The default IPCC EF was previously used to evaluate the overall climate impacts on peatland restoration using a modeling approach (Gunther et al., 2020). Our results thus suggest those studies might significantly underestimate the climate cooling effects for Canadian bog sites that have been restored using MLTT.

Our climate change simulations show the regulating effect of air temperature on the CO₂ uptake of restored peatlands, where future global warming is predicted to moderately weaken the sink strength (Fig. S3). However, it should be noted that future changes in seasonal patterns and extremes were not accounted for in our climate change scenarios. Helbig et al. (2022) analyzed flux measurements from northern peatlands and showed earlier summer warming increases NEE uptake, while late summer warming decreases it. The seasonal patterns and particularly extremes of climate can be additional factors controlling the CO₂ fluxes. For mosses and peatland vegetation to develop, a stable water table is required. However, this can be challenging under the altered ecohydrological and climate conditions, especially in areas where drainage has lowered groundwater levels, causing less resilience of the ecosystem. Our results show that under dry years (e.g., 3 CO₂ source years: 1995, 1997, and 2015 for BDB), even temporary water level drops can lead to poor establishment and even die-off of mosses, reducing the restoration success and consequently the ecosystem uptake of CO₂. In addition, there is also the possibility of fire structurally altering the peatlands. Our simulations do not include fire, which is much less common in eastern Canadian peatlands than in the west (Zoltai et al., 1998; Lavoie and Pellerin, 2007). Thus, our climate change simulations probably overall represent a conservative prediction, which might in turn explain the moderate reduction in sink strength. As our extended simulations show, it is possible that over extreme years the site can switch to a small CO₂ source and that potentially the number of source years could increase in the future.

5 Conclusions

This study applied the CoupModel to a peatland site restored 14–16 years previously. We conclude the following:

• The CoupModel can describe the measured subdaily CO₂ fluxes, hydrology, and heat for the restored bog system.

• Restored bogs have less resilience to climate variability than pristine bog systems.

• The CoupModel simulation results in an emission factor of $-\mathrm{1.01}\pm \mathrm{0.64}\mathrm{t}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for Canadian bogs that have been restored for 14–16 years with the moss layer transfer technique, ∼ 5 times larger than the IPCC default emission factor and much smaller than the current emission factor used in the Canadian NIR.

• Moderate reduction in CO₂ uptake is predicted for restored bogs with full vegetation cover under future climate change conditions.

• The CoupModel now simulates all stages of peat extraction and restoration, and it can be used for exploring land use change issues, suggesting climate-smart management practices, and reporting Tier 3 emissions.