Coupled eco-hydrology and biogeochemistry algorithms enable the simulation of water table depth effects on boreal peatland net CO 2 exchange

Water table depth (WTD) effects on net ecosystem CO2 exchange of boreal peatlands are largely mediated by hydrological effects on peat biogeochemistry and the ecophysiology of peatland vegetation. The lack of representation of these effects in carbon models currently limits our predictive capacity for changes in boreal peatland carbon deposits under potential future drier and warmer climates. We examined whether a process-level coupling of a prognostic WTD with (1) oxygen transport, which controls energy yields from microbial and root oxidation–reduction reactions, and (2) vascular and nonvascular plant water relations could explain mechanisms that control variations in net CO2 exchange of a boreal fen under contrasting WTD conditions, i.e., shallow vs. deep WTD. Such coupling of eco-hydrology and biogeochemistry algorithms in a process-based ecosystem model, ecosys, was tested against net ecosystem CO2 exchange measurements in a western Canadian boreal fen peatland over a period of drier-weather-driven gradual WTD drawdown. A May–October WTD drawdown of ∼ 0.25 m from 2004 to 2009 hastened oxygen transport to microbial and root surfaces, enabling greater microbial and root energy yields and peat and litter decomposition, which raised modeled ecosystem respiration (Re) by 0.26 μmol CO2 m−2 s−1 per 0.1 m of WTD drawdown. It also augmented nutrient mineralization, and hence root nutrient availability and uptake, which resulted in improved leaf nutrient (nitrogen) status that facilitated carboxylation and raised modeled vascular gross primary productivity (GPP) and plant growth. The increase in modeled vascular GPP exceeded declines in modeled nonvascular (moss) GPP due to greater shading from increased vascular plant growth and moss drying from near-surface peat desiccation, thereby causing a net increase in modeled growing season GPP by 0.39 μmol CO2 m−2 s−1 per 0.1 m of WTD drawdown. Similar increases in GPP and Re caused no significant WTD effects on modeled seasonal and interannual variations in net ecosystem productivity (NEP). These modeled trends were corroborated well by eddy covariance measured hourly net CO2 fluxes (modeled vs. measured: R2∼ 0.8, slopes ∼ 1± 0.1, intercepts ∼ 0.05 μmol m−2 s−1), hourly measured automated chamber net CO2 fluxes (modeled vs. measured: R2∼ 0.7, slopes ∼ 1± 0.1, intercepts ∼ 0.4 μmol m−2 s−1), and other biometric and laboratory measurements. Modeled drainage as an analog for WTD drawdown induced by climate-changedriven drying showed that this boreal peatland would switch from a large carbon sink (NEP∼ 160 g C m−2 yr−1) to carbon neutrality (NEP∼ 10 g C m−2 yr−1) should the water table deepen by a further ∼ 0.5 m. This decline in projected NEP indicated that a further WTD drawdown at this fen would eventually lead to a decline in GPP due to water limitation. Therefore, representing the effects of interactions among hydrology, biogeochemistry and plant physiological ecology on ecosystem carbon, water, and nutrient cycling in global carbon models would improve our predictive capacity for changes in boreal peatland carbon sequestration under changing climates. Published by Copernicus Publications on behalf of the European Geosciences Union. 5508 M. Mezbahuddin et al.: Coupled eco-hydrology and biogeochemistry algorithms

Abstract.Water table depth (WTD) effects on net ecosystem CO 2 exchange of boreal peatlands are largely mediated by hydrological effects on peat biogeochemistry and the ecophysiology of peatland vegetation.The lack of representation of these effects in carbon models currently limits our predictive capacity for changes in boreal peatland carbon deposits under potential future drier and warmer climates.We examined whether a process-level coupling of a prognostic WTD with (1) oxygen transport, which controls energy yields from microbial and root oxidation-reduction reactions, and (2) vascular and nonvascular plant water relations could explain mechanisms that control variations in net CO 2 exchange of a boreal fen under contrasting WTD conditions, i.e., shallow vs. deep WTD.Such coupling of eco-hydrology and biogeochemistry algorithms in a process-based ecosystem model, ecosys, was tested against net ecosystem CO 2 exchange measurements in a western Canadian boreal fen peatland over a period of drier-weather-driven gradual WTD drawdown.A May-October WTD drawdown of ∼ 0.25 m from 2004 to 2009 hastened oxygen transport to microbial and root surfaces, enabling greater microbial and root energy yields and peat and litter decomposition, which raised modeled ecosystem respiration (R e ) by 0.26 µmol CO 2 m −2 s −1 per 0.1 m of WTD drawdown.It also augmented nutrient mineralization, and hence root nutrient availability and uptake, which resulted in improved leaf nutrient (nitrogen) status that facilitated carboxylation and raised modeled vascular gross primary productivity (GPP) and plant growth.The increase in modeled vascular GPP exceeded declines in modeled nonvascular (moss) GPP due to greater shading from increased vascular plant growth and moss drying from near-surface peat desiccation, thereby causing a net increase in modeled growing season GPP by 0.39 µmol CO 2 m −2 s −1 per 0.1 m of WTD drawdown.Similar increases in GPP and R e caused no significant WTD effects on modeled seasonal and interannual variations in net ecosystem productivity (NEP).These modeled trends were corroborated well by eddy covariance measured hourly net CO 2 fluxes (modeled vs. measured: R 2 ∼ 0.8, slopes ∼ 1 ± 0.1, intercepts ∼ 0.05 µmol m −2 s −1 ), hourly measured automated chamber net CO 2 fluxes (modeled vs. measured: R 2 ∼ 0.7, slopes ∼ 1 ± 0.1, intercepts ∼ 0.4 µmol m −2 s −1 ), and other biometric and laboratory measurements.Modeled drainage as an analog for WTD drawdown induced by climate-changedriven drying showed that this boreal peatland would switch from a large carbon sink (NEP ∼ 160 g C m −2 yr −1 ) to carbon neutrality (NEP ∼ 10 g C m −2 yr −1 ) should the water table deepen by a further ∼ 0.5 m.This decline in projected NEP indicated that a further WTD drawdown at this fen would eventually lead to a decline in GPP due to water limitation.Therefore, representing the effects of interactions among hydrology, biogeochemistry and plant physiological ecology on ecosystem carbon, water, and nutrient cycling in global carbon models would improve our predictive capacity for changes in boreal peatland carbon sequestration under changing climates.

Introduction
Northern boreal peatlands have been accumulating carbon (C) at a rate of about 20-30 g m −2 yr −1 over several thousand years (Gorham, 1991;Turunen et al., 2002).Drier and warmer future climates can affect the resilience of long-term boreal peatland C stocks by lowering the water table (WT) that can halt or even reverse the C accumulation in boreal peatlands (Limpens et al., 2008;Dise, 2009;Frolking et al., 2011).To maintain and protect the C sequestration potentials of boreal peatlands, we need an improved predictive capacity of how these C stocks would behave under future drier and warmer climates.However, boreal peatland C processes are currently underrepresented in global C models largely due to an inadequate simulation of hydrologic feedbacks to C cycles (St-Hilaire et al., 2010;Sulman et al., 2012).This can be overcome by integrating interactions between eco-hydrology of peatland vegetation and peat biogeochemistry into finerresolution process models that can eventually be scaled up into larger spatial-and temporal-scale C models (Waddington et al., 2015).
The hydrologic feedbacks to boreal peatland C processes are largely mediated by water table depth (WTD) variation and its effects on peat-microbe-plant-atmosphere exchanges of C, energy, water and nutrients (Grant et al., 2012).WTD drawdown can affect net ecosystem productivity (NEP) of boreal peatlands through its effects on ecosystem respiration (R e ) and gross primary productivity (GPP).A receding WT can cause peat pore drainage that enhances microbial O 2 availability, energy yields, growth and decomposition and hence increases R e (Sulman et al., 2009(Sulman et al., , 2010;;Cai et al., 2010;Flanagan and Syed, 2011;Peichl et al., 2014).The rate of increase in R e due to the WTD drawdown may vary with peat moisture retention and quality of peat-forming substrates (Preston et al., 2012).For instance, peats with low moisture retention capacities exhibit more rapid pore drainage than those with high moisture retention, thus causing more increase in R e for similar WTD drawdowns (Parmentier et al., 2009;Sulman et al., 2009Sulman et al., , 2010;;Cai et al., 2010).Peats formed from Sphagnum mosses degrade at rates slower than those formed from remains of vascular plants (Moore and Basiliko, 2006).So for similar WTD drawdowns, moss peats would experience less increase in microbial decomposition and hence R e than would sedge, reed, or woody peats (Updegraff et al., 1995).Continued WTD drawdown can also cause near-surface peat desiccation from inadequate recharge through capillary rise from deeper WT.The desiccation of near-surface shallow peat layers can cause a reduction in microbial decomposition that can partially or fully offset the increased decomposition in the deeper peat layers, thereby yielding indistinct net effects of WTD drawdown on R e (Dimitrov et al., 2010a).
The interactions between WTD and GPP vary across peatlands depending upon peat-forming vegetation.For instance, increased aeration due to WTD drawdown enhances root O 2 availability and growth in vascular plants (Lieffers and Rothwell, 1987;Murphy et al., 2009).Enhanced root growth is also associated with greater root nutrient availability and uptake from more rapid mineralization facilitated by greater microbial energy yields, growth and decomposition under a deeper WT (Choi et al., 2007).Greater root nutrient uptake increases the rate of vascular CO 2 fixation and hence GPP (Sulman et al., 2009(Sulman et al., , 2010;;Cai et al., 2010;Flanagan and Syed, 2011;Peichl et al., 2014).WTD drawdown, however, does not affect the nonvascular (e.g., moss) GPP in the same way it does the vascular GPP (Lafleur et al., 2005).Nonvascular plants mostly depend upon the water available for uptake in the near-surface shallow peat layers (Dimitrov et al., 2011).These layers can drain quickly with receding WT and thus have to depend on moisture supply through capillary rise from deeper WT (Dimitrov et al., 2011;Peichl et al., 2014).If recharge through the capillary rise is not adequate, near-surface peat desiccation occurs which slows moss water uptake, causes eventual drying of mosses and reduces moss GPP (Lafleur et al., 2005;Riutta 2008;Sonnentag et al., 2010;Sulman et al., 2010;Dimitrov et al., 2011;Kuiper et al., 2014;Peichl et al., 2014).Near-surface peat desiccation also suppresses vascular root water uptake from the desiccated layers (Lafleur et al., 2005;Dimitrov et al., 2011).But enhanced root growth and elongation facilitated by improved O 2 status in the newly aerated deeper peat layers under a deeper WT enables vascular roots to take up water from wetter deeper layers (Dimitrov et al., 2011).If deeper root water uptake offsets the reduction in water uptake from desiccated near-surface layers, vascular transpiration (T ), canopy stomatal conductance (g c ), and hence GPP are sustained under a deeper WT (Dimitrov et al., 2011).But if the WT falls below a certain threshold level under which deeper root water uptake can no longer sustain vascular T , g c and hence vascular GPP declines (Lafleur et al., 2005;Wu et al., 2010).
WTD variation can thus affect boreal peatland NEP through its effects on peat moisture and aeration and consequent root and microbial oxidation-reduction reactions and energy yields.To predict how boreal peatlands would behave under future drier and warmer climates, a peatland C model thus needs to simulate WTD dynamics that determine the boundary between aerobic and anaerobic zones and control peat biogeochemistry.However, most of the current processbased peatland C models either do not simulate a continuous anaerobic zone below a prognostic WT (e.g., Baker et al., 2008;Schaefer et al., 2008;Tian et al., 2010) or do not simulate peat saturation since any water in excess of field capacity is drained in those models (e.g., Gerten et al., 2004;Krinner et al., 2005;Weng and Luo, 2008).Moreover, instead of explicitly simulating the above-described hydrological and biological interactions between peat aeration and biogeochemistry, most of those models use scalar functions of soil moisture contents to inhibit R e and GPP under low or high moisture conditions (e.g., Frolking et al., 2002;Zhang et al., 2002;Bond-Lamberty et al., 2007;St-Hilaire et al., 2010;Sulman et al., 2012).Consequently, those peatland C models could not simulate declines in GPP and R e due to shallow WT while simulating WTD effects on CO 2 exchange of peatlands across the northern US and Canada (Sulman et al., 2012).Furthermore, the approach of using scalar functions to simulate moisture limitations to GPP and R e requires site-specific parameterization of model algorithms, which reduces the scalability of those peatland C models.

Objective and rationale
In this study, we tested a process-based ecosystem model, ecosys, against eddy covariance (EC) net CO 2 fluxes measured over a drying period from 2004 to 2009 in a western Canadian boreal fen peatland in Alberta, Canada (termed as WPL hereafter) (Syed et al., 2006;Flanagan and Syed, 2011).The objective was to test whether the coupling of a dynamic WTD that arises from vertical and lateral water fluxes as a function of a soil moisture retention scheme with (1) oxygen transport, (2) microbial and root oxidation-reduction reactions and energy yields, and (3) root, microbial and plant growth and uptake within a soil-plant-microbe-atmosphere water, C and nutrient (nitrogen, phosphorus) scheme in an ecosystem process model could explain underlying processes that govern hydrological effects on net CO 2 exchange of a northern boreal fen under contrasting WTD conditions (e.g., shallower vs. deeper WTD).This study would reconcile our knowledge on the feedback mechanisms among hydrology, ecophysiology and biogeochemistry of peatlands, which is predominantly based upon inferences drawn from EC gapfilled values that include empirically modeled estimates.It would also provide us with a better insight into -and an improved predictive capacity of -how carbon deposits in northern boreal peatlands would behave under changing climates.Rigorous site-scale testing of coupled eco-hydrology and biogeochemistry algorithms in ecosystem process models such as ecosys would also provide us with important insights on how to improve large-scale representation of these processes into next-generation land surface models.

Hypotheses
In an EC study, Flanagan and Syed (2011) found no net effect of a drier-weather-driven WTD drawdown on NEP at the WPL over 2004-2009.From the regressions of EC-derived GPP and R e on site-measured WTD, they inferred that the absence of a net WTD effect on NEP was caused by similar increases in GPP and R e with WTD drawdown.We hypothesized that coupled eco-hydrology and biogeochemistry algorithms in ecosys would be able to simulate and explain underlying mechanisms of these effects of WTD drawdown on GPP and R e and hence NEP at the WPL.We tested the following four central hypotheses: 1. WTD drawdown would increase R e of the northern fen at the WPL.This effect of WTD drawdown on R e would be modeled by the simulation of peat pore drainage and improved peat aeration that would increase the energy yields from aerobic microbial decomposition and hence would increase R e .
2. WTD drawdown would increase GPP of the northern fen at the WPL.This effect of WTD drawdown on GPP would be modeled by simulating enhanced microbial activity due to WTD drawdown that would cause more rapid nutrient mineralization and greater root nutrient availability and uptake, greater leaf nutrient concentrations and hence increased GPP.
3. An increase in R e with WTD drawdown (hypothesis 1) would cease should WTD fall below a threshold depth.This threshold WTD effect on R e would be modeled by simulating inadequate recharge of the near-surface peat layers through capillary rise from the deeper WTD below the threshold level that would cause desiccation of those layers.Drying of near-surface peat layers and the surface residue would reduce near-surface and surface peat respiration that would partially offset the increase in deeper peat respiration due to aeration.
4. The net effect of threshold WTD on GPP would be driven by the balance between how WTD would affect vascular vs. nonvascular GPP.This threshold WTD effect on GPP would be modeled by simulating vascular vs. nonvascular water relations under deeper WTD below the threshold level.Near-surface peat desiccation in hypothesis 3 would reduce peat water potential and hydraulic conductivity and hence vascular and nonvascular water uptake from desiccated near-surface layers.Since nonvascular mosses depend mainly on nearsurface peat layers for moisture supply, reduction in moss water uptake would cause a reduction in moss water potential and hence moss GPP.By contrast, the suppression of vascular root water uptake from desiccated near-surface layers under a deeper WT would be offset by increased deeper root water uptake from newly aerated deeper peat layers with higher water potentials that would sustain vascular canopy water potential (ψ c ), g c and GPP.

Model development
Ecosys is a process-based ecosystem model which simulated 3-D water, energy, carbon and nutrient (nitrogen, phosphorus) cycles in different peatlands (Dimitrov et al., 2011;Grant et al., 2012;Sulman et al., 2012;Mezbahuddin et al  2014 , 2015, 2016).Ecosys algorithms that govern the modeled effects of WTD variations on peatland net CO 2 exchange are described below.These algorithms in ecosys are derived from published independent basic research which describes eco-hydrological and biogeochemical mechanisms that govern the carbon, nutrient (N, P), water and energy balance of a typical boreal peatland ecosystem (Fig. 1).Ecosys algorithms which are related to our hypotheses are depicted as a flowchart (Fig. 1) and cited as equations within the text.These equations are also listed in appendices A to D in the Supplement.These site-independent basic ecosystem process algorithms in ecosys are thus not parameterized for each peatland site.Instead the coupled algorithms are fed with peatland-specific measurable soil, weather, vegetation and management inputs to simulate the C, nutrient (N, P), water and energy balance of a particular peatland ecosystem.

Water table depth
The WTD in ecosys is calculated at the end of each time step as the depth to the top of the saturated zone, below which air-filled porosity is 0 (Eq.SD32).It is the depth at which lateral water flux is in equilibrium with the difference between vertical influxes (precipitation) and effluxes (evapotranspiration). Lateral water transfer between modeled grid cells in ecosys and the adjacent ecosystem occurs to and from a set external WTD (WTD x ) over a set distance (L t ) (Fig. 2).The WTD x represents average watershed WTD with reference to average hummock surface.The WTD in ecosys is thus not prescribed but rather controls, and is controlled by, lateral and vertical surface and subsurface water fluxes (Eqs.SD1-SD31).More detail about how peatland WTD, vertical and lateral soil water flow and soil moisture retention are modeled in ecosys can be found in Dimitrov et al. (2010b) and Mezbahuddin et al. (2015Mezbahuddin et al. ( , 2016)).

Heterotrophic respiration and WTD
WTD fluctuation in ecosys determines the boundary between and the extent of aerobic vs. anaerobic soil zones.So, WTD fluctuation affects ecosys's algorithms of organic oxidation-reduction transformations and microbial energy  yields, which drive microbial growth and substrate decomposition and uptake (Fig. 1) (Eqs.SA1-SA30).Organic transformations in ecosys occur in a residue layer and in each of the user-defined soil layers within five organicmatter-microbe complexes, i.e., coarse woody litter, fine nonwoody litter, animal manure, particulate organic C and humus (Fig. 1).Each of the complexes has three decomposition substrates, i.e., solid organic C, sorbed organic C and microbial residue C; the decomposition agent, i.e., microbial biomass; and the decomposition product, i.e., dissolved organic C (DOC) (Fig. 1).The rates of the decomposition and resulting DOC production in each of the complexes are a first-order function of the fraction of substrate colonized by active biomasses (M) of diverse microbial functional types (MFTs).The MFTs in ecosys are obligate aerobes (bacteria and fungi), facultative anaerobes (denitrifiers), obligate anaerobes (fermenters), heterotrophic (acetotrophic) and autotrophic (hydrogenotrophic) methanogens, and aerobic and anaerobic heterotrophic diazotrophs (nonsymbiotic N 2 fixers) (Fig. 1) (Eqs.SA1-SA2, SA4).Biomass (M) growth of each of the MFTs (Eq.SA25a) is calculated from its DOC uptake (Fig. 1) (Eq.SA21).The rate of M growth is driven by energy yield from growth respiration (R g ) (Eq.SA20) that is calculated by subtracting maintenance respiration (R m ) (Eq.SA18) from heterotrophic respiration (R h ) (Eq.SA11 ] slows O 2 uptake (Eq.SA17) and hence R h (Eq.SA14), R g (Eq.SA20) and the growth of M (Eq.SA25).
Lower M slows decomposition of organic C (Eqs.SA1-SA2) and production of DOC, which further slows R h (Eq.SA13), R g and the growth of M (Fig. 1).Although some MFTs can sustain DOC oxidation by reducing alternative electron acceptors (e.g., methanogens reducing acetate or CO 2 to CH 4 and denitrifiers reducing NO x to N 2 O or N 2 ), lower energy yields from these reactions reduce R g (Eq.SA21) and hence M growth, organic C decomposition and subsequent DOC production (Fig. 1).Slower decomposition of organic C under low [O 2s ] also causes slower decomposition of organic nitrogen (N) and phosphorus (P) (Eq.SA7) and the production of dissolved organic nitrogen (DON) and phosphorus (DOP), which causes slower uptake of microbial N and P (Eq.SA22) and growth of M (Eq.SA29) (Fig. 1).Slower M growth causes slower mineralization (Eq.SA26) and hence lowers aqueous concentrations of NH + 4 , NO − 3 and H 2 PO − 4 (Fig. 1).WTD drawdown can increase θ g that results in greater D g (Eq.SD44) and more rapid gaseous O 2 transport.A consequent rise in [O 2s ] increases O 2 uptake (Eq.SA17) and R h (Eq.SA14), R g (Eq.SA20) and the growth of M(Eq.SA25).Larger M hastens the decomposition of organic C (Eqs.SA1-SA2) and the production of DOC which further hastens R h (Eq.SA13), R g and the growth of M.More rapid decomposition of organic C under adequate [O 2s ] in this period also causes more rapid decomposition of organic N and P (Eq.SA7) and production of DON and DOP, which increases the uptake of microbial N and P (Eq.SA22) and the growth of M (Eq.SA29) (Fig. 1).Rapid M growth causes rapid mineralization (Eq.SA26) and hence greater aqueous concentrations of NH + 4 , NO − 3 and H 2 PO − 4 (Fig. 1).When WTD recedes below a certain threshold level, capillary rise from the WT can no longer support adequate recharge of the near-surface peat layers and the surface litter (Eqs.SD9, SD12).It causes desiccation of the residue and the near-surface peat layers, thereby causing a reduction in water potential (ψ s ) and an increase in aqueous microbial concentrations ([M]) in each of these layers (Eq.SA15).Increased [M] caused by the peat desiccation reduces microbial access to the substrate for decomposition in each of the desiccated layers and reduces R h (Eq.SA13).The reduction in R h is calculated in ecosys from competitive inhibition of microbial exoenzymes with increasing concentrations (Eq.SA4) (Lizama and Suzuki, 1991).

WTD effects on vascular gross primary productivity
Ecosys simulates effects of WTD variation on vascular GPP from WTD variation effects on root O 2 and nutrient availability and root growth and uptake.Root O 2 and nutrient up-take in ecosys are coupled with a hydraulically driven soilplant-atmosphere water scheme.Root growth in each vascular plant population in ecosys is calculated from its assimilation of the nonstructural C product of CO 2 fixation (σ C ) (Eq.SC20).Assimilation is driven by R g (Eq.SC17) remaining after subtracting R m (Eq.SC16) from autotrophic respiration (R a ) (Eq.SC13) driven by oxidation of σ C (Eq. SC14 ] slows root O 2 uptake (Eq.SC14c-d) and hence R a (Eq.SC14b), R g (Eq.SC17) and root growth (Eq.SC20b) and root N and P uptake (Eqs.SC23b, d, f) (Fig. 1).Root N and P uptake under a shallow WT is further slowed by reductions in aqueous concentrations of NH + 4 , NO − 3 and H 2 PO − 4 (Eqs.SC23a, c, e) from slower mineralization of organic N and P (Fig. 1).Slower root N and P uptake reduces concentrations of nonstructural N and P products of root uptake (σ N and σ P ) with respect to that of σ C in leaves (Eq.SC11), thereby slowing CO 2 fixation (Eq.SC6) and GPP.
WTD drawdown facilitates rapid D g which allows root O 2 demand to be almost entirely met from [O 2s ] (Eqs.SC14c-d) and so enables more rapid root growth and N and P uptake (Eqs.SC23b, d, f).Increased root growth and nutrient uptake is further stimulated by increased aqueous concentrations of NH + 4 , NO − 3 and H 2 PO − 4 (Eqs.SC23a, c, e) from more rapid mineralization of organic N and P during deeper WT (Fig. 1).Greater root N and P uptake increases concentrations of σ N and σ P with respect to σ C in leaves (Eq.SC11), thereby facilitating rapid CO 2 fixation (Eq.SC6) and GPP. when the WT falls below a certain threshold, inadequate capillary rise (Eqs.SA9, SA12) from deeper WT causes near-surface peat desiccation that reduces ψ s and raises soil hydraulic resistance ( s ) (Eq.SB9), thereby forcing lower root water uptake (U w ) from desiccated layers (Fig. 1) (Eq.SB6).However, deeper rooting facilitated by increased [O 2s ] under a deeper WT can sustain U w (Eq.SB6) from wetter deeper peat layers with higher ψ s and lower s (Eq.SB9).If U w from the deeper wetter layers cannot offset the suppression in U w from desiccated near-surface layers, the resultant net decrease in U w causes a reduction in root (ψ r ), canopy and turgor potentials (ψ t ) (Eq.SB4) and hence g c (Eq. SB2b) in ecosys when equilibrating U w with transpiration (T ) (Eq.SB14).Lower g c reduces CO 2 diffusion into the leaves, thereby reducing CO 2 fixation (Eq.SC6) and GPP (Eq.SC1) (Fig. 1).

WTD effects on nonvascular gross primary productivity
Ecosys simulates nonvascular plants (e.g., mosses) as tiny plants with no stomatal regulations that grow on modeled hummock and hollow grid cells (Dimitrov et al., 2011).Model input for moss population is usually larger and hence intraspecific competition for lights and nutrients is greater so that individual moss plant and moss belowground growth (i.e., root like structures for water and nutrient uptake) are smaller (Eq.SC21b).Shallower belowground growth of simulated mosses in ecosys means that the water uptake of mosses is mostly confined to the near-surface peat layers.when the WT deepens past a threshold level, inadequate capillary rise (Eqs.SD9, SD12) causes near-surface peat desiccation, thereby reducing ψ s and increasing s (Eq.SB9) of those layers (Fig. 1).Reduced ψ s along with increased s causes a reduction in moss ψ c while equilibrating moss evaporation with moss U w (Eq.SB6).Reduced moss ψ c causes a reduction in moss carboxylation rate (Eqs.SC3, SC6a) and moss GPP (Fig. 1) (Eq.SC1).

Study site
Eco-hydrology and biogeochemistry algorithms in ecosys were tested in this study against measurements of WTD and ecosystem net CO 2 fluxes at a flux station of the FLUXNET-Canada Research Network established at the WPL (latitude: 54.95 • N; longitude: 112.47 • W).The study site is a moderately nutrient-rich treed fen peatland within the central mixed-wood subregion of boreal alberta, Canada.Peat depth around the flux station was about 2 m.This peatland is dominated by stunted trees of black spruce (Picea mariana) and tamarack (Larix laricina) with an average canopy height of 3 m.A high abundance of a shrub species Betula pumila (dwarf birch) and the presence of a wide range of mosses, e.g., Sphagnum spp., feather moss and brown moss characterize the understory vegetation of WPL.The topographic, climatic, edaphic and vegetative characteristics of this site were described in more details by Syed et al. (2006).

Field data sets
Ecosys model inputs of half-hourly weather variables, i.e., incoming shortwave and longwave radiation, air temperature (T a ), wind speed, precipitation and relative humidity during 2003-2009, were measured by Syed et al. (2006) and Flanagan and Syed (2011) at the micrometeorological station installed at the WPL.To test the adequacy of WTD simulation in ecosys, modeled outputs of hourly WTD were tested against WTD measured at the WPL with respect to average hummock surface by Flanagan andSyed (2011) (Mezbahuddin et al., 2016).To examine how well ecosys simulated the net ecosystem CO 2 exchange at the WPL, we tested hourly modeled net ecosystem CO 2 fluxes against hourly averaged measurements (average of two half-hourly measurements) collected by Syed et al. (2006) and Flanagan and Syed (2011) by using an EC micrometeorological approach.Each of these EC-measured net CO 2 fluxes consisted of an eddy flux and a storage flux (Syed et al., 2006).Erroneous flux measurements due to stable air conditions were screened out with the use of a minimum friction velocity (u * ) threshold of 0.15 m s −1 (Syed et al., 2006).The net CO 2 fluxes that survived the quality control were used to derive half-hourly R e (nighttime net CO 2 fluxes) and GPP (daytime net CO 2 fluxes -R e ) (Barr et al., 2004;Syed et al., 2006).To derive daily, seasonal and annual estimates of GPP and R e , the data gaps resulting from the quality control were filled based on empirical relationships between soil temperature at a shallow depth (0.05 m) and measured half-hourly R e and between incoming shortwave radiation and measured half-hourly GPP using 15day moving windows.The gap-filled R e and GPP were then summed up for each half hour to fill NEP data gaps to derive daily, seasonal and annual estimates of EC gap-filled NEP (Barr et al., 2004;Syed et al., 2006).Soil CO 2 fluxes measured by automated chambers can provide a valuable supplement to EC CO 2 fluxes in testing modeled respiration by providing more continuous measurements than EC.So, we tested our modeled outputs against half-hourly automated chamber measurements by Cai et al. (2010) at the WPL.These CO 2 flux measurements were carried out during ice-free periods (May-October) of 2005 and 2006 over both hummocks and hollows by using a total of nine non-steady-state automatic transparent chambers (Cai et al., 2010).Along with soil respiration these chamber CO 2 fluxes included fixation and autotrophic respiration from dwarf shrubs, herbs and mosses (Cai et al., 2010).
Modeled WTD effects on peatland biogeochemistry and hence on peatland nutrient and carbon cycling were also corroborated against leaf nitrogen concentrations, foliar N-to-P ratios, N mineralization, and rooting depths biometrically measured at either our site or at sites that had similar peat substrates, hydrology and/or plant functional types.Needles of black spruce and tamarack and leaves of dwarf birch were sampled over our study site during midsummer of 2004 for foliar nutrient content analyses (Syed et al., 2006).Leaf nitrogen contents were analyzed on N 2 gas that was generated from the reduction of dried leaf tissues in an elemental analyzer and quantified using a gas isotope ratio mass spectrometer (Syed et al., 2006).Leaf phosphorus contents were analyzed on black spruce and tamarack needle leaf tissues that were dry-ashed and then digested using a dilute HNO 3 and HCl mixture and then quantified using an inductively coupled plasma (ICP) spectroscopic analysis technique (Syed et al., 2006).

Model run
An ecosys model run to simulate WTD effects on net CO 2 exchange of WPL had a hummock and a hollow grid cell that exchanged water, heat, carbon and nutrients (N, P) between them and with surrounding vertical and lateral boundaries (Fig. 2).The hollow grid cell had a near-surface peat layer that was 0.3 m thinner than the hummock cell representing a hummock-hollow surface difference of 0.3 m observed in the field (Long, 2008) (Fig. 2).Any depth with respect to the modeled hollow surface would thus be 0.3 m shallower than the depth with respect to the modeled hummock surface.
Peat organic and chemical properties at different depths of the WPL were represented in ecosys by inputs from measurements either at the site (e.g., Syed et al., 2006;Flanagan and Syed, 2011) or at similar nearby sites (e.g., Rippy and Nelson, 2007) (Fig. 2).Ecosys was run for a spin-up period of 1961-2002 under repeating 7-year sequences of hourly weather data (shortwave and longwave radiation, air temperature, wind speed, humidity and precipitation) recorded at the site from 2003 to 2009.There was a drying trend observed from 2003 to 2009 due to diminishing precipitation that caused WTD drawdown in the watershed in which WPL is located, which lowered the WT of this fen peatland (Flanagan and Syed, 2011).To accommodate the gradual drying effects of catchment hydrology on modeled fen peatland WTD, we set the WTD x at different levels based on the annual wetness of weather, e.g., shallow, intermediate, and deep (WTD x = 0.19, 0.35 and 0.72 m below the hummock surface or 0.11 m above and 0.05 and 0.42 m below the hollow surface) (Fig. 2).There was no exchange of water through the lower model boundary to represent the presence of nearly impermeable clay sediment underlying the peat (Syed et al., 2006) (Fig. 2).Variations in peat surface with WTD variations, which are an important hydrologic self-regulation of boreal peatlands (Dise, 2009), were not represented in this version of ecosys.
At the start of the spin-up run, the hummock grid cell was seeded with an evergreen needle leaf and a deciduous needle leaf overstory plant functional type (PFT) to represent the black spruce and tamarack trees at the WPL.The hollow grid cell was seeded with only the deciduous needle leaf overstory PFT since the black spruce trees at the WPL only grew on the raised areas.Each of the modeled hummock and the hollow was also seeded with a deciduous broadleaved vascular (to represent dwarf birch) and a nonvascular (to represent mosses) understory PFTs.The planting densities were such that the population densities of the black spruce, tamarack, dwarf birch and moss PFTs were 0.16, 0.14, 0.3, and 500 m −2 , respectively, at the end of the spin-up run so as to best represent field vegetation (Syed et al., 2006;Mezbahuddin et al., 2016).To include wetland adaptation, we used a root porosity (θ pr ) value of 0.1 for the two overstory PFTs and a higher θ pr value of 0.3 for the understory vascular PFT to represent better wetland adaptation in the understory than the overstory PFTs.These θ pr values were used in calculating root O 2 transport through aerenchyma (Eq.SD45) and did not change with waterlogging throughout the model run.These θ pr values were representatives of root porosities measured for various northern boreal peatland plant species (Cronk and Fennessy, 2001).Nonsymbiotic N 2 fixation through the association of cyanobacteria and mosses are also reported for Canadian boreal forests (Markham, 2009).This was represented in ecosys as N 2 fixation by nonsymbiotic heterotrophic diazotrophs (Eq.SA27) in the moss canopy.Further details about ecosys model setup to represent the hydrological, physical and ecological characteristics of WPL can be found in Mezbahuddin et al. (2016).
When the modeled ecosystem attained stable values of net ecosystem CO 2 exchange at the end of the spin-up run, we continued the spin-up run into a simulation run from 2003 to 2009 by using a real-time weather sequence.We tested our outputs from 2004 to 2009 of the simulation run against the available site measurements of WTD, net EC CO 2 fluxes and net chamber CO 2 fluxes over those years.

Model validation
To examine the adequacy of modeling WTD effects on canopy, root and soil CO 2 fluxes, which were summed for net ecosystem CO 2 exchange at the WPL, we spatially averaged hourly net CO 2 fluxes modeled over the hummock and the hollow to represent a 50 : 50 hummock-hollow ratio; they were regressed against hourly EC-measured net ecosystem CO 2 fluxes for each year from 2004 to 2009 with varying WTD.Each of these hourly EC-measured net ecosystem CO 2 fluxes used in these regressions is an average of two half-hourly net CO 2 fluxes measured at a friction velocity (u * ) greater than 0.15 m s −1 that survived quality control procedure (Sect.3.2.2).Model performance was evaluated from regression intercepts (a → 0), slopes (b → 1), coefficients of determination (R 2 → 1), and root means squares for errors (RMSE → 0) for each study year to test whether there was any systematic divergence between the modeled and EC-measured CO 2 fluxes.
Similar regressions were performed between modeled and automated chamber-measured net CO 2 fluxes for ice-free periods (May-October) of 2005 and 2006 to further test the robustness of modeled soil respiration under contrasting WTD conditions.Each of the half-hourly measured chamber net CO 2 fluxes included soil respiration and fixation and autotrophic respiration from understory vegetation (e.g., shrubs, herbs and mosses).So, we combined modeled soil respiration with modeled fixation and autotrophic respiration from understory PFTs for comparison against these chambermeasured net CO 2 fluxes.We also averaged net CO 2 flux measurements from all of the nine chambers for each half hour to accommodate the variations in those fluxes due to microtopography (e.g., hummock vs. hollow).Two halfhourly averaged values of net CO 2 fluxes were then aver-aged again to obtain hourly mean net chamber CO 2 fluxes for comparison against modeled hourly sums of soil and understory fluxes averaged over modeled hummock and hollow.Model performance was evaluated from regression intercepts (a → 0), slopes (b → 1), coefficients of determination (R 2 → 1), and root means squares for errors (RMSE → 0) for each of 2005 and 2006.

Sensitivity of modeled peatland CO 2 exchange to artificial drainage
Large areas of northern boreal peatlands in Canada have been drained primarily for increased forest and agricultural production since plant productivity in pristine peatlands are known to be constrained by shallow WTD (Choi et al., 2007).
Drainage and resultant WTD drawdown can affect both GPP and R e on a short-term basis and the vegetation composition on a longer timescale, thereby changing overall net CO 2 exchange trajectories of a peatland.To predict the short-term effects of drainage on WTD and hence ecosystem net CO 2 exchange of WPL, we extended our simulation run into a projection run consisting of two 7-year cycles by using repeated weather sequences of 2003-2009.While doing so, we forced a stepwise drawdown in WTD x by 1.0 and 2.0 m from that used in spin-up and simulation runs (Fig. 2) in the first (drainage cycle 1) and the second (drainage cycle 2) cycle, respectively.This projection run would give us a further insight into how the northern boreal peatland of western Canada would be affected by further WTD drawdown as a result of drier and warmer weather as well as a disturbance such as drainage.It would also provide us with a test of how sensitive the modeled C processes were to the changes in model lateral boundary condition as defined by WTD x in ecosys.

Model performance in simulating diurnal variations in ecosystem net CO 2 fluxes
Variations in precipitation can cause change in WTD and consequent variation in diurnal net CO 2 exchange across years.Ecosys simulated hourly EC-measured net CO 2 fluxes well from 2004 to 2009 with varying precipitation (Table 1).On a year-to-year basis, regressions of hourly modeled vs. EC-measured net ecosystem CO 2 fluxes gave intercepts within 0.1µmol m −2 s −1 of 0 and slopes within 0.1 of 1, indicating minimal bias in modeled outputs during each year from 2004 to 2008 (Table 1a).On a growing season (May-August) basis, regressions of modeled on measured hourly net CO 2 fluxes yielded larger positive intercepts from 2004 to 2009 (Table 1b).The larger intercepts were predominantly caused by model overestimation of growing season daytime CO 2 fluxes.This overestimation was offset by model overestimation of nighttime CO 2 fluxes during the winter, thus yielding smaller intercepts from throughout-theyear regressions of modeled vs. EC-measured fluxes (Table 1a vs. b).Values for coefficients of determination (R 2 ) were ∼ 0.8 (P < 0.001) for all years from both throughoutthe-year and growing season regressions (Table 1a, b).RM-SEs were < 2.0 and ∼ 2.5 µmol m −2 s −1 for whole-year regressions from 2004 to 2008 (Table 1a) and for growing season regressions from 2004 to 2009 (Table 1b), respectively.Much of the variations in EC-measured CO 2 fluxes that were not explained by the modeled fluxes could be attributed to a random error of ∼ 20 % in EC methodology (Wesely and Hart, 1985).This attribution was further corroborated by root mean squares for random errors (RMSRE) in EC measurements, calculated for forests with similar CO 2 fluxes from Richardson et al. (2006), which were similar to RMSE (Table 1a, b).The similar values of RMSE and RMSRE also indicated that further constraint in model testing could not be achieved without further precision in EC measurements.
Regressions of modeled vs. chamber-measured net CO 2 fluxes gave R 2 of ∼ 0.7 for ice-free periods (May-October) of 2005, and 2006 indicated that the variations in soil respiration and the fixation and aboveground autotrophic respiration due to WTD drawdown were modeled well (Table 1c).Smaller intercepts from those regressions meant lower model biases in simulating soil and understory CO 2 fluxes under a deepening WT (Table 1c).Although the slope was within 0.1 of one in 2005, it was a bit smaller in 2006.It indicated that the modeled soil and understory net CO 2 fluxes were smaller than the chamber-measured fluxes in 2006 (Table 1c).This was because some of the nighttime chamber fluxes in warmer nights of summer 2006 were as large as corresponding ECmeasured net ecosystem CO 2 fluxes, which could not be modeled to their full extent.An RMSE lower than RMSRE meant that the errors in modeling soil and understory CO 2 fluxes were within the limit of random errors due to chamber measurements (Table 1c).It further indicated the robustness of modeled outputs for soil and understory CO 2 fluxes under different WTD conditions (e.g., shallower in 2005 vs. deeper in 2006) (Table 1c).

Seasonality in WTD and net ecosystem CO 2 exchange
Ecosys simulated the seasonal and interannual variations in WTD well from 2004 to 2009 at the WPL (Fig. 3b, d, f, h, j, l) (Mezbahuddin et al., 2016).Seasonality in net CO 2 exchange at the WPL was predominantly governed by that in temperature, which controlled the seasonality in phenology and GPP as well as that in R e .Ecosys simulated the seasonality in phenology and hence GPP and R e well during a gradual growing season WTD drawdown from 2004 to 2009, which was apparent by good agreements between modeled vs. EC gap-filled daily NEP (Fig. 3a, c, e, g, i, k) and modeled vs. EC-measured hourly net CO 2 fluxes ( more negative than the EC gap-filled NEP, indicating larger modeled CO 2 effluxes than EC gap-filled fluxes during the winter (Fig. 3).The onset of photosynthesis at the WPL varied interannually depending upon spring temperature, which was also modeled well by ecosys.For instance, ecosys modeled a smaller early growing season (May) GPP and hence NEP in 2004, with a cooler spring than 2005, which was also apparent in daily EC gap-filled NEP (Fig. 3a vs. c).a, c, e, g, i, k) Three-day moving averages of modeled and EC gap-filled (Flanagan and Syed, 2011) NEP and (b, d, f, h, j, l) hourly modeled and half-hourly measured (Syed et al., 2006;Cai et al., 2010;Long et al., 2010;Flanagan and Syed, 2011) WTD from 2004 to 2009 at a western Canadian fen peatland.A positive NEP means that the ecosystem is a carbon sink, and a negative NEP means that the ecosystem is a carbon source.A negative WTD represents a depth below hummock/hollow surface, and a positive WTD represents a depth above hummock/hollow surface.(Syed et al., 2006;Cai et al.;2010, Long et al.;2010, Flanagan andSyed, 2011)   Beside WTD, temperature variation also profoundly affected ecosystem net CO 2 exchange at the WPL.For a given WTD condition, warmer weather raised R e at the WPL as represented by larger modeled and measured ecosystem, soil and understory nighttime CO 2 fluxes in warmer nights than those in cooler nights (Figs.4b and 5a-b).However, for a given temperature modeled and EC-measured nighttime ecosystem CO 2 fluxes and modeled and chamber-measured nighttime soil and understory CO 2 fluxes were larger under deeper WT conditions in 2006 and 2008 than under shallower WT conditions in 2005 (denoted by the grey arrows in Figs.4b and 5a-b).The nighttime CO 2 fluxes that were larger under deeper WTD than under shallower WTD with similar temperatures showed net a WTD drawdown effect on R e (separated from the temperature effect) and hence on NEP.
The degree of R e stimulation due to warming was also influenced by WTD at the WPL.The warming events in early to mid-August of 2006 and 2008, when the WT was deeper than in late July of 2005, caused gradual increases in modeled and EC-and chamber-measured nighttime ecosystem, soil and understory CO 2 effluxes (Fig. 6h, i, k, l).This R e stimulation due to warming under a deeper WT contributed to declines in modeled and EC gap-filled July-August NEP in 2006 and2008 (Fig. 3e, i).A lack of similar stimulation in R e with warming under a shallower WT in 2005 did not yield a similarly evident stimulation of either modeled or measured ecosystem, soil and understory nighttime CO 2 effluxes (Fig. 6g, j vs. h, i, k, l), which resulted in the absence of a decline in July-August NEP in 2005 as occurred in 2006 and 2008 (Fig. 3c vs. e, i).Greater warming-driven R e stim-  (d-f) hourly modeled and half-hourly observed (Syed et al., 2006;Cai et al., 2010;Long et al., 2010;Flanagan and Syed, 2011) WTD; (g-i) half-hourly EC-measured and gap-filled (Flanagan and Syed, 2011) and hourly modeled ecosystem net CO 2 fluxes; (j-l) half-hourly automated chamber-measured (Cai et al., 2010)

Interannual variations in WTD and net ecosystem productivity
The effects of WTD drawdown on modeled and EC gapfilled diurnal net ecosystem CO 2 exchange also contributed to the effects of interannual variation in WTD on that of NEP.Ecosys simulated the observed gradual drawdown of average growing season (May-August) WTD well from 2004 to 2009 (Fig. 7d) (Mezbahuddin et al., 2016).A small WTD drawdown caused a large increase in both modeled and ECderived growing season GPP from 2004 to 2005 (Fig. 7b).This increase in GPP was also contributed by a larger GPP in May 2005, which was warmer than May 2004.This small WTD drawdown, however, did not raise either modeled or EC-derived growing season R e from 2004 to 2005 (Fig. 7c).June and July of 2005 were cooler than 2004 by over 2 • C, which caused cooler soil that reduced R e in 2005.The reduction in R e in June-July 2005 due to cooler soil more than fully offset the increase in R e due to the small WTD drawdown and resulted in modeled and EC-derived growing season R e values that were smaller in 2005 than in 2004 (Fig. 7c).Larger GPP along with smaller R e in 2005 con-tributed to larger modeled and EC gap-filled growing season NEP in 2005 than in 2004 (Fig. 7a).WTD drawdown from 2005 to 2006 raised both modeled and EC-derived growing season GPP and R e (Fig. 7b, c, d).A warmer growing season in 2006 caused warmer soil that further contributed to the increase in modeled and EC-derived growing season R e from 2005 with a shallower WT to 2006 with a deeper WT (Figs. 5,6 and 7c,d).An increase in growing season R e that was greater than the increase in growing season GPP caused a decline in modeled and EC-derived growing season NEP from 2005 to 2006 (Fig. 7a).Continued growing season WTD drawdown from 2006 to 2008 caused similar increases in modeled growing season GPP and R e and hence caused no significant change in modeled growing season NEP (Fig. 7a-d).With this continued WTD drawdown from 2006 to 2008, however, EC-derived growing season GPP increased more than EC-derived growing season R e , which resulted in a larger EC-derived NEP in the growing season of 2008 than in 2006 (Fig. 7a-d).A further growing season WTD drawdown from 2008 to 2009 caused reductions in both modeled and EC-derived growing season GPP and R e (Fig. 7a-d).Reductions in GPP and R e from 2008 to 2009 were also contributed by lower T a in 2009 than in 2008, which caused cooler canopies and soil (Fig. 7b-d).The reduction in EC-derived growing season GPP was larger than that in EC-derived growing season R e , thereby causing www.biogeosciences.net/14/5507/2017/Biogeosciences, 14, 5507-5531, 2017 a decrease in growing season EC gap-filled NEP from 2008 to 2009 (Fig. 7a-c).By contrast, the reduction in modeled growing season GPP was smaller than the reduction in modeled R e that yielded an increase in modeled growing season NEP from 2008 to 2009 (Fig. 7a-c).
Modeled and EC-derived growing season GPP and R e in 2009 were larger than those in 2004 despite similar mean T a in those years (Fig. 7a-d).It suggested that increases in growing season GPP and R e from 2004 to 2009 was a net effect of the deepening of average growing season WT (Fig. 7a-d).It was further corroborated by polynomial regressions of modeled growing season GPP and R e against modeled average growing season WTD and by similar regressions of EC-derived growing season GPP and R e against site-measured average growing season WTD (Fig. 8b, c).These relationships showed that there were increases in modeled and EC-derived growing season GPP and R e with a deepening of the growing season WT from 2004 to 2008, after which further WTD drawdown in 2009 started to cause slight declines in both GPP and R e (Fig. 8b, c).Neither modeled nor EC gap-filled growing season NEP yielded significant regressions when regressed against modeled and measured growing season WTD, respectively (Fig. 8a).It indicated that similar increases in modeled and EC-derived growing season GPP and R e with a deepening of the WT had no net effects of WTD drawdown on either modeled or ECderived growing season NEP (Figs. 7a-d and 8a).
Similar to the growing season trend, drawdown of both measured and modeled WTD averaged over the ice-free periods (May-October) from 2004 to 2008 stimulated modeled and EC-derived annual GPP and R e (Figs. 7f,g,h and 8e,f).Similar increases in both modeled and EC-derived annual GPP and R e with WTD drawdown had no net WTD effects on modeled and EC gap-filled annual NEP (Figs. 7e,  8d).Although modeled WTD effects on GPP, R e and hence NEP were corroborated well by EC-derived GPP, R e and EC gap-filled NEP, the modeled values for growing season and annual GPP and R e were consistently larger than the ECderived GPP and R e throughout the study period (Figs.7b, c,  f, g and 8b, c, e, f).
Increased GPP with WTD drawdown (Figs.7b, f and 8b, e) was modeled predominantly through increased root growth and uptake of nutrients and consequently improved leaf nutrient status and hence more rapid CO 2 fixation in vascular PFTs.Under a shallow WT during the growing season of 2004, roots in modeled black spruce and tamarack PFTs hardly grew below 0.35 m from the hummock surface.Modeled root densities of both black spruce and tamarack were higher by 2-3 orders of magnitude in the top 0.19 m of the hummock (data not shown).A WTD drawdown by 0.35 m  ) EC derived carbon balance components vs. observed WTD Modeled carbon balance components vs. modeled WTD from the growing season of 2004 to that of 2009 caused an increase in maximum modeled rooting depth in both PFTs (Table 2).Increased root growth in modeled vascular PFTs augmented root surface area for nutrient uptake under a deeper WT in the growing season of 2009 than in 2004.Increased root surface area along with increased nutrient availability due to more rapid mineralization with improved aeration as a result of WTD drawdown from 2004 to 2009 caused improved root nutrient uptake in modeled vascular PFTs.Increased root growth, nutrient availability and hence uptake due to WTD drawdown from the growing season of 2004 to that of 2009 caused an increase in modeled foliar N concentrations in black spruce, tamarack and dwarf birch PFTs driving the increases in GPP modeled over this period (Figs.7b, f and 8b, e) (Table 2).

Simulated drainage effects on WTD and NEP
Artificial drainage can drastically alter WTD in a peatland that can cause dramatic changes in peatland NEP by shifting the balance between GPP and R e .Projected growing season WT was deeper by ∼ 0.5 m and ∼ 0.55 m, respectively, from those in the real-time simulation in drainage cycles 1 and 2 in all the years from 2004 to 2009 (Fig. 9a).Modeled growing season GPP increased with drainage-induced WTD drawdown up to ∼ 0.5 m below the hollow surface (∼ 0.8 m below the hummock surface), below which GPP decreased (Fig. 9c, f).The WTD drawdown affected modeled vascular and nonvascular growing season GPP quite differently.Modeled growing season vascular GPP increased with WTD drawdown before it plateaued and eventually decreased when WTD fell below ∼ 0.6 m from the hollow surface (∼ 0.9 m below the hummock surface) (Fig. 10a, c, e).By contrast, modeled nonvascular growing season GPP continued to decrease with WTD drawdown below ∼ 0.1 m from the hollow surface (∼ 0.4 m below the hummock surface) (Fig. 10a, b,  d).
WTD drawdown due to simulated drainage not only affected modeled growing season GPP but also affected, and was affected by, the associated change in transpiration from vascular canopies.Deeper WTD x in drainage cycle 1 caused larger hydraulic gradients and greater lateral discharge, thereby deepening the WT with respect to that in the real-time simulation (Fig. 9a).Larger GPP in the drainage cycle 1 than in the real-time simulation throughout the growing seasons of 2004-2007 caused a greater vertical water loss through rapid transpiration from vascular canopies that  further contributed to this deepening of the WT (Fig. 9c).However, greater lateral water discharge in drainage cycle 2 caused by deeper WTD x did not deepen the modeled growing season WT much below that in cycle 1 (Fig. 9a).The larger lateral water loss through discharge in drainage cycle 2 than in cycle 1 was mostly offset by slower vertical wa-ter losses due to vascular plant water stress as indicated by smaller GPP in the drainage cycle 2 (Fig. 9c).The changing feedbacks between WTD, and GPP and plant water relations in ecosys also indicated the ability of the model to simulate hydrological self-regulation which is an important characteristic of peatland eco-hydrology (Dise, 2009).
Modeled growing season R e continued to increase with projected drainage-driven WTD drawdown (Fig. 9d, g).Reductions in modeled growing season R e from drainage cycle 1 to 2 during 2006-2009 indicated R e inhibition due to desiccation of near-surface peat layers and surface residues (Fig. 9d).Overall GPP increased more than R e with drainage-driven initial WTD drawdown that caused a small increase in modeled growing season NEP (Fig. 9 b, e).Continued drainage-driven WTD drawdown, however, caused declines in GPP particularly in the model years of 2008 and 2009 (Fig. 9c).This decrease in GPP was also accompanied by increased R e , thereby causing a decrease in NEP when the WT fell below a threshold of about ∼ 0.45 m from the hollow surface, particularly during the drier years (Fig. 9bg).This projected drainage effect on WTD and NEP may be transient.Long-term manipulation of WTD may produce different trajectories of WTD effects on C processes and plant water relations in northern boreal peatlands via vegetation adaptation and succession (Strack et al., 2006;Munir et al., 2014).

Modeling WTD effects on northern boreal peatland NEP
Hourly modeled and measured net ecosystem, soil and understory CO 2 fluxes, and modeled and EC-derived seasonal and annual NEP, GPP and R e showed that WTD drawdown raised both GPP and R e at the WPL .Similar increases in GPP and R e with WTD drawdown yielded no net effect of WTD drawdown on NEP at the WPL during 2004-2009 (Figs. 7-8) (Figs. 7-8).Our hypotheses, which outlined how coupled eco-hydrology and biogeochemistry algorithms in ecosys would simulate and explain the mechanisms of these WTD effects on R e , GPP and NEP at this boreal fen, are examined in details in the following Sect. of 5.1.1 to 5.1.4.

Hypothesis 1: increase in R e with WTD drawdown
Shallow WTD in ecosys caused a shallow aerobic zone above the WT and a thicker anaerobic zone below the WT.In the shallow aerobic zone, peat O 2 concentration [O 2s ] was well above the Michaelis-Menten constant for O 2 reduction (K m = 0.064 g m −3 ), and hence DOC oxidation and consequent microbial uptake and growth in ecosys were not much limited by [O 2s ] (Eqs.SA17a, SC14c).On the contrary, [O 2s ] in the thicker anaerobic zone below the WT was well below K m so that DOC oxidation was coupled with DOC reduction by anaerobic heterotrophic fermenters, which yielded much less energy (4.4 kJ g −1 C) than did DOC oxidation coupled with O 2 reduction (37.5 kJ g −1 C) (Fig. 1) (Eq.SA21).Lower energy yields in the thicker anaerobic zone caused slower microbial growth (Eq.SA25) and R h (Eq.SA13).Since the anaerobic zone in ecosys was thicker than the aerobic zone under a shallow WT, smaller R h modeled in the anaerobic zone contributed to reduced modeled soil respiration and hence R e which was corroborated by EC and chamber measurements at the WPL (Figs. 5-8) (Tables 1 and 2).WTD drawdown in ecosys caused peat pore drainage and increased θ g , thereby deepening the aerobic zone.It raised D g (Eq.SD44) and increased O 2 influxes into the peat (Fig. 5c) (Eqs.SD42-SD43).Increased O 2 influxes enhanced [O 2s ] and stimulated R h (Eqs.SA13, SA20), soil respiration and hence R e .Rapid mineralization of DON and DOP due to improved [O 2s ] under a deeper WT also raised aqueous concentrations of NH + 4 , NO − 3 and H 2 PO − 4 (Eqs.SC23a, c, e), which increased microbial nutrient availability, uptake (Eq.SA22) and growth (Eq.SA29) and enhanced R e further (Fig. 1) .Modeled R e stimulation by improved peat oxygenation due to WTD drawdown was corroborated well by EC and chamber measurements at the site (Figs.5-8) (Tables 1 and 2).Although the modeled rate of increase in R e with each 0.1 m of WTD drawdown was larger than the EC-derived rate, it is still comparable to rates reported for other similar peatlands (Table 2).Kotowska ( 2013) carried out chamber-based field measurements and laboratory incubation experiments in a moderately rich fen very close to our study site; she reported that a WTD drawdown-driven stimulation of aerobic microbial decomposition contributed to increased R e .Mäkiranta et al. (2009) also found rapid microbial decomposition in a Finish peatland due to a thicker aerobic zone and consequently larger amounts of decomposable organic matter exposed to aerobic oxidation.
Apart from WTD, peat warming in ecosys also increased rates of decomposition (Eq.SA1) through an Arrhenius function (Eq.SA6) and hence raised R h and R e (Figs. 4-7).The warming effect on decomposition in ecosys was also modified by WTD.For a similar warming, greater thermal diffusivity in peat with a deeper WT and consequent smaller water contents caused greater peat warming (Eqs.SD34, SD36).It enabled larger simulation of R e during warming periods in 2006 and 2008 with a deeper WT than in 2005 (Fig. 6).Increased stimulation of peat decomposition by warming under a deeper WT was also modeled by Grant et al. (2012) using the same model, ecosys, over a northern fen peatland at Wisconsin, USA, and by Ise et al. (2008) using a land surface scheme named ED-RAMS (Ecosystem Demography Model version 2 integrated with the Regional Atmospheric Modeling System) coupled with a soil biogeochemical model across several shallow and deep peat deposits in Manitoba, Canada.

Hypothesis 2: increase in GPP with WTD drawdown
Modeled WTD variations influenced GPP by controlling root and microbial O 2 availability, energy yields, root and mi-crobial growth and decomposition, rates of mineralization, and hence root nutrient (nitrogen) availability and uptake (Fig. 1).Wet soils under a shallow WT caused low O 2 diffusion (Fig. 5c) (Eqs.SD42-SD44) into the peat, and consequent low [O 2s ] meant that root O 2 demand had to be mostly met by [O 2r ].Ecosys inputs for root porosity (θ pr = 0.1) that governed O 2 transport through aerenchyma (Eq.SD45) and hence maintained [O 2r ] were not enough to meet the root O 2 demand in saturated soil by the two overstory tree PFTs, i.e., black spruce and tamarack, causing shallow root systems to be simulated in these two tree PFTs under shallow WTD (Sect.4.4).The understory shrub PFT (dwarf birch) had a higher root porosity (θ pr = 0.3) and hence had deeper rooting under a shallow WT than the two tree PFTs.Shallow rooting in the tree PFTs reduced root surface area for nutrient (nitrogen) uptake.Root nutrient (nitrogen) uptake (Eqs.SC23b, d, f) in all the PFTs was also constrained by low nitrogen availability due to smaller aqueous concentrations of NH + 4 and NO − 3 (Eqs.SC23a, c, e) resulting from slower mineralization (Eq.SA26) of DON (Eq.SA7) because of low [O 2s ] in the wet soils under a shallow WT (Fig. 1).Slower root growth and nutrient (nitrogen) uptake caused lower foliar σ N with respect to foliar σ C (Eq. SC11) that slowed the rates of carboxylation (Eq.SC6) and hence reduced vascular GPP (Eq.SC1) under a shallow WT.
WTD drawdown enhanced O 2 diffusion (Fig. 5c) (Eqs.SD42-SD44) and raised [O 2s ] so that root O 2 demand in all the three vascular PFTs was almost entirely met by [O 2s ].Consequently roots in all the PFTs could grow deeper, which increased the root surface for nutrient (nitrogen) uptake (Table 2) (Sect. 4.4).An increase in modeled rooting depth due to WTD drawdown was corroborated well by studies on same PFTs grown on similar peatlands very close to our study site (Table 2).Murphy et al. (2009) found a significant increase in tree fine-root production with WTD drawdown by 0.15-0.2m during a WTD manipulation study in a Finish peatland.Beside improved root growth, greater [O 2s ] under a deeper WT also enhanced rates of mineralization (Eq.SA26) of DON (Eq.SA7), which raised aqueous concentrations of NH + 4 and NO − 3 and hence facilitated root nutrient (nitrogen) availability and uptake (Fig. 1).Enhanced root nutrient (nitrogen) uptake increased foliar σ N with respect to foliar σ C (Eq. SC11), which hastened the rates of carboxylation (Eq.SC6) and hence raised vascular GPP (Eq.SC1) under a deeper WT.
The three modeled vascular PFT were predominantly N limited as indicated by mass-based modeled foliar N-to-P ratios that matched well with site-measured mass-based foliar N-to-P ratios (Table 2).Mass-based modeled and measured foliar N-to-P ratios in all the PFTs were less than 16 : 1, indicating that the vegetation at the WPL was N limited (Aerts and Chapin III, 1999).Since the modeled PFTs were predominantly N limited, increases in foliar N concentrations as a result of improved root growth and root nitrogen availability and uptake with WTD drawdown enhanced modeled carboxylation rates and hence modeled GPP.In a similar fen peatland close to our study site, Choi et al. (2007) found an increase in peat NO − 3 -N due to enhanced nitrogen mineralization and nitrification stimulated by a WTD drawdown, which improved foliar N status and hence increased radial tree growth of black spruce and tamarack trees (Table 2).Macdonald and Lieffers (1990) also found improved foliar N concentrations in black spruce and tamarack trees that enhanced net photosynthetic C assimilation rates by those tree species in a northern Alberta moderately rich fen (Table 2).The rates of increases in foliar N concentrations in black spruce and tamarack trees due to WTD drawdown as reported in those studies are comparable to those in our modeled outputs (Table 2).Although the modeled rate of increase in GPP with each 0.1 m of WTD drawdown was larger than the ECderived rate, it is still comparable to rates reported for other similar peatlands (Table 2).
5.1.3Hypothesis 3: microbial water stress on R e due to WT deepening below a threshold WTD When modeled WTD fell below a threshold of ∼ 0.3 m from the hollow surface (∼ 0.6 m below the hummock surface), the desiccation of the surface residue layer and near-surface shallow peat layers reduced microbial access to substrate for decomposition (Eq.SA15), which enabled the simulation of reduced R h in those layers.When reduction in surface residue and near-surface R h more than fully offset the increase in deeper R h , net ecosystem R h decreased.The offsetting effect on R h partly contributed to a simulated decrease in growing season R e (= R h + R a ) from 2008 to 2009 with a WTD drawdown that was corroborated by a similar decrease in ECderived R e (Fig. 7c).Greater reductions in R h in desiccated surface residue and near-surface peat layers also caused the reductions in growing season R e in drainage cycle 2 from those in cycle 1 during 2007-2009 in the projected drainage simulation (Fig. 9d).Similar to our study, Peichl et al. (2014) found reductions in R e when WTD fell below a threshold of 0.25 m from the peat surface in a Swedish fen, which could be partially attributed to a reduction in near-surface R h due to desiccation.Mettrop et al. (2014) in a controlled incubation experiment found that the rates of microbial respiration in a nutrient-rich Dutch fen had initially increased with peat drying and consequent improved aeration before excessive drying and consequent peat desiccation reduced microbial respiration efficiency, growth and biomass.While modeling the CO 2 exchange of a northern temperate bog using ecosys, Dimitrov et al. (2010a) found that a decrease in desiccated near-surface peat respiration partially offset increased deeper peat respiration when the WT deepened below a threshold of 0.6-0.7 m from the hummock surface.
www.biogeosciences.net/14/5507/2017/Biogeosciences, 14, 5507-5531, 2017 5.1.4Hypothesis 4: plant water stress on GPP due to WT deepening below a threshold WTD The desiccation of near-surface peat layers with WTD drawdown below a threshold level also caused reduced plant water uptake from those layers which were colonized by most of the vascular root systems and all of the belowground biomasses of nonvascular mosses (Eqs.SD9-SD29).When WTD fell below ∼ 0.1 m from the hollow surface (∼ 0.4 m below the hummock surface), vertical recharge through capillary rise from the WT was not adequate to maintain nearsurface peat moisture.It reduced peat ψ s and raised peat s (Eq.SB9), which suppressed root and moss water uptake (U w ) (Eq.SB6) from desiccated near-surface peat layers (Fig. 1).Since moss U w entirely depended upon moisture supply from the near-surface layers, a reduction in U w from the desiccation of these layers caused a reduction in moss ψ c and hence moss GPP (Fig. 1) (Eqs.SC1, SC4) (Mezbahuddin et al., 2016).A reduction in root U w from desiccated near-surface layers, however, was offset by increased root U w (Eq.SB6) from deeper wetter layers, which had higher ψ s and lower s , due to deeper root growth facilitated by enhanced aeration.It enabled the vascular PFTs in ecosys to sustain ψ c , canopy ψ t (Eq.SB4), g c (Eqs.SB2, SC4) and hence to sustain increased GPP (Eq.SC1) due to higher root nutrient (nitrogen) availability and uptake (Fig. 1) (Mezbahuddin et al., 2016).Increased vascular GPP and consequent greater vascular plant growth further imposed limitations of water, nutrient and light on the modeled nonvascular PFTs due to interspecific competition and greater shading from the overstory vascular PFTs.However, an increases in vascular GPP due to enhanced plant nutrient (nitrogen) status more than fully offset the suppression in moss GPP due to moss drying and greater shading and competition from the overstory, thereby causing a net increase in modeled GPP with WTD drawdown (Figs.7b and 10b, c).This simulation of increased vascular dominance over moss with a deepening of the WT was corroborated by several WTD manipulation studies (e.g., Moore et al., 2006, Munir et al., 2014) in similar peatlands in Alberta that reported increased tree, shrub and herb growths over mosses with WTD drawdown.However, an increase in projected vascular GPP eventually plateaued, and GPP started to decline when the WT fell below ∼ 0.6 m from the hollow surface (∼ 0.9 m below the hummock surface) (Figs.10c, e).This was because deeper root U w (Eq.SB6) could no longer offset the suppression of near-surface root U w when the WT fell below threshold WTD, thereby causing lower ψ c , ψ t (Eq.SB4), g c (Eqs.SB2, SC4) and slower CO 2 fixation (Eq.SC6) (Fig. 1).
These threshold WTD effects on modeled vascular and nonvascular plant water relations were validated well by testing modeled vs. site-measured hourly energy fluxes (latent and sensible heat) and Bowen ratios and modeled vs. sitemeasured daily peat moisture at different depths throughout 2004-2009 as described in Mezbahuddin et al. (2016).Ri-utta et al. (2007) measured a reduction in moss productivity due to water limitation when WTD fell below ∼ 0.15 m from the surface in a Finish fen peatland.However, they reported a sustained vascular GPP during that period, indicating no vascular water stress (Riutta et al., 2007).Peichl et al. (2014) reported a reduction in moss GPP due to moss drying caused by an insufficient moisture supply through capillary rise when WTD fell below 0.25 m from the surface in a Swedish fen.Reductions in moss GPP due to decreased moss canopy water potentials were also modeled by Dimitrov et al. (2011) using the same model, ecosys, when WTD fell below 0.3 m from the hummock surface of a Canadian temperate bog.They, however, found no vascular plant water stress and hence no reduction in vascular GPP during that period.Similarly, Kuiper et al. (2014) found reductions in moss productivity with peat drying while vascular productivity was sustained in a simulated drought experiment on a Danish peat.
Continued deepening of the WT can also cause vascular plant water stress and hence reductions in vascular GPP as projected in our drainage simulation (Fig. 10c, e).This can also be corroborated by field measurements across various northern boreal fen peatlands in Canada and Sweden.Sonnentag et al. ( 2010) found a reduction in g c of a canopy that included tamarack and dwarf birch and a consequent decline in GPP when the WT fell below 0.3 m from the ridge surface at a fen peatland in Saskatchewan.Peichl et al. (2014) also found a reduction in vascular GPP due to plant water stress when WTD fell below 0.25 m from the surface in a Swedish fen.The WTD threshold for reductions in vascular GPP in those two field studies were shallower than that in our modeled projection, i.e., ∼ 0.6 m from the hollow surface (∼ 0.9 m below the hummock surface) (Fig. 10a, c, e), thereby indicating different vertical rooting patterns determined by specific interactions between hydrologic properties and rooting.Lafleur et al. (2005) and Schwärzel et al. (2006) found much deeper WTD thresholds for reductions in vascular transpiration that negatively affected vascular GPP over a Canadian pristine bog and a German drained fen, respectively.Those WTD thresholds were ∼ 0.65 and ∼ 0.9 m below the surface for the pristine and drained peatland, respectively, further indicating the importance of root-hydrology interactions and the resultant root adaptations, growth and uptake in determining WTD effects on vascular GPP across peatlands.

Divergences between modeled and EC-derived
annual GPP, R e and NEP Modeled seasonal and annual GPP and R e were consistently larger than EC-derived estimates of GPP and R e during 2004-2009 (Fig. 7b, c, f, g) (Fig. 7b, c, f, g).Although modeled annual GPP and R e were larger than the EC-derived estimates, modeled annual NEP was consistently lower than the EC gap-filled annual NEP (Fig. 7e).Modeled annual NEP was smaller than EC gap-filled NEP since the divergence between modeled and EC-derived R e was larger than the deviation between modeled and EC-derived GPP (Fig. 7f, g).Modeled R e were larger than EC-derived R e mainly due to the presence of gapfilled nighttime CO 2 fluxes (= R e ) in EC-derived R e , which were smaller than corresponding modeled values.This was apparent in negative intercepts that resulted from regressions of modeled vs. gap-filled net CO 2 fluxes (Table S1 in the Supplement).Gap-filled R e fluxes were calculated from soil temperature (T s ) at a shallow depth (0.05 m) (Sect.3.2.2).During nighttime and in the winter, peat at this shallow depth rapidly cooled down, which produced smaller nighttime gapfilled CO 2 fluxes empirically derived from lower T s (Figs. 3,5 and 6).By contrast, corresponding modeled CO 2 effluxes were affected by not only the cooler shallow peat layers but also the warmer deeper peat layers and thus were larger than the gap-filled fluxes (e.g., Figs. 5a,6h,i).Like modeled CO 2 effluxes, chamber-measured soil and understory CO 2 effluxes in cooler nights did not decline as rapidly as did the corresponding gap-filled net ecosystem CO 2 fluxes as night progressed, which further indicated the likely contribution of a gap-filling artifact to CO 2 effluxes that were smaller than corresponding modeled fluxes (e.g., Figs.5b vs. a, 6j-k vs. g-h).Systematic uncertainties embedded in EC methodology could also have contributed to EC-derived annual and growing season R e which was smaller than the modeled values (Fig. 7c, g).The major uncertainty in the EC methodology is the possible underestimation of nighttime EC CO 2 flux measurements due to poor turbulent mixing under stable air conditions (Goulden et al., 1997;Miller et al., 2004).By contrast, modeled biological production of CO 2 by plant and microbial respiration was independent of turbulent mixing, which could thus contribute to modeled R e that was larger than EC-derived estimates.
Complete energy balance closure in the model as opposed to incomplete (∼ 75 %) energy balance closure in EC measurements would also give rise to modeled evapotranspiration and GPP values that were larger than EC-derived estimates (Fig. 7b, f) (Mezbahuddin et al., 2016).Modeled GPP influenced modeled R e through root exudation and litter fall (Fig. 1).Therefore, modeled GPP that was larger than ECderived estimates would have further contributed to modeled growing season and annual R e that was larger than ECderived estimates (Fig. 7c, g).
Modeled GPP and hence R e larger than EC-derived estimates could also be a result of uncertainties in model inputs for soil organic N, N deposition, N 2 fixation and any other sources of N inputs into the modeled ecosystem.In ecosys, plant productivity is governed by foliar N status, which is constrained by root N availability and uptake.Our input for organic N into each modeled peat layer was measured for corresponding depth at the site (Fig. 2).To simulate N deposition, background wet deposition rates of 0.5 mg ammonium-N, and 0.25 mg nitrate-N per liter of precipi-tation were used as model inputs.However, from visual field observations, it was evident that there was a significant amount of nutrient inflow with the lateral water influxes into this fen peatland from the surrounding upland forests; this inflow was not quantified.To mimic this lateral nutrient inflow, we doubled the background wet deposition of NH + 4 and NO − 3 as reported for the area and used these as surrogates of lateral nutrient inflows into the modeled ecosystem.Ecosys also included an N 2 fixing algorithm which simulated symbiotic N 2 fixation in moss canopies that was reported for the boreal forests (Sect.3.2.3).We tested the adequacy of these N inputs into the model by comparing modeled leaf N concentrations against those measured in the field.The modeled foliar N concentrations for all the vascular PFTs compared well to site measurements (Table 2).To further examine the contribution of uncertainty due to model inputs towards the divergence between modeled and EC-derived seasonal and annual GPP and R e , we performed a sensitivity test where we had a parallel run without doubling the background N wet deposition rates in the model, hence simulating no lateral N influx into the modeled ecosystem.Unlike the run with lateral N inflow, the parallel run without lateral N inflow simulated GPP and R e , which were very close to the EC-derived estimates.However, the regressions between hourly modeled net CO 2 fluxes from the parallel run and EC-measured hourly net CO 2 fluxes gave slopes of ∼ 0.8, indicating under-simulation of the EC-measured fluxes in the parallel run with no lateral N inflow.

Conclusions
Our modeling study showed that, when adequately coupled to algorithms of a process-based ecosystem model, the existing knowledge of peatland eco-hydrological and biogeochemical processes could explain underlying mechanisms that governed WTD effects on net ecosystem CO 2 exchange of a boreal fen.Testing of our hypotheses against ECmeasured net CO 2 fluxes, automated chamber fluxes and other biometric measurements at the site revealed that a drierweather-driven WTD drawdown at this boreal fen raised both R e and GPP due to improved aeration that facilitated (1) microbial and root O 2 availability, energy yields, growth and decomposition, which raised microbial and root respiration; and (2) rapid nitrogen mineralization and consequently increased root nitrogen availability and uptake that improved leaf nitrogen status and hence raised carboxylation (Figs. 1, 7, 8) (Table 2).Similar increases in R e and GPP with WTD drawdown to a certain depth caused no net WTD drawdown effect on NEP (Fig. 8).Modeled drainage projection, however, showed that further WTD drawdown caused by either drainage or climate-change-induced drying would produce plant water stress and reduce GPP and hence NEP of this boreal fen (Figs. 9,10).This study further reconciled and mechanistically explained the WTD effects on seasonal and www.biogeosciences.net/14/5507/2017/Biogeosciences, 14, 5507-5531, 2017 annual GPP, R e and hence NEP of this boreal fen, on which there was previously speculation on the basis of EC-derived estimates.However, although modeled CO 2 fluxes were validated well by the EC and chamber-measured net CO 2 fluxes, modeled values of annual and seasonal GPP and R e were consistently larger than the EC-derived estimates (Table 1) (Figs. 7, 8) (Sect.5.2).These discrepancies between modeled and EC-derived GPP and R e estimates also created an opportunity of value addition by using more robust processbased estimates while forecasting WTD effects on these peatland C balance components in addition to empirically modeled EC-derived estimates that might not include some of the above-discussed offsetting feedbacks in peatland ecohydrology and biogeochemistry.
The model algorithms that were used in this study represented coupled feedbacks among ecosystem processes that governed carbon, water, energy and nutrient (N, P) cycling in a peatland ecosystem.These feedbacks were thus not parameterized for this particular boreal fen peatland site.Instead, the modeled boreal fen was simulated from peatlandspecific model inputs for weather, soil and vegetation properties that had physical meaning and were quantifiable at the site (Fig. 2).These modeled process-level interactions were also validated by corroborating modeled outputs against site measurements.By contrast, most of the current peatland C models use scalar functions to represent these feedbacks, and so those model algorithms have to be parameterized for each peatland site.Therefore, the modeling approach as described in this study should be more robust than the scalar feedback approach while assessing WTD effects on peatland C balance under contrasting peat types, climates and hydrology or under unknown future climates.These process-level feedbacks are also scalable once the peatlands of interest are defined within the modeled landscapes by scalable model inputs for weather; peat hydrological, physical, and biological properties; and plant functional types (Sect. 3.3.3).Current global land surface models either lack or have very poor representation of these feedbacks, which is thus far limiting our largescale predictive capacity of WTD effects on boreal peatland C stocks.This modeling exercise would thus provide valuable information to improve the representation of these feedbacks into next-generation land surface models.Therefore, the insights gained from this modeling study should be a significant contribution to our understanding and apprehension of how peatlands would behave with changing hydrology under future drier and warmer climates.

Figure 1 .
Figure 1.Schematic diagram of ecosys algorithms representing coupled key eco-physiological and biogeochemical (aerobic and anaerobic) processes and plant water relations of a typical boreal fen peatland ecosystem that are affected by WTD fluctuation.a , c and s : atmospheric, canopy and soil water potentials; r : vascular root or nonvascular belowground water potential; r c : canopy stomatal resistance (1/g c ); s : soil hydraulic resistance; r : hydraulic resistance to water flow through plants; OM: organic matter; DOC, DON, DOP: dissolved organic carbon, nitrogen and phosphorus; POC: particulate organic C; and POM: particulate organic matter.

Figure 2 .
Figure 2. Layout for ecosys model run to represent biological, chemical and hydrological characteristics of a western Canadian fen peatland.Figure is not drawn to scale.D humm : depth to the bottom of a layer from the hummock surface; D holl : depth to the bottom of a layer from the hollow surface; TOC: total organic C (Flanagan and Syed, 2011); TN: total nitrogen (Flanagan and Syed, 2011); TP: total phosphorus (Flanagan and Syed, 2011); CEC: cation exchange capacity (Rippy and Nelson, 2007) -the value for pH was obtained from Syed et al. (2006); WTD x : external reference water table depth representing average water table depth of the adjacent ecosystem; L t : distance from modeled grid cells to the adjacent watershed over which lateral discharge/recharge occurs.
Figure 2. Layout for ecosys model run to represent biological, chemical and hydrological characteristics of a western Canadian fen peatland.Figure is not drawn to scale.D humm : depth to the bottom of a layer from the hummock surface; D holl : depth to the bottom of a layer from the hollow surface; TOC: total organic C (Flanagan and Syed, 2011); TN: total nitrogen (Flanagan and Syed, 2011); TP: total phosphorus (Flanagan and Syed, 2011); CEC: cation exchange capacity (Rippy and Nelson, 2007) -the value for pH was obtained from Syed et al. (2006); WTD x : external reference water table depth representing average water table depth of the adjacent ecosystem; L t : distance from modeled grid cells to the adjacent watershed over which lateral discharge/recharge occurs.

4. 3
WTD effects on diurnal net ecosystem CO 2 exchange WTD variation can affect diurnal net CO 2 exchange by affecting peat O 2 status and consequently root and microbial O 2 and nutrient availability, growth and uptake, thereby in-fluencing CO 2 fixation and respiration.Ecosys simulated WTD effects on diurnal net CO 2 exchange at the WPL well over three 10-day periods with comparable weather conditions (radiation and air temperature) during the late growing seasons (August) of2005 , 2006  and 2008 (Fig. 4) (Fig. 4).A WTD drawdown from August 2005 to August 2006 in ecosys caused a reduction in peat water contents and a consequent increase in O 2 influxes from atmosphere into the peat that caused an increase in modeled soil CO 2 effluxes(Figs.4c,  5c).This stimulation of soil respiration was corroborated by modeled vs. chamber-measured(Cai et al., 2010) nighttime soil CO 2 fluxes and understory aboveground autotrophic respiration (R a ) in August 2006 with a deeper WTD which were larger than those in 2005 with a shallower WTD (Fig.5b).
Figure3.(a, c, e, g, i, k) Three-day moving averages of modeled and EC gap-filled(Flanagan and Syed, 2011) NEP and (b, d, f, h, j, l)   hourly modeled and half-hourly measured(Syed et al., 2006;Cai et al., 2010;Long et al., 2010;Flanagan and Syed, 2011) WTD from 2004 to 2009 at a western Canadian fen peatland.A positive NEP means that the ecosystem is a carbon sink, and a negative NEP means that the ecosystem is a carbon source.A negative WTD represents a depth below hummock/hollow surface, and a positive WTD represents a depth above hummock/hollow surface.

Figure 4 .
Figure 4. Half-hourly measured (a) incoming shortwave radiation and (b) T a and (c) hourly modeled and half-hourly measured(Syed et al., 2006;Cai et al.;2010, Long et al.; 2010, Flanagan andSyed, 2011) WTD during August of 2005, 2006 and 2008 at a western Canadian fen peatland.A negative WTD represents a depth below hummock/hollow surface, and a positive WTD represents a depth above hummock/hollow surface.Grey arrows indicate nights with similar temperatures.

Figure 5 .
Figure 5. (a) Half-hourly EC-measured and gap-filled(Flanagan and Syed, 2011) and hourly modeled ecosystem net CO 2 fluxes, (b) halfhourly automated chamber-measured(Cai et al., 2010) and hourly modeled understory and soil CO 2 fluxes, and (c) hourly modeled soil CO 2 and O 2 fluxes during August of 2005, 2006 and 2008 at a western Canadian fen peatland.No chamber CO 2 flux measurement was available for 2008.Bars represent standard errors of means of chamber CO 2 fluxes (n = 9).A negative flux represents an upward flux or a flux out of the ecosystem, and a positive flux represents a downward flux or a flux into the ecosystem.Grey arrows indicate nights with similar temperatures in the same way as in Fig. 4.
Figure 6.(a-c) Half-hourly observed T a ;(d-f) hourly modeled and half-hourly observed(Syed et al., 2006;Cai et al., 2010;Long et al., 2010;Flanagan and Syed, 2011) WTD; (g-i) half-hourly EC-measured and gap-filled(Flanagan and Syed, 2011) and hourly modeled ecosystem net CO 2 fluxes; (j-l) half-hourly automated chamber-measured(Cai et al., 2010) and hourly modeled understory and soil CO 2 fluxes during July-August of 2005, 2006 and 2008 at a western Canadian fen peatland.No chamber CO 2 flux measurement was available for 2008.Bars represent standard errors of means of chamber CO 2 fluxes (n = 9).A negative flux represents an upward flux or a flux out of the ecosystem, and a positive flux represents a downward flux or a flux into the ecosystem.A negative WTD represents a depth below hummock/hollow surface, and a positive WTD represents a depth above hummock/hollow surface Figure 6.(a-c) Half-hourly observed T a ;(d-f) hourly modeled and half-hourly observed(Syed et al., 2006;Cai et al., 2010;Long et al., 2010;Flanagan and Syed, 2011) WTD; (g-i) half-hourly EC-measured and gap-filled(Flanagan and Syed, 2011) and hourly modeled ecosystem net CO 2 fluxes; (j-l) half-hourly automated chamber-measured(Cai et al., 2010) and hourly modeled understory and soil CO 2 fluxes during July-August of 2005, 2006 and 2008 at a western Canadian fen peatland.No chamber CO 2 flux measurement was available for 2008.Bars represent standard errors of means of chamber CO 2 fluxes (n = 9).A negative flux represents an upward flux or a flux out of the ecosystem, and a positive flux represents a downward flux or a flux into the ecosystem.A negative WTD represents a depth below hummock/hollow surface, and a positive WTD represents a depth above hummock/hollow surface

Figure 7 .
Figure 7. Modeled and EC-derived (Flanagan and Syed, 2011) growing season (May-August) sums of (a) NEP, (b) GPP and (c) R e during 2004-2009; (d) observed mean growing season T a and measured and modeled average growing season WTD during 2004-2009; modeled and EC-derived(Flanagan and Syed, 2011) annual sums of (e) NEP, (f) GPP and (g) R e during 2004-2008 and (h) observed mean annual T a and measured and modeled average WTD during ice-free periods (May-October) of 2004-2008 at a western Canadian fen peatland.A negative WTD represents a depth below hollow surface, and a positive WTD represents a depth above hollow surface.A positive NEP means that the ecosystem is a carbon sink.Annual modeled vs. EC gap-filled NEP, GPP and R e estimates for 2009 were not compared due to the lack of flux measurements from September to December in that year.
Modeled Eddy-covariance-derived/ biometrically measured at the site a Values from other studies in similar peatlands Growing season (May to August) mean WTD drawdown from 2004 to 2009 from 0.3 m below the hummock surface (at the hollow surface) in 2004 to 0.65 m below the hummock surface (0.35 below the hollow surface) in 2009 from 0.32 m below the hummock surface (0.02 m below the hollow surface) in 2004 to 0.62 m below the hummock surface (0.32 m below the hollow surface) in 2009 Rate of increase in R 45.1 g N kg −1 C with the growing season WTD drawdown from 0.3 m below the hummock surface in 2004 to 0.65 m below the hummock surface in 2009 Increase in maximum rooting depth with WTD drawdown black spruce from 0.35 to 0.65 m below the hummock surface with the growing season WTD drawdown from 0.3 m below the hummock surface in 2004 to 0.65 m below the hummock surface in 2009 from 0.1 to 0.6 m below hummock surface with a WTD drawdown from 0.1 to 0.7 m below hummock surface f tamarack from 0.35 to 0.65 m below the hummock surface with the growing season WTD drawdown from 0.3 m below the hummock surface in 2004 to 0.65 m below the hummock surface in 2009 from 0.2 to 0.6 m below hummock surface with a WTD drawdown from 0.1 to 0.6 m below hummock surface f a Syed et al. (2006), Flanagan and Syed (2011).b Peichl et al. (2012) for a Swedish fen.c Ballantyne et al. (2014) for a Michigan fen peatland complex that has similar peat and plant functional types as our study site.d Choi et al. (2007) for a central Alberta fen peatland located ∼ 350 km to the southwest of the study site.e Macdonald and Lieffers (1990) for a northern Alberta fen peatland located ∼ 250 km to the northwest of the study site.f Lieffers and Rothwell (1987) for a northern Alberta fen peatland located ∼ 250 km to the northwest of the study site.a, d, e unit conversions for leaf nutrient values were done by assuming organic C as 50 % of total leaf dry matter.

m - 2 )Figure 8 .
Figure 8. Regressions (P < 0.001) of growing season (May-August) sums of modeled and EC-derived (Flanagan and Syed, 2011) (a) NEP, (b) GPP and (c) R e on growing season averages of modeled and observed WTD during 2004-2009; regressions (P < 0.001) of annual sums of modeled and EC-derived (Flanagan and Syed, 2011) (d) NEP, (e) GPP and (f) R e on average modeled and measured WTD during ice-free periods (May-October) of 2004-2008 at a western Canadian fen peatland.A negative WTD represents a depth below hollow surface, and a positive WTD represents a depth above hollow surface.A positive NEP means that the ecosystem is a carbon sink.

Figure 9 .Figure 10 .
Figure 9. (a) Observed, simulated (real-time simulation) and projected (drainage cycles 1 and 2) average growing season (May-August) WTD; EC-derived, simulated and projected growing season sums of (b) NEP, (c) GPP and (d) R e ; and regressions (P < 0.001) of simulated and projected sums of (e) NEP, (f) GPP and (g) R e on simulated and projected average growing season WTD during 2004-2009 at a western Canadian fen peatland.A negative WTD represents a depth below hollow surface, and a positive WTD represents a depth above hollow surface.A positive NEP means that the ecosystem is a C sink.
) in the wetter peat layers above the WT.Lower θ g reduces O 2 diffusivity in the gaseous phase (D g ) (Eq.SD44) and gaseous O 2 transport (Eqs.SD41-SD42) in these layers.Peat layers below the WT have zero θ g that prevents gaseous O 2 transport in these layers.So, under a shallow

Table 1 .
Statistics from regressions between hourly modeled and measured net CO 2 fluxes from 2004 to 2009 at a western Canadian fen peatland.(a) Regressions of modeled vs. EC-measured (recorded at u * > 0.15 m s −1 ) net ecosystem CO 2 fluxes over the whole years of 2004-2008 * Regressions of modeled vs. EC-measured (recorded at u * > 0.15 m s −1 ) net ecosystem CO 2 fluxes over the growing seasons (May-August) of 2004-2009 Regressions of modeled vs. measured chamber net CO 2 fluxes (understory vegetation and soil CO 2 fluxes) over ice-free periods (May-October) of 2005-2006 Table 1).Modeled NEPs throughout the winters of most of the years were www.biogeosciences.net/14/5507/2017/Biogeosciences, 14, 5507-5531, 2017 WTD: water table depth below the hummock surface; (a, b) from simple linear regressions of modeled on measured (± standard errors) and R 2 : coefficient of determination; RMSE: root mean square for errors from simple linear regressions of measured on simulated; RMSRE: root mean square for random errors in measurements -RMSRE for eddy covariance (EC) measurements was estimated by inputting EC CO 2 fluxes recorded at u * (friction velocity) > 0.15 m s −1 into algorithms for the estimation of random errors due to EC CO 2 measurements developed for forests by Richardson et al. (2006); * whole-year modeled vs. EC net CO 2 flux regression for 2009 could not be done due to the lack of flux measurements from September to December in that year.
WTD during August of 2005, 2006 and 2008 at a western Canadian fen peatland.A negative WTD represents a depth below hummock/hollow surface, and a positive WTD represents a depth above hummock/hollow surface.Grey arrows indicate nights with similar temperatures.

Table 2 .
Effects of WTD drawdown on components of ecosystem carbon and nutrient cycles of a western Canadian fen peatland.