ORCHIDEE-PEAT (revision 4596), a model for northern peatland CO2, water, and energy fluxes on daily to annual scales

Abstract

Peatlands store substantial amounts of carbon and are vulnerable to climate change. We present a modified version of the Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) land surface model for simulating the hydrology, surface energy, and CO${}_{\mathrm{2}}$ fluxes of peatlands on daily to annual timescales. The model includes a separate soil tile in each 0.5${}^{\circ }$ grid cell, defined from a global peatland map and identified with peat-specific soil hydraulic properties. Runoff from non-peat vegetation within a grid cell containing a fraction of peat is routed to this peat soil tile, which maintains shallow water tables. The water table position separates oxic from anoxic decomposition. The model was evaluated against eddy-covariance (EC) observations from 30 northern peatland sites, with the maximum rate of carboxylation ($V_{\mathrm{cmax}})$ being optimized at each site. Regarding short-term day-to-day variations, the model performance was good for gross primary production (GPP) ($r^{\mathrm{2}}=$ 0.76; Nash–Sutcliffe modeling efficiency, MEF $=$ 0.76) and ecosystem respiration (ER, $r^{\mathrm{2}}=$ 0.78, MEF $=$ 0.75), with lesser accuracy for latent heat fluxes (LE, $r^{\mathrm{2}}=$ 0.42, MEF $=$ 0.14) and and net ecosystem CO${}_{\mathrm{2}}$ exchange (NEE, $r^{\mathrm{2}}=$ 0.38, MEF $=$ 0.26). Seasonal variations in GPP, ER, NEE, and energy fluxes on monthly scales showed moderate to high $r^{\mathrm{2}}$ values (0.57–0.86). For spatial across-site gradients of annual mean GPP, ER, NEE, and LE, $r^{\mathrm{2}}$ values of 0.93, 0.89, 0.27, and 0.71 were achieved, respectively. Water table (WT) variation was not well predicted ($r^{\mathrm{2}}$ < 0.1), likely due to the uncertain water input to the peat from surrounding areas. However, the poor performance of WT simulation did not greatly affect predictions of ER and NEE. We found a significant relationship between optimized $V_{\mathrm{cmax}}$ and latitude (temperature), which better reflects the spatial gradients of annual NEE than using an average $V_{\mathrm{cmax}}$ value.

Introduction

Peatlands cover only 3–5 % of the Earth's land area but store large amounts of soil organic carbon (SOC). This carbon is primarily located in the boreal and subarctic regions (75–80 %), while about 15 % is located in tropical regions (Frolking et al., 2011; Page et al., 2011). Current estimates of the northern peatland SOC vary from 270 to 450 Pg C (Gorham, 1991; Turunen et al., 2002; Yu et al., 2010). Northern peat accumulation occurred mainly during the Holocene, originating from plant litter production that exceeds decomposition in water-logged soil conditions, with low pH and low temperatures (Parish et al., 2008). The future of the carbon stored in these peatlands under a warmer environment and altered hydrological regimes is very uncertain. Logically, higher CO${}_{\mathrm{2}}$ concentrations and elevated temperatures will stimulate higher carbon uptake because of longer growing seasons and higher photosynthetic rates (Aurela et al., 2004; Adkinson et al., 2011). However, the accumulation is also coupled with a high evaporative demand that will lower the groundwater table, resulting in increased heterotrophic respiration rates (i.e., carbon loss; Mertens et al., 2001; Sulman et al., 2009; Adkinson et al., 2011). In addition to these potential climatic influences, other natural and anthropogenic disturbances (permafrost thaw, drainage, fires, etc.) further play a role in determining the future carbon balance of these vulnerable ecosystems (Turetsky et al., 2002; Parish et al., 2008). Drainage and fires have particularly important impacts on the carbon balance of the tropical peatlands (Page et al., 2002; Hooijer et al., 2010).

A number of peat carbon models have been reported in the literature. For example, Frolking et al. (2010) developed the Holocene Peat Model (HPM), which includes feedbacks between plant communities, water table, peat properties, and peat decomposition. This model was applied at Mer Bleue Bog in southern Canada and validated with data from peat-core observations. HPM is a long-term peat accumulation model that works at an annual time step but cannot simulate seasonal variations of key water processes in peatlands. Wania et al. (2009a, b) integrated peatlands and permafrost into the Lund–Potsdam–Jena model (LPJ-WHy), where the upper 0.3 m of peatland soils (the acrotelm) experience a fluctuating water table and the underlying layer (the catotelm) is permanently inundated. A constant soil moisture modifier (0.35) was used to reduce acrotelm decomposition. Spahni et al. (2013) adopted and improved LPJ-WHy by considering the effects of varying water table depth on acrotelm decomposition rates, using a weighted average of the aerobic and anaerobic respiration modifier, and implementation of a dynamic nitrogen cycle. In the dynamic global vegetation model (DGVM) CLIMBER2-LPJ, Kleinen et al. (2012) quantified the fraction of oxic decomposition in the acrotelm by comparing the water table position and the acrotelm height. Chaudhary et al. (2017a, b) included a dynamic multilayer peat accumulation functionality in a customized Arctic version of the Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS). In their approach, new layers of litter were added at the top of the soil every year, and the remaining litter mass, after decomposition, was treated as a new individual peat layer from the first day of the following year. The decomposition rate of peat, modulated by temperature and moisture, declined over time. In these four peatland models, the water table depth was calculated from a bucket model. In the context of Earth system modeling, the land surface processes are better represented by multilayer schemes, such as multilayer plant canopy and root, multilayer snow, multilevel soil carbon, and energy budgets (Best et al., 2011; Mcgrath et al., 2016; Zhu et al., 2016). To model peatlands consistently in land surface models, a multilayer soil hydrology scheme is needed. Meanwhile, a more physically based multilayer scheme can provide more prognostic power in predicting peatland water table dynamics.

In this study, we present the development of a multilayer peat hydrology and carbon model in the Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) land surface scheme, with a focus on the water table dynamics and its effects on the energy budgets, and on carbon decomposition occurring within the oxic and the water-saturated parts of the peat profile. CH${}_{\mathrm{4}}$ fluxes and DOC loss through runoff are important components of the carbon balance of a peatland (Chu et al., 2014; Olefeldt et al., 2012) but are not included in this study. This new peat model is incorporated consistently into the land surface scheme in order to conserve water, carbon, and energy at scales from local sites to grid-based large-scale applications in an Earth system modeling context.

Model description

General structure of the model

The ORCHIDEE land surface model simulates biophysical processes of rainfall interception, soil water transport, latent (LE) and sensible ($H$) heat fluxes, heat diffusion in the soil, and photosynthesis on a 30 min time step (Ducoudré et al., 1993). Carbon cycle processes (e.g., carbon allocation, respiration, mortality, litter, and soil carbon dynamics) are simulated on a daily time step (Krinner et al., 2005).

ORCHIDEE discretizes the vegetation into plant functional types (PFTs): eight for trees, two for natural C${}_{\mathrm{3}}$ and C${}_{\mathrm{4}}$ grasses, two for C${}_{\mathrm{3}}$ and C${}_{\mathrm{4}}$ crops, and one for bare-soil type. Across the PFTs, plants are described with the same equations but different parameter values, except for leaf onset and senescence that follow PFT-specific equations (Botta et al., 2000). In grid-based simulations, PFTs are grouped into three soil tiles: one with bare soil, one with all tree PFTs, and one with all short vegetation. The water budget of each soil tile is calculated independently. The version of ORCHIDEE implemented in this study uses the same (dominant) soil texture for all the soil tiles of a grid cell to define the reference saturated hydraulic conductivity ($K_{\text{s-ref}})$ and the saturated and residual volumetric water contents (${\mathit{\theta }}_{\mathrm{s}}$, ${\mathit{\theta }}_{\mathrm{r}})$. Dominant soil textural classes are taken from the Zobler's soil texture map (Zobler, 1986) at 1${}^{\circ }$ resolution. The original five soil textures (fine, medium-fine, medium, medium-coarse, coarse) in Zobler's map are reduced to three (fine, medium, coarse) by grouping the medium-fine, medium, and medium-coarse textures into a single class. Hydrological parameters of the three dominant soil textures are taken from Carsel and Parrish (1988) (Table 1).

Schematic of the hydrology module in ORCHIDEE. Water balance components: (a) a soil tile with either trees or grasses; (b) a peatland soil tile. Black dashed lines indicate the position of nodes in the 11 soil layers of the model. Blue lines: vertical profiles of saturated hydraulic conductivity for different soil textures. Green lines: diffusivity for different soil textures. The vertical axis indicates soil depth, the horizontal axis indicates values of saturated hydraulic conductivity ($K$, mm day${}^{-\mathrm{1}}$) and diffusivity ($D$, mm${}^{\mathrm{2}}$ day${}^{-\mathrm{1}}$). Note that the horizontal axis is on a base-10 logarithmic scale.

[Figure omitted. See PDF]

Each soil tile in ORCHIDEE has 11 vertical layers (up to 2.0 m) with exponentially coarser vertical resolution (Fig. 1). The Fokker–Planck equation is used to describe the vertical diffusion of water in the soil. The Mualem (1976) and Van Genuchten (1980) model (Eqs. 1 and 2) is used to define the hydraulic conductivity ($K$, m s${}^{-\mathrm{1}})$ and diffusivity ($D$, m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}})$ as a function of volumetric water content ($\mathit{\theta }$, m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$: $$\begin{array}{}& {\displaystyle }& {\displaystyle }K\left(\mathit{\theta }\right)=K_{\mathrm{s}}\sqrt{{\mathit{\theta }}_{\mathrm{f}}}{\left(\mathrm{1},-,{\left(\mathrm{1},-,{\mathit{\theta }}_{\mathrm{f}}^{\mathrm{1}/m}\right)}^m\right)}^{\mathrm{2}},& {\displaystyle }& {\displaystyle }D\left(\mathit{\theta }\right)={\displaystyle \frac{\left(\mathrm{1},-,m\right)K\left(\mathit{\theta }\right)}{\mathit{\alpha }m}}{\displaystyle \frac{\mathrm{1}}{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}}{\mathit{\theta }}_{\mathrm{f}}^{-\frac{\mathrm{1}}{m}}{\left({\mathit{\theta }}_{\mathrm{f}}^{-\frac{\mathrm{1}}{m}},-,\mathrm{1}\right)}^{-m},\end{array}$$ where $\mathit{\theta }$ is the volumetric water content (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$, ${\mathit{\theta }}_{\mathrm{s}}$ is the saturated water content (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$, ${\mathit{\theta }}_{\mathrm{r}}$ is the residual water content (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$, ${\mathit{\theta }}_{\mathrm{f}}$ is the relative water content and is calculated as ${\mathit{\theta }}_{\mathrm{f}}=\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}$, $K_{\mathrm{s}}$ is the saturated hydraulic conductivity (m s${}^{-\mathrm{1}})$, $\mathit{\alpha }$ is the inverse of the air entry suction (m${}^{-\mathrm{1}})$, and $m$ is a dimensionless parameter.

Van Genuchten parameters used for different soil texture classes for non-peat soils (coarse, medium, fine) and for peat. ${\mathit{\theta }}_{\mathrm{s}}$ is the saturated water content (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$; ${\mathit{\theta }}_{\mathrm{r}}$ is the residual water content (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$; $K_{\text{s-ref}}$ is the reference saturated hydraulic conductivity (m s${}^{-\mathrm{1}})$; $\mathit{\alpha }$ is the inverse of the air entry suction (m${}^{-\mathrm{1}})$; $n$ is a dimensionless parameter. In Eqs. (1) and (2), $m=\mathrm{1}-\mathrm{1}/n$.

	$K_{\text{s-ref}}$	$n$	$\mathit{\alpha }$	${\mathit{\theta }}_{\mathrm{s}}$	${\mathit{\theta }}_{\mathrm{r}}$
	(m s${}^{-\mathrm{1}})$		(m${}^{-\mathrm{1}})$	(m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$	(m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$
Coarse	1.23 $\times$ 10${}^{-\mathrm{5}}$	1.89	7.5	0.41	0.065
Medium	2.89 $\times$ 10${}^{-\mathrm{6}}$	1.56	3.6	0.43	0.078
Fine	7.22 $\times$ 10${}^{-\mathrm{7}}$	1.31	1.9	0.41	0.095
Peat	2.45 $\times$ 10${}^{-\mathrm{5}}$	1.38	5.07	0.90	0.15

Following d'Orgeval (2006) and d'Orgeval et al. (2008), $K_{\mathrm{s}}$ exponentially decreases with soil depth ($z$) below $z_{\lim }=$ 30 cm ($F_{\mathrm{d}}$ ($z$)), while a root-fracturing factor increases $K_{\mathrm{s}}$ where roots are denser ($F_{\mathrm{root}}$ ($z$)): $$K_{\mathrm{s}}\left(z\right)=K_{\mathrm{s}-\mathrm{ref}}\times F_{\mathrm{d}}\left(z\right)\times F_{\mathrm{root}}\left(z\right),$$ with $F_{\mathrm{d}}\left(z\right)=min\left(\max ,\left(\exp ,\left(-,f,\left(z,-,z_{\lim }\right)\right),,,\mathrm{0.1}\right),,,\mathrm{1}\right)$ and $F_{\mathrm{root}}\left(z\right)=\prod _{j\mathit{ϵ}c}\max {\left(\mathrm{1},,,{\left(\frac{K_{\mathrm{s}}^{\max }}{K_{\text{s-ref}}}\right)}^{\frac{\mathrm{1}-{\mathit{\alpha }}_{j}z}{\mathrm{2}}}\right)}^{{\mathrm{f}}_j}$, where $K_{\text{s-ref}}$ is the reference top-soil saturated hydraulic conductivity determined by soil texture (m s${}^{-\mathrm{1}})$, ${\mathrm{K}}_{\mathrm{s}}^{\max }$ is the value of the coarser (sandy) texture and equals 8.25 $\times$ 10${}^{-\mathrm{5}}$ m s${}^{-\mathrm{1}}$, ${\mathit{\alpha }}_j$ is a root profile decay factor for PFT $j$ with a coverage fraction $f_j$, and $c$ is the soil tile to which PFT $j$ was assigned.

Modifications in ORCHIDEE-PEAT

To simulate peat, we (1) modified the parameters of plants growing on peat, (2) added a new peat soil tile with specific peat soil hydraulic properties, and (3) changed the decomposition of peat carbon as being controlled by saturated conditions, through the modeled water table (WT).

Modified peat plant parameters

As a response to the unique stress conditions in peatlands (i.e., oxygen deficit, nutrient limitation), peatland vegetation has shallow and extensive root systems (Boutin and Keddy, 1993; Iversen et al., 2015). Previous peatland models have incorporated more than one PFT to represent peatland plants and dynamically simulate fractional vegetation cover. For example, Wania et al. (2009b) separated flood-tolerant C${}_{\mathrm{3}}$ graminoids and Sphagnum moss in LPJ-WHy to represent peatland-specific vegetation, with peatland extent defined from an organic soil map and the fractional cover of PFTs determined by bioclimatic conditions including temperature, water table depth, inundation stress, etc. Stocker et al. (2014) applied a version of this model but removed the upper temperature limitation of the peatland-specific PFTs and further included three additional PFTs: flood-tolerant C${}_{\mathrm{4}}$ grasses, tropical evergreen, and tropical raingreen tree PFTs, with peatland extent diagnosed by the TOPMODEL scheme. At present, however, ORCHIDEE-PEAT lacks representation of dynamic moss and shrub covers, and we do not know the fractional coverage of different vegetation types at each site in grid-based simulations. Previous studies have shown that there are considerable overlaps for the plant trait ranges among different plant functional types, while variations in plant traits within a PFT can be larger than the differences in means of different PFTs (Verheijen et al., 2013; Wright et al., 2005; Laughlin et al., 2010). Therefore, for simplicity, we applied only the PFT of C${}_{\mathrm{3}}$ grass with a shallower rooting depth to represent the average of vegetation growing in northern peatlands.

Only one key photosynthetic parameter ($V_{\mathrm{cmax}}$) of this PFT has been tuned to match with observations at each site. This simplification may cause discrepancies between model output and observations. Druel et al. (2017) added non-vascular plants (bryophytes and lichens), boreal grasses, and shrubs into ORC-HL-VEGv1.0. Their work is parallel to our model and will be incorporated into the model in the future. It will then be possible to verify how many plant functional types are needed by the model to reliably simulate the peatlands at site level and larger scale. The maximum rate of carboxylation ($V_{\mathrm{cmax}})$ typically varies across peat sites (Rennermalm et al., 2005; Bubier et al., 2011) and further varies with leaf nitrogen, phosphorus content, and specific leaf area (Wright et al., 2004; Walker et al., 2014). For instance, $V_{\mathrm{cmax}}$ for Sphagnum at the Old Black Spruce site (53.985${}^{\circ }$ N, 105.12${}^{\circ }$ W) in Canada was 5, 14, and 6 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ during spring, summer, and autumn, respectively, while $V_{\mathrm{cmax}}$ for Pleurozium was 7, 5, and 7 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ during the three seasons (Williams and Flanagan, 1998). Bui (2013) conducted a fertilization experiment at the Mer Bleue Bog (Canada; 45.41${}^{\circ }$ N, 75.52${}^{\circ }$ W) on the dominant ericaceous shrub and reported that $V_{\mathrm{cmax}}$ values ranged between 6 and 179 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, with significantly higher $V_{\mathrm{cmax}}$ values after addition of nitrogen (6.4 g N m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$ at 20 times the growing season ambient wet N deposition rate with or without phosphorus (P) and potassium (K). In this study (Sect. 4.1), we calibrated $V_{\mathrm{cmax}}$ at each site so that modeled peak gross primary production (GPP) matched peak values derived from direct EC measurements, and then regressed this adjusted $V_{\mathrm{cmax}}$ value with environmental and climate variables. We note that this adjustment of $V_{\mathrm{cmax}}$ may over- or undercompensate for biases in other model parameters that impact maximum GPP, such as leaf area index (LAI), specific leaf area (SLA), canopy light absorption parameters, water, and temperature stresses (Fig. S1 in the Supplement).

Peat-specific soils hydraulics

Peatlands generally occur in flat areas that are poorly drained and/or receive runoff and subsurface water from the surrounding landscape (Graniero and Price, 1999). The low permeability catotelm peat layer is permanently saturated. In ORCHIDEE-PEAT, the new soil tile added in a grid cell to represent peatland as a landscape element was assumed to receive surface runoff from the other three soil tiles (bare soil, trees, grasses) and has a drainage flux reduced to zero (Largeron et al., 2017). Further, considering that the water table of a peatland can rise above the ground surface, an above-surface water reservoir with a maximum height of 10 cm was added (Fig. 1b). In the model, the partitioning between water infiltration and surface runoff is computed through a time-splitting procedure, with the maximum infiltration rates described as an exponential probability density distribution (d'Orgeval, 2006). The infiltration-excess water of peatland first fills the above-surface water reservoir, then leaves the grid cell as runoff. Water in this above-surface reservoir re-infiltrates into the peat soil on the next time step (Largeron et al., 2017). We verified that the measured standing water remained below 10 cm above the soil surface at 16 out of 20 northern peat sites where water table depth was recorded in this study (Table S1 in the Supplement). The four exceptions were Winous Point North Marsh (US-WPT), Himmelmoor (DE-Hmm), an Alaska fen (US-Fen), and an Alaskan bog (US-Bog), where observed water tables reached up to 77, 39, 46, and 34 cm above the soil surface, respectively.

Peat soils cannot be described with any of the mineral soil textures used for other tiles (Table 1) because the low bulk density and high porosity increase the downward water percolation (Rezanezhad et al., 2016). Observed peat-saturated hydraulic conductivity ($K$) and diffusivity ($D$) strongly vary in space, depth, and time. This is partly related to the degree of decomposition and compression of organic matter (Gnatowski et al., 2010). Morris et al. (2015) reported near-surface saturated hydraulic conductivities ($K$) of 2.69 $\times$ 10${}^{-\mathrm{2}}$ to 7.16 $\times$ 10${}^{-\mathrm{6}}$ m s${}^{-\mathrm{1}}$ in bogs. Gnatowski et al. (2010) measured values of 5 $\times$ 10${}^{-\mathrm{6}}$ m s${}^{-\mathrm{1}}$ in a moss-covered peat, which was 2 orders of magnitude larger than for a woody peat (5.56 $\times$ 10${}^{-\mathrm{8}}$ m s${}^{-\mathrm{1}})$. Peat hydraulic parameters values used in this study were applied after Largeron et al. (2017), based on Letts et al. (2000) and Dawson (2006) (Table 1). The peat-saturated hydraulic conductivity value of 2.45 $\times$ 10${}^{-\mathrm{5}}$ m s${}^{-\mathrm{1}}$ is comparable to the harmonic mean value (6 $\times$ 10${}^{-\mathrm{5}}$ m s${}^{-\mathrm{1}})$ of Morris et al. (2015). The values of the other Van Genuchten parameters for peat (Table 1) are similar to those employed in other peatland models (Wania et al., 2009a; Wu et al., 2016).

The peatland water table depth (cm) is diagnosed by summing water heights in the 11 soil layers, calculated from the relative water content (Largeron et al., 2017): $$\begin{array}{cc}{\displaystyle }\mathrm{WT}& {\displaystyle }=H_{\mathrm{tot}}-{\sum }_{i=\mathrm{1}}^{\mathrm{11}}\left({\mathit{\theta }}_{fi},\times ,{\mathrm{dz}}_i\right)-H_{\mathrm{ab}},\\ & {\displaystyle }& {\displaystyle }\mathrm{with}{\mathit{\theta }}_{fi}={\displaystyle \frac{{\mathit{\theta }}_i-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}},\end{array}$$ where ${\mathit{\theta }}_{fi}$ is the relative volumetric water content of the $i$th soil layer, ${\mathit{\theta }}_{\mathrm{s}}$ is the saturated water content (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$, ${\mathit{\theta }}_{\mathrm{r}}$ is the residual water content (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}})$, dz${}_i$ is the distance between node $i-\mathrm{1}$ and node $i$ (Fig. 1; m), $H_{\mathrm{tot}}$ is the total soil column height being fixed to 2.0 m, and $H_{\mathrm{ab}}$ is the height of the water reservoir above soil surface (m). Thus, when the water table is above the surface, the modeled WT takes negative values.

Decomposition of peat carbon controlled by water saturation

In the standard version of ORCHIDEE, plant litter carbon is added to two litter pools: the metabolic and the structural pool. Decomposed litter carbon from these two pools is then distributed into three soil carbon pools: the active, slow, and passive pools, similar to the CENTURY model (Parton et al., 1988). Both temperature and moisture functions are used to control soil carbon decomposition rates (Text S1 in the Supplement). In ORCHIDEE-PEAT, these standard processes are kept the same as in Krinner et al. (2005) for non-peatland vegetation (Fig. S2, black dashed box). For the peatland vegetation, we added a peat carbon module, in which the three soil carbon pools (active, slow, and passive) are replaced by two pools forming distinct layers, following Kleinen et al. (2012) (Fig. S2, red dashed box). Specifically, carbon from decomposed litter pools is added to the acrotelm carbon pool where it is decomposed aerobically above the simulated water table and anaerobically below it. The permanently saturated deep catotelm carbon pool receives a prescribed fraction of the acrotelm carbon, and is decomposed only anaerobically at a very slow rate. While the acrotelm depth is fixed to 30 cm in some peat decomposition models (Yurova et al., 2007; Wania et al., 2009a; Spahni et al., 2013), we used the average of simulated minimum summer water table position (WT${}_{\min })$ over the observational period to demarcate the boundary between the acrotelm and the catotelm at each site to take into account local site conditions. We conducted a “preparation run (S0)”, in which the model was run at each site using the same protocol (Sect. 3.3) but with the peat carbon module deactivated. WT${}_{\min }$ was diagnosed from the output of S0 before feeding into the peat carbon module in S1 and S2 (Sect. 3.3). Soil carbon exerts no feedback effects on the soil temperature and hydraulic in the structure of our model; thus, S0 and S1 produce the same simulated water table. WT${}_{\min }$ values were estimated based on current climate due to the lack of knowledge of initiation histories of these sites. For the long-term carbon accumulation estimations, the Holocene climate may be a better proxy since northern peatlands show peak initiation in the early Holocene (Yu et al., 2010). By comparing the height of the acrotelm (Fig. S2, Eq. 9) with the WT depth, we derived the fraction of the acrotelm where carbon decomposes under oxic ($\mathit{\beta })$ vs. anoxic conditions (1$\mathit{\beta })$. Acrotelm height ($H_{\mathrm{A}}$, Eq. 10) was calculated from acrotelm carbon stock (C${}_{\mathrm{A}}$ in Eqs. 5–7), acrotelm carbon fraction (C${}_{\mathrm{f},\mathrm{A}})$ and acrotelm bulk density (${\mathit{\rho }}_{\mathrm{A}})$. Decomposition of peat carbon is controlled by temperature ($f_{\mathrm{T}})$ and parameterized as an exponential function: $f_{\mathrm{T}}=Q_{\mathrm{10}}$ exp(($T-T_{\mathrm{ref}})/$10 ${}^{\circ }$C) with $Q_{\mathrm{10}}=$ 2.0 and $T_{\mathrm{ref}}=$ 30 ${}^{\circ }$C (Text S1). Soil carbon fluxes are given by $$\begin{array}{}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{AC}}=k_{\mathrm{p}}{f}_{\mathrm{T}}{\mathrm{C}}_{\mathrm{A}},& {\displaystyle }& {\displaystyle }{R}_{\mathrm{A},\mathrm{o}}=\mathit{\beta }{k}_{\mathrm{A}}{f}_{\mathrm{T}}{\mathrm{C}}_{\mathrm{A}},& {\displaystyle }& {\displaystyle }{R}_{\mathrm{A},\mathrm{a}}=(\mathrm{1}-\mathit{\beta }){\mathrm{vk}}_{\mathrm{A}}{f}_{\mathrm{T}}{\mathrm{C}}_{\mathrm{A}},& {\displaystyle }& {\displaystyle }{R}_{\mathrm{C}}=k_{\mathrm{C}}{f}_{\mathrm{T}}{\mathrm{C}}_{\mathrm{C}},& {\displaystyle }& {\displaystyle }\mathit{\beta }=\left\{\begin{array}{l}\mathit{\beta }=\mathrm{1},{\mathrm{WT}}_{\min }-\mathrm{WT}\le \mathrm{0}\\ \mathit{\beta }={\displaystyle \frac{H_{\mathrm{A}}-({\mathrm{WT}}_{\min }-\mathrm{WT})}{H_{\mathrm{A}}}},\mathrm{0}<{\mathrm{WT}}_{\min }-\mathrm{WT}<H_{\mathrm{A}}\\ \mathit{\beta }=\mathrm{0},{\mathrm{WT}}_{\min }-\mathrm{WT}\ge H_{\mathrm{A}}\end{array}\right.,& {\displaystyle }& {\displaystyle }{H}_{\mathrm{A}}={\displaystyle \frac{{\mathrm{C}}_{\mathrm{A}}}{{\mathit{\rho }}_{\mathrm{A}}\cdot {\mathrm{C}}_{\mathrm{f},\mathrm{A}}}},\end{array}$$ where $F_{\mathrm{AC}}$ is the carbon flux from acrotelm to catotelm; $R_{\mathrm{A},\mathrm{o}}$ is aerobically decomposed acrotelm carbon; $R_{\mathrm{A},\mathrm{a}}$ is anaerobically decomposed acrotelm carbon; $R_{\mathrm{C}}$ is decomposed carbon in catotelm; C${}_{\mathrm{A}}$ is carbon stored in the acrotelm; C${}_{\mathrm{C}}$ is carbon stored in the catotelm; and $\mathit{\beta }$ is the fraction of acrotelm under oxic conditions. A 10 100-year spinup was conducted to initialize peat depth at each site (Sect. 3.3). Following the study of Kleinen et al. (2012), the catotelm formation rate $k_{\mathrm{p}}=$ 1.91 $\times$ 10${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$, the acrotelm decomposition rate $k_{\mathrm{A}}=$ 0.067 yr${}^{-\mathrm{1}}$, the catotelm decomposition rate $k_{\mathrm{C}}=$ 3.35 $\times$ 10${}^{-\mathrm{5}}$ yr${}^{-\mathrm{1}}$, the ratio of anaerobic to aerobic CO${}_{\mathrm{2}}$ production $\mathit{\mu }=$ 0.35, carbon fraction in the acrotelm peat C${}_{\mathrm{f},\mathrm{A}}=$ 0.50, the acrotelm density ${\mathit{\rho }}_{\mathrm{A}}=$ 35.0 kg m${}^{-\mathrm{3}}$, carbon fraction in the catotelm peat C${}_{\mathrm{f},\mathrm{C}}=$ 0.52, and the catotelm density ${\mathit{\rho }}_{\mathrm{C}}=$ 91.0 kg m${}^{-\mathrm{3}}$.

The distribution of 30 peatland sites used in this study. Triangles are bogs; circles are fens; squares are tundra and marsh. Colors of the markers indicate peatland fractions in the 0.5${}^{\circ }$ grid cell. Mean air temperatures is the annual mean from 1999 to 2015, based on the 6-hourly CRU-NCEP 0.5${}^{\circ }$ global database.

[Figure omitted. See PDF]

Site characteristics of the 30 peatlands (sites are sorted by latitude from south to north). The first column denotes if the site is used in the second set of simulation (S2, with water table prescribed in the model equal to observed values): y – YES; n – NO. Lat: latitude; Long: longitude; MAT: long-term mean annual air temperature; MAP: long-term mean annual precipitation; peatland fraction (%): fraction of peatland in the 0.5${}^{\circ }$ grid cell which is read from the map of Yu et al. (2010); for cells where there is no peatland, mean fraction (22 %) is used. Note that, at US-Bog and US-Fen, the precipitation is the growing season (from 16 May to 31 August) mean value, thus clarified as “GS” in the table. Details of S2 and peatland fraction are provided in Sect. 3.3.

S2	Code	Lat	Long	Climatic	Type	MAP	MAT	Elevation	Peatland	Period	Dominant	LAI	Aboveground	Citation
				zone		(mm)	(${}^{\circ }$C)	(m a.s.l.)	fraction		vegetation	(m${}^{\mathrm{2}}$ m${}^{-\mathrm{2}})$	biomass
											type		(kg m${}^{-\mathrm{2}})$
n	US-WPT	41.5	$-$83.0	temperate	marsh	840	9.2	175	mean	2011–2013	grasses	area average: 2.3;	area average: 1.94;	Chu et al. (2014,
												emergent ;	emergent vegetation	2015)
												vegetation: 3.3	area: 3.04;
												open water: 1.0	open water area: 0.44
n	CA-Mer	45.4	$-$75.5	temperate	bog	944	6	70	mean	1999–2012	shrubs, mosses	1.5	moss: 0.144 $\pm$ 0.03;	Lafleur et al. (2005)
													vascular: 0.356 $\pm$ 0.1
y	US-Los	46.1	$-$90.0	temperate	fen	666	3.8	470	mean	2000–2010	trees, shrubs,	4.24	1.336	Sulman et al. (2009)
											grasses
n	LA-GUE	47.3	2.3	temperate	fen	880	11	145	mean	2011–2013	grasses			D'Angelo et al. (2016);
														Laggoun-Défarge et al. (2016)
y	DE-Sfn	47.8	11.3	temperate	bog	1127	8.6	590	3.01 %	2012–2014	trees, shrubs,			Hommeltenberg et al. (2014)
											grasses, mosses
y	CZ-Wet	49.0	14.8	temperate	fen	614	7.4	426.5	mean	2007–2013	grasses	2.45	0.57	Dušek et al. (2009)
n	DE-Spw	51.9	14.0	temperate	fen	559	9.5	61	11.01 %	2010–2014	trees	3.6		Petrescu et al. (2015)
y	IE-Kil	52.0	$-$9.9	temperate	blanket	2467	10.5	150	28.97 %	2002–2012	shrubs, grasses,	from 0.4 to 0.6 in		Sottocornola et al. (2009);
					bog						mosses	different years		McVeigh et al. (2014)
y	DE-Bou	52.7	7.2	temperate	bog	799	10	19	63.98 %	2011–2014	grasses,	0.7	grass dominated:	Hurkuck et al. (2016)
											mosses		0.577 $\pm$ 0.029;
													heather and moss
													dominated:
													0.517.0 $\pm$ 0.026;
													mixed:
													0.303 $\pm$ 0.015
n	PL-Wet	52.5	16.2	temperate	fen	526	8.5	54	4.01 %	2006–2013	shrubs, grasses,			Chojnicki et al. (2007);
											mosses			Barabach (2012);
														Milecka et al. (2017)
n	PL-Kpt	53.6	22.9	temperate	fen	600	7.1	109	mean	2013–2015	grasses, reeds,	sedges: 4.3;		Fortuniak et al. (2017)
											and ferns	reeds and ferns: 4.8
n	DE-Hmm	53.7	9.9	temperate	bog	838	9	12	15.99 %	2012–2014	90 % bare peat,			Vanselow-Algan et al. (2015)
											10 % vegetation
											cover: trees, grasses
n	DE-Zrk	53.9	12.9	temperate	fen	584	8.7	< 0.5	23.16 %	2013–2014	grasses			Franz et al. (2016)
n	CA-Wp3	54.5	$-$113.3	boreal	fen	504	2.1	670	29.77 %	2004–2006	grasses, mosses	1.1	0.157	Adkinson et al. (2011)
n	CA-Wp1	55.0	$-$112.5	boreal	fen	504	2.1	540	0.20 %	2003–2009	trees, shrubs,	2.6	1.08	Flanagan and Syed (2011)
											mosses
n	CA-Wp2	55.5	$-$112.3	boreal	fen	504	2.1	730	8.07 %	2004–2006	shrubs, grasses,	1.5	0.231	Adkinson et al. (2011)
											mosses
y	SE-Faj	56.3	13.6	temperate	bog	700	6.2	140	mean	2005–2009	shrubs, grasses,		dwarf shrub: 0.153;	Lund et al. (2007, 2012)
											mosses		Sphagnum: 0.192;
													graminoid: 0.077
n	FI-Sii	61.8	24.2	boreal	fen	713	3.3	162	mean	2005–2014	shrubs, grasses,	0.55		Aurela et al. (2007);
											mosses	(maximum value,		Riutta et al. (2007)
												occurs in June–July)
n	DK-NuF	64.1	$-$51.4	arctic	fen	750	$-$1.4	40	mean	2008–2014	grasses, mosses	0.7		Westergaard-Nielsen
														et al. (2013)
y	SE-Deg	64.2	19.6	boreal	fen	523	1.2	270	mean	2001–2005	shrubs, grasses,	0.47	moss: 0.065;	Sagerfors et al. (2008);
											mosses		vascular: 0.049	Nilsson et al. (2008);
														Peichl et al. (2014)
n	US-Bog	64.7	$-$148.3	boreal,	bog	146 (GS)	$-$2.2	100	28.01 %	2011–2015	trees, mosses			Euskirchen et al. (2014)
				thermokarst
n	US-Fen	64.7	$-$148.3	boreal	fen	146 (GS)	$-$2.2	100	28.01 %	2011–2015	grasses,			Euskirchen et al. (2014)
											forbs
y	FI-Lom	68.0	24.2	boreal	fen	521	$-$1	269	5.08 %	2007–2009	shrubs, grasses,	1.3		Aurela et al. (2009)
											mosses
n	SE-Sto	68.4	19.1	boreal,	bog	322	$-$0.14	360	mean	2014–2015	shrubs, grasses,			Malmer et al. (2005);
				permafrost							mosses			Olefeldt et al. (2012)
n	US-Ics	68.6	$-$149.3	arctic,	fen	318	$-$7.4	920	mean	2007–2011	shrubs,			Euskirchen et al. (2012,
				permafrost							grasses			2016)
n	RU-Che	68.6	161.3	arctic,	tundra	200–215	$-$12.5	4	64.09 %	2002–2005	shrubs,	0.3–0.4		Corradi et al. (2005);
				permafrost							grasses			Merbold et al. (2009)
n	NO-And	69.1	16.0	boreal	bog	1060	3.6	17	mean	2008–2014	shrubs, grasses,			Lund et al. (2015)
											mosses
n	US-Bes	71.3	$-$156.6	arctic,	tundra	173	$-$12	4	mean	2005–2008	grasses,			Zona et al. (2009)
				permafrost							mosses
n	DK-ZaF	74.5	$-$20.6	arctic,	fen	211	$-$9	35	mean	2008–2011	grasses,	0.65	0.471	Stiegler et al. (2016)
				permafrost							mosses
n	NO-Adv	78.2	15.9	arctic,	fen	190	$-$6.7	17	mean	2011–2014	shrubs, grasses,	0.41 $\pm$ 0.12	0.85 $\pm$ 0.28	Pirk et al. (2017)
				permafrost							mosses

${}^{*}$ For most of the sites, NEE was partitioned into GPP and ecosystem respiration following the nighttime partitioning method of Reichstein et al. (2005), except that NO-And used a light response curve approach following Lund et al. (2015); CA-Wp1 used the FLUXNET Canada Research Network (FCRN) standard NEE partitioning procedure following Barr et al. (2004); and DE-Spw used the online gap filling and flux partitioning tool (http://www.bgc-jena.mpg.de/~MDIwork/eddyproc/) which uses the method proposed by Lloyd and Taylor (1994). Note that the we grouped sedges, grasses, and herbaceous plants into one class (grasses) in the table.

Validation of ORCHIDEE-PEAT at Northern Hemisphere peatland eddy-covariance sites

Sites description

To evaluate the performance of ORCHIDEE-PEAT in simulating CO${}_{\mathrm{2}}$, water, and energy fluxes on daily to annual timescales, we compiled data from 30 northern peatland sites where eddy-covariance data and physical variables (water table, snow depth, soil temperature) were collected (Fig. 2, Table 2). These sites are spread between the temperate and the arctic climate zones, and include nine bogs and 18 fens. A marsh and two wet tundra sites (note that these two wet tundra sites are neither a fen nor a bog; hereafter, they are referred to as “tundra”) with a $\sim$ 30–50 cm thick organic layer are also included in this study. Among them, six sites are underlain by permafrost and one site is in a thermokarst area. The peatland fractional cover in the 0.5${}^{\circ }$ grid cell containing each site is from the Yu et al. (2010) map (Fig. 2, Table 2). A short description of all sites can be found in the Supplement.

Meteorological forcing data

We ran the model for 30 different 0.5${}^{\circ }$ grid cells corresponding to each peatland site (US-Fen and US-Bog are in the same grid cell but their local meteorological data were different). Peatland fraction in each grid cell was prescribed from Yu et al. (2010), adapted by Largeron et al. (2017) to be matched with a high-resolution land cover map. For the 16 out of 30 cells without peatland (Fig. 2, Table 2) in the large-scale map from Yu et al. (2010), a mean peatland fraction of 22 % was assigned.

Time series of half-hourly air temperature, wind speed, wind direction, longwave incoming radiation, shortwave incoming radiation, specific humidity, atmospheric pressure, and precipitation were used to drive ORCHIDEE-PEAT. All variables were from measurements made at each flux tower where CO${}_{\mathrm{2}}$ and energy (latent heat: LE; sensible heat: $H$) fluxes, water table position, soil temperature, and snow depth were recorded on a half-hourly time step. The linearly interpolated 6-hourly CRU-NCEP 0.5${}^{\circ }$ global climate forcing dataset was used to fill the gaps in the driving variables. A linear correction was applied to meteorological forcing variables (except precipitation) in the CRU-NCEP dataset to match observations before gap filling. For precipitation, no correction was applied. At CA-Wp2 and CA-Wp3, meteorological forcing data were measured only during the growing season, so CRU-NCEP data were linearly corrected using relationships derived from the available data. For some sites, several meteorological variables were not measured, such as longwave incoming radiation at NO-And, atmospheric pressure, shortwave incoming radiation, and longwave incoming radiation at CZ-Wet. In these cases, uncorrected CRU-NCEP data were used.

Model setup

ORCHIDEE-PEAT was first spun up for 10 100 years, forced by the pre-industrial atmospheric CO${}_{\mathrm{2}}$ concentration of 285 ppm, with repeated site-specific observational meteorological fields, and present-day vegetation fractions for each site. In reality, the climate changed through the Holocene, but since the initiation and climate history of each site are unknown, we assumed a constant present-day climate condition and peatland area. Thus, this model is only suitable for simulating water, energy and CO${}_{\mathrm{2}}$ fluxes from peat on timescales ranging from days to decades. To accelerate the spinup, ORCHIDEE-PEAT was first run for 100 years to reach the equilibrium for hydrology and soil thermal conditions, fast carbon pools, and soil carbon input from dead plants. Then, a submodel simulating only soil carbon dynamics (with fixed daily litter input from the previous simulation) was run for 10 000 years to accumulate soil carbon. Peatlands can reach equilibrium only when the addition of carbon equals carbon lost, which is attained on timescales of 10${}^{\mathrm{4}}$ years (Clymo, 1984; Wania et al., 2009b). The catotelm carbon pool in this study was still not fully equilibrated even after 10 100 years due to the low carbon decomposition rate in this reservoir (3.35 $\times$ 10${}^{-\mathrm{5}}$ yr${}^{-\mathrm{1}}$, Kleinen et al., 2012). The modeled peat carbon pool thus depends on the time length of spinup, which was fixed at 10 100 years, while in the real world, peat age at some sites can be younger. For example, the sample from the second last 10 cm peat segment at CA-Wp1 has an uncalibrated radiocarbon date of $\sim$ 2200 years (Flanagan and Syed, 2011). Since we focus on carbon and water fluxes on daily to annual scales in this study, rather than on the simulation of peat carbon stocks, we conducted a sensitivity analysis of modeled heterotrophic respiration to the length of the spinup, which shows only a slight increase of catotelm respiration with increasing simulation time (Fig. S3). After the spinup, transient simulations were conducted for each site, forced by repeated site-specific climates and rising atmospheric CO${}_{\mathrm{2}}$ concentration during the period 1901–2015. Finally, the model outputs corresponding to the respective measurement periods (all during 1999–2015) were compared to observed time series for each site.

Two sets of simulations were conducted. In the first one (S1), soil water content and WT position were modeled by ORCHIDEE-PEAT, and the WT was used in the carbon module to define the fraction of oxic and anoxic decomposition in the acrotelm. S1 was performed for all the 30 sites. In the second set (S2) of simulations, we prescribed water table in the model to equal the observed values (WT${}_{\mathrm{obs}})$. That is, soil moisture at layers below the measured water table was prescribed as saturated ($\mathit{\theta }(z$ > WT${}_{\mathrm{obs}})$ =${\mathit{\theta }}_{\mathrm{s}})$, while soil moisture above WT${}_{\mathrm{obs}}$ was simulated. WT${}_{\mathrm{obs}}$ was further used in the carbon module in S2. S2 was performed only for a subset of eight sites where at least 2 years of water table measurements were available and where there were sufficient observations to gap fill the WT${}_{\mathrm{obs}}$ time series (Table 2). For these sites, the gaps of WT${}_{\mathrm{obs}}$ were filled with the mean value of the same period from other years of measurement (Table S2). The simulation S2 was designed to check if the model performance will improve (or deteriorate) when prescribing WT exactly to its observed value, since WT is known to be a critical variable impacting peat water, CO${}_{\mathrm{2}}$, and CH4 fluxes (Dušek et al., 2009; Parmentier et al., 2011; Strack et al., 2006). Fixing the simulated water table to WT${}_{\mathrm{obs}}$ in S2 violated the water mass conservation of the model but allowed us to evaluate the carbon module independently from the hydrological module biases.

Measures for evaluating model performance

Following Jung et al. (2011) and Tramontana et al. (2016), we used site-specific daily means, annual means, seasonal variations, and daily anomalies to evaluate the model performance. For each site, seasonal variations are calculated by removing the annual mean value from the mean seasonal cycle (averaged value for each month across all available years). Anomalies are calculated as the deviation of a daily flux value from the corresponding mean seasonal cycle.

A series of measures was used to assess the model performance (Kobayashi and Salam, 2000; Jung et al., 2011; Tramontana et al., 2016).

The root mean square deviation (RMSD) reports the model accuracy by measuring the differences between simulation and observation. $$\mathrm{RMSD}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(x_i-y_i{)}^{\mathrm{2}}},$$ where $x_i$ is simulated variable, $y_i$ is measured variable, and $n$ is the number of observations.

Two signals (SDSD and LCS) are discriminated from the mean squared deviation (Kobayashi and Salam, 2000). The squared difference (SDSD) between the standard deviation of the simulation (SD${}_{\mathrm{s}})$ and the measurement (SD${}_{\mathrm{m}})$ shows if the model can reproduce the magnitude of fluctuation among the $n$ measurements. $$\begin{array}{cc}{\displaystyle }\mathrm{SDSD}& {\displaystyle }=({\mathrm{SD}}_{\mathrm{s}}-{\mathrm{SD}}_{\mathrm{m}}{)}^{\mathrm{2}};\\ {\displaystyle }\mathrm{with}{\mathrm{SD}}_{\mathrm{s}}& {\displaystyle }=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(x_i-\overline{x}{)}^{\mathrm{2}}},\\ & {\displaystyle }{\mathrm{SD}}_{\mathrm{m}}& {\displaystyle }=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(y_i-\overline{y}{)}^{\mathrm{2}}},\end{array}$$ where $\overline{x}$ is simulated mean value; $\overline{y}$ is measured mean value.

The lack of correlation weighted by the standard deviations (LCS) is a measure to examine if the model reproduces the observed phase of variability. $$\begin{array}{cc}{\displaystyle }\mathrm{LCS}& {\displaystyle }=\mathrm{2}{\mathrm{SD}}_{\mathrm{s}}{\mathrm{SD}}_{\mathrm{m}}\left(\mathrm{1},-,r\right);\\ & {\displaystyle }& {\displaystyle }\mathrm{with}r=\left[{\displaystyle \frac{\mathrm{1}}{n}},\sum _{i=\mathrm{1}}^n,\left(x_i,-,\overline{x}\right),\left(y_i,-,\overline{y}\right)\right]/\left({\mathrm{SD}}_{\mathrm{s}}{\mathrm{SD}}_{\mathrm{m}}\right),\end{array}$$ where $r$ is Pearson's correlation coefficient.

The Nash–Sutcliffe modeling efficiency (MEF) is used to indicate the predictive accuracy of the model. MEF varies between negative infinity ($-$inf) and 1: an efficiency of 1 indicates a perfect fit between simulations and observations; an efficiency of 0 indicates the simulations are as accurate as the mean value of observations; a negative MEF indicates that the mean value of observations has greater predictive power than the model. The modeling efficiency is defined as $$\mathrm{MEF}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^n(x_i-y_i{)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n(y_i-\overline{y}{)}^{\mathrm{2}}}}.$$

Results

Site-specific $V_{\mathrm{cmax}}$ reduces errors in carbon flux simulations

Out of the 30 sites, 22 sites provided observed daily GPP (based on measured NEE). The values of optimized $V_{\mathrm{cmax}}$ at each site were listed in Table 3. The optimized $V_{\mathrm{cmax}}$ varied from 19 to 89 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ (Table 3), with a mean value of 40 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. The calibration of $V_{\mathrm{cmax}}$ may compensate for biases in other model parameters. A brief comparison between simulated and reported (measured/estimated) LAI and aboveground biomass showed that there are no systematic errors (Fig. S1).

Taylor diagrams were used to evaluate model results at these 22 sites (Fig. 3). The model had the best performance for GPP, with the correlation coefficient between simulated and observed GPP varying between 0.66 and 0.93, and all data points fell within the 0.9 root mean square difference circle. Simulated water table depth had a larger spread in correlation (0.16–0.82) and root mean square difference (0.4–4.0). We found no significant patterns of model–data misfits among different peatland types (fen, bog, and others) or climate zones (temperate, boreal, and arctic; Fig. 3).

Optimized $V_{\mathrm{cmax}}$ ($\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$ at each site.

Site	$V_{\mathrm{cmax}}$	Site	$V_{\mathrm{cmax}}$
US-WPT	80	FI-Sii	19
CA-Mer	25	DK-NuF	31
US-Los	65	SE-Deg	23
DE-Sfn	45	US-Bog	42
CZ-Wet	54	US-Fen	56
DE-Spw	89	FI-Lom	28
IE-Kil	28	RU-Che	35
DE-Bou	34	NO-And	21
DE-Zrk	33	DK-ZaF	37
CA-Wp1	38	NO-Adv	28
SE-faj	21	PL-Kpt	52

Taylor diagrams of (a) GPP (g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$); (b) ER (g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$); (c) NEE (g C m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$); (d) LE (W m${}^{-\mathrm{2}}$); (e) $H$ (W m${}^{-\mathrm{2}}$); and (f) water table depth (cm). All statistics were calculated using daily averaged data. All points were normalized by dividing the standard deviation of model results by the standard deviation of the corresponding measurement; thus, the reference point is 1.0. Light green markers – temperate sites; dark green markers – boreal sites; blue markers – arctic sites.

[Figure omitted. See PDF]

For the 22 sites where NEE and ER measurements were available, the errors in the three carbon fluxes (GPP, ER, and NEE) were significantly reduced by optimizing $V_{\mathrm{cmax}}$ at each site (Table 4, Figs. 4, S4). With site-specific $V_{\mathrm{cmax}}$ values (site-by-site model performances are shown in Figs. S5 to S10), the overall (all the daily data from all the 22 sites) performance of the model was high for GPP ($r^{\mathrm{2}}=$ 0.76, MEF $=$ 0.76) and ER ($r^{\mathrm{2}}=$ 0.78, MEF $=$ 0.75), and lower for NEE ($r^{\mathrm{2}}=$ 0.38, MEF $=$ 0.26; Fig. 4, Table 4). Seasonal variations in carbon fluxes were well captured by the model ($r^{\mathrm{2}}=$ 0.61 to 0.86). The spatial across-site gradients of annual mean GPP and ER were generally good, with $r^{\mathrm{2}}$ of 0.93 and 0.89, and lower for NEE ($r^{\mathrm{2}}=$ 0.27). Compared to simulations with a fixed $V_{\mathrm{cmax}}$ (the mean of the optimized values of 40 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$, there were large improvements in capturing spatial gradients of carbon fluxes with a site-specific $V_{\mathrm{cmax}}$ (e.g., $r^{\mathrm{2}}$ increased from 0.20 to 0.93, from 0.27 to 0.89, and from 0.16 to 0.27 for GPP, ER, and NEE, respectively, while the RMSD was reduced by 63, 48, and 9 %). This result indicates that model–data disagreement can be largely reduced by using site-specific $V_{\mathrm{cmax}}$ instead of a fixed (mean) value. In future regional simulations, spatial variations in $V_{\mathrm{cmax}}$ should be taken into account. There was, however, no significant improvement in LE, H, and WT by using site-specific $V_{\mathrm{cmax}}$ values (Table 4). The model performance was poor for predicting daily anomalies of all fluxes, with $r^{\mathrm{2}}$ < 0.20. For both temporal and spatial variation, the MEF values of the WT were negative, and $r^{\mathrm{2}}$ smaller than 0.10, indicating that the model had a low predictive capability for the WT. Possible reasons for this could be (1) peat disturbance was not parameterized; i.e., the removal of beaver dams resulted in a decline of water level at US-Los; water levels at US-WPT, CZ-Wet, and RU-Che were manipulated. (2) The model diagnosed all peatland sites as fens by routing runoff from non-peatland areas into the peatland soil tile, whereas in reality, bogs receive water and nutrients only through precipitation. In other words, we included an extra water source for bogs other than rainfall. However, the model did not perform better for fens (Fig. 3f), possibly because the amount of water that was routed into the fen was in error. (3) WT depends on water input from surrounding non-peatland areas: the greater the peatland fraction in the grid cell, the smaller the runoff input from other soils to the peatland, hence resulting in a deeper water table in the peatland (Fig. S11). The peatland area fraction derived from the map of Yu et al. (2010) cannot represent the local area providing water for fens. (4) For global applications, the effects of micro-relief were not represented in the model, although they have been shown to be an important regulator of the local hydrology cycle (Gong et al., 2012; Shi et al., 2015).

Observed ($x$ axis) vs. simulated ($y$ axis) fluxes (GPP, ER, NEE, LE, H, and WT) at the 22 sites where GPP derived from EC measurements were available. Fluxes were simulated using site-specific optimized $V_{\mathrm{cmax}}$. The colors of points indicate the amount of data in each bin; in panel (b), each data point represents one peatland site. The red line shows the observations equal to the simulations.

[Figure omitted. See PDF]

Model performance measures for GPP, ER, NEE, LE, $H$, and WT. The left-hand column shows results with site-specific optimized $V_{\mathrm{cmax}}$ at each site; the right-hand column shows results with the fixed $V_{\mathrm{cmax}}$ (40 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$ at all sites.

	Site-specific optimized $V_{\mathrm{cmax}}$					Mean $V_{\mathrm{cmax}}$ (constant value, 40 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$
Flux	RMSD	SDSD	LCS	$r^{\mathrm{2}}$	MEF	RMSD	SDSD	LCS	$r^{\mathrm{2}}$	MEF
	Overall (daily variability)					Overall (daily variability)
GPP	1.39	0.11	1.80	0.76	0.76	2.17	0.06	4.60	0.47	0.41
ER	0.83	0.09	0.52	0.78	0.75	1.09	0.14	1.04	0.57	0.56
NEE	1.30	0.02	1.56	0.38	0.26	1.48	0.00	2.01	0.29	0.03
LE	31.67	21.65	932.76	0.42	0.14	31.67	21.19	933.95	0.42	0.14
$H$	35.40	96.59	1151.28	0.24	$-$0.50	35.40	97.21	1150.59	0.24	$-$0.50
WT	25.93	10.26	661.80	0.01	$-$0.56	26.14	7.63	675.51	0.01	$-$0.59
	Across-site variability					Across-site variability
GPP	0.41	0.03	0.10	0.93	0.89	1.11	0.42	0.80	0.20	0.19
ER	0.38	0.01	0.06	0.89	0.79	0.72	0.16	0.33	0.27	0.23
NEE	0.60	0.06	0.20	0.27	$-$0.01	0.66	0.17	0.13	0.16	$-$0.21
LE	9.85	1.13	65.49	0.71	0.50	9.80	1.04	65.21	0.71	0.50
$H$	14.31	2.67	155.85	0.01	$-$1.04	14.28	2.83	154.38	0.01	$-$1.03
WT	24.40	15.20	444.83	0.02	$-$0.82	25.10	4.65	478.84	0.03	$-$0.92
	Mean seasonal variability					Mean seasonal variability
GPP	0.92	0.03	0.81	0.86	0.86	1.36	0.02	1.83	0.70	0.69
ER	0.51	0.05	0.22	0.86	0.86	0.65	0.05	0.37	0.77	0.77
NEE	0.80	0.00	0.64	0.61	0.54	0.95	0.01	0.88	0.50	0.35
LE	11.49	7.75	124.23	0.83	0.78	11.47	7.46	124.02	0.83	0.78
$H$	17.85	65.77	252.65	0.57	0.11	17.85	66.40	252.30	0.57	0.11
WT	9.87	8.32	88.88	0.06	$-$1.38	9.77	12.73	82.69	0.12	$-$1.33
	Anomalies					Anomalies
GPP	1.03	0.03	1.02	0.18	0.01	1.10	0.02	1.19	0.13	$-$0.13
ER	0.61	0.08	0.29	0.19	0.17	0.64	0.07	0.34	0.16	0.10
NEE	0.96	0.12	0.81	0.07	$-$0.07	0.99	0.12	0.85	0.04	$-$0.14
LE	27.43	26.14	726.25	0.07	$-$0.94	27.46	26.19	727.76	0.07	$-$0.94
$H$	28.09	81.43	707.43	0.12	$-$1.12	28.10	82.12	707.49	0.12	$-$1.12
WT	13.25	0.40	174.69	0.10	$-$0.47	13.43	0.47	179.41	0.09	$-$0.51

To better understand the influence of the water table dynamics on ER and NEE in the model, we compared the second set of simulations (S2, with observed water table used in the carbon module to define the fraction of oxic and anoxic decomposition in the acrotelm) with the first set (S1, water table calculated by the model). ORCHIDEE-PEAT showed only a small improvement in reproducing ER and NEE when WT${}_{\mathrm{obs}}$ was used (Tables 5 and 6). To illustrate this effect, we took the Lompolojänkkä (FI-Lom) fen site as an example, in which WT was most severely underestimated among the 22 sites where NEE and ER measurements were available (Fig. S8). While modeled WT varied between 5 and 54 cm below the surface, WT${}_{\mathrm{obs}}$ was always above the soil surface. Figure 5a showed that in comparison to S1, there was no aerobic respiration and larger anaerobic respiration in the acrotelm in S2. Due to the smaller acrotelm respiration (aerobic plus anaerobic) in S2, carbon input from acrotelm to catotelm was larger and consequently, more carbon accumulated in the catotelm in S2. Thus, the catotelm respiration in S2 was higher than that in S1 (Fig. 5c), even though the catotelm respiration rate was very small. Because the growth of the peatland vegetation was not constrained by water in the model, the simulated GPP values were similar between S1 and S2 (Fig. 5a). With similar GPP but smaller soil respiration (sum of the acrotelm and the catotelm respiration), S2 simulations thus resulted in more negative NEE values than S1 (higher net CO${}_{\mathrm{2}}$ uptake). Simulated leaf onset occurred earlier than observed at the Lompolojänkkä site, causing the ecosystem to switch from carbon source to carbon sink in May, while the start of the carbon uptake was observed to occur later (Fig. 5b). Although the modeled NEE was similar in amplitude to the observations, the day-to-day variations of this flux were not captured (Fig. 6), causing an overestimation (more negative values) of NEE in the warm period (May–September).

Monthly mean (averaged over 2007–2009) of (a) GPP and ecosystem respiration (ER); (b) NEE; (c) catotelm respiration at the Lompolojänkkä fen site (FI-Lom). S1: simulated WT was used in the carbon module; S2: observed WT values (WT${}_{\mathrm{obs}})$ were used; ob: measured NEE. The graph inserted shows catotelm respiration. By convention, a source of CO${}_{\mathrm{2}}$ to the atmosphere is a positive number.

[Figure omitted. See PDF]

The influence of WT on respiration was parameterized as the separation of oxic ($\mathit{\beta }$ in Eq. 6) vs. anoxic (1$\mathit{\beta }$ in Eq. 7) decomposition in the acrotelm. Although absolute values of simulated WT in S1 and WT${}_{\mathrm{obs}}$ in S2 were quite different (Fig. S8), the values of $\mathit{\beta }$ were not very different (Fig. S12). Therefore, the simulated WT was good enough to properly replicate ER (Fig. S13). An additional simulation (S3) performed at FI-Lom showed that if WT was more severely underestimated, i.e., WT in S3 was consistently 20 cm deeper than in S1, the acrotelm was exposed to oxygen for a longer time, resulting in larger ER and hence smaller carbon sequestration in S3 (Figs. S12, S13).

Model performance measures of ER simulations for the site-by-site comparison, comparison across sites, mean seasonal cycle, and anomalies, using modeled (S1) and observed (S2) WT.

	Modeled WT used (S1)					Observed WT used (S2)
Site	RMSD	SDSD	LCS	$r^{\mathrm{2}}$	MEF	RMSD	SDSD	LCS	$r^{\mathrm{2}}$	MEF
CZ-Wet	1.45	0.86	0.87	0.81	0.68	1.51	1.05	0.79	0.81	0.66
DE-Bou	0.78	0.03	0.50	0.69	0.64	0.77	0.03	0.50	0.69	0.65
DE-Sfn	0.96	0.10	0.79	0.61	0.59	0.97	0.09	0.82	0.60	0.58
FI-Lom	0.46	0.00	0.19	0.85	0.84	0.45	0.02	0.18	0.85	0.84
IE-Kil	0.44	0.01	0.01	0.09	0.51	0.42	0.01	0.01	0.13	0.48
SE-Deg	0.69	0.26	0.19	0.75	0.62	0.64	0.16	0.23	0.75	0.68
SE-Faj	0.58	0.07	0.08	0.87	0.60	0.59	0.08	0.07	0.88	0.59
US-Los	0.63	0.01	0.39	0.85	0.85	0.60	0.00	0.35	0.87	0.87
Overall	0.79	0.09	0.51	0.78	0.76	0.79	0.09	0.51	0.78	0.76
Across sites	0.31	0.01	0.06	0.82	0.76	0.32	0.01	0.06	0.82	0.74
Seasonal	0.45	0.06	0.15	0.91	0.89	0.44	0.07	0.13	0.92	0.89
Anomalies	0.62	0.07	0.31	0.21	0.19	0.63	0.08	0.31	0.20	0.17

Model performance measures of NEE simulations for the site-by-site comparison, comparison across sites, mean seasonal cycle, and anomalies, using modeled (S1) and observed (S2) WT.

	Modeled WT used (S1)					Observed WT used (S2)
Site	RMSD	SDSD	LCS	$r^{\mathrm{2}}$	MEF	RMSD	SDSD	LCS	$r^{\mathrm{2}}$	MEF
CZ-Wet	2.97	3.61	4.38	0.46	0.37	2.86	3.22	4.27	0.50	0.41
DE-Bou	1.30	0.02	1.40	0.31	$-$0.21	1.31	0.03	1.41	0.31	$-$0.23
DE-Sfn	2.98	2.98	4.27	0.20	0.02	2.98	3.08	4.15	0.21	0.02
FI-Lom	1.05	0.01	0.94	0.46	0.21	1.08	0.02	0.95	0.49	0.16
IE-Kil	0.48	0.000	0.16	0.29	$-$0.37	0.49	0.002	0.16	0.32	$-$0.44
SE-Deg	0.64	0.03	0.33	0.51	0.09	0.57	0.01	0.29	0.51	0.26
SE-Faj	0.65	0.01	0.33	0.31	$-$0.36	0.65	0.02	0.33	0.32	$-$0.39
US-Los	3.15	0.05	8.78	0.47	$-$3.37	3.10	0.06	8.57	0.39	$-$3.23
Overall	1.95	0.20	3.52	0.02	$-$0.35	1.92	0.18	3.42	0.04	$-$0.31
Across sites	0.67	0.27	0.16	0.40	0.29	0.65	0.26	0.14	0.46	0.32
Seasonal	1.30	0.05	1.64	0.25	0.13	1.27	0.03	1.58	0.28	0.17
Anomalies	1.18	0.22	1.17	0.003	$-$0.34	1.17	0.21	1.17	0.001	$-$0.33

Relationship between optimized $V_{\mathrm{cmax}}$ and meteorological variables

Several univariate ANOVA models were used to explain the spatial gradient of optimized $V_{\mathrm{cmax}}$, explanatory variables including air temperature ($T$), precipitation ($P$), net radiation (NET_RAD), water use efficiency (WUE), water balance (WB), and latitude (LAT). All explanatory variables were calculated as daily mean values during the growing season. Water use efficiency (g C m${}^{-\mathrm{2}}$ mm${}^{-\mathrm{1}}$ H${}_{\mathrm{2}}$O) was calculated as the ratio of GPP and evapotranspiration (ET). Water balance (mm day${}^{-\mathrm{1}})$ was calculated as the difference between precipitation and ET.

There was no significant difference between optimized $V_{\mathrm{cmax}}$ among peatland types (fen vs. bog, $p$ $=$ 0.16), climate zones (temperate vs. boreal vs. arctic, $p$ $=$ 0.17), or dominant vegetation types (grasses and/or mosses dominated vs. shrubs and/or trees dominated, $p$ $=$ 0.67; Fig. S14). However, we found a significant positive relationship between $V_{\mathrm{cmax}}$ and the growing season mean air temperature (Fig. S15, Table 6, $V_{\mathrm{cmax}}=\mathrm{2.78}T+\mathrm{8.74}$, with $r^{\mathrm{2}}=$ 0.19, $p$ < 0.05) and a significant negative relationship between $V_{\mathrm{cmax}}$ and the latitude (Fig. S15, Table 6, $V_{\mathrm{cmax}}=$ $-$0.92LAT $+$ 93.56, with $r^{\mathrm{2}}=$ 0.23, $p$ < 0.05).

To verify the applicability of the empirical relationship found across sites between optimized $V_{\mathrm{cmax}}$ and the latitude (Fig. S15), we used the seven sites where there were no GPP observations available (US-Bes, DE-Hmm, US-Ics, PL-wet, SE-Sto, CA-Wp2, and CA-Wp3) as cross-validated sites. We compared model performance in simulating NEE, with $V_{\mathrm{cmax}}$ being calculated according to the empirical relationship, and with $V_{\mathrm{cmax}}$ being fixed to its mean value of all 22 sites from Table 3 (40 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$. The model performance in reproducing spatial gradients of NEE was improved when the $V_{\mathrm{cmax}}$ values derived from the empirical relationship were used (Fig. S16b, with RMSD reduced by 11 %, $r^{\mathrm{2}}$ increased from 0.20 to 0.38, and MEF increased from $-$0.04 to 0.17). This implies that, compared to a fixed $V_{\mathrm{cmax}}$, the usage of $V_{\mathrm{cmax}}$ value from the empirical relationship can better capture spatial gradients of NEE. It is worth mentioning that the empirical relationship was built on climate conditions from the last two decades (1999–2015) and thus may change in the future when the climate changes.

The results of the ANOVA analysis – the variance of optimized $V_{\mathrm{cmax}}$ in relation to chosen variables.

Variable	$F$-ratio	$p$-value	$r^{\mathrm{2}}$ (%)
$T$	4.67	0.04${}^{*}$	18.95
$P$	0.95	0.34	4.52
NET_RAD	0.22	0.64	1.11
WUE	0.39	0.54	1.91
WB	1.35	0.26	6.32
LAT	6.08	0.023${}^{*}$	23.30

${}^{*}$ Indicates statistical significance at a significance level of 0.05.

Observed and simulated daily mean NEE at the FI-Lom fen site in (a) S1 (simulated WT was used in the carbon module) and (b) S2 (modeled water table was assimilated to WT${}_{\mathrm{obs}}$ and was used in the carbon module).

[Figure omitted. See PDF]

Soil temperature and a snow depth underestimation in the model

For most of the sites, soil temperature was underestimated in winter and overestimated in summer by our model (Figs. 7 and 8; results from sites DK-NuF and CA-Wp1 are shown as illustrative examples). One possible reason for the underestimation of soil temperature in winter is the underestimation of snow depth (Fig. 9), since snow insulates the soil-changing thermal conditions in comparison to a snow-free surface. The underestimation of the snow depth can be caused by the bias in snow processes of the model, such as underestimation of snow mass, and/or overestimation of snow density and subsequently overestimation of snow compaction, and/or overestimation of sublimation. The insulation effect of the moss layer and the top organic layer is not included in this study, which may explain why soil temperature was overestimated in summer but underestimated in winter. ORCHIDEE-PEAT calculates one energy budget for the vegetation and soil columns in one grid cell. Key parameters used for solving the heat diffusion equations in the soil, such as soil heat capacity and thermal conductivity, were prescribed by the dominant soil texture in the grid cell (Gouttevin et al., 2012). Nevertheless, similarly to the case of the hydrology module, the three default (coarse, medium, fine) soil textures cannot represent thermal properties of a peat soil (Paavilainen and Päivänen, 1995; Abu-Hamdeh and Reeder, 2000).

Measured (a) and simulated (b) soil temperature, and their differences (c) at the DK-NuF (64.13${}^{\circ }$ N, 51.39${}^{\circ }$ W) fen site. Soil temperature was measured at 2, 10, 20, 50, and 70 cm below soil surface. To compare simulated soil temperatures with the measurements, we linearly interpolated simulated soil temperature in different layers to the depths of the measurements.

[Figure omitted. See PDF]

Discussion

ORCHIDEE-PEAT groups various peatland vegetation into one plant functional type (PFT). This PFT cannot represent the true range in vegetation composition (shrubs, sedges, mosses, etc.) of peatlands. However, by optimizing the value of $V_{\mathrm{cmax}}$ at each site, simulated GPP well represented observations and yielded reasonable soil carbon input. The $V_{\mathrm{cmax}}$ values estimated in this study ranged from 19 to 89 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, with a mean value of 40 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. These values were not fully comparable with values reported for a specific vegetation type, as they are averages for all plants growing in the peatland ecosystem. As stated in Sect. 2.2, observed $V_{\mathrm{cmax}}$ varies strongly among different species and sites. $V_{\mathrm{cmax}}$ of mosses at the Old Black Spruce site (Canada) ranged from 5 to 14 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ (Williams and Flanagan, 1998). In a nutrient addition experiment conducted by Bubier et al. (2011), $V_{\mathrm{cmax}}$ for ericaceous shrubs in a temperate bog ranged from 67 to 137 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, with $V_{\mathrm{cmax}}$ for Vaccinium myrtilloides, Ledum groenlandicum, and Chamaedaphne calyculata valued at 84.6 $\pm$ 13.5, 78.1 $\pm$ 13.4, and 132.1 $\pm$ 31.2 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ in the plots with no nutrient addition. The optimized model $V_{\mathrm{cmax}}$ in our study was within the range of these observations. Meanwhile, the values we inferred from sites to match peak GPP are comparable to those used in other land surface models: the McGill wetland model used a value of 17 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for evergreen shrubs (St-Hilaire et al., 2010); the CLASS-CTEM model (Wu et al., 2016) used 60, 50, and 40 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for evergreen shrubs, deciduous shrubs, and sedges, respectively; the values for mosses in these two models were adapted from the study of Williams and Flanagan (1998). Here, we found that optimized $V_{\mathrm{cmax}}$ has a significant positive relationship with temperature and a significant negative relationship with latitude of chosen peatland sites. A decrease of $V_{\mathrm{cmax}}$ with latitude in the Northern Hemisphere, like the one inferred from optimized site values, has also been documented by Walker et al. (2017), who assumed that $V_{\mathrm{cmax}}$ was constrained by the rate of N uptake, with the rate of N uptake calculated as a function of soil C, N, and mean annual air temperature. We speculate that the dependence of optimized $V_{\mathrm{cmax}}$ on latitude found in Sect. 4.2 can be attributed to two effects. First, there is an increase of the length of the growing season as latitude decreases. Simultaneously, temperature and incoming solar radiation increase. The longer growing season may enhance vegetation productivity (Fang et al., 2003; Nemani et al., 2003; Piao et al., 2007). Second, temperature influences the nutrient availability for plants. The decomposition of plant litter and the release of nitrogen can be enhanced by high temperature, although litter decomposition is also driven by soil moisture, vegetation composition, litter quality, and their interactions with temperature (Aerts, 2006; Cornelissen et al., 2007; Gogo et al., 2016). Because nitrogen (N) is one key element in proteins that are involved in the photosynthesis process, photosynthesis capacity is highly correlated to N availability (Evans, 1989; Takashima et al., 2004; Walker et al., 2014). Since the N cycle is not explicitly included in ORCHIDEE-PEAT, the relationship between $V_{\mathrm{cmax}}$ and the latitude (and temperature) possibly reflected the impact of N on photosynthesis rates.

Previous studies have shown that peatlands can have contrasting responses to variations in water table depth. Concerning sites analyzed in our study, Aurela et al. (2007) reported that at the nutrient-poor fen FI-Sii site, drought increased respiration and thus diminished carbon uptake; Adkinson et al. (2011) reported that reduced water availability constrained photosynthesis capacity at the rich fen CA-Wp3 and consequently suppressed NEE, while the poor fen CA-Wp2 did not show a significant response to the lower water table. At the moderately rich treed fen CA-Wp1 site, Flanagan and Syed (2011) reported that both photosynthesis and respiration increased in response to the warmer and drier conditions; Hurkuck et al. (2016) stated that temperature and light played a more important role than water table depth in controlling respiration and photosynthesis at the DE-Bou bog. Based on the field observations, the timing, duration, and intensity of drought have a major impact on the responses of peatland ecosystems. Lund et al. (2012) demonstrated that at the raised bog SE-Faj, a relatively short but severe drought that occurred in the middle of growing season of 2006 amplified respiration while a long-lasting drought that occurred at the beginning of growing season of 2008 reduced GPP. Lafleur et al. (2005) and Sulman et al. (2009) concluded from their studies at the CA-Mer bog and US-Los fen that wetter peatlands would show a stronger relationship between respiration and water table than drier peatlands because in a narrow range of the upper soils, small increases in WT (shallower WT) can result in a large increase in soil water content and therefore respiration decrease, while below a critical level, soil water content shows only small increase with increasing WT, and respiration changes are not so pronounced. Sulman et al. (2010) found that wetter conditions decreased respiration at fens but increased respiration at bogs, mainly due to different vegetation composition at these two types of peatlands: the fen sites had more shrubs and sedges while the bog sites had more mosses. In this study, we did not distinguish between fens and bogs, and growth of peatland vegetation was not constrained by water table depth in the model. Therefore, the sensitivity of GPP to WT fluctuations in observations was not included in the model. As a consequence, the model captured neither the reported decrease of photosynthesis due to drought at CA-Wp3 (Adkinson et al., 2011) and SE-Faj (Lund et al., 2012) nor the increase of photosynthesis as a result of lower water table at CA-Wp1 (Flanagan and Syed, 2011). However, the model can reproduce the pattern where, above a critical level (acrotelm depth), peat respiration decreases with increasing WT (Figs. 5, S13), as reported at CA-Mer and US-Los (Lafleur et al., 2005; Sulman et al., 2009). ORCHIDEE-PEAT adequately captured the daily, seasonal, and across-site annual variations in GPP (with $r^{\mathrm{2}}=$ 0.75, 0.86, and 0.93, respectively) and ER (with $r^{\mathrm{2}}=$ 0.78, 0.86, and 0.89, respectively) but did not perform as well in reproducing NEE variations (with $r^{\mathrm{2}}=$ 0.38, 0.61, and 0.27, respectively). Note that in the two-layer soil carbon scheme, the dependence of soil respiration on temperature was parameterized as an exponential function of the soil layer-weighted average temperature (Text S1); the vertical temperature gradient in the soil profile was ignored by the model. However, field studies have shown that soil temperature is one of the most important predictors of respiration, and values of $Q_{\mathrm{10}}$ coefficient depend on the soil depth (Lafleur et al., 2005; D'Angelo et al., 2016).

Measured (a) and simulated (b) soil temperature, and their difference (c) at the CA-Wp1 (54.95${}^{\circ }$ N, 112.47${}^{\circ }$ W) fen site. The measured soil temperature (a) is the mean of a hummock and a hollow. Soil temperature was measured at 2, 10, 20, 50, and 100 cm below soil surface. To compare simulated soil temperatures with the measurements, we linearly interpolated simulated soil temperature in different layers to the depths of the measurements.

[Figure omitted. See PDF]

Simulated vs. measured snow depth (m) at the (a) DK-NuF and (b) CA-Wp1 fen sites.

[Figure omitted. See PDF]

Correct representation of peatland hydrology is a challenging problem in large-scale land surface models (Wania et al., 2009a; Wu et al., 2016). The simulated water table by ORCHIDEE-PEAT depends on water inflows from the surrounding non-peatland areas, and a water-routing analysis on subgrid scales can be included to improve the model performance for water table in the future (Ringeval et al., 2012; Stocker et al., 2014). Other studies have shown that microtopography exerts important influences on hydrological dynamics of peatlands; however, to capture the influence of microtopography on water table, high-resolution microtopographic feature and vegetation information are needed (Gong et al., 2013; Shi et al., 2015).

The poor correspondence between simulated and observed energy fluxes was not completely unexpected, since ORCHIDEE-PEAT only calculates one energy budget for the whole grid cell and not for each soil tile/PFT present in the same grid cell. A site-varied and/or time-varied correction of LE and H measurements to force energy balance closure, and parameterizations of an independent energy budget in peatlands would be helpful for better comparison of simulated and observed energy fluxes in peatlands.

Conclusions

We developed ORCHIDEE-PEAT to simulate soil hydrology and carbon dynamics in peatlands. The model was evaluated at 30 northern peatland sites (Europe, USA, Canada, and Russia). The optimization of $V_{\mathrm{cmax}}$ reduced the errors in the simulated carbon budget. The model, generally, reproduced the spatial gradient and temporal variations in GPP, ER, and NEE well. Water table depth was poorly simulated, possibly due to uncertainties in water input from non-peatland areas in the grid cell, and to a lack of representation of micro-relief, as well as the lack of consideration of peat disturbance. A significant relationship between $V_{\mathrm{cmax}}$ and latitude was found. This may be attributed to the influence of temperature on growing season length and nutrient availability. For ER and NEE, the improvement brought by forcing the carbon module to use observed WT values (WT${}_{\mathrm{obs}})$, instead of being calculated by the model, was small, indicating that the simulated WT was reliable to predict ER and NEE properly.

Our study shows that in order to reproduce spatial gradients of NEE for northern peatlands, an average $V_{\mathrm{cmax}}$ value is not sufficient. To represent a spatial gradient of carbon fluxes in large-scale simulations of northern peatlands, incorporating the peatland nitrogen cycle would be helpful. Alternatively, an empirical relationship between $V_{\mathrm{cmax}}$ and the latitude (temperature) may be used as a proxy of nitrogen availability. Effects of water table variations on soil carbon decomposition are modeled as the partitioning of the acrotelm layer into oxic and anoxic zones, but effects of water table changes on GPP were not modeled in this study. Future priorities for improving ORCHIDEE-PEAT include better representing the influence of the water table on photosynthesis and depth-dependent influence of soil temperature on soil respiration, as well as including an independent subgrid energy budget for peatland areas.

The access to the source code is available online via (http://forge.ipsl.jussieu.fr/orchidee/browser/perso/chunjing.qiu/ORCHIDEE, but its access is restricted. Readers interested in running the model should follow the instructions at http://orchidee.ipsl.fr/index.php/you-orchidee and contact the corresponding author for a username and password.

Measured eddy-covariance fluxes and related meteorological data can be obtained from the FLUXNET database (http://fluxnet.ornl.gov/), the AmeriFlux database (http://ameriflux.lbl.gov/), and from investigators upon request. Model outputs are available at https://files.lsce.ipsl.fr/public.php?service=files&t=0f319ede335dc37d43edf617c94f83d0.

The Supplement related to this article is available online at https://doi.org/10.5194/gmd-11-497-2018-supplement.

The authors declare that they have no conflict of interest.

Acknowledgements

This study was supported by the European Research Council Synergy grant ERC-2013-SyG-610028 IMBALANCE-P. We would like to thank all the PIs for giving us permission to use the flux and ancillary data and for all the help and advice they provided while we were preparing the manuscript. We thank the Polish National Science Centre which provided funds for site Kopytkowo (PL-Kpt) under projects UMO-2011/01/B/ST10/07550 and UMO-2015/17/B/ST10/02187, and the Department of Energy for supporting measurements at Lost Creek fen (US-Los) through the AmeriFlux Network Management Project. We gratefully acknowledge the financial support provided for the La Guette site under the Labex VOLTAIRE (ANR-10-LABX-100-01) and the PIVOTS project of the Région Centre – Val de Loire (ARD 2020 program and CPER 2015–2020). Data from the Greenlandic sites (DK-ZaF and DK-NuF) were provided by the Greenland Ecosystem Monitoring Programme. The US-Bes tower is funded by NSF (award nos. 1204263 and 1702797), NASA ABoVE (NNX15AT74A; NNX16AF94A), EU Horizon 2020 INTAROS (under grant agreement no. 727890), and the NERC UAMS grant (NE/P002552/1). Edited by: Carlos Sierra Reviewed by: two anonymous referees