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

. Peatlands store substantial amounts of carbon and are vulnerable to climate change. We present a modiﬁed version of the Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) land surface model for simulating the hydrology, surface energy, and CO 2 ﬂuxes of peatlands on daily to annual timescales. The model includes a separate soil tile in each 0.5 ◦ grid cell, deﬁned from a global peatland map and identiﬁed with peat-speciﬁc 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 cmax ) being optimized at each site. Regarding short-term day-to-day variations, the model performance was good for gross primary production (GPP) ( r 2 = 0.76; Nash– Sutcliffe modeling efﬁciency, MEF = 0.76) and ecosystem respiration (ER, MEF

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 2 fluxes of peatlands on daily to annual timescales. The model includes a separate soil tile in each 0.5 • 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 cmax ) being optimized at each site. Regarding short-term day-to-day variations, the model performance was good for gross primary production (GPP) (r 2 = 0.76; Nash-Sutcliffe modeling efficiency, MEF = 0.76) and ecosystem respiration (ER, r 2 = 0.78, MEF = 0.75), with lesser accuracy for latent heat fluxes (LE, r 2 = 0.42, MEF = 0.14) and and net ecosystem CO 2 exchange (NEE, r 2 = 0.38, MEF = 0.26). Seasonal variations in GPP, ER, NEE, and energy fluxes on monthly scales showed moderate to high r 2 values (0.57-0.86). For spatial across-site gradients of annual mean GPP, ER, NEE, and LE, r 2 values of 0.93, 0.89, 0.27, and 0.71 were achieved, respectively. Water table (WT) variation was not well predicted (r 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 cmax and latitude (temperature), which better reflects the spatial gradients of annual NEE than using an average V 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 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 feed-Geosci. Model Dev., 11, 497-519, 2018 www.geosci-model-dev.net/11/497/2018/ backs 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 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 largescale applications in an Earth system modeling context. 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 3 and C 4 grasses, two for C 3 and C 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 s-ref ) and the saturated and residual volumetric water contents (θ s , θ r ). Dominant soil textural classes are taken from the Zobler's soil texture map (Zobler, 1986) at 1 • 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).
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 −1 ) and diffusivity (D, m 2 s −1 ) as a function of volumetric water content (θ, m 3 m −3 ): Figure 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 −1 ) and diffusivity (D, mm 2 day −1 ). Note that the horizontal axis is on a base-10 logarithmic scale. Table 1. Van Genuchten parameters used for different soil texture classes for non-peat soils (coarse, medium, fine) and for peat. θ s is the saturated water content (m 3 m −3 ); θ r is the residual water content (m 3 m −3 ); K s-ref is the reference saturated hydraulic conductivity (m s −1 ); α is the inverse of the air entry suction (m −1 ); n is a dimensionless parameter. In Eqs. (1) and (2), m = 1 − 1/n.
, while a root-fracturing factor increases K s where roots are denser (F root (z)): with F d (z) = min (max (exp (−f (z − z lim )) , 0.1) , 1) and reference top-soil saturated hydraulic conductivity determined by soil texture (m s −1 ), K max s is the value of the coarser (sandy) texture and equals 8.25 × 10 −5 m s −1 , α 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 floodtolerant C 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 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 3 grass with a shallower rooting depth to represent the average of vegetation growing in northern peatlands.
Only one key photosynthetic parameter (V 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 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 cmax for Sphagnum at the Old Black Spruce site (53.985 • N, 105.12 • W) in Canada was 5, 14, and 6 µmol m −2 s −1 during spring, summer, and autumn, respectively, while V cmax for Pleurozium was 7, 5, and 7 µmol m −2 s −1 during the three seasons (Williams and Flanagan, 1998). Bui (2013) conducted a fertilization experiment at the Mer Bleue Bog (Canada; 45.41 • N, 75.52 • W) on the dominant ericaceous shrub and reported that V cmax values ranged between 6 and 179 µmol m −2 s −1 , with significantly higher V cmax values after addition of nitrogen (6.4 g N m −2 yr −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 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 cmax value with environmental and climate variables. We note that this adjustment of V 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 abovesurface 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 reinfiltrates 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 peatsaturated 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 × 10 −2 to 7.16 × 10 −6 m s −1 in bogs. Gnatowski et al. (2010) measured values of 5 × 10 −6 m s −1 in a mosscovered peat, which was 2 orders of magnitude larger than for a woody peat (5.56 × 10 −8 m s −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 × 10 −5 m s −1 is comparable to the harmonic mean value (6 × 10 −5 m s −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): where θ f i is the relative volumetric water content of the ith soil layer, θ s is the saturated water content (m 3 m −3 ), θ r is the residual water content (m 3 m −3 ), dz i is the distance between node i −1 and node i ( Fig. 1; m), H tot is the total soil column height being fixed to 2.0 m, and H 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 (β) vs. anoxic conditions (1β). Acrotelm height (H A , Eq. 10) was calculated from acrotelm carbon stock (C A in Eqs. 5-7), acrotelm carbon fraction (C f,A ) and acrotelm bulk density (ρ A ). Decomposition of peat carbon is controlled by temperature (f T ) and parameterized as an exponential function: f T = Q 10 exp((T − T ref )/10 • C) with Q 10 = 2.0 and T ref = 30 • C (Text S1). Soil carbon fluxes are given by where F AC is the carbon flux from acrotelm to catotelm; R A,o is aerobically decomposed acrotelm carbon; R A,a is anaerobically decomposed acrotelm carbon; R C is decomposed carbon in catotelm; C A is carbon stored in the acrotelm; C C is carbon stored in the catotelm; and β 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 p = 1.91 × 10 −2 yr −1 , the acrotelm decomposition rate k A = 0.067 yr −1 , the catotelm decomposition rate k C = 3.35 × 10 −5 yr −1 , the ratio of anaerobic to aerobic CO 2 production µ = 0.35, carbon fraction in the acrotelm peat C f,A = 0.50, the acrotelm density ρ A = 35.0 kg m −3 , carbon fraction in the catotelm peat C f,C = 0.52, and the catotelm density ρ C = 91.0 kg m −3 .

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

Sites description
To evaluate the performance of ORCHIDEE-PEAT in simulating CO 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 ∼ 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 • 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 • 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 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 • 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 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 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 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 × 10 −5 yr −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 ∼ 2200 years . 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 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 obs ). That is, soil moisture at layers below the measured water table was prescribed as saturated (θ (z > WT obs ) = θ s ), while soil moisture above WT obs was simulated. WT 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 obs time series (Table 2). For these sites, the gaps of WT 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 2 , and CH4 fluxes (Dušek et al., 2009;Parmentier et al., 2011;Strack et al., 2006). Fixing the simulated water table to WT 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.
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 s ) and the measurement (SD m ) shows if the model can reproduce the magnitude of fluctuation among the n measurements.
wherex is simulated mean value;ȳ 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.
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 4 Results

Site-specific V 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 cmax at each site were listed in Table 3. The optimized V cmax varied from 19 to 89 µmol m −2 s −1 (Table 3), with a mean value of 40 µmol m −2 s −1 . The calibration of V 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).
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 cmax at each site ( ; 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. variations in carbon fluxes were well captured by the model (r 2 = 0.61 to 0.86). The spatial across-site gradients of annual mean GPP and ER were generally good, with r 2 of 0.93 and 0.89, and lower for NEE (r 2 = 0.27). Compared to simulations with a fixed V cmax (the mean of the optimized values of 40 µmol m −2 s −1 ), there were large improvements in capturing spatial gradients of carbon fluxes with a site-specific V cmax (e.g., r 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 cmax instead of a fixed (mean) value. In future regional simulations, spatial variations in V cmax should be taken into account. There was, however, no significant improvement in LE, H, and WT by using site-specific V cmax values (Table 4). The model performance was poor for predicting daily anomalies of all fluxes, with r 2 < 0.20. For both temporal and spatial variation, the MEF values of the WT were negative, and r 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 Geosci. Model Dev., 11, 497-519, 2018 www.geosci-model-dev.net/11/497/2018/ 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).
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 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 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 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).
The influence of WT on respiration was parameterized as the separation of oxic (β in Eq. 6) vs. anoxic (1β in Eq. 7) decomposition in the acrotelm. Although absolute values of simulated WT in S1 and WT obs in S2 were quite different (Fig. S8), the values of β 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).

Relationship between optimized V cmax and meteorological variables
Several univariate ANOVA models were used to explain the spatial gradient of optimized V 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 −2 mm −1 H 2 O) was calculated as the ratio of GPP and evapotranspiration (ET). Water balance (mm day −1 ) was calculated as the difference between precipitation and ET. There was no significant difference between optimized V 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 cmax and the growing season mean air temperature (Fig. S15, Table 6, V cmax = 2.78T + 8.74, with r 2 = 0.19, p < 0.05) and a significant negative relationship between V cmax and the latitude (Fig. S15, Table 6, V cmax = −0.92LAT + 93.56, with r 2 = 0.23, p < 0.05).
To verify the applicability of the empirical relationship found across sites between optimized V 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 cmax being calculated according to the empirical relationship, and with V cmax being fixed to its mean value of all 22 sites from Table 3 (40 µmol m −2 s −1 ). The model performance in reproducing spatial gradients of NEE was improved when the V cmax values derived from the empirical relationship were used (Fig. S16b, with RMSD reduced by 11 %, r 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 cmax , the usage of V 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)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and thus may change in the future when the climate changes.  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 . 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).

ORCHIDEE-PEAT groups various peatland vegetation into one plant functional type (PFT). This PFT cannot represent
Geosci. Model Dev., 11, 497-519, 2018 www.geosci-model-dev.net/11/497/2018/  (Williams and Flanagan, 1998). In a nutrient addition experiment conducted by Bubier et al. (2011), V cmax for ericaceous shrubs in a temperate bog ranged from 67 to 137 µmol m −2 s −1 , with V cmax for Vaccinium myrtilloides, Ledum groenlandicum, and Chamaedaphne calyculata valued at 84.6 ± 13.5, 78.1 ± 13.4, and 132.1 ± 31.2 µmol m −2 s −1 in the plots with no nutrient addition. The optimized model V 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 µmol m −2 s −1 for evergreen shrubs (St-Hilaire et al., 2010); the CLASS-CTEM model (Wu et al., 2016) used 60, 50, and 40 µmol m −2 s −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 cmax has a significant positive relationship with temperature and a significant negative relationship with latitude of chosen peatland sites. A decrease of V 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 cmax (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 obs ) were used; ob: measured NEE. The graph inserted shows catotelm respiration. By convention, a source of CO 2 to the atmosphere is a positive number. 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 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 . 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 obs and was used in the carbon module).  (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 re-lationship between V 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 se- vere 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 . However, the model can reproduce the pattern where, above a critical level (acrotelm depth), peat 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. 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 2 = 0.75, 0.86, and 0.93, respectively) and ER (with r 2 = 0.78, 0.86, and 0.89, respectively) but did not perform as well in reproducing NEE variations (with r 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 10 coefficient depend on the soil depth (Lafleur et al., 2005;D'Angelo et al., 2016).
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 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 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 obs ), instead of being calcu-lated 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 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 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.
Code availability. 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. Data availability. 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.